tdma.m

温度场的数值模拟程序

% Twodimensional heat conduction
% Finite Volume Method
% Line by Line TDMA with SOR
clear all;
x=[];y=[];T=[];Told=[];Su=[];Sp=[];ap=[];ae=[];aw=[];as=[];an=[];
Q=[];P=[];Fe=[];Fw=[];Fn=[];Fs=[];De=[];Dw=[];Dn=[];Ds=[];
great = 1e20;
lambda = 10; % thermal conductivity
alfa = 10; % heat transfer coefficient
dt = great; % Time step. If great steady state
density = 6000;% density
cp = 500;% heat capacity
Lx = 0.12; % length x-direction
Ly = 0.12; % length y -direction
Tfluid = 20; % Fluid temperature
Tinit = 50; % Initial guess and top- and bottom tempearature
u=0.2; %convection velocity of T in x direction
v=0.0; %convection velocity of T in y direction
%cv_x = input('Number of CVs in x-direction = ')
%cv_y = input('Number of CVs in y-direction = ')
cv_x=10;cv_y=10;
ni = cv_x+2; % grid points x-direction
nj = cv_y+2; % grid points y-direction
dx = Lx/cv_x;
dy = Ly/cv_y;
x(1) = 0;
x(2)=dx/2;
for i = 3:ni-1
   x(i)=x(i-1)+dx;
end;
x(ni)=Lx;
y(1) = 0;
y(2)=dy/2;
for j = 3:nj-1
   y(j)=y(j-1)+dy;
end
y(nj)=Ly;
% Initial values and coefficients
for i = 1:ni
  for j = 1:nj
    T(i,j) = Tinit;  %Initial temperature
    Told(i,j) = Tinit;
    T(i,1) = 50;
    T(i,nj) = 50
    Su(i,j)=0;       %Initial indendendent source term
    Sp(i,j)=0;       %Initial dependent source term
    Fe = density*u*dy;
    Fw = density*u*dy;
    Fn = density*v*dx;
    Fs = density*v*dx;
    De = lambda*dy/dx;
    Dw = lambda*dy/dx;
    Dn = lambda*dx/dy;
    Ds = lambda*dx/dy;
    ae(i,j) = Dw-Fw/2;
    aw(i,j) = Dw+Fw/2;
    an(i,j) = Dn-Fn/2;
    as(i,j) = Ds+Fs/2;
    dV = dx*dy;
    ap0 = density*cp*dV/dt;
    if i==2  % convective heat transfer boundary
       Su(i,j) = -Tfluid/(1/alfa+dx/(2*lambda))*dy/dV;
       Sp(i,j) = 1/(1/alfa+dx/(2*lambda))*dy/dV+density*u/2*dy/dV;
       aw(i,j) = 0;
    end;
    if i==ni-1 % insulated boundary 热流密度为0
       Su(i,j) = 0;
       Sp(i,j) = density*u/2*dy/dV;
       ae(i,j) = 0;
    end
    if j==2 % bottom boundary, given temperature
       Su(i,j) = (2*lambda*dx/dy+density*v*dx)*Tinit/dV;
       Sp(i,j) = -(2*lambda*dx/dy+density*v*dx)/dV;
       as(i,j) = 0;
    end   
    if j==nj-1 % top boundary, given temperature
       Su(i,j) = (2*lambda*dx/dy+density*v*dx)*Tinit/dV;
       Sp(i,j) = -(2*lambda*dx/dy+density*v*dx)/dV;
       an(i,j) = 0;
    end   
    ap(i,j) = ae(i,j)+aw(i,j)+an(i,j)+as(i,j)-Sp(i,j)*dV+ap0;
  end;
end;
%%%%%%%%%%%
maxerr = 1.0e-6;
maxit = 1000;
time=0;
maxtime =100;
omega=0.1;
while (time < (maxtime+dt/2))
 Told=T;  
 sumerr = 1;
 counter = 0;
 while (sumerr>maxerr&counter<maxit)
  %Sweep in j-direction  
  for i = 2:ni-1
    P(1) = 0;
    Q(1) = T(1); 
    for j = 2:nj-1
      a=ap(i,j)/omega;b=an(i,j);c=as(i,j);
      d=ae(i,j)*T(i+1,j)+aw(i,j)*T(i-1,j)+Su(i,j)*dV +ap0*Told(i,j)...
          +(1-omega)*a*T(i,j);
      P(j) = b/(a-c*P(j-1));
      Q(j) = (c*Q(j-1)+d)/(a-c*P(j-1));
    end;
    for j = nj-1:-1:2
      T(i,j) = P(j)*T(i,j+1)+Q(j);
    end;
  end;
  %Sweep in i-direction
  for j = 2:nj-1  
    P(1) = (lambda/(dx/2))/(lambda/(dx/2)-alfa);
    Q(1) = (-alfa*Tfluid)/(lambda/(dx/2)-alfa);
    for i = 2:ni-1
      a=ap(i,j)/omega;b=ae(i,j);c=aw(i,j);
      d=an(i,j)*T(i,j+1)+as(i,j)*T(i,j-1)+Su(i,j)*dV+ap0*Told(i,j)...
           +(1-omega)*a*T(i,j);
      P(i) = b/(a-c*P(i-1));
      Q(i) = (c*Q(i-1)+d)/(a-c*P(i-1));
    end;
    for i = ni-1:-1:2
      T(i,j) = P(i)*T(i+1,j)+Q(i);
    end;
  end;
  % Calculate boundary values
  for j = 2:nj-1
%%%%第三类边界条件Given heat convection
%%%%给定表面传热系数alfa
     T(1,j)=(lambda/(dx/2)*T(2,j)-alfa*Tfluid)/(lambda/(dx/2)-alfa);
%%%%第一类边界条件Given Temperature
     T(ni,j) = T(ni-1,j);
  end;
  sumres = 0;
  % Calculate residual
  for i = 2:ni-1
   for j = 2:nj-1    
    res(i,j) =  abs(ap(i,j)*T(i,j)-(ae(i,j)*T(i+1,j)+aw(i,j)*T(i-1,j)+...
      an(i,j)*T(i,j+1)+as(i,j)*T(i,j-1)+Su(i,j)*dV+ap0*Told(i,j)));
%%%%残差的定义还没有搞清楚
    sumres=sumres+res(i,j);
   end;
  end
  sumerr=sumres
  counter = counter + 1
 end;%end while
 time = time+dt;
end; %end time loop 
%%%% 全隐式离散化方程(不引入松弛因子时)的格式为 
%%%% ap*Tp=ae*Te+aw*Tw+b
%%%% b=Su*dV+ap0*Tp0
%%%% 引入松弛因之a后的格式为
%%%% (ap/a)*Tp=ae*Te+aw*Tw+(b+(1-a)*(ap/a)Tp0)
%%%% TDMA迭代式格式为
%%%% ai*T(i)=bi*T(i+1)＋ci*T(i-1)+di
%%%% 对应系数,有
%%%% ai=ap/a;bi=ae;ci=aw;di=(b+(1-a)*(ap/a)Tp0)
%%%% 扩展到二维情况，并采用交替方向线迭代方法(ADI)
%%%% 将y方向的 an(i,j)*T(i,j+1)+as(i,j)*T(i,j-1) 并入b中
%%%% P(i)=bi/(ai-ci*P(i-1)) Q(i)=(di+ci*Q(i-1))/(ai-ci*P(i-1))


%
pcolor(x,y,T');shading interp;xlabel('x');ylabel('y');title('Temperature distribution');colorbar;
%