Code of the Week #4: solve(UEqn == -fvc::grad(p))

今天，我们要讲的是

solve(UEqn == -fvc::grad(p));

这一行程序来自于 icoFoam.C ，从字面上理解是求解方程组。

首先，来看一下括号里面各项的定义：

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

以类 fvVectorMatrix 定义了一个以 UEqn 为名的变量，括号内的各项已经分别在 code of the week #1: $\nabla(UU)$ div(phi, U), #2: $\dfrac{\partial U}{\partial t}$ ddt(U), #3 $\nabla\cdot(\nu\nabla U)$ laplacian(nu, U) 分别进行了解读，其分别代表了方程等号左边的各项 $$ \begin{equation} \frac{\partial U}{\partial t} + \nabla\cdot(UU) + \nabla\cdot(\nu\nabla U) = -\nabla p \end{equation} \tag{1} $$

对方程左边部分离散化后可以得到以下矩阵方程。 $$ \begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N} \mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=0\end{equation} \tag{2} $$

其中 $U^r_P$ 指待求网格的速度， $U^n_N$ 指相邻网格已知速度值 (前一时间步)， $E_P^n$ 为该点已知速度计算的值。这一方程加上 $\nabla p$ 就成为动量预测方程。 $$ \begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N} \mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=-\nabla p \end{equation} \tag{3} $$ 可简写为 $$ A U = B $$

该行程序双等号右边的 -fvc::grad(p) 是动量方程右边的部分，而使用 fvc:: 表示利用前一时间步的值来获得此梯度。

这一行就是求解 U 方程获取预测速度，不能满足连续性方程。因此，需要进行修正。

如果我们定义 $\mathbf U_\mathrm{P}^{n+1}$ 为 PISO 循环后最终收敛的速度，那么，进行修正后的 $\mathbf U_\mathrm{P}^{n+1}$ 就应当满足动量方程和连续性方程。

速度 $\mathbf U_\mathrm{P}^{n+1}$ 符合动量方程

$$ \begin{equation} 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} \end{equation} \tag{4} $$

其中， $p^{n+1}$ 为 PISO 循环内最后求得的压力。

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

式中， $\mathbf{U}^{*,n+1}$ 为预测速度值。要使可求得的速度满足连续性方程，只需要对 $\mathbf U_\mathrm{P}^{n+1} $ 进行求散度操作。 $$ \begin{equation} \nabla \cdot \mathbf U_\mathrm{P}^{n+1}=0 \end{equation}\tag{6} $$ 离散后得

$$ \begin{equation} \sum \left(\mathbf U_\mathrm{P}^{n+1} \cdot \mathbf S_f \right)=0 \end{equation}\tag{7} $$ 而此时，该网格面上的中心值也未求得，因此，需要另辟他法，Rhie-Chow 插值建议这样进行获取面上的 $\mathbf U_{\mathrm{P,}f}^{n+1}$：

$$ \begin{equation} \mathbf U_{\mathrm{P},f}^{n+1} = \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}} \nabla_f p^{n+1} \end{equation}\tag{8} $$

将获取的面心上的速度 $\mathbf U_{\mathrm{P,}f}^{n+1}$ 代入连续方程，

$$ \begin{equation} \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 \end{equation}\tag{9} $$

写成张量形式，有

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

PISO 循环

方程的解 $\mathbf U_\mathrm{P}^{n+1}$ 和压力 $p^{n+1}$ 应当满足连续性方程和动量方程。而到目前为止，仅有通过动量预测方程求出来的 $\mathbf U_\mathrm{P}^r$。动量预测方程 和 替代速度连续性方程 无法自行求解。Issa 于 1986 年提出 PISO (Pressure Implicit with Splitting of Operators) 技术来解决这个问题，通过对当前时间步内的多次修正来获得最终的 $\mathbf U_\mathrm{P}^{n+1}$ 和 $p^{n+1}$ 。

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

1. 根据初始条件（压力场、速度场等）求解预测速度场 $\mathbf U^r$

2. 根据预测速度或第 4 步求得的速度，计算出无压力梯度项的速度场 $\mathbf U^n$

3. 根据无压力梯度项的速度场 $\mathbf U^n$，求解压力（泊松）方程 (9) 或 (10)，得到压力场 $p^n$

4. 根据压力场 $p^n$，使用方程 (8) 修正预测速度 $\mathbf U^{n+1}$ 以满足连续性方程，获得一个新的速度值

5. 接第 2 步

社区讨论请移步 Code of the Week #4: solve(UEqn == -fvc::grad(p))