example1b.m

来自「基于matlab计算流体力学(CFD)源代码。」· M 代码 · 共 42 行

%Temperature distribution in a rod
%See example 9
%Discretization method: Finite difference
% Solution method: SOR

clear all;
T = [];x=[];P=[];Q=[];
nn = input('Number of increment = ')
n = nn+1;% number of  grid points
L = 0.6;
dx = L/nn; % size of increment
Qprim = 50000; % heat source
lambda = 20; % thermal conductivity
x = 0:dx:L;
T(n) = 30;
% Solver TDMA
% first order appr. at boundary to get a tri-diagonal matrix
% dT/dx = 0 gives  T(1)=T(2) results in P(1)=1
P(1) = 1; 
Q(1) = 0;
a = 2; b = 1; c = 1;
d = Qprim*dx^2/lambda;
for k = 2:n-1
   P(k) = b/(a-c*P(k-1));
   Q(k) = (c*Q(k-1)+d)/(a-c*P(k-1));
end;
for k = n-1:-1:2
   T(k) = P(k)*T(k+1)+Q(k);
end;
T(1)=T(2);% Insulated end
% Analytical solution
for k = 1:n
   analytical(k)=Qprim*(L^2-x(k)^2)/(lambda*2)+T(n);
end;   
% Plot
hold on;
plot(x,T,'*');
plot(x,analytical,'o');
hold off;
legend('Numerical','Analytical',0);
title('Temperature distribution');