【计算流体力学】附录E 有限体积法求解库塔流(E)

1、-E. -附录E二维不可压缩黏性流体Couette流动问题的有限体积算法与计算程序二维Couette流动是一个不可压缩黏性流动的典型流动，并有解析解，可以用来检验数值算法计算精度和可靠性。对它采用有限体积算法一阶迎风型离散格式进行数值求解。同时，为了初学者入门和练习方便，这里给出了用C语言和Fortran7语言编写的计算二维不可压缩黏性流体Couette流动问题计算程序，供大家学习参考。E-1利用有限体积算法一阶迎风型离散格式求解二维不可压缩黏性流体Couette流动问题二维不可压缩黏性流体Couette流动问题的提法二维不可压缩黏性流体Couette流动：有两个无限长平板，间距为L，两板之间

2、充满密度为1、静止的不可压缩黏性流体。无限长平板组成的通道两端压力相等，下板固定不动，上板以量纲为一的速度1自左向右平移运动（图E.1）。w=0=0图E.1二维不可压缩黏性流体Couette流动问题示意图基本方程组、初始条件和边界条件设流体是黏性流体。二维Couette流动问题在数学上可以由二维不可压缩黏性流动N-S方程组来表示，把它写成通用变量的微分方程组形式，有:(E.1)其中u为变量在水平x方向的流速，v为0在垂直y方向的流速，卩为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。初始条件：上板以量纲为一的速度1沿着上壁面方向自左向右运动。边界条件：流动速度u、v在上下壁面可采用

3、无滑移边界条件，在左右两端采用自由输出边界条件；压强p采用自由输出边界条件。计算网格划分和控制体单元与节点定义采用交错网格。图E.2和图E.3是计算网格、控制体单元和节点示意图。图E.2方腔流动计算网格、控制体单元和节点示意图图E.3计算采用的交错网格示意图节点P所在主控制体单元如图E.2中有阴影部分所示。在x方向与节点P相邻的节点为W和E，在y方向与节点P相邻的节点为S和N，主控制单元界面分别为w、s、e、n。压力p和速度u、v分别在三套不同网格中，如图E.3中有阴影部分所示。有限体积算法离散格式对方程（E.1）在图E.2所示节点P主控制单元内积分，有:JSdV(E.2)AVJ2(pu)dV

4、+(pvdV=J纠dV+卩绡dV+CxoyuxICx丿oy(oyAVAVAVAV由于不可压缩黏性流体Couette流动是二维问题，因此，控制体单元体积AV仅是面积，而它的边界是长度。设A=A=Ay,A=A=Ax，利用Gauss定理，可wesn将方程（E.2）改写成如下有限体积离散格式：r(pua)-(puA)+(pvA)-(pvA)Iewns(人QIQy丿n0八TITtQuA-1uA0X丿w、Qy丿0 x丿Ie00y丿lQx丿0、0y丿E.3)+SAxAy采用一阶向前差分近似，则有：s-=EP,Ax-NP,Ay0、l0X丿w0、0,uO,(F0,F0),主控制单元界面取值：ewew HYPER

5、LINK l bookmark72 =,=(E.8)wWeP则方程(E.7)离散为： HYPERLINK l bookmark74 (D+F)+D+(F-F)=(D+F)+D(E.9)wweewPwwWeE当流动为负向时，即u0,u0,(F0,F#includevstdlib.h#include#defineMX100最大网格数#defineMY20最大网格数#defineLambda0.002#defineRe100.00/Chorin压力迭代常数，直接影响收敛性雷诺数#definedt0.0005时间步长doubleuMX+lMY+2,vMX+2MY+l,pMX+2MY+2,unMX+lM

6、Y+2,vnMX+2MY+l;全局变量doublemax(doublea,doubleb)定义max函数doublec;if(ab)c=b;elsec=a;returnc;/*初始化输入:无；输出：u、v、p、dx、dy，初始速度、压强和空间步长。*/voidinit(doubleuMX+1MY+2,doublevMX+2MY+1,doublepMX+2MY+2,double&dx,double&dy)inti,j;dx=1.0/MX;dy=0.2/MY;for(i=0;i=MX;i+)for(j=0;j=MY+1;j+)uij=0.0;if(j=MY+1)uij=2.0;for(i=0;i=

7、MX+1;i+)for(j=0;j=MY;j+)vij=0.0;for(i=0;i=MX+1;i+)for(j=0;j=MY+1;j+)pij=1.0;/*Chorin迭代更新压强输入：u、v、p、dx、dy，n时刻速度、压强、空间步长；输出：p，n+1时刻压强。*/voidgetp(doubleuMX+1MY+2,doublevMX+2MY+1,doublepMX+2MY+2,doubledx,doubledy)doubledMX+2MY+2;inti,j;for(i=1;i=MX;i+)for(j=1;j=MY;j+)dij=(uij-ui-1j)/dx+(vij-vij-1)/dy;fo

8、r(i=1;i=MX;i+)for(j=1;j=MY;j+)pij=pij-Lambda*dij;for(i=1;i=MX;i+)pi0=pi1;piMY+1=piMY;for(j=1;j=MY;j+)p0j=p1j;pMX+1j=pMXj;/*一阶迎风格式输入：u、v、p、dx、dy,n时刻速度、n+1时刻压强，空间步长；输出：un、vn，n+1时刻速度。*/voidupwind(doubleuMX+1MY+2，doublevMX+2MY+1，doublepMX+2MY+2，doubleunMX+1MY+2，doublevnMX+2MY+1，doubledx，doubledy)inti，j;

9、doubleaw，ae，as，an，df，ap，miu;miu=1.0/Re;for(i=1;i=MX-1;i+)for(j=1;j=MY;j+)aw=miu+max(0.5*(ui-1j+uij)*dy，0.0);ae=miu+max(0.0，-0.5*(uij+ui+1j)*dy);as=miu+max(0.5*(vij-1+vi+1j-1)*dx，0.0);an=miu+max(0.0，-0.5*(vij+vi+1j)*dx);df=0.5*(ui+1j-ui-1j)*dy+0.5*(vij+vi+1j-vij-1-vi+1j-1)*dx;ap=aw+ae+as+an+df;unij=u

10、ij+dt/dx/dy*(-ap*uij+aw*ui-1j+ae*ui+1j+as*uij-1+an*uij+1)-dt*(pi+1j-pij)/dx;/边界条件for(i=1;i=MX-1;i+)uni0=-uni1;uniMY+1=2.0-uniMY;for(j=0;j=MY+1;j+)un0j=un1j;unMXj=unMX-1j;for(i=1;i=MX;i+)for(j=1;j=MY-1;j+)aw=miu+max(0.5*(ui-1j+ui-1j+1)*dy,0.0);ae=miu+max(0.0,-0.5*(uij+uij+1)*dy);as=miu+max(0.5*(vij-1

11、+vij)*dx,0.0);an=miu+max(0.0,-0.5*(vij+vij+1)*dx);df=0.5*(uij+uij+1-ui-1j-ui-1j+1)*dy+0.5*(vij+1-vij-1)*dx;ap=aw+ae+as+an+df;vnij=vij+dt/dx/dy*(-ap*vij+aw*vi-1j+ae*vi+1j+as*vij-1+an*vij+1)-dt*(pij+1-pij)/dy;/边界条件for(i=1;i=MX;i+)vni0=0.0;vniMY=0.0;for(j=0;j=MY;j+)vn0j=-vn1j;vnMX+1j=-vnMXj;/*格式输出，采用Te

12、cplot数据格式画图输入：u、v、p、dx、dy，最后结果；输出：无。*/voidoutput(doubleuMX+1MY+2,doublevMX+2MY+1,doublepMX+2MY+2,doubledx,doubledy)doubleuo,vo,po;inti,j;FILE*fp;fp=fopen(Result.plt,w);fprintf(fp,variables=x,y,u,v,pnzonei=%d,j=%d,f=pointrn,MX+1,MY+1);for(j=0;j=MY;j+)for(i=0;i1e-6&step1e6)err=0.0;getp(u,v,p,dx,dy);up

13、wind(u,v,p,un,vn,dx,dy);for(i=0;i=MX;i+)for(j=0;jerr)err=value;uij=unij;for(i=0;i=MX+1;i+)for(j=0;jerr)err=value;vij=vnij;-E.ll-printf(err=%len,err);output(u,v,p,dx,dy);/输出结果2.Fortran7语言源程序programfvm_upwind_MAC_couette!以一阶迎风型离散格式和Chorin压力迭代求解二维Couette流动问题(Fortran77吾言版本)!parameter(mx=101,my=21)implic

14、itdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my)doubleprecisionlambdaintegerstepre=100.!雷诺数dt=0.0005!时间步长lambda=0.002!Chorin压力迭代常数，直接影响收敛性dx=1.0d0/(mx-1)dy=0.2d0/(my-1)step=0err=100.callinitial(u,v,p)!初始化dowhile(errle-6.and.stepvle6)!err,收敛限制；step，迭代次数限制e

15、rr=0.0step=step+l!write(*,*)stepcallcalp(u,v,p,mx,my,dx,dy,lambda)!Chorin压力迭代callupwind(u,v,p,un,vn,mx,my,dx,dy,dt,re)!一阶迎风型离散格式求解动量方程doi=l,mxdoj=l,my+ltemp=abs(u(i,j)-un(i,j)/dt-E. -if(temperr)err=tempu(i,j)=un(i,j)enddoenddodoi=1,mx+1doj=1,mytemp=abs(v(i,j)-vn(i,j)/dtif(temperr)err=tempv(i,j)=vn(i

16、,j)enddoenddowrite(*,*)err=,errenddocalloutput(u,v,p,mx,my,dx,dy)!输出End!初始化!输入：无；!输出：u、v、p,初始速度、压强。!subroutineinitial(u,v,p)parameter(mx=101,my=21)doubleprecisionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)doi=1,mx+1doj=1,my+1p(i,j)=1.0enddoenddodoi=1,mxdoj=1,my+1u(i,j)=0.0if(j=my+1)u(i,j)=2.0enddoenddodoi=1,

17、mx+1doj=1,myv(i,j)=0.0enddoenddoendsubroutine!Chorin压力迭代!输入：u、v、p、mx、my、dx、dy、lambda,n时刻速、压强、其他各量;!输出：p,n+1时刻压强。!subroutinecalp(u,v,p,mx,my,dx,dy,lambda)implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),d(mx+1,my+1)doubleprecisionlambdadoi=2,mxdoj=2,myd(i,j)=(u(i,j)-u(i-1,j

18、)/dx+(v(i,j)-v(i,j-1)/dyenddoenddodoi=2,mxdoj=2,myp(i,j)=p(i,j)-lambda*d(i,j)enddoenddodoi=2,mxp(i,1)=p(i,2)p(i,my+1)=p(i,my)enddodoj=1,my+1p(1,j)=p(2,j)p(mx+1,j)=p(mx,j)enddoendsubroutine!一阶迎风型离散格式!输入：u、v、pn、mx、my、dx、dy、dt、re,n时刻速度、n+1时刻压强、其他各量;!输出：un,vn，n+1时刻速度。!subroutineupwind(u,v,pn,un,vn,mx,my

19、,dx,dy,dt,re)implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),pn(mx+1,my+1),un(mx,my+1),vn(mx+1,my)doubleprecisionmiumiu=1.0/redoi=2,mx-1doj=2,myaw=miu+max(0.5*(u(i-1,j)+u(i,j)*dy,0.0)ae=miu+max(0.0,-0.5*(u(i,j)+u(i+1,j)*dy)as=miu+max(0.5*(v(i,j-1)+v(i+1,j-1)*dx,0.0)an=miu+max(0.0,-0.5

20、*(v(i,j)+v(i+1,j)*dx)df=0.5*(u(i+1,j)-u(i-1,j)*dy+0.5*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)*dxap=aw+ae+as+an+dfun(i,j)=u(i,j)+dt/dx/dy*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1)&-dt*(pn(i+1,j)-pn(i,j)/dxenddoenddo!doi=2,mx-1un(i,1)=-un(i,2)un(i,my+1)=2.0-un(i,my)enddodoj=1,my+1un(1,j

21、)=un(2,j)un(mx,j)=un(mx-1,j)enddo!u边界条件!doi=2,mxdoj=2,my-1aw=miu+max(0.5*(u(i-1,j)+u(i-1,j+1)*dy,0.0)ae=miu+max(0.0,-0.5*(u(i,j)+u(i,j+1)*dy)as=miu+max(0.5*(v(i,j-1)+v(i,j)*dx,0.0)an=miu+max(0.0,-0.5*(v(i,j)+v(i,j+1)*dx)df=0.5*(u(i,j)+u(i,j+1)-u(i-1,j)-u(i-1,j+1)*dy+0.5*(v(i,j+1)-v(i,j-1)*dxap=aw+ae+as+an+dfvn(i,j)=v(i,j)+dt/dx/dy*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1)&-dt*(pn(i,j+1)-pn(i,j)/dyenddoenddo!d