第 6 章：一些简单的 CFD 技术入门

本章在解决什么问题

前面已经完成：

控制方程→PDE 数学性质→离散化基本理论→网格\text{控制方程} \rightarrow \text{PDE 数学性质} \rightarrow \text{离散化基本理论} \rightarrow \text{网格}

现在的问题变成：

给定一个具体 PDE，到底怎样真正把离散方程一步一步算出来？

这一章介绍的各种方法并不是彼此孤立的算法。

它们主要分别解决：

  • 双曲型传播问题；
  • 抛物型扩散问题；
  • 椭圆型边值问题；
  • 多维隐式计算；
  • 不可压缩流压力—速度耦合；
  • 高频振荡与稳定性。

时间推进与空间推进

大多数非定常 CFD 问题采用：

tn→tn+1t^n \rightarrow t^{n+1}

进行时间推进。

但某些稳态双曲型问题，也可以把一个空间方向当成类似“时间”的推进方向，例如超声速流。

教材中某些方法会同时被描述为：

  • time marching；
  • space marching。

这与 PDE 的信息传播方向直接相关。

Lax-Wendroff 方法

简单例子

考虑某个时间项的 Taylor 展开：

un+1=un+Δt ut+Δt22utt+⋯u^{n+1} = u^n + \Delta t\,u_t + \frac{\Delta t^2}{2}u_{tt} +\cdots

然后利用 PDE 本身，把：

ut,utt,…u_t,\qquad u_{tt},\dots

转换成空间导数。

以线性对流为例：

ut+aux=0u_t+a u_x=0

得出：

utt=a2uxx.u_{tt}=a^2u_{xx}.

代回 Taylor 展开以后，就得到一个二阶时间精度的显式格式。

主要特征

Lax-Wendroff 的目标是：

利用原本的 PDE 构造和消去时间高阶导数，从而构造高阶显式推进

优点：

  • 二阶精度；
  • 结构清晰；
  • 对光滑波传播性能较好。

缺点：

  • 在激波和间断附近容易产生非物理振荡；
  • 对 CFL 条件敏感；
  • 对复杂非线性方程直接求高阶导数比较麻烦。

为什么会出现振荡

Lax-Wendroff 的主导误差往往具有明显的色散性质。

这意味着不同波数具有不同的数值传播速度。

因此在陡峭波前附近可能出现波前前后非物理振铃

这也是后来需要：

  • artificial viscosity；
  • flux limiter；
  • TVD；
  • ENO/WENO

等方法的重要原因。

MacCormack 方法

相比 Lax-Wendroff 的优化

Lax-Wendroff 直接对复杂非线性系统反复求：

utt,uttt,…u_{tt},u_{ttt},\dots

会很麻烦，书上仅在二维的系统上进行推导就进行了大量的反复求导、方程之间的消去等操作。

MacCormack 使用预测—校正思想避免显式计算这些高阶导数。

基本结构

第一步：预测。

例如使用前向空间差分：

Ui∗=Uin−ΔtΔx(Fi+1n−Fin).U_i^* = U_i^n - \frac{\Delta t}{\Delta x} \left( F_{i+1}^n-F_i^n \right).

第二步：校正。

使用反向空间差分：

Uin+1=12[Uin+Ui∗−ΔtΔx(Fi∗−Fi−1∗)].U_i^{n+1} = \frac12 \left[ U_i^n+U_i^* - \frac{\Delta t}{\Delta x} \left( F_i^*-F_{i-1}^* \right) \right].

也可以把前向、后向顺序反过来，也可以在不同的时间步之间交替使用，均为二阶精度。

MacCormack 与 Lax-Wendroff 的关系

两者对于简单线性方程具有非常接近的二阶行为。

但 MacCormack：

  • 编程更直接；
  • 不需要显式写出高阶时间导数；
  • 更容易推广到 Euler/Navier-Stokes 方程组。

因此曾经是经典 CFD 中非常常见的方法。

缺点仍然是：

在激波附近容易出现色散振荡。

所以通常需要额外 人工耗散。

Crank-Nicolson 方法

基本思想

对于：

ut=L(u)u_t=L(u)

Crank-Nicolson 不只使用 nn 时刻，也不只使用 n+1n+1 时刻，而是取两者平均：

un+1−unΔt=12[L(un)+L(un+1)]\boxed{ \frac{u^{n+1}-u^n}{\Delta t} = \frac12 \left[ L(u^n)+L(u^{n+1}) \right] }

可以理解为在时间方向使用梯形积分。

主要特点

Crank-Nicolson：

  • 时间二阶精度；
  • 对扩散方程具有很好的稳定性；
  • 可以使用比显式扩散方法大得多的时间步。

但：

稳定并不等于任何时间步都准确。

即使理论稳定，时间步过大仍然会造成严重相位和振幅误差。

对非常刚性的高频模态，Crank-Nicolson 还可能出现轻微非物理时间振荡。

松弛方法

主要解决的问题

Laplace/Poisson 等椭圆型问题：

∇2u=f\nabla^2u=f

通常不能像双曲型方程那样沿某个方向直接推进。

因为整个区域是全局耦合的，需要迭代求解整个离散代数系统。

基本思想

例如二维 Laplace：

ui,j=14(ui+1,j+ui−1,j+ui,j+1+ui,j−1).u_{i,j} = \frac14 \left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} \right).

可以不断利用邻居更新当前节点：

u(0)→u(1)→u(2)→…u^{(0)} \rightarrow u^{(1)} \rightarrow u^{(2)} \rightarrow\dots

直到残差足够小。

Jacobi、Gauss-Seidel 与 SOR

Jacobi：一次迭代全部使用上一轮数据。

Gauss-Seidel：已经更新的新数据立即投入本轮后续计算。通常收敛更快。

SOR 则进一步引入松弛系数：

unew=uold+ω(uGS−uold).u^{new} = u^{old} + \omega \left( u^{GS}-u^{old} \right).

1<ω<21<\omega<2 时常可以加速收敛，称为超松弛。ω<1\omega<1 称作亚松驰。

ADI 方法

ADI：

Alternating Direction Implicit\boxed{ \text{Alternating Direction Implicit} }

即交替方向隐式法。

为什么需要 ADI

二维隐式方程往往会让：

(i,j)(i,j)

节点和上下左右节点全部耦合。

直接求解完整二维矩阵代价较高。

ADI 的想法是：

把一次二维隐式计算拆成两个一维隐式计算。

一维隐式计算时，系数阵为 3 对角阵，可以更简便地用Thomas算法解决。

一个典型过程

第一半步：

x 方向隐式,y 方向显式.x\text{ 方向隐式}, \qquad y\text{ 方向显式}.

第二半步：

y 方向隐式,x 方向显式.y\text{ 方向隐式}, \qquad x\text{ 方向显式}.

这样每一次只需求解许多条一维三对角系统。

因此：

二维大系统→大量一维三对角系统\boxed{ \text{二维大系统} \rightarrow \text{大量一维三对角系统} }

计算成本显著下降。

实现代价

方向分裂不是完全免费的。

拆分：

Lx+LyL_x+L_y

会引入一定的 splitting error。

但对于许多扩散型和弱耦合问题，其效率优势非常明显。

压力修正技术

不可压缩 NS 问题

不可压缩方程：

∇⋅u=0.\nabla\cdot\mathbf u=0.

动量方程：

∂u∂t+u⋅∇u=−1ρ∇p+ν∇2u.\frac{\partial\mathbf u}{\partial t} + \mathbf u\cdot\nabla\mathbf u = -\frac1\rho\nabla p + \nu\nabla^2\mathbf u.

这里压力并没有类似：

pt=⋯p_t=\cdots

这样的独立演化方程。

因此：

压力的作用主要是强制速度场满足不可压缩约束。

所以压力更像一个约束变量。

预测—压力修正—校正

一个典型 projection / MAC 思路是：

先忽略新压力预测：

u∗\mathbf u^*

然后令：

un+1=u∗−Δtρ∇pn+1.\mathbf u^{n+1} = \mathbf u^* - \frac{\Delta t}{\rho} \nabla p^{n+1}.

要求：

∇⋅un+1=0.\nabla\cdot\mathbf u^{n+1}=0.

得到：

∇2pn+1=ρΔt∇⋅u∗\boxed{ \nabla^2p^{n+1} = \frac{\rho}{\Delta t} \nabla\cdot\mathbf u^* }

所以一个时间步为：

预测速度→压力 Poisson→压力修正速度\boxed{ \text{预测速度} \rightarrow \text{压力 Poisson} \rightarrow \text{压力修正速度} }

最后速度重新满足：

∇⋅un+1=0.\nabla\cdot\mathbf u^{n+1}=0.

交错网格与不可压缩流

经典交错网格中：

压力：

pi,jp_{i,j}

放在单元中心。

横向速度：

ui+1/2,ju_{i+1/2,j}

放在左右面。

纵向速度：

vi,j+1/2v_{i,j+1/2}

放在上下面。

于是连续方程：

ui+1/2,j−ui−1/2,jΔx+vi,j+1/2−vi,j−1/2Δy=0\frac{ u_{i+1/2,j} - u_{i-1/2,j} }{\Delta x} + \frac{ v_{i,j+1/2} - v_{i,j-1/2} }{\Delta y} =0

正好就是控制体四个面的净流量为零

而压力梯度直接作用在面速度上：

pi+1,j−pi,jΔx.\frac{p_{i+1,j}-p_{i,j}}{\Delta x}.

这构成：

pcell→∇uface→∇⋅cell\boxed{ p_{\text{cell}} \xrightarrow{\nabla} u_{\text{face}} \xrightarrow{\nabla\cdot} \text{cell} }

非常自然的离散结构。

数值耗散、色散，人工粘性

一个差分格式不仅可以被看成，原 PDE 的近似。

也可以反过来理解成，某个稍微被修改过的 PDE 的精确离散。

例如一阶迎风离散线性对流：

ut+aux=0u_t+a u_x=0

通过 Taylor 展开以后，可得到类似：

ut+aux=C2uxx+C3uxxx+⋯u_t+a u_x = C_2u_{xx} + C_3u_{xxx} +\cdots

右侧这些额外项完全由数值离散产生。

数值耗散

修正方程中的偶数阶导数，例如：

uxx,uxxxxu_{xx}, \qquad u_{xxxx}

往往产生类似扩散的效果。

例如：

ut+aux=νnumuxx.u_t+a u_x = \nu_{\text{num}}u_{xx}.

其中：

νnum\nu_{\text{num}}

并不是物理黏性，而是数值黏性。

它会：

  • 抑制高频误差；
  • 平滑激波；
  • 把尖锐波前抹宽。

因此称为：

数值耗散\boxed{ \text{数值耗散} }

对 Fourier 模态：

u=eikxu=e^{ikx}

有：

uxx=−k2u.u_{xx}=-k^2u.

这是一个实负系数。

所以高波数模态会快速衰减。

而且：

k 越大⇒k2 越大,k\text{ 越大} \Rightarrow k^2\text{ 越大},

因此耗散首先杀死的是短波和网格尺度附近的振荡。

数值色散

奇数阶导数，例如：

uxxxu_{xxx}

具有：

uxxx=−ik3u.u_{xxx} = -ik^3u.

它主要改变 Fourier 模态的相位，而不是直接降低其幅值。

所以会造成：

不同波长以错误的速度传播\boxed{ \text{不同波长以错误的速度传播} }

这叫数值色散。

其表现通常是：

  • 波前附近振铃；
  • 波形拉伸；
  • 波峰和波谷位置偏移。

耗散和色散的统一理解

如果放大因子写成：

G=∣G∣eiϕ,G=|G|e^{i\phi},

那么：

∣G∣|G|

描述振幅变化，对应数值耗散；

ϕ\phi

描述相位变化，对应数值色散。

所以 修正方程分析 和 von Neumann / Fourier 分析 实际上是在从两个角度研究同一组数值误差。

人工黏性

数值耗散有时是格式无意产生的。

但在双曲型激波计算中，人们也会故意加入耗散(artificial viscosity)。

因为激波具有近似间断结构，没有足够耗散时，高阶中心型格式容易出现强烈数值振荡。

适量人工黏性可以：

  • 稳定激波；
  • 消除高频振荡；
  • 保持计算可用。

但人工黏性不是越大越好

过强的人工黏性会：

  • 把激波抹得很宽；
  • 消除真实的小尺度流动；
  • 降低有效分辨率。

因此理想的数值耗散应该，只在真正需要稳定的区域足够强，而在光滑区域尽可能弱。

现代：

  • TVD；
  • flux limiter；
  • Roe；
  • HLLC；
  • ENO/WENO

很多思想都可以理解为在更精细地控制这种耗散。

本章几个方法的总结

双曲型传播问题

重点关注：

  • 信息方向；
  • CFL；
  • 数值色散；
  • 激波稳定性。

典型方法：

  • Lax-Wendroff；
  • MacCormack；
  • 迎风方法。

抛物型问题

重点关注：

Δt∼Δx2\Delta t \sim \Delta x^2

造成的显式稳定性限制。

因此常考虑：

  • Crank-Nicolson；
  • ADI；
  • 其他隐式方法。

椭圆型问题

没有简单时间推进方向。

主要需求是求解大型代数系统。

典型方法：

  • relaxation；
  • Gauss-Seidel；
  • SOR。

现代实际程序中还会大量使用：

  • Krylov 方法；
  • multigrid。

不可压缩流

核心问题是压力—速度耦合

因此要使用：

  • pressure correction；
  • projection；
  • MAC；
  • SIMPLE 类方法。

HPC 视角

不同 CFD 方法的计算结构差别很大。

显式 stencil

例如 Lax-Wendroff、MacCormack：

  • 局部通信；
  • stencil 规则；
  • 并行度高；
  • GPU 友好。

但：

Δt\Delta t

往往受 CFL 严格限制。

隐式算法

例如 Crank-Nicolson：

  • 时间步可以较大；
  • 但每一步要解线性系统。

因此主要成本从：

stencil evaluation\text{stencil evaluation}

转向：

linear solver.\text{linear solver}.

ADI

把二维/三维系统转化成许多线求解。

在传统 CPU 上效率很高。

但不同方向之间存在较强的数据依赖，在现代大规模并行架构上并不一定像显式 stencil 那样容易扩展。

Poisson 求解

不可压缩流每一步往往都要解：

∇2p=b.\nabla^2p=b.

压力 Poisson 方程具有全局耦合性。

因此在大规模不可压 CFD 中：

压力 Poisson 求解器通常是最主要的性能瓶颈之一\boxed{ \text{压力 Poisson 求解器通常是最主要的性能瓶颈之一} }

这也是 multigrid、domain decomposition、Krylov solver 对 HPC CFD 如此重要的原因。

总结

不同 PDE 类型→不同数值困难→不同算法设计\boxed{ \text{不同 PDE 类型} \rightarrow \text{不同数值困难} \rightarrow \text{不同算法设计} }

例如：

双曲型→传播、CFL、色散、激波\text{双曲型} \rightarrow \text{传播、CFL、色散、激波} 抛物型→扩散、时间步限制、隐式求解\text{抛物型} \rightarrow \text{扩散、时间步限制、隐式求解} 椭圆型→全局耦合、迭代求解\text{椭圆型} \rightarrow \text{全局耦合、迭代求解} 不可压缩 NS→压力约束、Poisson 方程\text{不可压缩 NS} \rightarrow \text{压力约束、Poisson 方程}

与此同时，所有离散方法还必须始终考虑：

精度+稳定性+耗散+色散+计算成本\boxed{ \text{精度} + \text{稳定性} + \text{耗散} + \text{色散} + \text{计算成本} }

真正的 CFD 方法设计，本质上就是在这些要求之间寻找平衡。