example2b.m

基于matlab计算流体力学(CFD)源代码。

% Transient temperature distribution in an
% infinite plate with the thickness 2L.
% Initially the plate has the temperature
% Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Implicite scheme
% Solver Gauss-Siedel
clear all;
n = 20; % number of increments
nk = n+1; % number of node points
L = 0.05; % half thickness
dx = L/n; % increment
dt = 1.0 % time step
lambda = 40; % thermal conductivity
alfa = 1600; % heat transfer coefficient
diffusivity = 1.0e-6; % thermal diffusivity
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
x=[];a=[];b=[];c=[];d=[];t=[];told=[];
x = 0:dx:L;

% Fully Implicit formulation
for k = 1:nk % coefficients
   t(k) = Tinit;
   told(k) = Tinit;
   a(k) = 1+2*diffusivity*dt/dx^2;
   b(k) = diffusivity*dt/dx^2;
   c(k) = diffusivity*dt/dx^2;
   d(k) = 1;
end;
% second order approximation at boundary
c(1) = 0;
b(1) = 4;
a(1) = 3 + 2*alfa*dx/lambda;
a(nk) = 3;
b(nk) = 0;
c(nk) = 4;
time = 0;
maxres = 1.e-5;
maxit = 10;
while (time < (maxtime+dt/2))
   told = t;
   sumres = 1;
   counter = 0;
   while ((sumres > maxres)&(counter < maxit) )
     d(1) = 2*Tfluid*alfa*dx/lambda - t(3);
     t(1)=(b(1)*t(2)+d(1))/a(1);
     for k = 2:n
        d(k)=told(k);
        t(k) = (b(k)*t(k+1)+c(k)*t(k-1)+d(k))/a(k);
     end;
     d(nk)=-t(nk-2);
     t(nk) =(c(nk)*t(nk-1)+d(nk))/a(nk);
     sumres = 0.0;
     for k = 2:n % residual
       res=abs(a(k)*t(k)-(b(k)*t(k+1)+c(k)*t(k-1)+told(k)));
       sumres = sumres+res;
     end % residual
     summa = sumres
     counter = counter +1 % counts number of iterations    
   end; %while convergence loop   
   time = time +dt;
end; %while  time loop
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution');