流体力学核心贯通课 II

2025–2026 学年第二学期

北京理工大学 空天科学与技术学院

任杰

第3章 计算策略

参考书

1. John D. Anderson Computational Fluid Dynamics: The Basics with Applications
   McGraw-Hill, 1995

2. 任玉新 陈海昕
   计算流体力学基础 清华大学出版社，2006

3. 张德良
   计算流体力学教程 高等教育出版社，2010

4. 蔡庆东
   计算流体动力学 北京大学本科生研究生专业基础课讲义，2018

真正理解 CFD 的最好方式，是自己实现数值算法并运行计算 - John D. Anderson

3.1 Lax–Wendroff 方法 ❶

Lax–Wendroff 方法是一种显式有限差分方法。考虑二维非定常无粘流，非守恒形式方程为：

连续方程

\[ \frac{\partial \rho}{\partial t}= -\left(\rho \frac{\partial u}{\partial x}+ u \frac{\partial \rho}{\partial x}+ \rho \frac{\partial v}{\partial y}+ v \frac{\partial \rho}{\partial y}\right) \]

x方向动量方程

\[ \frac{\partial u}{\partial t}= -\left(u \frac{\partial u}{\partial x}+ v \frac{\partial u}{\partial y}+ \frac{1}{\rho}\frac{\partial p}{\partial x}\right) \]

y方向动量方程

\[ \frac{\partial v}{\partial t}= -\left(u \frac{\partial v}{\partial x}+ v \frac{\partial v}{\partial y}+ \frac{1}{\rho}\frac{\partial p}{\partial y}\right) \]

能量方程

\[ \frac{\partial e}{\partial t}= -\left(u \frac{\partial e}{\partial x}+ v \frac{\partial e}{\partial y}+ \frac{p}{\rho}\frac{\partial u}{\partial x}+ \frac{p}{\rho}\frac{\partial v}{\partial y}\right) \]

3.1 Lax–Wendroff 方法 ❷

Lax–Wendroff 方法基于时间方向的泰勒展开。以密度为例：

\[ \rho_{i,j}^{t+\Delta t}=\rho_{i,j}^t+ \left(\frac{\partial \rho}{\partial t}\right)_{i,j}^t \Delta t+ \frac{1}{2}\left(\frac{\partial^2 \rho}{\partial t^2}\right)_{i,j}^t(\Delta t)^2+ \cdots \]

类似地：

\[ u_{i,j}^{t+\Delta t}=u_{i,j}^t+ \left(\frac{\partial u}{\partial t}\right)_{i,j}^t \Delta t+ \frac{1}{2}\left(\frac{\partial^2 u}{\partial t^2}\right)_{i,j}^t(\Delta t)^2 \]

\[ v_{i,j}^{t+\Delta t}=v_{i,j}^t+ \left(\frac{\partial v}{\partial t}\right)_{i,j}^t \Delta t+ \frac{1}{2}\left(\frac{\partial^2 v}{\partial t^2}\right)_{i,j}^t(\Delta t)^2 \]

\[ e_{i,j}^{t+\Delta t}=e_{i,j}^t+ \left(\frac{\partial e}{\partial t}\right)_{i,j}^t \Delta t+ \frac{1}{2}\left(\frac{\partial^2 e}{\partial t^2}\right)_{i,j}^t(\Delta t)^2 \]

3.1 Lax–Wendroff 方法 ❸

一阶时间导数

由连续方程可得：

\[ \left(\frac{\partial \rho}{\partial t}\right)_{i,j}=-\left[\rho_{i,j}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+ u_{i,j}\frac{\rho_{i+1,j}-\rho_{i-1,j}}{2\Delta x}+ \rho_{i,j}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}+ v_{i,j}\frac{\rho_{i,j+1}-\rho_{i,j-1}}{2\Delta y}\right] \]

二阶时间导数

对连续方程再次对时间求导：

\[ \frac{\partial^2 \rho}{\partial t^2}=-\rho \frac{\partial^2 u}{\partial x \partial t}- \frac{\partial u}{\partial t}\frac{\partial \rho}{\partial x}- u \frac{\partial^2 \rho}{\partial x \partial t}- \frac{\partial \rho}{\partial t}\frac{\partial u}{\partial x}+ \rho \frac{\partial^2 v}{\partial y \partial t}+ \frac{\partial v}{\partial t}\frac{\partial \rho}{\partial y}+ v \frac{\partial^2 \rho}{\partial y \partial t} \]

这些混合导数通过对动量方程进一步求导得到。例如右端第一项

\[ \frac{\partial^2 u}{\partial x \partial t}=- u \frac{\partial^2 u}{\partial x^2}+ \left(\frac{\partial u}{\partial x}\right)^2+ v \frac{\partial^2 u}{\partial x \partial y}+ \frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+ \frac{1}{\rho} \frac{\partial^2 p}{\partial x^2}- \frac{1}{\rho^2} \frac{\partial p}{\partial x}\frac{\partial \rho}{\partial x} \]

3.1 Lax–Wendroff 方法 ❹

右端第一项

\[ \frac{\partial^2 u}{\partial x \partial t}=- u \frac{\partial^2 u}{\partial x^2}+ \left(\frac{\partial u}{\partial x}\right)^2+ v \frac{\partial^2 u}{\partial x \partial y}+ \frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+ \frac{1}{\rho} \frac{\partial^2 p}{\partial x^2}- \frac{1}{\rho^2} \frac{\partial p}{\partial x}\frac{\partial \rho}{\partial x} \]

所有项都在时间 \(t\) 处采用二阶中心差分格式表示，即：

\[ \begin{aligned} \left(\frac{\partial^2 u}{\partial x \partial t}\right)^t_{i,j}=& - u^t_{i,j} \frac{u^t_{i+1,j} - 2u^t_{i,j} + u^t_{i-1,j}}{(\Delta x)^2}+ \left(\frac{u^t_{i+1,j} - u^t_{i-1,j}}{2\Delta x}\right)^2 \\&+ v^t_{i,j}\frac{u^t_{i+1,j+1} + u^t_{i-1,j-1} - u^t_{i-1,j+1} - u^t_{i+1,j-1}}{4(\Delta x)(\Delta y)}+ \frac{u^t_{i,j+1} - u^t_{i,j-1}}{2\Delta y}\frac{v^t_{i+1,j} - v^t_{i-1,j}}{2\Delta x} \\ &+ \frac{1}{\rho^t_{i,j}}\frac{p^t_{i+1,j} - 2p^t_{i,j} + p^t_{i-1,j}}{(\Delta x)^2}- \frac{1}{(\rho^t_{i,j})^2}\frac{p^t_{i+1,j} - p^t_{i-1,j}}{2\Delta x}\frac{\rho^t_{i+1,j} - \rho^t_{i-1,j}}{2\Delta x} \end{aligned} \]

右端的所有项都来自于 \(t\) 时刻已知的流场，因此可以计算出左端的数值

3.1 Lax–Wendroff 方法 ❺

通过连续方程关于 \(x\) 求导，并将右端所有导数用二阶中心差分替代即可得到：

\[ \frac{\partial^2 \rho}{\partial x \partial t} \]

进一步，连续方程关于 \(y\) 求导，可得到
\[ \frac{\partial^2 \rho}{\partial y \partial t} \]

法向动量方程关于 \(y\) 求导，可得到 \[ \frac{\partial^2 v}{\partial y \partial t} \]

所有这些导数同样用二阶中心差分表示。

3.1 Lax–Wendroff 方法 ❻

\[ \frac{\partial^2 \rho}{\partial t^2}=-\rho \frac{\partial^2 u}{\partial x \partial t}- \frac{\partial u}{\partial t}\frac{\partial \rho}{\partial x}- u \frac{\partial^2 \rho}{\partial x \partial t}- \frac{\partial \rho}{\partial t}\frac{\partial u}{\partial x}+ \rho \frac{\partial^2 v}{\partial y \partial t}+ \frac{\partial v}{\partial t}\frac{\partial \rho}{\partial y}+ v \frac{\partial^2 \rho}{\partial y \partial t} \]

此时，上式右端唯一剩下的是一阶空间导数：

\[ \frac{\partial u}{\partial x}, \quad \frac{\partial v}{\partial y}, \quad \frac{\partial \rho}{\partial x}, \quad \frac{\partial \rho}{\partial y} \]

这些同样用二阶中心差分表示，例如：

\[ \left(\frac{\partial u}{\partial x}\right)^t_{i,j}=\frac{u_{i+1,j} - u_{i-1,j}}{2\Delta x} \]

3.1 Lax–Wendroff 方法 ❼

此外，还需要一阶时间导数：

\[ \frac{\partial \rho}{\partial t}, \quad \frac{\partial u}{\partial t}, \quad \frac{\partial v}{\partial t} \]

其中：

• \(\frac{\partial \rho}{\partial t}\) 已得到
  
• \(\frac{\partial u}{\partial t}\)、\(\frac{\partial v}{\partial t}\) 可分别流向和法向动量方程，并采用二阶中心差分计算得到

通过上述步骤，我们最终得到：

\[ \frac{\partial^2 \rho}{\partial t^2} \]

因此可以通过如下展开得到下一时刻的密度

\[ \rho_{i,j}^{t+\Delta t}=\rho_{i,j}^t+ \left(\frac{\partial \rho}{\partial t}\right)_{i,j}^t \Delta t+ \frac{1}{2}\left(\frac{\partial^2 \rho}{\partial t^2}\right)_{i,j}^t(\Delta t)^2+ \cdots \]

3.1 Lax–Wendroff 方法 ❽

核心思想总结

1. 用控制方程计算时间导数
   
2. 用空间差分替代空间导数
   
3. 逐步构造高阶时间导数
   
4. 代入泰勒展开完成时间推进

👉 从而实现：已知时刻 ( t ) → 显式推进到 ( t + t )

方法特点

• 显式方法（无需解线性方程组）
• 二阶时间精度
• 二阶空间精度
• 计算简单但代数推导较长

3.2 MacCormack 方法

MacCormack 方法是 Lax–Wendroff 方法的一种变体，但在实际应用中要简单得多。与 Lax–Wendroff 方法类似，MacCormack 方法也是一种显式有限差分方法，在时间和空间上均具有二阶精度。

该方法最早于1969年提出，在随后15年中成为求解流体流动问题最常用的显式差分方法。如今，MacCormack 方法在很大程度上已被更复杂的方法所取代。然而，它具有很强的“教学友好性”，是最容易理解和编程实现的方法之一。此外，其计算结果在许多流动问题中已经足够准确。

基本思想

时间推进求解欧拉方程，假设在时刻 \(t\)，所有网格点上的流场变量已知，目标是计算时刻 \(t + \Delta t\) 的流场变量。以密度为例：

\[ \rho_{i,j}^{t+\Delta t}=\rho_{i,j}^t+\left(\frac{\partial \rho}{\partial t}\right)_{av} \Delta t \]

其中：\(\left(\frac{\partial \rho}{\partial t}\right)_{av}\) 表示时间区间 \([t, t+\Delta t]\) 内的平均时间导数。

3.2 MacCormack 方法

与 Lax–Wendroff 方法的对比

• Lax–Wendroff 方法需要计算二阶时间导数
  
• MacCormack 方法通过平均时间导数实现二阶精度
  
• 从而避免复杂推导

👉 本质优势：避免计算二阶时间导数

类似地：

\[ u_{i,j}^{t+\Delta t}=u_{i,j}^t+\left(\frac{\partial u}{\partial t}\right)_{av} \Delta t \]

\[ v_{i,j}^{t+\Delta t}=v_{i,j}^t+\left(\frac{\partial v}{\partial t}\right)_{av} \Delta t \]

\[ e_{i,j}^{t+\Delta t}=e_{i,j}^t+\left(\frac{\partial e}{\partial t}\right)_{av} \Delta t \]

3.2 MacCormack 方法

预测-校正（Predictor–Corrector）方法

预测步（Predictor Step）

在连续方程中，将空间导数用前向差分（forward difference）表示：

\[ \left(\frac{\partial \rho}{\partial t}\right)^t_{i,j}=-\left[\rho^t_{i,j}\frac{u^t_{i+1,j}-u^t_{i,j}}{\Delta x}+u^t_{i,j}\frac{\rho^t_{i+1,j}-\rho^t_{i,j}}{\Delta x}\right. \]

\[ \left.+\rho^t_{i,j}\frac{v^t_{i,j+1}-v^t_{i,j}}{\Delta y}+v^t_{i,j}\frac{\rho^t_{i,j+1}-\rho^t_{i,j}}{\Delta y}\right] \]

然后利用泰勒展开得到预测值：

\[ \bar{\rho}_{i,j}^{t+\Delta t}=\rho^t_{i,j}+\left(\frac{\partial \rho}{\partial t}\right)^t_{i,j}\Delta t \]

3.2 MacCormack 方法

预测-校正（Predictor–Corrector）方法

预测步（Predictor Step）

👉 注意：

• 该值仅为预测值（predicted value）
• 仅具有一阶精度

类似地：

\[ \bar{u}_{i,j}^{t+\Delta t}=u^t_{i,j}+\left(\frac{\partial u}{\partial t}\right)^t_{i,j}\Delta t \]

\[ \bar{v}_{i,j}^{t+\Delta t}=v^t_{i,j}+\left(\frac{\partial v}{\partial t}\right)^t_{i,j}\Delta t \]

\[ \bar{e}_{i,j}^{t+\Delta t}=e^t_{i,j}+\left(\frac{\partial e}{\partial t}\right)^t_{i,j}\Delta t \]

3.2 MacCormack 方法

预测-校正（Predictor–Corrector）方法

校正步（Corrector Step）

将预测值代入连续方程，并将空间导数改为后向差分（backward difference）：

\[ \left(\frac{\partial \rho}{\partial t}\right)^{t+\Delta t}_{i,j}=-\left[\bar{\rho}_{i,j}^{t+\Delta t}\frac{\bar{u}_{i,j}^{t+\Delta t}-\bar{u}_{i-1,j}^{t+\Delta t}}{\Delta x}+\bar{u}_{i,j}^{t+\Delta t}\frac{\bar{\rho}_{i,j}^{t+\Delta t}-\bar{\rho}_{i-1,j}^{t+\Delta t}}{\Delta x}\right. \]

\[ \left.+\bar{\rho}_{i,j}^{t+\Delta t}\frac{\bar{v}_{i,j}^{t+\Delta t}-\bar{v}_{i,j-1}^{t+\Delta t}}{\Delta y}+\bar{v}_{i,j}^{t+\Delta t}\frac{\bar{\rho}_{i,j}^{t+\Delta t}-\bar{\rho}_{i,j-1}^{t+\Delta t}}{\Delta y}\right] \]

平均时间导数

\[ \left(\frac{\partial \rho}{\partial t}\right)_{av}=\frac{1}{2}\left[\left(\frac{\partial \rho}{\partial t}\right)^t_{i,j}+\left(\frac{\partial \rho}{\partial t}\right)^{t+\Delta t}_{i,j}\right] \]

3.2 MacCormack 方法

预测-校正（Predictor–Corrector）方法

最终更新

\[ \rho_{i,j}^{t+\Delta t}=\rho_{i,j}^t+\left(\frac{\partial \rho}{\partial t}\right)_{av} \Delta t \]

方法总结

MacCormack 方法的关键特点：

• 使用两步法（预测 + 校正）
• 预测：前向差分
  
• 校正：后向差分
  
• 自动达到二阶精度
  
• 避免复杂二阶时间导数

方法优势

✔ 实现简单 ✔ 精度高（与 Lax–Wendroff 相同） ✔ 计算量更低 ✔ 易于编程

3.2 MacCormack 方法

👉 前向预测 + 后向修正 + 平均

本质是一种：对时间导数进行对称近似的策略

• 预测步和校正步可以交换（forward/backward 对调）
• 仍保持二阶精度
  
• 可用于构造交替时间推进格式

👉 MacCormack 方法 = “不用二阶导数的 Lax–Wendroff”

3.2 MacCormack 方法

黏性流、守恒形式与空间推进

在前面的讨论中，我们通过假设无粘流、采用欧拉方程的非守恒形式，并考虑时间推进过程，来介绍 Lax–Wendroff（第6.2节）和 MacCormack（第6.3节）方法。事实上，这些假设并非必要；上述方法同样适用于：

• 黏性流
• 控制方程的守恒形式
• 空间推进（space marching）

下面分别加以说明。

3.2 MacCormack 方法

黏性流

黏性流由 Navier–Stokes 方程控制。在稳态情况下，这些方程在数学上表现为部分椭圆型，因此 Lax–Wendroff 和 MacCormack 方法不适用于椭圆型方程。

然而，对于非定常 Navier–Stokes 方程，其性质为抛物型与椭圆型的混合，因此上述方法是适用的。事实上，MacCormack 方法已被广泛用于通过时间推进求解非定常 Navier–Stokes 方程。

其基本思想相同：

• 时间导数放在方程左侧
  
• 空间导数放在方程右侧
  
• 预测步使用前向差分
  
• 校正步使用后向差分

唯一的区别在于：Navier–Stokes 方程中的空间导数项比欧拉方程更多。

注：对流项可以采用上述差分方式，而黏性项通常建议在预测和校正步中均采用中心差分。

3.2 MacCormack 方法

守恒形式

为简单起见，仍以欧拉方程为例。适用于CFD计算的守恒形式可写为：

\[ \frac{\partial U}{\partial t}=- \frac{\partial F}{\partial x}- \frac{\partial G}{\partial y}+ J \]

其中：

• \(U, F, G, J\) 为向量
• \(U\) 的分量包括：\(\rho, \rho u, \rho v, \rho(e + V^2/2)\)

这些量可以通过 Lax–Wendroff 或 MacCormack 方法进行时间推进计算。

⚠️ 注意：由于守恒变量是通量变量，在每个时间步结束后，需要将其转换为原始变量（primitive variables）。

3.2 MacCormack 方法

空间推进（Space Marching）

为了说明空间推进思想，考虑图6.3中的二维流动。假设：

• 流动方向为 \(x\) 方向（从左到右）
• 流动为无粘流（欧拉方程）

守恒形式下，方程可写为：

\[ \frac{\partial F}{\partial x}=J - \frac{\partial G}{\partial y} \]

适用条件

• 亚声速流：方程为椭圆型 ❌ 不适用空间推进
  
• 超声速流：方程为双曲型 ✅ 可采用空间推进

3.2 MacCormack 方法

空间推进（Space Marching）

基本思想

假设在某一条垂直线（初始数据线）上，流场变量已知：

\[ x = x_0 \]

则可以沿 \(x\) 方向逐步推进计算。

👉 类比：

• 时间推进：\(t \rightarrow t+\Delta t\)
  
• 空间推进：\(x \rightarrow x+\Delta x\)

3.2 MacCormack 方法

MacCormack 空间推进格式

更新公式

\[ F_j^{i+1}=F_j^i+\left(\frac{\partial F}{\partial x}\right)_{av} \Delta x \]

预测步

将式 (6.24) 中的 \(y\) 导数用前向差分表示：

\[ \left(\frac{\partial F}{\partial x}\right)^i_j=J^i_j-\frac{G^i_{j+1} - G^i_j}{\Delta y} \]

然后得到预测值：

\[ \bar{F}^{i+1}_j=F^i_j+\left(\frac{\partial F}{\partial x}\right)^i_j \Delta x \]

3.2 MacCormack 方法

MacCormack 空间推进格式

校正步

利用预测值，采用后向差分：

\[ \left(\frac{\partial F}{\partial x}\right)^{i+1}_j=J^{i+1}_j-\frac{G^{i+1}_j - G^{i+1}_{j-1}}{\Delta y} \]

平均导数

\[ \left(\frac{\partial F}{\partial x}\right)_{av}=\frac{1}{2}\left[\left(\frac{\partial F}{\partial x}\right)^i_j+\left(\frac{\partial F}{\partial x}\right)^{i+1}_j\right] \]

最终更新

\[ F_j^{i+1}=F_j^i+\left(\frac{\partial F}{\partial x}\right)_{av} \Delta x \]

3.2 MacCormack 方法

空间推进本质上是：👉 把空间方向当作“时间”来推进

重要区别

空间推进 vs 时间推进：

1. 需要使用守恒形式
2. 需要进行变量解码（flux → primitive）
3. 仅适用于双曲型问题（如超声速流）

总结

• MacCormack 方法不仅适用于时间推进
  
• 也可以用于空间推进
  
• 关键在于方程类型（双曲性）

👉 空间推进 = “沿流向做时间推进”

3.3 松弛方法及其在低速无粘流中的应用

松弛方法是一种有限差分方法，特别适用于求解椭圆型偏微分方程。低速亚声速无粘流通常由椭圆型方程控制（见第3.4.3节），因此松弛方法常用于此类问题的求解。

松弛方法可以是显式或隐式形式。本节介绍一种显式松弛方法，也称为：逐点迭代方法（point-iterative method）

控制方程

考虑一个二维、不可压、无粘、无旋流动，其控制方程为：

• 速度势函数： \[ \mathbf{V} = \nabla \Phi \]

• 控制方程（Laplace方程）： \[ \frac{\partial^2 \Phi}{\partial x^2} + \frac{\partial^2 \Phi}{\partial y^2} = 0 \]

3.3 松弛方法及其在低速无粘流中的应用

差分离散

3.3 松弛方法及其在低速无粘流中的应用

数值方法

理论上可以用：

• 克拉默法则（Cramer’s rule） ❌（计算量极大）
• 高斯消元 ✔（可行但复杂）

实际 CFD 中：

• 网格数量可达上千甚至更多
  
• 👉 直接法不可行

因此采用：

• 迭代方法（松弛法）

3.3 松弛方法及其在低速无粘流中的应用

松弛方法基本思想（Jacobi方法）

在第 \(n\) 次迭代中：

• 4个邻点 → 已知，当前点 → 未知，从而解出 \(\Phi_{i,j}\)：

\[ \Phi_{i,j}^{n+1}=\frac{(\Delta x)^2 (\Delta y)^2}{2(\Delta y)^2 + 2(\Delta x)^2}\left[\frac{\Phi_{i+1,j}^n + \Phi_{i-1,j}^n}{(\Delta x)^2}+\frac{\Phi_{i,j+1}^n + \Phi_{i,j-1}^n}{(\Delta y)^2}\right] \]

• 所有右端项使用 上一轮迭代值
• 同步更新

步骤：

1. 假设初值
2. 用上式计算新值
3. 全部更新后进入下一轮

3.3 松弛方法及其在低速无粘流中的应用

Gauss-Seidel 方法（更重要）

3.3 松弛方法及其在低速无粘流中的应用

扫描策略（Sweeping）

通常采用：

• 从左到右
• 从下到上

也可以：

• 右 → 左
• 上 → 下

收敛判据

当满足：

\[ |\Phi_{i,j}^{n+1} - \Phi_{i,j}^{n}| < \epsilon \]

👉 认为收敛

3.3 松弛方法及其在低速无粘流中的应用

10. 超松弛（SOR）

定义中间值：

\[ \tilde{\Phi}_{i,j}^{n+1}=\frac{(\Delta x)^2 (\Delta y)^2}{2(\Delta y)^2 + 2(\Delta x)^2}\left[\frac{\Phi_{i+1,j}^n + \Phi_{i-1,j}^{n+1}}{(\Delta x)^2}+\frac{\Phi_{i,j+1}^n + \Phi_{i,j-1}^{n+1}}{(\Delta y)^2}\right] \]

然后：

\[ \Phi_{i,j}^{n+1}=\Phi_{i,j}^{n}+\omega \left(\tilde{\Phi}_{i,j}^{n+1} - \Phi_{i,j}^{n}\right) \]

松弛因子 \(\omega\)

• \(\omega = 1\) → Gauss-Seidel
• \(\omega > 1\) → 超松弛（加速）
• \(\omega < 1\) → 欠松弛

3.3 松弛方法及其在低速无粘流中的应用

评注

• 松弛法本质是： 局部更新 + 全局收敛

• Gauss-Seidel 比 Jacobi 更快

• SOR 可大幅加速：👉 收敛步数可减少一个数量级（甚至30倍）

本节核心理解

👉 本质问题：

• 椭圆方程 → 全局耦合
• 不能逐点显式推进

👉 解决方案：

• 用“假时间”（迭代）逼近稳态解

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

现实往往并不像表面看起来那样简单——CFD 亦是如此。我们已经讨论了多种用于求解控制方程的数值方法。我们一直采用这样一种观点：

数值解是带有截断误差和舍入误差的近似解

但实际上，还有另一种理解方式。

考虑一维波动方程：

\[ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 \]

其中 \(a>0\)。时间：一阶前向差分；空间：一阶后向差分。得到差分方程：

\[ \frac{u_i^{t+\Delta t} - u_i^t}{\Delta t}+a \frac{u_i^t - u_{i-1}^t}{\Delta x}= 0 \]

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

Taylor展开分析

对差分项进行展开：

\[ u_i^{t+\Delta t}=u_i^t+\left(\frac{\partial u}{\partial t}\right)_i \Delta t+\left(\frac{\partial^2 u}{\partial t^2}\right)_i \frac{(\Delta t)^2}{2}+\left(\frac{\partial^3 u}{\partial t^3}\right)_i \frac{(\Delta t)^3}{6}+ \cdots \]

\[ u_{i-1}^t=u_i^t-\left(\frac{\partial u}{\partial x}\right)_i \Delta x+\left(\frac{\partial^2 u}{\partial x^2}\right)_i \frac{(\Delta x)^2}{2}-\left(\frac{\partial^3 u}{\partial x^3}\right)_i \frac{(\Delta x)^3}{6}+ \cdots \]

代入差分方程， \[ \frac{u_i^{t+\Delta t} - u_i^t}{\Delta t}+a \frac{u_i^t - u_{i-1}^t}{\Delta x}= 0 \]

得到：

\[ \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=-\frac{\partial^2 u}{\partial t^2}\frac{\Delta t}{2}-\frac{\partial^3 u}{\partial t^3}\frac{(\Delta t)^2}{6}+a\frac{\partial^2 u}{\partial x^2}\frac{\Delta x}{2}-a\frac{\partial^3 u}{\partial x^3}\frac{(\Delta x)^2}{6}+\cdots \]

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

修正方程（Modified Equation）

\[ \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=-\frac{\partial^2 u}{\partial t^2}\frac{\Delta t}{2}-\frac{\partial^3 u}{\partial t^3}\frac{(\Delta t)^2}{6}+a\frac{\partial^2 u}{\partial x^2}\frac{\Delta x}{2}-a\frac{\partial^3 u}{\partial x^3}\frac{(\Delta x)^2}{6}+\cdots \]

通过进一步推导（用时间导数替换为空间导数）最终得到：

\[ \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=\frac{a \Delta x}{2}(1-\nu)\frac{\partial^2 u}{\partial x^2}+\frac{a(\Delta x)^2}{6}(3\nu - 2\nu^2 -1)\frac{\partial^3 u}{\partial x^3}+ \cdots \]

其中：

\[ \nu = \frac{a \Delta t}{\Delta x} \]

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

关键思想：数值解 ≠ 原方程解

重要结论：差分方程的“精确解”，实际上是某个不同偏微分方程（修正方程）的解

即：

• 数值方法并不是在解原始 PDE，而是在解：带有额外项的“修正方程”

数值耗散（Numerical Dissipation）

观察

\[ \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=\frac{a \Delta x}{2}(1-\nu)\frac{\partial^2 u}{\partial x^2}+\frac{a(\Delta x)^2}{6}(3\nu - 2\nu^2 -1)\frac{\partial^3 u}{\partial x^3}+ \cdots \]

其中

\[ \frac{a \Delta x}{2}(1-\nu)\frac{\partial^2 u}{\partial x^2} \]

类似粘性项（\(\mu \partial^2 u/\partial x^2\)）👉 结论：数值方法会引入“类粘性效应”

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

物理意义：

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

数值色散（Numerical Dispersion）

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

人工粘性（Artificial Viscosity）

在CFD中：数值耗散 ≈ 人工粘性

定义：人为加入项

\[ S_{i,j}=C_x \frac{|p_{i+1,j}-2p_{i,j}+p_{i-1,j}|}{p_{i+1,j}+2p_{i,j}+p_{i-1,j}}(U_{i+1,j}-2U_{i,j}+U_{i-1,j})+ \cdots \]

特点：

• 类似四阶耗散
• 控制震荡
• 提高稳定性

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

人工粘性的作用

优点：

• 提高稳定性
• 抑制振荡
• 捕捉激波

缺点：

• 降低精度
• 过度平滑

实际计算现象（Fig. 6.8–6.10）

随着人工粘性增大：

✔ 好处：

• 解更稳定
• 振荡减少

✘ 坏处：

• 激波被“抹平”
• 结构变模糊

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

关键认知

👉 CFD的核心矛盾：稳定性 vs 精度

评注

1. 数值解对应“改进方程”
2. 数值误差具有物理含义：
   • 耗散 → 类粘性
   • 色散 → 相速度误差
3. 人工粘性是工程手段：
   • 必要但危险

工程经验

• 激波问题 → 必须加人工粘性
• 光滑流 → 尽量减少
• 参数需要调节（经验性强）

现代发展

近年来：自适应人工粘性（如 TVD 方法）

特点：

• 只在需要的地方加入
• 自动控制强度

3.5 交替方向隐式（ADI）方法

让我们回到隐式格式的讨论，以第2章介绍的 Crank–Nicolson 方法为例。本节中，我们考虑一个推进求解（marching）的例子；将 \(t\) 作为推进变量。

👉 Crank–Nicolson格式

\[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \]

\[ \frac{T_i^{n+1} - T_i^n}{\Delta t}=\alpha \frac{1}{2}\left[\frac{T_{i+1}^{n} - 2T_i^{n} + T_{i-1}^{n}}{(\Delta x)^2}+\frac{T_{i+1}^{n+1} - 2T_i^{n+1} + T_{i-1}^{n+1}}{(\Delta x)^2}\right] \]

该方程只有一个额外的自变量，即 \(x\)。由于我们处理的是线性方程，因此采用 Crank–Nicolson 方法得到的隐式解可以直接通过 Thomas 算法（见教材附录A）求解。

3.5 交替方向隐式（ADI）方法

👉 多维问题带来的困难

当问题涉及多个空间变量时（三维或二维问题），差分方程将不再保持三对角结构。

考虑二维非定常导热方程：

\[ \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \]

采用 Crank–Nicolson 离散，可得：

\[ \frac{T_{i,j}^{n+1} - T_{i,j}^{n}}{\Delta t}=\alpha \frac{1}{2} \left[\frac{T_{i+1,j}^{n+1} -2T_{i,j}^{n+1} + T_{i-1,j}^{n+1}}{(\Delta x)^2}+\frac{T_{i+1,j}^{n} -2T_{i,j}^{n} + T_{i-1,j}^{n}}{(\Delta x)^2} \right] \]

\[ + \alpha \frac{1}{2} \left[\frac{T_{i,j+1}^{n+1} -2T_{i,j}^{n+1} + T_{i,j-1}^{n+1}}{(\Delta y)^2}+\frac{T_{i,j+1}^{n} -2T_{i,j}^{n} + T_{i,j-1}^{n}}{(\Delta y)^2}\right] \]

3.5 交替方向隐式（ADI）方法

关键问题

\[ \frac{T_{i,j}^{n+1} - T_{i,j}^{n}}{\Delta t}=\alpha \frac{1}{2} \left[\frac{T_{i+1,j}^{n+1} -2T_{i,j}^{n+1} + T_{i-1,j}^{n+1}}{(\Delta x)^2}+\frac{T_{i+1,j}^{n} -2T_{i,j}^{n} + T_{i-1,j}^{n}}{(\Delta x)^2} \right] \]

\[ + \alpha \frac{1}{2} \left[\frac{T_{i,j+1}^{n+1} -2T_{i,j}^{n+1} + T_{i,j-1}^{n+1}}{(\Delta y)^2}+\frac{T_{i,j+1}^{n} -2T_{i,j}^{n} + T_{i,j-1}^{n}}{(\Delta y)^2}\right] \]

该方程包含 5 个未知量：

\[ T_{i,j}^{n+1},\; T_{i+1,j}^{n+1},\; T_{i-1,j}^{n+1},\; T_{i,j+1}^{n+1},\; T_{i,j-1}^{n+1} \]

👉 因此：

• 不再是三对角矩阵、Thomas算法无法直接使用
• 需要求解大型稀疏线性系统（计算代价高）

3.5 交替方向隐式（ADI）方法

ADI 方法思想

为了恢复三对角结构，引入： 交替方向隐式（Alternating-Direction Implicit, ADI）方法

核心思想：

• 将二维问题拆成两个一维问题
• 分两步推进时间

第一步（x方向隐式）

在时间区间 \(\Delta t/2\) 内：x方向 👉 隐式、y方向 👉 显式

\[ \frac{T_{i,j}^{n+1/2} - T_{i,j}^{n}}{\Delta t/2}=\alpha \frac{T_{i+1,j}^{n+1/2} -2T_{i,j}^{n+1/2} + T_{i-1,j}^{n+1/2}}{(\Delta x)^2}+\alpha \frac{T_{i,j+1}^{n} -2T_{i,j}^{n} + T_{i,j-1}^{n}}{(\Delta y)^2} \]

3.5 交替方向隐式（ADI）方法

整理得到三对角形式：

\[ A T_{i-1,j}^{n+1/2} - B T_{i,j}^{n+1/2} + A T_{i+1,j}^{n+1/2} = K_i \]

其中：

\[ A = \frac{\alpha \Delta t}{2(\Delta x)^2} \]

\[ B = 1 + \frac{\alpha \Delta t}{(\Delta x)^2} \]

\[ K_i =- T_{i,j}^{n}- \frac{\alpha \Delta t}{2(\Delta y)^2}\left(T_{i,j+1}^{n} -2T_{i,j}^{n} + T_{i,j-1}^{n} \right) \]

求解方法（x方向扫掠）

固定 \(j\)：对 \(i = 1 \rightarrow N\) 扫掠（每一行使用 Thomas 算法）

总共执行：\(M\) 次扫掠（每个 \(j\) 一次）👉 得到：\(T_{i,j}^{n+1/2}\)

3.5 交替方向隐式（ADI）方法

第二步（y方向隐式）

从 \(t + \Delta t/2\) 推进到 \(t + \Delta t\)

此时：

• y方向 👉 隐式、x方向 👉 显式

离散格式：

\[ \frac{T_{i,j}^{n+1} - T_{i,j}^{n+1/2}}{\Delta t/2}=\alpha \frac{T_{i+1,j}^{n+1/2} -2T_{i,j}^{n+1/2} + T_{i-1,j}^{n+1/2}}{(\Delta x)^2}+\alpha \frac{T_{i,j+1}^{n+1} -2T_{i,j}^{n+1} + T_{i,j-1}^{n+1}}{(\Delta y)^2} \]

化为三对角形式

\[ C T_{i,j+1}^{n+1} - D T_{i,j}^{n+1} + C T_{i,j-1}^{n+1} = L_j \]

3.5 交替方向隐式（ADI）方法

第二步（y方向隐式）

化为三对角形式

\[ C T_{i,j+1}^{n+1} - D T_{i,j}^{n+1} + C T_{i,j-1}^{n+1} = L_j \]

其中：

\[ C = \frac{\alpha \Delta t}{2(\Delta y)^2} \]

\[ D = 1 + \frac{\alpha \Delta t}{(\Delta y)^2} \]

\[ L_j =- T_{i,j}^{n+1/2}- \frac{\alpha \Delta t}{2(\Delta x)^2}\left(T_{i+1,j}^{n+1/2} -2T_{i,j}^{n+1/2} + T_{i-1,j}^{n+1/2}\right) \]

3.5 交替方向隐式（ADI）方法

3.5 交替方向隐式（ADI）方法

计算流程

1. x方向隐式（半步）
2. y方向隐式（完整步）

评注

• 将二维问题拆为两个一维问题
• 时间：二阶精度、空间：二阶精度。截断误差：\(O\left( (\Delta t)^2, (\Delta x)^2, (\Delta y)^2 \right)\)
• 每一步都是三对角系统
• 可高效使用 Thomas 算法
• ADI方法广泛应用于：导热问题、粘性流动、抛物型偏微分方程

3.6 不可压粘性流的压力修正方法

在本章，我们已经讨论了求解无粘不可压流动的一种数值方法，即松弛法。
无粘不可压流动由椭圆型偏微分方程控制，而松弛法本质上是一种迭代方法，适用于求解椭圆问题。

然而：

黏性不可压流动由Navier–Stokes方程控制，其具有混合的椭圆-抛物特性

因此，松弛方法并不适用于此类问题。

本节介绍一种新的方法：压力修正方法（pressure correction technique）

该方法：

• 广泛应用于不可压Navier–Stokes方程
• 由 Patankar 和 Spalding 提出
• 核心算法：SIMPLE方法

3.6 不可压粘性流的压力修正方法

\[ \text{连续方程:}\quad \nabla \cdot \mathbf{V} = 0 \]

\[ \text{x动量:}\quad\rho \frac{Du}{Dt}=-\frac{\partial p}{\partial x}+ \mu \nabla^2 u + \rho f_x \]

\[ \text{y动量:}\quad\rho \frac{Dv}{Dt}=-\frac{\partial p}{\partial y}+ \mu \nabla^2 v + \rho f_y \]

\[ \text{z动量:}\quad\rho \frac{Dw}{Dt}=-\frac{\partial p}{\partial z}+ \mu \nabla^2 w + \rho f_z \]

👉 数值困难

对于可压缩流：声速有限 → 时间步长有限

对于不可压流：声速 \(\to \infty\)。显式方法要求 \(\Delta t \to 0\)

👉 不可行！

3.6 不可压粘性流的压力修正方法

中心差分问题：奇偶失联

3.6 不可压粘性流的压力修正方法

中心差分问题：奇偶失联

3.6 不可压粘性流的压力修正方法

解决方法：交错网格（Staggered Grid）

• 压力：单元中心
• 速度：单元面

连续方程变为：

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

👉 使用相邻点 → 消除奇偶失联（棋盘格）

3.6 不可压粘性流的压力修正方法

压力修正思想

核心步骤：

1. 假设压力 \(p^*\)
2. 解动量方程 → 得 \(u^*, v^*\)
3. 不满足连续方程
4. 引入修正：

\[ p = p^* + p', \qquad u = u^* + u', \qquad v = v^* + v' \]

压力修正方程推导

从动量方程离散：

\[ (\rho u)^{n+1}_{i+1/2,j}=(\rho u)^n + A\Delta t- \frac{\Delta t}{\Delta x}(p'_{i+1,j} - p'_{i,j}) \]

3.6 不可压粘性流的压力修正方法

压力修正方程推导

\[ (\rho u)^{n+1}_{i+1/2,j}=(\rho u)^n + A\Delta t- \frac{\Delta t}{\Delta x}(p'_{i+1,j} - p'_{i,j}) \]

\[ (\rho v)^{n+1}_{i,j+1/2}=(\rho v)^n + B\Delta t- \frac{\Delta t}{\Delta y}(p'_{i,j+1} - p'_{i,j}) \tag{6.93} \]

代入连续方程

\[ \frac{(\rho u)_{i+1/2,j} - (\rho u)_{i-1/2,j}}{\Delta x}+\frac{(\rho v)_{i,j+1/2} - (\rho v)_{i,j-1/2}}{\Delta y}= 0 \]

得到压力修正方程

\[ a p'_{i,j}+ b p'_{i+1,j}+ b p'_{i-1,j}+ c p'_{i,j+1}+ c p'_{i,j-1}+ d = 0 \]

3.6 不可压粘性流的压力修正方法

压力修正方程推导

\[ a p'_{i,j}+ b p'_{i+1,j}+ b p'_{i-1,j}+ c p'_{i,j+1}+ c p'_{i,j-1}+ d = 0 \]

其中：

\[ a = 2\left(\frac{\Delta t}{\Delta x^2} + \frac{\Delta t}{\Delta y^2}\right) \]

\[ b = -\frac{\Delta t}{\Delta x^2},\quad c = -\frac{\Delta t}{\Delta y^2} \]

\[ d =\frac{1}{\Delta x}\left[(\rho u^*)_{i+1/2,j} - (\rho u^*)_{i-1/2,j}\right]+\frac{1}{\Delta y}\left[(\rho v^*)_{i,j+1/2} - (\rho v^*)_{i,j-1/2}\right] \]

3.6 不可压粘性流的压力修正方法

该方程本质是：

\[ \nabla^2 p' = Q \]

👉 即：Poisson方程

3.6 不可压粘性流的压力修正方法

SIMPLE算法流程

1. 猜测 \(p^*\)
2. 解动量方程 → \(u^*, v^*\)
3. 解压力修正方程 → \(p'\)
4. 更新压力：

\[ p = p^* + p' \]

5. 修正速度
6. 重复迭代直到收敛

3.6 不可压粘性流的压力修正方法

SIMPLE算法流程

边界条件

• 入口：\(p'=0\)
• 出口：\(p'=0\)
• 壁面：

\[ \frac{\partial p}{\partial y} = 0 \]

3.6 不可压粘性流的压力修正方法

SIMPLE算法流程

核心思想 👉 压力修正 = 用Poisson方程强制满足连续方程

方法本质

• 将不可压问题转化为椭圆问题
• 通过迭代逼近质量守恒

广泛应用于：

• CFD工业软件
• 不可压流动模拟
• SIMPLE / SIMPLER / PISO 方法

👉 压力修正方法 = 用压力调整速度，使流场满足连续方程