附录E-有限体积法求解库塔流(E).doc

附录E 二维不可压缩黏性流体流动问题
的有限体积算法与计算程序
二维流动是一个不可压缩黏性流动的典型流动，并有解析解，可以用来检验数值算法计算精度和可靠性。对它采用有限体积算法一阶迎风型离散格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和语言编写的计算二维不可压缩黏性流体流动问题计算程序，供大家学习参考。
E-1利用有限体积算法一阶迎风型离散格式求解二维
不可压缩黏性流体流动问题
1.二维不可压缩黏性流体流动问题的提法
二维不可压缩黏性流体流动：有两个无限长平板，间距为，两板之间充满密度为1、静止的不可压缩黏性流体。无限长平板组成的通道两端压力相等，下板固定不动，上板以量纲为一的速度1自左向右平移运动（图E.1）。

图E.1二维不可压缩黏性流体流动问题示意图

2. 基本方程组、初始条件和边界条件
设流体是黏性流体。二维流动问题在数学上可以由二维不可压缩黏性流动方程组来表示，把它写成通用变量的微分方程组形式，有：
（E.1）
其中为变量在水平方向的流速，为在垂直方向的流速，为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。
初始条件：上板以量纲为一的速度1沿着上壁面方向自左向右运动。
边界条件：流动速度在上下壁面可采用无滑移边界条件，在左右两端采用自由输出边界条件；压强采用自由输出边界条件。
3. 计算网格划分和控制体单元与节点定义
采用交错网格。图E.2和图E.3是计算网格、控制体单元和节点示意图。

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

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

节点所在主控制体单元如图E.2中有阴影部分所示。在方向与节点相邻的节点为和，在方向与节点相邻的节点为和，主控制单元界面分别为。压力和速度分别在三套不同网格中，如图E.3中有阴影部分所示。
4. 有限体积算法离散格式
对方程（E.1）在图E.2所示节点主控制单元内积分，有：
（E.2）
由于不可压缩黏性流体流动是二维问题，因此，控制体单元体积仅是面积，而它的边界是长度。设 ，利用定理，可将方程（E.2）改写成如下有限体积离散格式：
（E.3）
对上式中采用一阶向前差分近似，则有：
（E.4）
同时记：
（E.5）
（E.6）
则式（E.2）写成：
（E.7）
式中都是主控制单元内节点上的已知量，如果利用差分计算得到主控制单元边界上的流通量，就可以求出节点上的未知量。、
为了便于讨论，现对一维对流扩散方程的一阶迎风型离散格式进行分析。首先讨论无源项一维对流扩散方程的一阶迎风型离散格式。
当流动为正向时，即，主控制单元界面取值：
（E.8）
则方程（E.7）离散为：
（E.9）
当流动为负向时，即，控制体单元界面取值为：
（E.10）
则方程（E.7）离散为：
（E.11）
将两种流动方向离散格式（E.9）和（E.11）合并后，得到统一的一维对流扩散方程一阶迎风型离散格式表达式：
（E.12）
式中
（E.12a）
同理，可以得到带有源项的二维对流扩散方程的一阶迎风型离散格式：
（E.13）
式中
（E.13a）
源项为：
（E.14）
若把表示为时刻的动量，表示时刻的动量，则可以得到源项离散格式为：
（E.15）
最后，得到带有源项的二维对流扩散方程有限体积算法一阶迎风型显式离散格式：

（E.16）
式中系数为一阶迎风格式中各对应系数。
5. 计算结果分析
利用上述有限体积算法一阶迎风型离散格式和相应的初始条件和边界条件，采用MAC算法中压力耦合方程，求解二维不可压缩黏性流体流动问题。图E.4 给出了二维不可压缩黏性流体流动沿y方向的速度分布，并和精确解进行了比较，十分吻合。图E.5是二维不可压缩黏性流体流动水平x方向的速度云图。图E.6是二维不可压缩黏性流体流动速度矢量分布图。

图E.4 二维不可压缩黏性流体流动沿y方向的速度分布

图E.5 二维不可压缩黏性流体流动水平x方向速度云图

图E.6 二维不可压缩黏性流体流动速度矢量分布图

E-2 二维不可压缩黏性流体流动问题的数值计算源程序

1. 语言源程序
// fvm_upwind_MAC_couette.cpp
/*-------------------------------------------------------------------------------------------------------------------
以一阶迎风型格式和压力迭代求解二维流动问题（C语言版本）
-----------------------------------------------------------/
#include
#include
#include
#define MX 100 //最大网格数
#define MY 20 //最大网格数
#define Lambda 0.002 //Chorin压力迭代常数，直接影响收敛性
#define Re 100.00 //雷诺数
#define dt 0.0005 //时间步长

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];
//全局变量

double max(double a,double b) //定义max函数
{
double c;

if(a ，初始速度、压强和空间步长。
----------------------------------------------------------------------*/
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=0.2/MY;
for(i=0;i<=MX;i++)
{
for(j=0;j<=MY+1;j++)
{
u[i][j]=00;
if(j==MY+1) u[i][j]=2.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;
}

/*------------------------------------------------------------------------
迭代更新压强
输入：u、v、p、dx、dy ，n时刻速度、压强、空间步长；
输出：p n+1时刻压强。
-----------------------------------------------------------------------------*/
void getp(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+2],double dx,double dy)
{
double d[MX+2][MY+2];
int i,j;
for(i=1;i<=MX;i++)
for(j=1;j<=MY;j++)
d[i][j]=(u[i][j]-u[i-1][j])/dx + (v[i][j]-v[i][j-1])/dy;
for(i=1;i<=MX;i++)
for(j=1;j<=MY;j++)
p[i][j]=p[i][j]-Lambda*d[i][j];
for(i=1;i<=MX;i++)
{
p[i][0]=p[i][1];
p[i][MY+1]=p[i][MY];
}
for(j=1;j<=MY;j++)
{
p[0][j]=p[1][j];
p[MX+1][j]=p[MX][j];
}
}

/*---------------------------------
一阶迎风格式
输入：u、v、p、dx、dy ，n时刻速度、n+1时刻压强,空间步长；
输出：un、vn ，n+1时刻速度。
------------------------------------------------------------------------------------------*/
void upwind(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)
{
int i,j;
double aw,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*(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;
}
}
//边界条件
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]=un[1][j];
un[MX][j]=un[MX-1][j];
}
for(i=1;i<=MX;i++)
{
for(j=1;j<=MY-1;j++)
{
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[i1][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;
}
}
//边界条件
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];
}
}

/*------------------------------------------------
格式输出，采用数据格式画图
输入：u、v、p、dx、dy ，最后结果；
输出：无。
--------------------------------------------------------------------------------*/
void output(double u[MX+1][MY+2],double v[MX+2][MY+1],double p[MX+2][MY+2],double dx,double dy)
{
double uo,vo,po;
int i,j;
FILE *fp;
fp=fopen("Result.plt","w");
fprintf(fp,"variables=x,y,u,v,p\nzone i=%d,j=%d,f=point\r\n",MX+1,MY+1);
for(j=0;j<=MY;j++)
{
for(i=0;i<=MX;i++)
{
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]);
fprintffp,"%20.10e %20.10e %20.10e %20.10e %20.10e\r\n",i*dx,j*dy,uo,vo,po);
}
}
fclose(fp);
}

/*------------------------------------------------------
主程序
-------------------------------------------------------*/
void main()
{
double dx,dy,err,value;
int i,j,step;
init(u,v,p,dx,dy);
err=100.0;
step=0;
while(err>1e-6 && step<1e6) //err，收敛限制；step，迭代次数限制
{
err=0.0;
getp(u,v,p,dx,dy); //更新压强
upwind(u,v,p,un,vn,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];
}
}
printf("err=%le\n",err);
}
output(u,v,p,dx,dy); //输出结果

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

2. 语言源程序
program fvm_upwind_MAC_couette
!-------------------------------------------------------------------
以一阶迎风型离散格式和压力迭代求解二维流动问题
（语言版本）
!-------------------------------------------------------------------------------------------------------------------
parameter(mx=101,my=21)
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 lambda
integer step
re=100. !雷诺数
dt=0.0005 !时间步长
lambda=0.002 !Chorin压力迭代常数，直接影响收敛性
dx=1.0d0/(mx-1)
dy=0.2d0/(my-1)
step=0
err=100.
call initial(u,v,p)
!初始化
do while(err>1e-6.and.step<1e6) !err，收敛限制；step，迭代次数限制
err=0.0
step=step+1
! write(*,*) step
call calp(u,v,p,mx,my,dx,dy,lambda)
!Chorin压力迭代
call upwind(u,v,p,un,vn,mx,my,dx,dy,dt,re)
!一阶迎风型离散格式求解动量方程
do i=1,mx
do j=1,my+1
temp=abs(u(i,j)-un(i,j))/dt
if(temp>err) err=temp
u(i,j)=un(i,j)
enddo
enddo
do i=1,mx+1
do j=1,my
temp=abs(v(i,j)-vn(i,j))/dt
if(temp>err) err=temp
v(i,j)=vn(i,j)
enddo
enddo
write(*,*) 'err=',err
enddo
call output(u,v,p,mx,my,dx,dy) !输出
End
!------------------------------------------------------------------------------
!初始化
!输入：无；
!输出：u、v、p，初始速度、压强。
!-------------------------------------------------------------------------------
subroutine initial(u,v,p)
parameter(mx=101,my=21)
double precision u(mx,my+1),v(mx+1,my),p(mx+1,my+1)
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==my+1) u(i,j)=2.0
enddo
enddo
do i=1,mx+1
do j=1,my
v(i,j)=0.0

enddo
enddo
endsubroutine

!-----------------------------------------------------------------------------------------------
! Chorin压力迭代
!输入：u、v、p、mx、my、dx、dy、lambda ，n时刻速、压强、其他各量；
!输出：p ，n+1时刻压强。
!--------------------------------------------------------------------------------

subroutine calp(u,v,p,mx,my,dx,dy,lambda)
implicit double precision(a-h,o-z)
dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1),d(mx+1,my+1)
double precision lambda
do i=2,mx
do j=2,my
d(i,j)=(u(i,j)-u(i-1,j))/dx+(v(i,j)-v(i,j-1))/dy
enddo
enddo
do i=2,mx
do j=2,my
p(i,j)=p(i,j)-lambda*d(i,j)
enddo
enddo
do i=2,mx
p(i,1)=p(i,2)
p(i,my+1)=p(i,my)
enddo
do j=1,my+1
p(1,j)=p(2,j)
p(mx+1,j)=p(mx,j)
enddo
endsubroutine

!-----------------------------------------------------------------
!一阶迎风型离散格式
!输入：u、v、pn、mx、my、dx、dy、dt、re ，n时刻速度、n+1时刻压强、其他各量；
!输出：un,vn，n+1时刻速度。
!----------------------------------------------------------------------------------------------------------------
subroutine upwind(u,v,pn,un,vn,mx,my,dx,dy,dt,re)

implicit double precision(a-h,o-z)
dimension u(mx,my+1),v(mx+1,my),pn(mx+1,my+1),un(mx,my+1),vn(mx+1,my)
double precision miu
miu=1.0/re
do i=2,mx-1
do j=2,my
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*(pn(i+1,j)-pn(i,j))/dx
enddo
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)=un(2,j)
un(mx,j)=un(mx-1,j)
enddo
!u边界条件
!--------------------------------
!---------------------
do i=2,mx
do j=2,my-1
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)
df0.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*(pn(i,j+1)-pn(i,j))/dy
enddo
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
!v边界条件
!--------------------------------
Endsubroutine

!--------------------------------------------------------------------------------
!格式输出，采用Tecplot数据格式画图
!输入：u、v、p、mx、my、dx、dy，最后结果；
!输出：无。
!--------------------------------------------------------------------------------
subroutine output(u,v,p,mx,my,dx,dy)
implicit double precision(a-h,o-z)
dimension u(mx,my+1),v(mx+1,my),p(mx+1,my+1)
dimension uc(mx,my),vc(mx,my),pc(mx,my),x(mx),y(my)
do i=1,mx
x(i)=(i-1)*dx
enddo
do j=1,my
y(j)=(j-1)*dy
enddo
do i=1,mx
do j=1,my
uc(i,j)=0.5*(u(i,j)+u(i,j+1))
vc(i,j)=0.5*(v(i,j)+v(i+1,j))
pc(i,j)=0.25*(p(i,j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1))
enddo
enddo
open(1,file='result.plt')
write(1,*) 'TITLE = "2D Results"'
write(1,*) 'VARIABLES = "X", "Y", "U", "V", "P"'
write(1,*) 'ZONE I=',mx, 'J=',my, 'F=POINT'
write(1,10) ((x(i),y(j),uc(i,j),vc(i,j),pc(i,j),i=1,mx),j=1,my)
close(1)
10 format(5(e15.8,5x))
endsubroutine

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

-E.15-