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);

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.