[转帖]matlab编写的流体计算和传热程序
%Temperature distribution in a rod%See example 9
%Discretization method: Finite difference
% Solution method: SOR
clear all;
a=[];b=[];c=[];d=[];x=[];T=[];analytical=[];
nn = input('Number of increment = ')
n = nn+1;% number ofgrid points
L = 0.6;
dx = L/nn; % size of increment
Qprim = 50000; % heat source
lambda = 20; % thermal conductivity
for k = 1:n % Set coefficients
a(k)=2.0;
b(k)=1.0;
c(k)=1.0;
d(k)=Qprim*dx^2/lambda;
T(k) = 30.0; % start value
x(k) = (k-1)*dx;
end; %
a(1) = 3;
b(1) = 4;
c(1) = 0;
% Solver SOR if omega = 1 Gauss-Siedel
omega=2/(1+sin(pi/nn));
counter = 0;
maxit = 200;
sumres = 1.;
maxres = 1.0e-5;
while ((sumres > maxres)&(counter < maxit) )
for k = 2:n-1 % SOR
T(k)=omega*(b(k)*T(k+1)+c(k)*T(k-1)+d(k))/a(k)+(1-omega)*T(k);
end; % SOR
d(1) =-T(3);
T(1)=(b(1)*T(2)+d(1))/a(1);% Insulated end
sumres = 0.0
for k = 2:n-1 % residual
res=abs(a(k)*T(k)-(b(k)*T(k+1)+c(k)*T(k-1)+d(k)));
sumres = sumres+res;
end % residual
summa = sumres
counter = counter +1 % counts number of iterations
end; %while
%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');
[ 本帖最后由 suffer 于 2006-10-9 19:49 编辑 ] %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 ofgrid 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 givesT(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'); %Temperature distribution in a rod
%See example 9
%Finite difference methof
%TDMA
clear all;
T = [];P=[];Q=[];analytical=[];
nn = input('Number of increment = ')
n = nn+1;% number ofgrid 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;
for k =1:n
T(k)=30;
end;
% Solver TDMA
% second order appr. dT/dx = 0 gives3T(1)=4T(2)-T(3)
% to get a tri-diagonal matrix the eqn. T(1) = T(1) is used
% For the results in P(1)=0 and Q(1)=T(1)
a = 2; b = 1; c = 1;
d = Qprim*dx^2/lambda;
counter = 0;
maxit = 200;
sumres = 1.;
maxres = 1.0e-5;
P(1) = 0;
while ((sumres > maxres)&(counter < maxit) )
Q(1) = T(1);
sumres = 0;
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));
sumres=sumres+abs(a*T(k)-(b*T(k+1)+c*T(k-1)+d));
end;
for k = n-1:-1:1
T(k) = P(k)*T(k+1)+Q(k);
end;
T(1)=(4*T(2)-T(3))/3;
counter=counter +1
end;
for k = 1:n
analytical(k)=Qprim*(L^2-x(k)^2)/(lambda*2)+T(n);
end;
hold on;
plot(x,T,'*');
plot(x,analytical,'o');
hold off;
legend('Numerical','Analytical',0);
title('Temperature distribution'); %Temperature distribution in a rod
%See example 9
%Finite volume method
%TDMA
clear all;
T = [];P=[];Q=[];ae=[];aw=[];ap=[];analytical=[];
nn = input('Number of control volumes = ');
n = nn+2;% number ofgrid points
L = 0.6;
dx = L/nn; % size of increment
Qprim = 50000; % heat source
lambda = 20; % thermal conductivity
x(1)=0;x(2)=dx/2;
for k = 3:n-1
x(k)=x(k-1)+dx;
end
x(n)=L;
T(n) = 30;
for k =1:n
T(k)=30;
end;
% Solver TDMA
% P(1)=0 and Q(1)=T(1)
for k = 2:n-1
d(k)=Qprim*dx;
ae(k)=lambda/dx;
aw(k)=lambda/dx;
if (k==2);aw(k)=0;end;
if (k==n-1);ae(k)=2*lambda/dx;end;
ap(k)=ae(k)+aw(k);
end;
P(1) = 0;
Q(1) = T(1);
sumres = 0;
for k = 2:n-1
P(k) = ae(k)/(ap(k)-aw(k)*P(k-1));
Q(k) = (aw(k)*Q(k-1)+d(k))/(ap(k)-aw(k)*P(k-1));
end;
for k = n-1:-1:1
T(k) = P(k)*T(k+1)+Q(k);
end;
T(1)=T(2);
for k = 1:n
analytical(k)=Qprim*(L^2-x(k)^2)/(lambda*2)+T(n);
end;
hold on;
plot(x,T,'*');
plot(x,analytical,'o');
hold off;
legend('Numerical','Analytical',0);
title('Temperature distribution'); % Transient temperature distribution in aninfinite plate with the thickness 2L.
% Initially temperature is Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Explicit method
%
clear all;
n = 20; % number of increments
nk = n+1; % number of node points
L = 0.05; % half thickness
dx = L/n % increment
dt =0.1 % time step
conductivity = 40; % thermal conductivity
density = 8000; % density
cp = 500 % Heat capacity
alfa = 1600; % heat transfer coefficient
diffusivity = conductivity/(density*cp); % thermal diffusivity
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
for k = 1:nk
x(k) = (k-1)*dx;
end;
% Explicit formulation
for k = 1:nk % coefficients
t(k) = Tinit;
told(k) = Tinit;
a(k) = 1;
b(k) = diffusivity*dt/dx^2;
c(k) = diffusivity*dt/dx^2;
d(k) = 1 - 2*diffusivity*dt/dx^2;
end;
% second order approximation at boundary
c(1) = 0;
b(1) = 4;
a(1) = 3 + 2*alfa*dx/conductivity;
a(nk) = 3;
b(nk) = 0;
c(nk) = 4;
time = 0
while (time < (maxtime+dt/2))
told = t;
d(1) = 2*Tfluid*alfa*dx/conductivity - t(3);
t(1)=(b(1)*t(2)+d(1))/a(1);
for k = 2:n
t(k) = (b(k)*told(k+1)+c(k)*told(k-1)+d(k)*told(k))/a(k);
end;
d(nk)=-t(nk-2);
t(nk) =(c(nk)*t(nk-1)+d(nk))/a(nk);
time = time +dt;
end; %while
%
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution'); % Transient temperature distribution in aninfinite plate with the thickness 2L.
% Initially temperature is Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Explicit method
%
clear all;
t=[];ap=[];aw=[];ae=[];d=[];
n_cv = 30; % number of control volumes
ni = n_cv+2; % number of node points
L = 0.05; % half thickness
dx = L/n_cv % increment
dt =0.1 % time step
conductivity = 40; % thermal conductivity
density = 8000; % density
cp = 500 % Heat capacity
alfa = 1600; % heat transfer coefficient
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
x(1) = 0;
x(2) = dx/2;
for k = 3:ni-1
x(k) = (k-1)*dx;
end;
x(ni) = L;
% Explicit formulation
for k = 1:ni % coefficients
t(k) = Tinit;
told(k) = Tinit;
ap(k) = density*cp*dx/dt;
ae(k) = conductivity/dx;
aw(k) = conductivity/dx;
if (k==2);aw(k)=0;end;
if (k==ni-1);ae(k)=0;end;
d(k) = ap(k)- aw(k)-ae(k);
end;
time = 0;
while (time < (maxtime+dt/2))
told = t;
for k = 2:ni-1
dd = d(k)*told(k);
if k == 2;
dd = d(k)*told(k)+(Tfluid-told(k))/(1/alfa+2*dx/conductivity);
end;
t(k) = (ae(k)*told(k+1)+aw(k)*told(k-1)+dd)/ap(k);
end;
time = time +dt;
end; %while
% Boundaries
t(ni)=t(ni-1)
t(1)=(alfa*Tfluid+2*conductivity/dx*t(2))/(alfa+2*conductivity/dx);
%
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution'); % 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; %whiletime loop
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution'); % Transient temperature distribution in an infinite plate with the thickness 2L.
% Initially the plate is at Tinit when it is exposed to a fluid with
% the temperature Tfluid.
% Implicit scheme
% Finite volume method
% Solver TDMA
t=[];ap=[];aw=[];ae=[];d=[];Q=[];P=[];
n_cv = 30; % number of control volumes
ni = n_cv+2; % number of node points
L = 0.05; % half thickness
dx = L/n_cv % increment
dt =1 % time step
conductivity = 40; % thermal conductivity
density = 8000; % density
cp = 500 % Heat capacity
alfa = 1600; % heat transfer coefficient
maxtime = 100; % total time
Tinit = 1; % initial temperature (dimensionless)
Tfluid = 0; % fluid temperature (dimensionless)
x(1) = 0;
x(2) = dx/2;
for k = 3:ni-1
x(k) = (k-1)*dx;
end;
x(ni) = L;
% Implicit formulation
for k = 1:ni % coefficients
t(k) = Tinit;
told(k) = Tinit;
ae(k) = conductivity/dx;
aw(k) = conductivity/dx;
if (k==2);aw(k)=0;end;
if (k==ni-1);ae(k)=0;end;
d(k) = density*cp*dx/dt;
ap(k)=aw(k)+ae(k)+density*cp*dx/dt;
end;
time = 0;
maxres = 1.e-5;
maxit = 10;
while (time < (maxtime+dt/2))
told = t;
P(1) = 0;
Q(1) = t(1);
for k = 2:ni-1
dd = d(k)*told(k);
a = ap(k);
if k == 2;
a = ap(k)+1/(1/alfa+2*dx/conductivity);
dd = d(k)*told(k)+Tfluid/(1/alfa+2*dx/conductivity);
end;
P(k) = ae(k)/(a-aw(k)*P(k-1));
Q(k) = (aw(k)*Q(k-1)+dd)/(a-aw(k)*P(k-1));
end;
for k = ni-1:-1:2
t(k) = P(k)*t(k+1)+Q(k);
end;
time = time +dt;
end; %whiletime loop
% Boundaries
t(ni)=t(ni-1)
t(1)=(alfa*Tfluid+2*conductivity/dx*t(2))/(alfa+2*conductivity/dx);
%
plot(x,t,'*');ylabel('Temperature');xlabel('x-pos');
title('Transient temperature distribution'); % Assignment 6 MMV042A 2003
% Boundary layer eqn.
% Tangential flow along a flat plate
% Marching procedure in x-direction
% using TDMA in y-direction
clear all;
x=[];y=[];P=[];Q=[];U=[];V=[];T=[];
viscosity = 18.1e-6; % viscosity
density = 1.19 ;
kvisc = viscosity/density; % kinematic viscosity
Pr = 0.72; % Prandtl number
Umax = 5; % free stream velocity
maxres = 1e-5; % max residual
maxit = 3;
nx = input('Number of increment in x-dir= ')
ny = input('Number of increment in y-dir= ')
ni = nx +1;
nj = ny +1;
% Increment length in y-dir 0.05 - 0.2 mm
% Increment length in x-dir 0.1 - 1.0 mm
dx = input('Increment length in x-dir (mm)= ')
dy = input('Increment length in y-dir (mm)=')
dx = dx/1000;
dy = dy/1000;
x=0:dx:((ni-1)*dx)
y=0:dy:(nj-1)*dy
% Boundary and initial values
for i=1:ni
for j=1:nj
U(i,j)=Umax;
if j==1; U(i,j)=0; end;
V(i,j) = 0;
T(i,j) = 0;
if j==1; T(i,j)=1.0; end;
end;
end;
counter = 0;
for i = 2:ni
counter = 0;
sumres = 1;
while ((sumres > maxres)&(counter < maxit) )
P(1)=0;Q(1)=0;
sumres=0;
%% momentuum eqn
for j = 2:nj-1
a = U(i,j)/dx+2*kvisc/dy^2;
b = kvisc/dy^2-V(i,j)/(2*dy);
c = kvisc/dy^2+V(i,j)/(2*dy);
d = U(i,j)*U(i-1,j)/dx;
P(j)=b/(a-c*P(j-1));
Q(j)=(c*Q(j-1)+d)/(a-c*P(j-1));
sumres=sumres+abs(a*U(i,j)-(b*U(i,j+1)+c*U(i,j-1)+d));
end;
for j=nj-1:-1:2
U(i,j) = P(j)*U(i,j+1)+Q(j);
end;
P(1)=0;Q(1)=0;
%%continuity eqn.
for j = 2:nj-1
a = 1/dy;
b = 0;
c = 1/dy;
d = -((U(i,j)-U(i-1,j))+(U(i,j-1)-U(i-1,j-1)))/(2*dx);
P(j)=b/(a-c*P(j-1));
Q(j)=(c*Q(j-1)+d)/(a-c*P(j-1));
sumres=sumres+abs(a*V(i,j)-(b*V(i,j+1)+c*V(i,j-1)+d));
end;
for j=nj-1:-1:2
V(i,j) = P(j)*V(i,j+1)+Q(j);
end;
counter = counter +1; % counts number of iterations
end;
TotalResidual = sumres
end;
for i = 2:ni % Temperature field
P(1)=0;Q(1)=T(i,1);% TDMA
% for j = 2:nj-1
% a =???
% b =???
% c =???
% d =???
% 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;
for i=2:ni
% calc. Nusselt number and friction factor
end;
subplot(2,1,1);pcolor(x,y,U');shading interp;title('U-velocity');colorbar;
% subplot(2,1,2);pcolor(x,y,T');shading interp;title('Temperature');colorbar; % Twodimensional heat conduction
% Finite Volume Method
% SOR
clear all;
x=[];y=[];T=[];Told=[];Su=[];Sp=[];ap=[];ae=[];aw=[];as=[];an=[];
great = 1.e20;
lambda = 10; % thermal conductivity
alfa = 10; % heat transfer coefficient
dt = great; % Time step. If great stedy 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
%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
ae(i,j) = lambda*dy/dx;
aw(i,j) = lambda*dy/dx;
an(i,j) = lambda*dx/dy;
as(i,j) = lambda*dx/dy;
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;
aw(i,j) = 0;
end;
if i==ni-1 % insulated boundary
ae(i,j) = 0;
end
if j==2 % bottom boundary, given temperature
as(i,j)=2*lambda*dx/dy;
end
if j==nj-1 % top boundary, given temperature
an(i,j)=2*lambda*dx/dy;
end
ap(i,j) = ae(i,j)+aw(i,j)+an(i,j)+as(i,j)-Sp(i,j)*dV+ap0;
end;
end;
%%%%%%%%%%%
maxres = 1.0e-6;
maxit = 100;
time=0;
maxtime=100;
s=(cos(pi/cv_x)+(dx/dy)^2*cos(pi/cv_y))/(1+(dx/dy)^2);
omega =2/(1+sqrt(1-s^2));omega=1;
while (time < (maxtime+dt/2))
Told=T;
sumres = 1;
counter = 0;
while (sumres>maxres&counter<maxit)
sumres = 0;
for i = 2:ni-1
for j = 2:nj-1
T(i,j)=omega*(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))/ap(i,j)+(1-omega)*T(i,j);
res =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;
end;
end;
for i = 2:ni-1
for j = 2:nj-1
res =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;
end;
end
sumerr=sumres
counter = counter + 1
end;
time = time +dt;
end;
% Calculate boundary values
for j = 2:nj-1
T(1,j)=(alfa*Tfluid+lambda/(dx/2)*T(2,j))/(alfa+lambda/(dx/2));
T(ni,j) = T(ni-1,j);
end;
%
pcolor(x,y,T');shading interp;xlabel('x');ylabel('y');title('Temperature distribution');colorbar;
% % 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=[];
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
%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
ae(i,j) = lambda*dy/dx;
aw(i,j) = lambda*dy/dx;
an(i,j) = lambda*dx/dy;
as(i,j) = lambda*dx/dy;
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;
aw(i,j) = 0;
end;
if i==ni-1 % insulated boundary
ae(i,j) = 0;
end
if j==2 % bottom boundary, given temperature
as(i,j)=2*lambda*dx/dy;
end
if j==nj-1 % top boundary, given temperature
an(i,j)=2*lambda*dx/dy;
end
ap(i,j) = ae(i,j)+aw(i,j)+an(i,j)+as(i,j)-Sp(i,j)*dV+ap0;
end;
end;
%%%%%%%%%%%
maxres = 1.0e-6;
maxit = 100;
time=0;
maxtime = 100;
omega=1;
while (time < (maxtime+dt/2))
Told=T;
sumres = 1;
counter = 0;
while (sumres>maxres&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) = 0;
Q(1) = T(1);
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;
sumres = 0;
% Calculate residual
for i = 2:ni-1
for j = 2:nj-1
res =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;
end;
end
sumerr=sumres
counter = counter + 1
end;%end while
time = time+dt;
end; %end time loop
% Calculate boundary values
for j = 2:nj-1
T(1,j)=(alfa*Tfluid+lambda/(dx/2)*T(2,j))/(alfa+lambda/(dx/2));
T(ni,j) = T(ni-1,j);
end;
%
pcolor(x,y,T');shading interp;xlabel('x');ylabel('y');title('Temperature distribution');colorbar;
% 这个有些看不明白。 搂主能否给些简单的说明呢 原帖由 无水1324 于 2006-10-12 09:09 发表
搂主能否给些简单的说明呢
程序开头的注解都有简单的说明 太长了
有简单的没有？
页:
[1]
2