**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 load femlib; % 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.

**Grade **

**11.26.04** *Joe, your objective statement is for the steady Pe problem??? The unsteady heaat conduction problem has NO fluid convection term??! Problem statement equations are right, but the text remains specified for Pe? In the discussion and results section you did not follow the directions in conducting the experiments, hence your results defy prediction of the accurate temperature. Conclusions again refer to pe problem, but at least you do comment that non-uniform mesh has some advantages. Make sure to look into the archive when it is opened - the lights will turn on! Score is 4/6 AJB*