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 linear basis convergence study. This study includes verification of the template construction that forms the [DIFF] matrix including modifications for the convection (HBC) and Dirichlet boundary conditions (BC).

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:

The following Matlab code carries out the discussed matrix construction and solves the given conduction problem for a range of nodal elements:

%                                                              %
% Joseph Tipton                                                %
% ES551, LAB #2                                                %
% 1-D, steady-state, conduction in a homogeneous slab          %
%                                                              %

global X1;        % array for node coordinates

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 = 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]
  DIFF = asjac1D(1,[],k,-1,A3011L,[]);

  % 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;
  Q1(i) = Q(1);
format long				% to see all the digits
deltaQ = -(Q1(1:end-1) - Q1(2:end));
eh2 = deltaQ/3;
Texact = eh2 + Q1(2:end);
slope = log10(eh2(1:end-1) ./ eh2(2:end)) / log10(2);
format short		% put Matlab back to default format
answer(:,1) = M';
answer(:,2) = Q1;
answer(:,3) = [0;deltaQ];
answer(:,4) = [0;eh2];
answer(:,5) = [0;Texact];
answer(:,6) = [0;0;slope];
answer(:,7) = 1000-Q1;

save c:/lab2conv answer -ascii -double -tabs



In order to better clarify the construction of the [DIFF] matrix, an example is given for the case with four nodal elements. The initial [DIFF] matrix construction begins with the Jacobian assembly routine, asjac1D. This yields the proper hypermatrix with distributed conductivity:


    45   -45     0     0     0
   -45   100   -55     0     0
     0   -55   120   -65     0
     0     0   -65   140   -75
     0     0     0   -75    75

Next, the appropriate convection boundary condition is added via the [HBC] matrix. Since this boundary condition only affects the left end of the given slab, we expect only the first element of the [DIFF] matrix to change. The [DIFFHBC] matrix confirms this:


    65   -45     0     0     0
   -45   100   -55     0     0
     0   -55   120   -65     0
     0     0   -65   140   -75
     0     0     0   -75    75

The given problem specifies a Dirichlet-type boundary condition on the right end of the slab that results in a fixed temperature (Tb) at x=L. Thus, a further modification to the new [DIFFHBC] matrix involves forcing the final linear equations to represent the known temperature (i.e. Q5 = Tb). The new matrix designation is [Dhdiri] which yields:

DHdiri =

    65   -45     0     0     0
   -45   100   -55     0     0
     0   -55   120   -65     0
     0     0   -65   140   -75
     0     0     0     0     1

Finally, the {b} matrix contains the fixed temperature at the right end of the slab (Tb) and the constant part of the convection boundary condition at the right end of the slab (h*Tr):

b =


The solution to the given conduction problem then reduces to a simple algebraic solution as seen in the code via “Q = DHdiri \ b”.

Finally, it is desirable to incorporate an error estimate due to quantization in this convergence study. Assuming smooth data, which holds true for this study, a Taylor series truncation error estimate combined with the quantization of error from a uniform mesh refinement yields the following relation:

This study uses only linear interpolation between the element solutions and thus k=1. A validation of this error estimate can be made via a slope comparison. The log of the slope of the error versus the mesh refinement (here a factor of 2) yields:

Thus, a slope value near 2 between all meshes confirms the validity of the Taylor series truncation error estimate.


The referenced Matlab program, when run, yields the temperature and error results for mesh sizes of 1 to 32 by factors of 2. The results are seen below in Table 1:

Table 1: Convergence results for 1-D slab conduction with varying conductivity and no sources.

Here, “le“ represents the element length, “Q1” represents the temperature solution for the first element, and “eh/2” represents the Taylor series truncation error estimate. Addition of the leftmost element temperature with the error estimate yields an estimation of the exact temperature, “Texact(est.)”, at x = 0.

The slope calculation allows for validation of the Taylor series truncation error estimate. As seen, the data exhibits a slope very close to 2 with the smaller mesh size beginning to depart. Knowledge of the true analytical solution of 1000.0000oF at the left end allows for an exact error calculation (eexact). In Figure 2 gives a plot of both the Taylor series error estimate and the exact errors versus mesh size. This figure graphically details the slope behavior as well as the departure from the truncation error theory at the lower mesh size.

Figure 2: Actual and estimated truncation results versus element length showing proper error convergence and departure from truncation error assumptions at lower mesh sizes.

Finally, Figure 3 shows the temperature distribution across the length of the conducting slab for the coursest and finest meshes. As seen, the coursest mesh of 1 node results in two nodal temperatures that yields a straight line. The finest mesh of 32 nodes results in 33 nodal temperatures and a decidedly nonlinear temperature profile.

Figure 3: Temperature distribution across the given slab length for the course (M=1) and fine (M=32) mesh sizes.