SIMPLE 算法全称是压力耦合方程组的半隐式方法(Semi-Implicit Method for Pressure Linked Equations), 作为经典的压力修正算法在 Fluent 等工业软件中被广泛应用。

SIMPLE 算法介绍

考虑无量纲化的二维不可压缩方程组

\[ \begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot \nabla) \mathbf{u} &= -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u}\\ \nabla \cdot \mathbf{u} &= 0\\ \end{aligned} \]

现在我们已有 \(\mathbf{u}^n\)，希望下一步的\(\mathbf{u}^{n+1},p^{n+1}\)满足

\[ \begin{aligned} \frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u})^{n+1} + \nabla p^{n+1} &= 0\\ \nabla \cdot \mathbf{u}^{n+1} &= 0 \end{aligned} \]

我们把动量方程(不等价地)拆成两部分

\[ \begin{aligned} \frac{\mathbf{u}^{*}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u})^{*} + \nabla p^{*} = 0\\ \frac{\mathbf{u}^{n+1}-\mathbf{u}^*}{\Delta t} + \nabla (p^{n+1}-p^*) &= 0 \end{aligned} \]

其中记\(p^{n+1}=p^*+p'\)，\(p^*\)是预测压强，\(p'\)是修正压强；记\(\mathbf{u}^{n+1}=\mathbf{u}^*+\mathbf{u}'\)，\(\mathbf{u}^*\)是预测速度，\(\mathbf{u}'\)是修正速度。

这里的过程与投影法相似，但是投影法只使用了预测速度，SIMPLE 算法同时使用了预测速度和预测压强。

第二个式子可以写成

\[ \frac{\mathbf{u}'}{\Delta t} + \nabla p' = 0 \]

可以理解为：在修正速度的计算中，只考虑了当前位置压强的影响，没有考虑相邻单元相互之间的影响，这是 SIMPLE 算法的假设，也是它的一大特点。

假设我们先给出了预测压强\(p^*\)，可以通过下式计算预测速度\(\mathbf{u}^*\)

\[ \frac{\mathbf{u}^*-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u})^* + \nabla p^* = 0 \]

这里对\(\mathbf{u}^*\)的求解仍然是隐式的，甚至关于\(\mathbf{u}^*\)不是线性的，离散化后得到的问题可以使用迭代法求解直到收敛，对于定常流问题，由于时间推进是无意义的，在这里的求解中也可以不用迭代至收敛，只需要保证最终迭代的解是足够准确的即可。

由于\(\mathbf{u}^{n+1}\)满足连续性方程，因此有

\[ \begin{aligned} 0 &=\nabla \cdot \mathbf{u}^{n+1} = \nabla \cdot (\mathbf{u}^* + \mathbf{u}')\\ &= \nabla \cdot \mathbf{u}^* -\Delta t \nabla^2 p' \end{aligned} \]

其中使用了\(\frac{\mathbf{u}'}{\Delta t} + \nabla p' = 0\)进行代换，得到修正压强需要满足的 Poisson 方程

\[ \nabla^2 p' = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^* \]

这里离散化后得到关于\(p'\)的线性方程组的求解问题，也可以使用迭代法求解直到收敛，对于定常流问题，同样地，可以不用迭代至收敛，只需要保证最终迭代的解是足够准确的即可。

最终得到新一步的速度和压强

\[ \begin{aligned} \mathbf{u}^{n+1}& =\mathbf{u}^*+\mathbf{u}' = \mathbf{u}^* - \Delta t \nabla p'\\ p^{n+1} &= p^* + p' \end{aligned} \]

计算流程

一个粗略的计算流程如下：

  • 假设我们已有速度\(\mathbf{u}^n\)，取一个合适的预测压强 \(p^*\) (例如直接取 \(p^*=p^n\))
  • 基于预测压强 \(p^*\)，利用下式隐式求解 \(\mathbf{u}^{*}\) ：(离散后是关于\(\mathbf{u}^*\)的非线性方程组) \[ \frac{\mathbf{u}^{*}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u})^{*} + \nabla p^{*} = 0 \]
  • 基于预测速度\(\mathbf{u}^{*}\)，利用下式求解\(p'\)：(离散后是关于\(p'\)的线性方程组) \[ \nabla^2 p' = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^{*} \]
  • 更新速度和压强 \[ \begin{aligned} \mathbf{u}^{n+1}& =\mathbf{u}^*+\mathbf{u}' = \mathbf{u}^* - \Delta t \nabla p'\\ p^{n+1} &= p^* + p' \end{aligned} \]

定常流

对于定常流，迭代对应的时间推进是没有意义的，因此我们在两个离散求解过程中都可以很大程度地偷懒，只需要通过时间推进，最终达到需要的精度即可，尤其是求解\(\mathbf{u}^*\)涉及的非线性方程组，它的原始形式是

\[ \begin{aligned} \frac{\mathbf{u}^{*}-\mathbf{u}^n}{\Delta t} + (\mathbf{u}\cdot\nabla \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u})^{*} + \nabla p^{*} = 0\\ \end{aligned} \]

可以把全隐改成半隐或者显式计算\(\mathbf{u}^*\)：

  • 可以把对流项中的一个用\(\mathbf{u}^n\)代替，转化为关于\(\mathbf{u}^*\)的线性方程组求解： \[ \frac{\mathbf{u}^*-\mathbf{u}^n}{\Delta t} + (\mathbf{u}^n\cdot\nabla \mathbf{u}^* - \frac{1}{Re}\nabla^2 \mathbf{u}^*) + \nabla p^* = 0 \]
  • 可以把时间项和对流项中的一个用\(\mathbf{u}^n\)代替，相当于没有取时间项的离散方程组： \[ (\mathbf{u}^n\cdot\nabla \mathbf{u}^* - \frac{1}{Re}\nabla^2 \mathbf{u}^*) + \nabla p^* = 0 \]
  • 甚至可以把对流项和黏性项完全使用\(\mathbf{u}^n\)代替，显式计算\(\mathbf{u}^*\)： \[ \frac{\mathbf{u}^*-\mathbf{u}^n}{\Delta t} + (\mathbf{u}^n\cdot\nabla \mathbf{u}^n - \frac{1}{Re}\nabla^2 \mathbf{u}^n) + \nabla p^* = 0 \]

对于压强修正值的求解，是一个 Poisson 方程

\[ \nabla^2 p' = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^* \]

离散得到的线性方程组，也可以通过 G-S 迭代等进行固定步数的迭代，近似求解。

非定常流

对于非定常流，迭代的每一步都是有意义的时间推进，两个离散求解过程都需要保持足够的精度，这里我们可以在预测速度和修正压强的隐式求解中分别进行内迭代，达到精确求解的目的，例如一种可选的计算流程如下：

  • 假设我们已有速度\(\mathbf{u}^n\)，取一个合适的预测压强 \(p^{*}\) (例如直接取 \(p^{*}=p^n\))
  • 基于预测压强\(p^{*}\)，进入\(\mathbf{u}^*\)求解的内迭代：
    • 获取合适的预测速度启动值\(\mathbf{u}^{*,0}\)；
    • 迭代计算预测速度：通过\(\mathbf{u}^{*,k}\)计算\(\mathbf{u}^{*,k+1}\)，直至收敛；
    • 取\(\mathbf{u}^* = \mathbf{u}^{*,k_M}\)结束内迭代。
  • 基于预测速度\(\mathbf{u}^*\)，进入\(p'\)求解的内迭代：
    • 获取合适的修正压强启动值\(p^{',0}\)；
    • 迭代计算预测速度：通过\(p^{',k}\)计算\(p^{',k+1}，\)直至收敛；
    • 取\(p' = p^{',k_N}\)结束内迭代。
  • 更新速度和压强 \[ \begin{aligned} \mathbf{u}^{n+1}& =\mathbf{u}^*+\mathbf{u}' = \mathbf{u}^* - \Delta t \nabla p'\\ p^{n+1} &= p^* + p' \end{aligned} \]

更高效的作法是把两个内迭代过程合在一起，组成一个内迭代过程，计算流程如下：

  • 假设我们已有速度\(\mathbf{u}^n\)；
  • 进入内迭代，启动值\(\mathbf{u}^{n+1,0},p^{n+1,0}\)
    • 获取预测压强\(p^{*,k}=p^{n+1,k-1}\)；
    • 基于预测压强\(p^{*,k}\)和\(\mathbf{u}^{*,k-1}\)，近似求解\(\mathbf{u}^{*,k}\)；
    • 基于预测速度\(\mathbf{u}^{*,k}\)和\(p^{',k-1}\)，近似求解\(p^{',k}\)；
    • 更新 \[ \begin{aligned} \mathbf{u}^{n+1,k} &=\mathbf{u}^{*,k}+\mathbf{u}^{',k} = \mathbf{u}^{*,k} - \Delta t \nabla p^{',k}\\ p^{n+1,k} &= p^{*,k} + p^{',k}\\ \end{aligned} \]
    • 直到 \(\mathbf{u}^{n+1,k_N},p^{n+1,k_N}\) 满足收敛条件，退出内迭代；
  • 更新速度和压强 \[ \begin{aligned} \mathbf{u}^{n+1}& = \mathbf{u}^{n+1,k_N}\\ p^{n+1} &= p^{n+1,k_N} \end{aligned} \]

其他

强烈推荐李新亮的计算流体力学视频，讲得 SIMPLE 比网上的其他各种教程都透彻详细。 SIMPLE 算法的改进包括 SIMPLEC 算法，PISO 算法等，它们的改进主要是在迭代过程，以及系数上的小改动，以期望更快地达到收敛。