Objectives

This lab exercise explores a benchmark problem of 2-D steady-state heat conduction through a square-shaped domain heated uniformly by a source. By removing one quarter of the square, an "L" shaped domain is created. The concave corner then introduces a singularity that degrades the solution convergence rate. The uniform mesh refinement study for this problem impliments the natural coordinate FE basis (triangular elements) and the tensor product FE basis (quadrilateral elements without quadrature).

Problem Statement

Problem Statement

This lab considers 2-D steady-state heat conduction through an "L" shaped region. This region has a thermal conductivity of unity (k = 1) as well as a unity source term (s = 1). Dirichlet conditions (Tb = 0) are imposed along the entire boundary. The domain consists of a 2x2 unit square with coordinates running from -1 to 1 on each axis. Then, to create the "L" shape, one quarter of the domain is forced to a fixed temperature of zero, equivalent to the Dirichlet BC implimentation.

The resulting characteristic equation (Laplace's equation) is well known and has been mentioned in previous labs. Also, as previously discussed, the error convergence rate is known to have a slope dominated by 2*min[1, r-1] for a linear FE basis. Here, the "r" term depends upon data non-smoothness. For the "L" shaped domain, the concave corner is expected to introduce a boundary singularity that will control the error convergence rate. This will be benchmarked against an equivalent full square domain which should exhibit the expected convergence slope of 2.

Part A utilizes the natural coordinate system to create the domain mesh. Thus, a linear triangle mesh is created via the Delauney method. The MATLAB program code shown below impliments this mesh proceedure as based upon the MATLAB codes given in labs 7 and 8. Part B impliments the tensor product coordinate system to create the linear domain mesh. Thus, quadrilaterals are created (without quadrature) using an instructor supplied MATLAB program.

# PART A

## Natural Coordinate System(Linear Triangle Elements)

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #9a                                                                         %
% 2-D, steady-state, conduction using triangular, elements (k=1).  A square and "L"      %
% shape region are created by calling either the do_boxshape2D_dlist or                  %
% do_Lshape2D_dlist functions.  The domain is defined on a 2x2 unit square, with         %
% coordinates running from -1 to 1 on each axis.  k*Laplacian(q) - s = 0 is solved,      %
% with k=1 and s=1 on the interior.  A Dirichlet condition (q=0) is imposed along the    %
% entire boundary.                                                                       %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

Mesh = [8 16 32 64 128];

for i = 1:size(Mesh,2)

% Declare global variables.  These variables must be made available as
% global variables for the assembly routine to work.

elwide = Mesh(i);	% number of elements along the width of the grid
eltall = Mesh(i);	% number of elements along the height of the grid

% Construct the grid seeds.  Be sure that X1seed and X2seed both include 0.
X1seed=linspace(-1,1,elwide+1);
X2seed=linspace(1,-1,eltall+1);

% construct the nodal coordinate lists
[X1 X2] = meshgrid(X1seed, X2seed);
nnodes = prod(size(X1));

% construct the element connectivity table
elcon = delaunay(X2,X1);

% save element mesh size
elements(i,1) = size(elcon,1);

% Scan through all nodes and select the ones for which Dirichlet
% conditions will be enforced. Also, select the nodes outside the
% computational domain for blanking.  This will be enforced in exactly
% the same manner as Dirichlet B.C's.
diri_list = do_Lshape2D_dlist(X1,X2);

% calculate element zeta's and determinants.
[ZETAe,DETe]=do_ZETAe_2Dt(X1,X2,elcon);

% generate element diffusion matrices, [DIFF]e
DIFFe = do_DIFFe_2Dt(ZETAe);

% assemble global diffusion matrix, [DIFF].
DIFF = asjac([],[],[],-1,DIFFe,[],elcon,DETe,nnodes);

% assemble the residual
B200L = 1/24 * [2 1 1; 1 2 1; 1 1 2];
SRC = ones(size(X1));
b = asres([],[],[],1,B200L,SRC,elcon,DETe,nnodes);

% apply the BC's/blanking
[DIFF_diri b_diri] = do_diri(DIFF,b,diri_list);

% solve the linear system
Q = DIFF_diri\b_diri;

% save Qmax and Enorm data
Qmax(i,1) = max(Q);
enorm(i,1) = 1/2 * (Q' * DIFF * Q);

end

delta_Qmax = Qmax(2:end) - Qmax(1:end-1);
slope_Qmax = log10(delta_Qmax(1:end-1) ./ delta_Qmax(2:end)) / log10(2);
delta_enorm =  enorm(2:end) - enorm(1:end-1);
slope_enorm = log10(delta_enorm(1:end-1) ./ delta_enorm(2:end)) / log10(2);

fprintf('Elements   Qmax  (Qmax)change  (Qmax)slope    ||Q||      change||Q||     Slope(Q)\n\n')
fprintf('%6.0f  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g\n',ANSWER');

% plot the solution
figure(1)
trisurf(elcon,X1,X2,Q); view(0,0);
figure(2)
trisurf(elcon,X1,X2,Q); view(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

# PART B

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #9b                                                                         %
% A square and "L" shape region are created by calling either the do_boxshape2D_dlist or %
% do_Lshape2D_dlist functions.  The domain is defined on a 2x2 unit square, with         %
% coordinates running from -1 to 1 on each axis.  k*Laplacian(q) - s = 0 is solved,      %
% with k=1 and s=1 on the interior.  A Dirichlet condition (q=0) is imposed along the    %
% entire boundary.                                                                       %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

Mesh = [8 16 32 64 128];

for i = 1:size(Mesh,2)

% Declare global variables.  These variables must be made available as
% global variables for the assembly routine to work.

elwide = Mesh(i);	% number of elements along the width of the grid
eltall = Mesh(i);	% number of elements along the height of the grid

% Construct the grid seeds.  Be sure that X1seed and X2seed both include 0.
X1seed=linspace(-1,1,elwide+1);
X2seed=linspace(1,-1,eltall+1);

% construct the nodal coordinate lists
[X1 X2] = meshgrid(X1seed, X2seed);
nnodes = prod(size(X1));

% construct the element connectivity table

% save element mesh size
elements(i,1) = size(elcon,1);

% Scan through all nodes and select the ones for which Dirichlet
% conditions will be enforced. Also, select the nodes outside the
% computational domain for blanking.  This will be enforced in exactly
% the same manner as Dirichlet B.C's.
diri_list = do_Lshape2D_dlist(X1,X2);

% calculate d(etai)/d(xj) and determinants for each element
[DETe,DETeI,ETAe]=do_ETAe_q0(X1,X2,elcon);

% calculate diffusion matrix for each element
DIFFe = do_DIFFe_q0(ETAe);

% assemble the system diffusion matrix
DIFF = asjac([],[],[],-1,DIFFe,[],elcon,DETe,nnodes);

% assemble the residual
B200L = 1/9 * [4 2 1 2; 2 4 2 1; 1 2 4 2; 2 1 2 4];
SRC = ones(size(X1));
b = asres([],[],[],1,B200L,SRC,elcon,DETe,nnodes);

% apply the BC's/blanking
[DIFF_diri b_diri] = do_diri(DIFF,b,diri_list);

% solve the linear system
Q = DIFF_diri\b_diri;

% save Qmax and Enorm data
Qmax(i,1) = max(Q);
enorm(i,1) = 1/2 * (Q' * DIFF * Q);

end

delta_Qmax = Qmax(2:end) - Qmax(1:end-1);
slope_Qmax = log10(delta_Qmax(1:end-1) ./ delta_Qmax(2:end)) / log10(2);
delta_enorm =  enorm(2:end) - enorm(1:end-1);
slope_enorm = log10(delta_enorm(1:end-1) ./ delta_enorm(2:end)) / log10(2);

fprintf('Elements   Qmax  (Qmax)change  (Qmax)slope    ||Q||      change||Q||     Slope(Q)\n\n')
fprintf('%6.0f  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g\n',ANSWER');

% reshape Q from a column matrix to an array of the same shape as X1 for
% plotting purposes.
Q_filled = zeros(size(X1));
Q_filled(:) = Q;

surfc(X1,X2,Q_filled);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```
Discussion of Results

Discussion of Results

# PART A

## Natural Coordinate System(Linear Triangle Elements) Figure 1A: Temperature solution planform view for square-shaped region using triangular mesh (k=1). Figure 1B: Temperature solution angle view for square-shaped region using triangular mesh (k=1).

Figures 1A-B detail the steady-state temperature solution for the full square domain using triangular elements. Table 1 lists the convergence results for the uniform mesh refinement study. As expected the results indicate a convergence slope of 2, indicating that the "k+1-1" FE basis term dominates the convergence rate.

Table 1: Error convergence results for square-shaped region using triangular mesh (k=1).
```==================================================================================
Elements      Qmax   (Qmax)change  (Qmax)slope     ||Q||   change||Q||  Slope(Q)
__________________________________________________________________________________
128    0.289331              0            0  0.267858             0         0
512    0.293794     0.00446326            0  0.277844    0.00998513         0
2048    0.294461    0.000667315      2.74166  0.280303    0.00245964   2.02133
8192    0.294705     0.00024389      1.45214  0.280941   0.000637795   1.94728
32768    0.294668  -3.71168e-005      2.71608  0.281101   0.000159556   1.99903
==================================================================================
``` Figure 2A: Temperature solution planform view for "L"-shaped region using triangular mesh (k=1). Figure 2B: Temperature solution angle view for "L"-shaped region using triangular mesh (k=1).

Figures 2A-B detail the steady-state temperature solution for the "L" shaped domain. Table 2 lists the convergence results for the uniform mesh refinement study. Now, a degredation of the convergence rate occurs with an indicated slope near 5/3. This would indicate that the singularity induced by the concave corner produces non-smoothness that degrades the solution process. Thus, the "r-1" term controls the convergence slope.

Table 2: Error convergence results for "L"-shaped region using triangular mesh (k=1).
```==================================================================================
Elements      Qmax  (Qmax)change  (Qmax)slope      ||Q||   change||Q||  Slope(Q)
__________________________________________________________________________________
128    0.138975             0            0   0.095661             0         0
512    0.146308    0.00733293            0   0.103536   0.007874550         0
2048    0.148214    0.00190621      1.94368   0.105930   0.002394370   1.71755
8192    0.148998   0.000783458      1.28278   0.106680   0.000749491   1.67566
32768    0.149254   0.000255973      1.61386   0.106919   0.000239033   1.64870
==================================================================================
```

# PART B Figure 3A: Temperature solution planform view for square-shaped region using quadrilateral mesh (k=1). Figure 3B: Temperature solution angle view for square-shaped region using quadrilateral mesh (k=1).

Figures 3A-B detail the steady-state temperature solution for the full square domain using quadilateral elements. Table 3 lists the convergence results for the uniform mesh refinement study. As expected the results indicate a convergence slope of 2, indicating that the "k+1-1" FE basis term dominates the convergence rate. Upon comparison of Table 3 with Table 1, it is evident that the tensor product FE basis implimenation has produced equivalent solution accuracies using half the number of elements as required in the natural coordinate, triangular mesh.

Table 3: Error convergence results for square-shaped region using quadrilateral mesh (k=1).
```====================================================================================
Elements      Qmax   (Qmax)change  (Qmax)slope      ||Q||   change||Q||   Slope(Q)
____________________________________________________________________________________
64    0.298393             0             0   0.274669             0          0
256    0.295597   -0.00279598             0   0.279521    0.00485257          0
1024    0.294912  -0.000684757       2.02969   0.280745    0.00122365    1.98756
4096    0.294742  -0.000170347       2.00712   0.281052   0.000306698     1.9963
16384      0.2947  -4.25345e-005      2.00177   0.281128  7.67318e-005    1.99892
====================================================================================
``` Figure 4A: Temperature solution planform view for "L"-shaped region using quadrilateral mesh (k=1). Figure 4B: Temperature solution angle view for "L"-shaped region using quadrilateral mesh (k=1).

Figures 4A-B detail the steady-state temperature solution for the "L" shaped domain. Table 4 lists the convergence results for the uniform mesh refinement study. Again, a degredation of the convergence rate occurs with an indicated slope near 5/3. Also, the solution again shows equivalent results compared to the natural coordinate implimentation but with half the number of elements.

Table 4: Error convergence results for "L"-shaped region using quadrilateral mesh (k=1).
```===================================================================================
Elements      Qmax  (Qmax)change  (Qmax)slope       ||Q||   change||Q||  Slope(Q)
___________________________________________________________________________________
64    0.141772             0            0   0.0995121             0         0
256    0.147248    0.00547607            0    0.104916    0.00540436         0
1024    0.148688    0.00143958      1.92749    0.106413     0.0014967   1.85234
4096    0.149172    0.00048379       1.5732    0.106844   0.000431161   1.79548
16384    0.149336   0.000164377      1.55737    0.106975   0.000130358   1.72575
===================================================================================
```

Conclusion

Conclusion

In conclusion, a benchmark study has been completed for the case of 2-D steady-state heat conduction through a square-shaped domain with uniform heat generation. As expected, the linear FE basis dominates the converence rate with a slope of 2. When one quarter of the square is removed to create an "L" shaped domain, the concave corner creates a singularity resulting in convergence solution degredation. Now, the data non-smoothness term (r-1) is seen to drive the convergence with an indicated slope near 5/3. Both triangular elements (natural coordinate basis) and quadrilateral elements (tensor product basis) have been implimented for the problem solutions.