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

实用文案附录F 二维不可压缩黏性流体方腔流动问题的有限体积算法与计算程序二维方腔流动问题是一个不可压缩黏性流动中典型流动。 虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。 在本算例中采用有限体积算法三阶迎风型QUICK离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用C语言和Fortran77语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。F-1利用有限体积算法三阶迎风型 QUICK离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（ cavity flow）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为p=1.0、=1.0它周围壁面（左右壁面和

图F.1二维不可压缩黏性方腔流问题示意图底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。基本方程组、初始条件和边界条件设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N-S方程组来表示，把它写成通用变量的微分方程组形式，有：图F.1二维不可压缩黏性方uv腔流动问题示意图（F.1）Sxyxxyy标准实用文案其中u为变量 在水平x方向的流速，v为 在垂直y方向的流速， 为黏度，S为源项。源项中不仅包含压力梯度项，也包含时间导数项。初始条件：方腔上壁面以量纲为一的速度 1.0沿着上壁面方向自左向右运动。边界条件：流动速度u、v均可采用无滑移边界条件， 压力p采用自由输出边界条件。3.计算网格划分和控制体单元与节点定义采用交错网格，图 F.2和图F.3是计算网格、控制体单元和节点示意图。图F.2方腔流动计算网格、图F.3计算采用的交错网格示意图控制体单元和节点示意图节点P所在主控制单元如图 F.2中有阴影部分所示。在 x方向与节点P相邻的节点为W和E，在y方向与节点P相邻的节点为S和N，主控制单元界面分别为w、s、e、n。压力p和速度u、v分别在三套不同网格中如图 F.3中有阴影部分所示。4.有限体积算法三阶迎风型 QUICK离散格式标准实用文案对方程（F.1）在图F.2所示节点P所在控制体单元内积分，有：udVVyvdVVxdVVydVSdV（F.2）VxxyV由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积V仅是面积，而它的边界是长度。设AwAey,AsAnx，利用Gauss定理，可将方程（F.2）改写成如下有限体积算法离散格式：uAeuAwvAnvAsAAAAS（F.3）xyxexwynys对上式中、、、采用一阶向前差分近似，则有：xexwynysxy

E P,e x xN P,n y y

Wwx（F.4）Ssy同时记：FeueAe,FwuwAw（F.5）FnvnAn,FsvsAsDeeAe,DwwAw,xPExPW（F.6）nAn,DssAsDnyPNyPS则可由式（F.2）写成：FeeFwwFnnFssDeEPDwPW（F.7）DnNPDsPSSxy式中P、E、W、N、S、De、Dw、Dn、Ds都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量Fee、Fww、Fnn、Fss，就可以求出节点上未知量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计算公式如下：当uw0时，界面参数w计算公式为：636（F.8a）ww8wW8wP8wWW当ue0时，界面参数w计算公式为：636（F.8b）e8P8E8W对于一维无源项一维对流扩散方程三阶迎风型QUICK离散格式：当ue0,uw0时，三阶迎风型QUICK离散格式为：标准实用文案aP P aWWaEEaWWWW其中aWDw3Fw1Fe88aD3FEe8eaWW1Fw8aPaWaEaWWFeFw同理，若ue0,uw0，三阶迎风型QUICK离散格式为：aPPaWWaEEaEEEE其中aWDw3Fw,8aEDe6Fe1Fw,88aWW1Fe8aPaWaEaEEFeFw

F.9）F.9a）F.10）F.10a）将两种流动方向离散方程（ F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型 QUICK离散格式：aPPaWWaEEaEEEEaEEEE（F.11）其中aWDw6wFw1eFe31wFw888aWW1wFw8aEDe3eFe61eFe11wFw（F.11a）888aEE11eFe8aPaWaEaWWaEEFeFw标准实用文案式中当Fw时，w；01当Fw时，w；00（F.11b）当Fe时，e；01当Fe时，e。00同理，可以得到带有源项的二维对流扩散方程三阶迎风型QUICK离散格式为：aPPaWWaEEaEEEEaEEEE（F.12）aSSaNNaSSSSaNNNNSV其中S为有限体积算法中源项平均值。式中各个系数为：aWDw6wFw1eFe31wFw888aWW1wFw8361aEDe8eFe81eFe81wFw1aEE1eFe86131aSDssFsnFnsFs888aSS1sFs8aNDn3nFn61nFn11sFs888aNN11nFn8aPaWaEaWWaEEaSaNaSSaNNFeFwFnFs（F.12a）式中当Fk0时，当Fk0时，

kk

1；0。 （F.12b）k w,e,s,n源项S为：up）S（F.13tx标准实用文案nn1S离散若把u表示tn时刻动量，u表示tn1时刻动量，则可以得到源项格式为：n1nSdVuPuPxypepwy（F.14）tV最后，得到有限体积算法二维对流扩散方程三阶迎风型QUICK离散格式：aPuPnaWuWnaEuEnaSuSnaNuNnn1n（F.15）uPuPxynnypepwt式中系数ak为一阶迎风格式中各对应系数。5.计算结果分析利用三阶迎风型QUICK离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。图 F.5是不同雷诺数Re条件下采用三阶迎风型QUICK离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。 这表明有限体积算法三阶迎风型 QUICK离散格式具有相当高的计算精度。Re=100 Re=1000标准实用文案Re=5000 Re=10000图F.5不同雷诺数 Re条件下采用三阶迎风型 QUICK离散格式计算二维不可压缩黏性方腔流动的计算结果从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。这是因为，方腔上壁面以量纲为一的速度 1.0沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件，这是一个带奇性的方腔流动。计算结果表明，方腔流动和雷诺数有关，随着雷诺数 Re增加，计算精度在降低。当雷诺数 Re较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数Re的增加，二次小涡变得越来越模糊。由于三阶迎风型 QUICK离散格式计算精度较高，因此三阶迎风型 QUICK离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。标准实用文案F-2 二维不可压缩黏性方腔流动问题数值计算源程序C语言源程序//fvm_quick_cavity.cpp/*----------------------------------------------------------------------------!利用有限体积算法三阶迎风型 QUICK离散格式和!人工压缩算法求解方腔流动问题（ C语言版本）-------------------------------------------------------------------------------*/#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMX100#defineMY100#defineRe100.00#definedt0.0005#definec22.25Doubleu[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];标准实用文案全局变量/*-------------------------------------------------------定义两个要用到的函数--------------------------------------------------------*/doublemax(doublea,doubleb){if(a<b)returnb;elsereturna;}doublealfa(doublex){if(x>=0)return1.0;elsereturn0.0;}/*------------------------------------------------------------------------标准实用文案初始化入口:无；出口：u、v、p、dx、dy，初始速度、压强和空间步长。---------------------------------------------------------------------------*/voidinit(doubleu[MX+1][MY+2],doublev[MX+2][MY+1],doublep[MX+2][MY+2],double&dx,double&dy){inti,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],doublep[MX+2][MY+2],doubleun[MX+1][MY+2],doubledx,doubledy,inti,intj){doubleaw,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],doublep[MX+2][MY+2],doublevn[MX+2][MY+1],doubledx,doubledy,inti,intj){doubleaw,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],doublep[MX+2][MY+2],doubleun[MX+1][MY+2],doublevn[MX+2][MY+1],doubledx,doubledy){doublemiu,fw,fe,fs,fn,df,aw,ae,as,an,aww,aee,ass,ann,ap;inti,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],doublep[MX+2][MY+2],doublepn[MX+2][MY+2],doubledx,doubledy){inti,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];}标准实用文案}/*------------------------------------------------------主程序-------------------------------------------------------*/voidmain(){doubledx,dy,err,value,uo,vo,po;inti,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\nzonei=%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);标准实用文案----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------Fortran77语言源程序————————————————————————————————————!利用有限体积算法三阶迎风型 QUICK离散格式和人工压缩算法求解方腔流动问题（Fortran77语言版本）————————————————————————————————————programQUICK_cavityparameter(mx=101,my=101)implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)dimensionun(mx,my+1),vn(mx+1,my),pn(mx+1,my+1)common/ini/u,v,pc2=2.25标准实用文案re=100.0dt=0.0005dx=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=0err=100.00!nun，计数器；err，判断人工压缩法求解收敛的标准callinitial调入初始条件，以下为人工压缩算法求解err=0.0callquick(u,v,p,un,vn,mx,my,dx,dy,dt,re)!QUICK离散格式求解动量方程，得到 un、vncallcalp(p,un,vn,pn,mx,my,dx,dy,dt,c2)求压强pncallcheck(u,v,p,un,vn,pn,mx,my,dx,dy,dt,c2,err)校验流场信息，判断是否收敛，同时更新u、v、pwrite(*,*)'error=',err标准实用文案num=num+1write(*,*)num屏幕跟踪输出enddocalloutput(u,v,p,mx,my,dx,dy)输出结果文件End!!subroutineinitial初始化流场parameter(mx=101,my=101)doubleprecisionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)common/ini/u,v,pdoi=1,mx+1doj=1,my+1p(i,j)=1.0enddoenddodoi=1,mxdoj=1,my+1u(i,j)=0.0标准实用文案enddoenddodoi=1,mx+1doj=1,myv(i,j)=0.0enddoenddoendsubroutine!!subroutinequick(u,v,p,un,vn,mx,my,dx,dy,dt,re)以QUICK格式离散动量方程implicitdoubleprecision(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)doubleprecisionmiumiu=1.0/re! 以 下 求 解 x 方 向 速 度un----------------------------------------------------------------------------doi=3,mx-2标准实用文案doj=3,my-1fw=0.5*(u(i-1,j)+u(i,j))*dyfe=0.5*(u(i,j)+u(i+1,j))*dyfs=0.5*(v(i,j-1)+v(i+1,j-1))*dxfn=0.5*(v(i,j)+v(i+1,j))*dxdf=fe-fw+fn-fsaw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=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))/dxenddoenddo标准实用文案!-------------------------------------------------------------------------------------------j=2doi=3,mx-2callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoj=mydoi=3,mx-2callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoi=2doj=2,mycallupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoi=mx-1doj=2,mycallupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddo!内层边界由一阶迎风型离散格式得到 ----------------------------------------------------!-------------------------------------------------------------------------------------------doi=2,mx-1标准实用文案un(i,1)=-un(i,2)un(i,my+1)=2.0-un(i,my)enddodoj=1,my+1un(1,j)=0.0un(mx,j)=0.0enddo外层边界条件按物理边界条件给出!-----------------------------------------------------------------------------以下求解y方向速度vn-----------------------------------------------doi=3,mx-1doj=3,my-2fw=0.5*(u(i-1,j)+u(i-1,j+1))*dyfe=0.5*(u(i,j)+u(i,j+1))*dyfs=0.5*(v(i,j-1)+v(i,j))*dxfn=0.5*(v(i,j)+v(i,j+1))*dxdf=fe-fw+fn-fsaw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=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))*fsass=-0.125*alfa(fs)*fsan=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+ae+as+an+aww+aee+ass+ann+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)+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))/dyenddoenddo!-----------------------------------------------------------------------------j=2doi=3,mx-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddoj=my-1doi=3,mx-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddoi=2doj=2,my-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)标准实用文案enddoi=mxdoj=2,my-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddo!----------------------------------------------------------------------------doi=2,mxvn(i,1)=0.0vn(i,my)=0.0enddodoj=1,myvn(1,j)=-vn(2,j)vn(mx+1,j)=-vn(mx,j)enddo!----------------------------------------------------------------------------Endsubroutine!!functionalfa(x)函数,正1负0doubleprecisionalfa,x标准实用文案alfa=1.0elsealfa=0.0endifend!!subroutineupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)以一阶迎风型离散格式得到内层边界un的值implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1)doubleprecisionmiumiu=1.0/reaw=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))*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*(p(i+1,j)-p(i,j))/dxEndsubroutine标准实用文案!!subroutineupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)!以一阶迎风型离散格式得到内层边界 vn值implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),vn(mx+1,my)doubleprecisionmiumiu=1.0/reaw=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