有限体积法求解二维方腔流（一）——动量方程和连续性方程的离散

CFD_Tyro

于 2022-05-03 15:26:58 发布

阅读量3.2k

点赞数 5

分类专栏：理论基础文章标签：算法线性代数

本文链接：https://blog.csdn.net/CFD_Tyro/article/details/124555888

版权

1. 有限体积法求解二维方腔流–理论手册

1.1. 不可压缩流体控制方程

连续性方程
$\nabla \cdot U=0 \tag{1}$
动量方程
$\frac{\partial U}{\partial t}+\nabla\cdot (UU)=-\nabla p + \nabla \cdot (\nu \nabla U) \tag{2}$
其中 $p$ 为真实压力与流体密度之比，量纲为 $m^2/s^2$ 。

上式中的扩散项为空间的二阶导数项，为了保证精度，对上式离散的空间精度也应至少满足二阶精度（Jasak, 1996, phd）

有限体积法对上式离散的第一步是积分：
$\int_t^{t+\Delta t}\left[ \frac{\partial }{\partial t}\int_{V_P} U dV + \int_{V_P}\nabla\cdot (UU) dV\right]dt=\int_t^{t+\Delta t}\left[ -\int_{V_P}\nabla p dV + \int_{V_P}\nabla \cdot (\nu \nabla U) dV \right]dt$

1.2. 逐项离散

使用网格质心值代替网格平均值，有二阶精度，详情参考
$\int_{V_p} U dV = U_P V_P + \varepsilon(\Delta x^2)$

1.2.1. 非稳态项

$\int_t^{t+\Delta t} \frac{\partial }{\partial t}\int_{V_P} U dVdt\approx \frac{U^{t+\Delta t}_P-U^t_P}{\Delta t} \Delta t V_P = (U_P^{n+1}-U_P^n) V_P$

下文中以上标 $(n + 1)$ 替代上标 $(t+\Delta t)$ ，以 $n$ 替代 $t$

1.2.2. 对流项

结合高斯散度定理，
$\int_{V_P}\nabla\cdot (UU) dV=\int_{\partial V_P} (UU) dS\approx\sum(UU)_f S$

注意上式中 $UU)_f$ 为张量在面心上的平均值，二阶精度； $S$ 为矢量，将其乘积做线性化处理
$\sum(UU)_f S \approx \sum F^n_f U^{n+1}_f$

其中 $F_f=S \cdot U_f$ ，为标量。本文中的面法向方向由内指向外，这与OpenFOAM中根据owner和neighbor来定义面法向的方式不同。对于二维情况， $U = (u, v)$ ：
$\sum F_f^n U_f^{n+1}=F_wU_w^{n+1}+ F_eU_e^{n+1}+ F_sU_s^{n+1}+ F_nU_n^{n+1}$

$F_w=U_w\cdot S_w=(u_w,v_w)\cdot (-h,0)=-u_wh=-\frac{u_P+u_W}{2}h$

$F_e=U_e\cdot S_e=(u_e,v_e)\cdot (h,0)=u_eh=\frac{u_P+u_E}{2}h$

$F_s=U_s\cdot S_s=(u_s,v_s)\cdot (0,-h)=-v_sh=-\frac{v_P+v_S}{2}h$

$F_n=U_n\cdot S_n=(u_n,v_n)\cdot (0,h)=v_nh=\frac{v_P+v_N}{2}h$

上式中 $U_f^{n+1}$ 的获取是有限体积法中的关键。介绍几种常用格式：

一阶迎风
$\phi_e=\left\{ \begin{matrix} \phi_P & F_e>0 \\ \phi_E & F_e<0 \end{matrix} \right.$
二阶迎风
$\phi_e=\left\{ \begin{matrix} 0.5(3\phi_P-\phi_W) & F_e>0 \\ 0.5(3\phi_E-\phi_{EE}) & F_e<0 \end{matrix} \right.$
中心差分
$\phi_e=\frac{\phi_P+\phi_E}{2}$

上式中 $\phi$ 为标量，下标 $W, E$ 分别为 $P$ 网格的左和右邻居网格， $e$ 表示 $P$ 网格的右面。如果不是六面体网格（三维）或矩形网格（二维），则上述格式需要修正。

OpenFOAM自带cavity算例中，system/fvSchemes-divSchemes-div(phi,U)使用的离散格式是Gauss linear，其中phi是已知值，所以线性插值仅是应用于变量U，换句话说，cavity默认的对流项离散方式是中心差分

按照中心差分， $U_f=(U_N+U_P)/2$
$U_w^{n+1}=\frac{U_P^{n+1}+U_W^{n+1}}{2}$

$U_e^{n+1}=\frac{U_P^{n+1}+U_E^{n+1}}{2}$

$U_s^{n+1}=\frac{U_P^{n+1}+U_S^{n+1}}{2}$

$U_n^{n+1}=\frac{U_P^{n+1}+U_N^{n+1}}{2}$

$\sum F^n_fU^{n+1}_f = (-u_w^n h)U_w^{n+1} + (u_e^n h)U_e^{n+1} + (-v_s^n h)U_s^{n+1} + (v_n^n h)U_n^{n+1} \\ =\frac{h}{4}\left[ -(u_P^n+u_W^n)(U^{n+1}_P+U^{n+1}_W)+(u_P^n+u_E^n)(U^{n+1}_P+U^{n+1}_E)-(v_P^n+v_S^n)(U^{n+1}_P+U^{n+1}_S)+(v_P^n+v_N^n)(U^{n+1}_P+U^{n+1}_N)\right] \\ =\frac{h}{4}\left[ (u_E^n-u_W^n+v_N^n-v_S^n)U^{n+1}_P -(u_P^n+u_W^n)U^{n+1}_W+(u_P^n+u_E^n)U^{n+1}_E-(v_P^n+v_S^n)U^{n+1}_S+(v_P^n+v_N^n)U^{n+1}_N \right]$