流体力学核心贯通课 II

2025–2026 学年第二学期

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

任杰

第2章 数值方法与网格变换

参考书

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

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

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

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

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

2.1 离散化 ❶

定义：

• 离散化（Discretization）：用有限信息逼近无限连续
• 将连续问题（continuous problem）转化为离散问题（discrete problem）

\[ \frac{\partial u}{\partial x} \rightarrow \frac{u_{i+1}-u_i}{\Delta x} \]

➡ 微分算子 → 差分算子
➡ PDE → 代数方程组

2.1 离散化 ❷

结构化网格与非结构化网格 Structured vs Unstructured

2.1 离散化 ❸

结构化网格

• 坐标索引： \[ (i,j) \]

• 相邻点：

  • 右：\((i+1,j)\)
  • 左：\((i-1,j)\)
  • 上：\((i,j+1)\)
  • 下：\((i,j-1)\)

网格间距

• x方向：\(\Delta x\)，y方向：\(\Delta y\)
• 均匀网格（uniform grid），非均匀网格（non-uniform grid）

2.1 离散化 ❹

设函数 \(u(x)\)，在 \(x_i\) 处Taylor展开：

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

示例

\[ f(x) = \sin(2\pi x) \]

已知 \(f(0.2)\) 求 \(f(0.2+0.02)\)

• 仅用第一项 \(f(0.22) \approx f(0.2)\)，误差 ≈ 3.176%
• 加入一阶导数 \(f(0.22) \approx f(0.2) + f'(0.2)\times0.02\)，误差 ≈ 0.775%
• 加入二阶导数 \(f(0.22) \approx f(0.2) + f'(0.2)\times0.02 + \frac{f''(0.2)}{2}\times(0.02)^2\)，误差 ≈ 0.01%

2.1 离散化 ❺

从泰勒展开到有限差分

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

前向差分（Forward Difference）

\[ \left(\frac{\partial u}{\partial x}\right)_i = \frac{u_{i+1}-u_i}{\Delta x} +O(\Delta x) \]

后向差分（Backward Difference）

\[ \left(\frac{\partial u}{\partial x}\right)_i = \frac{u_i-u_{i-1}}{\Delta x}+O(\Delta x) \]

中心差分（Central Difference）

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

2.1 离散化 ❻

差分格式推导：以中心差分为例

在点 \(x_i\) 处，对 \(u_{i+1}\) 和 \(u_{i-1}\) 分别做泰勒展开：

\[ u_{i+1} = u_i + u'_i \Delta x + \frac{1}{2}u''_i (\Delta x)^2 + \frac{1}{6}u'''_i (\Delta x)^3 + O(\Delta x^4) \]

\[ u_{i-1} = u_i - u'_i \Delta x + \frac{1}{2}u''_i (\Delta x)^2 - \frac{1}{6}u'''_i (\Delta x)^3 + O(\Delta x^4) \]

将两式相减（消去了全部偶次项）：

\[ u_{i+1} - u_{i-1} = 2 u'_i \Delta x + \frac{1}{3}u'''_i (\Delta x)^3 + O(\Delta x^5) \]

整理得：

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

2.1 离散化 ❼

二阶前向差分

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

四阶中心差分

\[ \left(\frac{\partial u}{\partial x}\right)_i= \frac{-u_{i+2} + 8u_{i+1} - 8u_{i-1} + u_{i-2}}{12\Delta x}+ O(\Delta x^4) \]

五阶迎风偏置差分

\[ \left(\frac{\partial u}{\partial x}\right)_i=\frac{-2u_{i-3}+15u_{i-2}-60u_{i-1}+20u_i+30u_{i+1}-3u_{i+2}}{60\Delta x}+ O(\Delta x^5) \]

六阶中心差分？

\[ \left(\frac{\partial u}{\partial x}\right)_i= \frac{-u_{i+3} + 9u_{i+2} - 45u_{i+1} + 45u_{i-1} - 9u_{i-2} + u_{i-3}}{60\Delta x}+ O(\Delta x^6) \]

2.1 离散化 ❽

六阶中心差分？

\[ \left(\frac{\partial u}{\partial x}\right)_i= \frac{-u_{i+3} + 9u_{i+2} - 45u_{i+1} + 45u_{i-1} - 9u_{i-2} + u_{i-3}}{60\Delta x}+ O(\Delta x^6) \]

• 用线性函数（\(u = x\)）检验差分格式

设：\(u(x) = x\)，则：\(\frac{du}{dx} = 1\)

• 网格表示

设：\(x_i = i\Delta x\)，则：\(u_i = x_i = i\Delta x\)

• 代入：

\[ \frac{-(i+3)\Delta x+ 9(i+2)\Delta x- 45(i+1)\Delta x + 45(i-1)\Delta x- 9(i-2)\Delta x+ (i-3)\Delta x}{60\Delta x}=-1 \]

2.1 离散化 ❾

六阶中心差分 ❌ 错误形式

\[ \left(\frac{\partial u}{\partial x}\right)_i= \frac{-u_{i+3} + 9u_{i+2} - 45u_{i+1} + 45u_{i-1} - 9u_{i-2} + u_{i-3}}{60\Delta x}+ O(\Delta x^6) \]

六阶中心差分

\[ \left(\frac{\partial u}{\partial x}\right)_i=\frac{-u_{i-3}+9u_{i-2}-45u_{i-1}+45u_{i+1}-9u_{i+2}+u_{i+3}}{60\Delta x}+ O(\Delta x^6) \]

精度 ↑ ⇒ stencil变宽 ⇒ 需要更多点 ⇒ 对边界更敏感

2.1 离散化 ❿

二阶导数差分格式

二阶精度（3点中心差分）

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

四阶精度（5点中心差分）

\[ \left(\frac{\partial^2 u}{\partial x^2}\right)_i= \frac{-u_{i+2} + 16u_{i+1} - 30u_i + 16u_{i-1} - u_{i-2}}{12(\Delta x)^2}+ O(\Delta x^4) \]

六阶精度（7点中心差分）

\[ \left(\frac{\partial^2 u}{\partial x^2}\right)_i= \frac{2u_{i+3} - 27u_{i+2} + 270u_{i+1} - 490u_i + 270u_{i-1} - 27u_{i-2} + 2u_{i-3}}{180(\Delta x)^2}+ O(\Delta x^6) \]

2.1 离散化 ⓫

二阶导数差分格式

二阶精度前向差分（4点）

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

四阶精度前向差分（6点）

\[ \left(\frac{\partial^2 u}{\partial x^2}\right)_i= \frac{45u_i - 154u_{i+1} + 214u_{i+2} - 156u_{i+3} + 61u_{i+4} - 10u_{i+5}}{12(\Delta x)^2}+ O(\Delta x^4) \]

混合导数

\[ \frac{\partial^2 u}{\partial x \partial y}\approx \frac{u_{i+1,j+1}-u_{i+1,j-1}-u_{i-1,j+1}+u_{i-1,j-1}}{4\Delta x \Delta y}+O(\Delta x^2,\Delta y^2) \]

2.1 离散化 ⓬

紧致差分格式（Compact Scheme）

紧致格式 = 用更少的网格点，获得更高精度的差分格式 👉 “用隐式关系换精度”

为什么需要紧致格式？ 回顾普通中心差分：

问题：

• 精度越高 → stencil越宽
  
• 边界难处理
  
• 并行效率下降

2.1 离散化 ⓭

紧致差分格式（Compact Scheme）

普通差分（显式）：

\[ u'_i = \frac{u_{i+1} - u_{i-1}}{2\Delta x} \]

紧致格式（隐式）：

\[ a u'_{i-1} + u'_i + a u'_{i+1} = \frac{b(u_{i+1} - u_{i-1})}{2\Delta x} \]

👉 导数之间也建立关系！ 例：四阶紧致格式

\[ \frac{1}{4}u'_{i-1} + u'_i + \frac{1}{4}u'_{i+1}= \frac{3}{4}\frac{u_{i+1} - u_{i-1}}{\Delta x}+O(\Delta x^4) \]

• 只用 3个点（i-1, i, i+1） 却达到 四阶精度

2.1 离散化 ⓮

紧致差分格式（Compact Scheme）

👉 对比：

本质：需要解线性方程组

把所有点写出来：

\[ \begin{cases} \frac{1}{4}u'_0 + u'_1 + \frac{1}{4}u'_2 = RHS_1 \\ \frac{1}{4}u'_1 + u'_2 + \frac{1}{4}u'_3 = RHS_2 \\ \cdots \end{cases} \]

👉 形成：

\[ A \mathbf{u'} = \mathbf{b} \]

2.1 离散化 ⓯

紧致差分格式（Compact Scheme）

👉 物理意义（可以这样讲）

紧致格式 = 更接近“谱方法”的差分方法

• 对波动问题（CFD、声学、稳定性）特别好
• 能更准确捕捉高频信息

2.1 离散化 ⓰

紧致差分格式（Compact Scheme）

⚠️ 常见误区

❌ 认为紧致格式更简单
→ 实际上需要解线性系统

❌ 认为一定更快
→ 小规模可能慢，但精度优势明显

❌ 忽略边界处理
→ 边界需要专门设计紧致格式

紧致格式通过引入导数之间的耦合关系，在保持紧 stencil 的同时显著提高精度，是连接有限差分与谱方法的重要桥梁。

2.1 离散化 ⓱

需要特别注意的是，当控制方程以如下形式书写时，即使对于粘性流动，所需的也仅是一阶导数。

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

• 进一步地，\(F\)、\(G\) 和 \(H\) 中的一些项包含粘性应力以及热传导项。这些物理量依赖于速度或温度的梯度，而这些梯度本身也是一阶导数。因此，对于这些粘性项，同样可以使用一阶导数的有限差分格式进行离散。通过这种方式，可以避免使用二阶导数的有限差分格式！

• 我们这里要强调的是：可以构造出几乎无限多种有限差分表达式，并且其精度可以不断提高。过去，在大多数计算流体力学（CFD）应用中，二阶精度通常被认为已经足够，因此本节中推导的这些差分格式一直是最常用的形式。关于高阶精度的优缺点如下：

• 缺点（con）：更高阶精度的差分格式通常需要更多的网格点，因此在每一个时间步或空间步中，会增加计算机的计算量。

• 优点（pro）：相比之下，高阶差分格式可能只需要更少的总体网格点，就可以获得相当的整体精度。高阶差分格式通常可以得到“质量更高”的数值解，例如捕捉到更清晰、更尖锐的激波结构。这一方面目前仍是CFD领域的研究热点。

基于以上原因，对于不同CFD问题究竟应采用何种精度，并不存在一个绝对明确的答案。二阶精度在绝大多数CFD应用中已被广泛接受。

2.1 离散化 ⓲

边界如何处理？

在内部点，我们可以用中心差分：

\[ \frac{\partial u}{\partial y}\Big|_i \approx \frac{u_{i+1}-u_{i-1}}{2\Delta y} \]

但在边界：没有 \(i-1\)，无法使用中心差分。👉 思路：单边差分（One-sided difference）

一阶精度（前向差分）

\[ \left(\frac{\partial u}{\partial y}\right)_1= \frac{u_2 - u_1}{\Delta y}+ O(\Delta y) \]

二阶精度单边差分（重点）

\[ \left(\frac{\partial u}{\partial y}\right)_1= \frac{-3u_1 + 4u_2 - u_3}{2\Delta y}+ O(\Delta y^2) \]

边界处理质量 = CFD精度关键来源之一

2.1 离散化 ⓳

边界如何处理？

在粘性流动中：

壁面剪切应力：

\[ \tau_w = \mu \left(\frac{\partial u}{\partial y}\right)_w \]

壁面热流：

\[ q_w = k \left(\frac{\partial T}{\partial y}\right)_w \]

👉 导数越准确 → 剪切应力 / 热流越准确

2.1 离散化 ⓴

边界如何处理？

2.2 显式/隐式格式 ❶

一致性问题

考虑一维热传导方程 \(\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}\)

时间（前向差分）： \[ \frac{\partial T}{\partial t} =\frac{T_i^{n+1}-T_i^n}{\Delta t}+O(\Delta t) \]

空间（中心差分）： \[ \frac{\partial^2 T}{\partial x^2}=\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta x)^2}+O((\Delta x)^2) \]

得到差分方程：

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

2.2 显式/隐式格式 ❷

一致性问题

⚠️ 重要概念 差分方程 ≠ 原微分方程

• 是近似
• 存在截断误差

定义一致性（Consistency）

当： \[ \Delta x \to 0,\quad \Delta t \to 0 \]

👉 差分方程 → 原方程

2.2 显式/隐式格式 ❸

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

显式方法（Explicit）

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

隐式方法（Implicit）

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

👉 Crank–Nicolson格式

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

2.2 显式/隐式格式 ❹

👉 Crank–Nicolson格式

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

隐式方法未知量：\(T_{i-1}^{n+1}\)、\(T_i^{n+1}\)、\(T_{i+1}^{n+1}\)

👉 一个方程多个未知数、不能逐点计算、必须：同时求解整个网格系统

\[ A T_{i-1}^{n+1} - B T_i^{n+1} + A T_{i+1}^{n+1} = K_i \]

其中

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

2.2 显式/隐式格式 ❺

\[ A T_{i-1}^{n+1} - B T_i^{n+1} + A T_{i+1}^{n+1} = K_i \]

若考虑一维区间上共有 \(7\) 个网格点，其中边界点 \(1\) 和 \(7\) 已知，则未知量为

\[ T_2^{n+1},\; T_3^{n+1},\; T_4^{n+1},\; T_5^{n+1},\; T_6^{n+1} \]

于是线性方程组可写成：

\[ \begin{aligned} -BT_2^{n+1}+AT_3^{n+1} &= K_2-AT_1^{n+1},\\ AT_2^{n+1}-BT_3^{n+1}+AT_4^{n+1} &= K_3,\\ AT_3^{n+1}-BT_4^{n+1}+AT_5^{n+1} &= K_4,\\ AT_4^{n+1}-BT_5^{n+1}+AT_6^{n+1} &= K_5,\\ AT_5^{n+1}-BT_6^{n+1} &= K_6-AT_7^{n+1}. \end{aligned} \]

2.2 显式/隐式格式 ❻

写成矩阵形式就是

\[ \begin{bmatrix} -B & A & 0 & 0 & 0 \\ A & -B & A & 0 & 0 \\ 0 & A & -B & A & 0 \\ 0 & 0 & A & -B & A \\ 0 & 0 & 0 & A & -B \end{bmatrix} \begin{bmatrix} T_2^{n+1}\\ T_3^{n+1}\\ T_4^{n+1}\\ T_5^{n+1}\\ T_6^{n+1} \end{bmatrix} = \begin{bmatrix} K_2-AT_1^{n+1}\\ K_3\\ K_4\\ K_5\\ K_6-AT_7^{n+1} \end{bmatrix} \]

也可将右端向量记为

\[ \mathbf{K}= \begin{bmatrix} K_2'=K_2-AT_1^{n+1}\\ K_3\\ K_4\\ K_5\\ K_6'=K_6-AT_7^{n+1} \end{bmatrix} \]

2.2 显式/隐式格式 ❼

因此最终得到标准形式

\[ \mathbf{A}\mathbf{T}=\mathbf{K} \]

其中

\[ \mathbf{A}= \begin{bmatrix} -B & A & 0 & 0 & 0 \\ A & -B & A & 0 & 0 \\ 0 & A & -B & A & 0 \\ 0 & 0 & A & -B & A \\ 0 & 0 & 0 & A & -B \end{bmatrix}, \qquad \mathbf{T}= \begin{bmatrix} T_2^{n+1}\\ T_3^{n+1}\\ T_4^{n+1}\\ T_5^{n+1}\\ T_6^{n+1} \end{bmatrix}, \qquad \mathbf{K}= \begin{bmatrix} K_2'\\ K_3\\ K_4\\ K_5\\ K_6' \end{bmatrix} \]

说明

• 系数矩阵 \(\mathbf{A}\) 只在主对角线和两条相邻对角线上非零，因此称为三对角矩阵。
• 边界点 \(T_1^{n+1}\) 和 \(T_7^{n+1}\) 已知，因此被移入右端向量。
• 这类方程组通常可用 Thomas 算法（追赶法） 高效求解。

2.2 显式/隐式格式 ❽

显式 vs 隐式

时间推进（Time Marching）：从已知时间层 \(n\) 推进到 \(n+1\)

1. 稳态问题
   • 通过时间推进逼近
2. 非稳态问题
   • 必须时间推进

2.2 显式/隐式格式 ❾

显式 vs 隐式

既然显式方法更简单，为什么不一直使用显式方法？

👉 答案：稳定性约束

在差分格式中，会出现：

• 空间步长：\(\Delta x\)、时间步长：\(\Delta t\)
• 一旦 \(\Delta x\) 给定、\(\Delta t\) 不能任意选择

必须满足稳定性条件（如 CFL 条件），否则会发生：

• 数值变为无穷大
• 出现非法运算（如负数开方）

程序崩溃

2.2 显式/隐式格式 ❿

显式 vs 隐式

一些隐式方法是 无条件稳定（unconditionally stable）

👉 不存在“最优方法”，只有：针对问题选择合适方法

2.2 显式/隐式格式 ⓫

显式 vs 隐式

1969 年至约 1979 年期间：

• 绝大多数涉及“推进计算（marching solutions）”的实际 CFD 计算
  
• 都采用 显式方法（explicit methods）

1980年代：

• 隐式方法成为 CFD 的主流研究方向

如今：

• 计算机体系结构发生了重要变化，对于大规模并行计算（Massively Parallel Computing）可以在同一时刻处理：成千上万个网格点
• 👉 因此：显式方法非常适合并行计算

2.2 误差与稳定性分析 ❶

误差的来源

1️⃣ 截断误差（Discretization Error）👉 一致性问题

\[ \text{Discretization error} = A - D \]

• \(A\)：微分方程精确解
• \(D\)：差分方程精确解

2️⃣ 舍入误差（Round-off Error）👉 稳定性问题

\[ \epsilon = N - D \]

• \(N\)：计算机数值解
• 计算机解 = 理论离散解 + 数值误差

👉 舍入误差如何随时间演化 = 稳定性本质

2.2 误差与稳定性分析 ❷

\[ \epsilon = N - D \]

• \(N\)：计算机数值解
• \(D\)：差分方程精确解
  • 若误差：不增长 / 衰减 → ✅ 稳定
  • 若误差：增长 → ❌ 不稳定

对于热传导方程：

\[ \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{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta x)^2} \]

2.2 误差与稳定性分析 ❸

显式格式：

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

考虑数值解与精确解的分解：

\[ T_i^n = D_i^n + \epsilon_i^n \]

其中：

• \(D_i^n\)：差分方程的精确解
• \(\epsilon_i^n\)：数值误差

代入差分格式

\[ \frac{D_i^{n+1} + \epsilon_i^{n+1} - D_i^n - \epsilon_i^n}{\alpha \Delta t}=\frac{D_{i+1}^n + \epsilon_{i+1}^n - 2D_i^n - 2\epsilon_i^n + D_{i-1}^n + \epsilon_{i-1}^n}{(\Delta x)^2} \]

\(D\) 是差分方程的精确解，因此满足：

\[ \frac{D_i^{n+1} - D_i^n}{\alpha \Delta t}=\frac{D_{i+1}^n - 2D_i^n + D_{i-1}^n}{(\Delta x)^2} \]

2.2 误差与稳定性分析 ❹

\[ \frac{D_i^{n+1} + \epsilon_i^{n+1} - D_i^n - \epsilon_i^n}{\alpha \Delta t}=\frac{D_{i+1}^n + \epsilon_{i+1}^n - 2D_i^n - 2\epsilon_i^n + D_{i-1}^n + \epsilon_{i-1}^n}{(\Delta x)^2} \]

\[ \frac{D_i^{n+1} - D_i^n}{\alpha \Delta t}=\frac{D_{i+1}^n - 2D_i^n + D_{i-1}^n}{(\Delta x)^2} \]

两式相减

\[ \frac{\epsilon_i^{n+1} - \epsilon_i^n}{\alpha \Delta t}=\frac{\epsilon_{i+1}^n - 2\epsilon_i^n + \epsilon_{i-1}^n}{(\Delta x)^2} \]

👉 误差 \(\epsilon\)：满足 与原方程完全相同形式的差分方程

👉 意味着：数值误差的传播规律 ≈ 物理量的传播规律

👉 von Neumann 稳定性分析的起点

2.2 误差与稳定性分析 ❺

von Neumann 稳定性分析（Fourier模态分析）

我们从误差满足的差分方程出发：

\[ \frac{\epsilon_i^{n+1} - \epsilon_i^n}{\Delta t}=\alpha \frac{\epsilon_{i+1}^n - 2\epsilon_i^n + \epsilon_{i-1}^n}{(\Delta x)^2} \]

误差的 Fourier 展开：

\[ \epsilon(x,t) = \sum_m A_m(t) e^{i k_m x} \]

由于方程是线性的：👉 每个模态可以独立分析

2.2 误差与稳定性分析 ❻

von Neumann 稳定性分析（Fourier模态分析）

单一波数：

\[ \epsilon_i^n = \hat{\epsilon}^n e^{i k x_i} \]

其中：

• \(x_i = i \Delta x\)
• \(\hat{\epsilon}^n\)：振幅（随时间变化）

进一步假设时间指数增长：

\[ \hat{\epsilon}^n = e^{a t_n} \]

因此：

\[ \epsilon_i^n = e^{a t_n} e^{i k x_i} \]

2.2 误差与稳定性分析 ❼

\[\epsilon_i^n = e^{a t_n} e^{i k x_i}\]

代入差分方程

\[ \frac{\epsilon_i^{n+1} - \epsilon_i^n}{\Delta t}=\alpha \frac{\epsilon_{i+1}^n - 2\epsilon_i^n + \epsilon_{i-1}^n}{(\Delta x)^2} \]

时间推进项

\[ \epsilon_i^{n+1} = e^{a (t_n+\Delta t)} e^{i k x_i}= e^{a \Delta t} \epsilon_i^n \]

空间项

\[ \epsilon_{i+1}^n = \epsilon_i^n e^{i k \Delta x} \]

\[ \epsilon_{i-1}^n = \epsilon_i^n e^{-i k \Delta x} \]

2.2 误差与稳定性分析 ❽

2.2 误差与稳定性分析 ❾

放大因子定义：

\[ G = \left|\frac{\epsilon^{n+1}}{\epsilon^n}\right|=\left|e^{a \Delta t}\right|=\left|1 - \frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2\left(\frac{k\Delta x}{2}\right)\right| \]

稳定性条件

\[ G \le 1 \]

条件 1

\[ 1 - \frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k_m \Delta x}{2} \right)\le 1\Longrightarrow\frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k_m \Delta x}{2} \right)\ge 0 (自动满足) \]

条件 2

\[ 1 - \frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k_m \Delta x}{2} \right)\ge -1\Longrightarrow\frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k_m \Delta x}{2} \right)\le 2 \]

2.2 误差与稳定性分析 ❿

\[ \frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k_m \Delta x}{2} \right)\le 2 \]

由于：

\[ \sin^2 \left( \frac{k_m \Delta x}{2} \right) \le 1 \]

👉 最不利情况：

\[ \sin^2(\cdot) = 1 \]

因此得到最终稳定性条件：

\[ \frac{4\alpha \Delta t}{(\Delta x)^2} \le 2 \Longrightarrow\boxed{\frac{\alpha \Delta t}{(\Delta x)^2} \le \frac{1}{2}} \]

2.2 误差与稳定性分析 ⓫

物理与数值意义

• 对于给定 \(\Delta x\)
• 时间步长 \(\Delta t\) 必须足够小 👉 CFL条件（Courant-Friedrichs-Lewy条件）

上述分析属于：

von Neumann 稳定性分析方法

特点：

• 基于 Fourier 模态（“误差波”）；稳定性 = 这些波是否被放大
• 适用于线性差分方程
• CFD 中最常用稳定性分析工具

2.2 误差与稳定性分析 ⓬

von Neumann 稳定性分析方法

👉 方法核心：

1. 误差满足同一差分方程
2. 用 Fourier 分解误差
3. 研究单一波数传播
4. 得到放大因子 \(G\)
5. 判断 \(G\le1\)

2.2 误差与稳定性分析 ⓭

双曲方程

考虑一阶波动方程：

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

时间离散（前向差分）、空间离散（中心差分）👉 Euler显式格式

\[ \frac{u_i^{n+1} - u_i^n}{\Delta t}=- c \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x} \]

对该格式进行 von Neumann 分析：

无论 \(\Delta t\) 取何值，解都是不稳定的

2.2 误差与稳定性分析 ⓮

引入 Fourier 模态

\[ u_i^n = \hat{u}^n e^{ikx_i} \]

其中：

• \(x_i = i\Delta x\)
• \(k\)：波数

时间推进项

\[ u_i^{n+1} = \hat{u}^{n+1} e^{ikx_i} \]

空间项

\[ u_{i+1}^n = \hat{u}^n e^{ik(x_i+\Delta x)} = u_i^n e^{ik\Delta x} \]

\[ u_{i-1}^n = u_i^n e^{-ik\Delta x} \]

2.2 误差与稳定性分析 ⓯

\[ \frac{u_i^{n+1} - u_i^n}{\Delta t}=- c \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x} \]

左边：

\[ \frac{u_i^{n+1} - u_i^n}{\Delta t}=\frac{u_i^n}{\Delta t}\left(\frac{\hat{u}^{n+1}}{\hat{u}^n} - 1\right) \]

右边：

\[ -c \frac{u_i^n(e^{ik\Delta x} - e^{-ik\Delta x})}{2\Delta x} \]

约去公共因子 \(u_i^n\)

\[ \frac{\hat{u}^{n+1}}{\hat{u}^n} =1- \frac{c\Delta t}{2\Delta x}\left(e^{ik\Delta x} - e^{-ik\Delta x}\right) \]

2.2 误差与稳定性分析 ⓰

放大因子

\[ G = \left|\frac{\hat{u}^{n+1}}{\hat{u}^n}\right|= \left|1- \frac{c\Delta t}{2\Delta x}\left(e^{ik\Delta x} - e^{-ik\Delta x}\right)\right| \]

利用指数恒等式

\[ e^{ik\Delta x} - e^{-ik\Delta x} = 2i\sin(k\Delta x) \]

得到：

\[ |G|^2 = 1 + C^2 \sin^2(k\Delta x)>1 （无条件不稳定） \]

其中：

\[ C = \frac{c\Delta t}{\Delta x} \]

👉 FTCS（前向时间 + 中心空间）格式对对流方程是无条件不稳定的

2.2 误差与稳定性分析 ⓱

双曲方程

FTCS格式：无条件不稳定

\[ \frac{u_i^{n+1} - u_i^n}{\Delta t}=- c \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x} \]

Lax 方法：用平均值替代

\[ u_i^n \rightarrow \frac{1}{2}(u_{i+1}^n + u_{i-1}^n) \]

代入后：

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

2.2 误差与稳定性分析 ⓲

双曲方程

设误差：

\[ \epsilon = e^{a t} e^{ikx} \]

代入可得：

\[ G=\left|e^{a\Delta t}\right|=\left|\cos(k\Delta x) - iC\sin(k\Delta x)\right| \]

其中：

\[ C = \frac{c\Delta t}{\Delta x} \]

稳定性条件要求：

\[ |e^{a\Delta t}| \le 1 \Longrightarrow C \le 1 \]

2.2 误差与稳定性分析 ⓳

对于二阶波动方程： \[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]

稳定性条件同样要求：

\[ C = \frac{c\Delta t}{\Delta x} \le 1 \]

两条特征线 \[ x = ct,\quad x = -ct \]

👉 \(C < 1\)（稳定）：数值域包含物理域，信息传播范围被覆盖，计算可靠

👉 \(C > 1\)（不稳定）：数值域小于物理域，信息传播超出计算能力，解不稳定

\[ \boxed{ \text{数值依赖域必须包含物理依赖域} } \]

2.2 误差与稳定性分析 ⓴

👉 CFD核心思想：

• 误差不可避免
• 关键是控制误差传播

思考：如果我们拥有一台没有舍入误差的理想计算机，那么数值计算中就不会出现不稳定性？

实际上，数值稳定性的基本概念是建立在解本身随时间的演化行为之上的；它并不本质依赖于舍入误差的行为。将解本身表示为傅里叶级数：

\[ U(x) = \sum_m V_m e^{ik_m x} \]

其中，\(V_m\) 是第 \(m\) 个谐波的振幅。相应地，放大因子可写为：

\[ G = \frac{V_m^{n+1}}{V_m^n} \]

2.2 误差与稳定性分析 ㉑

👉 回忆：以下内容是否掌握：

1. 误差的两种来源：数值误差 = 截断误差 + 舍入误差
2. 误差传播 = 稳定性本质
3. von Neumann 方法：Fourier 模态思想
4. 放大因子 \(G\)
5. CFL 条件的物理意义

2.3 网格与坐标变换 ❶

👉 如何构造“合适的网格”？

• 需要：
  • 工程直觉（物理行为）
  • 数学理解（函数性质）
  • 一定想象力
• 网格直接影响：
  • 数值精度
  • 收敛性
  • 物理边界刻画能力

2.3 网格与坐标变换 ❷

网格生成（Grid Generation）

两个本质问题

① 网格生成 👉 几何问题

② 方程求解 👉 数值问题

⚠️ 二者本质不同！

标准有限差分法要求：👉 均匀矩形网格

但现实问题：

• 复杂几何（如机翼）
• 非规则区域

👉 无法直接使用矩形网格

2.3 网格与坐标变换 ❸

示例：机翼流动问题

• 使用矩形网格的问题：
  1. 网格点落在机翼内部 ❌
     
  2. 机翼表面缺乏网格点 ❌
     
  3. 无法正确施加边界条件

思路：网格变换

👉 将物理空间中的非均匀网格 → 转换为计算空间中的规则网格

👉 使用 贴体网格（boundary-fitted grid）

• 网格线贴合物体边界
• 机翼表面成为一条坐标线

2.3 网格与坐标变换 ❹

曲线网格（物理空间）

• 非均匀
• 非矩形
• 难以直接差分计算

问题：如何进行数值离散？

👉 解决方法：坐标变换

\[ (x, y) \rightarrow (\xi, \eta) \]

在 \((\xi, \eta)\) 空间中：

• 网格变为：
  • 规则、矩形、均匀

👉 可直接使用有限差分

2.3 网格与坐标变换 ❺

关键要求：一一映射（One-to-one mapping）

• 物理空间 ↔︎ 计算空间
• 每个点唯一对应，例如：物理点 \(a,b,c\) 对应计算空间 \(a,b,c\)

数值求解流程

1. 物理空间坐标变换 → 计算空间
2. 控制方程从 \((x,y)\) 转换到 \((\xi,\eta)\)，在计算空间求解 PDE
3. 结果映射回物理空间

👉需要掌握

1. 方程坐标变换形式
2. 网格生成方法
3. 理解CFD中的网格哲学

2.3 网格与坐标变换 ❻

考虑二维非定常流：

\[ (x,y,t) \rightarrow (\xi,\eta,\tau) \]

定义：

\[ \begin{cases} \xi = \xi(x,y,t) \\ \eta = \eta(x,y,t) \\ \tau = \tau(t) \end{cases},~~~通常：\tau = t \]

链式法则

\[ \boxed{\frac{\partial}{\partial x}=\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial}{\partial \eta}\frac{\partial \eta}{\partial x}} ~~~~~~\boxed{ \frac{\partial}{\partial y}=\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial y}+\frac{\partial}{\partial \eta}\frac{\partial \eta}{\partial y}} ~~~~~~ \boxed{ \frac{\partial}{\partial t}=\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial t}+\frac{\partial}{\partial \eta}\frac{\partial \eta}{\partial t}+\frac{\partial}{\partial \tau} } \]

以下项称为度量系数（metrics）（反映了网格的几何特征）：

\[ \frac{\partial \xi}{\partial x},\; \frac{\partial \eta}{\partial x},\; \frac{\partial \xi}{\partial y},\; \frac{\partial \eta}{\partial y} \]

2.3 网格与坐标变换 ❼

链式法则（二阶导数）

设：

\[ A = \frac{\partial}{\partial x}=\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial}{\partial \eta}\frac{\partial \eta}{\partial x} \]

\[ \frac{\partial^2}{\partial x^2}=\frac{\partial}{\partial x}(A)=\frac{\partial}{\partial x}\left[\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial}{\partial \eta}\frac{\partial \eta}{\partial x}\right] \]

使用乘积法则与链式法则：

\[ \boxed{\begin{aligned}\frac{\partial^2}{\partial x^2}=&\frac{\partial}{\partial \xi}\frac{\partial^2 \xi}{\partial x^2}+\frac{\partial}{\partial \eta}\frac{\partial^2 \eta}{\partial x^2}\\&+\frac{\partial^2}{\partial \xi^2}\left(\frac{\partial \xi}{\partial x}\right)^2+\frac{\partial^2}{\partial \eta^2}\left(\frac{\partial \eta}{\partial x}\right)^2\\&+2\frac{\partial^2}{\partial \xi \partial \eta}\frac{\partial \xi}{\partial x}\frac{\partial \eta}{\partial x}\end{aligned}}~~~~~~~~~\boxed{\begin{aligned}\frac{\partial^2}{\partial y^2}=&\frac{\partial}{\partial \xi}\frac{\partial^2 \xi}{\partial y^2}+\frac{\partial}{\partial \eta}\frac{\partial^2 \eta}{\partial y^2}\\&+\frac{\partial^2}{\partial \xi^2}\left(\frac{\partial \xi}{\partial y}\right)^2+\frac{\partial^2}{\partial \eta^2}\left(\frac{\partial \eta}{\partial y}\right)^2\\&+2\frac{\partial^2}{\partial \xi \partial \eta}\frac{\partial \xi}{\partial y}\frac{\partial \eta}{\partial y}\end{aligned}} \]

2.3 网格与坐标变换 ❽

混合导数

\[ \boxed{\begin{aligned}\frac{\partial^2}{\partial x \partial y}=\frac{\partial}{\partial \xi}\frac{\partial^2 \xi}{\partial x \partial y}+\frac{\partial}{\partial \eta}\frac{\partial^2 \eta}{\partial x \partial y}+\frac{\partial^2}{\partial \xi^2}\frac{\partial \xi}{\partial x}\frac{\partial \xi}{\partial y}+\frac{\partial^2}{\partial \eta^2}\frac{\partial \eta}{\partial x}\frac{\partial \eta}{\partial y}+\frac{\partial^2}{\partial \xi \partial \eta}\left(\frac{\partial \xi}{\partial x}\frac{\partial \eta}{\partial y}+\frac{\partial \eta}{\partial x}\frac{\partial \xi}{\partial y}\right)\end{aligned}} \]

Laplace方程变换示例，原方程：

\[ \frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=0 \]

变换后形式

\[ \begin{aligned} \frac{\partial^2 \phi}{\partial \xi^2} \left[ \left(\frac{\partial \xi}{\partial x}\right)^2 + \left(\frac{\partial \xi}{\partial y}\right)^2 \right] + \frac{\partial^2 \phi}{\partial \eta^2} \left[ \left(\frac{\partial \eta}{\partial x}\right)^2 + \left(\frac{\partial \eta}{\partial y}\right)^2 \right] \\+ 2\frac{\partial^2 \phi}{\partial \xi \partial \eta} \left[ \frac{\partial \xi}{\partial x}\frac{\partial \eta}{\partial x} + \frac{\partial \xi}{\partial y}\frac{\partial \eta}{\partial y} \right] + \frac{\partial \phi}{\partial \xi} \left( \frac{\partial^2 \xi}{\partial x^2} + \frac{\partial^2 \xi}{\partial y^2} \right) + \frac{\partial \phi}{\partial \eta} \left( \frac{\partial^2 \eta}{\partial x^2} + \frac{\partial^2 \eta}{\partial y^2} \right) =0 \end{aligned} \]

2.3 网格与坐标变换 ❾

评注

• 计算在规则网格进行
• 方程复杂性增加
• 几何复杂性转移到metrics
• 若使用守恒形式：可避免复杂二阶导变换

\[ \frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}+\frac{\partial \mathbf{H}}{\partial z}= \mathbf{RHS} \]

2.3 网格与坐标变换 ❿

度量项系数（Metrics）

\[ \frac{\partial \xi}{\partial x},\ \frac{\partial \xi}{\partial y},\ \frac{\partial \eta}{\partial x},\ \frac{\partial \eta}{\partial y} \]

如果变换\((x,y,t) \rightarrow (\xi,\eta,\tau)\)是解析给定的，则可以直接得到这些度量项的解析表达式。然而，在许多 CFD 应用中，该变换是数值给定的，因此度量项通常通过有限差分来计算。

在许多应用中，变换更方便写为其逆变换形式：

\[ x = x(\xi, \eta, \tau),~~y = y(\xi, \eta, \tau),~~ t = t(\tau) \]

\(\xi, \eta, \tau\) 是自变量。 要从逆变换得到度量系数，需要建立如下关系： \[ \frac{\partial \xi}{\partial x},\ \frac{\partial \eta}{\partial y} \quad \Longleftrightarrow \quad \frac{\partial x}{\partial \xi},\ \frac{\partial y}{\partial \eta} \]

2.3 网格与坐标变换 ⓫

2.3 网格与坐标变换 ⓬

考虑流动变量 \(u = u(x, y)\)，\(x = x(\xi, \eta), \quad y = y(\xi, \eta)\)

对 \(\xi\) 求导：

\[\frac{\partial u}{\partial \xi} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \xi} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \xi}\]

对 \(\eta\) 求导：

\[\frac{\partial u}{\partial \eta} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \eta} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \eta}\]

看作关于未知量 \(\frac{\partial u}{\partial x}, \quad \frac{\partial u}{\partial y}\) 的线性方程组：

\[\begin{cases} \frac{\partial u}{\partial \xi} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \xi} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \xi} \\[6pt] \frac{\partial u}{\partial \eta} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \eta} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \eta} \end{cases}\]

2.3 网格与坐标变换 ⓭

用克拉默法则求解

\(\frac{\partial u}{\partial x} = \frac{ \begin{vmatrix} \frac{\partial u}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial u}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{vmatrix} }{ \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{vmatrix} }\)

分母定义为雅可比行列式（Jacobian）：

\[ J = \frac{\partial(x,y)}{\partial(\xi,\eta)} = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{vmatrix} \]

展开得：

\(\frac{\partial u}{\partial x} = \frac{1}{J} \left[ \left(\frac{\partial u}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right) - \left(\frac{\partial u}{\partial \eta}\right)\left(\frac{\partial y}{\partial \xi}\right) \right]\)

\(\frac{\partial u}{\partial y} = \frac{1}{J} \left[ \left(\frac{\partial u}{\partial \eta}\right)\left(\frac{\partial x}{\partial \xi}\right) - \left(\frac{\partial u}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right) \right]\)

2.3 网格与坐标变换 ⓮

将上述表达推广为算子形式：

\[ \frac{\partial}{\partial x} = \frac{1}{J} \left[ \left(\frac{\partial}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right) - \left(\frac{\partial}{\partial \eta}\right)\left(\frac{\partial y}{\partial \xi}\right) \right] \]

\[ \frac{\partial}{\partial y} = \frac{1}{J} \left[ \left(\frac{\partial}{\partial \eta}\right)\left(\frac{\partial x}{\partial \xi}\right) - \left(\frac{\partial}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right) \right] \]

推导度量系数表达式

\[ \frac{\partial \xi}{\partial x} = \frac{1}{J}\frac{\partial y}{\partial \eta},~~~~~~~~~ \frac{\partial \eta}{\partial x} = -\frac{1}{J}\frac{\partial y}{\partial \xi} \]

\[ \frac{\partial \xi}{\partial y} = -\frac{1}{J}\frac{\partial x}{\partial \eta},~~~~~~~~~ \frac{\partial \eta}{\partial y} = \frac{1}{J}\frac{\partial x}{\partial \xi} \]

2.3 网格与坐标变换 ⓯

评注

• 直接度量（direct metrics）

\[\frac{\partial \xi}{\partial x}, \frac{\partial \eta}{\partial x}\]

• 逆度量（inverse metrics）

\[\frac{\partial x}{\partial \xi}, \frac{\partial y}{\partial \eta}\]

• 该推导是 CFD 中网格使用的核心基础。

• 两者通过 Jacobian \(J\) 联系：

\[ J = \frac{\partial(x,y)}{\partial(\xi,\eta)} = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{vmatrix} = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta} - \frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi} \]

2.3 网格与坐标变换 ⓰

CFD控制方程的变换形式

对于二维非定常流动且无源项的情况，守恒型方程可写为：

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

（此处仅考虑二维\(x,y\)，是为了简化说明，推广到三维是直接的。）

我们希望得到如下形式：

\[\frac{\partial U_1}{\partial t}+\frac{\partial F_1}{\partial \xi}+\frac{\partial G_1}{\partial \eta}=0 \]

其中\(F_1,G_1\)是由原始通量\(F,G\)组合得到的新通量。

根据链式法则：

\[\frac{\partial U}{\partial t}+\frac{\partial F}{\partial \xi}\left(\frac{\partial \xi}{\partial x}\right)+\frac{\partial F}{\partial \eta}\left(\frac{\partial \eta}{\partial x}\right)+\frac{\partial G}{\partial \xi}\left(\frac{\partial \xi}{\partial y}\right)+\frac{\partial G}{\partial \eta}\left(\frac{\partial \eta}{\partial y}\right)=0 \]

2.3 网格与坐标变换 ⓱

两边乘以\(J\)：

\[J\frac{\partial U}{\partial t}+J\left(\frac{\partial F}{\partial \xi}\right)\left(\frac{\partial \xi}{\partial x}\right)+J\left(\frac{\partial F}{\partial \eta}\right)\left(\frac{\partial \eta}{\partial x}\right)+J\left(\frac{\partial G}{\partial \xi}\right)\left(\frac{\partial \xi}{\partial y}\right)+J\left(\frac{\partial G}{\partial \eta}\right)\left(\frac{\partial \eta}{\partial y}\right)=0 \]

考虑如下项：

\[J\left(\frac{\partial F}{\partial \xi}\right)\left(\frac{\partial \xi}{\partial x}\right)=\frac{\partial}{\partial \xi}\left[J F \left(\frac{\partial \xi}{\partial x}\right)\right]-F\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial x}\right)\]

类似地：

\[ \begin{aligned} J\left(\frac{\partial F}{\partial \eta}\right)\left(\frac{\partial \eta}{\partial x}\right) &= \frac{\partial}{\partial \eta}\left[J F \left(\frac{\partial \eta}{\partial x}\right)\right]-F\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial x}\right) \\ J\left(\frac{\partial G}{\partial \xi}\right)\left(\frac{\partial \xi}{\partial y}\right) &= \frac{\partial}{\partial \xi}\left[J G \left(\frac{\partial \xi}{\partial y}\right)\right]-G\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial y}\right) \\ J\left(\frac{\partial G}{\partial \eta}\right)\left(\frac{\partial \eta}{\partial y}\right) &= \frac{\partial}{\partial \eta}\left[J G \left(\frac{\partial \eta}{\partial y}\right)\right]-G\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial y}\right) \end{aligned} \]

2.3 网格与坐标变换 ⓲

导入到守恒型方程，整理得到：

\[ \begin{aligned} J\frac{\partial U}{\partial t}+\frac{\partial}{\partial \xi}\left(JF\frac{\partial \xi}{\partial x}+JG\frac{\partial \xi}{\partial y}\right)+\frac{\partial}{\partial \eta}\left(JF\frac{\partial \eta}{\partial x}+JG\frac{\partial \eta}{\partial y}\right)\\-F\left[\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial x}\right)+\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial x}\right)\right]-G\left[\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial y}\right)+\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial y}\right)\right]&=0 \end{aligned} \]

注意恒等式

\[\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial x}\right)+\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial x}\right)=\frac{\partial}{\partial \xi}\left(\frac{\partial y}{\partial \eta}\right)-\frac{\partial}{\partial \eta}\left(\frac{\partial y}{\partial \xi}\right)=\frac{\partial^2 y}{\partial \xi \partial \eta}-\frac{\partial^2 y}{\partial \eta \partial \xi}=0\]

\[\frac{\partial}{\partial \xi}\left(J\frac{\partial \xi}{\partial y}\right)+\frac{\partial}{\partial \eta}\left(J\frac{\partial \eta}{\partial y}\right)=\frac{\partial}{\partial \xi}\left(-\frac{\partial x}{\partial \eta}\right)+\frac{\partial}{\partial \eta}\left(\frac{\partial x}{\partial \xi}\right)=-\frac{\partial^2 x}{\partial \xi \partial \eta}+\frac{\partial^2 x}{\partial \eta \partial \xi}=0\]

2.3 网格与坐标变换 ⓳

最终守恒形式化简为：

\[ \frac{\partial U_1}{\partial t}+\frac{\partial F_1}{\partial \xi}+\frac{\partial G_1}{\partial \eta}=0 \]

其中：

\[\begin{aligned} U_1 &= JU \\ F_1 &= JF\frac{\partial \xi}{\partial x}+JG\frac{\partial \xi}{\partial y} = F\frac{\partial y}{\partial \eta}-G\frac{\partial x}{\partial \eta}\\ G_1 &= \underbrace{JF\frac{\partial \eta}{\partial x}+JG\frac{\partial \eta}{\partial y}}_{\text{度量系数表示}} = \underbrace{-F\frac{\partial y}{\partial \xi}+G\frac{\partial x}{\partial \xi}}_{\text{逆度量系数表示}} \end{aligned}\]

2.3 网格与坐标变换 ⓴

小结

这一节的关键是：

✔ 将控制方程保持为强守恒形式
✔ 引入新的守恒变量：\(U_1=JU\)
✔ 构造新的通量：\(F_1,G_1\)
✔ 利用Jacobian实现几何映射

本质上：

坐标变换 ≠ 改变物理规律
而是改变“表达方式”以适应计算

2.4 网格变换实例 ❶

例：平板附近粘性流动的网格需求

考虑图示的物理平面与计算平面。

设我们研究的是：

• 平直壁面上的粘性流动；
• 且如图左侧速度分布示意所示，靠近壁面处速度变化很快。

因此，为了精细计算近壁区域的流动细节：

• 在物理平面中，
• 应在 \(y\) 方向布置更密的网格。

但是，在离壁面较远处，流场变化较缓，网格则可以较粗。

另一方面，我们又希望在计算平面 \((\xi,\eta)\) 中处理的是：均匀网格

2.4 网格变换实例 ❷

例：平板附近粘性流动的网格需求

实现这种网格拉伸的一个简单解析变换为

\[ \xi = x,~~~~~~\eta = \ln(y+1) \]

对应的逆变换（inverse transformation）为

\[ x=\xi,~~~~~~y=e^\eta-1 \]

可得

\[ \frac{dy}{d\eta}=e^\eta \]

• 当 \(\eta\) 变大时，即当我们在图中从平板向上远离壁面时，对于同一个固定的 \(\Delta\eta\)，物理空间中的 \(\Delta y\) 会越来越大。
• 虽然计算平面中是均匀网格，但在物理平面中，网格会在竖直方向逐渐变疏。

2.4 网格变换实例 ❸

例：控制方程在物理空间和计算空间中的形式

• 假定为定常流
• 仅以连续性方程为例说明

物理空间中

\[ \frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}=0 \]

利用链式法则，得到

\[ \frac{\partial (\rho u)}{\partial \xi}\left(\frac{\partial \xi}{\partial x}\right)+\frac{\partial (\rho u)}{\partial \eta}\left(\frac{\partial \eta}{\partial x}\right)+\frac{\partial (\rho v)}{\partial \xi}\left(\frac{\partial \xi}{\partial y}\right)+\frac{\partial (\rho v)}{\partial \eta}\left(\frac{\partial \eta}{\partial y}\right)=0 \]

由直接变换式

\[ \xi=x,\qquad \eta=\ln(y+1) \Rightarrow \frac{\partial \xi}{\partial x}=1, \qquad \frac{\partial \xi}{\partial y}=0, \qquad \frac{\partial \eta}{\partial x}=0, \qquad \frac{\partial \eta}{\partial y}=\frac{1}{y+1} \]

2.4 网格变换实例 ❹

代入后得到的连续性方程

\[ \frac{\partial (\rho u)}{\partial \xi}+\frac{1}{y+1}\frac{\partial (\rho v)}{\partial \eta}=0 \]

又因为根据式 \(y+1=e^\eta\)，计算空间的方程可写为

\[ e^\eta \frac{\partial (\rho u)}{\partial \xi}+\frac{\partial (\rho v)}{\partial \eta}=0 \]

此时应该能更清楚地看到：

• 网格变换并不是只改变点的位置；
• 它还会改变控制方程的数学形式；
• 这种改变正是通过度量系数体现出来的。

2.4 网格变换实例 ❺

从逆变换角度推导

\[ \frac{1}{J}\left[\frac{\partial (\rho u)}{\partial \xi}\left(\frac{\partial y}{\partial \eta}\right)-\frac{\partial (\rho u)}{\partial \eta}\left(\frac{\partial y}{\partial \xi}\right)\right]+\frac{1}{J}\left[\frac{\partial (\rho v)}{\partial \eta}\left(\frac{\partial x}{\partial \xi}\right)-\frac{\partial (\rho v)}{\partial \xi}\left(\frac{\partial x}{\partial \eta}\right)\right]=0 \]

其中 \(J\) 是雅可比。由逆变换 \(x=\xi,\qquad y=e^\eta-1\) 得到逆度量系数

\[ \frac{\partial x}{\partial \xi}=1, \qquad \frac{\partial x}{\partial \eta}=0, \qquad \frac{\partial y}{\partial \xi}=0, \qquad \frac{\partial y}{\partial \eta}=e^\eta \]

代入逆度量系数后的结果

\[ e^\eta \frac{\partial (\rho u)}{\partial \xi}+\frac{\partial (\rho v)}{\partial \eta}=0 \]

结果完全一致。

2.4 网格变换实例 ❻

评注

注意前面的推导过程：

• 连续性方程首先通过导数变换被改写；
• 在这一步，式仍然是一般形式；
• 只有当具体变换对应的度量系数代入之后，变换才真正变成“某一个具体网格”。

因此可以认识到：

控制方程变换中的具体信息，全部由度量系数来承载。

也就是说：

• 你不一定需要知道网格是怎样“画出来的”；
• 但你一定需要知道该网格对应的度量系数。

2.4 网格变换实例 ❼

评注

一个很形象的想象：假设你负责在某个物体外部进行数值计算，而另一个人负责网格生成。

那么：

• 你去向网格生成的人索取的，
• 其实就是这个坐标变换的度量系数（metrics）。

因为对于在计算平面中求解流动问题而言，

度量项已经包含了你所需的全部几何变换信息。

2.5 贴体坐标系：椭圆型方程用于网格生成 ❶

先看一个简单问题。考虑图示的扩张管道（divergent duct）内流动。

其中：

• 曲线 (de) 是管道上壁；
• 线段 (fg) 是中心线。

对于这个流动问题，在物理平面中采用简单矩形网格并不合适。

• 边界是曲线；
• 若仍用简单矩形网格，
• 网格线就无法准确贴合边界。

因此，需要构造曲线网格（curvilinear grid）。

2.5 贴体坐标系：椭圆型方程用于网格生成 ❷

绘制曲线网格，使得：

• 上壁 (de) 本身就是一条坐标线；
• 中心线 (fg) 也是一条坐标线。

这样，网格就能精确贴合边界。

一个简单的贴体变换：记图中上边界的纵坐标为

\[ y_s=f(x) \]

则如下变换可将其映射为 \((\xi,\eta)\) 空间中的矩形网格：

\[ \xi=x, \qquad \eta=\frac{y}{y_s},\qquad \text{其中 } y_s=f(x) \]

2.5 贴体坐标系：椭圆型方程用于网格生成 ❸

为什么这个变换是“贴体”的？

例如，在物理平面中考虑点 \(d\)，其位置满足

\[ y=y_d=y_s(x_d) \]

将其代入式 (5.66)，得

\[ \eta_d=\frac{y_d}{y_s}=\frac{y_s(x_d)}{y_s(x_d)}=1 \]

因此，在计算平面中，点 \(d\) 位于

\[ \eta=\eta_d=1 \]

这说明：上边界上的这个点被映射到了计算平面中的一条水平直线 \(\eta=1\) 上。

2.5 贴体坐标系：椭圆型方程用于网格生成 ❹

为什么这个变换是“贴体”的？

再看物理平面中的点 \(c\)，其纵坐标为

\[ y=y_c=y_s(x_c) \]

显然它在物理平面中的纵坐标与点 \(d\) 不同，也即

\[ y_c>y_d \]

但将 \(y_c\) 代入可得

\[ \eta_c=\frac{y_c}{y_s}=\frac{y_s(x_c)}{y_s(x_c)}=1 \]

所以，在计算平面中，点 \(c\) 也位于

\[ \eta=\eta_c=1 \]

2.5 贴体坐标系：椭圆型方程用于网格生成 ❺

2.5 贴体坐标系：椭圆型方程用于网格生成 ❻

2.5 贴体坐标系：椭圆型方程用于网格生成 ❼

\(\xi\) 方向坐标线

2.5 贴体坐标系：椭圆型方程用于网格生成 ❽

如何将物理平面中的曲线网格映射为计算平面中的均匀网格？

基本思路 将物理网格画在笛卡尔坐标 \((x,y)\) 上：

内边界：

\[ (x,y)\ \text{在}\ \Gamma_1\ \text{已知} \]

外边界：

\[ (x,y)\ \text{在}\ \Gamma_2\ \text{已知} \]

👉 这构成一个边值问题（boundary-value problem）

• 边界上所有点的 \((x,y)\) 都已知
• 需要求解内部点

2.5 贴体坐标系：椭圆型方程用于网格生成 ❾

采用椭圆型偏微分方程定义坐标变换：

Laplace 方程

\[ \frac{\partial^2 \xi}{\partial x^2} + \frac{\partial^2 \xi}{\partial y^2} = 0 \]

\[ \frac{\partial^2 \eta}{\partial x^2} + \frac{\partial^2 \eta}{\partial y^2} = 0 \]

反变换（求 \(x,y\)）\(x = x(\xi,\eta), \quad y = y(\xi,\eta)\)

得到：

\[ \alpha \frac{\partial^2 x}{\partial \xi^2}- 2\beta \frac{\partial^2 x}{\partial \xi \partial \eta}+ \gamma \frac{\partial^2 x}{\partial \eta^2} = 0 \]

\[ \alpha \frac{\partial^2 y}{\partial \xi^2}- 2\beta \frac{\partial^2 y}{\partial \xi \partial \eta}+ \gamma \frac{\partial^2 y}{\partial \eta^2} = 0 \]

2.5 贴体坐标系：椭圆型方程用于网格生成 ❿

系数定义

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

\[ \beta = \frac{\partial x}{\partial \xi}\frac{\partial x}{\partial \eta}+ \frac{\partial y}{\partial \xi}\frac{\partial y}{\partial \eta} \]

\[ \gamma = \left(\frac{\partial x}{\partial \xi}\right)^2 + \left(\frac{\partial y}{\partial \xi}\right)^2 \]

为了保证问题可解：👉 必须在整个边界上给定条件

包括：

• \(\Gamma_1\)（内边界）
• \(\Gamma_2\)（外边界）
• \(\Gamma_3, \Gamma_4\)（切割产生的侧边界）

2.5 贴体坐标系：椭圆型方程用于网格生成 ⓫

C型网格切割

通过在尾缘切开：

• 引入两条新边界：\(\Gamma_3\), \(\Gamma_4\)
• 物理上重合，计算上分离

👉 形成矩形计算域

• 计算平面：规则矩形网格（均匀）
• 物理平面：非均匀曲线网格

对应关系：

\[ (\xi_i, \eta_j) \leftrightarrow (x_i, y_j) \]

2.5 贴体坐标系：椭圆型方程用于网格生成 ⓬

数值求解方法

椭圆型 PDE常用方法：

• 有限差分
• 松弛法（relaxation）

度量系数：👉 用中心差分计算

\[ \left(\frac{\partial \xi}{\partial x}\right)_{i,j}=\frac{\xi_{i+1,j} - \xi_{i-1,j}}{x_{i+1,j} - x_{i-1,j}} \]

2.5 贴体坐标系：椭圆型方程用于网格生成 ⓭

评注

1. 椭圆方程用于生成平滑网格
   
2. 需要完整边界条件
   
3. 计算平面规则 → 物理平面非规则
   
4. 网格质量由 PDE 控制
   
5. 远场边界需足够远（特别是亚声速流）
6. 近壁区域需加密网格
7. 网格生成与流动物理无关、仅是坐标映射、Laplace 方程只是工具

2.6 自适应网格 ❶

在流场中：

• 梯度大的区域 → 需要密集网格
• 梯度小 → 可稀疏

平板边界层示例

边界层厚度：

\[ \delta = \delta(x) \]

若网格均匀：

• 第一层网格点满足 \(\Delta y > \delta\)
• ❌ 无法解析边界层

2.6 自适应网格 ❷

优势

1. 在固定网格点数下提高精度
   
2. 在给定精度下降低网格数量

2.6 自适应网格 ❸

一个典型自适应变换，网格步长定义为：

\[ \Delta x = \frac{B\,\Delta \xi}{1 + b\left(\frac{\partial g}{\partial x}\right)} \]

\[ \Delta y = \frac{C\,\Delta \eta}{1 + c\left(\frac{\partial g}{\partial y}\right)} \]

其中：

• \(g\)：流场变量（如 \(p,\ \rho,\ T\)）
• \(b,c\)：调节参数
• \(B,C\)：尺度因子
• 梯度大 → 分母大 → \(\Delta x,\Delta y\) 小 → 网格加密
• 梯度小 → 网格变稀

2.6 自适应网格 ❹

网格随时间演化

由于：

\[ \frac{\partial g}{\partial x},\quad \frac{\partial g}{\partial y} \]

随时间变化：👉 网格点位置也随时间变化

设第 \(N\) 个网格点：

\[ x_N^t = \sum_{i=1}^{N} (\Delta x)_i^t \]

下一时刻：

\[ x_N^{t+\Delta t} = \sum_{i=1}^{N} (\Delta x)_i^{t+\Delta t} \]

2.6 自适应网格 ❺

网格随时间演化

类似地：

\[ y_M^t = \sum_{i=1}^{M} (\Delta y)_i^t \]

\[ y_M^{t+\Delta t} = \sum_{i=1}^{M} (\Delta y)_i^{t+\Delta t} \]

时间导数（动网格）

\[ \left(\frac{\partial x}{\partial t}\right)_{\xi,\eta}=\frac{x_N^{t+\Delta t} - x_N^t}{\Delta t} \]

\[ \left(\frac{\partial y}{\partial t}\right)_{\xi,\eta}=\frac{y_M^{t+\Delta t} - y_M^t}{\Delta t} \]

2.6 自适应网格 ❻

坐标变换关系

设：

\[ x = x(\xi,\eta,t), \qquad y = y(\xi,\eta,t) \]

则：

\[ dx = \left(\frac{\partial x}{\partial \xi}\right)d\xi+ \left(\frac{\partial x}{\partial \eta}\right)d\eta+ \left(\frac{\partial x}{\partial t}\right)dt, \qquad dy = \left(\frac{\partial y}{\partial \xi}\right)d\xi+ \left(\frac{\partial y}{\partial \eta}\right)d\eta+ \left(\frac{\partial y}{\partial t}\right)dt \]

时间度量项

\[ \frac{\partial \xi}{\partial t} \neq 0, \quad \frac{\partial \eta}{\partial t} \neq 0 \]

👉 与固定网格不同！

2.6 自适应网格 ❼

推导结果

\[ \frac{\partial \xi}{\partial t}=-\frac{1}{J}\left[\left(\frac{\partial x}{\partial t}\right)\left(\frac{\partial y}{\partial \eta}\right)-\left(\frac{\partial y}{\partial t}\right)\left(\frac{\partial x}{\partial \eta}\right)\right] \]

\[ \frac{\partial \eta}{\partial t}=\frac{1}{J}\left[\left(\frac{\partial x}{\partial t}\right)\left(\frac{\partial y}{\partial \xi}\right)-\left(\frac{\partial y}{\partial t}\right)\left(\frac{\partial x}{\partial \xi}\right)\right] \]

其中：

\[ J = \frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta}-\frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi} \]

👉 这意味着：控制方程中必须包含网格运动项

2.6 自适应网格 ❽

示例：后向台阶流

• 网格自动聚集在：
  • 激波
  • 剪切层
  • 分离区
• 自适应网格仍在发展中
  • 有多种不同实现方法
    
  • 现代 CFD 广泛使用

2.7 网格生成的现代发展与网格生成软件介绍 ❶

网格生成是 CFD 中：非常活跃的研究方向 例如：

• F-20 战斗机：几何复杂、需要多块网格（multi-block）

• 现代方法通常结合：椭圆网格、自适应网格、多块结构

• 典型特征：网格在高梯度区域高度聚集、表面附近分辨率高、远场较稀疏

网格技术发展方向：

1. 自适应
   
2. 高质量
   
3. 多块结构
   
4. 自动化

2.7 网格生成的现代发展与网格生成软件介绍 ❷

在计算流体力学（CFD）中，网格质量直接决定计算精度与稳定性

ANSYS Meshing

• 集成在 ANSYS Workbench 中
• 操作简单，适合入门
• 支持：
  • 非结构网格（四面体、六面体混合）
  • 边界层加密（Inflation）

👉 适合： - 初学者 - 常规工程问题

2.7 网格生成的现代发展与网格生成软件介绍 ❸

ICEM CFD

• 工业级网格工具
• 支持：
  • 结构网格（Hexa）
  • 非结构网格（Tet）
  • O-grid / C-grid
• Block结构非常强大（适合翼型、叶片）

👉 优点： - 网格质量高 - 可控性强

👉 缺点： - 学习曲线陡

👉 适合： - 高质量CFD计算 - 转捩 / 边界层问题

2.7 网格生成的现代发展与网格生成软件介绍 ❹

Pointwise

• 被称为“网格生成界的 CAD”
• 主要特点：
  • 高质量结构/混合网格
  • 边界拟合能力强
  • 网格光滑度极高

👉 优点： - 非常适合： - 高超声速 - 转捩问题 - DNS / LES

👉 缺点： - 商业授权较贵

2.7 网格生成的现代发展与网格生成软件介绍 ❺

Star-CCM+（Simcenter）

• Siemens 出品
• 一体化 CFD 平台（含网格）

👉 特点： - 自动网格（polyhedral） - 对复杂几何非常友好

👉 优点： - 自动化程度高 - 工程效率高

👉 缺点： - 网格控制不如 ICEM 精细

2.7 网格生成的现代发展与网格生成软件介绍 ❻

Fluent Meshing

• ANSYS Fluent 内置网格模块

👉 特点： - 自动化网格生成 - Mosaic / Poly 网格

👉 适合： - 工程快速计算 - 复杂几何

2.7 网格生成的现代发展与网格生成软件介绍 ❼

Gmsh（最推荐开源工具）

• 轻量级
• 支持：
  • 结构 / 非结构网格
  • 脚本化建模（.geo）

👉 优点： - 免费 - 可编程 - 适合科研

👉 缺点： - GUI较简单 - 大规模问题略弱

2.7 网格生成的现代发展与网格生成软件介绍 ❽

OpenFOAM 自带工具

blockMesh

• 结构网格
• 完全可控

snappyHexMesh

• 非结构网格
• 支持 STL 几何

👉 优点： - 完全开源 - 与求解器无缝耦合

👉 缺点： - 调参复杂 - 学习成本高

2.7 网格生成的现代发展与网格生成软件介绍 ❾

GridPro

• 专注结构网格
• 航空航天常用

👉 特点： - 高质量 O/C 网格 - 多块结构

NUMECA (Autogrid)

• 主要用于：
  • 涡轮机械
  • 叶片网格

2.7 网格生成的现代发展与网格生成软件介绍 ❿