fvm_quick_cavity.zip_FVM_FVM c_Quick_二维方腔_二维方腔问题

// fvm_quick_cavity.cpp /*---------------------------------------------------------------------------- !利用有限体积算法三阶迎风型 离散格式和 !人工压缩算法求解方腔流动问题（ 语言版本） -------------------------------------------------------------------------------*/ #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; } /*----------------------------------------------------------------------------------------------- 一阶迎风型离散格式 二维的三阶迎风型 离散格式为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-