Objectives

This lab exercise performs a uniform mesh refinement convergence study in the energy norm for the steady "Peclet problem". Both linear and quadratic finite element bases are used over a range of Peclet numbers (1, 10, 100, and 1000). Finally, a non-uniform adaptive mesh following a geometric progression is used in an attempt to achieve improved solution accuracy for Peclet numbers of 100 and 1000.

Problem Statement

This lab involves a convergence study for the classic "Peclet problem". The Peclet problem is a 1-D, source-free, steady-state, heat convection/conduction problem. The non-dimensionalized form is

where represents the non-dimensional temperature,

and "Pe" repsresents the non-dimensional Peclet number which is a function of the Prandtl and Reynolds numbers,

With Dirichlet boundary conditions

and

the analytical solution for this problem is well-known

and results in a sharp temperature "spike" near the right boundary that becomes more pronounced at higher Peclet values as shown in Figure 1 below. At this scale, a Peclet number of 1000 results in a curve almost undistinguishable from the right axis.

Figure 1: Analytical solutions for the Peclet Problem.

# PART A

## Uniform Mesh Refinement

It is first desired to perform a uniform mesh refinement convergence study for Peclet numbers 1, 10, 100, and 1000. The MATLAB template code below impliments the Dirichlet boundary conditions and performs a convergence study using either a linear or quadratic Finite Element basis.

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #11                                                          %
% 1-D, steady-state, convection/diffusion in a homogeneous slab           %
% (The Peclet Problem)                                                    %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;

global X1;                              % array for node coordinates

Pe = input('Enter the Peclet Number: \n\n');
K = input('Enter the basis degree (1 or 2): \n\n');
M = [80 160 320 640];               % number of elements for each run
PEI = 1/Pe;

for i=1:size(M,2)

% uniform discretization,
nnodes = K*M(i) + 1;
X1 = linspace(0,1,nnodes);
nodes(i,1) = nnodes;

% load the FE matrix library

% assemble Matrix for WSe term [DIFF]
if K == 1
DIFF = asjac1D(PEI,[],[],-1,A211L,[]);
DIFF = DIFF + asjac1D([],[],[],0,A201L,[]);
end

if K == 2
DIFF = asjac1D(PEI,[],[],-1,A211Q,[]);
DIFF = DIFF + asjac1D([],[],[],0,A201Q,[]);
end

% set up data matrix (b)
b = zeros(nnodes,1);                  % an nnodes by 1 matrix of zeros

% modify b for the Dirichlet boundary conditions
b(1,1) = 0;
b(nnodes,1) = 1;

% modify DIFFHBC for Dirichlet BC direct solve
DHdiri = DIFF;                        % copy DIFF
DHdiri(1,:) = zeros(1,nnodes);        % zero out the first row
DHdiri(1,1) = 1;                      % i.e., 1*QM = 0
DHdiri(nnodes,:) = zeros(1,nnodes);   % zero out the last row
DHdiri(nnodes,nnodes) = 1;            % i.e., 1*QM = 1

% solve the linear system
Q =  DHdiri \ b;

% compute Energy Norm using [DIFF] (without Dirichlet BC)
enorm(i,1) = 0.5*Q'*DIFF*Q;

subplot(size(M,2),1,i)
plot(X1,Q,'bs-')
axis([0 1 -1 1])

end

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  Nodes   ||Q||      change||Q||     Slope(Q)\n\n')

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

# PART B

## Geometric Mesh Refinement

Finally, it is desired to impliment a non-uniform mesh study for Peclet numbers 100 and 1000. The MATLAB template code below creates a 20 element mesh with decreasing element length according to a geometric progression ratio. The solution utilizes only a linear FE basis.

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #11                                                          %
% 1-D, steady-state, convection/diffusion in a homogeneous slab           %
% (The Peclet Problem)                                                    %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;

global X1;                              % array for node coordinates

Pe = input('Enter the Peclet Number: \n\n');
PEI = 1/Pe;

M = 20;                                 % Mesh size
nnodes = M + 1;

Pr = input('Enter the Geometric Progression Rate: \n\n');

% Create the geometric progression across the mesh
Seed(1) = 1;
for i = 2:M
Seed(i) = Seed(i-1) * Pr;
end

% Convert the Seed array to an element length array
Elength = Seed .* (1/sum(Seed));

% Convert the element length array to the X-axis
X1(1) = 0;
for i = 1:size(Elength,2)
X1(i+1) = X1(i) + Elength(i);
end

% load the FE matrix library

% assemble Matrix for WSe term [DIFF]
DIFF = asjac1D(PEI,[],[],-1,A211L,[]);
DIFF = DIFF + asjac1D([],[],[],0,A201L,[]);

% set up data matrix (b)
b = zeros(nnodes,1);                  % an nnodes by 1 matrix of zeros

% modify b for the Dirichlet boundary conditions
b(1,1) = 0;
b(nnodes,1) = 1;

% modify DIFFHBC for Dirichlet BC direct solve
DHdiri = DIFF;                        % copy DIFF
DHdiri(1,:) = zeros(1,nnodes);        % zero out the first row
DHdiri(1,1) = 1;                      % i.e., 1*QM = 0
DHdiri(nnodes,:) = zeros(1,nnodes);   % zero out the last row
DHdiri(nnodes,nnodes) = 1;            % i.e., 1*QM = 1

% solve the linear system
Q =  DHdiri \ b;

% compute Energy Norm using [DIFF] (without Dirichlet BC)
enorm = 0.5*Q'*DIFF*Q

plot(X1,Q,'bs-')

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

Discussion of Results

# PART A

## Uniform Mesh Refinement

Figures 1A-B and Table 1 detail the convergence results for Pe = 1 using both a linear and quadratic FE basis. The results as seen in Table 1 indicate the baseline 20 element mesh is sufficient to model the temperature solution. The slope values are not trusted due to probable computer round-off errors.

Figure 1A: Temperature solution for the Peclet Problem (Pe = 1) using a linear FE basis.

Figure 1B: Temperature solution for the Peclet Problem (Pe = 1) using a quadratic FE basis.

Table 1: Error convergence results for the Peclet Problem (Pe = 1).
```=======================================================
Elements  FE Basis     ||Q||    change||Q||  Slope(Q)
_______________________________________________________
20       1      0.790892              0         0
40       1      0.790964   7.19411e-005         0
60       1      0.790978   1.33207e-005   2.43315
80       1      0.790982   4.66210e-006   1.51461
_______________________________________________________
20       2      0.790988              0         0
40       2      0.790988  -3.74684e-009         0
60       2      0.790988  -2.00729e-010   4.22235
80       2      0.790988  -3.36140e-011   2.57812
=======================================================
```

Figures 2A-B and Table 2 detail the convergence results for Pe = 10 using both a linear and quadratic FE basis. The results as seen in Table 2 indicate the baseline 20 element mesh is sufficient to model the temperature solution. Again, the slope values are not trusted due to probable computer round-off errors.

Figure 2A: Temperature solution for the Peclet Problem (Pe = 10) using a linear FE basis.

Figure 2B: Temperature solution for the Peclet Problem (Pe = 10) using a quadratic FE basis.

Table 2: Error convergence results for the Peclet Problem (Pe = 10).
```=======================================================
Elements  FE Basis     ||Q||    change||Q||  Slope(Q)
_______________________________________________________
20       1      0.500018              0         0
40       1      0.500022   3.25677e-006         0
60       1      0.500022   6.41145e-007   2.34472
80       1      0.500022   2.27255e-007   1.49634
_______________________________________________________
20       2      0.500023              0         0
40       2      0.500023  -1.87723e-008         0
60       2      0.500023  -9.92590e-010   4.24127
80       2      0.500023  -1.66619e-010   2.57465
=======================================================
```

Figures 3A-B and Table 3 detail the convergence results for Pe = 100 using both a linear and quadratic FE basis. As one might expect, the pronounced wall effect as seen in the analytical solution introduces problems in the numerical solutions. For a linear FE basis, Figure 3A indicates that the oscillating temperature solutions persist until an 80 element is used. For the quadratic FE basis, Figure 3B indicates oscillating temperature solutions until the use of a 60 element mesh. It is interesting to note that the energy norm does not appear to be affected by this oscillatory phenomenon. As seen in Table 3, all meshes yield the same energy norm value of 0.5.

Figure 3A: Temperature solution for the Peclet Problem (Pe = 100) using a linear FE basis.

Figure 3B: Temperature solution for the Peclet Problem (Pe = 100) using a quadratic FE basis.

Table 3: Error convergence results for the Peclet Problem (Pe = 100).
```=====================================================
Elements  FE Basis  ||Q||    change||Q||   Slope(Q)
_____________________________________________________
40       1        0.5              0          0
80       1        0.5   1.11022e-016          0
160       1        0.5   1.11022e-016          0
320       1        0.5  -3.88578e-016   -1.80735
_____________________________________________________
40       2        0.5              0          0
80       2        0.5   1.11022e-016          0
160       2        0.5  -1.16573e-015   -3.39232
320       2        0.5   5.38458e-015    -2.2076
=====================================================
```

Finally, Figures 4A-B and Table 4 detail the convergence results for Pe = 1000 using both a linear and quadratic FE basis. Now, the appropriate monotone temperature solution requires a 640 element mesh for a linear FE basis while the quadratic FE basis requires a 320 element mesh. Again, Table 4 illustrates the energy norm's failure to indicate the oscillatory phenomenon. Instead, all meshes yield the same energy norm of 0.5. Since this energy norm matches that of the Pe = 100 case, it would seem that a limiting situation has been reached for large values of the Peclet number (Pe => 100).

Figure 4A: Temperature solution for the Peclet Problem (Pe = 1000) using a linear FE basis.

Figure 4B: Temperature solution for the Peclet Problem (Pe = 1000) using a quadratic FE basis.

Table 4: Error convergence results for the Peclet Problem (Pe = 1000).
```=====================================================
Elements  FE Basis  ||Q||    change||Q||   Slope(Q)
_____________________________________________________
80       1        0.5              0          0
160       1        0.5  -3.05267e-012          0
320       1        0.5  -5.55112e-017    15.7469
640       1        0.5  -5.55112e-017          0
_____________________________________________________
80       2        0.5              0          0
160       2        0.5              0          0
320       2        0.5   5.55112e-017       -Inf
640       2        0.5  -5.55112e-017          0
=====================================================
```

# PART B

## Geometric Mesh Refinement

Finally, a non-uniform adaptive mesh is implimented in an attempt to mitigate the oscillatory solution phenomenon for a lower element mesh. Figure 5 gives the "correct" temperature solution for the case of Pe = 100 using only 20 elements with a geometric progession ratio of 0.9. For the more severe wall-effect case of Pe = 1000, again a "correct" temperature solution is obtained with a 20 element mesh using a geometric progression ratio of 0.65. Both solutions used a linear FE basis.

Figure 5: Temperature solution for the Peclet Problem (Pe = 100) using a linear FE basis with adaptive, non-uniform meshing.

Figure 6: Temperature solution for the Peclet Problem (Pe = 1000) using a linear FE basis with adaptive, non-uniform meshing.

Conclusion

In conclusion, a uniform mesh refinement convergence study in the energy norm has been completed for the steady "Peclet problem". Results indicate an oscillatory instability in the temperature solution at high Peclet numbers (Pe => 100). This oscillatory phenomenon is seen only via visual inspection since the energy norm does not appear to be affected. In these situations, the correct monotone solutions seem to require meshes on the same order as the Peclet number. Finally, it was seen that the implimentation of a non-uniform, adaptive mesh that follows a decreasing geometric progression provides superior solution accuracy with only a 20 element, computer time-friendly, mesh.