**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 GWS^{h} solution. Finally, the convergence of the solution is estimated using the point norm (T_{max}) 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 (T_{b} = 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) load femlib; % load library matrix data % 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 answer(:,1) = M'; answer(:,2) = Qmax; answer(:,3) = [0;delta_Q]; answer(:,4) = [0;eh2]; answer(:,5) = [0;Texact]; answer(:,6) = [0;0;slope_Taylor]; answer(:,7) = Enorm; answer(:,8) = [0;error_est]; answer(:,9) = [0;Enorm_exact_est]; answer(:,10) = [0;0;slope_Enorm]; format long g answer 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 (T_{b} = 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 l_{e}:

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 % load library matrix data load femlib; % 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 answer(:,1) = M'; answer(:,2) = Qmax; answer(:,3) = [0;delta_Q]; answer(:,4) = [0;eh2]; answer(:,5) = [0;Texact]; answer(:,6) = [0;0;slope_Taylor]; answer(:,7) = Enorm; answer(:,8) = [0;error_est]; answer(:,9) = [0;Enorm_exact_est]; answer(:,10) = [0;0;slope_Enorm]; format long g answer 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: T_{max} convergence results with a linear (*k*=1) Finite Element basis.

Table 2: T_{max} 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: T_{max} convergence results with a linear (*k*=1) Finite Element basis.

Table 2: T_{max} 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 (T_{max}) and energy norm principles and compared to theory. In addition, for the second case, the boundary fluxes have been determined via the GWS^{h} solution.

**Grade **

**Grade **

**9.27.04** *Joe, objectives and problem statement are fine for lab 4. However, you have some serious bugs in your computation of energy norm for this problem. We discussed this on the phone today. Your convergence data is not in agreement with what is anticipated for either case a or b. However, your data resentations are good and concise. Conclusions cannot make sense out of bad data. You need to check your Enorm computation. Score is 5/6 AJB*