有限体积法求解二维方腔流（三）——代码以及与icoFoam结果对比

理论手册内容请参考理论一与理论二

2. 有限体积法求解二维方腔流–求解

二维模拟区域尺寸为$x\in [0,0.1],y\in [0,0.1]$，每一维各有$N$个网格，网格长度为$h=1/N$，网格体心间距$d=h$，物理量均定义在网格体心。流体粘度为$\nu$。待求物理量为速度矢量$U=[u,v]$以及压力标量$p$，上边界$x$方向速度$u=1$，$y$方向速度$v=0$，其他三边$u=v=0$，四边的压力边界条件均为第二类边界条件即$\partial p/\partial x=\partial p/\partial y=0$

离散格式：对流项中心差分，扩散项中心差分，压力项中心差分，非稳态项一阶欧拉

模拟区域左下角坐标为$(0,0)$，右上角坐标为$(1,1)$，网格全局编号从0开始，令网格$[i,j]$表示该网格为$x$方向的第$i$个网格、$y$方向的第$j$个网格。对于二维，标量$S_f$为表面积，等于$h$

2.1. 代码链接

python代码链接

2.2. 求解流程

1. 给定$t_0$时刻网格值，所有网格$u=v=p=0$
2. 根据$t_0$时刻的$uv$计算系数$A_P,A_W,A_E,A_S,A_N$，对应update_coefficient函数
3. 根据$t_0$时刻的$p$计算$uv$，记为$u^{\ast },v^{\ast }$，对应update_velo_GS函数
4. 进入内循环，循环次数$nCorr$
   • 由$u^{\ast },v^{\ast }$计算$HbyA{|}^{\ast }$，由t_0时刻速度更新$A_P$，对应update_HbyA_AP函数
   • 组建压力泊松方程，由$HbyA{|}^{\ast }$计算$p$，记为$p^{\ast }$，对应update_pressure_GS函数
   • 由$p^{\ast }$更新速度，对应update_velo_final函数
   • 计算连续性方程误差，对应update_continuity_error函数。开始下一次内循环
5. 时间推进，重复以上步骤

可以看出，以上步骤对应PISO算法

2.3. 注意事项

• 由$t_0$时刻推进到$t_1$时刻的过程中，假设求出了中间速度$U^{\ast }$

  • 系数$A_P,A_N$均由$t_0$时刻的速度求解
  • HbyA由最新速度$U^{\ast }$更新

• 压力泊松方程中，不要忘记$V_P$，推导过程中很容易忽视$V_P\nabla p=\sum p_{f}S$，有了