附录F-有限体积法计算方腔流(F)

《附录F-有限体积法计算方腔流(F)》由会员分享，可在线阅读，更多相关《附录F-有限体积法计算方腔流(F)（23页珍藏版）》请在装配图网上搜索。

1、精选优质文档-倾情为你奉上附录F 二维不可压缩黏性流体方腔流动问题的有限体积算法与计算程序二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。F-1利用有限体积算法三阶迎风型离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（cavity

2、 flow图F.1二维不可压缩黏性方腔流问题示意图）：有一正方形腔室，其量纲为一的宽度为，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度沿着上壁面方向自左向右运动（图F.1）。2. 基本方程组、初始条件和边界条件图F.1 二维不可压缩黏性方腔流动问题示意图设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S方程组来表示，把它写成通用变量的微分方程组形式，有： （F.1）其中为变量在水平方向的流速，为在垂直方向的流速，为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。初始条件：方腔上壁面以

3、量纲为一的速度沿着上壁面方向自左向右运动。边界条件： 流动速度均可采用无滑移边界条件，压力采用自由输出边界条件。3.计算网格划分和控制体单元与节点定义采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。图F.2方腔流动计算网格、控制体单元和节点示意图图F.3计算采用的交错网格示意图节点所在主控制单元如图F.2中有阴影部分所示。在方向与节点相邻的节点为和，在方向与节点相邻的节点为和，主控制单元界面分别为。压力和速度分别在三套不同网格中如图F.3中有阴影部分所示。4.有限体积算法三阶迎风型离散格式对方程（F.1）在图F.2所示节点所在控制体单元内积分，有：（F.2）由于二维不可压缩

4、黏性流体方腔流动是二维问题，因此控制体单元体积仅是面积，而它的边界是长度。设 ，利用定理，可将方程（F.2）改写成如下有限体积算法离散格式： （F.3）对上式中采用一阶向前差分近似，则有： （F.4）同时记： （F.5） （F.6）则可由式（F.2）写成： （F.7）式中都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量，就可以求出节点上未知量。图F.4三阶迎风型离散格式示意图为了便于讨论，现对一维对流扩散方程的三阶迎风型离散格式进行分析：在三阶迎风型离散格式中，计算主控制单元界面上流动量需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界

5、面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。当时，通过、和三个节点值拟合曲线来计算主控制单元左侧界面参数。通过节点、和三个节点值拟合曲线来计算主控制单元右侧界面参数。当，则分别通过节点、和、三个节点值计算主控制单元左、右两侧界面参数和。根据上述计算原则，可以得到界面参数计算公式如下：当时，界面参数计算公式为： （F.8a）当时，界面参数计算公式为： （F.8b）对于一维无源项一维对流扩散方程三阶迎风型离散格式：当时，三阶迎风型离散格式为： （F.9）其中 （F.9）同理，若，三阶迎风型离散格式为： （F.10）其中 （F.10a）将两种流动方向离散方程（F.9）和（F.10）

6、合并后，可得到统一的一维对流扩散方程三阶迎风型离散格式： （F.11）其中 （F.11a）式中 （F.11b）同理，可以得到带有源项的二维对流扩散方程三阶迎风型离散格式为： （F.12）其中为有限体积算法中源项平均值。式中各个系数为： （F.12a）式中 （F.12b）源项为： （F.13）若把表示时刻动量，表示时刻动量，则可以得到源项离散格式为： （F.14）最后，得到有限体积算法二维对流扩散方程三阶迎风型离散格式： （F.15）式中系数为一阶迎风格式中各对应系数。5.计算结果分析利用三阶迎风型离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。图F.5是不同雷诺数条件

7、下采用三阶迎风型离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。这表明有限体积算法三阶迎风型离散格式具有相当高的计算精度。=1 000=100 =10 000=5 000 图F.5不同雷诺数条件下采用三阶迎风型离散格式计算二维不可压缩黏性方腔流动的计算结果 从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。这是因为，方腔上壁面以量纲为一的速度沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件

8、，这是一个带奇性的方腔流动。计算结果表明，方腔流动和雷诺数有关，随着雷诺数增加，计算精度在降低。当雷诺数较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数的增加，二次小涡变得越来越模糊。由于三阶迎风型离散格式计算精度较高，因此三阶迎风型离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。F-2 二维不可压缩黏性方腔流动问题数值计算源程序1. 语言源程序/ fvm_quick_cavity.cpp/*-!利用有限体积算法三阶迎风型离散格式和!人工压缩算法求解方腔流动问题（语言版本）-*/#include #include #include #define MX 100#defin

9、e MY 100#define Re 100.00#define dt 0.0005#define c2 2.25Double uMX+1MY+2,vMX+2MY+1,pMX+2MY+2,unMX+1MY+2,vnMX+2MY+1,pnMX+2MY+2;/全局变量/*- 定义两个要用到的函数-*/double max(double a,double b)if(a=0)return 1.0;else return 0.0;/*- 初始化 入口:无； 出口：u、v、p、dx、dy，初始速度、压强和空间步长。-*/void init(double uMX+1MY+2,double vMX+2MY+1

10、,double pMX+2MY+2,double& dx,double& dy) int i,j;dx=1.0/MX;dy=1.0/MY; for(i=0;i=MX;i+)for(j=0;j=MY+1;j+)uij=0.0;if(j=MY+1) uij=4.0/3.0;if(j=MY) uij=2.0/3.0;for(i=0;i=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;/*- 一阶迎风型离散格式 二维的三阶迎风型离散格式为9点格式，因此有两层边界网格需要处理，本程序采用一阶迎风型

11、离散格式处理内层，用物理边界条件处理外层。 入口:u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号； 出口：un，新的x方向速度。-*/void upwind_u(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double unMX+1MY+2,double dx,double dy,int i,int j) double aw,ae,as,an,df,ap,miu;miu=1.0/Re;aw=miu+max(0.5*(ui-1j+uij)*dy,0.0);ae=miu+max(0.0,-0.5*(uij+ui+1j

12、)*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=uij + dt/dx/dy*(-ap*uij+aw*ui-1j+ae*ui+1j+as*uij-1+an*uij+1) - dt*(pi+1j-pij)/dx;/*- 入口: u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号； 出口：vn，新的y方向速

13、度。-*/void upwind_v(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double vnMX+2MY+1,double dx,double dy,int i,int j)double aw,ae,as,an,df,ap,miu;miu=1.0/Re; 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+vij)*dx,0.0);an=miu+max(0.0,-0.5*(vij+vij+

14、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;/*- 三阶迎风型离散格式 入口:u、v、p、dx、dy，当前速度、压强，空间步长； 出口：un、vn，新的速度。-*/void quick(double uMX+1MY+2,double vMX+2MY+1,double pMX+2MY+2,double un

15、MX+1MY+2,double vnMX+2MY+1,double dx,double dy)double miu,fw,fe,fs,fn,df,aw,ae,as,an,aww,aee,ass,ann,ap;int i,j; miu=1.0/Re; for(i=2;i=MX-2;i+)for(j=2;j=MY-1;j+) fw=0.5*(ui-1j+uij)*dy; fe=0.5*(uij+ui+1j)*dy; fs=0.5*(vij-1+vi+1j-1)*dx; fn=0.5*(vij+vi+1j)*dx; df=fe-fw+fn-fs; aw=miu+0.750*alfa(fw)*fw+0

16、.125*alfa(fe)*fe+0.375*(1.0-alfa(fw)*fw; aww=-0.125*alfa(fw)*fw; ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe)*fe-0.125*(1.0-alfa(fw)*fw; aee=0.125*(1.0-alfa(fe)*fe; as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs)*fs; ass=-0.125*alfa(fs)*fs; an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn)

17、*fn-0.125*(1.0-alfa(fs)*fs; ann=0.125*(1.0-alfa(fn)*fn; ap=aw+ae+as+an+aww+aee+ass+ann+df;/aw、ae、as、an.均为有限体积算法中各项系数，详见前文三阶迎风型QUICK离散格式。 unij=uij+ dt/dx/dy*(-ap*uij+aw*ui-1j+ae*ui+1j+as*uij-1+an*uij+1+aww*ui-2j+aee*ui+2j+ass*uij-2+ann*uij+2)-dt*(pi+1j-pij)/dx;/- j=1; for(i=2;i=MX-2;i+) upwind_u(u,v,

18、p,un,dx,dy,i,j); j=MY; for(i=2;i=MX-2;i+) upwind_u(u,v,p,un,dx,dy,i,j); i=1; for(j=1;j=MY;j+) upwind_u(u,v,p,un,dx,dy,i,j); i=MX-1; for(j=1;j=MY;j+) upwind_u(u,v,p,un,dx,dy,i,j);/内层边界由一阶迎风型离散格式得到-/- for(i=1;i=MX-1;i+) uni0=-uni1; uniMY+1=2.0-uniMY; for(j=0;j=MY+1;j+) un0j=0.0; unMXj=0.0; /外层边界条件按物理边

19、界条件给出/- for(i=2;i=MX-1;i+) for(j=2;j=MY-2;j+) fw=0.5*(ui-1j+ui-1j+1)*dy; fe=0.5*(uij+uij+1)*dy; fs=0.5*(vij-1+vij)*dx; fn=0.5*(vij+vij+1)*dx; df=fe-fw+fn-fs; aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw)*fw; aww=-0.125*alfa(fw)*fw; ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe)*fe-0.1

20、25*(1.0-alfa(fw)*fw; aee=0.125*(1.0-alfa(fe)*fe; as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs)*fs; ass=-0.125*alfa(fs)*fs; an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn)*fn-0.125*(1.0-alfa(fs)*fs; ann=0.125*(1.0-alfa(fn)*fn; ap=aw+ae+as+an+aww+aee+ass+ann+df; vnij=vij + dt/dx/dy*(-ap*

21、vij+aw*vi-1j+ae*vi+1j+as*vij-1+an*vij+1+aww*vi-2j+aee*vi+2j+ass*vij-2+ann*vij+2) - dt*(pij+1-pij)/dy; /- j=1; for(i=2;i=MX-1;i+) upwind_v(u,v,p,vn,dx,dy,i,j); j=MY-1; for(i=2;i=MX-1;i+) upwind_v(u,v,p,vn,dx,dy,i,j); i=1; for(j=1;j=MY-1;j+) upwind_v(u,v,p,vn,dx,dy,i,j); i=MX; for(j=1;j=MY-1;j+) upwin

22、d_v(u,v,p,vn,dx,dy,i,j);/- for(i=1;i=MX;i+) vni0=0.0; vniMY=0.0; for(j=0;j=MY;j+) vn0j=-vn1j; vnMX+1j=-vnMXj; /-/*- 更新压强 入口: un、vn、p、dx、dy，新的速度，当前压强，空间步长； 出口: pn，新的压强。-*/void getp(double unMX+1MY+2,double vnMX+2MY+1,double pMX+2MY+2,double pnMX+2MY+2,double dx,double dy)int i,j;for(i=1;i=MX;i+)for(j

23、=1;j=MY;j+)pnij=pij-dt*c2/dx*(unij-uni-1j+vnij-vnij-1);for(i=1;i=MX;i+)pni0=pni1; pniMY+1=pniMY;for(j=0;j1e-4 & step1e6) /err1e-5，定常解判据；step，限制迭代次数printf(step=%dn,step);step+; err=0.0;quick(u,v,p,un,vn,dx,dy);getp(un,vn,p,pn,dx,dy);/-for(i=0;i=MX;i+)for(j=0;jerr) err=value;uij=unij;for(i=0;i=MX+1;i+

24、)for(j=0;jerr) err=value;vij=vnij;for(i=0;i=MX+1;i+)for(j=0;jerr)err=value;pij=pnij;printf(err=%len,err);/-/输出结果，用数据格式画图FILE *fp;fp=fopen(Result.plt,w);fprintf(fp,variables=x,y,u,v,pn zone i=%d,j=%d,f=pointn,MX+1,MY+1);for(i=0;i=MX;i+)for(j=0;j=MY;j+)uo=0.5*(uij+uij+1); vo=0.5*(vij+vi+1j); po=0.25*(

25、pij+pi+1j+pij+1+pi+1j+1);fprintf(fp,%20.10e %20.10e %20.10e %20.10e %20.10en,i*dx,j*dy,uo,vo,po);fclose(fp);-2. 语言源程序! !利用有限体积算法三阶迎风型离散格式和!人工压缩算法求解方腔流动问题（语言版本）! program QUICK_cavity parameter(mx=101,my=101) implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1) dimension un(

26、mx,my+1),vn(mx+1,my),pn(mx+1,my+1) common /ini/u,v,p c2=2.25 re=100.0 dt=0.0005 dx=1.0/float(mx-1) dy=1.0/float(my-1)!-! u、v、p为t时刻值，un、vn、pn为t+1时刻值，! mx、my为最大网格数，c2为虚拟压缩系数的平方，re为雷诺数。!- num=0 err=100.00!nun，计数器；err，判断人工压缩法求解收敛的标准 call initial!调入初始条件，以下为人工压缩算法求解 do while(err.gt.1e-4.and.num.lt.1e6) er

27、r=0.0 call quick(u,v,p,un,vn,mx,my,dx,dy,dt,re)!QUICK离散格式求解动量方程，得到un、vn call calp(p,un,vn,pn,mx,my,dx,dy,dt,c2)!求压强pn call check(u,v,p,un,vn,pn,mx,my,dx,dy,dt,c2,err)!校验流场信息，判断是否收敛，同时更新u、v、p write(*,*) error=,err num=num+1 write(*,*) num!屏幕跟踪输出 enddo call output(u,v,p,mx,my,dx,dy)!输出结果文件 End! subrou

28、tine initial!初始化流场 parameter(mx=101,my=101) double precision u(mx,my+1),v(mx+1,my),p(mx+1,my+1) common /ini/u,v,p do i=1,mx+1 do j=1,my+1 p(i,j)=1.0 enddo enddo do i=1,mx do j=1,my+1 u(i,j)=0.0 if(j.eq.my+1) u(i,j)=4.0/3.0 if(j.eq.my) u(i,j)=2.0/3.0 enddo enddo do i=1,mx+1 do j=1,my v(i,j)=0.0 enddo

29、 enddo endsubroutine! subroutine quick(u,v,p,un,vn,mx,my,dx,dy,dt,re)!以QUICK格式离散动量方程 implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my) double precision miu miu=1.0/re!以下求解x方向速度un- do i=3,mx-2 do j=3,my-1 fw=0.5*(u(i-1,j)+u(i,j)*dy fe=0.5*(u(i,j)+

30、u(i+1,j)*dy fs=0.5*(v(i,j-1)+v(i+1,j-1)*dx fn=0.5*(v(i,j)+v(i+1,j)*dx df=fe-fw+fn-fs aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw)*fw aww=-0.125*alfa(fw)*fw ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe)*fe-0.125*(1.0-alfa(fw)*fw aee=0.125*(1.0-alfa(fe)*fe as=miu+0.750*alfa(fs)*fs+0.1

31、25*alfa(fn)*fn+0.375*(1.0-alfa(fs)*fs ass=-0.125*alfa(fs)*fs an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn)*fn-0.125*(1.0-alfa(fs)*fs ann=0.125*(1.0-alfa(fn)*fn ap=aw+ae+as+an+aww+aee+ass+ann+df!aw、ae、as、an.均为有限体积算法中各项系数，详见前文三阶迎风型QUICK离散格式 un(i,j)=u(i,j)+dt/dx/dy*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j) &+

32、as*u(i,j-1)+an*u(i,j+1)+aww*u(i-2,j)+aee*u(i+2,j) &+ass*u(i,j-2)+ann*u(i,j+2)-dt*(p(i+1,j)-p(i,j)/dx enddo enddo!- j=2 do i=3,mx-2 call upbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j) enddo j=my do i=3,mx-2 call upbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j) enddo i=2 do j=2,my call upbound_u(u,v,p,un,mx,my,dx

33、,dy,dt,re,i,j) enddo i=mx-1 do j=2,my call upbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j) enddo!内层边界由一阶迎风型离散格式得到-!- do i=2,mx-1 un(i,1)=-un(i,2) un(i,my+1)=2.0-un(i,my) enddo do j=1,my+1 un(1,j)=0.0 un(mx,j)=0.0 enddo!外层边界条件按物理边界条件给出!-!以下求解y方向速度vn- do i=3,mx-1 do j=3,my-2 fw=0.5*(u(i-1,j)+u(i-1,j+1)*dy f

34、e=0.5*(u(i,j)+u(i,j+1)*dy fs=0.5*(v(i,j-1)+v(i,j)*dx fn=0.5*(v(i,j)+v(i,j+1)*dx df=fe-fw+fn-fs aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw)*fw aww=-0.125*alfa(fw)*fw ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe)*fe-0.125*(1.0-alfa(fw)*fw aee=0.125*(1.0-alfa(fe)*fe as=miu+0.750*alfa(

35、fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs)*fs ass=-0.125*alfa(fs)*fs an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn)*fn-0.125*(1.0-alfa(fs)*fs ann=0.125*(1.0-alfa(fn)*fn ap=aw+ae+as+an+aww+aee+ass+ann+df vn(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)+aww*v(i-2,j

36、)+aee*v(i+2,j) &+ass*v(i,j-2)+ann*v(i,j+2)-dt*(p(i,j+1)-p(i,j)/dy enddo enddo!- j=2 do i=3,mx-1 call upbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j) enddo j=my-1 do i=3,mx-1 call upbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j) enddo i=2 do j=2,my-1 call upbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j) enddo i=mx do j=

37、2,my-1 call upbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j) enddo!- do i=2,mx vn(i,1)=0.0 vn(i,my)=0.0 enddo do j=1,my vn(1,j)=-vn(2,j) vn(mx+1,j)=-vn(mx,j) enddo!- Endsubroutine! function alfa(x)!函数,正1负0 double precision alfa, x if(x.gt.0.d0) then alfa=1.0 else alfa=0.0 endif end! subroutine upbound_u(u,

38、v,p,un,mx,my,dx,dy,dt,re,i,j)!以一阶迎风型离散格式得到内层边界un的值 implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1) double precision miu miu=1.0/re aw=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)

39、an=miu+max(0.0,-0.5*(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)*dx ap=aw+ae+as+an+df un(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*(p(i+1,j)-p(i,j)/dx Endsubroutine! subroutine upbound_v(u,v,p,vn,mx,my,dx,dy,dt

40、,re,i,j)!以一阶迎风型离散格式得到内层边界vn值 implicit double precision(a-h,o-z) dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),vn(mx+1,my) double precision miu miu=1.0/re aw=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

41、(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)*dx ap=aw+ae+as+an+df vn(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*(p(i,j+1)-p(i,j)/dy Endsubroutine! subroutine calp(p,un,vn,pn,mx,my,dx,dy,dt,c2)!根据un、vn求解压强pn implicit double precision(a-h,o-z) dimension p(mx+1,my+1),un(mx,my+1),vn(mx+1,my),pn(mx+1,my+1) do i=2,mx do j=2,my pn(i,j)=p(i,j)-dt*c2/dx*(un(i,j)-un(i-1,j)+vn(i,j)-vn(i,j-1) enddo enddo