计算流体力学附录F 有限体积法计算方腔流(F)附录F 二维不可压缩黏性流体方腔流动问题的
有限体积算法与计算程序
二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解,但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型
离散格式对它进行数值求解。同时,为了初学者入门和练习方便,这里给出了用
语言和
语言编写的计算二维不可压缩黏性方腔流动问题计算程序,供大家学习参考。
F-1利用有限体积算法三阶迎风型
离散格式求解
二维不可压缩黏性流体方腔流动问题
1.二维不可压...
附录F 二维不可压缩黏性流体方腔流动问题的
有限体积算法与计算程序
二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解,但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型
离散
对它进行数值求解。同时,为了初学者入门和练习方便,这里给出了用
语言和
语言编写的计算二维不可压缩黏性方腔流动问题计算程序,供大家学习参考。
F-1利用有限体积算法三阶迎风型
离散格式求解
二维不可压缩黏性流体方腔流动问题
1.二维不可压缩黏性流体方腔流动问题
二维不可压缩黏性流体方腔流动(cavity flow):有一正方形腔室,其量纲为一的宽度为
,里面充满静止的不可压缩黏性流体,方腔内初始时刻压力和密度为
它周围壁面(左右壁面和底面)固定不动,上壁面以量纲为一的速度
沿着上壁面方向自左向右运动(图F.1)。
2. 基本方程组、初始条件和边界条件
设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S方程组来表示,把它写成通用变量的微分方程组形式,有:
(F.1)
其中
为变量
在水平
方向的流速,
为
在垂直
方向的流速,
为黏度,
为源项。源项中不仅包含压力梯度项,也包含时间导数项。
初始条件:方腔上壁面以量纲为一的速度
沿着上壁面方向自左向右运动。
边界条件:
流动速度
均可采用无滑移边界条件,压力
采用自由输出边界条件。
3.计算网格划分和控制体单元与节点定义
采用交错网格,图F.2和图F.3是计算网格、控制体单元和节点示意图。
节点
所在主控制单元如图F.2中有阴影部分所示。在
方向与节点
相邻的节点为
和
,在
方向与节点
相邻的节点为
和
,主控制单元界面分别为
。压力
和速度
分别在三套不同网格中如图F.3中有阴影部分所示。
4.有限体积算法三阶迎风型
离散格式
对方程(F.1)在图F.2所示节点
所在控制体单元内积分,有:
(F.2)
由于二维不可压缩黏性流体方腔流动是二维问题,因此控制体单元体积
仅是面积,而它的边界是长度。设
,利用
定理,可将方程(F.2)改写成如下有限体积算法离散格式:
(F.3)
对上式中
采用一阶向前差分近似,则有:
(F.4)
同时记:
(F.5)
(F.6)
则可由式(F.2)写成:
(F.7)
式中
都是控制体单元内节点上的已知量,如果利用差分计算得到控制体单元边界上的流通量
,就可以求出节点上未知量
。
为了便于讨论,现对一维对流扩散方程的三阶迎风型
离散格式进行分析:在三阶迎风型
离散格式中,计算主控制单元界面上流动量
需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到,其中两个节点位于界面紧邻的两侧,第三个节点位于迎风一侧较远邻点,如图F.4所示。
当
时,通过
、
和
三个节点值拟合曲线来计算主控制单元左侧界面参数
。通过节点
、
和
三个节点值拟合曲线来计算主控制单元右侧界面参数
。当
,则分别通过节点
、
、
和
、
、
三个节点值计算主控制单元左、右两侧界面参数
和
。根据上述计算原则,可以得到界面参数
计算
如下:
当
时,界面参数
为:
(F.8a)
当
时,界面参数
计算公式为:
(F.8b)
对于一维无源项一维对流扩散方程三阶迎风型
离散格式:
当
时,三阶迎风型
离散格式为:
(F.9)
其中
(F.9
)
同理,若
,三阶迎风型
离散格式为:
(F.10)
其中
(F.10a)
将两种流动方向离散方程(F.9)和(F.10)合并后,可得到统一的一维对流扩散方程三阶迎风型
离散格式:
(F.11)
其中
(F.11a)
式中
(F.11b)
同理,可以得到带有源项的二维对流扩散方程三阶迎风型
离散格式为:
(F.12)
其中
为有限体积算法中源项平均值。式中各个系数为:
(F.12a)
式中
(F.12b)
源项
为:
(F.13)
若把
表示
时刻动量,
表示
时刻动量,则可以得到源项
离散格式为:
(F.14)
最后,得到有限体积算法二维对流扩散方程三阶迎风型
离散格式:
(F.15)
式中系数
为一阶迎风格式中各对应系数。
5.计算结果分析
利用三阶迎风型
离散格式和相应的初始条件和边界条件,求解二维不可压缩黏性流体方腔流动问题。图F.5是不同雷诺数
条件下采用三阶迎风型
离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。
计算结果和文献中其他高精度算法得到的计算结果进行了比较,两者计算结果十分吻合,能把方腔下壁面两个底角附近二次小涡清晰地计算出来。这表明有限体积算法三阶迎风型
离散格式具有相当高的计算精度。
从图F.5中可以看出:二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置,方腔内流动也并不对称。这是因为,方腔上壁面以量纲为一的速度
沿着上壁面方向自左向右运动时,在方腔上壁面两侧的两个顶角处不再满足边界条件,这是一个带奇性的方腔流动。
计算结果表明,方腔流动和雷诺数有关,随着雷诺数
增加,计算精度在降低。当雷诺数
较低时,方腔下壁面的两个底角附近的二次小涡十分清晰,随着雷诺数
的增加,二次小涡变得越来越模糊。由于三阶迎风型
离散格式计算精度较高,因此三阶迎风型
离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。
F-2 二维不可压缩黏性方腔流动问题数值计算源程序
1.
语言源程序
// fvm_quick_cavity.cpp
/*----------------------------------------------------------------------------
!利用有限体积算法三阶迎风型
离散格式和
!人工压缩算法求解方腔流动问题(
语言版本)
-------------------------------------------------------------------------------*/
#include
#include
#include
#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=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;
}
/*-----------------------------------------------------------------------------------------------
一阶迎风型离散格式
二维的三阶迎风型
离散格式为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;
}
/*----------------------------------------------------------------------
三阶迎风型
离散格式
入口: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);
//--------------------------------------------------------
}
//输出结果,用
数据格式画图
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.
语言源程序
!————————————————————————————————————
!利用有限体积算法三阶迎风型
离散格式和
!人工压缩算法求解方腔流动问题(
语言版本)
!————————————————————————————————————
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
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
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)+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
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,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
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
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=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,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)
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,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(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
do i=2,mx
pn(i,1)=pn(i,2)
pn(i,my+1)=pn(i,my)
enddo
do j=1,my+1
pn(1,j)=pn(2,j)
pn(mx+1,j)=pn(mx,j)
enddo
endsubroutine
!
!
subroutine check(u,v,p,un,vn,pn,mx,my,dx,dy,dt,c2,err)
!校验流场信息,得到收敛判断准则err,同时更新u、v、p
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)
do i=1,mx
do j=1,my+1
abc=abs(un(i,j)-u(i,j))/dt
if(abc.gt.err) err=abc
u(i,j)=un(i,j)
enddo
enddo
do i=1,mx+1
do j=1,my
abc=abs(vn(i,j)-v(i,j))/dt
if(abc.gt.err) err=abc
v(i,j)=vn(i,j)
enddo
enddo
do i=1,mx+1
do j=1,my+1
abc=abs(pn(i,j)-p(i,j))/c2/dt
if(abc.gt.err) err=abc
p(i,j)=pn(i,j)
enddo
enddo
endsubroutine
!
!
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),un(mx,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='out.plt')
write(1,*) 'TITLE = "result"'
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
!
!
-------------------------------------------------------------------
-------------------------------------------------------------------
-------------------------------------------------------------------
图F.1二维不可压缩黏性
方腔流问题示意图
图F.3计算采用的交错网格示意图
图F.2方腔流动计算网格、
控制体单元和节点示意图
图F.4三阶迎风型� EMBED Equation.DSMT4 ���离散格式示意图
图F.1 二维不可压缩黏性方腔流动问题示意图
� EMBED Equation.DSMT4 ���=1 000
� EMBED Equation.DSMT4 ���=100
� EMBED Equation.DSMT4 ���=10 000
� EMBED Equation.DSMT4 ���=5 000
图F.5不同雷诺数� EMBED Equation.DSMT4 ���条件下采用三阶迎风型� EMBED Equation.DSMT4 ���
离散格式计算二维不可压缩黏性方腔流动的计算结果
PAGE
-F.8-
_1341491435.unknown
_1341492992.unknown
_1345287821.unknown
_1345292736.unknown
_1345294576.unknown
_1345287856.unknown
_1341570295.unknown
_1343239008.unknown
_1341493857.unknown
_1341494546.unknown
_1341494890.unknown
_1341494942.unknown
_1341494990.unknown
_1341495058.unknown
_1341494965.unknown
_1341494924.unknown
_1341494743.unknown
_1341494858.unknown
_1341494597.unknown
_1341494067.unknown
_1341494201.unknown
_1341493987.unknown
_1341493485.unknown
_1341493765.unknown
_1341493835.unknown
_1341493687.unknown
_1341493314.unknown
_1341493431.unknown
_1341493206.unknown
_1341491821.unknown
_1341492777.unknown
_1341492859.unknown
_1341492661.unknown
_1341491716.unknown
_1341491742.unknown
_1341491625.unknown
_1262191536.unknown
_1263107967.unknown
_1263108053.unknown
_1263108720.unknown
_1263108923.unknown
_1263109814.unknown
_1263108754.unknown
_1263108668.unknown
_1263108017.unknown
_1263108037.unknown
_1263107985.unknown
_1263059671.unknown
_1263107909.unknown
_1263107933.unknown
_1263107188.unknown
_1263107769.unknown
_1263107203.unknown
_1263061168.unknown
_1263059531.unknown
_1263059550.unknown
_1262267354.unknown
_1263059519.unknown
_1262191663.unknown
_1262192719.unknown
_1262191629.unknown
_1262002541.unknown
_1262072924.unknown
_1262072969.unknown
_1262073230.unknown
_1262084892.unknown
_1262084734.unknown
_1262073091.unknown
_1262072940.unknown
_1262004207.unknown
_1262067436.unknown
_1262002694.unknown
_1262002503.unknown
_1262002528.unknown
_1261145244.unknown
_1261145365.unknown
_1262001401.unknown
_1261145379.unknown
_1261145358.unknown
_1257315109.unknown
_1261143476.unknown
_1261143494.unknown
_1257315263.unknown
_1261143168.unknown
_1257257718.unknown
_1257314795.unknown
_1256934151.unknown
_1256931457.unknown
本文档为【计算流体力学附录F 有限体积法计算方腔流(F)】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。