方腔内自然对流MATLAB程序数值传热学

%==========================================================================
un=0.0; us=0.0; vw=0.0; ve=0.0; Tw=100.0; Te=0.0; T(1,1:ny+1)=100; TT(1,1:ny+1)=100; %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vw-v(2,1:ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1);
%========================================================================== uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,1:ny+1)+u(1:nx+1,2:ny+2)); vv(1:nx+1,1:ny+1)=0.5*(v(1:nx+1,1:ny+1)+v(2:nx+2,1:ny+1));
自然对流传热
问题描述：
一个二维矩形腔体，高 H 1 ，宽 L 1 ，物理模型如图 1。
绝热壁面
热 壁 面 T=100
冷 壁 面 T=0
绝热壁面 图 1 物理模型
上下壁面为绝热壁面，左壁面为热壁面 T 100 ，右壁面为冷壁面 T 0 。腔体内为初始 温度 T0 0 的空气，计算由此引起的自然对流传热。
(6)
vit, j vit, j 1 vit, j vit, j 1 2 2
2
2

vit+1, j 2vit, j vit1, j h2
vit, j 1 2vit, j vit, j 1 h2
T u T 2T t
(2)
其中， 是热扩散系数。 边界条件： 由无滑移边界条件得，四个壁面上的速度均为零，即： un us uw ue 0 。 在热壁面上 T 100 ，在冷壁面上 T 0 ，在上下绝热壁面上处
T 0。 y
数值处理：
区域离散化如图 2 所示。 对于动量守恒方程，在不考虑压力的情况下先计算出一个临时速度
(4)
根据求得的临时速度，在考虑压力的情况下，得修正速度
1 uin, j ui , j
dtBiblioteka Baidu

1 pi 1, j pi , j dx
2
(5)
方程(3)对应的离散方程如下：
t t t t t t t t 1 ui , j ui 1, j ui , j ui 1, j vi , j vi 1, j ui , j 1 ui , j dt h 2 2 2 2 t t t t t t t t t t ui +1, j 2ui , j ui 1, j ui , j 1 2ui , j ui , j 1 vi , j 1 vi 1, j 1 ui , j ui , j 1 2 2 2 h h2 Ti , j Ti 1, j gx T0 2 t t t t t t t t vi , j vi , j 1 ui , j ui , j 1 vi , j vi 1, j ui 1, j uit1, j 1 vit, j vit1, j dt h 2 2 2 2 t t uit, j ui , j 2
u u n un un 2un g T T0 dt
(3)
用 SOR（超松弛迭代）法求压力场
1 1 pi, j1 c pi 1, j pi 1, j pi , j 1 pi , j 1
h ui , j ui1, j vi, j vi1, j 1 pi, j dt
%========================================================================== for i=2:nx for j=2:ny T(i,j)=TT(i,j)+dt*( -(0.25/h)*( (uu(i,j)+uu(i+1,j))*(TT(i,j)+... TT(i+1,j))-(uu(i-1,j)+uu(i,j))*(TT(i-1,j)+TT(i,j))+... (vv(i,j)+vv(i,j+1))*(TT(i,j)+TT(i,j+1))-(vv(i,j)+... vv(i,j-1))*(TT(i,j)+TT(i,j-1)) )+(alpha/h^2)*(TT(i+1,j)+... TT(i-1,j)+TT(i,j+1)+TT(i,j-1)-4*TT(i,j)) ); end end for i=2:nx T(i,1)=T(i,2); T(i,ny+1)=T(i,ny); end TT=T; %========================================================================== end %========================================================================== for i=1:nx+1 for j=1:ny+1 velocity(i,j)=sqrt(uu(i,j)*uu(i,j)+vv(i,j)*vv(i,j)); end end %========================================================================== startx=h:10*h:0.5; starty=zeros(size(startx))+0.5; streamline(xh,yh,flipud(rot90(uu)),flipud(rot90(vv)),startx,starty); clf(figure) quiver(xh,yh,flipud(rot90(uu)),flipud(rot90(vv)),'r'),title('Velocity'),axis equal,axis([0 L 0 H]); clf(figure) contour(xh,yh,flipud(rot90(T))),title('Temperature map'),axis equal,axis([0 L 0 H]); clf(figure) [X,Y,Z]=griddata(xh,yh,flipud(rot90(T)),linspace(0,1,nx+1)',linspace(0,1,ny+1),'v4'); pcolor(X,Y,Z);shading interp;title('Temperature contour'); colorbar,colormap(jet); %========================================================================== %求温度场
(7)
Ti , j Ti , j 1 gx T0 2
图 2 交错网格
数值结果：
速度场
图 3 速度场
温度等值线和温度云图
图 4 温度等值线
图 5 温度云图
由图 3-5 可知，
程序：
%========================================================================== H=1;L=1;nx=32;ny=32;h=1/nx;dt=0.002; gx=0;gy=-100;T0=0.0;kviscosity=0.01; alpha=0.01;Wallhot=100;Wallcold=0;nstep=2000; maxit=100;beta0=1.2; %========================================================================== u=zeros(nx+1,ny+2); v=zeros(nx+2,ny+1); p=zeros(nx+2,ny+2); tmp1=zeros(nx+2,ny+2); tmp2=zeros(nx+2,ny+2); c=zeros(nx+2,ny+2)+0.25; ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); T=zeros(nx+1,ny+1); TT=zeros(nx+1,ny+1); c(2,3:ny)=1/3; c(3:nx,2)=1/3; c(2,2)=1/2; c(nx+1,2)=1/2; c(nx+1,3:ny)=1/3; c(3:nx,ny+1)=1/3; c(2,ny+1)=1/2; c(nx+1,ny+1)=1/2; beta=1./((Tw+Te)/2+273); % 边界条件
数学模型：
控制方程： 使用 Boussinesq 近似来考虑温差所引起的自然对流， 在 N-S 方程中添加一个体积力项。
u 1 uu p 2u g T T0 t
(1)
其中， T0 是初始温度， 是运动粘度， 是流体的体胀系数。 能量守恒方程如下：
%==========================================================================
for i=1:nx+1 xh(i)=h*(i-1); end for j=1:ny+1 yh(j)=h*(j-1); end %========================================================================== for is=1:nstep for i=2:nx for j=2:ny+1 ut(i,j)=u(i,j)+dt*( -(0.25/h)*( (u(i,j)+u(i+1,j))^2-(u(i,j)+... u(i-1,j))^2+(v(i,j)+v(i+1,j))*(u(i,j+1)+u(i,j))-... (v(i,j-1)+v(i+1,j-1))*(u(i,j)+u(i,j-1)) )+... (kviscosity/h^2)*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-... 4*u(i,j))+beta*gx*(0.5*(T(i,j-1)+T(i,j))-T0) ); end end for i=2:nx+1 for j=2:ny vt(i,j)=v(i,j)+dt*( -(0.25/h)*( (u(i,j)+u(i,j+1))*(v(i,j)+... v(i+1,j))-(u(i-1,j)+u(i-1,j+1))*(v(i-1,j)+v(i,j))+... (v(i,j)+v(i,j+1))^2-(v(i,j)+v(i,j-1))^2 )+... (kviscosity/h^2)*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-... 4*v(i,j))+beta*gy*(0.5*(T(i-1,j)+T(i,j))-T0) ); end end for it=1:maxit for i=2:nx+1 for j=2:ny+1 p(i,j)=beta0*c(i,j)*(p(i+1,j)+p(i-1,j)+p(i,j+1)+p(i,j-1)-... (h/dt)*(ut(i,j)-ut(i-1,j)+vt(i,j)-vt(i,j-1)))+... (1-beta0)*p(i,j); end end end %========================================================================== u(2:nx,2:ny+1)=ut(2:nx,2:ny+1)-(dt/h)*(p(3:nx+1,2:ny+1)-p(2:nx,2:ny+1)); v(2:nx+1,2:ny)=vt(2:nx+1,2:ny)-(dt/h)*(p(2:nx+1,3:ny+1)-p(2:nx+1,2:ny)); %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); v(1,1:ny+1)=2*vw-v(2,1:ny+1); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1); % 求 修 正 速 度 % 求压力 % 临时速度v % 临时速度u