计算流体驱动方腔程序

驱动方腔（粘性不可压流），方腔示意图如左下：物理背景：如左图所示,装满黏性不可压流体的方腔上底面以速度U运动,其它壁面固定,其内部流体将作类似湍流的运动.流体运动复杂度将依赖于以下要素:(1)方腔形状；(2)U的大小；(3)流体的黏性,即雷诺数Re的大小。本题只讨论理想流体在二维正方形空腔中的流动，如左图，作无量纲化处理后,方腔长度和上底面运动速度均取为1,通过差分数值模拟计算出运动的流函数，涡量函数，画出Re=100和200时的流线图和涡量图。并标出窝心位置。采用流函数-涡方法：对于二维不可压缩的理想流体,流速),(uvu应满足如下方程组:：（1）连续方程：0u（1）(2)N-S方程：0)(2222xpyuxuuyuvxuutu（2）0)(2222ypyuxuuyuvxvutv（3）引入流函数：),,(tyx和),,(tyx：yu，xv，yuxv（4）（4）代入（1）的流函数的方程：）＋－（2222yx；（5）把（5）代入)2()3(xy，的涡量方程：)(Re12222yxyvxut（6）以上方法称之为流函数－涡方法。数值方法：将区域分为N×N个方格单元，只有（N+1）×（N+1）格网格点要计算，令h=1/N为空间步长，为时间步长。1：流函数与速度的边界条件：;0.........;0,0,1,1....,........2,1,0,,,x1,yx0,,,,,,ijiNiNiNiNiNijvuelsevuyifNjijhyih，2涡量的边界条件：固壁上的涡量：NjihhhHjNjNjjjiNi..........0,,2;2;2;/22,1,2,1,02,10,,计算：1，方腔内的流函数：采用中心五点差分格式对（5）进行离散并简化得：1...........1,;4/)(,21,1,,1,1,Njihjijijijijiji算出ji,后，对（4）采用中心差分离散可以求得u,v.2,方腔内得涡量：采用迎风格式对涡量方程（6）进行离散：);22(Re1212/1,1,11,21,11,1,112/1,12/1,1,12/11,2/11,,,hhhvhunjinjinjinjinjinjinjinjinjininjinjinjinji程序实现过程如下：clearallclcjishi=cputime;timestep=2000;[x,y]=meshgrid(0:1/100:1);Nx=length(x);Ny=length(y);dh=1/100;dt=0.6*dh;renold=200;fluidfun=zeros(Nx,Ny);vortex=zeros(Nx,Ny);velocityx=zeros(Nx,Ny);velocityy=zeros(Nx,Ny);%%%初始条件velocityx(:,Ny)=1;vortex(:,Ny)=-2/dh;%fluidfun(Ny,:)=0;%velocityy(Ny,:)=0;fort=1:timestepfori=2:Nx-1forj=2:Ny-1fluidfun(i,j)=(fluidfun(i+1,j)+fluidfun(i-1,j)+fluidfun(i,j+1)+fluidfun(i,j-1)+dh^2*vortex(i,j))/4;endendfori=2:Nx-1forj=2:Ny-1velocityx(i,j)=(fluidfun(i,j+1)-fluidfun(i,j-1))/(2*dh);velocityy(i,j)=(fluidfun(i-1,j)-fluidfun(i+1,j))/(2*dh);endendfori=2:Nx-1forj=2:Ny-1a=renold*(abs(velocityx(i,j))+abs(velocityy(i,j)))/dh+4/dh^2+renold/(dt);b=renold*abs(velocityx(i,j))/dh+1/dh^2;c=renold*abs(velocityy(i,j))/dh+1/dh^2;if(velocityx(i,j)=0)&&(velocityy(i,j)=0)vortex(i,j)=((vortex(i+1,j)+vortex(i,j+1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/(dt))/a;elseif(velocityx(i,j)=0)&&(velocityy(i,j)0)vortex(i,j)=((vortex(i+1,j)+vortex(i,j-1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt))/a;elseif(velocityx(i,j)0)&&(velocityy(i,j)=0)vortex(i,j)=((vortex(i-1,j)+vortex(i,j+1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/(dt))/a;elseif(velocityx(i,j)0)&&(velocityy(i,j)0)vortex(i,j)=((vortex(i-1,j)+vortex(i,j-1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt))/a;endendendendc1=max(abs(fluidfun));b1=find(c1==max(c1));%%%在第几列d1=find(abs(fluidfun(:,b1))==max(c1));%%%在第几行figurecontour(x,y,fluidfun,10)%text(0,0.2,['re=200时驱动方腔流线图']);holdonplot(d1*0.01,b1*0.01,'ks')title('re=200时驱动方腔流线图');figurecontour(x,y,vortex,10)%text(0,0.2,['re=200时驱动方腔等涡量图']);holdonplot(d1*0.01,b1*0.01,'ks')title('re=200时驱动方腔等涡量图');数值模拟结果：00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91re=200时驱动方腔等涡量图00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91re=200时驱动方腔流线图00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91re=100时驱动方腔流线图00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91re=100时驱动方腔等涡量图其中，小方块表示中心位置。