Code of the Week #7: fvm::laplacian(rUA, p) == fvc::div(phi)

如需讨论，请移步至 CFDwired Forum

fvm::laplacian(rUA, p) == fvc::div(phi)

这一行程序是 icoFoam 中的一行，写成张量形式为 $$ \frac{1}{A_p}\nabla^2 p =\frac{1}{A_p}\nabla\cdot(\nabla p)= \nabla \cdot U \tag{1} $$ 称为压力泊松方程。rUA 就是方程中的 $\dfrac{1}{A_p}$ 。

$A_p$ 是动量方程离散后，去掉压力梯度项，当前 Cell 的矩阵系数： $$ A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = 0 \tag{2} $$

不过，我们先来看一下 U 方程：

下面 icoFoam.C 片段：

        fvVectorMatrix UEqn       // 速度方程
        (
            fvm::ddt(U)           // 非稳态项，隐式离散
          + fvm::div(phi, U)      // 对流项，隐式离散
          - fvm::laplacian(nu, U) // 扩散项，隐式离散
        );

此为无压力梯度项的速度方程，张量形式为 $$ \text{UEqn} = \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})- \nabla \cdot(\nu \nabla \mathbf{U}) \tag{3} $$ 该方程在后台就形成离散形式 (2)。

速度方程的张量形式为 $$ \text{UEqn} =-\nabla \frac{p}{\rho}\tag{4} $$ 求解速度方程，得预测速度 $U^*$

        solve(UEqn == -fvc::grad(p)); // 不可压缩时, p=压力/密度

此时求得的预测速度 $U^*$ 并不能满足连续性方程。因此，需要一点特殊的方法来修正。假设修正后的速度为 $\mathbf U^{n+1}$

$$ A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = - \nabla p^{n+1} \tag{5} $$

修正前的速度 (去掉压力梯度项) 为

$$ \mathbf{U}^* = \mathbf U_\mathrm{P}^r = \frac{1}{A_\mathrm{P}}\left( - \sum A_\mathrm{N} \mathbf U_\mathrm{N}^r + E_\mathrm{P}^r \right) \tag{6} $$

方程 $(5)-A_\mathrm{P}(6)$ 为

$$ \mathbf U_\mathrm{P}^{n+1} -\mathbf{U}^{*,n+1}= - \frac{1}{A_\mathrm{P}}\nabla p^{n+1} \tag{7} $$

稍作变换

$$ \mathbf U_\mathrm{P}^{n+1} = \mathbf{U}^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla p^{n+1} \tag{7a} $$

面插值版本为

$$ \mathbf U_f^{n+1} = \mathbf U_f^{*,n+1} - \frac{1}{A_{\mathrm P}}\nabla_{\! f} \,p^{n+1} \tag{7b} $$

而速度 $\mathbf U_\mathrm{P}^{n+1}$ 符合连续性方程

$$ \nabla \cdot \mathbf U_\mathrm{P}^{n+1}=0 \tag{8} $$

将 (8) 代入 (7) 的散度，有

$$ \nabla\cdot\mathbf{U}^{*,n+1}= \nabla\cdot\left(\frac{1}{A_\mathrm{P}}\nabla p^{n+1}\right) \tag{9} $$

这就是压力泊松方程。将 (9) 离散，得

$$ \sum \left[ \left( \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}}\nabla_f p^{n+1} \right) \cdot \mathbf S_f \right] = 0 \tag{10} $$

求解该方程将用到修正前的速度 $U^*$ ，从而可求解 $p$ ，然后将 $p$ 用于修正速度场 (方程 (7a))。

以下为压力泊松方程的循环。

...     // 前面已经完成预测速度 Ｕ* 的计算

        for (int corr = 0; corr < nCorr; corr++)
        {
            volScalarField rUA = 1.0/UEqn.A();      // 取得 U 方程的 1/Ap

            U = rUA*UEqn.H();                       // 无压力梯度项之速度值 U'
            phi = (fvc::interpolate(U) & mesh.Sf()) // U' 使用线性插值，需要在算例中设定interpolate(U)的格式
                + fvc::ddtPhiCorr(rUA, U, phi);     // 在 ext-4.0 版本中，这一项已经去掉

            adjustPhi(phi, U, p);            // 修正压力边界（梯度为0）上的 phi，使其满足连续性方程
           //网格非正交循环
            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)  // 压力泊松方程，此处间接使用 Rhie-Chow 插值原理
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();               // 求解压力泊松方程

                if (nonOrth == nNonOrthCorr)
                {
                    phi -= pEqn.flux();   // 方程 (7b)，使用求解的压力场修正通量场，即在最后一次修正的时候使通量守恒，
                }// pEqn.flux()返回方程 (7) RHS，为 fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf())
            }    // 某些可压缩求解器中, 修正项为 phi += pEqn.flux(), 这是因为 pEqn 中的 laplacian 项为 '-' 号

#           include "continuityErrs.H"

            U -= rUA*fvc::grad(p);   // 公式 (7a) 
            U.correctBoundaryConditions();
        }

icoFoam 使用 PISO 算法解 N-S 方程。 其基本程序为

1. 根据初始条件（压力场、速度场等）求解预测速度场 $U^*$
2. 根据预测速度求无压力梯度项的速度场 $U'$
3. 根据无压力梯度项的速度场 $U'$ ，求解压力（泊松）方程，得到压力场 $p$
4. 根据压力场 $p$ ，修正预测速度 $U^*$ 以满足连续性方程
5. 接第 2 步，直到满足收敛要求

如需讨论，请移步至 CFDwired Forum