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

自然对流传热问题描述：一个二维矩形腔体，高H1，宽L1，物理模型如图1。绝热壁面绝热壁面图1物理模型上下壁面为绝热壁面，左壁面为热壁面100T，右壁面为冷壁面0T。腔体内为初始温度00T的空气，计算由此引起的自然对流传热。数学模型：控制方程：使用Boussinesq近似来考虑温差所引起的自然对流，在N-S方程中添加一个体积力项。201pgTTtuuuu(1)其中，0T是初始温度，是运动粘度，是流体的体胀系数。能量守恒方程如下：2TTTtu(2)其中，是热扩散系数。边界条件：由无滑移边界条件得，四个壁面上的速度均为零，即：0unusuwue。在热壁面上100T，在冷壁面上0T，在上下绝热壁面上处0Ty。数值处理：区域离散化如图2所示。对于动量守恒方程，在不考虑压力的情况下先计算出一个临时速度20nnnnuugTTdtuuu(3)用SOR（超松弛迭代）法求压力场111,1,1,,1,1,1,,1,,1ijijijijijijijijijijpcpppphuuvvpdt(4)根据求得的临时速度，在考虑压力的情况下，得修正速度冷壁面T=0热壁面T=1001,,1,,1nijijijijuuppdtdx(5)方程(3)对应的离散方程如下：22,,,1,,1,,1,,1,,11,1,,1+1,,1,212222222tttttttttttijijijijijijijijijijtttttttijijijijijijijuuuuuuvvuudthvvuuuuuh,1,,12,1,022tttijijijijijuuuhTTgxT(6),,,,1,1,1,1,1,1,22,,1,,1+1,,1,212222222tttttttttttijijijijijijijijijijtttttttijijijijijijijvvuuvvuuvvdthvvvvvvvh,1,,12,,1022tttijijijijijvvvhTTgxT(7)图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);ut=zeros(nx+1,ny+2);vt=zeros(nx+2,ny+1);tmp1=zeros(nx+2,ny+2);uu=zeros(nx+1,ny+1);vv=zeros(nx+1,ny+1);tmp2=zeros(nx+2,ny+2);T=zeros(nx+1,ny+1);TT=zeros(nx+1,ny+1);c=zeros(nx+2,ny+2)+0.25;%==========================================================================c(2,3:ny)=1/3;c(nx+1,3:ny)=1/3;c(3:nx,2)=1/3;c(3:nx,ny+1)=1/3;c(2,2)=1/2;c(2,ny+1)=1/2;c(nx+1,2)=1/2;c(nx+1,ny+1)=1/2;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;beta=1./((Tw+Te)/2+273);%==========================================================================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);%==========================================================================fori=1:nx+1xh(i)=h*(i-1);endforj=1:ny+1yh(j)=h*(j-1);end%==========================================================================foris=1:nstepfori=2:nxforj=2:ny+1%临时速度uut(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));endendfori=2:nx+1forj=2:ny%临时速度vvt(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));endendforit=1:maxit%求压力fori=2:nx+1forj=2:ny+1p(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);endendend%==========================================================================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);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));%==========================================================================fori=2:nxforj=2:nyT(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%求温度场endfori=2:nxT(i,1)=T(i,2);T(i,ny+1)=T(i,ny);endTT=T;%==========================================================================end%==========================================================================fori=1:nx+1forj=1:ny+1velocity(i,j)=sqrt(uu(i,j)*uu(i,j)+vv(i,j)*vv(i,j));endend%==========================================================================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'),axisequal,axis([0L0H]);clf(figure)contour(xh,yh,flipud(rot90(T))),title('Temperaturemap'),axisequal,axis([0L0H]);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)