FDM之二维静态热传导--含有不同传导系数K

function head2d_K

% 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 = 1000;
BCright = 0;
BCleft  = 0;


nodes_x = 40;       %width x-direction
nodes_y = 80;       %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)
KK=ones(nodes_y,nodes_x);
K1 = 1.; K2 = 5.; K3= 1.5;
j = 1;k=1;
for j = 1:nodes_y
    for k = 1:nodes_x
        if (j > nodes_y/3&& j <  nodes_y*2/3&&k<nodes_x*2/3)
            KK(j,k) = K2;
        else
            KK(j,k) = K1;
        end
         
    end
   
    
end
% 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  
        x=mod(i-1,nodes_x);
        y=floor((i-1)/nodes_x);
        A(i,i)=-(KK(y,x)+KK(y+1,x)+KK(y,x+1)+KK(y+1,x+1));
        A(i,i-1)=0.5*(KK(y,x)+KK(y+1,x));
        A(i,i+1)=0.5*(KK(y,x+1)+KK(y+1,x+1));
        A(i,i-nodes_x)=0.5*(KK(y,x)+KK(y,x+1));
        A(i,i+nodes_x)=0.5*(KK(y+1,x)+KK(y+1,x+1));
    end
    %left boundary and right
    if (mod(i-1,nodes_x)==0)
        A(i,i) = -1;
        A(i,i+1)= 1;
    end
    if (mod(i,nodes_x)==0)
        A(i,i) = 1;
        A(i,i-1)= -1;
    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) =0;     
B(2*nodes_x:nodes_x:nodes-1)         =0;     


%solve
T = inv(A)*B;


dx=1:nodes_x;
dy=1:nodes_y;
[X,Y]=meshgrid(dx,dy);
ii=1;
for i=1:nodes_y
    for j=1:nodes_x
        TT(i,j)=T(ii);
        ii=ii+1;
    end
end


figure(111)
pcolor(X,Y,TT)
shading flat
colorbar
axis ij
figure(222)
[C,h]=contourf(TT);
colorbar
%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)


figure;
contourf(KK);

title('Map of Conductivities');

结果如下：