资源描述：
《计算流体驱动方腔程序.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、驱动方腔（粘性不可压流），方腔示意图如左下：物理背景：如左图所示,装满黏性不可压流体的方腔上底面以速度U运动,其它壁面固定,其内部流体将作类似湍流的运动.流体运动复杂度将依赖于以下要素:(1)方腔形状；(2)U的大小；(3)流体的黏性,即雷诺数Re的大小。本题只讨论理想流体在二维正方形空腔中的流动，如左图，作无量纲化处理后,方腔长度和上底面运动速度均取为1,通过差分数值模拟计算出运动的流函数，涡量函数，画出Re=100和200时的流线图和涡量图。并标出窝心位置。采用流函数-涡方法：对于二维不可压
2、缩的理想流体,流速应满足如下方程组:：（1）连续方程：（1）(2)N-S方程：（2）（3）引入流函数：和：，，（4）（4）代入（1）的流函数的方程：；（5）把（5）代入，的涡量方程：（6）以上方法称之为流函数－涡方法。数值方法：将区域分为N×N个方格单元，只有（N+1）×（N+1）格网格点要计算，令h=1/N为空间步长，为时间步长。1：流函数与速度的边界条件：2涡量的边界条件：固壁上的涡量：计算：1，方腔内的流函数：采用中心五点差分格式对（5）进行离散并简化得：算出后，对（4）采用中心差分离散可
3、以求得u,v.2,方腔内得涡量：采用迎风格式对涡量方程（6）进行离散：程序实现过程如下：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,N
4、y);%%%初始条件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-1ve
5、locityx(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;
6、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)=((vor
7、tex(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(v
8、elocityx(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));%%%在第几行
