计算流体驱动方腔程序

驱动方腔（粘性不可压流），方腔示意图如左下：物理背景：如左图所示,装满黏性不可压流体的方腔上底面以速度U运动,其它壁面固定,其内部流体将作类似湍流的运动.流体运动复杂度将依赖于以下要素: (1)方腔形状 ；(2)U的大小；(3)流体的黏性,即雷诺数Re的大小。本题只讨论理想流体在二维正方形空腔中的流动，如左图，作无量纲化处理后,方腔长度和上底面运动速度均取为1,通过差分数值模拟计算出运动的流函数，涡量函数，画出Re=100和200时的流线图和涡量图。并标出窝心位置。采用流函数-涡方法：对于二维不可压缩的理想流体,流速 应满足如下方程组:：（1）连续方程： （1）(2) N-S方程： （2） （3）引入流函数：和：， （4）（4）代入（1）的流函数的方程：； （5）把（5）代入，的涡量方程： （6）以上方法称之为流函数涡方法。数值方法：将区域分为NN个方格单元，只有（N+1）（N+1）格网格点要计算，令h=1/N为空间步长，为时间步长。1：流函数与速度的边界条件：2涡量的边界条件：固壁上的涡量：计算：1，方腔内的流函数：采用中心五点差分格式对（5）进行离散并简化得：算出后，对（4）采用中心差分离散可以求得u,v.2,方腔内得涡量：采用迎风格式对涡量方程（6）进行离散：程序实现过程如下：clear allclcjishi=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;for t=1:timestep for i=2:Nx-1 for j=2:Ny-1 fluidfun(i,j)=(fluidfun(i+1,j)+fluidfun(i-1,j)+fluidfun(i,j+1)+fluidfun(i,j-1)+dh2*vortex(i,j)/4; endend for i=2:Nx-1 for j=2:Ny-1 velocityx(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); end end for i=2:Nx-1 for j=2:Ny-1 a=renold*(abs(velocityx(i,j)+abs(velocityy(i,j)/dh+4/dh2+renold/(dt); b=renold*abs(velocityx(i,j)/dh+1/dh2; c=renold*abs(velocityy(i,j)/dh+1/dh2; if (velocityx(i,j)=0)&(velocityy(i,j)=0) vortex(i,j)=(vortex(i+1,j)+vortex(i,j+1)/dh2+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)/dh2+b*vortex(i-1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt)/a; elseif (velocityx(i,j)=0) vortex(i,j)=(vortex(i-1,j)+vortex(i,j+1)/dh2+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)/dh2+b*vortex(i+1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt)/a; end end endendc1=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时驱动方腔流线图); hold onplot(d1*0.01,b1*0.01,ks)title(re=200时驱动方腔流线图);figure contour(x,y,vort