Objectives

This lab exercise aims to perform a convergence study for 1-D conduction in a slab of material over a domain with a distributed source. The first problem consists of a uniform slab with a sinusoidal source. The second problem consists of two slabs of uniform dimensions but different thermal conductivities that are fused together. Here a separate sinusoidal source is applied to each slab. The analysis of both problems utilizes both linear and quadratic Finite Element bases. The boundary fluxes for both problems are determined via the GWSh solution. Finally, the convergence of the solution is estimated using the point norm (Tmax) and energy norm principles and compared to theory.

Problem Statement

PART A:

The first part of this lab requests an accuracy/convergence assessment for 1-D conduction across a homogeneous domain with a distributed source. Figure 1 details the given problem. As seen, the conductivity is assumed constant and Dirichlet boundary conditions of a fixed temperature (Tb = 0) are applied at both endpoints of the domain. The source distribution is a sinusoidal function.

Figure 1: 1-D Heat conduction in a slab of material with a distributed source.

The FE solution to Laplace's Heat Equation has already been well defined through Labs 2 and 3 using the FemPSE MATLAB Toolkit. Now, however, a distributed source appears. The MATLAB “linspace” function provides an easy method for the interpolation of the source data to fit within the well-defined template syntax. Over a course mesh, this interpolation is expected to generate considerable data error thus degrading the convergence performance of the FE solution. In the energy norm, the contribution of data interpolation error to the total error estimate is:

where

The MATLAB code given below describes the analysis of this problem:

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #4a                                                                         %
% 1-D, steady-state, conduction in a homogeneous slab with source                        %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;

global X1;	                            %  node coordinates

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

% set up problem parameters
Tb = 0;                                     % fixed temperature applied to both ends
A = 100;                                    % source strength parameter
L = 1;                                      % length of the domain
cond = 0.1;                                 % constant conductivity

for i=1:7

% set up grid and source data (uniform discretization)
nnodes = K*M(i) + 1;
X1 = linspace(0,L,nnodes);
SRC = (A * sin(pi*X1/L))';              % nodal source data (column matrix)

% set up diffusion matrix and column matrix

if K == 1
DIFF = asjac1D(cond,[],[],-1,A211L,[]);
b = asres1D([],[],[],1,A200L,SRC);
end

if K == 2
DIFF = asjac1D(cond,[],[],-1,A211Q,[]);
b = asres1D([],[],[],1,A200Q,SRC);
end

% enforce the left Dirichlet B.C.
DIFFdiri = DIFF;                        % copy DIFF
DIFFdiri(1,:) = zeros(1,nnodes);
DIFFdiri(1,1) = 1;
b(1,1) = Tb;

%  enforce the right Dirichlet B.C.
DIFFdiri(nnodes,:) = zeros(1,nnodes);
DIFFdiri(nnodes,nnodes) = 1;
b(nnodes,1) = Tb;

%  Solve for Q
Q = DIFFdiri \ b;

% compute Energy Norm using [DIFF] (without Dirichlet BC)
Enorm(i,1) = 0.5*Q'*DIFF*Q;

% save max temperature
Qmax(i,1) = max(Q);

end
format long				              % to see all the digits
delta_Q = Qmax(2:end) - Qmax(1:end-1);
eh2 = delta_Q / (2^(2*K) - 1);
Texact = eh2 + Qmax(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

format long g

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

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

PART B:

The second part of this lab analyzes heat conduction through a multi-material composite slab. Figure 2 represents the given problem pictorally. Here, the two slabs are assumed to be fused together. The two materials have two different, but constant, thermal conductivities. Again, Dirichlet boundary conditions (Tb = 0) are applied to the model. Also, the source is distributed as two separate sine functions, one for each material, with differing magnitudes.

Figure 1: 1-D Heat conduction in two slabs of material with equal dimensions but different thermal conductivities and distributed sources.

Now, the temperature solution will be piecewise continuous with the discontinuity occuring at the joining of the two materials (x = L). Because the Taylor Series truncation error estimate depends upon the derivative of the temprature solution, the error estimate is expected to become problematic and depart from theory. The piecewise continuity does not present a problem to the energy norm error estimate, however, which should remain an accurate measure of error convergence.

Finally, the flux at the endpoints of the domain space is determined via Fourier's conduction law:

From the FE solution, the flux becomes a function of the finite temperature difference across the FE of length le:

When evaluated at x = 0 and at x = 2L, this equation should yield a solution for the flux at the boundaries.

The MATLAB code given below describes the analysis of this second problem. It should be noted that this code ensures that an even number of nodes is generated over the domain space which allows one element node to correspond to the boundary between the two materials. Thus the discontinuity present in the conductivities is handled directly by the FE method. Also, the conductivity is modeled as an element property instead of a nodal property. This avoids a potential problem at the nodal interface between the two materials.

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #4b                                                                         %
% 1-D, steady-state, conduction in a homogeneous slab of two disparate materials         %
% each of different thermal conductivity and varying source                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;

global X1;                                  % node coordinates

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

% set up problem parameters
Tb = 0;                                     % fixed temperature applied to both ends
Aleft = 100;                                % source strength in left half-domain
Aright = 10;                                % source strength in right half-domain
L = 1;                                      % half-length of the domain
kleft = 0.1;                                % conductivity in left half-domain
kright = 0.001;                             % conductivity in left half-domain

for i=1:7

% set up grid data (uniform discretization)
nnodes = K*M(i) + 1;
nelements = M(i);
X1 = linspace(0,2*L,nnodes);

% set up nodal source data
SRC = (abs(sin(pi*X1/L)))';             %  uses symmetry of sine for right half-domain
nnhd = (nnodes-1)/2;                    % number of nodes in a half-domain
SRC(1:nnhd,1) = SRC(1:nnhd,1)*Aleft;    % apply Aleft to left half-domain
SRC(nnhd+2:nnodes,1) = SRC(nnhd+2:nnodes,1)*Aright;        % apply Aright

% set up element conductivities
nel2 = nelements/2;                     % number of elements in each half-domain
elcond = ones(nelements,1)*kright;      %create elcond with all values set to kright
elcond(1:nel2,1) = ones(nel2,1)*kleft;  % put in kleft

% set up diffusion matrix and column matrix

if K == 1
DIFF = asjac1D([],elcond,[],-1,A211L,[]);
S = asjac1D([],[],[],1,A200L,[]);
b = S * SRC;
end

if K == 2
DIFF = asjac1D([],elcond,[],-1,A211Q,[]);
S = asjac1D([],[],[],1,A200Q,[]);
b = S * SRC;
end

% enforce the left Dirichlet B.C.
DIFFdiri = DIFF;                           % copy DIFF
DIFFdiri(1,:) = zeros(1,nnodes);
DIFFdiri(1,1) = 1;
b(1,1) = Tb;

%  enforce the right Dirichlet B.C.
DIFFdiri(nnodes,:) = zeros(1,nnodes);
DIFFdiri(nnodes,nnodes) = 1;
b(nnodes,1) = Tb;

%  Solve for Q
Q = DIFFdiri \ b;

% compute Energy Norm using [DIFF] (without Dirichlet BC)
Enorm(i,1) = 0.5*Q'*DIFF*Q;

% save max temperature
Qmax(i,1) = max(Q);

end
format long				                    % to see all the digits
delta_Q = Qmax(2:end) - Qmax(1:end-1);
eh2 = delta_Q / (2^(2*K) - 1);
Texact = eh2 + Qmax(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

format long g

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

% calculate the boundary fluxes
q_left = -kleft * M(i) * (Q(2) - Q(1))
q_right = -kright * M(i) * (Q(end) - Q(end-1))

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

Discussion of Results

PART A:

Figure 1: Steady-state temperature distribution across the given slab length at a 64 element mesh size.

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: Energy norm convergence results with a linear (k=1) Finite Element basis.

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

PART B:

Figure 2: Steady-state temperature distribution across the given composite slab length at a 64 element mesh size.

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: Energy norm convergence results with a linear (k=1) Finite Element basis.

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

Finally, the linear FE basis yields an outward left end flux of 132.7879 and an outward right end flux of 7.0437. For the quadratic FE basis, these fluxes change to 66.4856 and 3.5285 respectively. These values are in units of energy per cross-sectional area.

Conclusion

In conclusion, this a convergence study has been completed for two cases. In the first, a solution has been generated for 1-D conduction in a slab of material over a domain with a distributed source. In the second case, a solution has been generated for the more complex problem with two joined materials of different conductivity and with two different distributed sources. The analysis has been completed for both linear and quadratic Finite Element bases. Error estimations using the point norm (Tmax) and energy norm principles and compared to theory. In addition, for the second case, the boundary fluxes have been determined via the GWSh solution.