fdm之二维静态热传导

function heat_2d

% Geol 5200, Humphrey 2013

% 2D, steady state temperature.  Simplest approach
% full (not sparse) matrices, with BCnodes included in matrix
clear
% solve temperature field in crust with a constant temperatures
BCtop = 0;
BCbottom = 5;
BCright = 3;
BCleft  = 4;


nodes_x = 36;       %width
nodes_y = 25;       %depth


nodes   = nodes_y*nodes_x;


% all BCs at boundaries are fixed
% steady state and constant K
% values of the diagonal and off diagonal coeff
alpha = -4;
beta  =  1;


% make a diagonal matrix of correct size, with '1's on diagonal
v(1:nodes)=1;
A = diag(v);
% better would be to use sparse matrices as follows
% A = speye(nodes,nodes)


% construct the non-trivial bands in the matrix, using a 5 point FD kernel
for i=nodes_x+1:nodes-nodes_x;
    % check if we are at a boundary node, if not then do:
    if((mod(i-1,nodes_x)~=0) && (mod(i,nodes_x)~=0))    %leave BC nodes as '1's  
        A(i,i)=alpha;
        A(i,i-1)=beta;
        A(i,i+1)=beta;
        A(i,i-nodes_x)=beta;
        A(i,i+nodes_x)=beta;
    end
end


% make the known vector, with BC's in correct locations
B = zeros(nodes,1);
B(1:nodes_x)                         =BCtop;     %
B(nodes-nodes_x:nodes)               =BCbottom;    
B(nodes_x+1:nodes_x:nodes-nodes_x-1) =BCleft;     
B(2*nodes_x:nodes_x:nodes-1)         =BCright;     %bottom BC


%solve
T = A\B;


%Plot a contour plot of the solution, note the solution is a vector, and needs to be 
% reshaped into a matrix to plot
B = reshape(T,nodes_x,nodes_y);
%[C, h]=contour(B,[ 1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 3.5 3.75]);
[C, h]=contour(B');
clabel(C,h);
axis ij;
xlabel('Bottom');
ylabel('Left');
title('Steady Temperature field, contoured');


% make a new figure and make a surface plot of the temperature field
figure;
surf(B);
axis ij;
title('Temperature as a surface, note upside down');
view(140,30)