fdm之一维静态热传导不同导热系数K

function heat_neil_varK
% program to solve temperature in the earth
% note this does not use realistic values, it is mainly to 
% demonstrate correct usage of variable conductivity (K) layers
% surface BC = 0, with constant flux throughout, applied as a BC at the top
% radioactive decay not included
% this version calcs heat input at center of element
clear


% Using a Matrix implicit finite difference approach
q=.065;             % heat flux
K=2.;               % conductivity
Z = 100000;         % thickness of problem
n = 11;             %number of nodes in Z
delz = Z/(n-1);     % delta z, length of problem / number of nodes
z = 0:delz:Z;       % 'z' is a vector for plotting


% construct the K vector
K1 = 2.; K2 = 3.; K3= 1.5;
zK12 = 20001; zK23 = 50001;  % transition depths


j = 1;
for zz = 0:delz:Z
    if    ( zz < zK12) Kvec(j) = K1;
    elseif (zz < zK23) Kvec(j) = K2;
    else               Kvec(j) = K3;
    end
    j = j+1;
end


% matrix construction
a = ones(1,n).*Kvec;      % make the A matrix, 'a' is the diagonal, 'b' is the off-diagonal
b = ones(1,n-1)*-1.*(Kvec(1:n-1)+Kvec(2:n));
c = ones(1,n-2).*Kvec(2:n-1);
A= diag(a,0) + diag(b,-1) + diag(c,-2);


%boundary conditions
A(1,1) = 1;         % adjust the A matrix for the top BC and for the gradient BC at the Bottom
A(1,2) = 0;
A(2,1) = -1;         % gradient BC
A(2,2) = 1;




C=zeros(n,1);       % make the BC vector
C(1) = 0;       %BC
C(2) = delz*q/K1;






%T=inv(A)*C;         % multi by the inverse to solve
T = A\C;


depth = z/1000;
plot(T,depth,'r*')
l1=line(T,depth)       % put a line through the output of a different color
set(l1,'color','green');
axis ij             % put the zero at the top of the plot
xlabel('Temperature')
ylabel('Depth [kilometers]')
handlet=title('Geothermal Temperature Profile, radioactive decay')
line([0 T(n)],[zK12 zK12]/1000);
line([0 T(n)],[zK23 zK23]/1000);