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
  load femlib;

  % 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);

ANSWER(:,1) = M';
ANSWER(:,2) = nodes;
ANSWER(:,3) = enorm;
ANSWER(:,4) = [0;delta_enorm];
ANSWER(:,5) = [0;0;slope_enorm];

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
load femlib;

% 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.

Grade

11.18.04 Dr. Baker, I briefly checked my Lab 11 report and, indeed, I am calculating the energy norm as: Enorm = 0.5 * Q' * DIFF * Q When you grade my lab report, you will see that my results seem to indicate a limiting value of 0.5 as the Peclet number increases from 1 to 1000. I'm currently at a loss to explain how this is happening. - Joseph

Troops, a problem assignment for section CD1 asked you to verify that the E norm for the Pe problem approached 0.25 as Pe approached a large number. Your Lab 11 reports all indicate that the E norm is approaching 1/2. Anyone have an idea as to the cause of the discrepancy??

Joe, your lab 11 report is very complete. Your lack of firm convergence data for the two small Pe cases for k = 1,2 is because you did not follow the "rule" of doubling the mesh! You went 20, 40, 60, 80 when you should have gone 20, 40, 80, 160!! Look at Scott Lusted's report to see these firm data. Non uniform mesh solution graphs look right. i do not understand your comments about the E norms though. Dispersion error adds huge amounts of false energy to the E norm - you cannot do convergence study until the solution is approximately monotone. I note in all reports that the E norm is approaching 0.5 as Pe gets large. that means something is wrogn, either with the coding or with the problem assignment that yielded 0.25. Will ask the class to figure this out. lab score is 6/6 AJB