程序

附录F 二维不可压缩黏性流体方腔流动问题的

有限体积算法与计算程序

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

F-1利用有限体积算法三阶迎风型QUICK离散格式求解

二维不可压缩黏性流体方腔流动问题

1.二维不可压缩黏性流体方腔流动问题

二维不可压缩黏性流体方腔流动（cavity flow）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度

=1.0它周围壁面（左右壁面和为p=1.0、

底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。

2. 基本方程组、初始条件和边界条件

设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动

N - S

图F.1二维不可压缩黏性

方腔流问题示意图

） u v x y x x

其中u为变量 在水平x方向的流速，v为 在垂直y方向的流速， 为黏度，S 为源项。源项中不仅包含压力梯度项，也包含时间导数项。

初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。

程序

边界条件：

v均可采用无滑移边界条件，压力p采用自由输出边界条件。 流动速度u、

3.计算网格划分和控制体单元与节点定义

采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。

图F.2方腔流动计算网格、 控制体单元和节点示意图

图F.3计算采用的交错网格示意图

节点P所在主控制单元如图F.2中有阴影部分所示。在x方向与节点P相邻的节点为W和E，在y方向与节点P相邻的节点为S和N，主控制单元界面分

sen。压力p和速度u、v分别在三套不同网格中如图F.3中有阴影部分别为w、、、

所示。

4.有限体积算法三阶迎风型QUICK离散格式

对方程（F.1）在图F.2所示节点P所在控制体单元内积分，有：

u dV v dV dV dV S dV（F.2） x y x x y y V V V V V由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积 V仅是面积，而它的边界是长度。设 Aw Ae y,As An x，利用Gauss定理，可将方程（F.2）改写成如下有限体积算法离散格式：

程序

u A e u A w v A n v A s

（F.3）

A A A A S x y

x e x w y n y s

对上式中 、 、 采用一阶向前差分近似，则有： 、

x x y e w n y s

P W E P

,

x x x e x w

N P P S

,

y y y n y s

Fe u eAe,Fw u wAwFn v nAn,Fs v sAs

（F.4）

同时记：

（F.5）

eAe A

,Dw ww, xPE xPW

（F.6）

nAn sAs

Dn ,Ds

yPN yPS

De

则可由式（F.2）写成：

Fe e Fw w Fn n Fs s De E P Dw P W

Dn N P Ds P S S x y

（F.7）

式中 P、 E、 W、 N、 S、De、Dw、Dn、Ds都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量Fe e、Fw w、Fn n、Fs s，就可以求出节点上未知量 P、 E、 W、 N、 S。

图F.4三阶迎风型QUICK离散格式示意图

程序

为了便于讨论，现对一维对流扩散方程的三阶迎风型QUICK离散格式进行分析：在三阶迎风型QUICK离散格式中，计算主控制单元界面上流动量 需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。

当ue 0,uw 0时，通过WW、W和P三个节点值拟合曲线来计算主控制单元左侧界面参数 w。通过节点W、P和E三个节点值拟合曲线来计算主控制单元右侧界面参数 e。当ue 0,uw 0，则分别通过节点W、P、E和P、E、EE三个节点值计算主控制单元左、右两侧界面参数 w和 e。根据上述计算原则，可以得到界面参数 w计算公式如下：

当uw 0时，界面参数 w计算公式为：

ww wW wP wWW （F.8a）

当ue 0时，界面参数 w计算公式为：

636

e P E W （F.8b）

888

683868

对于一维无源项一维对流扩散方程三阶迎风型QUICK离散格式： 当ue 0,uw 0时，三阶迎风型QUICK离散格式为：

aP P aW W aE E aWW WW （F.9） 其中

31

aW Dw Fw Fe

883

aE De Fe

8 （F.9a） 1

aWW Fw

8

aP aW aE aWW Fe Fw

同理，若ue 0,uw 0，三阶迎风型QUICK离散格式为：

aP P aW W aE E aEE EE （F.10）

程序

其中

3

aW Dw Fw,

861

aE De Fe Fw,

88 （F.10a） 1aWW Fe

8

aP aW aE aEE Fe Fw

将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型QUICK离散格式：

（F.11） aP P aW W aE E aEE EE aEE EE

其中

613

aW Dw wFw eFe 1 w Fw

8881

aWW wFw

8361

aE De eFe 1 e Fe 1 w Fw （F.11a）

8881

aEE 1 e Fe

8

aP aW aE aWW aEE Fe Fw

式中

当Fw 0时， w 1；

当Fw 0时， w 0；

（F.11b）

当Fe 0时， e 1；当Fe 0时， e 0。

同理，可以得到带有源项的二维对流扩散方程三阶迎风型QUICK离散格式为：

aP P aW W aE E aEE EE aEE EE

aS S aN N aSS SS aNN NN S V

其中 为有限体积算法中源项平均值。式中各个系数为：

（F.12）

程序

613

aW Dw wFw eFe 1 w Fw

8881

aWW wFw

8361

aE De eFe 1 e Fe 1 w Fw

888

1

1 e Fe8

613

aS Ds sFs nFn 1 s Fs

8881

aSS sFs

8361

aN Dn nFn 1 n Fn 1 s Fs

8881

aNN 1 n Fn

8aEE

aP aW aE aWW aEE aS aN aSS aNN Fe Fw Fn Fs （F.12a）

式中

当Fk 0时， k 1；

当Fk 0时， k 0。 （F.12b）

k w,e,s,n

源项 为：

若把 u 表示tn时刻动量， u 格式为：

n

u p

（F.13） t x

n 1

表示tn 1时刻动量，则可以得到源项 离散

V

dV

u P

n 1

u P

n

t

x y pe pw y （F.14）

最后，得到有限体积算法二维对流扩散方程三阶迎风型QUICK离散格式：

nnnnn

aPuP aWuW aEuE aSuS aNuN

u P

n 1

u P

n

t

x y p p

n

e

nw

y

（F.15）

式中系数ak为一阶迎风格式中各对应系数。

程序

5.计算结果分析

利用三阶迎风型QUICK离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。图F.5是不同雷诺数Re条件下采用三阶迎风型QUICK离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。

计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。这表明有限体积算法三阶迎风型QUICK离散格式具有相当高的计算精度。

Re=100 Re=1 000

图F.5不同雷诺数Re条件下采用三阶迎风型QUICK

离散格式计算二维不可压缩黏性方腔流动的计算结果

Re=5 000 Re=10 000

程序

从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。这是因为，方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件，这是一个带奇性的方腔流动。

计算结果表明，方腔流动和雷诺数有关，随着雷诺数Re增加，计算精度在降低。当雷诺数Re较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数Re的增加，二次小涡变得越来越模糊。由于三阶迎风型QUICK离散格式计算精度较高，因此三阶迎风型QUICK离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。

F-2 二维不可压缩黏性方腔流动问题数值计算源程序

1. C语言源程序

// fvm_quick_cavity.cpp

/*---------------------------------------------------------------------------- !利用有限体积算法三阶迎风型QUICK离散格式和

!人工压缩算法求解方腔流动问题（C语言版本）

-------------------------------------------------------------------------------*/ #include <stdio.h> #include <stdlib.h> #include <math.h> #define MX 100 #define MY 100 #define Re 100.00 #define dt 0.0005 #define c2 2.25

Double u[MX+1][MY+2],v[MX+2][MY+1],p[MX+2][MY+2],

un[MX+1][MY+2],vn[MX+2][MY+1],pn[MX+2][MY+2];

//全局变量

程序

/*------------------------------------------------------- 定义两个要用到的函数

--------------------------------------------------------*/ double max(double a,double b) { if(a<b) return b; else return a; }

double alfa(double x) { if(x>=0) return 1.0; else return 0.0; }

/*------------------------------------------------------------------------ 初始化 入口:无；

出口：u、v、p、dx、dy，初始速度、压强和空间步长。 ---------------------------------------------------------------------------*/

void init(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+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++) { u[i][j]=0.0; if(j==MY+1) u[i][j]=4.0/3.0; if(j==MY) u[i][j]=2.0/3.0; } } for(i=0;i<=MX+1;i++) for(j=0;j<=MY;j++) v[i][j]=0.0; for(i=0;i<=MX+1;i++)

程序

for(j=0;j<=MY+1;j++) p[i][j]=1.0; }

/*----------------------------------------------------------------------------------------------- 一阶迎风型离散格式

二维的三阶迎风型QUICK离散格式为9点格式，因此有两层边界网格需要

处理，本程序采用一阶迎风型离散格式处理内层，用物理边界条件处理外层。 入口:u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号； 出口：un，新的x方向速度。

-----------------------------------------------------------------------------------------------*/

void upwind_u(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+2], double un[MX+1][MY+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*(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*(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; }

/*------------------------------------------------------------------------------------------------ 入口: u、v、p、dx、dy、i、j，当前速度、压强，空间步长和网格节点编号； 出口：vn，新的y方向速度。

-----------------------------------------------------------------------------------------------------*/

void upwind_v(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+2], double vn[MX+2][MY+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*(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])*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;

程序

}

/*---------------------------------------------------------------------- 三阶迎风型QUICK离散格式

入口:u、v、p、dx、dy，当前速度、压强，空间步长； 出口：un、vn，新的速度。

-----------------------------------------------------------------------*/

void quick(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+2], double un[MX+1][MY+2],double vn[MX+2][MY+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*(u[i-1][j]+u[i][j])*dy; fe=0.5*(u[i][j]+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.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; //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]+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; } }

//------------------------------------------------------------------------------ j=1; for(i=2;i<=MX-2;i++) upwind_u(u,v,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++) { un[i][0]=-un[i][1]; un[i][MY+1]=2.0-un[i][MY]; } for(j=0;j<=MY+1;j++) { un[0][j]=0.0; un[MX][j]=0.0; }

//外层边界条件按物理边界条件给出

//------------------------------------------------------------------------------- for(i=2;i<=MX-1;i++) { for(j=2;j<=MY-2;j++) { fw=0.5*(u[i-1][j]+u[i-1][j+1])*dy; fe=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(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]+aee*v[i+2][j]+ass*v[i][j-2]+ann*v[i][j+2])

- dt*(p[i][j+1]-p[i][j])/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++) upwind_v(u,v,p,vn,dx,dy,i,j);

//---------------------------------------------------------------------------- for(i=1;i<=MX;i++) { vn[i][0]=0.0; vn[i][MY]=0.0; } for(j=0;j<=MY;j++) { vn[0][j]=-vn[1][j]; vn[MX+1][j]=-vn[MX][j]; } }

//----------------------------------------------------------------------------

/*---------------------------------------------------------------------------- 更新压强

入口: un、vn、p、dx、dy，新的速度，当前压强，空间步长； 出口: pn，新的压强。

-----------------------------------------------------------------------------*/

void getp(double un[MX+1][MY+2],double vn[MX+2][MY+1],double p[MX+2][MY+2], double pn[MX+2][MY+2],double dx,double dy) { int i,j; for(i=1;i<=MX;i++) for(j=1;j<=MY;j++) pn[i][j]=p[i][j]-dt*c2/dx*(un[i][j]-un[i-1][j]+vn[i][j]-vn[i][j-1]); for(i=1;i<=MX;i++) { pn[i][0]=pn[i][1]; pn[i][MY+1]=pn[i][MY];

程序

} for(j=0;j<=MY+1;j++) { pn[0][j]=pn[1][j]; pn[MX+1][j]=pn[MX][j]; } }

/*------------------------------------------------------ 主程序

-------------------------------------------------------*/ void main() { double dx,dy,err,value,uo,vo,po; int i,j,step; init(u,v,p,dx,dy); err=100.0; step=0; while(err>1e-4 && step<1e6) //err<1e-5，定常解判据；step，限制迭代次数 { printf("step=%d\n",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;j<=MY+1;j++) { value=fabs(un[i][j]-u[i][j])/dt; if(value>err) err=value; u[i][j]=un[i][j]; } } for(i=0;i<=MX+1;i++) { for(j=0;j<=MY;j++) { value=fabs(vn[i][j]-v[i][j])/dt; if(value>err) err=value; v[i][j]=vn[i][j]; } }

程序

for(i=0;i<=MX+1;i++) { for(j=0;j<=MY+1;j++) { value=fabs(pn[i][j]-p[i][j])/c2/dt; if(value>err) err=value; p[i][j]=pn[i][j]; } } printf("err=%le\n",err);

//-------------------------------------------------------- } //输出结果，用Tecplot数据格式画图

FILE *fp;

fp=fopen("Result.plt","w");

fprintf(fp,"variables=x,y,u,v,p\n zone i=%d,j=%d,f=point\n",MX+1,MY+1); for(i=0;i<=MX;i++) { for(j=0;j<=MY;j++) { uo=0.5*(u[i][j]+u[i][j+1]); vo=0.5*(v[i][j]+v[i+1][j]);

po=0.25*(p[i][j]+p[i+1][j]+p[i][j+1]+p[i+1][j+1]); fprintf(fp,"%20.10e %20.10e %20.10e %20.10e %20.10e\n",i*dx,j*dy,uo,vo,po); } }

fclose(fp);

----------------------------------------------------------------------------------------------------- -----------------------------------------------------------------------------------------------------

2. Fortran77语言源程序

!———————————————————————————————————— !利用有限体积算法三阶迎风型QUICK离散格式和

!人工压缩算法求解方腔流动问题（Fortran77语言版本）

!———————————————————————————————————— 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(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) err=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 ! !

subroutine 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