标记网格法(Marker-and-Cell method，MAC)，是一种求解压力 Poisson 方程算法，最初被设计用于求解含自由面的不可压缩流动，下文的算法并不涉及自由面，只是关注对不可压缩条件的处理。

求解压力 Poisson 方程算法

求解压力 Poisson 方程算法是一大类 N-S 方程的求解算法，它的特点是：通过求解一个 Poisson 方程获取压强\(p\)，并且使得速度散度为零间接地得到满足。(张德良，计算流体力学教程)

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

\[ \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} \]

记\(D = \nabla \cdot \mathbf{u}\)，对动量方程求散度，得到

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

其中记\(S_p\)如下

\[ S_p = - \nabla \cdot \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right] \]

从而压强满足

\[ \nabla^2 p = S_p - D_t \]

通过如上的 Poisson 方程求解\(p\)，就是此类方法的主要特点。求解压力 Poisson 方程算法主要包括以下几种：

  • MAC 方法，SOLA 方法(相当于 MAC 方法的迎风版)
  • 投影法(又称为时间分裂法，分步法)

本文中我们主要介绍 MAC 方法，但是严格来说，投影法也属于求解压力 Poisson 方程算法的一种，并且和 MAC 方法相似。

MAC 方法

MAC 方法又称为标记网格法(Marker-and-Cell method)，在求解自由面的不可压缩流动时，使用了 Euler-Lagrange 混合方法，下文没有涉及自由面，因此没有体现处这一特点。MAC 方法是最早提出使用交错网格的方法，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} \]

首先，时间离散得到下式

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

直接对动量方程取散度，得到

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

这里对两个时间层的速度散度，记\(D^n = \nabla \cdot \mathbf{u}^n\)，同理 $ D^{n+1}$ ，

  • \(D^{n+1}=\nabla \cdot \mathbf{u}^{n+1}=0\)是方程组的要求，强制成立，因此直接消去；
  • \(D^{n}=\nabla \cdot \mathbf{u}^{n}\) 是上一步的求解结果，在求解中不能严格成立，因此我们允许这种余量的存在，并体现在下一步计算压强的过程中。

计算流程

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

  • 已知第 n 步的速度\(\mathbf{u}^n\)；
  • 代入下式，计算第 n+1 步的压强\(p^{n+1}\) \[ -\frac{\nabla \cdot \mathbf{u}^n}{\Delta t} + \nabla \cdot \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right]^{n} = -\nabla^2 p^{n+1} \]
  • 代入下式，计算第 n+1 步的速度\(\mathbf{u}^{n+1}\) \[ \frac{\mathbf{u}^{n+1}-\mathbf{u}^n}{\Delta t} + \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right]^{n} = -\nabla p^{n+1} \]

注意到，上述单步流程实际上不能保证\(D^{n+1}=0\)严格成立，只是在计算压强时利用了\(D^{n+1}=0\)，因此可以使用 \(\max(D^{n+1})\le \epsilon\) 作为内迭代终止的判断条件；在计算压强时也可以使用迭代法，例如使用 G-S 迭代法，SOR 迭代法等，有以下几种典型的迭代计算过程。

第一种计算流程：

  • 已有第 n 步的速度场\(\mathbf{u}^n\)；
  • 进入求解内迭代：
    • 取\(\mathbf{u}^{n+1,0} = \mathbf{u}^n\)，取\(p^{n+1,0}=p^n\)
    • 利用 G-S 迭代法近似求解泊松方程，基于\(p^{n+1，k}\)和\(\mathbf{u}^{n+1,k}\)，计算\(p^{n+1,k+1}\)； \[ \frac{\nabla \cdot \mathbf{u}^{n+1,k} - \nabla \cdot \mathbf{u}^n}{\Delta t} + \nabla \cdot \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right]^{n} = -\nabla^2 p^{n+1,k+1} \]
    • 利用动量方程，基于\(p^{n+1,k+1}\)计算\(\mathbf{u}^{n+1,k+1}\)； \[ \frac{\mathbf{u}^{n+1,k+1}-\mathbf{u}^n}{\Delta t} + \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right]^{n} = -\nabla p^{n+1,k+1} \]
    • 计算\(D^{n+1,k+1}= \nabla \cdot \mathbf{u}^{n+1,k+1}\)，判断\(\max(D^{n+1,k+1})\le \epsilon\)是否成立？
    • 假设\(k=k_N\)时满足，退出内迭代；
  • 取\(\mathbf{u}^{n+1}=\mathbf{u}^{n+1,k_N},p^{n+1}=p^{n+1,k_N}\)。

第二种计算流程(Chorin 迭代法)： 不使用 Poisson 方程进行迭代法求解，而是对压强进行直接的修正。

  • 已有第 n 步的速度场\(\mathbf{u}^n\)；
  • 进入求解内迭代：
    • 取\(\mathbf{u}^{n+1,0} = \mathbf{u}^n\)，取\(p^{n+1,0}=p^n\)
    • 计算\(D^{n+1,k} = \nabla \cdot \mathbf{u}^{n+1,k}\)，通过\(p^{n+1，k}\)直接计算\(p^{n+1,k+1}\)，其中\(\lambda\)是可选的松弛因子； \[ p^{n+1,k+1} = p^{n+1,k} - \lambda D^{n+1,k} \]
    • 利用动量方程，基于\(p^{n+1,k+1}\)计算\(\mathbf{u}^{n+1,k+1}\)； \[ \frac{\mathbf{u}^{n+1,k+1}-\mathbf{u}^n}{\Delta t} + \left[(\mathbf{u}\cdot \nabla) \mathbf{u} - \frac{1}{Re}\nabla^2 \mathbf{u}\right]^{n} = -\nabla p^{n+1,k+1} \]
    • 计算\(D^{n+1,k+1}= \nabla \cdot \mathbf{u}^{n+1,k+1}\)，判断\(\max(D^{n+1,k+1})\le \epsilon\)是否成立？
    • 假设\(k=k_N\)时满足，退出内迭代；
  • 取\(\mathbf{u}^{n+1}=\mathbf{u}^{n+1,k_N},p^{n+1}=p^{n+1,k_N}\)。

要求\(\lambda\)满足

\[ \lambda \le \lambda_{\text{max}} = \frac{1}{\Delta t \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} \right)} \]

例如可以取\(\lambda = \frac12 \lambda_{\text{max}}\)。