• 返回主页
  • Rhie-Chow插值在OpenFOAM中的实现

    2017.03.06,Update:增加部分内容、符号统一;
    Rhie-Chow插值是C.M. Rhie和W.L. Chow在1983年对面速度进行速度插值进而引入$1-\delta$压差的一种方法$^{[1],[2]}$。本文先对棋盘型压力场进行分析,然后看Rhie-Chow如何在同位网格上解决压力棋盘分布问题,最后分析Rhie-Chow插值在OpenFOAM中的实现。

    棋盘型压力分布(Checker-board pressure field)

    ● 动量方程:考虑粘度值为$1$Pas,无时间项,线性化后的一维不可压缩动量方程: \begin{equation} u_{\mathrm{P}}^n \frac{\mathrm{d} u}{\mathrm{d}x}-\frac{\mathrm{d}^2 u}{\mathrm{d}x^2}=-\frac{\mathrm{d}p}{\mathrm{d}x} \end{equation} 其中$u$为$x$方向速度,$u_{\mathrm{P}}^n$为采用当前时间步已知的速度,$p$为压力,$_\mathrm{P}$为节点。如果我们对面速度进行中心插值,有限体积离散后有: \begin{equation} \int^e _w u_{\mathrm{P}}^n \frac{\mathrm{d} u}{\mathrm{d}x} \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} \begin{equation} \int^e _w \frac{\mathrm{d}^2 u}{\mathrm{d}x^2} \mathrm{d}x = \left(\frac{\mathrm{d}u}{\mathrm{d}x}\right)_e-\left(\frac{\mathrm{d}u}{\mathrm{d}x}\right)_w = \frac{u_{\mathrm{E}}-u_{\mathrm{P}}}{\Delta x} - \frac{u_{\mathrm{P}}-u_{\mathrm{W}}}{\Delta x} = \frac{u_{\mathrm{E}}- 2u_{\mathrm{P}}+u_{\mathrm{W}}}{\Delta x} \end{equation} \begin{equation} -\int^e _w \frac{\mathrm{d}p}{\mathrm{d}x}\mathrm{d}x = -(p_e - p_w ) \end{equation} 其中$_\mathrm{E,W}$为相邻节点,$\Delta x$为节点的距离。把公式(2),(3),(4)代入公式(1)中有下面的紧凑形式: \begin{equation} u_{\mathrm{P}}=H-(p_e-p_w) \end{equation} 同样的对压力使用线性插值有: \begin{equation} u_{\mathrm{P}}=H-\left(\frac{p_\mathrm{E}+p_\mathrm{P}}{2} - \frac{p_\mathrm{P}+p_\mathrm{W}}{2} \right) = H-\frac{p_\mathrm{E}-p_\mathrm{W}}{2} \end{equation} 其中$H$表示和临点速度有关的项。公式(4)的压力插值我们采用了二阶精度的中心差分,这导致了公式(6)引入的压力为$\frac{p_\mathrm{E}-p_\mathrm{W}}{2}$,我们称$p_\mathrm{E}-p_\mathrm{W}$为$2-\delta$压差。相反的,如果我们对公式(4)使用一阶向前或向后差分,公式(6)并不会出现$2-\delta$压差$^{[4]}$,但是压力作为驱动流体的一种驱动力,最好使用高阶精度$^{[4]}$。$2-\delta$压差不可避免。

    ● 压力泊松方程:考虑一维不可压缩流体连续性方程为: \begin{equation} \frac{\mathrm{d} u}{\mathrm{d}x} =0 \end{equation} 对其离散有: \begin{equation} u_e-u_w =0 \end{equation} 如果我们把公式(8)中的面速度采用常规的中心插值处理(对比与Rhie-Chow插值),有$u_e=\frac{u_\mathrm{E}+u_\mathrm{P}}{2}$,$u_w=\frac{u_\mathrm{W}+u_\mathrm{P}}{2}$并代入到公式(8),我们有: \begin{equation} \frac{1}{2}(u_\mathrm{E}-u_\mathrm{W})=0 \end{equation} 基于上文动量方程导出的公式(6),同样的我们有: \begin{equation} u_{\mathrm{E}}=H-\frac{p_\mathrm{EE}-p_\mathrm{P}}{2} \end{equation} \begin{equation} u_{\mathrm{W}}=H-\frac{p_\mathrm{P}-p_\mathrm{WW}}{2} \end{equation} 将公式(9),(10),(11)代入到公式(8),即为离散后的压力泊松方程: \begin{equation} \frac{1}{4}(2p_\mathrm{P}-p_\mathrm{WW}-p_\mathrm{EE}) = 0 \end{equation} 可见,即使$_\mathrm{P}$节点的压力值被引入但公式(12)中引入的也为$2-\delta$压差$^{[4]}$。早期棋盘型压力分布的解决办法是使用“错位网格”来在动量方程中引入$1-\delta$压差。更多的棋盘型分布推荐阅读相关资料$^{[1,2,3]}$。

    如何引入$1-\delta$压差

    ● 交错网格:在交错网格中,定义了针对速度、压力不同的有限体积。通过对速度在有限体积上的积分,可以把公式(1)中压力梯度项很容易的转化为$1-\delta$压差。有关错位网格如何解决棋盘压力分布请参考陶文栓的专著$^{[2]}$。

    ● 同位网格上的Rhie-Chow插值原理:Rhie-Chow插值的精髓即为不同于$u_e=\frac{u_\mathrm{E}+u_\mathrm{P}}{2}$,如果按照公式(6)的形式来对速度进行插值$^{[1,2]}$: \begin{equation} u_e=H-\frac{p_\mathrm{E}-p_\mathrm{P}}{2} \end{equation} 类似的整理代入到公式(8)有最终的: \begin{equation} \frac{1}{4}(2p_\mathrm{P}-p_\mathrm{W}-p_\mathrm{E}) = 0 \end{equation} 上述公式中成功的引入了引入$1-\delta$压差。其可在同位网格中针对压力修正方程的方法来消除棋盘型压力分布。

    OpenFOAM实现

    在OpenFOAM中,并没有直接采用Rhie-Chow插值的原始步骤。相反的,其通过对拉普拉斯项的离散巧妙的获得一种近似Rhie-Chow插值的原理。参考icoFoam解析,我们有: \begin{equation} \nabla \cdot (\mathbf{HbyA}^{n+1}) = \nabla \cdot(\frac{1}{A_{\mathrm{P},f}} \nabla p^{n+1}) \label{poss} \end{equation} 其通过求解压力泊松方程来组建离散的压力方程:

    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );
    

    正是这一步,引入了关键的$1-\delta$压差得以消除压力波。参考《OpenFOAM编程指南》中对fvm::laplacian()函数的定义。对于一维问题: \begin{equation} \mathrm{fvm::laplacian(rAU, p)} = \int_w^e \mathrm{rAU} \left( \frac{\mathrm{d}p}{\mathrm{d}x} \right) {\mathrm{d}x}= \mathrm{rAU}\left(\frac{p_{\mathrm{E}}-p_{\mathrm{P}}}{\delta x}+\frac{p_{\mathrm{W}}-p_{\mathrm{P}}}{\delta x}\right)=\mathrm{div}(\mathrm{phiHbyA}) \end{equation} 其可进一步化简为: \begin{equation} \mathrm{rAU}\frac{p_{\mathrm{E}}+p_{\mathrm{W}}-2p_{\mathrm{P}}}{\delta x}=\mathrm{div}(\mathrm{phiHbyA}) \end{equation} 这和公式(14)是相符的。可见,对于消除压力波,其中的fvm::laplacian()函数起了关键作用。如果我们只是在网格体心计算压力梯度并将其插值到面上,压力波将依然存在。

    [1]. C.M. Rhie, W.L. Chow, A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation, AIAA J. 21 (1983) 1525–1532.
    [2]. 陶文铨. 数值传热学[M]. 西安交通大学出版社, 2001.
    [3]. Versteeg H. K., Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Pearson Education, 2007.
    [4]. Ferziger J.H., Peric M. Computational methods for fluid dynamics[M]. Springer Science & Business Media, 2012.

    东岳流体
    勘误、讨论、补充内容请前往CFD中国