Objectives

In Lab #9, a uniform mesh refinement study was performed for the case of 2-D steady-state heat conduction through a square domain heated uniformly by a source. Removal of one quarter of the square (via the imposition of Dirichlet-style conditions) created an "L" shaped domain was created, and the temperature solution was solved directly using the resulting diffusion and data matrices. For the current lab, we desire to modify the Lab #9 MATLAB program code to solve for the steady-state temperature distribution via Newton's iterative proceedure. This is motivated by the desire to develop the capabilities to handle non-linear problems as in the case of an unsaturated aquifer.

Problem Statement

This lab considers 2-D steady-state heat conduction through an "L" shaped region. This region has a thermal conductivity of unity (k = 1) as well as a unity source term (s = 1). Dirichlet conditions (Tb = 0) are imposed along the entire boundary. The domain consists of a 2x2 unit square with coordinates running from -1 to 1 on each axis. Then, to create the "L" shape, one quarter of the domain is forced to a fixed temperature of zero, equivalent to the Dirichlet BC implimentation.

The resulting characteristic equation (Laplace's equation) is well known and has been mentioned in previous labs

with Dirichlet boundary conditions imposed at all outer surfaces.

either the natural coordinate system or the tensor product system could be used to create the element mesh resulting in triangular or quadrilateral (without quadriture) elements.

It is now desired to solve for the steady-state temperature distribution using Newton's iterative method. The Newton Method solution format is

where {FQ} signifies the element-wise assembly of the Galerkin Weak Statement, Se({WS}h). The Jacobian is the partial derivative of WSh with respect to the temperature, Q.

Since Dirichlet conditions yield fixed temperatures around the edge of the domain, the respective nodes of the Jacobian matrix should be 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.

Discussion of Results

The MATLAB program code below impliments the given problem and solves for the steady-state temperature distribution using Newton's iterative method.

 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #10                                                                         %
% 2-D, steady-state, conduction using triangular, elements (k=1).  A "L			 %
% shape region is created by calling either the do_Lshape2D_dlist			 %
% subroutine.  The domain is defined on a 2x2 unit square, with coordinates	         %
% running from -1 to 1 on each axis.  k*Laplacian(q) - s = 0 is solved			 %
% using Newton's iterative method, with k=1 and s=1 on the interior.  A			 %
% Dirichlet condition (q=0) is imposed along the entire boundary.			 %
% 											 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

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

Mesh = [8 16 32 64 128];

for i = 1:size(Mesh,2)

	% Declare global variables.  These variables must be made available as
	% global variables for the assembly routine to work.
	
	elwide = Mesh(i);	% number of elements along the width of the grid
	eltall = Mesh(i);	% number of elements along the height of the grid
						
	% Construct the grid seeds.  Be sure that X1seed and X2seed both include 0.
	X1seed=linspace(-1,1,elwide+1);
	X2seed=linspace(1,-1,eltall+1);
	
	% construct the nodal coordinate lists
	[X1 X2] = meshgrid(X1seed, X2seed);
	nnodes = prod(size(X1));
	
	% construct the element connectivity table
	elcon = delaunay(X2,X1);
    
    % save element mesh size
    elements(i,1) = size(elcon,1);
	
	% Scan through all nodes and select the ones for which Dirichlet
	% conditions will be enforced. Also, select the nodes outside the
	% computational domain for blanking.  This will be enforced in exactly
	% the same manner as Dirichlet B.C's.
	diri_list = do_Lshape2D_dlist(X1,X2);
	
	% calculate element zeta's and determinants.
	[ZETAe,DETe]=do_ZETAe_2Dt(X1,X2,elcon);
	
    % generate element diffusion matrices, [DIFF]e
    DIFFe = do_DIFFe_2Dt(ZETAe);

    % initial temperature quess
        Q = 0.25 * ones((Mesh+1)^2,1);
        
        % include Dirichlet conditions
        for k=1:size(Q,1)
            if diri_list(k,1) == 0
                % this col is not diri, skip it.
            else
                Q(k,1) = 0;
            end
        end
    
    for j=1:maxit
        % assemble residual matrix [RES]

            % assemble global diffusion matrix, [DIFF]
            DIFF = asres([],[],[],-1,DIFFe,Q,elcon,DETe,nnodes);
            
            % assemble the residual
            B200L = 1/24 * [2 1 1; 1 2 1; 1 1 2];
            SRC = ones(size(X1));
            b = asres([],[],[],1,B200L,SRC,elcon,DETe,nnodes);
            
            % combine
            RES = DIFF - b;
    
        % assemble the Jacobian [JAC]
            JAC = asjac([],[],[],-1,DIFFe,[],elcon,DETe,nnodes);
            
            % enforce Dirichlet conditions to keep dQ~0
                % THE DESIRED SUBROUTINE WOULD USE THE MATRIX diri_list TO
                % DETERMINE WHICH ENTRIES IN THE JACOBIAN MATRIX REPRESENTED
                % NODES WHERE THE DIRICHLET CONDITIONS SHOULD BE ENFORCED.
                % THEN THOSE NODES WOULD BE SET TO A HIGH NUMBER, ~1E15, TO
                % 'INFORM' THE ITERATION PROCESS THAT dQ IS APPROXIMATELY ZERO.
                %  THUS, ENTRIES IN {Q} THAT ARE RESTRAINED BY THE DIRICHLET
                %  BCs ARE NOT UPDATED.
            
        % compute solution
        dQ = JAC\(-1*(RES));   % linear solve for dQ
        Q = Q + dQ;            % update Q
        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;
	
	% save max temperatures
	Qmax(i,1) = max(Q);
    
    % calculate and display energy norm. 
    % res must be recalculated using the converged solution      
            DIFF = asres([],[],[],-1,DIFFe,Q,elcon,DETe,nnodes);
            B200L = 1/24 * [2 1 1; 1 2 1; 1 1 2];
            SRC = ones(size(X1));
            b = asres([],[],[],1,B200L,SRC,elcon,DETe,nnodes);
            RES = DIFF - b;
	
    enorm(i,1) = 0.5 * Q' * RES;
	
end

delta_Qmax = Qmax(2:end) - Qmax(1:end-1);
slope_Qmax = log10(delta_Qmax(1:end-1) ./ delta_Qmax(2:end)) / log10(2);
delta_enorm =  enorm(2:end) - enorm(1:end-1);
slope_enorm = log10(delta_enorm(1:end-1) ./ delta_enorm(2:end)) / log10(2);

ANSWER(:,1) = elements;
ANSWER(:,2) = Qmax;
ANSWER(:,3) = [0;delta_Qmax];
ANSWER(:,4) = [0;0;slope_Qmax];
ANSWER(:,5) = enorm;
ANSWER(:,6) = [0;delta_enorm];
ANSWER(:,7) = [0;0;slope_enorm];

fprintf('Elements   Qmax  (Qmax)change  (Qmax)slope    ||Q||      change||Q||     Slope(Q)\n\n')
fprintf('%6.0f  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g\n',ANSWER');

% plot the solution
figure(1)
trisurf(elcon,X1,X2,Q); view(0,0);
figure(2)
trisurf(elcon,X1,X2,Q); view(3);

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

As seen, the only thing lacking to the code is a subroutine to apply a very large number at every node on the Jacobian matrix that exists on the domain boundary. An implimentation of this is currently beyond the capability of the author. Were the appropriate subroutine supplied, however, the problem could be solved and the solution validated against the direct solve results from Lab #9.

With the iterative process of Newton's method defined and implimented, the above code can now handle non-linear problem definitions. A prime example would be that of an unsaturated aquifer.

In this case, what has previously been defined as the coefficient of thermal conductivity becomes the porepressure which is also the dependent variable. The above code could then easy accomodate this kind of problem via a rather simple modification to the residual and jacobian matrices.

Conclusion

In conclusion, the MATLAB program code from Lab #9 has been modified to solve for the steady-state temperature distribution via Newton's iterative proceedure. This cannot be run until a subroutine is developed to handle the necessary manipulations to the boundary nodes in the Jacobian matrix due to the Dirichlet BCs. Finally, it is recognized that this new solution code can be utilized in the solution of a non-linear problem such as the case of an unsaturated aquifer.

Grade

11.11.04 Joe, I did not see the detail of your jacobian for this lab, which of course is THE main issue. Did you recognize that matrices with indices alternated, e.g., [B3101] are generated? You did recognize that the BC link had to be added to handle Dirichlet nodes - GOOD! Score is 3/3 AJB

Dr. Baker, You might want to reword the lab assignment as listed on the course website. After talking with some other students, I think the wording confused us. You could construe the assignment as asking the student to just rewrite Lab 9 as a newton iterative statement as a motivation for a non-linear problem statement similar to an unconfined aquifer. Regards, Joseph

Joe, indeed I will totally reword this to request the action using template essence nomenclature, hence modify what is presented on courseware page HTn.29. Thanks for the comments. AJB

Troops - some of you got close on this, but no one really got the non-linearity into the jacobian. Look back at the temperature dependent conductivity example for guidance on this issue. You should have had asjac contributions that employed the matrix set {B3K0J] as well as the standard [B30KJ]. Only Joe noticed that you also had to hit the diagonals of the JAC with an E10 to enforce the fixed BCs. AJB