fvm_quick_cavity.zip_FVM_FVMc_Quick_二维方腔_二维方腔问题资源-CSDN文库

共1个文件

cpp：1个

版权申诉

5星 · 超过95%的资源 175 浏览量 2022-09-22 22:59:15 上传评论收藏 3KB ZIP 举报

标题 "fvm_quick_cavity.zip" 提供了一个关于FVM（有限体积方法）在解决二维不可压缩黏性方腔流动问题的应用。FVM是一种广泛用于流体动力学、电磁学等领域的数值计算方法，它通过将连续域离散化为一系列互不重叠的控制体积，并基于这些体积来近似偏微分方程的解。描述中提到的"二维不可压缩黏性方腔流动问题"是指在二维空间内，研究一个充满流体的方形腔体内，由于黏性效应引起的流体流动情况。在这种问题中，流体被视为不可压缩的，意味着它的密度在整个过程中保持恒定。黏性则考虑了流体内部的摩擦力，是纳维-斯托克斯方程中的关键因素。 "二维方腔问题"通常用作测试案例，用于验证数值求解方法的准确性和稳定性。在这个问题中，腔体的边界条件通常设定为周期性的，即两侧的流动速度相等，上下面为无滑移边界，这意味着流体与壁面之间没有相对运动。这样的设定使得流体在腔体内形成一个典型的对流结构，即著名的“Poiseuille流动”或“四方腔流动”。标签 "fvm_c" 指的是使用C语言实现的FVM算法。C语言是一种高效且灵活的编程语言，常被用于科学计算和工程应用。"quick"可能指的是快速算法或快速求解策略，这在数值模拟中通常是为了提高计算效率。压缩包内的文件 "fvm_quick_cavity.cpp" 是C++源代码文件，它包含了实现上述FVM算法的程序。这个程序可能包括以下部分： 1. **初始化**：设置网格、边界条件和初始条件。 2. **离散化**：将连续的Navier-Stokes方程转化为离散形式，例如有限体积形式。 3. **时间推进**：利用时间步进方法如欧拉法或龙格-库塔法来更新流场信息。 4. **边界处理**：处理边界条件，如周期性边界和无滑移边界。 5. **求解器**：设计高效的迭代求解器来求解得到的代数系统，如Gauss-Seidel法或雅可比法。 6. **输出和可视化**：可能包含数据输出和结果的后处理，如使用ParaView或MayaVi进行可视化。这个压缩包提供的资源是一个用C++实现的FVM代码，用于解决经典的二维不可压缩黏性方腔流动问题。对于学习和理解FVM、流体力学数值模拟以及C/C++编程实践的人员来说，这是一个非常有价值的实例。

资源详情

资源评论

资源推荐

收起资源包目录

fvm_quick_cavity.zip （1个子文件）

fvm_quick_cavity.cpp 9KB

// 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-