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

附录F 二维不可压缩黏性流体方腔流动问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 的 有限体积算法与计算程序 二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型 离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用 语言和 语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。 F-1利用有限体积算法三阶迎风型 离散格式求解 二维不可压缩黏性流体方腔流动问题 1.二维不可压缩黏性流体方腔流动问题 二维不可压缩黏性流体方腔流动（cavity flow）：有一正方形腔室，其量纲为一的宽度为 ，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为 它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度 沿着上壁面方向自左向右运动（图F.1）。 2. 基本方程组、初始条件和边界条件 设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S方程组来 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示，把它写成通用变量的微分方程组形式，有： （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