同位网格Rhie-Chow插值


李东岳


1. 动量方程

Rhie-Chow插值是C.M. Rhie和W.L. Chow在1983年对面速度进行速度插值进而消除压力波的一种方法。在Rhie-Chow插值之前,主要采用的是错位网格技术。在Rhie-Chow插值之后至今,大量的CFD求解器采用同位网格。考虑常粘度,无时间项,线性化后的一维不可压缩动量方程:

\begin{equation}\label{00} \frac{\rd u_1u_1}{\rd x}-\frac{\rd }{\rd x}\left(\nu\frac{\rd u_1}{\rd x}\right)=-\frac{\rd p}{\rd x} \end{equation}
其中$u$为$x$方向速度,$p$表示除掉密度的压力,$\nu$表示粘度。其离散形式为:
\begin{equation}\label{0} u_{\mathrm{P}}^n \left(\frac{\mathrm{d} u}{\mathrm{d}x}\right)_\rP-\left(\nu\frac{\mathrm{d}^2 u}{\mathrm{d}x^2}\right)_\rP=-\left(\frac{\mathrm{d}p}{\mathrm{d}x}\right)_\rP \end{equation}
$u_{\mathrm{P}}^n$为当前网格点、当前时间步已知的速度,$_\mathrm{P}$为当前网格点。如果我们对面速度进行中心插值,对对流项进行有限体积离散有:
\begin{equation}\label{1} \int^e _w u_{\mathrm{P}}^n \left(\frac{\mathrm{d} u}{\mathrm{d}x}\right)_\rP \mathrm{d}x = u_{\mathrm{P}}^n u_e - u_{\mathrm{P}}^n u_w = u_{\mathrm{P}}^n \frac{u_{\mathrm{E}}+u_{\mathrm{P}}}{2} - u_{\mathrm{P}}^n \frac{u_{\mathrm{P}}+u_{\mathrm{W}}}{2} = u_{\mathrm{P}}^n \frac{u_{\mathrm{E}}-u_{\mathrm{W}}}{2} \end{equation}
可以看出中心格式当前网格点中心的对角系数为0。对扩散项进行离散有:
\begin{equation}\label{2} \int^e _w \left(\nu\frac{\mathrm{d}^2 u}{\mathrm{d}x^2}\right)_\rP \rd x= \left(\nu\frac{\mathrm{d}u}{\mathrm{d}x}\right)_e-\left(\nu\frac{\mathrm{d}u}{\mathrm{d}x}\right)_w = \nu\left(\frac{u_{\mathrm{E}}-u_{\mathrm{P}}}{\Delta x} - \frac{u_{\mathrm{P}}-u_{\mathrm{W}}}{\Delta x}\right) = \nu\left(\frac{u_{\mathrm{E}}- 2u_{\mathrm{P}}+u_{\mathrm{W}}}{\Delta x}\right) \end{equation}
其离散后为一个对称的矩阵结构。对压力项进行离散有:
\begin{equation}\label{3} -\int^e _w \left(\frac{\mathrm{d}p}{\mathrm{d}x}\right)_\rP \mathrm{d}x= -(p_e - p_w ) \end{equation}
其中$_\mathrm{E,W}$为相邻节点,$\Delta x$为节点的距离并默认为1。在这里假定$\nu=1$,同时将方程\eqref{1}、\eqref{2}、\eqref{3}代入到\eqref{0}中有下面的紧凑形式:
\begin{equation}\label{tri} u_{\mathrm{P}}=\frac{H_\rP}{A_\rP}-\frac{1}{A_\rP}(p_e-p_w) \end{equation}
其中
\begin{equation}\label{H} H_\rP=-\left(\frac{u_\rP^n}{2}-\frac{1}{\Delta x}\right)u_\mathrm{E}-\left(-\frac{u_\rP^n}{2}-\frac{1}{\Delta x}\right)u_\mathrm{W} \end{equation}
\begin{equation} A_\rP=\frac{2}{\Delta x} \end{equation}
若对压力使用线性插值 \begin{equation}\label{p1} p_e=\frac{p_\mathrm{E}+p_\mathrm{P}}{2} \end{equation} \begin{equation}\label{p2} p_w=\frac{p_\mathrm{W}+p_\mathrm{P}}{2} \end{equation} 有:
\begin{equation}\label{pdiff} u_{\mathrm{P}}=\frac{H_\rP}{A_\rP}-\frac{1}{A_\rP}\frac{p_\mathrm{E}-p_\mathrm{W}}{2} \end{equation}
方程\eqref{pdiff}中的压力插值我们采用了二阶精度的中心格式,这导致了\eqref{pdiff}中引入的压力为$\frac{p_\mathrm{E}-p_\mathrm{W}}{2}$,我们称$p_\mathrm{E}-p_\mathrm{W}$为$2-\delta$压差。相反的,如果我们对\eqref{pdiff}使用一阶向前或向后差分,\eqref{pdiff}并不会出现$2-\delta$压差,但是压力作为驱动流体的一种驱动力,且压力的传递不具备方向性。迎风格式不符合物理且精度较低,传统方式动量方程中$2-\delta$压差不可避免。

2. 压力方程

现在考虑压力泊松方程。对于一维不可压缩流体,有连续性方程:

\begin{equation} \left(\frac{\mathrm{d} u}{\mathrm{d}x} \right)_\rP=0 \end{equation}
对其离散有:
\begin{equation}\label{contd} u_e-u_w =0 \end{equation}
如果把方程\eqref{contd}中的面速度采用常规的中心插值处理,有
\begin{equation}\label{u1} u_e=\frac{u_\mathrm{E}+u_\mathrm{P}}{2} \end{equation} \begin{equation}\label{u2} u_w=\frac{u_\mathrm{W}+u_\mathrm{P}}{2} \end{equation}
并代入到方程\eqref{contd}中有:
\begin{equation}\label{pf} \frac{1}{2}(u_\mathrm{E}-u_\mathrm{W})=0 \end{equation}
基于上文动量方程导出的方程\eqref{pdiff},同样的我们有:
\begin{equation}\label{E} u_{\mathrm{E}}=\frac{H_\rE}{A_\rE}-\frac{1}{A_\rE}\frac{p_\mathrm{EE}-p_\mathrm{P}}{2} \end{equation}
\begin{equation}\label{W} u_{\mathrm{W}}=\frac{H_\rW}{A_\rW}-\frac{1}{A_\rW}\frac{p_\mathrm{P}-p_\mathrm{WW}}{2} \end{equation}
将方程\eqref{E}、\eqref{W}代入到\eqref{pf},即为离散后的压力泊松方程:
\begin{equation}\label{pfinal} \frac{H_\rE}{A_\rE}-\frac{H_\rW}{A_\rW}=\frac{1}{2A_\rE}p_\mathrm{EE}-\left(\frac{1}{2A_\rE}+\frac{1}{2A_\rW}\right)p_\mathrm{P}+\frac{1}{2A_\rW}p_\mathrm{WW} \end{equation}
可见,即使$_\mathrm{P}$节点的压力值被引入,但方程\eqref{pfinal}中引入压力也并非相邻网格点的压力值。早期棋盘型压力分布的解决办法是使用错位网格来在动量方程中引入$1-\delta$压差。

3. Rhie-Chow插值

Rhie-Chow插值的精髓在于处理压力插值,即摒弃传统插值方式,如方程\eqref{u1}、\eqref{u2}。类似的,Rhie-Chow插值认为面上的值可以按照方程\eqref{tri}的形式来进行: \begin{equation}\label{rc} u_e=\frac{H_e}{A_e}-\frac{1}{A_e}(p_\rE-p_\rP) \end{equation} 其中$\frac{H_e}{A_e}$不会影响棋盘分布。类似的,将方程\eqref{rc}代入到\eqref{contd}有

\begin{equation}\label{rcfinal} \frac{H_e}{A_e}-\frac{H_w}{A_w}=\frac{1}{A_e}p_\mathrm{E}-\left(\frac{1}{A_e}+\frac{1}{A_w}\right)p_\mathrm{P}+\frac{1}{A_w}p_\mathrm{W} \end{equation}
对比方程\eqref{pfinal},方程\eqref{rcfinal}成功的引入了引入相邻网格点的压力分布。其可在同位网格中针对压力修正方程的方法来消除棋盘型压力分布。

4. OpenFOAM实现

在OpenFOAM中,其通过对拉普拉斯项的离散巧妙的获得一种Rhie-Chow插值的方法。参考icoFoam解析有:

\begin{equation} \nabla \cdot (\mathbf{HbyA}) = \nabla \cdot\left(\frac{1}{A_{\mathrm{P},f}} \nabla p\right) \label{poss} \end{equation}
其中如何对右端的压力梯度项进行离散是问题的关键。在OpenFOAM中,对右端进行积分并结合散度定理有(略去$\frac{1}{A_{\mathrm{P},f}}$):
\begin{equation} \int \nabla \cdot\left( \nabla p\right)\rd V=\sum \left( \nabla p\right)_f\cdot\bfS_f \label{poss2} \end{equation}
其中$\sum \left( \nabla p\right)_f\cdot\bfS_f$定义为面法向梯度,直接和相邻网格点相连接,因此引入相邻网格点的压力梯度。对方程\eqref{poss2}离散后反推,可以获得icoFoam解析中方程(23)的插值公式。类似的,参考方程\eqref{pdiff},在动量方程中,若直接对体积力进行离散,也会产生压力振荡问题。另一种方式是首先将梯度项转换为面法向梯度(直接和相邻网格点的数值相连接),然后再重组到体心,同样可以消除震荡。即对应OpenFOAM中对体积力处理的fvc::reconstruct

2020.02.03CloudBird:修正方程\eqref{H}

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网