1 N-S 方程

无量纲化的不可压 N-S 方程组如下：

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

下文将不可压 N-S 方程组简写为 INSE。与可压缩流动相比，不可压缩流动有如下特点:

  • INSE 的压强：只出现在动量方程中，从数学上只是起到了一个 Lagrange 乘子的作用，失去了热力学意义；
  • INSE 不含压强的时间导数项 \(p_t\) ，压强\(p\)的显式计算比较困难，往往需要求解一个全局的 Poisson 方程；
  • 从数学上，方程组呈现抛物-椭圆型的特点，而非可压缩情形的双曲型；
  • 从物理意义上，不可压缩条件意味着声速无穷大，在一点的扰动会瞬间传遍整个区域，从而边界的小误差也会造成很大的影响。

对于不可压方程的处理中，对流项和扩散项分别发挥着不同的角色：

  • 在高雷诺数时，对流占优，扩散项的影响通常较小(除了边界层等区域)；(双曲型)
  • 在低雷诺数时，扩散占优，此时为了避免时间步长过小，对扩散项的离散必须使用隐式。(抛物-椭圆型)

在雷诺数极小或较小的定常流情形下，不可压流体还有如下的简化模型：

  1. Stokes 方程组： \[ \begin{aligned} -\nu \nabla^2 \mathbf{u} + \nabla p &= \mathbf{f}\\ \nabla \cdot \mathbf{u} &= 0 \end{aligned} \]
  2. Oseen 方程组：（速度 \(\mathbf{w}\)已知） \[ \begin{aligned} -\nu \nabla^2 \mathbf{u} + (\mathbf{w}\cdot \nabla)\mathbf{u} + \nabla p &= \mathbf{f}\\ \nabla \cdot \mathbf{u} &= 0 \end{aligned} \]

2 常见算法

对于 INSE 的求解，可以分成两类算法：

  • 耦合式算法(coupled method)：所有变量联立求解，方法直观，但是离散得到的代数方程组具有很强的刚性，计算效率低；
  • 分离式算法(segregated method)：逐个变量进行求解。

从求解变量的角度，算法又可以分成如下两类：

  • 原始变量法：求解变量为速度\(\mathbf{u}\)和压强\(p\)，它们具有明确的物理含义，边界处理较容易；
  • 非原始变量法：求解变量为导出变量，如使用导出变量流函数\(\phi\)和涡量\(w\)，对应为涡量—流函数方法，还有流函数方法、涡量-速度方法等。它们面临的方程组形式上比原始方程组更容易处理，但是导出变量的边界处理更加困难。

2.1 非原始变量法

涡量-流函数方法

通过变量代换，可以把二维不可压缩 N-S 方程组变换为关于涡量和流函数的方程组，进而求解。

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

其中速度\(\mathbf{u}=(u,v)^T\)和流函数\(\phi\)，涡量\(w\)的关系为：

\[ \begin{aligned} & u = \phi_y \\ &v = -\phi_x \\ &w = u_y - v_x \end{aligned} \]

详细内容见后续笔记。

对于涡量—流函数方法为代表的非原始变量法，优缺点都非常明显：

  • 优点：消去了压强项，不可压缩条件自动成立，对方程组的处理更加容易；
  • 缺点：导出变量的边界条件难以确定，不易从二维推广到三维。

2.2 原始变量法

常见的使用原始变量的分离式算法包括以下几大类：

  • 人工压缩性方法(ACM)
  • 求解压力 Poisson 方程算法
    • MAC 方法
    • 投影法
  • 压力修正算法
    • SIMPLE 算法

下文是对几个主流算法的简要介绍。

人工压缩性方法(ACM)

将连续性方程替换为下式进行求解

\[ \frac{\partial p}{\partial t} + c^2 \nabla \cdot \mathbf{u} =0 \]

对于定常流的求解，在时间推进到收敛时，压强场不变，速度场仍然处处满足\(\nabla \cdot\mathbf{u}=0\)；对于非定常流的求解，则需要加入内迭代，计算效率低，详细内容见后续笔记。

投影法

最早的无压力增量投影法如下，引入预测速度\(\mathbf{u}^*\)，将动量方程拆分为两个式子

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

对第二个式子取散度得到

\[ \begin{align} \nabla^2 p^{n+1} = \frac{1}{\Delta t} \nabla \cdot \mathbf{u}^* \end{align} \]

计算流程即：

  • 通过(1)式计算\(\mathbf{u}^*\)；
  • 通过(2)式计算\(p^{n+1}\)；
  • 通过(3)式计算\(\mathbf{u}^{n+1}\)。

最后一步相当于把 \(\mathbf{u}^*\) 向散度为零的子空间做投影得到 \(\mathbf{u}^{n+1}=P(\mathbf{u}^*)\)，这也是此类方法称为投影法的原因。

后续的改进包括考虑压力增量，以及速度的旋度的修正，引入中间层等，详细内容见后续笔记。

SIMPLE 算法

SIMPLE 算法全称是压力耦合方程组的半隐式方法 (Semi-Implicit Method for Pressure Linked Equations)，是一种经典的压力修正算法。

SIMPLE 算法在商业软件 Fluent 中被广泛使用，在 Fluent 中默认提供了四种算法选择，包括三种分离式算法和一种耦合式算法：

  • SIMPLE 算法
  • SIMPLEC 算法
  • PISO 算法
  • Coupled 算法(耦合式算法)

其中 SIMPLEC 和 PISO 都是基于 SIMPLE 算法进行的改进。Fluent 建议，在非定常流动(瞬态流动)问题中选择 PISO 算法，在定常流动问题中选择 SIMPLE 或 SIMPLEC 算法；在高速可压缩流动中使用耦合式算法。

SIMPLE 算法对每一步计算速度和压强都使用了预测值和修正值，一个粗略的计算流程如下：

  • 假设我们已有速度\(\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}^{*} \]
  • 更新速度\(\mathbf{u}^{n+1}\)和压强 \(p^{n+1}\) \[ \begin{aligned} \mathbf{u}^{n+1}& =\mathbf{u}^*+\mathbf{u}' = \mathbf{u}^* - \Delta t \nabla p'\\ p^{n+1} &= p^* + p' \end{aligned} \]

另见后续笔记。

3 空间离散——交错网格

在使用有限差分/有限体积法的空间离散中，往往会采用交错网格来避免压力场的奇偶失联：因为在通常的均匀矩形网格中，对压强使用中心差分，如下高度振荡的压强场，压强梯度的计算却全都是零:

\[ P_{IJ} = \left\{ \begin{aligned} &P_1 & I+J\text{为奇数}\\ &P_2 & I+J\text{为偶数}\\ \end{aligned} \right. \]

这就是奇偶失联现象，交错网格的设置可以有效解决这一问题： 对于有限差分法，把标量(压强，温度等)定义在整格点例如\((I,J)\)，把速度向量的 u 分量定义在水平方向的半格点例如\((I+1/2,J)\)，把 v 分量定义在竖直方向的半格点例如\((I,J+1/2)\)。对于有限体积法把相应的控制体定义为以上述整格点或半格点为中心的矩形，如下图。