流体力学核心贯通课 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 方法基于 时间方向的泰勒展开：

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

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

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

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

\[ t^{n+1}=t^n+\Delta t \]

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

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

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

\(e\) 表示 单位质量内能， 若采用理想气体状态方程，还需要补充：

\[ p=(\gamma-1)\rho e \]

3.1 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) \]

在均匀网格上，若空间导数采用 二阶中心差分，有：

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

二阶时间导数可由连续方程再次对时间求导得到：

\[ \rho_{tt} =-\left[ \rho_t u_x +\rho u_{xt} +u_t\rho_x +u\rho_{xt} +\rho_t v_y +\rho v_{yt} +v_t\rho_y +v\rho_{yt} \right] \]

混合导数可以通过对动量方程继续求空间导数得到。由 \(x\) 方向动量方程：

\[ u_t =-\left( u u_x +v u_y +\frac{1}{\rho}p_x \right) \]

3.1 Lax–Wendroff 方法

\[ u_t =-\left( u u_x +v u_y +\frac{1}{\rho}p_x \right) \]

对 \(x\) 求导：

\[ u_{xt} =-\frac{\partial}{\partial x} \left( u u_x +v u_y +\frac{1}{\rho}p_x \right) \]

逐项展开：

\[ \frac{\partial}{\partial x}(u u_x) =u_x^2+u u_{xx} \]

\[ \frac{\partial}{\partial x}(v u_y) =v_xu_y+v u_{xy} \]

\[ \frac{\partial}{\partial x} \left( \frac{1}{\rho}p_x \right) =\frac{1}{\rho}p_{xx} -\frac{\rho_x}{\rho^2}p_x \]

3.1 Lax–Wendroff 方法

类似地，由 \(y\) 方向动量方程：

\[ v_t =-\left( u v_x +v v_y +\frac{1}{\rho}p_y \right) \]

对 \(y\) 求导：

\[ v_{yt} =-\frac{\partial}{\partial y} \left( u v_x +v v_y +\frac{1}{\rho}p_y \right) \]

展开得到：

\[ v_{yt} =-u_yv_x -u v_{xy} -v_y^2 -v v_{yy} -\frac{1}{\rho}p_{yy} +\frac{\rho_y}{\rho^2}p_y \]

3.1 Lax–Wendroff 方法

\[ \rho_{tt} =-\left[ \rho_t u_x +\rho u_{xt} +u_t\rho_x +u\rho_{xt} +\rho_t v_y +\rho v_{yt} +v_t\rho_y +v\rho_{yt} \right] \]

为了计算 \(\rho_{tt}\)， 还需要 \(\rho_{xt}\) 与 \(\rho_{yt}\)。

由连续方程：

\[ \rho_t =-\left( \rho u_x +u\rho_x +\rho v_y +v\rho_y \right) \]

对 \(x\) 求导：

\[ \rho_{xt} =-\frac{\partial}{\partial x} \left( \rho u_x +u\rho_x +\rho v_y +v\rho_y \right) \]

展开：

\[ \rho_{xt} =-2\rho_xu_x -\rho u_{xx} -u\rho_{xx} -\rho_xv_y -\rho v_{xy} -v_x\rho_y -v\rho_{xy} \]

同理，对 \(y\) 求导：

\[ \rho_{yt} =-\rho_yu_x -\rho u_{xy} -u_y\rho_x -u\rho_{xy} -2\rho_yv_y -\rho v_{yy} -v\rho_{yy} \]

3.1 Lax–Wendroff 方法

因此，密度的二阶时间导数可计算：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \rho_{tt} =-\rho_t u_x -\rho u_{xt} -u_t\rho_x -u\rho_{xt} -\rho_t v_y -\rho v_{yt} -v_t\rho_y -v\rho_{yt} } $} \]

其中：

\[ \rho_t =-\left( \rho u_x +u\rho_x +\rho v_y +v\rho_y \right) \]

\[ u_t =-\left( u u_x +v u_y +\frac{1}{\rho}p_x \right) \]

\[ v_t =-\left( u v_x +v v_y +\frac{1}{\rho}p_y \right) \]

\[ u_{xt} =-u_x^2 -u u_{xx} -v_xu_y -vu_{xy} -\frac{1}{\rho}p_{xx} +\frac{\rho_x}{\rho^2}p_x \]

\[ v_{yt} =-u_yv_x -u v_{xy} -v_y^2 -v v_{yy} -\frac{1}{\rho}p_{yy} +\frac{\rho_y}{\rho^2}p_y \]

3.1 Lax–Wendroff 方法

空间离散均采用 中心差分，例如

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

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

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

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

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

方法总结

1. 将 \(n+1\) 时刻的基本物理量在 \(n\) 时刻泰勒展开

2. 用控制方程将时间导数转为空间导数

3. 用空间差分替代空间导数

4. 代入泰勒展开完成时间推进

👉 从而实现： 已知时刻 \(\color{black}{t}\) → 显式推进到 \(\color{black}{t+\Delta t}\)

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

课堂互动：Lax–Wendroff 方法的本质是什么？

从数值离散的角度看，Lax–Wendroff 方法最核心的思想是什么？

3.2 MacCormack 方法

与Lax-Wendroff方法类似，以密度为例：

\[ \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]\) 内的 平均时间导数。

• MacCormack 方法是 Lax–Wendroff 方法的一种变体， 但在实际应用中 要简单得多。

• 与 Lax–Wendroff 方法类似，MacCormack 方法也是一种 显式有限差分方法， 在 时间 和 空间 上均具有 二阶精度。

• 该方法最早于 1969 年 提出，在随后约 15 年 中成为求解流体流动问题最常用的显式差分方法之一。

• 如今，MacCormack 方法在很大程度上已被更复杂的方法所取代，但它仍然是理解 预测-校正思想 和 显式时间推进策略 的重要入门格式。

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）

\[ \rho_t =-\left( \rho u_x +u\rho_x +\rho v_y +v\rho_y \right) \]

在连续方程中，将空间导数用 前向差分（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}+\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(\overline{\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}+\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(\overline{\frac{\partial \rho}{\partial t}}\right)^{t+\Delta t}_{i,j}\right] \]

最终更新

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

3.2 MacCormack 方法（分预测-校正（Predictor–Corrector）两步）

方法总结

MacCormack 方法的关键特点：

• 使用 两步法（预测 + 校正）

• 预测： 前向差分

• 校正： 后向差分

• 自动达到 二阶精度 （避免复杂 二阶时间导数）

方法优势

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

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

• 预测步和校正步可以交换（forward/backward 对调）

• 仍保持 二阶精度

3.2 MacCormack 方法

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

在前面的讨论中，我们通过假设 无粘流、 采用欧拉方程的 非守恒形式， 并考虑 时间推进过程， 来介绍 Lax–Wendroff和 MacCormack方法。

事实上，这些假设并非必要；上述方法同样适用于：

• 黏性流

• 控制方程的守恒形式

• 空间推进（space marching）

下面分别加以说明。

3.2 MacCormack 方法

黏性流

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 方法进行时间推进计算。

3.2 MacCormack 方法

空间推进（Space Marching）

3.2 MacCormack 方法

空间推进（Space Marching）

基本思想

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

\[ x=x_0 \]

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

👉 类比：

• 时间推进： \(\color{black}{t\rightarrow t+\Delta t}\)

• 空间推进： \(\color{black}{x\rightarrow x+\Delta x}\)

3.2 MacCormack 方法

MacCormack 空间推进格式

控制方程形式

考虑如下空间推进形式：

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

其中，\(F\)、\(G\) 和 \(J\) 均为流场变量的函数。

3.2 MacCormack 方法

MacCormack 空间推进格式

预测步（前向差分）

在第 \(i\) 条网格线上，采用 前向差分 近似 \(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 \]

这里， \(\color{black}{\bar{F}^{i+1}_j}\) 表示在第 \(i+1\) 条网格线上的预测值。

3.2 MacCormack 方法

MacCormack 空间推进格式

由预测值计算通量和源项

由于 \(G\) 和 \(J\) 通常依赖于流场变量，因此校正步中应使用由预测值 \(\color{black}{\bar{F}^{i+1}_j}\) 计算得到的预测通量和预测源项：

\[ \bar{G}^{i+1}_j=G\left(\bar{F}^{i+1}_j\right) \]

\[ \bar{J}^{i+1}_j=J\left(\bar{F}^{i+1}_j\right) \]

3.2 MacCormack 方法

MacCormack 空间推进格式

校正步（后向差分）

在第 \(i+1\) 条网格线上，采用 后向差分 近似 \(y\) 方向导数：

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

注意：这里必须使用预测量 \(\color{black}{\bar{G}^{i+1}}\) 和 \(\color{black}{\bar{J}^{i+1}}\)， 而不是最终未知量 \(\color{black}{G^{i+1}}\) 和 \(\color{black}{J^{i+1}}\)。

3.2 MacCormack 方法

MacCormack 空间推进格式

平均导数

将预测步和校正步得到的空间导数取平均：

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

3.2 MacCormack 方法

MacCormack 空间推进格式

最终更新

最终得到第 \(i+1\) 条网格线上的校正值：

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

也可以写成：

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

代入预测步和校正步的具体表达式，可得：

\[ F^{i+1}_j=F^i_j+\frac{\Delta x}{2}\left[J^i_j-\frac{G^i_{j+1}-G^i_j}{\Delta y}+\bar{J}^{i+1}_j-\frac{\bar{G}^{i+1}_j-\bar{G}^{i+1}_{j-1}}{\Delta y}\right] \]

3.2 MacCormack 方法

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

重要区别

空间推进 vs 时间推进：

1. 需要使用 守恒形式

2. 需要进行 变量解码（flux → primitive）

3. 仅适用于 双曲型问题（如超声速流）

总结

• MacCormack 方法不仅适用于时间推进

• 也可以用于空间推进

• 关键在于方程类型（双曲性）

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

另一种看待数值误差的方式 ❶

我们已经讨论了多种求解控制方程的数值方法。通常会这样理解：

但也可以换一个角度：

这个被修改后的方程称为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \text{Modified Equation} } $} \]

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

模型方程与差分格式 ❷

考虑一维线性波动方程：

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

采用如下差分格式：

得到差分方程：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{u_i^{t+\Delta t}-u_i^t}{\Delta t}+a\frac{u_i^t-u_{i-1}^t}{\Delta x}=0 } $} \]

一阶迎风格式（FTBS）。

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

Taylor展开分析 ❸

对时间项作 Taylor 展开：\(u_i^{t+\Delta t}=u_i^t+\left(\frac{\partial u}{\partial t}\right)_i^t\Delta t+\left(\frac{\partial^2u}{\partial t^2}\right)_i^t\frac{(\Delta t)^2}{2}+\left(\frac{\partial^3u}{\partial t^3}\right)_i^t\frac{(\Delta t)^3}{6}+\cdots\)

对空间项作 Taylor 展开：\(u_{i-1}^t=u_i^t-\left(\frac{\partial u}{\partial x}\right)_i^t\Delta x+\left(\frac{\partial^2u}{\partial x^2}\right)_i^t\frac{(\Delta x)^2}{2}-\left(\frac{\partial^3u}{\partial x^3}\right)_i^t\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 \]

得到：

\[ u_t+\frac{\Delta t}{2}u_{tt}+\frac{(\Delta t)^2}{6}u_{ttt}+a\left(u_x-\frac{\Delta x}{2}u_{xx}+\frac{(\Delta x)^2}{6}u_{xxx}\right)+\cdots=0 \]

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

初步整理 ❹

将上一页结果整理为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+au_x=-\frac{\Delta t}{2}u_{tt} -\frac{(\Delta t)^2}{6}u_{ttt} +\frac{a\Delta x}{2}u_{xx} -\frac{a(\Delta x)^2}{6}u_{xxx} +\cdots } $} \]

首先对初步整理式关于 \(t\) 求导：

\[ u_{tt}+au_{xt}=-\frac{\Delta t}{2}u_{ttt}-\frac{(\Delta t)^2}{6}u_{tttt}+\frac{a\Delta x}{2}u_{xxt}-\frac{a(\Delta x)^2}{6}u_{xxxt}+\cdots \]

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

处理二阶时间导数 ❺

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+au_x=-\frac{\Delta t}{2}u_{tt} -\frac{(\Delta t)^2}{6}u_{ttt} +\frac{a\Delta x}{2}u_{xx} -\frac{a(\Delta x)^2}{6}u_{xxx} +\cdots } $} \]

再对初步整理式关于 \(x\) 求导，并乘以 \(a\)：

\[ a u_{tx}+a^2u_{xx}=-\frac{a\Delta t}{2}u_{ttx}-\frac{a(\Delta t)^2}{6}u_{tttx}+\frac{a^2\Delta x}{2}u_{xxx}-\frac{a^2(\Delta x)^2}{6}u_{xxxx}+\cdots \]

\[ u_{tt}+au_{xt}=-\frac{\Delta t}{2}u_{ttt}-\frac{(\Delta t)^2}{6}u_{tttt}+\frac{a\Delta x}{2}u_{xxt}-\frac{a(\Delta x)^2}{6}u_{xxxt}+\cdots \]

两式相减，利用 \(u_{xt}=u_{tx}\)，得到：

\[ \colorbox{#9ACD32}{$\displaystyle \color{black}{ u_{tt}=a^2u_{xx} +\frac{\Delta t}{2}\left(-u_{ttt}+a u_{ttx}+O(\Delta t)\right) +\frac{\Delta x}{2}\left(a u_{xxt}-a^2u_{xxx}+O(\Delta x)\right) } $} \]

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

处理三阶时间导数 ❻

由上一页的结果，保留最低阶关系可得：

\[ u_{tt}=a^2u_{xx}+O(\Delta t,\Delta x) \]

对其关于 \(t\) 求导：

\[ u_{ttt}=a^2u_{xxt}+O(\Delta t,\Delta x) \]

另一方面，由：

\[ u_t+au_x=O(\Delta t,\Delta x) \Longrightarrow u_{xxt}=-au_{xxx}+O(\Delta t,\Delta x) \]

因此：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_{ttt}=-a^3u_{xxx}+O(\Delta t,\Delta x) } $} \]

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

处理混合导数 ❼

由

\[ u_{tt}=a^2u_{xx}+O(\Delta t,\Delta x)\quad\Longrightarrow\quad u_{ttx}=a^2u_{xxx}+O(\Delta t,\Delta x) \]

\[ \left. \begin{aligned} u_{ttt}&=a^2u_{xxt}+O(\Delta t,\Delta x)\\ u_{ttt}&=-a^3u_{xxx}+O(\Delta t,\Delta x) \end{aligned} \right\} \quad\Longrightarrow\quad u_{xxt}=-au_{xxx}+O(\Delta t,\Delta x) \]

于是有：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_{ttx}=a^2u_{xxx}+O(\Delta t,\Delta x),\qquad u_{xxt}=-au_{xxx}+O(\Delta t,\Delta x) } $} \]

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

得到 \(u_{tt}\) 的高阶近似 ❽

将

\[ u_{ttt}=-a^3u_{xxx}+O(\Delta t,\Delta x),\qquad u_{ttx}=a^2u_{xxx}+O(\Delta t,\Delta x),\qquad u_{xxt}=-au_{xxx}+O(\Delta t,\Delta x) \]

代入：

\[ \colorbox{#9ACD32}{$\displaystyle \color{black}{ u_{tt}=a^2u_{xx} +\frac{\Delta t}{2}\left(-u_{ttt}+a u_{ttx}+O(\Delta t)\right) +\frac{\Delta x}{2}\left(a u_{xxt}-a^2u_{xxx}+O(\Delta x)\right) } $} \]

得到：

\[ u_{tt}=a^2u_{xx}+\frac{\Delta t}{2}\left(a^3u_{xxx}+a^3u_{xxx}+O(\Delta t,\Delta x)\right)+\frac{\Delta x}{2}\left(-a^2u_{xxx}-a^2u_{xxx}+O(\Delta t,\Delta x)\right)+\cdots \]

即：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_{tt}=a^2u_{xx}+a^3\Delta t\,u_{xxx}-a^2\Delta x\,u_{xxx}+\cdots } $} \]

这就是教材推导中修正三阶项的关键结果。

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

代回初步整理式 ❾

初步整理式为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+au_x=-\frac{\Delta t}{2}u_{tt} -\frac{(\Delta t)^2}{6}u_{ttt} +\frac{a\Delta x}{2}u_{xx} -\frac{a(\Delta x)^2}{6}u_{xxx} +\cdots } $} \]

代入：

\[ u_{tt}=a^2u_{xx}+a^3\Delta t\,u_{xxx}-a^2\Delta x\,u_{xxx}+\cdots \]

以及：

\[ u_{ttt}=-a^3u_{xxx}+\cdots \]

得到：

\[ u_t+au_x=-\frac{a^2\Delta t}{2}u_{xx}-\frac{a^3(\Delta t)^2}{2}u_{xxx}+\frac{a^2\Delta x\Delta t}{2}u_{xxx}+\frac{a^3(\Delta t)^2}{6}u_{xxx}+\frac{a\Delta x}{2}u_{xx}-\frac{a(\Delta x)^2}{6}u_{xxx}+\cdots \]

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

引入 CFL 数 ❿

将同类项合并：

\[ u_t+au_x=\left(\frac{a\Delta x}{2}-\frac{a^2\Delta t}{2}\right)u_{xx}+\left(\frac{a^2\Delta x\Delta t}{2}-\frac{a^3(\Delta t)^2}{2}+\frac{a^3(\Delta t)^2}{6}-\frac{a(\Delta x)^2}{6}\right)u_{xxx}+\cdots \]

定义 CFL 数：\(\nu=\frac{a\Delta t}{\Delta x}\)

于是：

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

并且：

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

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

一阶迎风格式的修正方程 ⓫

最终得到修正方程：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+au_x=\frac{a\Delta x}{2}(1-\nu)u_{xx}+\frac{a(\Delta x)^2}{6}(-2\nu^2+3\nu-1)u_{xxx}+\cdots } $} \]

其中：

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

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

数值耗散与人工粘性 ⓬

修正方程中的第一项为：

\[ \frac{a\Delta x}{2}(1-\nu)u_{xx} \]

它与粘性扩散项形式相似。因此，可以把

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \mu_{\mathrm{num}}=\frac{a\Delta x}{2}(1-\nu) } $} \]

理解为一阶迎风格式引入的数值粘性。

这就是数值耗散，也常被称为人工粘性。

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

数值色散 ⓭

修正方程中的第二项为：

\[ \frac{a(\Delta x)^2}{6}(3\nu-2\nu^2-1)u_{xxx} \]

这一项为色散项。

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

物理意义：

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

数值色散（Numerical Dispersion）

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

Fourier空间描述

考虑一维线性对流方程：

\[ u_t+au_x=0 \]

其物理含义很简单：

将扰动写成 Fourier 模态：

\[ u(x,t)=\hat{u}e^{i(kx-\omega t)} \]

其中 \(k\) 为波数，\(\omega\) 为角频率。

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

原始方程：所有波数同速传播

对 Fourier 模态，有：

\[ u_t=-i\omega u,\qquad u_x=iku \]

代入原始方程：

\[ u_t+au_x=0 \]

得到：

\[ -i\omega u+aiku=0 \]

因此：

\[ \omega=ak \]

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

原始方程：所有波数同速传播

相速度为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{c_p=\frac{\omega}{k}=a}$} \]

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

二阶导数项：改变振幅

若修正方程中含有二阶导数项：

\[ u_t+au_x=\mu u_{xx} \]

对 Fourier 模态，有：

\[ u_{xx}=(ik)^2u=-k^2u \]

代入可得：

\[ -i\omega u+aiku=-\mu k^2u \Longrightarrow \omega=ak-i\mu k^2 \]

\(\omega\)不再是实数，因此解\(u(x,t)=\hat{u}e^{i(kx-\omega t)}\)中包含衰减因子：

\[ \colorbox{yellow}{$\displaystyle \color{black}{e^{-\mu k^2t}}$} \]

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

二阶导数项对应数值耗散

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

三阶导数项：改变相速度

若修正方程中含有三阶导数项：

\[ u_t+au_x=\beta u_{xxx} \]

对 Fourier 模态，有：

\[ u_{xxx}=(ik)^3u=-ik^3u \]

代入可得：

\[ -i\omega u+aiku=-i\beta k^3u \Longrightarrow \omega=ak+\beta k^3 \]

因此相速度为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{c_p=\frac{\omega}{k}=a+\beta k^2}$} \]

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

三阶导数项对应数值色散

原始方程中：

\[ c_p=a \]

所有波数传播速度相同。

而加入三阶导数项后：

\[ c_p=a+\beta k^2 \]

相速度依赖于波数 \(k\)。

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

偶数阶与奇数阶导数的区别

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

放回修正方程中理解

例如迎风格式的修正方程：

\[ u_t+au_x =\frac{a\Delta x}{2}(1-\nu)u_{xx} +\frac{a(\Delta x)^2}{6}(3\nu-2\nu^2-1)u_{xxx} +\cdots \]

其中：

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

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

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

教学算例：一维线性对流方程 ❶

考虑最简单的一维对流方程：

\[ u_t+cu_x=0,\qquad c>0 \]

其精确解表示：

\[ u(x,t)=u_0(x-ct) \]

本算例用不同差分格式求解同一个方程，比较它们的数值误差特征。

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

为什么选择这个方程？❷

一维线性对流方程虽然简单，但非常适合展示数值格式的本质差异。

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

本代码比较的格式 ❸

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

初始条件设计 ❹

代码中提供三类初始条件：

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

诊断量 ❺

代码中不仅画波形，还计算三个诊断量。

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

FTCS 格式 ❻

FTCS 格式为：

\[ u_i^{n+1}=u_i^n-\frac{\nu}{2}\left(u_{i+1}^n-u_{i-1}^n\right), \qquad \nu=\frac{c\Delta t}{\Delta x} \]

修正方程：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+cu_x =-\frac{c\Delta x}{2}\nu u_{xx} +\frac{c(\Delta x)^2}{6}(\nu^2-1)u_{xxx} +\cdots } $} \]

课堂观察：

\[ \text{波形逐渐失真，误差快速增长} \]

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

Upwind 格式 ❼

当 \(c>0\) 时，采用后向差分：

\[ u_i^{n+1}=u_i^n-\nu\left(u_i^n-u_{i-1}^n\right) \]

其修正方程主导项为：

\[ u_t+cu_x =\frac{c\Delta x}{2}(1-\nu)u_{xx} +\frac{c(\Delta x)^2}{6}(3\nu-2\nu^2-1)u_{xxx} +\cdots \]

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

Upwind 的数值特征 ❽

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

Lax 格式 ❾

Lax 格式为：

\[ u_i^{n+1} =\frac{1}{2}\left(u_{i+1}^n+u_{i-1}^n\right) -\frac{\nu}{2}\left(u_{i+1}^n-u_{i-1}^n\right) \]

它与 FTCS 的主要区别是：

\[ u_i^n \quad\Longrightarrow\quad \frac{1}{2}\left(u_{i+1}^n+u_{i-1}^n\right) \]

修正方程：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_t+cu_x= \frac{c\Delta x}{2\nu}\left(1-\nu^2\right)u_{xx} -\frac{c(\Delta x)^2}{6}\left(1-\nu^2\right)u_{xxx} +\cdots } $} \]

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

Lax 与 Upwind 的对比 ❿

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

Lax-Wendroff 格式 ⓫

Lax-Wendroff 格式为：

\[ u_i^{n+1}=u_i^n -\frac{\nu}{2}\left(u_{i+1}^n-u_{i-1}^n\right) +\frac{\nu^2}{2}\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right) \]

其修正方程主导项为：

\[ u_t+cu_x =\frac{c(\Delta x)^2}{6}(\nu^2-1)u_{xxx} +\cdots \]

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

Lax-Wendroff 的数值特征 ⓬

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

人工粘性格式 ⓭

在 FTCS 基础上加入人工粘性项：

\[ u_i^{n+1}=u_i^n -\frac{\nu}{2}\left(u_{i+1}^n-u_{i-1}^n\right) +\mu_{av}\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right) \]

其中：

\[ \mu_{av} \]

为无量纲人工粘性强度。

对应原始方程可理解为：

\[ u_t+cu_x=\epsilon_{num}u_{xx}+\cdots \]

其中：

\[ \epsilon_{num}=\mu_{av}\frac{(\Delta x)^2}{\Delta t} \]

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

人工粘性的作用 ⓮

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

几种格式的对比 ⓯

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

课堂互动：哪种格式耗散最大？

从本算例的波形和振幅曲线看，哪种格式通常表现出最强的数值耗散？

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

课堂互动：哪种格式最容易振荡？

对于方波或陡峭梯度，哪种格式更容易在间断附近产生前后振荡？

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

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

人工粘性 ❶

对于非定常二维流动，守恒型控制方程可写为：

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

其中：

\[ U=\left[\rho,\rho u,\rho v,\rho\left(e+\frac{V^2}{2}\right)\right] \]

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

人工粘性的基本形式 ❷

在每一个时间步，可向控制方程右端加入人工粘性项：

\[ U_{i,j}^{n+1}=U_{i,j}^{n}+\left(\frac{\partial U}{\partial t}\right)_{i,j}^{n}\Delta t+S_{i,j}^{n} \]

其中 \(S_{i,j}^{n}\) 是人工粘性项。

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

压力传感器：哪里需要加粘性？ ❸

人工粘性通常不希望在全场均匀加入，而应集中作用于激波、间断和强梯度区域。

一种常用做法是利用压力的二阶差分构造传感器：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \phi_{i,j}^{x}= \frac{\left|p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}\right|} {p_{i+1,j}^{n}+2p_{i,j}^{n}+p_{i-1,j}^{n}} } $} \]

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \phi_{i,j}^{y}= \frac{\left|p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}\right|} {p_{i,j+1}^{n}+2p_{i,j}^{n}+p_{i,j-1}^{n}} } $} \]

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

二维人工粘性项 ❹

人工粘性项可写为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \begin{aligned} S_{i,j}^{n}=&C_x\phi_{i,j}^{x}\left(U_{i+1,j}^{n}-2U_{i,j}^{n}+U_{i-1,j}^{n}\right)\\ &+C_y\phi_{i,j}^{y}\left(U_{i,j+1}^{n}-2U_{i,j}^{n}+U_{i,j-1}^{n}\right) \end{aligned} } $} \]

其中：

\[ \phi_{i,j}^{x}= \frac{\left|p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}\right|} {p_{i+1,j}^{n}+2p_{i,j}^{n}+p_{i-1,j}^{n}}, \qquad \phi_{i,j}^{y}= \frac{\left|p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}\right|} {p_{i,j+1}^{n}+2p_{i,j}^{n}+p_{i,j-1}^{n}} \]

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

为什么说它是“四阶耗散”？ ❺

人工粘性项中包含两个二阶差分结构：

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

MacCormack格式中的加入方式 ❻

在预测步中，人工粘性根据已知时刻 \(t\) 的变量计算：

\[ \bar{U}_{i,j}^{n+1} =U_{i,j}^{n} +\left(\frac{\partial U}{\partial t}\right)_{i,j}^{n}\Delta t +S_{i,j}^{n} \]

在校正步中，人工粘性根据预测后的变量重新计算：

\[ U_{i,j}^{n+1} =U_{i,j}^{n} +\left(\frac{\partial U}{\partial t}\right)_{i,j}^{av}\Delta t +\bar{S}_{i,j}^{n+1} \]

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

系数 \(C_x\)、\(C_y\) 如何理解？ ❼

人工粘性中有两个可调参数：

\[ C_x,\qquad C_y \]

它们控制 \(x\) 方向和 \(y\) 方向人工耗散的强弱。

典型范围为（实际计算中通常需要通过数值试验来选择）：

\[ 0.01\le C_x,C_y\le0.3 \]

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

人工粘性的本质 ❽

从修正方程角度看，人工粘性相当于在原方程右端加入额外耗散项：

\[ u_t+cu_x=\epsilon_{\mathrm{num}}u_{xx}+\cdots \]

其中 \(\epsilon_{\mathrm{num}}\) 表示等效数值粘性。

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

传统人工粘性的局限 ❾

人工粘性的一个重要问题是：
它往往具有较强的经验性和人为调参色彩。

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

小结：人工粘性

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

松弛方法是一种有限差分方法，特别适用于求解椭圆型偏微分方程。低速亚声速无粘流通常由椭圆型方程控制（见第3.3.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.4 松弛方法及其在低速无粘流中的应用

差分离散

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

数值方法

理论上可以用：

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

实际 CFD 中：

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

因此采用：

• 迭代方法（松弛法）

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

松弛方法基本思想（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.4 松弛方法及其在低速无粘流中的应用

Gauss-Seidel 方法（更重要）

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

扫描策略（Sweeping）

通常采用：

• 从左到右
• 从下到上

也可以：

• 右 → 左
• 上 → 下

收敛判据

当满足：

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

👉 认为收敛

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

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.4 松弛方法及其在低速无粘流中的应用

评注

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

• Gauss-Seidel 比 Jacobi 更快

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

本节核心理解

👉 本质问题：

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

👉 解决方案：

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

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 方法

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