流体力学核心贯通课 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)\)

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 离散化 ❷ 有限差分格式

从泰勒展开到有限差分

向前差分（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 离散化 ❸ 紧致差分格式

回顾普通中心差分：

2.1 离散化 ❸ 紧致差分格式

2.1 离散化 ❸ 紧致差分格式

下面给出其推导过程

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

在 \(x_i\) 处对两侧各项作 Taylor 展开，记步长为 \(\Delta x\)。

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

\[ u'_{i+1}= u'_i + \Delta x\,u''_i + \frac{\Delta x^2}{2}u'''_i+ \frac{\Delta x^3}{6}u^{(4)}_i+ \frac{\Delta x^4}{24}u^{(5)}_i + O(\Delta x^5) \]

因此，

\[ u'_{i-1}+u'_{i+1} = 2u'_i + \Delta x^2 u'''_i + \frac{\Delta x^4}{12}u^{(5)}_i + O(\Delta x^6) \]

2.1 离散化 ❸ 紧致差分格式

再乘上 \(\frac14\)，得到

\[ \frac14 u'_{i-1} + \frac14 u'_{i+1} = \frac12 u'_i + \frac14 \Delta x^2 u'''_i + \frac{1}{48}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

于是左端为

\[ \frac14 u'_{i-1} + u'_i + \frac14 u'_{i+1} = \frac32 u'_i + \frac14 \Delta x^2 u'''_i + \frac{1}{48}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

展开右端项 \(u_{i+1}\) 与 \(u_{i-1}\)：

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

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

2.1 离散化 ❸ 紧致差分格式

两式相减得

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

两边除以 \(\Delta x\)：

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

再乘上 \(\frac34\)，得到右端

\[ \frac34 \frac{u_{i+1}-u_{i-1}}{\Delta x} = \frac32 u'_i + \frac14 \Delta x^2 u'''_i + \frac{1}{80}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

2.1 离散化 ❸ 紧致差分格式

比较左右两端。左端展开为

\[ \frac14 u'_{i-1} + u'_i + \frac14 u'_{i+1}= \frac32 u'_i + \frac14 \Delta x^2 u'''_i + \frac{1}{48}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

右端展开为

\[ \frac34 \frac{u_{i+1}-u_{i-1}}{\Delta x}= \frac32 u'_i + \frac14 \Delta x^2 u'''_i + \frac{1}{80}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

可见两端在 \(u'_i\) 和 \(\Delta x^2 u'''_i\) 项上完全一致，故其差为

\[ \left(\frac{1}{48}-\frac{1}{80}\right)\Delta x^4 u^{(5)}_i + O(\Delta x^6)= \frac{1}{120}\Delta x^4 u^{(5)}_i + O(\Delta x^6) \]

因此该格式的截断误差为

\[ O(\Delta x^4) \]

即得到四阶精度格式

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

2.1 离散化 ❸ 紧致差分格式

评注

紧致差分格式（compact finite difference scheme）的特点：

1. 左端包含未知导数 \(u'_{i-1},u'_i,u'_{i+1}\)；
2. 右端只用到相邻点函数值 \(u_{i-1},u_{i+1}\)；
3. 只用三点模板就能达到四阶精度，因此比普通显式差分格式更“紧致”。

若写成求解形式，则对所有内部点可组成一个三对角线性方程组，再求得全部 \(u'_i\)。

2.1 离散化 ❸ 紧致差分格式

本质：需要解线性方程组

把所有点写出来：

\[ \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 离散化 ❸ 紧致差分格式

给定四阶紧致差分格式：

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

\[ \begin{cases} \frac{1}{4}u'_0 + u'_1 + \frac{1}{4}u'_2= \frac{3}{4}\frac{u_2 - u_0}{\Delta x} \\[8pt] \frac{1}{4}u'_1 + u'_2 + \frac{1}{4}u'_3= \frac{3}{4}\frac{u_3 - u_1}{\Delta x} \\[8pt] \frac{1}{4}u'_2 + u'_3 + \frac{1}{4}u'_4= \frac{3}{4}\frac{u_4 - u_2}{\Delta x} \\[6pt] \qquad \vdots \\[6pt] \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} \\[6pt] \qquad \vdots \end{cases} \]

2.1 离散化 ❸ 紧致差分格式

写成统一矩阵形式（便于数值实现）

\[ \begin{bmatrix} 1 & \tfrac{1}{4} & 0 & \cdots & 0 \\ \tfrac{1}{4} & 1 & \tfrac{1}{4} & \cdots & 0 \\ 0 & \tfrac{1}{4} & 1 & \tfrac{1}{4} & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \tfrac{1}{4} & 1 \end{bmatrix} \begin{bmatrix} u'_1 \\ u'_2 \\ u'_3 \\ \vdots \\ u'_{N-1} \end{bmatrix}= \begin{bmatrix} \frac{3}{4}\frac{u_2-u_0}{\Delta x} \\ \frac{3}{4}\frac{u_3-u_1}{\Delta x} \\ \frac{3}{4}\frac{u_4-u_2}{\Delta x} \\ \vdots \\ \frac{3}{4}\frac{u_N-u_{N-2}}{\Delta x} \end{bmatrix} \]

• RHS 是显式中心差分（但带权系数）
• LHS 是隐式三对角耦合
• 常用于DNS/LES

2.1 离散化 ❸ 紧致差分格式

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

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

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

2.1 离散化 ❸ 紧致差分格式

⚠️ 常见误区

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

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

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

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

2.1 离散化 ❸ 紧致差分格式

四阶精度的二阶导数中心紧致格式

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

其中：

• 左侧包含相邻节点的二阶导数项，因此这是一个隐式紧致格式
• 右侧采用三点中心差分形式
• 该格式整体达到 四阶精度

对内部节点 \(i=2,\dots,N-1\)，可写成

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

2.1 离散化 ❸ 紧致差分格式

其中矩阵 \(A\) 为三对角矩阵：

\[ A= \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ \frac{1}{10} & 1 & \frac{1}{10} & \cdots & 0 \\ 0 & \frac{1}{10} & 1 & \frac{1}{10} & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 1 \end{pmatrix} \]

右端项内部节点为

\[ \mathbf{rhs}_i= \frac{6}{5}\frac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2}, \qquad i=2,\dots,N-1 \]

边界点通常需要结合边界格式单独处理。

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

2.1 离散化 ❹ 讨论与小结

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.1 离散化 ❹ 讨论与小结

本节应掌握的知识点

1. Taylor展开
2. 基本差分格式推导 （所有模板点Taylor展开的组合）
3. 紧致格式（用隐式关系换精度）
4. 边界处理（在边界采用单边差分）

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

2.2 显式/隐式格式 ❶ 差分方程与微分方程的一致性

2.2 显式/隐式格式 ❷ 求解思路

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

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

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 隐式

2.2 显式/隐式格式 ❷ 求解思路

显式 vs 隐式

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

1. 稳态问题 通过时间推进逼近

2. 非稳态问题 必须时间推进

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

2.2 显式/隐式格式 ❷ 求解思路

显式 vs 隐式

1969 年至约 1979 年期间：

• 绝大多数实际 CFD 计算都采用 显式方法（explicit methods）

1980年代：

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

如今：

• 计算机体系结构发生了重要变化，对于大规模并行计算（Massively Parallel Computing）可以在同一时刻处理：成千上万个网格点
• 👉 显式方法重新受到重视：
  • ✔ 局部计算（local stencil）
  • ✔ 无需求解大型线性方程组
  • ✔ 极易并行化

2.3 误差与稳定性分析 ❶ 误差的来源与发展

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

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

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

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

\[ \epsilon = N - D \]

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

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

2.3 误差与稳定性分析 ❶ 误差的来源与发展

\[ \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.3 误差与稳定性分析 ❶ 误差的来源与发展

显式格式：

\[ \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.3 误差与稳定性分析 ❶ 误差的来源与发展

\[ \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.3 误差与稳定性分析 ❷ 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_i^n=\epsilon(x_i,t_n) = \sum_m A_m(t_n) e^{i k_m x_i} \]

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

\[ \epsilon_i^n = \sum_m e^{a_m t_n} e^{i k_m x_i} \]

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

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

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

放大因子定义：

\[ 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 - \frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k \Delta x}{2} \right)\le 1\Longrightarrow\frac{4\alpha \Delta t}{(\Delta x)^2}\sin^2 \left( \frac{k \Delta x}{2} \right)\ge 0 (自动满足) \]

同时需要满足

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

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

\[ \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 \colorbox{yellow}{$\displaystyle \textcolor{black}{\frac{\alpha \Delta t}{(\Delta x)^2} \le \frac{1}{2}}$} \]

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

物理与数值意义

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

von Neumann 稳定性分析方法

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

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

von Neumann 稳定性分析方法

👉 方法梳理：

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

2.3 误差与稳定性分析 ❷ von Neumann 稳定性分析（Fourier模态分析）

📘 作业二：采用von Neumann方法分析Crank–Nicolson格式的稳定性

一维热传导方程

\[\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^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] \]

提交时间：5月9日（课代表收齐后统一提交）

2.3 误差与稳定性分析 ❸ 示例：线性对流方程

FTCS格式

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

我们采用时间前向差分、空间中心差分 👉 FTCS （Forward Time, Central Space）格式

\[ \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 稳定性分析

设误差：

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

2.3 误差与稳定性分析 ❸ 示例：线性对流方程

FTCS格式

引入 Fourier 误差模态

\[ \epsilon_i^n = e^{a t_n} e^{ikx_i} \]

其中： \(x_i = i\Delta x\)， \(k\) 表示 \(x\) 方向波数

时间推进项

\[ \epsilon_i^{n+1} = e^{a t_{n+1}} e^{ikx_i} \]

空间项

\[ \epsilon_{i+1}^n = e^{a t_n} e^{ik(x_i+\Delta x)} = \epsilon_i^n e^{ik\Delta x} \]

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

2.3 误差与稳定性分析 ❸ 示例：线性对流方程

FTCS格式

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

左边：

\[ \frac{\epsilon_i^{n+1} - \epsilon_i^n}{\Delta t} =\frac{\epsilon_i^n}{\Delta t}\left(e^{a\Delta t} - 1\right) \]

右边：

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

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

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

2.3 误差与稳定性分析 ❸ 示例：线性对流方程

FTCS格式

放大因子

\[ G = \left|e^{a\Delta t}\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 + \left(\frac{c\Delta t}{\Delta x}\right)^2 \sin^2(k\Delta x)>1 （无条件不稳定） \]

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

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

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

当 \(c>0\) 时，采用 一阶迎风格式（Upwind） 👉 FTBS（Forward Time, Backward Space）

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

对该格式进行 von Neumann 稳定性分析

设误差：

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

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

引入 Fourier 误差模态

\[ \epsilon_i^n = e^{a t_n} e^{ikx_i} \]

其中： \(x_i = i\Delta x\)

时间推进项：

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

空间项：

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

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

代入差分格式：

\[ \frac{\epsilon_i^{n+1} - \epsilon_i^n}{\Delta t} = -c \frac{\epsilon_i^n - \epsilon_{i-1}^n}{\Delta x} \]

代入Fourier分解：

\[ \frac{\epsilon_i^n}{\Delta t}(e^{a\Delta t} - 1)=-c \frac{\epsilon_i^n(1 - e^{-ik\Delta x})}{\Delta x} \]

约去 \(\epsilon_i^n\)：

\[ e^{a\Delta t} = 1 - \frac{c\Delta t}{\Delta x}(1 - e^{-ik\Delta x}) \]

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

定义 CFL 数：

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

放大因子：

\[ G = 1 - \nu(1 - e^{-ik\Delta x}) \]

展开指数：

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

得到：

\[ G = 1 - \nu(1 - \cos\theta) - i\nu\sin\theta,\qquad \theta = k\Delta x \]

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

模平方：

\[ |G|^2 = [1 - \nu(1-\cos\theta)]^2 + \nu^2 \sin^2\theta \]

整理：

\[ |G|^2 = 1 - 2\nu(1-\cos\theta) + \nu^2[(1-\cos\theta)^2 + \sin^2\theta] \]

利用恒等式：

\[ (1-\cos\theta)^2 + \sin^2\theta = 2(1-\cos\theta) \]

得到：

\[ |G|^2 = 1 - 2\nu(1-\nu)(1-\cos\theta) \]

2.3 误差与稳定性分析 ❹ 示例：线性对流方程

FTBS格式

由于：

\[ 1-\cos\theta \ge 0 \]

要使

\[ |G| \le 1 \]

必须满足

\[ 0 \le \nu \le 1 \]

👉 结论：

• 当 \(0 \le \nu \le 1\) 时，FTBS格式稳定
• 当 \(\nu > 1\) 时，FTBS格式不稳定

2.3 误差与稳定性分析 ❺ 示例：线性对流方程

Lax 格式

仍考虑线性对流方程：

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

Lax 格式将FTCS时间推进中的 \(u_i^n\) 用相邻点平均替代：

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

也可写成

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

其中

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

2.3 误差与稳定性分析 ❺ 示例：线性对流方程

Lax 格式

设 Fourier 误差模态为

\[ \epsilon_i^n = e^{a t_n} e^{ikx_i} \]

则：

\[ \epsilon_{i+1}^n = \epsilon_i^n e^{ik\Delta x},\qquad \epsilon_{i-1}^n = \epsilon_i^n e^{-ik\Delta x} \]

以及

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

代入 Lax 格式：

\[ \epsilon_i^{n+1} = \frac{1}{2}(\epsilon_{i+1}^n + \epsilon_{i-1}^n) - \frac{\nu}{2}(\epsilon_{i+1}^n - \epsilon_{i-1}^n) \]

2.3 误差与稳定性分析 ❺ 示例：线性对流方程

Lax 格式

约去公共因子 \(\epsilon_i^n\)：

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

利用恒等式：

\[ \frac{1}{2}(e^{i\theta}+e^{-i\theta}) = \cos\theta \]

\[ \frac{1}{2}(e^{i\theta}-e^{-i\theta}) = i\sin\theta \]

得到放大因子：

\[ G = \cos\theta - i\nu\sin\theta, \qquad \theta = k\Delta x \]

2.3 误差与稳定性分析 ❺ 示例：线性对流方程

Lax 格式

模平方：

\[ |G|^2 = \cos^2\theta + \nu^2\sin^2\theta \]

要使格式稳定，需满足

\[ |G|^2 \le 1 \]

即

\[ \cos^2\theta + \nu^2\sin^2\theta \le 1 \]

整理可得：

\[ \nu^2 \le 1 \]

2.3 误差与稳定性分析 ❺ 示例：线性对流方程

Lax 格式

👉 结论：

Lax 格式的稳定条件为

\[ |\nu| \le 1 \]

即：

\[ \left|\frac{c\Delta t}{\Delta x}\right| \le 1 \]

👉 物理理解：

• Lax 格式通过“邻点平均”引入了人工耗散
• 因而比 FTCS 稳定得多（但会使解出现明显平滑）

2.3 误差与稳定性分析 ❻ 示例：线性对流方程

BTCS 格式（Backward Time, Central Space）

同样考虑线性对流方程：

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

采用 BTCS 格式 👉（后向时间 + 中心空间）

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

设误差模态：

\[ \epsilon_i^n = e^{a t_n} e^{ikx_i} \]

则：

\[ \epsilon_i^{n+1} = \epsilon_i^n e^{a\Delta t}, \qquad \epsilon_{i\pm1}^{n+1}= \epsilon_i^{n+1} e^{\pm ik\Delta x} \]

2.3 误差与稳定性分析 ❻ 示例：线性对流方程

BTCS 格式

代入格式：

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

得到

\[ \frac{\epsilon_i^n(e^{a\Delta t}-1)}{\Delta t}=-c \frac{\epsilon_i^{n+1}(e^{ik\Delta x}-e^{-ik\Delta x})}{2\Delta x} \]

约去 \(\epsilon_i^n\) 后有：

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

2.3 误差与稳定性分析 ❻ 示例：线性对流方程

BTCS 格式

利用：

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

得到：

\[ e^{a\Delta t} = \frac{1}{1 + i\nu\sin(k\Delta x)} \]

因此放大因子为：

\[ G = \frac{1}{1 + i\nu\sin\theta}, \qquad \theta = k\Delta x, \qquad |G| = \frac{1}{\sqrt{1+\nu^2\sin^2\theta}} \le 1 \]

对任意 \(\nu\) 都成立

2.3 误差与稳定性分析 ❻ 示例：线性对流方程

BTCS 格式

👉 结论：

BTCS 对线性对流方程无条件稳定

👉 说明：

虽然稳定，但该格式对纯对流问题会引入较强数值耗散，甚至相位误差也较明显。

2.3 误差与稳定性分析 ❼ 示例：线性对流方程

格式稳定性对比

其中：

\[ \theta = k\Delta x, \qquad \nu = c \Delta t /\Delta x \]

👉 对于双曲型对流问题，稳定往往需要以“耗散”为代价；

算例实践：Couette 流动

参考书

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

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

\[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]

这里：

• \(c\) 为波速
• 方程具有两条特征线

\[ x - ct = \text{const}, \qquad x + ct = \text{const} \]

也即

\[ \frac{dx}{dt} = c, \qquad \frac{dx}{dt} = -c \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

采用二阶中心差分离散：

• 时间二阶中心差分
• 空间二阶中心差分

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

整理得

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

定义 Courant 数

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

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

对该格式进行 von Neumann 稳定性分析

设误差模态为

\[ \epsilon_i^n = G^n e^{ikx_i} \]

其中：\(G\) 为放大因子、\(k\) 为波数、\(x_i = i\Delta x\)

\[ \epsilon_i^{n+1} = G^{n+1} e^{ikx_i},\qquad \epsilon_i^{n-1} = G^{n-1} e^{ikx_i} \]

以及

\[ \epsilon_{i+1}^n = G^n e^{ik(x_i+\Delta x)} = \epsilon_i^n e^{ik\Delta x} \]

\[ \epsilon_{i-1}^n = G^n e^{ik(x_i-\Delta x)} = \epsilon_i^n e^{-ik\Delta x} \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

将误差模态代入差分格式：

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

代入各项后得到：

\[ \frac{G^{n+1}e^{ikx_i} - 2G^n e^{ikx_i} + G^{n-1}e^{ikx_i}}{\Delta t^2}=c^2 \frac{G^n e^{ikx_i}(e^{ik\Delta x}-2+e^{-ik\Delta x})}{\Delta x^2} \]

约去公共因子 \(G^{\,n-1}e^{ikx_i}\)：

\[ \frac{G^2 - 2G + 1}{\Delta t^2}=c^2 \frac{G(e^{ik\Delta x}-2+e^{-ik\Delta x})}{\Delta x^2} \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

利用恒等式

\[ e^{ik\Delta x}-2+e^{-ik\Delta x}=2\cos\theta -2=-4\sin^2\frac{\theta}{2}, \qquad \theta = k\Delta x \]

代回可得

\[ G^2 - 2G + 1=-4C^2 G \sin^2\frac{\theta}{2} \]

整理为关于 \(G\) 的二次方程：

\[ G^2 + \left(4C^2\sin^2\frac{\theta}{2}-2\right)G + 1 = 0 \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

记

\[ \mu = 2C\sin\frac{\theta}{2} \]

则上式写成

\[ G^2 + (\mu^2 - 2)G + 1 = 0 \]

要使格式稳定，必须要求所有波数下误差不增长，即

\[ |G|\le 1 \]

对于这种二次方程，稳定要求根位于单位圆上或单位圆内。

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

令

\[ G = e^{i\varphi} \Rightarrow G + \frac{1}{G} = 2\cos\varphi \]

将方程两边除以 \(G\)：

\[ G + \frac{1}{G} + \mu^2 -2 = 0 \]

即

\[ 2\cos\varphi = 2 - \mu^2 \]

\[ -1 \le \cos\varphi \le 1 \Longrightarrow -2 \le 2-\mu^2 \le 2\Longrightarrow 0 \le \mu^2 \le 4 \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

因此

\[ \mu^2 = 4C^2 \sin^2\frac{\theta}{2} \le 4 \]

即

\[ C^2 \sin^2\frac{\theta}{2} \le 1 \]

由于对任意波数都有

\[ \sin^2\frac{\theta}{2} \le 1 \]

为了使所有模态都稳定，必须满足最严格条件

\[ C^2 \le 1 \Longrightarrow \colorbox{yellow}{$\displaystyle \textcolor{black}{C = \frac{c\Delta t}{\Delta x} \le 1}$} \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

👉 结论：

对于二阶中心差分离散的波动方程，

\[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]

其稳定条件（CFL 条件）为

\[ \colorbox{yellow}{$\displaystyle \textcolor{black}{C = \frac{c\Delta t}{\Delta x} \le 1}$} \]

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

物理解释

波动方程有两条物理特征线：

\[ x - ct = \text{const}, \qquad x + ct = \text{const} \]

也可写为

\[ \frac{dx}{dt} = c,\qquad \frac{dx}{dt} = -c \]

在一个时间步 \(\Delta t\) 内，物理信息传播距离为

\[ c\Delta t \]

而数值格式在空间上只能通过相邻网格点传递信息，其作用范围约为

\[ \Delta x \]

因此必须要求：\(c\Delta t \le \Delta x\Longrightarrow C = \frac{c\Delta t}{\Delta x} \le 1\)

2.3 误差与稳定性分析 ❽ 示例：一维波动方程

物理解释

2.3 误差与稳定性分析 ❾

👉 本节核心思想：

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

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

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

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

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

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

2.3 误差与稳定性分析 ❾

👉 检查：以下内容是否掌握：

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

2.4 网格与坐标变换 ❶ 基本思想

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

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

2.4 网格与坐标变换 ❶ 基本思想

网格生成（Grid Generation）

两个本质问题

① 网格生成 👉 几何问题

② 方程求解 👉 数值问题

⚠️ 二者本质不同！

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

但现实问题：

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

👉 无法直接使用矩形网格

2.4 网格与坐标变换 ❶ 基本思想

2.4 网格与坐标变换 ❷ 坐标变换

曲线网格（物理空间）

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

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

👉 解决方法：坐标变换

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

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

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

👉 可直接使用有限差分

2.4 网格与坐标变换 ❷ 坐标变换

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

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

数值求解流程

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

👉需要掌握

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

2.4 网格与坐标变换 ❷ 坐标变换

考虑二维非定常流：

\[ (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.4 网格与坐标变换 ❷ 坐标变换

链式法则（二阶导数）

\[ \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}\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.4 网格与坐标变换 ❷ 坐标变换

混合导数

\[ \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.4 网格与坐标变换 ❷ 坐标变换

评注

• 计算在规则网格进行
• 方程复杂性增加
• 几何复杂性转移到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.4 网格与坐标变换 ❷ 坐标变换

度量系数（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.4 网格与坐标变换 ❷ 坐标变换

2.4 网格与坐标变换 ❷ 坐标变换

考虑流动变量 \(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.4 网格与坐标变换 ❷ 坐标变换

用克拉默法则求解

\(\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.4 网格与坐标变换 ❷ 坐标变换

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

\[ \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.4 网格与坐标变换 ❷ 坐标变换

评注

• 直接度量（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.4 网格与坐标变换 ❷ 坐标变换

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.4 网格与坐标变换 ❷ 坐标变换

两边乘以\(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.4 网格与坐标变换 ❷ 坐标变换

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

\[ \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.4 网格与坐标变换 ❷ 坐标变换

最终守恒形式化简为：

\[ \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.4 网格与坐标变换 ❷ 坐标变换

小结

这一节的关键是：

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

本质上：

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

2.4 网格与坐标变换 ❸ 网格变换实例

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

2.4 网格与坐标变换 ❸ 网格变换实例

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

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

\[ \xi = x,\quad \eta = \ln(y+1) \Longleftrightarrow x=\xi,\quad 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) \]

将其代入网格变换式，得

\[ \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 贴体坐标系：椭圆型方程用于网格生成 ❶

由上面的讨论可知：

• 物理平面中曲线边界上的所有点，经过变换后，都落在计算平面中的同一条水平线 \(\eta=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 贴体坐标系：椭圆型方程用于网格生成 ❷

保角变换与椭圆型网格生成

对于二维问题，保角变换（conformal mapping）提供了一种能够同时满足正交性和光滑性要求的网格生成方法。

我们把物理平面和计算平面中的网格点坐标分别表示为复数形式：

\[ Z=x+iy \]

\[ \zeta=\xi+i\eta \]

其中：

• \(Z=x+iy\) 表示物理平面；
• \(\zeta=\xi+i\eta\) 表示计算平面；
• \(i\) 为虚数单位。

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

保角变换与椭圆型网格生成

我们把物理平面和计算平面中的网格点坐标分别表示为复数形式：

\[ Z=x+iy \]

\[ \zeta=\xi+i\eta \]

如果

\[ \zeta=\zeta(Z) \]

是解析函数（函数不仅在某一点可导，而且在该点附近一小片区域内处处可导），则从 \((x,y)\) 到 \((\xi,\eta)\) 的变换称为保角变换。

保角变换的优良性质主要体现在：保角变换能够保持局部夹角不变。

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

保角变换与椭圆型网格生成

也就是说，如果在计算平面 \((\xi,\eta)\) 中两组网格线是正交的，那么经过保角变换后，在物理平面 \((x,y)\) 中对应的两组网格线仍然是正交的。

因此，保角变换具有两个重要特点：

• 网格线之间的夹角保持不变；
• 变换是光滑的；
• 若计算平面中网格正交，则物理平面中的网格也正交。

不过，保角变换主要适用于二维问题。对于一般复杂区域，构造合适的保角变换通常比较困难。

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

Cauchy–Riemann 条件

设

\[ \zeta=\xi+i\eta \]

是关于

\[ Z=x+iy \]

的解析函数，即

\[ \zeta=\zeta(Z). \]

则实部 \(\xi(x,y)\) 和虚部 \(\eta(x,y)\) 必须满足 Cauchy–Riemann 条件：

\[ \frac{\partial \xi}{\partial x}= \frac{\partial \eta}{\partial y}, \qquad \frac{\partial \xi}{\partial y}=-\frac{\partial \eta}{\partial x}. \]

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

Cauchy–Riemann 条件

\[ \frac{\partial \xi}{\partial x}= \frac{\partial \eta}{\partial y}, \qquad \frac{\partial \xi}{\partial y}=-\frac{\partial \eta}{\partial x}. \]

也可简写为：

\[ \xi_x=\eta_y, \qquad \xi_y=-\eta_x. \]

这两个关系是复变函数解析性的基本条件。

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

由 Cauchy–Riemann 条件：

\[ \xi_x=\eta_y, \qquad \xi_y=-\eta_x \]

对第一个式子关于 \(x\) 求导：

\[ \xi_{xx}=\eta_{yx} \]

对第二个式子关于 \(y\) 求导：

\[ \xi_{yy}=-\eta_{xy} \]

两式相加，得到：

\[ \xi_{xx}+\xi_{yy}= \eta_{yx}-\eta_{xy}=0 \]

同理，

\[ \eta_{xx}+\eta_{yy}=0. \]

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

因此，由 Cauchy–Riemann 条件可以推出：

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

也就是：

\[ \nabla^2 \xi=0, \qquad \nabla^2 \eta=0. \]

这说明：

如果 \(\zeta=\zeta(Z)\) 是解析函数，则 \(\xi\) 和 \(\eta\) 都是调和函数（指满足 Laplace 方程的函数）。

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

Laplace 方程是椭圆型方程，具有明显的平滑作用。

如果在边界上给定网格点的位置，那么求解 Laplace 方程可以把边界信息平滑地传递到区域内部。

其作用可以理解为：边界决定网格形状，Laplace 方程负责把这种形状平滑地扩散到内部。

因此，采用 Laplace 方程生成网格具有以下特点：

• 内部网格较光滑；
• 网格线不容易发生剧烈振荡；
• 适合处理封闭区域或复杂边界；
• 可以自然形成边界贴合网格。

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

与 CFD 网格生成的联系

在 CFD 中，网格生成的核心目标是：

• 边界要贴合物体；
• 近壁区域要足够密；
• 网格线要尽量光滑；
• 网格不能交叉；
• Jacobian 不能为零；
• 尽量减少非正交性和网格畸变。

保角变换可以很好地满足正交性和光滑性，但工程适用范围有限。

Laplace 方程虽然不能保证严格正交，但作为椭圆型方程，它具有天然的平滑作用，因此成为贴体网格生成的重要基础。

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 贴体坐标系：椭圆型方程用于网格生成 ❷

知识补充：二阶度量系数

前面已经得到一阶度量系数：

\[ \xi_x,\quad \xi_y,\quad \eta_x,\quad \eta_y \]

有：

\[ \begin{bmatrix} \xi_x & \xi_y\\ \eta_x & \eta_y \end{bmatrix}= \frac{1}{J} \begin{bmatrix} y_\eta & -x_\eta\\ -y_\xi & x_\xi \end{bmatrix} \]

其中 Jacobian 为：

\[ J= \frac{\partial(x,y)}{\partial(\xi,\eta)}= x_\xi y_\eta-x_\eta y_\xi \]

即：

\[ \boxed{ \xi_x=\frac{y_\eta}{J}, \qquad \xi_y=-\frac{x_\eta}{J}, \qquad \eta_x=-\frac{y_\xi}{J}, \qquad \eta_y=\frac{x_\xi}{J} } \]

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

知识补充：二阶度量系数

从恒等关系出发

\[ \frac{\partial \xi}{\partial \xi}=1 \]

根据链式法则：

\[ \xi_\xi=\xi_x x_\xi+\xi_y y_\xi=1 \]

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

\[ \xi_\eta=\xi_x x_\eta+\xi_y y_\eta=0 \]

对 \(\xi\) 求导：

\[ \left( \xi_x x_\xi+\xi_y y_\xi \right)_\xi=0 \]

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

知识补充：二阶度量系数

对 \(\xi\) 求导：

\[ \left( \xi_x x_\xi+\xi_y y_\xi \right)_\xi=0 \]

使用乘积法则：

\[ (\xi_x)_\xi x_\xi+\xi_x x_{\xi\xi}+ (\xi_y)_\xi y_\xi+\xi_y y_{\xi\xi}=0 \]

其中：

\[ (\xi_x)_\xi= \xi_{xx}x_\xi+\xi_{xy}y_\xi \]

\[ (\xi_y)_\xi= \xi_{xy}x_\xi+\xi_{yy}y_\xi \]

代入得：

\[ (\xi_{xx}x_\xi+\xi_{xy}y_\xi)x_\xi+ (\xi_{xy}x_\xi+\xi_{yy}y_\xi)y_\xi+ \xi_x x_{\xi\xi}+ \xi_y y_{\xi\xi}=0 \]

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

知识补充：二阶度量系数

\[ (\xi_{xx}x_\xi+\xi_{xy}y_\xi)x_\xi+ (\xi_{xy}x_\xi+\xi_{yy}y_\xi)y_\xi+ \xi_x x_{\xi\xi}+ \xi_y y_{\xi\xi}=0 \]

整理为：

\[ \boxed{ \xi_{xx}(x_\xi)^2+ 2\xi_{xy}x_\xi y_\xi+ \xi_{yy}(y_\xi)^2+ \xi_x x_{\xi\xi}+ \xi_y y_{\xi\xi}=0 } \]

仍从

\[ \xi_x x_\xi+\xi_y y_\xi=1 \]

出发，对 \(\eta\) 求导：

\[ \left( \xi_x x_\xi+\xi_y y_\xi \right)_\eta=0 \]

展开：

\[ (\xi_x)_\eta x_\xi+\xi_x x_{\xi\eta} + (\xi_y)_\eta y_\xi+\xi_y y_{\xi\eta}=0 \]

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

知识补充：二阶度量系数

\[ (\xi_x)_\eta x_\xi+\xi_x x_{\xi\eta} + (\xi_y)_\eta y_\xi+\xi_y y_{\xi\eta}=0 \]

其中：

\[ (\xi_x)_\eta= \xi_{xx}x_\eta+\xi_{xy}y_\eta \]

\[ (\xi_y)_\eta= \xi_{xy}x_\eta+\xi_{yy}y_\eta \]

代入得：

\[ (\xi_{xx}x_\eta+\xi_{xy}y_\eta)x_\xi+ (\xi_{xy}x_\eta+\xi_{yy}y_\eta)y_\xi+ \xi_x x_{\xi\eta}+ \xi_y y_{\xi\eta}=0 \]

整理为：

\[ \boxed{ \xi_{xx}x_\xi x_\eta+ \xi_{xy}(x_\xi y_\eta+x_\eta y_\xi)+ \xi_{yy}y_\xi y_\eta+ \xi_x x_{\xi\eta}+ \xi_y y_{\xi\eta}=0 } \]

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

知识补充：二阶度量系数

对第二个恒等式关于 \(\eta\) 求导

\[ \xi_x x_\eta+\xi_y y_\eta=0 \]

\[ \left( \xi_x x_\eta+\xi_y y_\eta \right)_\eta=0 \]

展开：

\[ (\xi_x)_\eta x_\eta+\xi_x x_{\eta\eta}+ (\xi_y)_\eta y_\eta+\xi_y y_{\eta\eta}=0 \]

代入

\[ (\xi_x)_\eta= \xi_{xx}x_\eta+\xi_{xy}y_\eta \]

\[ (\xi_y)_\eta= \xi_{xy}x_\eta+\xi_{yy}y_\eta \]

得到：

\[ (\xi_{xx}x_\eta+\xi_{xy}y_\eta)x_\eta+ (\xi_{xy}x_\eta+\xi_{yy}y_\eta)y_\eta+ \xi_x x_{\eta\eta}+ \xi_y y_{\eta\eta}=0 \]

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

知识补充：二阶度量系数

\[ (\xi_{xx}x_\eta+\xi_{xy}y_\eta)x_\eta+ (\xi_{xy}x_\eta+\xi_{yy}y_\eta)y_\eta+ \xi_x x_{\eta\eta}+ \xi_y y_{\eta\eta}=0 \]

整理为：

\[ \boxed{ \xi_{xx}(x_\eta)^2+ 2\xi_{xy}x_\eta y_\eta+ \xi_{yy}(y_\eta)^2+ \xi_x x_{\eta\eta}+ \xi_y y_{\eta\eta}=0 } \]

结合前面得到的

\[ \boxed{ \xi_{xx}x_\xi x_\eta+ \xi_{xy}(x_\xi y_\eta+x_\eta y_\xi)+ \xi_{yy}y_\xi y_\eta+ \xi_x x_{\xi\eta}+ \xi_y y_{\xi\eta}=0 } \]

\[ \boxed{ \xi_{xx}(x_\xi)^2+ 2\xi_{xy}x_\xi y_\xi+ \xi_{yy}(y_\xi)^2+ \xi_x x_{\xi\xi}+ \xi_y y_{\xi\xi}=0 } \]

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

知识补充：二阶度量系数

整理为矩阵形式：

\[ \begin{bmatrix} (x_\xi)^2 & 2x_\xi y_\xi & (y_\xi)^2 \\[4pt] x_\xi x_\eta & x_\xi y_\eta+x_\eta y_\xi & y_\xi y_\eta \\[4pt] (x_\eta)^2 & 2x_\eta y_\eta & (y_\eta)^2 \end{bmatrix} \begin{bmatrix} \xi_{xx}\\[4pt] \xi_{xy}\\[4pt] \xi_{yy} \end{bmatrix}=- \begin{bmatrix} \xi_x x_{\xi\xi}+\xi_y y_{\xi\xi}\\[4pt] \xi_x x_{\xi\eta}+\xi_y y_{\xi\eta}\\[4pt] \xi_x x_{\eta\eta}+\xi_y y_{\eta\eta} \end{bmatrix} \]

因此：

\[ \boxed{ \begin{bmatrix} \xi_{xx}\\[4pt] \xi_{xy}\\[4pt] \xi_{yy} \end{bmatrix}=- \begin{bmatrix} (x_\xi)^2 & 2x_\xi y_\xi & (y_\xi)^2 \\[4pt] x_\xi x_\eta & x_\xi y_\eta+x_\eta y_\xi & y_\xi y_\eta \\[4pt] (x_\eta)^2 & 2x_\eta y_\eta & (y_\eta)^2 \end{bmatrix}^{-1} \begin{bmatrix} \xi_x x_{\xi\xi}+\xi_y y_{\xi\xi}\\[4pt] \xi_x x_{\xi\eta}+\xi_y y_{\xi\eta}\\[4pt] \xi_x x_{\eta\eta}+\xi_y y_{\eta\eta} \end{bmatrix} } \]

这就给出了 \(\xi\) 的二阶度量系数。

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

知识补充：二阶度量系数

同理，

\[ \begin{bmatrix} (x_\xi)^2 & 2x_\xi y_\xi & (y_\xi)^2 \\[4pt] x_\xi x_\eta & x_\xi y_\eta+x_\eta y_\xi & y_\xi y_\eta \\[4pt] (x_\eta)^2 & 2x_\eta y_\eta & (y_\eta)^2 \end{bmatrix} \begin{bmatrix} \eta_{xx}\\[4pt] \eta_{xy}\\[4pt] \eta_{yy} \end{bmatrix}=- \begin{bmatrix} \eta_x x_{\xi\xi}+\eta_y y_{\xi\xi}\\[4pt] \eta_x x_{\xi\eta}+\eta_y y_{\xi\eta}\\[4pt] \eta_x x_{\eta\eta}+\eta_y y_{\eta\eta} \end{bmatrix} \]

因此：

\[ \boxed{ \begin{bmatrix} \eta_{xx}\\[4pt] \eta_{xy}\\[4pt] \eta_{yy} \end{bmatrix}=- \begin{bmatrix} (x_\xi)^2 & 2x_\xi y_\xi & (y_\xi)^2 \\[4pt] x_\xi x_\eta & x_\xi y_\eta+x_\eta y_\xi & y_\xi y_\eta \\[4pt] (x_\eta)^2 & 2x_\eta y_\eta & (y_\eta)^2 \end{bmatrix}^{-1} \begin{bmatrix} \eta_x x_{\xi\xi}+\eta_y y_{\xi\xi}\\[4pt] \eta_x x_{\xi\eta}+\eta_y y_{\xi\eta}\\[4pt] \eta_x x_{\eta\eta}+\eta_y y_{\eta\eta} \end{bmatrix} } \]

给出了 \(\eta\) 的二阶度量系数。

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

知识补充：二阶度量系数

求矩阵逆得到

\[ \begin{bmatrix} \xi_{xx}\\[4pt] \xi_{xy}\\[4pt] \xi_{yy} \end{bmatrix}= \frac{-1} { J\left[(x_\xi y_\eta)^2+(x_\eta y_\xi)^2\right] } \begin{bmatrix} J(y_\eta)^2 & -2J y_\xi y_\eta & J(y_\xi)^2 \\[4pt] -Jx_\eta y_\eta & J(x_\xi y_\eta+x_\eta y_\xi) & -Jx_\xi y_\xi \\[4pt] J(x_\eta)^2 & -2Jx_\xi x_\eta & J(x_\xi)^2 \end{bmatrix} \begin{bmatrix} \xi_x x_{\xi\xi}+\xi_y y_{\xi\xi}\\[4pt] \xi_x x_{\xi\eta}+\xi_y y_{\xi\eta}\\[4pt] \xi_x x_{\eta\eta}+\xi_y y_{\eta\eta} \end{bmatrix}. \]

\[ \begin{bmatrix} \eta_{xx}\\[4pt] \eta_{xy}\\[4pt] \eta_{yy} \end{bmatrix}= \frac{-1} { J\left[(x_\xi y_\eta)^2+(x_\eta y_\xi)^2\right] } \begin{bmatrix} J(y_\eta)^2 & -2J y_\xi y_\eta & J(y_\xi)^2 \\[4pt] -Jx_\eta y_\eta & J(x_\xi y_\eta+x_\eta y_\xi) & -Jx_\xi y_\xi \\[4pt] J(x_\eta)^2 & -2Jx_\xi x_\eta & J(x_\xi)^2 \end{bmatrix} \begin{bmatrix} \eta_x x_{\xi\xi}+\eta_y y_{\xi\xi}\\[4pt] \eta_x x_{\xi\eta}+\eta_y y_{\xi\eta}\\[4pt] \eta_x x_{\eta\eta}+\eta_y y_{\eta\eta} \end{bmatrix}. \]

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

回到椭圆型偏微分方程定义坐标变换：

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

代入二阶度量系数表达式，得到：

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

系数定义

\[ \alpha = \left(\frac{\partial x}{\partial \eta}\right)^2 + \left(\frac{\partial y}{\partial \eta}\right)^2,\qquad \beta = \frac{\partial x}{\partial \xi}\frac{\partial x}{\partial \eta}+ \frac{\partial y}{\partial \xi}\frac{\partial y}{\partial \eta},\qquad \gamma = \left(\frac{\partial x}{\partial \xi}\right)^2 + \left(\frac{\partial y}{\partial \xi}\right)^2 \]

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 网格生成的现代发展与网格生成软件介绍 ❹

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

Star-CCM+（Simcenter）

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

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

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

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

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

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

Fluent Meshing

• ANSYS Fluent 内置网格模块

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

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

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

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

Gmsh（开源工具）

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

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

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

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

OpenFOAM 自带工具

blockMesh

• 结构网格
• 完全可控

snappyHexMesh

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

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

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

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

GridPro

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

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

NUMECA (Autogrid)

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

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

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

一、作业内容

本次大作业要求同学们使用自编代码或开源 CFD 程序，完成一个二维或三维流动问题的数值模拟。流动可以是定常，也可以是非定常。

题目自选，例如：

• 二维圆形、二维方形、三维圆柱、三维球体、二维翼型或三维机翼绕流；
• 二维/三维方腔驱动流动；
• 激波管、激波反射流动；
• 喷管流动；
• 边界层流动；
• 后台阶流动；
• 射流、尾迹或剪切层流动；
• 其他 CFD 算例。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

二、基本要求

每个项目需至少体现 3 项 CFD 关键技术，包括但不限于：

1. 网格技术
   • 结构网格、非结构网格；
   • 贴体网格；
   • 网格拉伸；
   • 坐标变换；
   • 多块网格；
   • 动网格或自适应网格。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

二、基本要求

每个项目需至少体现 3 项 CFD 关键技术，包括但不限于：

2. 空间离散格式
   • 中心格式、迎风格式；
   • 高阶格式；
   • MUSCL、ENO、WENO；
   • 紧致差分格式。
3. 时间推进方法
   • 显式 Euler 方法；
   • Runge-Kutta 方法；
   • 隐式方法；
   • LU-SGS；
   • 双时间步方法。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

二、基本要求

每个项目需至少体现 3 项 CFD 关键技术，包括但不限于：

4. 激波与间断处理
   • Roe 格式；
   • HLL/HLLC 格式；
   • AUSM 格式；
   • TVD 限制器；
   • WENO 激波捕捉；
   • 人工黏性方法；
   • 激波装配方法。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

二、基本要求

每个项目需至少体现 3 项 CFD 关键技术，包括但不限于：

5. 物理模型与边界条件
   • 不可压/可压流动；
   • 层流、RANS、LES；
   • 无滑移壁面；
   • 入口/出口边界；
   • 远场边界；
   • 周期边界；
   • 对称边界。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

二、基本要求

每个项目需至少体现 3 项 CFD 关键技术，包括但不限于：

6. 计算实现与后处理
   • 收敛加速；
   • 并行计算；
   • 网格无关性分析；
   • 误差分析；
   • 流场可视化；
   • 与理论解、实验结果或文献结果对比。

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

三、完成方式

• 可个人完成，也可组队完成（但最多2人）；
• 可自编代码（鼓励借助AI）
• 也可使用 OpenFOAM、SU2、Nek5000 等开源代码；

四、提交与汇报要求

• 最后一节课（6月3日）进行课堂 PPT 汇报和代码演示，汇报时间3分钟（个人）或6分钟（组队）。
• PPT 内容建议包括：1. 计算问题描述；2. 网格与边界条件；3. 采用的 CFD 关键技术；4. 计算结果，须包含现场演示；
• 最终提交材料包括： PPT； 代码文件；

CFD 大作业：二维或三维流动数值模拟（6月3日随堂汇报）

五、评分方式（占总成绩15%）

评分重点：

• 是否实现二维或三维流动数值模拟；
• 是否体现至少 3 项 CFD 关键技术；
• 数值方法选择是否合理；
• 结果是否可信；
• 图像展示是否清晰；
• 汇报逻辑是否完整；
• 是否体现独立思考能力。

建议选择难度适中的算例，不要求追求复杂几何或大规模计算。