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

Figure 1: Blah, blah.

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #12                                                          %
% 1-D, unsteady conduction in a homogeneous, fin                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% INITIALIZE NON-UNIFORM MESH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global X1;                  % array for node coordinates
Pr = 1.3;                   % geometric progression ratio

% Create a geometric progression across the mesh
Seed(1) = 1;
for i = 2:8                 % where 8 is the initial mesh
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) = 1;
for i = 1:size(Elength,2)
X1(i+1) = X1(i) + Elength(i);
end

nnodes = size(X1,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```
```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #12                                                          %
% 1-D, unsteady conduction in a homogeneous, fin                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% REFINE NON-UNIFORM MESH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Elength_refine = cat(2,Elength./2,Elength./2);
Elength_refine = reshape(Elength_refine,[],2);
Elength_refine = rot90(Elength_refine);
Elength_refine = reshape(Elength_refine,1,[]);

Elength = Elength_refine;

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

nnodes = size(X1,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```
```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #12                                                          %
% 1-D, unsteady conduction in a homogeneous, fin                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% INPUT PROBLEM PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = 10;	                    % conductivity, Btu/(hr*ft*degF)
rho = 1;                    % density, ft^3/lbm
Cp = 1;                     % specific heat, Btu/(lbm*degF)
h = 20;			    % convection coefficient, Btu/(hr*ft^2*degF)
Tr = 1500;		    % convection reference temperature (F)
Tb = 306.8528;		    % Dirichlet BC (F)
delta_t = input('Enter the desired time step in hours: \n\n');
end_t = 0.001;              % final time in hours
theta = 0.5;

% physical constants (calculated)
last_timestep = end_t / delta_t;
L = X1(nnodes) - X1(1);             % length of the domain
Nu = h*L/k;                         % Nusselt number
alpha = rho*Cp/k;                   % thermal diffusivity
r_inner = X1(1);                    % inner radius
T_initial = Tb * ones(nnodes,1);    % initial condition

%
% NON-DIMENSIONAL CONVERSIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Tmax = Tr;                          % max reference dimensional temp
Tmin = Tb;                          % min reference dimensional temp
pot_Tr = (Tr-Tmin)/(Tmax-Tmin);     % convection reference temperature in
% potential form
r_inner_nd = r_inner/L;             % non-dimesnional inner radius
t_ref = L^2 * alpha;                % non-dimensional time
delta_t_nd = delta_t/t_ref;         % non-dimenstional time step
% transform initial condition to potential form
pot_T_initial = (T_initial - Tmin) ./ (Tmax-Tmin);

%
% CREATE GWS+TS MASTER MATRICES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load the FE matrix library

% The required matrices are constant, and need be computed only
% once, outside the time loop
MASS = asjac1D([],[],[X1],1,A3000L,[]);
DIFF = asjac1D([],[],[X1],-1,A3011L,[]);

HBC = zeros(nnodes,nnodes);
HBC(1,1) = r_inner_nd * Nu;

b = zeros(nnodes,1);
b(1) = r_inner_nd * Nu * pot_Tr;

% initialize Qn with the non-dimensional temperature profile
Qn = pot_T_initial;

% make a matrix to store the solutions at each time step
Q_history = zeros(nnodes,last_timestep+1);
Q_history(:,1) = Qn;     % store the initial condition

%
% TIME LOOP SOLUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for step = 1:last_timestep
SYS = MASS + theta*delta_t_nd*(DIFF+HBC);
RESn = (DIFF + HBC)*Qn - b;
RHS = -delta_t_nd * RESn;

% enforce the dirichlet bc on the last node by making the coefficient
% in the nnodes,nnodes position very large
SYS(nnodes,nnodes) = SYS(nnodes,nnodes) * 1E15;

% solve for deltaQ
deltaQ = SYS\RHS;

% update the solution, but don't use deltaQ on the diri node
deltaQ(nnodes) = 0;
Qn = Qn + deltaQ;

% store the solution for this time step
Q_history(:,step+1) = Qn;
end

%
% FORM & DISPLAY RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert the non-dimensional Q_history to dimensional form
T_history = Q_history*(Tmax-Tmin) + Tmin;

% plot the final temperature solution
figure(1)
plot(X1,T_history(:,last_timestep+1))
hold on

% plot the temperature history of the leftmost node
figure(2)
plot(0:delta_t:end_t,T_history(1,:))
hold on

% calculate the energy norm
enorm = 0.5*Qn'*DIFF*Qn

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

Discussion of Results

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 2A: Temperature solution for the Peclet Problem (Pe = 1) using a linear FE basis.

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

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

Figure 3B: 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  Time Step        Tmax  change(Tmax)     ||Q||  change||Q||   slope
(hours)         (F)           (F)
______________________________________________________________________________
8	   0.00001    524.735843                0.091791
16	   0.00001    525.842399     1.106556   0.092491  0.000699941
32	   0.000005   526.121550     0.279151   0.092792  0.000301851  1.2134
64	   0.000001   526.191497     0.069947   0.092937  0.000144459  1.0632
==============================================================================
```

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

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

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.