**Objectives**

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-ft^{2}-^{o}F and *T _{r} *=
1500

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.

where

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
WS_{e} 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); end 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 plot(1./M(2:end)',eh2,'bo',1./M(2:end)',1000-Q1(2:end),'r-') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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:

DIFF = 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:

DIFFHBC = 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 (T_{b})
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. Q_{5} = T_{b}). 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 (T** _{b}**) and the
constant part of the convection boundary condition at the right end
of the slab (h*T

b = 30000 0 0 0 306.85282

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.

**Conclusion **

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,
“l_{e}“ represents the element length, “Q1”
represents the temperature solution for the
first element, and “e^{h/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, “T_{exact}(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.0000^{o}F at the left
end allows for an exact error calculation (e_{exact}). 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.