In general, this lab aims to introduce the student to basic Finite Element (FE) model design, construction, and solution. Specifically, this lab introduces the FEmPSE Matlab toolbox which provides a general framework for the construction and solution of FE problems. These objectives are accomplished through a simple 1-D conduction heat transfer convergence study.

In Lab #2, this same problem was solved in a linear FE basis convergence study. For this lab, the convergence study is expanded to a quadratic and cubic FE basis. In addition, this lab introduces the Energy Norm concept to provide a more robust error estimate as compared to the point norm error estimate based upon Taylor series truncation.

Problem Statement

The problem statement as given by the course instructor states the following:

“Determine the steady-state temperature distribution in a slab of thickness L = 1ft subject to a convective thermal load at the left end with h = 20 Btu / hr-ft2-oF and Tr = 1500oF to an accuracy of 0.1oF. The thermal conductivity of the slab material is a linear function of x, such that k(x = 0) = 10 Btu / hr-ft-oF while k(x = L) = 20 Btu / hr-ft-oF. There is no heat source present. At x = L the slab temperature is held fixed at Tb = 306.85282 oF. The analytical solution at x = 0 is 1000.0000 oF.”

Analysis of this problem further assumes a slab of homogeneous material with 1-D conduction along the thickness dimension. In this case, and with no sources present, the conservation of energy principle along with Fourier's Law of Conduction yields Laplace's Heat Equation. The convective thermal load at the left end of the slab yields a Robin BC while the fixed temperature at the right end of the slab yields a Dirichlet BC. The equations below represent the mathematical formulation while Figure 1 offers a pictorial representation of the problem at hand.


and with boundary conditions

Figure 1: 1-D Heat conduction in a slab of material with length, L.

Discussion of Results

The Galerkin Weak Statement, when discretized according to the finite elements, assembles an algebraic system of equations:

The three main portions of this linear system equations, [DIFF]{Q}, [HBC]{Q}, and {b}, are then translated into the more generic Matlab FEMLib template:

where the template syntax amounts to:

This arrangement assumes, however, that the conductivity {COND} is taken as an average value over the element span. It is preferable, instead, to allow the conductivity to vary across the elements and thus across the slab length as specified in the problem statement. This can be represented in the following manner:

In doing so, the extra nodal matrix incorporates itself into the [DIFF] matrix and thus the number of FE bases in the integrand that form the WSe term increases from 2 to 3. The extra zero signifies that this extra FE base is not differentiated:

Thus, the conductivity matrix {COND} no longer represents average data. Instead, it becomes a distributed data matrix. The generic Matlab FEMlib template now changes accordingly with the conductivity in the distributed data matrix and the new “A3011L” hypermatrix. For a quadratic and cubic FE basis, the hypermatrix as called from the FEMLib toolbox becomes “A3011Q” and “A3011C” respectively.

The following Matlab code carries out the discussed matrix construction and solves the given conduction problem for a range of nodal elements. It has the capability to choose a linear, quadratic, or cubic FE basis depending upon user input.

%                                                                         %
% Joseph Tipton                                                           %
% ES551, LAB #3                                                           %
% 1-D, steady-state, conduction in a homogeneous slab                     %
% (for 1-3 degree Lagrange basis | using energy norm error estimation)    %
%                                                                         %

clear all;

global X1;        			% array for node coordinates

K = input('Enter the basis degree (for Lagrange interpretation): \n\n');
M = [1 2 4 8 16 32]; 			% number of elements for each run
Q1 = zeros(6,1);
for i=1:6

  % constant data
  h = 20;           			% convection coefficient
  Tr = 1500;        			% reference temperature
  Tb = 306.85282;   			% fixed temperature

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

  % distributed data at the nodes
  k = X1' * 10 + 10;

  % load the FE matrix library
  load femlib;

  % assemble Matrix for WSe term [DIFF]
  if K == 1
      DIFF = asjac1D(1,[],k,-1,A3011L,[]);
  if K == 2
      DIFF = asjac1D(1,[],k,-1,A3011Q,[]);
  if K == 3
      DIFF = asjac1D(1,[],k,-1,A3011C,[]);

  % complete Matrix for [HBC] on element 1
  DIFFHBC = DIFF;   			% copy DIFF
  DIFFHBC(1,1) = DIFFHBC(1,1) + h; 	% modify the 1,1 entry

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

  % modify b for convective BC at node 1
  b(1,1) = h * Tr;

  % modify b for the Dirichlet boundary condition
  b(nnodes,1) = Tb;

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

  % solve the linear system
  Q =  DHdiri \ b;
  % compute Energy Norm using [DIFF + HBC] (without Dirichlet BC)
  Enorm(i,1) = 0.5*Q'*DIFFHBC*Q;
  % save temperatures at the left end of the slab
  Q1(i) = Q(1);
format long				% to see all the digits
delta_Q = Q1(2:end) - Q1(1:end-1);
eh2 = delta_Q / (2^(2*K) - 1);
Texact = eh2 + Q1(2:end);
slope_Taylor = log10(eh2(1:end-1) ./ eh2(2:end)) / log10(2);
error_est =  (Enorm(2:end) - Enorm(1:end-1)) ./  (2^(2*K) - 1);
Enorm_exact_est = Enorm(2:end) + error_est;
slope_Enorm = log10(error_est(1:end-1) ./ error_est(2:end)) / log10(2);
format short				% put Matlab back to default format
answer(:,1) = M';
answer(:,2) = Q1;
answer(:,3) = [0;delta_Q];
answer(:,4) = [0;eh2];
answer(:,5) = [0;Texact];
answer(:,6) = [0;0;slope_Taylor];
answer(:,7) = 1000-Q1;
answer(:,8) = Enorm;
answer(:,9) = [0;error_est];
answer(:,10) = [0;Enorm_exact_est];
answer(:,11) = [0;0;slope_Enorm];

format long g


As seen, the above MATLAB program utilizes two methods to determine an estimate of the solution error. Assuming smooth data, which holds true for this study, a Taylor series truncation error estimate follows the relationship:

By taking the LOG10 of both sides of the equation, a convienient relationship arises for the error convergence. If plotted as log(le) versus log(eh), the asymptotic convergence error relationship yields a straight line of slope 2k where k is the order FE basis:

Now, for a uniform mesh refinement of factor two, a precise error quantification can be expressed in terms of a nodal temperature change after refinement (without knowledge of the exact solution):

The validity of this relationship depends upon actual data adherence to the theoretical convergence curve. Thus a simple slope computation verifies the founding error relationship. The calculated slope should match the theoretical slope of 2k for verification of the Taylor series truncation error estimate:

The energy norm concept allows for a final, more robust, error convergence model. While the Taylor series truncation error estimate uses only one nodal point for comparison, the energy norm uses the entire nodal solution:

In addition, the energy norm does not require a smooth data set as does the previous method. The energy norm follows the same previously discussed asymptotic error relationship between temperature solution and error. Thus the difference between temperature energy norms across a uniform mesh refinement yields an estimate of the error energy norm according to:

Finally, a similar slope calculation allows for verification of the asymptotic error estimate relationship that permitted calculation of the above energy norm error estimate:


Figure 2 displays the steady-state temperature distribution along the slab with the 32 element mesh. At this mesh size, the solution is virtually identical for the linear, quadratic, and cubic FE bases. The temperature profile is decidedly non-linear and properly reaches a maximum temperature at x=0 of 1000.0000oF. Figure 3 shows a comparison of the 1 element mesh temperature distribution across the slab for all three FE bases. As expected, the cubic FE basis provides the more robust approximation of the exact temperature distribution even at the coursest mesh.

Figure 2: Steady-state temperature distribution across the given slab length at a 32 element mesh size for all three FE bases (k=1,2,3).

Figure 2: Steady-state temperature distribution across the given slab length at a 1 element mesh size for all three FE bases (k=1,2,3).

The tables below list the temperature and point norm error results for increasing meshes of 1, 2, 4, 8, 16, and 32 elements as calculated via the referenced MATLAB program. Three cases were performed using linear, quadratic, and cubic FE bases:

Table 1: Tmax convergence results with a linear (k=1) Finite Element basis.

Table 2: Tmax convergence results with a quadratic (k=2) Finite Element basis.

Table 3: Tmax convergence results with a cubic (k=3) Finite Element basis.

Here, “le“ represents the element length, “Q1” represents the temperature solution for the first element (which is also the maximum temperature, Tmax), and “eh/2” represents the point norm error estimate at Tmax which is a Taylor series truncation error estimate. Addition of Tmax with the error estimate yields an estimation of the exact temperature, “Texact(est.)”, at x = 0. As seen in Figure 3 and confirmed in the tables, the greater order FE basis provides a better Tmax solution for a given mesh. Here, the linear FE basis solution reaches the known analytical solution to four decimal places (1000.0000oF) with a 32 element mesh. This compares to the quadratic FE basis solution at an 8 element mesh and the cubic FE basis solution at only a 4 element mesh.

Figure 4 gives a plot of the Tmax point norm error estimation for the three FE bases. As discussed, theory expects linear convergence trends of slope 2, 4, and 6 for the linear, quadratic, and cubic FE bases respectively. In reality, this does not hold for all points. As detailed in Tables 1-3, the slope at the course end of the mesh begins to deviate from the expected slope. This is attributed to a depedence of the constant, C, on mesh size for the courser meshes. In addition, Figure 4 and Table 3 indicate an extreme deviation from the expected linear behavior for the cubic FE basis. This is attributed to computer round-off errors.

Figure 4: Tmax point norm error estimation for the three FE bases (k=1,2,3).

Finally, Tables 4-6 relate the convergence of the error estimate in terms of the more robust energy norm. As expected, the energy norm error estimate follows the same patterns seen in the previous discussion. Figure 5 provides a Log-Log plot of the energy norm error estimate versus the element length. As expected, the same patterns are seen as in Figure 4. The previous paragraph also explains the observed deviations from the theoretical slopes.

Table 4: Energy norm convergence results with a linear (k=1) Finite Element basis.

Table 5: Energy norm convergence results with a quadratic (k=2) Finite Element basis.

Table 6: Energy norm convergence results with a cubic (k=3) Finite Element basis.

Figure 5: Energy norm error estimation for the three FE bases (k=1,2,3).


9.17.04 Joe, labs 2 and 3 both have the content required, hence the scores are 2/2 and 6/6. However, you are not following directions on content by subject. Specifically, ALL theory materials belong in the problem statement, ALL computed data belong in the discussion and results section, and the conclusions section should be very short in answering the objectives. You are pushing theory into discussion and data in conclusions, which is not appropriate. Please heed these comments for the remaining labs. AJB