Objectives

This lab exercise aims to perform a convergence study for 1-D conduction in a uniform slab of material over a given domain. For the first problem, the thermal conductivity varies linearly with the temperature. For the second problem, the thermal conductivity varies quadratically with the temperature. Both problems thus introduce non-linearities into the problem and require a new solution scheme via Newton's Method. The analysis of both problems utilizes both linear and quadratic Finite Element bases. The convergences of the solutions are estimated using the point norm (Tmid) and energy norm principles and compared to theory. This is then repeated for a modified Newton's Method solution process that omits the non-linear portions of the Jacobian matrix.

Problem Statement

This lab exercise concerns the conduction of heat through a 1 ft. thick slab with temperature dependent conductivity. The temperature at the left end of the slab (Tleft) is held fixed at 2000oF while the right end (Tright) is maintained at 1000oF. In Part A of this lab, a quadratic function describes the temperature dependent conductivity, k(T) = A + BT + CT2, where A = -1.00, B = 0.002, and C = 1.0e-5. For Part B, a linear function describes the conductivity, k(T) = kright + B(T – Tright), where kright = 1.0 and B = 0.002. In both situations, the units of conductivity remain in the English system (Btu / hr-ft).

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 fixed temperatures at both ends of the slab yield Dirichlet boundary conditions. The equations below represent the mathematical formulation while Figure 1 offers a pictorial representation of the problem at hand.

where:

with two different models of heat conduction as a function of temperature:


with boundary conditions:



Figure 1: 1-D heat conduction in a slab of material with no source, Dirichlet boundary conditions, and a temperature dependent thermal conductivity.

The MATLAB program file for the FE solution to this problem utilizes the FEMLib toolkit and a similar structure as with the other labs. Now, however, the temperature dependence in the conductivity creates a non-linear problem. An iterative process following the Newton Method allows for a solution with the non-linearities. The Newton Method solution format is

where {FQ} signifies the Galerkin Weak Statement (GWSh). The Jacobian is the partial derivative of WSh with respect to the temperature, Q.

Since Dirichlet conditions yield fixed temperatures at both ends of the slab, the respective nodes of the Jacobian matrix are rewritten as very large numbers, signifying an approximately non-existant temperature change at those nodes. With these changes defined, the relationship above allows for an iterative solution of the temperature until the change reaches some predefined threshold.

It should be noted that the non-linearities are contained in the Jacobian via the [A3110_] matrices. By simply “commenting-out” these lines to the code, the effect of the non-linearity may be ascertained in terms of the iteration convergence rate and the mesh refinement error convergence.

PART A


Conductivity Linear Temperature Variation:     k(T) = kright + B * (T - Tright)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #5a                                                                         %
% 1-D, steady-state, conduction in a homogeneous slab with Dirichlet boundary conditions %
% and a thermal conductivity that is a linear function of temperature                    %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;
global X1;
format compact                              % suppress extra line feeds

%  set problem parameters
K = input('Enter the basis degree (for Lagrange interpretation): \n\n');
M = [2 4 8 16 32 64 128];                   % number of elements for each run
kright = 1.0 ;                              % k(T) = kright + beta * (T - Tright)
beta = .002;
Tleft = 2000;                               % left end Dirichlet condition
Tright = 1000;                              % right end Dirichlet condition

%  set max iteration count and convergence criteria
maxit = 20;                                 % max iterations
eps = .001;                                 % convergence tolerance

%  load the element library matrices
load femlib

for i = 1:size(M,2)

    % set up the linearly spaced mesh between 0 and 1
    nnodes = K*M(i) + 1;
    X1 = linspace(0,1,nnodes);
	
    %  initial condition
    Q = (linspace(Tleft,Tright,nnodes))';   
	
    %  iteration loop
    for j = 1:maxit;
        
      % form residual
      if K == 1
          res =       asres1D(kright,[],[],-1,A211L,Q);
          res = res + asres1D(beta,[],Q,-1,A3011L,Q);
          res = res + asres1D(-1*beta*Tright,[],[],-1,A211L,Q);
      end
      
      if K == 2
          res =       asres1D(kright,[],[],-1,A211Q,Q);
          res = res + asres1D(beta,[],Q,-1,A3011Q,Q);
          res = res + asres1D(-1*beta*Tright,[],[],-1,A211Q,Q);
      end
	
      % form jacobian
      if K == 1
          jac =       asjac1D(kright,[],[],-1,A211L,[]);
          jac = jac + asjac1D(beta,[],Q,-1,A3011L,[]);
          jac = jac + asjac1D(beta,[],Q,-1,A3110L,[]);
          jac = jac + asjac1D(-1*beta*Tright,[],[],-1,A211L,[]);
      end
      
      if K == 2
          jac =       asjac1D(kright,[],[],-1,A211Q,[]);
          jac = jac + asjac1D(beta,[],Q,-1,A3011Q,[]);
          jac = jac + asjac1D(beta,[],Q,-1,A3110Q,[]);
          jac = jac + asjac1D(-1*beta*Tright,[],[],-1,A211Q,[]);
      end
	
      % modify the jacobian to enforce Dirichlet conditions at both ends
      jac(1,1) = 1E10;
      jac(nnodes,nnodes) = 1E10;
	
      dQ = jac\(-1*(res));   % linear solve for dQ
      Q = Q + dQ;            % update Q
      Q(1,1) = Tleft;        % undo dQ for Dirichlet node
      Q(nnodes,1) = Tright;  % undo dQ for Dirichlet node
      dQmax = max(abs(dQ));  % convergence criteria
      if dQmax < eps         % test for convergence against tolerance eps
        break                % quit if converged
      end
    end
    
    % save iteration convergence count
    iteration(i,1) = j;
    
    % calculate and display energy norm. 
    % res must be recalculated using the converged solution
      if K == 1
          res =       asres1D(kright,[],[],-1,A211L,Q);
          res = res + asres1D(beta,[],Q,-1,A3011L,Q);
          res = res + asres1D(-1*beta*Tright,[],[],-1,A211L,Q);
      end
      
      if K == 2
          res =       asres1D(kright,[],[],-1,A211Q,Q);
          res = res + asres1D(beta,[],Q,-1,A3011Q,Q);
          res = res + asres1D(-1*beta*Tright,[],[],-1,A211Q,Q);
      end
    Enorm(i,1) = 0.5 * Q' * res;
      
    % save midpoint temperatures
    Qmax(i,1) = Q((0.5*K*M(i)) + 1);
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) = iteration;
answer(:,3) = Qmax;
answer(:,4) = [0;delta_Q];
answer(:,5) = [0;eh2];
answer(:,6) = [0;Texact];
answer(:,7) = [0;0;slope_Taylor];
answer(:,8) = Enorm;
answer(:,9) = [0;error_est];
answer(:,10) = [0;Enorm_exact_est];
answer(:,11) = [0;0;slope_Enorm];

fprintf('Mesh  Iteration  Tmid  Tmid_change  TS_error  Tmid_exact  TS_Slope  Enorm ...
		... Enorm_error  Enorm_exact  Enorm_Slope\n\n')
fprintf('%6.0f  %2.0f  %12.4f  %12.6g  %12.6g  %12.4f  %12.6g  %12.6g  %12.6g  %12.6g ...
		... %12.6g\n',answer');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PART B


Conductivity Quadratic Temperature Variation:     k(T) = A + B*T + C*T2

Here, the thermal conductivity is a quadratic function of temperature.

The MATLAB toolkit template syntax for the {WS}e contribution of the thermal conductivity then follows.

The quadratic term yields a hypermatrix of order [A40110_]. Interpolation of the quadratic term via the element nodal product {Q,Q} avoids the use of this unnecessarily cumbersome matrix.

Aside from this change, the resulting MATLAB program code follows that of Part A.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #5b                                                                         %
% 1-D, steady-state, conduction in a homogeneous slab with Dirichlet boundary conditions %
% and a thermal conductivity that is a quadratic function of temperature                 %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;
global X1;
format compact                              % suppress extra line feeds

%  set problem parameters
K = input('Enter the basis degree (for Lagrange interpretation): \n\n');
M = [2 4 8 16 32 64 128];                   % number of elements for each run
A = -1;                                     % k = A + B*T + C*T^2
B = .002;
C = 1E-5;
Tleft = 2000;                               % left end Dirichlet condition
Tright = 1000;                              % right end Dirichlet condition

%  set max iteration count and convergence criteria
maxit = 20;                                 % max iterations
eps = .001;                                 % convergence tolerance

%  load the element library matrices
load femlib

for i = 1:size(M,2)

    % set up the linearly spaced mesh between 0 and 1
    nnodes = K*M(i) + 1;
    X1 = linspace(0,1,nnodes);
	
    %  initial condition
    Q = (linspace(Tleft,Tright,nnodes))';   
	
    %  iteration loop
    for j = 1:maxit;
	
      % form residual
      if K == 1
          res =       asres1D(A,[],[],-1,A211L,Q);
          res = res + asres1D(B,[],Q,-1,A3011L,Q);
          res = res + asres1D(C,[],Q.^2,-1,A3011L,Q);
      end
      
      if K == 2
          res =       asres1D(A,[],[],-1,A211Q,Q);
          res = res + asres1D(B,[],Q,-1,A3011Q,Q);
          res = res + asres1D(C,[],Q.^2,-1,A3011Q,Q);
      end
	
      % form jacobian
      if K == 1
          jac =       asjac1D(A,[],[],-1,A3011L,[]);
          jac = jac + asjac1D(B,[],Q,-1,A3011L,[]);
          jac = jac + asjac1D(B,[],Q,-1,A3110L,[]);
          jac = jac + asjac1D(C,[],Q.^2,-1,A3011L,[]);
          jac = jac + asjac1D(2*C,[],Q,-1,A3110L,Q);
      end
      
      if K == 2
          jac =       asjac1D(A,[],[],-1,A3011Q,[]);
          jac = jac + asjac1D(B,[],Q,-1,A3011Q,[]);
          jac = jac + asjac1D(B,[],Q,-1,A3110Q,[]);
          jac = jac + asjac1D(C,[],Q.^2,-1,A3011Q,[]);
          jac = jac + asjac1D(2*C,[],Q,-1,A3110Q,Q);
      end
	
      % modify the jacobian to enforce Dirichlet conditions at both ends
      jac(1,1) = 1E15;
      jac(nnodes,nnodes) = 1E15;
	
      dQ = jac\(-1*(res));   % linear solve for dQ
      Q = Q + dQ;            % update Q
      Q(1,1) = Tleft;        % undo dQ for Dirichlet node
      Q(nnodes,1) = Tright;  % undo dQ for Dirichlet node
      dQmax = max(abs(dQ));  % convergence criteria
      if dQmax < eps         % test for convergence against tolerance eps
        break                % quit if converged
      end
    end
    
    % save iteration convergence count
    iteration(i,1) = j;
    
    % calculate and display energy norm. 
    % res must be recalculated using the converged solution
      if K == 1
          res =       asres1D(A,[],[],-1,A211L,Q);
          res = res + asres1D(B,[],Q,-1,A3011L,Q);
          res = res + asres1D(C,[],Q.^2,-1,A3011L,Q);
      end
      
      if K == 2
          res =       asres1D(A,[],[],-1,A211Q,Q);
          res = res + asres1D(B,[],Q,-1,A3011Q,Q);
          res = res + asres1D(C,[],Q.^2,-1,A3011Q,Q);
      end
    Enorm(i,1) = 0.5 * Q' * res;
      
    % save midpoint temperatures
    Qmax(i,1) = Q((0.5*K*M(i)) + 1);
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) = iteration;
answer(:,3) = Qmax;
answer(:,4) = [0;delta_Q];
answer(:,5) = [0;eh2];
answer(:,6) = [0;Texact];
answer(:,7) = [0;0;slope_Taylor];
answer(:,8) = Enorm;
answer(:,9) = [0;error_est];
answer(:,10) = [0;Enorm_exact_est];
answer(:,11) = [0;0;slope_Enorm];

fprintf('Mesh  Iteration  Tmid  Tmid_change  TS_error  Tmid_exact  TS_Slope  Enorm ...
		... Enorm_error  Enorm_exact  Enorm_Slope\n\n')
fprintf('%6.0f  %2.0f  %12.4f  %12.6g  %12.6g  %12.4f  %12.6g  %12.6g  %12.6g  %12.6g ...
		... %12.6g\n',answer');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Discussion of Results

PART A


Conductivity Linear Temperature Variation:     k(T) = kright + B * (T - Tright)


Figure 1: Steady-state temperature distribution across the given slab at a 128 element mesh size with quadratic FE basis. Here, the thermal conductivity varies linearly with the temperature.

Figure 1 above shows the temperature distribution along the slab for the case of a thermal conductivity linear temperature dependence. This solution utilized a quadratic element basis with 128 mesh elements. The figure confirms that the solution meets the Dirichlet BCs at both ends of the slab.

Tables 1-4 list the convergence results for the case of a linear FE basis with and without the non-linear Jacobian contributions. As seen in Tables 1 and 3, the GWSh solved via Newton's Method yields a nodally exact solution at every mesh. The large flucuations in the error convergence slopes are attributed to computer round-off error. Tables 2 and 4 indicate that, with removal of the non-linear portions of the Jacobian matrix, the solutions are no longer exact. Thus, the error convergence slopes once again approach 2, which follows theory, and the number of iterations required to converge upon a solution using Newton's Method doubles.

Table 1: Taylor Series truncation error convergence results for a linear (k=1) Finite Element basis.
===========================================================================================
  Mesh  Iterations          Tmid    Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)            (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           4     1618.0340              0              0        0.0000             0 
     4           4     1618.0340              0              0     1618.0340             0 
     8           4     1618.0340              0              0     1618.0340           NaN 
    16           4     1618.0340   2.27374e-013   7.57912e-014     1618.0340          -Inf 
    32           4     1618.0340  -4.54747e-013  -1.51582e-013     1618.0340            -1 
    64           4     1618.0340  -1.13687e-012  -3.78956e-013     1618.0340      -1.32193 
   128           4     1618.0340   1.81899e-012    6.0633e-013     1618.0340     -0.678072 
===========================================================================================

Table 2: Taylor Series truncation error convergence results for a linear (k=1) Finite Element basis without the non-linear Jacobian contribution.
===========================================================================================
  Mesh  Iterations          Tmid    Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)            (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           6     1618.0340              0              0        0.0000             0 
     4           7     1618.0340  -1.26093e-005  -4.20309e-006     1618.0340             0 
     8           7     1618.0338   -0.000163938  -5.46461e-005     1618.0338       -3.7006 
    16           8     1618.0340     0.00021447   7.14901e-005     1618.0341     -0.387626 
    32           8     1618.0340   3.56257e-006   1.18752e-006     1618.0340       5.91172 
    64           8     1618.0340   9.05348e-007   3.01783e-007     1618.0340       1.97637 
   128           8     1618.0340   2.25515e-007   7.51716e-008     1618.0340       2.00525 
===========================================================================================

Table 3: Energy Norm convergence results for a linear (k=1) Finite Element basis.
===========================================================================
  Mesh  Iterations      Enorm   Enorm_error   Enorm_exact   Enorm_Slope    
  Size    Required                                                         
___________________________________________________________________________
     2           4     1e+006             0             0             0    
     4           4     1e+006  2.05783e-006        1e+006             0    
     8           4     1e+006  4.20473e-006        1e+006      -1.03089    
    16           4     1e+006  3.69359e-006        1e+006       0.18699    
    32           4     1e+006  2.53553e-006        1e+006      0.542737    
    64           4     1e+006  1.51577e-006        1e+006      0.742238    
   128           4     1e+006  8.36366e-007        1e+006       0.85784    
===========================================================================

Table 4: Energy Norm convergence results for a linear (k=1) Finite Element basis without the non-linear Jacobian contribution.
============================================================================
  Mesh  Iterations      Enorm    Enorm_error   Enorm_exact   Enorm_Slope    
  Size    Required                                                          
____________________________________________________________________________
     2           6     1e+006              0             0             0    
     4           7     1e+006    -0.00370151        1e+006             0    
     8           7     1e+006     -0.0359057        1e+006      -3.27803    
    16           8     1e+006      0.0412412        1e+006     -0.199873    
    32           8     1e+006   -0.000788305        1e+006       5.70919    
    64           8     1e+006   -0.000269675        1e+006       1.54753    
   128           8     1e+006  -7.30884e-005        1e+006       1.88351    
============================================================================

Tables 5-8 list the convergence results for the case of a quadratic FE basis with and without the non-linear Jacobian contributions. The results follow the same pattern as seen for the linear FE basis solutions. When the non-linear portions of the Jacobian matrix stay, the code generates an exact nodal solution. Thus the error convergence slopes become useless. With ommission of the non-linear portions of the Jacobian, however, the solutions are no longer exact and the error convergence slopes approach a value of 4 which agrees with theory. Finally, Table 8 adds one more "twist" to the results. Here, even with the degredation of the solution via omission of the non-linear Jacobian portions, the quadratic FE basis yields a robust solution with errors that approach the limit of computer round-off. This results once again in useless convergence rate slope calculations.

Table 5: Taylor Series truncation error convergence results for a quadratic (k=2) Finite Element basis.
===========================================================================================
  Mesh  Iterations          Tmid    Tmid_change     Truncation    Tmid_exact      TS_Slope 
  Size    Required           (F)            (F)     Error Est.    (Est.) (F)               
___________________________________________________________________________________________
     2           4     1618.0340              0              0        0.0000             0 
     4           4     1618.0340  -9.09495e-013   -6.0633e-014     1618.0340             0 
     8           4     1618.0340  -6.36646e-012  -4.24431e-013     1618.0340      -2.80735 
    16           4     1618.0340   2.95586e-011   1.97057e-012     1618.0340      -2.21501 
    32           4     1618.0340   5.77529e-011   3.85019e-012     1618.0340     -0.966317 
    64           4     1618.0340  -7.07132e-011  -4.71421e-012     1618.0340     -0.292086 
   128           4     1618.0340   4.54747e-013   3.03165e-014     1618.0340       7.28077 
===========================================================================================

Table 6: Taylor Series truncation error convergence results for a quadratic (k=2) Finite Element basis without the non-linear Jacobian contribution.
===========================================================================================
  Mesh  Iterations          Tmid    Tmid_change     Truncation    Tmid_exact      TS_Slope 
  Size    Required           (F)            (F)     Error Est.    (Est.) (F)               
___________________________________________________________________________________________
     2           7     1618.0339              0              0        0.0000             0 
     4           8     1618.0340    0.000102911   6.86072e-006     1618.0340             0 
     8           8     1618.0340   7.62007e-006   5.08005e-007     1618.0340       3.75545 
    16           8     1618.0340   9.56038e-007   6.37359e-008     1618.0340       2.99466 
    32           8     1618.0340   7.47518e-008   4.98345e-009     1618.0340       3.67689 
    64           8     1618.0340   4.48495e-009   2.98996e-010     1618.0340       4.05895 
   128           8     1618.0340  -5.00449e-010  -3.33633e-011     1618.0340       3.16379 
===========================================================================================

Table 7: Energy Norm convergence results for a quadratic (k=2) Finite Element basis.
==========================================================================
  Mesh  Iterations      Enorm    Enorm_error   Enorm_exact   Enorm_Slope  
  Size    Required                                                             
__________________________________________________________________________
     2           4     1e+006              0             0             0  
     4           4     1e+006   1.00429e-006        1e+006             0  
     8           4     1e+006   7.07347e-007        1e+006      0.505688  
    16           4     1e+006   4.41253e-007        1e+006      0.680813  
    32           4     1e+006   1.76858e-007        1e+006       1.31901  
    64           4     1e+006   2.74857e-007        1e+006     -0.636086  
   128           4     1e+006  -2.80662e-007        1e+006    -0.0301538  
==========================================================================

Table 8: Energy Norm convergence results for a quadratic (k=2) Finite Element basis without the non-linear Jacobian contribution.
===========================================================================
  Mesh  Iterations      Enorm    Enorm_error   Enorm_exact   Enorm_Slope   
  Size    Required                                                         
___________________________________________________________________________
     2           7     1e+006              0             0             0   
     4           8     1e+006     0.00327463        1e+006             0   
     8           8     1e+006   -0.000201019        1e+006       4.02593   
    16           8     1e+006  -5.25297e-005        1e+006       1.93613   
    32           8     1e+006  -5.61744e-006        1e+006       3.22515   
    64           8     1e+006  -9.16809e-008        1e+006       5.93715   
   128           8     1e+006  -2.58752e-007        1e+006      -1.49688   
===========================================================================

PART B


Conductivity Quadratic Temperature Variation:     k(T) = A + B*T + C*T2


Figure 2: Steady-state temperature distribution across the given slab at a 128 element mesh size with quadratic FE basis. Here, the thermal conductivity varies quadratically with the temperature.

Figure 2 above shows the temperature distribution along the slab for the case of a thermal conductivity quadratic temperature dependence. This solution utilized a quadratic element basis with 128 mesh elements. The figure confirms that the solution meets the Dirichlet BCs at both ends of the slab. It is of interest to note that this resulting temperature profile is virtually indistinguishable from the results of Part A with the thermal conductivity linear temperature dependence.

Tables 9-12 list the convergence results for the case of a linear FE basis with and without the non-linear Jacobian contributions. The convergence slopes, as calculated from the Taylor series truncation error estimate from the middle node temperature and from the energy norm, differ very little. They all indicate a slope of 2 which follows the theory. As Tables 10 and 12 show, omission of the non-linear portions of the Jacobian matrix yields basically the same results and convergence rates but requires twice as many Newton Method iterations.

Table 9: Taylor Series truncation error convergence results for a linear (k=1) Finite Element basis.
===========================================================================================
  Mesh  Iterations          Tmid   Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)           (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           4     1642.2144             0              0        0.0000             0  
     4           4     1646.5571       4.34274        1.44758     1648.0047             0  
     8           4     1647.9219       1.36484       0.454947     1648.3769       1.66987  
    16           4     1648.3007      0.378764       0.126255     1648.4270       1.84936  
    32           4     1648.3990     0.0982439       0.032748     1648.4317       1.94686  
    64           4     1648.4238     0.0248245     0.00827485     1648.4320        1.9846  
   128           4     1648.4300    0.00622354     0.00207451     1648.4321       1.99596  
===========================================================================================

Table 10: Taylor Series truncation error convergence results for a linear (k=1) Finite Element basis without the non-linear Jacobian contribution.
===========================================================================================
  Mesh  Iterations          Tmid   Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)           (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           6     1642.2143             0             0        0.0000             0   
     4           7     1646.5570       4.34268       1.44756     1648.0045             0   
     8           9     1647.9220       1.36501      0.455004     1648.3770       1.66967   
    16           9     1648.3008      0.378811       0.12627     1648.4271       1.84936   
    32           9     1648.3991     0.0982626     0.0327542     1648.4318       1.94676   
    64           9     1648.4239     0.0248299    0.00827662     1648.4322       1.98457   
   128           9     1648.4301    0.00622491    0.00207497     1648.4322       1.99595   
===========================================================================================

Table 11: Energy Norm convergence results for a linear (k=1) Finite Element basis.
===============================================================================
  Mesh  Iterations            Enorm   Enorm_error   Enorm_exact   Enorm_Slope  
  Size    Required                                                             
_______________________________________________________________________________
     2           4     1.29256e+007             0             0             0  
     4           4     1.27407e+007      -61607.4  1.26791e+007             0  
     8           4     1.26864e+007      -18115.1  1.26683e+007       1.76592  
    16           4     1.26717e+007      -4894.64  1.26668e+007       1.88791  
    32           4     1.26679e+007      -1258.32  1.26667e+007        1.9597  
    64           4      1.2667e+007      -317.154  1.26667e+007       1.98825  
   128           4     1.26667e+007      -79.4584  1.26667e+007       1.99691  
===============================================================================

Table 12: Energy Norm convergence results for a linear (k=1) Finite Element basis without the non-linear Jacobian contribution.
===============================================================================
  Mesh  Iterations            Enorm   Enorm_error   Enorm_exact   Enorm_Slope  
  Size    Required                                                             
_______________________________________________________________________________
     2           6     1.29256e+007             0             0             0  
     4           7     1.27407e+007      -61607.8  1.26791e+007             0  
     8           9     1.26864e+007      -18114.3  1.26683e+007       1.76598  
    16           9     1.26717e+007      -4894.52  1.26668e+007       1.88789  
    32           9     1.26679e+007      -1258.28  1.26667e+007       1.95971  
    64           9      1.2667e+007      -317.144  1.26667e+007       1.98825  
   128           9     1.26667e+007       -79.456  1.26667e+007       1.99691  
===============================================================================

Tables 13-16 list the convergence results for the case of a quadratic FE basis with and without the non-linear Jacobian contributions. Again, both methods for convergence slope calculation differ very little and indicate a slope of 4, corresponding to theory. Also, Tables 14 and 16 show that omission of the non-linear portions of the Jacobian matrix only serve to double the iterations needed for the Newton Method to converge upon its solution.

Table 13: Taylor Series truncation error convergence results for a quadratic (k=2) Finite Element basis.
===========================================================================================
  Mesh  Iterations          Tmid   Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)           (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           4     1648.0586             0             0        0.0000             0   
     4           4     1648.3854      0.326833     0.0217889     1648.4072             0   
     8           4     1648.4279     0.0424628    0.00283085     1648.4307       2.94428   
    16           4     1648.4318    0.00391238   0.000260825     1648.4320       3.44008   
    32           4     1648.4321   0.000286442  1.90961e-005     1648.4321       3.77173   
    64           4     1648.4321  1.88277e-005  1.25518e-006     1648.4321       3.92731   
   128           4     1648.4321  1.19073e-006   7.9382e-008     1648.4321       3.98294   
===========================================================================================

Table 14: Taylor Series truncation error convergence results for a quadratic (k=2) Finite Element basis without the non-linear Jacobian contribution.
===========================================================================================
  Mesh  Iterations          Tmid   Tmid_change     Truncation    Tmid_exact      TS_Slope  
  Size    Required           (F)           (F)     Error Est.    (Est.) (F)                
___________________________________________________________________________________________
     2           8     1648.0586             0             0        0.0000             0   
     4           9     1648.3855      0.326869     0.0217912     1648.4073             0   
     8           9     1648.4280     0.0424936    0.00283291     1648.4308       2.94339   
    16           9     1648.4319    0.00391706   0.000261137     1648.4321        3.4394   
    32           9     1648.4322   0.000286856  1.91237e-005     1648.4322       3.77137   
    64           9     1648.4322  1.88565e-005   1.2571e-006     1648.4322        3.9272   
   128           9     1648.4322    1.193e-006  7.95333e-008     1648.4322       3.98239   
===========================================================================================

Table 15: Energy Norm convergence results for a quadratic (k=2) Finite Element basis.
===============================================================================
  Mesh  Iterations            Enorm   Enorm_error   Enorm_exact   Enorm_Slope  
  Size    Required                                                             
_______________________________________________________________________________
     2           4     1.26783e+007             0             0             0
     4           4     1.26681e+007      -680.957  1.26674e+007             0
     8           4     1.26668e+007      -86.2331  1.26667e+007       2.98125
    16           4     1.26667e+007      -7.87018  1.26667e+007       3.45377
    32           4     1.26667e+007     -0.574382  1.26667e+007       3.77632
    64           4     1.26667e+007    -0.0377152  1.26667e+007       3.92879
   128           4     1.26667e+007    -0.0023557  1.26667e+007       4.00092
===============================================================================

Table 16: Energy Norm convergence results for a quadratic (k=2) Finite Element basis without the non-linear Jacobian contribution.
===============================================================================
  Mesh  Iterations            Enorm   Enorm_error   Enorm_exact   Enorm_Slope  
  Size    Required                                                             
_______________________________________________________________________________
     2           8     1.26783e+007             0             0             0  
     4           9     1.26681e+007      -680.918  1.26674e+007             0  
     8           9     1.26668e+007      -86.2192  1.26667e+007        2.9814  
    16           9     1.26667e+007      -7.86849  1.26667e+007       3.45385  
    32           9     1.26667e+007     -0.574245  1.26667e+007       3.77635  
    64           9     1.26667e+007    -0.0377049  1.26667e+007       3.92884  
   128           9     1.26667e+007   -0.00236235  1.26667e+007       3.99646  
 ===============================================================================

Conclusion

In conclusion, this a convergence study has been completed for two cases where the thermal conductivity varies linearly and quadratically with temperature. The resulting Galerkin Weak Statement (GWSh) introduces non-linearities that have been solved using an iterative scheme based upon Newton's Method. For the conductivity linear temperature dependence case, the GWSh is seen to yield nodally exact solutions for all meshes. For all cases, omission of the non-linear portions of the Jacobian matrix does not prohibit convergence upon a solution but does require approximately twice as many iterations in Newton's Method.

Grade

10.6.04 Joe, this is without doubt the BEST lab 5 report that I have ever read! Your data presentation with commentary are very well organized, and the interpretation of the data is clear and concise. Others generated similar data comparisons wherein the quasi-Newton procedure yielded data that appear to generate a measureable convergence rate for the linear temperature variation case. This is the first time I have seen this result. I would like to replace the current lab 5 archive with your report, with your permission. Score for this report is 8/6 (much extra credit!). AJB

Troops: this lab turns out having MANY interesting aspects to it. Some of your results were right one, while others for whatever reason missed generating correct data or getting the interpretation correct. I suggest you look into Joe Tipton's report discussion section to see an exceptionally well presented data set interpretation. The totally new issue to me is the quasi-Newton solution iteration error process yielding data that appears to be mesh convergence for the linear temperature variation case. Overall, you should all readily see that understanding numerical solutions to problems in the engineering sciences is a challenging process! AJB