2025–2026 学年第二学期
北京理工大学 空天科学与技术学院
任杰
John D. Anderson Computational Fluid
Dynamics: The Basics with Applications
McGraw-Hill, 1995
任玉新 陈海昕
计算流体力学基础 清华大学出版社，2006
张德良
计算流体力学教程 高等教育出版社，2010
蔡庆东
计算流体动力学
北京大学本科生研究生专业基础课讲义，2018
真正理解 CFD 的最好方式，是自己实现数值算法并运行计算 - John D. Anderson
前面章节已经介绍了 CFD 的基本思想和若干典型算法。
这些方法的特点是：
算法形式相对简单，便于入门理解；
同时又具有足够实用性，可以求解一批典型流动问题。
但是，现代 CFD 中大量算法更加复杂，也更高效。
它们往往来自：
CONSERVATION FORM REVISITED
重新认识守恒形式
守恒形式不仅是“写法问题”，
它直接关系到激波捕捉、方程类型、特征传播和现代数值格式构造。
现代 CFD 算法大量建立在控制方程的数学结构之上。
尤其在高速流和激波问题中，通常采用：
原因包括：
非守恒形式在光滑区域可以等价，但在含激波或间断的问题中，往往会导致错误的跳跃关系。
三维控制方程可写成一般形式：
\[ \frac{\partial U}{\partial t} +\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} +\frac{\partial H}{\partial z} =J \]
其中：
对于 Euler 或 Navier–Stokes 方程，通量通常不是 \(U\) 的线性函数，这里写做：
\[ F=F(U),\qquad G=G(U),\qquad H=H(U) \]
由于：
\[ F=F(U),\qquad G=G(U),\qquad H=H(U) \]
根据链式法则：
\[ \frac{\partial F}{\partial x}= \frac{\partial F}{\partial U} \frac{\partial U}{\partial x} \]
类似地：
\[ \frac{\partial G}{\partial y}= \frac{\partial G}{\partial U} \frac{\partial U}{\partial y}, \qquad \frac{\partial H}{\partial z}= \frac{\partial H}{\partial U} \frac{\partial U}{\partial z} \]
于是控制方程可写为：
\[ \frac{\partial U}{\partial t} +\frac{\partial F}{\partial U}\frac{\partial U}{\partial x} +\frac{\partial G}{\partial U}\frac{\partial U}{\partial y} +\frac{\partial H}{\partial U}\frac{\partial U}{\partial z} =J \]
定义：
\[ A\equiv\frac{\partial F}{\partial U}, \qquad B\equiv\frac{\partial G}{\partial U}, \qquad C\equiv\frac{\partial H}{\partial U} \]
则控制方程写成：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial U}{\partial t} +A\frac{\partial U}{\partial x} +B\frac{\partial U}{\partial y} +C\frac{\partial U}{\partial z} =J } $} \]
这里的 \(A,B,C\) 称为通量向量的 Jacobian 矩阵。
注意：这里的 Jacobian 不是坐标变换中的 Jacobian。
这里表示通量向量对守恒变量向量的导数。
Jacobian 矩阵的重要性体现在两个方面：
Jacobian 矩阵的特征值决定方程组的数学性质。
Jacobian 矩阵的特征值给出信息在流场中的传播速度和传播方向。
为简化讨论，考虑一维非定常无粘流：
\[ \frac{\partial \rho}{\partial t} +\frac{\partial(\rho u)}{\partial x} =0 \]
\[ \frac{\partial(\rho u)}{\partial t} +\frac{\partial(\rho u^2+p)}{\partial x} =0 \]
\[ \frac{\partial(\rho E)}{\partial t} +\frac{\partial(\rho uE+pu)}{\partial x} =0 \]
写成向量形式：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial U}{\partial t} +\frac{\partial F}{\partial x} =0 } $} \]
守恒变量为：
\[ U= \begin{Bmatrix} \rho\\ \rho u\\ \rho E \end{Bmatrix} \triangleq\begin{Bmatrix} \rho\\ m\\ \varepsilon \end{Bmatrix} \]
通量向量为：
\[ F= \begin{Bmatrix} \rho u\\ \rho u^2+p\\ \rho uE+pu \end{Bmatrix}= \begin{Bmatrix} m\\ \dfrac{m^2}{\rho}+p\\ \dfrac{m}{\rho}(\varepsilon+p) \end{Bmatrix} \]
其中：
\[ E=e+\frac{u^2}{2} \]
\(E\) 为单位质量总能量，\(e\) 为单位质量内能。
理想气体状态方程为：
\[ p=(\gamma-1)\rho e \]
又因为：
\[ \varepsilon=\rho E =\rho\left(e+\frac{u^2}{2}\right) =\rho e+\frac{\rho u^2}{2} \]
由 \(m=\rho u\) 得：
\[ \rho e=\varepsilon-\frac{m^2}{2\rho} \]
因此：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ p=(\gamma-1)\left(\varepsilon-\frac{m^2}{2\rho}\right) } $} \]
将压力表达式代入通量向量：
\[ F= \begin{Bmatrix} \rho u\\ \rho u^2+p\\ \rho uE+pu \end{Bmatrix}= \begin{Bmatrix} m\\ \dfrac{m^2}{\rho}+p\\ \dfrac{m}{\rho}(\varepsilon+p) \end{Bmatrix}= \begin{Bmatrix} m\\[6pt] \dfrac{m^2}{\rho} +(\gamma-1)\left(\varepsilon-\dfrac{m^2}{2\rho}\right)\\[8pt] \dfrac{m}{\rho} \left[ \varepsilon +(\gamma-1)\left(\varepsilon-\dfrac{m^2}{2\rho}\right) \right] \end{Bmatrix} \]
此时，通量已经完全写成 \(U=(\rho,m,\varepsilon)^T\) 的函数：
\[ F=F(U) \]
因此可以计算一维通量 Jacobian：
\[ A=\frac{\partial F}{\partial U} \]
因为：
\[ U= \begin{Bmatrix} \rho\\ m\\ \varepsilon \end{Bmatrix}, \qquad \frac{\partial U}{\partial t}= \begin{Bmatrix} \dfrac{\partial \rho}{\partial t}\\[6pt] \dfrac{\partial m}{\partial t}\\[6pt] \dfrac{\partial \varepsilon}{\partial t} \end{Bmatrix}, \qquad \frac{\partial U}{\partial x}= \begin{Bmatrix} \dfrac{\partial \rho}{\partial x}\\[6pt] \dfrac{\partial m}{\partial x}\\[6pt] \dfrac{\partial \varepsilon}{\partial x} \end{Bmatrix} \]
Jacobian 矩阵定义为：
\[ A=\frac{\partial F}{\partial U}= \begin{bmatrix} \dfrac{\partial F_1}{\partial \rho} & \dfrac{\partial F_1}{\partial m} & \dfrac{\partial F_1}{\partial \varepsilon}\\[8pt] \dfrac{\partial F_2}{\partial \rho} & \dfrac{\partial F_2}{\partial m} & \dfrac{\partial F_2}{\partial \varepsilon}\\[8pt] \dfrac{\partial F_3}{\partial \rho} & \dfrac{\partial F_3}{\partial m} & \dfrac{\partial F_3}{\partial \varepsilon} \end{bmatrix} \]
逐项求导后，可得：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ A=\frac{\partial F}{\partial U}= \begin{bmatrix} \dfrac{\partial F_1}{\partial \rho} & \dfrac{\partial F_1}{\partial m} & \dfrac{\partial F_1}{\partial \varepsilon}\\[8pt] \dfrac{\partial F_2}{\partial \rho} & \dfrac{\partial F_2}{\partial m} & \dfrac{\partial F_2}{\partial \varepsilon}\\[8pt] \dfrac{\partial F_3}{\partial \rho} & \dfrac{\partial F_3}{\partial m} & \dfrac{\partial F_3}{\partial \varepsilon} \end{bmatrix}= \begin{bmatrix} 0 & 1 & 0\\[6pt] (\gamma-3)\dfrac{u^2}{2} & (3-\gamma)u & \gamma-1\\[8pt] (\gamma-1)u^3-\gamma uE & \gamma E-\dfrac{3}{2}(\gamma-1)u^2 & \gamma u \end{bmatrix} } $} \]
其中：
\[ u=\frac{m}{\rho}, \qquad E=\frac{\varepsilon}{\rho} \]
这个矩阵并不是为了“多写一个复杂形式”，而是为了揭示一维 Euler 方程的信息传播结构。
一维 Euler 方程可以写成：
\[ \frac{\partial U}{\partial t} +A\frac{\partial U}{\partial x} =0 \]
展开为：
\[ \frac{\partial}{\partial t} \begin{Bmatrix} \rho\\ \rho u\\ \rho E \end{Bmatrix} + A \frac{\partial}{\partial x} \begin{Bmatrix} \rho\\ \rho u\\ \rho E \end{Bmatrix} =0 \]
虽然形式发生了变化，但并没有丢失原方程信息。
Jacobian 拟线性形式与原守恒形式是等价的。
因此：
特征值由下式确定：
\[ |A-\lambda I|=0 \]
对于一维 Euler 方程，特征值为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \lambda_1=u,\qquad \lambda_2=u+a,\qquad \lambda_3=u-a } $} \]
其中 \(a\) 为声速。
由于三个特征值均为实数，因此一维无粘 Euler 方程是：
三个特征值代表三类信息传播速度：
随流体质点传播
\[ \frac{dx}{dt}=u \]
对应流体微团的运动轨迹。
相对流体向右传播的声波
\[ \frac{dx}{dt}=u+a \]
对应右行 Mach 波。
相对流体向左传播的声波
\[ \frac{dx}{dt}=u-a \]
对应左行 Mach 波。
在 \(x-t\) 平面中，特征线斜率为：
\[ \frac{dt}{dx}=\frac{1}{\lambda} \]
因此三族特征线为：
\[ \frac{dt}{dx}=\frac{1}{u}, \qquad \frac{dt}{dx}=\frac{1}{u+a}, \qquad \frac{dt}{dx}=\frac{1}{u-a} \]
Jacobian 矩阵的特征值给出了信息在流场中的传播方向。
这正是现代 CFD 中构造迎风格式、特征分解和通量分裂的基础。
对一维 Euler 方程，通量 Jacobian 的三个特征值为：
\[ \lambda_1=u,\qquad \lambda_2=u+a,\qquad \lambda_3=u-a \]
它们主要说明什么？
A. 数值格式的截断误差大小
B. 信息在流场中的传播速度和方向
C. 网格生成中的坐标变换关系
D. 黏性项的离散精度
参考判断：B
Jacobian 的特征值决定特征传播速度，是理解双曲型流动方程和现代迎风格式的基础。
激波管初始时由隔膜分成左右两部分：
一个初始间断，不会一直保持一个间断，而会分裂成几种不同性质的波。
初始条件写成：
\[ (\rho,u,p)(x,0)= \begin{cases} (\rho_L,u_L,p_L),&x<x_0,\\ (\rho_R,u_R,p_R),&x>x_0. \end{cases} \]
典型 Sod 激波管：
\[ (\rho_L,u_L,p_L)=(1,0,1),\qquad (\rho_R,u_R,p_R)=(0.125,0,0.1) \]
💻 代码演示
shock_tube_fd_exact_demo.m
隔膜移除后，左侧高压气体向右运动，右侧低压气体被推动和压缩。
这一过程不是简单地“整体流过去”，而是会产生三种典型结构：
Sod 激波管的典型波系：左膨胀波、接触间断、右激波。
可以先用一句话理解：
左侧高压气体释放能量，因此产生膨胀波。
右侧低压气体被压缩，因此产生激波。
左右两团气体虽然被带动向同一方向运动，但密度一般不同，因此中间留下接触间断。
隔膜移除后，流场不再只有左、右两个初始状态。
三种波会把流场分成四个主要区域：
先把区域分清楚，再分别理解连接这些区域的波。
四个区域之间的连接关系：
\[ L \quad \longrightarrow \quad L^* \quad \longrightarrow \quad R^* \quad \longrightarrow \quad R \]
对于 Sod 激波管：
\[ L \xrightarrow{\text{膨胀波}} L^* \xrightarrow{\text{接触间断}} R^* \xrightarrow{\text{激波}} R \]
\(L^*\) 和 \(R^*\) 是隔膜撤去后形成的中间状态。
它们由接触间断分开。
接触间断两侧的速度和压力相同，但密度一般不同。
接触间断两侧满足：
\[ p_L^*=p_R^*=p^* \]
\[ u_L^*=u_R^*=u^* \]
但是：
\[ \rho_L^*\ne\rho_R^* \]
精确解的核心量是中间压力 \(p^*\) 和中间速度 \(u^*\)。本节课只需要知道它们代表中间区域的共同压力和速度。
守恒变量：
\[ \mathbf U= \begin{bmatrix} \rho\\ \rho u\\ \rho E \end{bmatrix} \]
通量：
\[ \mathbf F= \begin{bmatrix} \rho u\\ \rho u^2+p\\ u(\rho E+p) \end{bmatrix} \]
一维 Euler 方程守恒形式为：
\[ \frac{\partial \mathbf U}{\partial t} +\frac{\partial \mathbf F(\mathbf U)}{\partial x}=0 \]
其中：
\[ E=e+\frac{u^2}{2},\qquad p=(\gamma-1)\rho e \]
激波管问题必须重视守恒形式，因为激波处变量会发生跳跃，普通光滑函数意义下的微分方程不再足够。
一维 Euler 方程有三条特征速度：
\[ \lambda_1=u-a,\qquad \lambda_2=u,\qquad \lambda_3=u+a \]
其中：
\[ a=\sqrt{\gamma p/\rho} \]
为当地声速。
三条特征速度的物理意义：
| 特征速度 | 物理含义 | 对应波族 |
|---|---|---|
| (u-a) | 相对流体向左传播的声学信息 | 1-波 |
| (u) | 随流体一起运动的物质信息 | 2-波 |
| (u+a) | 相对流体向右传播的声学信息 | 3-波 |
三种波与三族特征不是简单的“一条线对应一种波”，而是：
第一族和第三族是声波族，可以形成膨胀波，也可以形成激波。
第二族是物质波族，对应接触间断。
对 Sod 激波管，通常对应关系为：
| 连接区域 | 波结构 | 对应特征族 | 典型表现 |
|---|---|---|---|
| \(L\to L^*\) | 左膨胀波 | 1-族，\(u-a\) | 连续展开 |
| \(L^*\to R^*\) | 接触间断 | 2-族，\(u\) | 密度跳跃 |
| \(R^*\to R\) | 右激波 | 3-族，\(u+a\) | 压缩间断 |
激波和膨胀波都属于声波族；区别在于一个是压缩，一个是释放。
移除隔膜后，左侧气体和右侧气体都会被加速到中间速度附近。
但是它们来自不同初始区域，密度一般不同。
因此，中间会保留一个物质分界面。
接触间断不是压力波，而是两团不同气体的分界面。
接触间断连接 \(L^*\) 和 \(R^*\)。
它沿流体速度运动：
\[ \frac{dx}{dt}=u^* \]
跨过接触间断：
\[ p_L^*=p_R^*,\qquad u_L^*=u_R^*,\qquad \rho_L^*\ne\rho_R^* \]
接触间断处，压力和速度必须连续。
如果压力不连续，界面两侧受力不同，界面会继续被加速。
如果速度不连续，两侧物质会相互穿透，这与接触面的定义不符。
| 变量 | 接触间断两侧 | 说明 |
|---|---|---|
| 压力 \(p\) | 连续 | 两侧力平衡 |
| 速度 \(u\) | 连续 | 两侧一起运动 |
| 密度 \(\rho\) | 可以跳跃 | 两团气体来源不同 |
| 温度 \(T\) | 可以跳跃 | 由状态方程决定 |
| 熵 | 可以跳跃 | 物质性质不同 |
判断接触间断最简单的方法：密度有跳跃，但压力和速度没有跳跃。
右侧低压气体受到中间高压区推动，发生压缩。
对于非线性可压缩流动，压缩扰动在传播过程中会逐渐陡化。
后面的压缩扰动可能追上前面的扰动，最终形成强间断。
压缩波不断汇聚，最后形成激波。
右激波连接 \(R^*\) 和 \(R\)。
激波跨越前后通常满足：
激波不是普通光滑解，而是守恒律意义下的弱解。
激波处变量不连续，因此不能只依赖光滑意义下的微分方程。
但是质量、动量和能量仍然必须守恒。
所以激波两侧状态必须满足 Rankine–Hugoniot 条件。
整体写法：
\[ s(\mathbf U_R-\mathbf U_L)=\mathbf F_R-\mathbf F_L \]
质量跳跃条件：
\[ s(\rho_R-\rho_L)=\rho_Ru_R-\rho_Lu_L \]
动量跳跃条件：
\[ s(\rho_Ru_R-\rho_Lu_L)= \rho_Ru_R^2+p_R-\rho_Lu_L^2-p_L \]
这里的 (s) 是激波传播速度，不是简单的 (u+a)。它由激波两侧守恒关系决定。
仅满足 Rankine–Hugoniot 条件还不够。
有些间断在数学上满足守恒跳跃条件，但物理上不允许。
典型例子是膨胀激波。
膨胀激波意味着熵降低，违反热力学第二定律。
物理激波要求：
\[ \Delta s>0 \]
也就是说，气体经过激波后熵增加。
守恒条件决定“能不能这样跳”；熵条件决定“物理上允不允许这样跳”。
左侧高压气体向右运动后，左侧区域压力降低。
这一过程不是压缩，而是释放和展开。
因此，左侧产生的是膨胀波，而不是激波。
膨胀波不是间断，而是一段连续展开的区域。
左膨胀波连接 \(L\) 和 \(L^*\)。
它对应第一族声学信息：
\[ \lambda_1=u-a \]
膨胀波内部变量连续变化：
\[ \rho,\quad u,\quad p \]
都不是突然跳跃，而是逐渐变化。
膨胀波可以理解为高压气体向外“松开”的过程。
在这个区域内：
压缩容易形成激波；释放通常形成连续膨胀波。
| 对比 | 膨胀波 | 激波 |
|---|---|---|
| 过程 | 释放 | 压缩 |
| 变量变化 | 连续变化 | 突然跳跃 |
| 熵变 | 近似等熵 | 熵增加 |
| 形态 | 一段区域 | 一个间断面 |
对某一侧初始状态 (K=L,R)，如果：
\[ p^*<p_K \]
说明这一侧气体压力降低，形成膨胀波。
如果：
\[ p^*>p_K \]
说明这一侧气体被压缩，形成激波。
对 Sod 激波管，通常有 \(p^*<p_L\)，所以左侧是膨胀波；同时 \(p^*>p_R\)，所以右侧是激波。
❶ 隔膜移除后，高压气体释放，低压气体被压缩。
❷ 三种波把流场分为 \(L,L^*,R^*,R\)
四个区域。
❸ 一维 Euler 方程有三族特征速度：\(u-a\)、\(u\)、\(u+a\)。
❹ \(L\to L^*\) 通常是膨胀波，\(L^*\to R^*\) 是接触间断，\(R^*\to R\) 通常是激波。
❺ 激波满足 Rankine–Hugoniot 条件，并且要满足熵条件。
❻ 膨胀波是连续变化区域，接触间断处压力和速度连续但密度可以跳跃。
掌握激波、接触间断和膨胀波的物理来源，再学习 Godunov、Roe、HLLC、WENO 等方法会自然很多。
延伸阅读：
建议有兴趣的同学课后继续阅读 E. F. Toro 的
Riemann Solvers and Numerical Methods for Fluid
Dynamics。
该书系统介绍了 Riemann 问题、Godunov 方法、Roe 格式、HLL/HLLC 格式以及
TVD/WENO 等高分辨率方法，是学习可压缩流动数值方法的重要参考书。
💻 代码演示
riemann2d_euler_contour_demo.m
Godunov 方法的核心思想是：
把每个网格界面附近看作一个局部 Riemann 问题， 通过界面处的波传播来构造数值通量。
在一个时间步内，数值解通常近似为分片常数：
设界面两侧状态为：\(U_L,\quad U_R\)，则界面附近可看作局部 Riemann 问题：
\[ U(x,0)= \begin{cases} U_L, & x<0\\ U_R, & x>0 \end{cases} \]
Godunov 方法就是利用这个局部问题的解，本质上是一种迎风思想，因为 Riemann 解自动包含了信息传播方向。
实际计算中，精确 Riemann 解通常代价较高，因此常用近似 Riemann 求解器：Roe、Osher、HLL、HLLC、AUSM 等。
UPWIND SCHEMES: FLUX VECTOR SPLITTING
迎风格式：Steger–Warming 通量向量分裂
数值格式不仅要逼近微分方程，
还要尊重流动中信息传播的方向。
Characteristics · Domain of dependence · Upwind differencing · Flux splitting
对于一维 Euler 方程：
\[ \frac{\partial U}{\partial t} +\frac{\partial F}{\partial x} =0 \]
通量 \(F\) 中同时包含不同传播方向的信息。
如果所有信息都用同一个差分方向处理，就无法体现流动中的上游影响关系。
通量分裂的目的，就是把向右传播和向左传播的信息分开。
因此可以写成：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ F=F^+ + F^- } $} \]
其中：
设通量 Jacobian 为：
\[ A=\frac{\partial F}{\partial U} \]
若 \(A\) 可对角化，则有：
\[ A=T\Lambda T^{-1} \]
其中：
\[ \Lambda= \begin{bmatrix} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3 \end{bmatrix} \]
对于一维 Euler 方程：
\[ \lambda_1=u,\qquad \lambda_2=u+c,\qquad \lambda_3=u-c \]
特征值的正负号，决定信息向哪个方向传播。
Steger–Warming 分裂的关键，是把每个特征值拆成正、负两部分：
\[ \lambda^+=\frac{\lambda+|\lambda|}{2}, \qquad \lambda^-=\frac{\lambda-|\lambda|}{2} \]
于是：
\[ \lambda=\lambda^++\lambda^-, \qquad \lambda^+\ge0, \qquad \lambda^-\le0 \]
对应关系为：
\[ \lambda^+=\lambda,\qquad \lambda^-=0 \]
信息向右传播。
\[ \lambda^+=0,\qquad \lambda^-=\lambda \]
信息向左传播。
将特征值矩阵分为：
\[ \Lambda=\Lambda^++\Lambda^- \]
其中：
\[ \Lambda^+=\operatorname{diag}(\lambda_1^+,\lambda_2^+,\lambda_3^+), \qquad \Lambda^-=\operatorname{diag}(\lambda_1^-,\lambda_2^-,\lambda_3^-) \]
由：
\[ A=T\Lambda T^{-1} \]
定义：
\[ A^+=T\Lambda^+T^{-1}, \qquad A^-=T\Lambda^-T^{-1} \]
于是：
\[ A=A^++A^- \]
Steger–Warming 分裂本质上是在特征空间中做迎风分裂。
对于齐次 Euler 方程，通量具有齐次性，可形式上写成：
\[ F=AU \]
于是定义：
\[ F^+=A^+U, \qquad F^-=A^-U \]
从而：
\[ F=F^++F^- \]
控制方程变为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial U}{\partial t} +\frac{\partial F^+}{\partial x} +\frac{\partial F^-}{\partial x} =0 } $} \]
先分裂通量，再分别迎风离散。
由于 \(F^+\) 对应向右传播的信息，因此使用后向差分：
\[ \left(\frac{\partial F^+}{\partial x}\right)_i \approx \frac{F_i^+-F_{i-1}^+}{\Delta x} \]
由于 \(F^-\) 对应向左传播的信息，因此使用前向差分：
\[ \left(\frac{\partial F^-}{\partial x}\right)_i \approx \frac{F_{i+1}^- - F_i^-}{\Delta x} \]
半离散格式为：
\[ \frac{dU_i}{dt}= -\frac{F_i^+-F_{i-1}^+}{\Delta x} -\frac{F_{i+1}^- - F_i^-}{\Delta x} \]
正向传播的通量用后向差分，负向传播的通量用前向差分。
Steger–Warming 通量向量分裂的关键步骤是什么？
A. 对压力项加入人工粘性
B. 将特征值按正负号分裂，并据此分裂通量
C. 在所有方向上统一使用中心差分
D. 只对密度方程进行迎风处理
参考判断：B
它的核心是：
先根据 Jacobian 的特征值判断信息传播方向，
再把通量分成 \(F^+\) 和 \(F^-\) 两部分进行迎风离散。
推荐阅读原始文献：J. L. Steger and R. F. Warming, “Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods,” Journal of Computational Physics, 40(2), 263–293, 1981.
本节讨论的迎风格式大多是一阶方法。
一阶迎风方法的主要优点是：
但它也有明显不足：
一阶迎风格式数值耗散较强，
会将激波、接触间断和剪切层过度抹宽。
因此，需要发展二阶甚至更高阶的迎风格式。
高分辨率格式的目标是：
光滑区保持高阶精度，间断附近避免非物理振荡。
为描述振荡程度，引入总变差：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ TV(u)=\sum_i |u_{i+1}-u_i| } $} \]
若数值格式满足：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ TV(u^{n+1})\le TV(u^n) } $} \]
则称为 TVD 格式：
\[ \colorbox{#9ACD32}{$\displaystyle \color{black}{ \text{Total Variation Diminishing} } $} \]
TVD 的物理含义是：数值解在推进过程中不应产生新的非物理振荡。
二阶格式的风险来自高阶修正项。通量限制器的作用是：
根据局部光滑程度，决定允许多少高阶修正。
一般可写成：
\[ F_{\text{limited}}=F_{\text{1st-order}}+\phi(r)\left(F_{\text{2nd-order}}-F_{\text{1st-order}}\right) \]
其中 \(r\) 表示相邻梯度比，\(\phi(r)\) 为限制器。
典型行为是：
常见限制器包括：Minmod、Superbee、Van Leer、Van Albada。
高分辨率格式的核心不是“越高阶越好”， 而是在精度和稳定性之间实现自适应平衡。
许多 CFD 计算采用迭代或时间推进方法。
在迭代过程中，误差可以分为：
常规迭代方法往往有一个特点：
高频误差衰减较快，低频误差衰减较慢。
因此，即使局部振荡已经被抹平，整体误差仍可能长时间存在，导致收敛缓慢。
在一维网格上，误差可以看成不同波长成分的叠加。
最短波长约为：
\[ \lambda_{\min}=2\Delta x \]
最长波长约为：
\[ \lambda_{\max}=L \]
其中：
在细网格上看起来是低频误差的成分，
转移到粗网格后可能变成较高频误差，因而更容易被迭代过程消除。
多重网格方法的流程可以概括为：
在细网格上进行若干步迭代，消除高频误差
将残差或误差转移到较粗网格
在粗网格上继续迭代，快速消除低频误差
将粗网格修正量插值回细网格
重复上述过程，直到细网格收敛
多重网格的关键不是“网格越粗越好”，
而是让不同尺度的误差在合适的网格上被快速消除。
多重网格方法常用于加速稳态 CFD 计算。
其优点包括：
但它也有一些实现难点：
二阶迎风格式比一阶迎风格式精度更高，为什么仍需要通量限制器？
A. 因为二阶格式一定不稳定
B. 因为二阶格式在间断附近可能产生非物理振荡
C. 因为一阶格式无法处理激波
D. 因为限制器可以消除所有数值误差
参考判断：B
限制器的作用不是简单增加耗散，
而是根据局部光滑程度自适应控制高阶修正项。
多重网格方法为什么可以提高稳态 CFD 计算的收敛速度？
A. 因为粗网格一定比细网格精确
B. 因为它完全不需要迭代
C. 因为不同尺度误差在不同网格上更容易被消除
D. 因为它可以避免求解控制方程
参考判断：C
细网格适合表示最终精细解，
粗网格则有助于更快消除低频误差。
CHAPTER 4
计算进阶：从格式到思想
下一步可以继续进入高分辨率格式、Riemann 求解器、隐式算法、湍流模拟与复杂几何 CFD。
Conservation · Characteristics · Upwinding · High resolution · Multigrid
PHYSICS-INFORMED NEURAL NETWORKS
物理信息神经网络：从数据拟合到方程约束
传统神经网络主要学习数据规律；
PINN 在学习数据的同时，把 控制方程、边界条件和初始条件
也写进训练目标。
普通神经网络通常学习：
\[ 输入\longrightarrow输出 \]
例如输入 \((x,t)\)，输出 \(u(x,t)\)、\(T(x,t)\) 或 \(p(x,t)\)。
但如果只依赖数据，可能出现两个问题：
PINN 的出发点：让神经网络不仅拟合数据，还要尽量满足物理方程。
考虑一个偏微分方程：
\[ \mathcal{N}[u](x,t)=0 \]
PINN 用神经网络近似未知解：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ u_\theta(x,t)\approx u(x,t)} $} \]
将网络输出代入控制方程，得到残差：
\[ r_\theta(x,t)=\mathcal{N}[u_\theta](x,t) \]
训练目标不只是：
\[ u_\theta\approx u_{data} \]
还要尽量满足：
\[ r_\theta\approx0 \]
PINN = 神经网络函数近似 + 物理方程残差约束。
考虑一维线性对流方程：
\[ u_t+a u_x=0 \]
它表示波形以速度 \(a\) 沿 \(x\) 方向传播。
令神经网络输出：
\[ u_\theta(x,t) \]
则方程残差为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ r_\theta=u_{\theta,t}+a u_{\theta,x} } $} \]
训练时希望在许多采样点上满足：
\[ r_\theta\approx0 \]
这里的导数不是用差分近似，而是通过神经网络的自动微分得到。
传统数值方法中，导数常通过差分近似：
\[ \frac{\partial u}{\partial x}\approx\frac{u_{i+1}-u_i}{\Delta x} \]
在 PINN 中，\(u_\theta(x,t)\) 是一个可微函数，因此可以直接计算：
\[ \frac{\partial u_\theta}{\partial x},\qquad \frac{\partial u_\theta}{\partial t},\qquad \frac{\partial^2u_\theta}{\partial x^2} \]
以MATLAB语言为例，通过 dlarray 和
dlgradient 实现自动微分。
PINN 的损失函数通常由几部分组成：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \mathcal{L}=\mathcal{L}_{data}+\lambda_f\mathcal{L}_{PDE}+\lambda_b\mathcal{L}_{BC}+\lambda_i\mathcal{L}_{IC} } $} \]
其中：
\[ \mathcal{L}_{PDE}=\frac{1}{N_f}\sum |r_\theta|^2 \]
\[ \mathcal{L}_{IC}=\frac{1}{N_i}\sum |u_\theta(x,0)-u_0(x)|^2 \]
\[ \mathcal{L}_{BC}=\frac{1}{N_b}\sum |u_\theta-u_b|^2 \]
PINN 的基本流程可以概括为：
输入：\((x,t)\)
输出：\(u_\theta(x,t)\)
自动微分计算导数，代入控制方程。
最小化总损失，得到近似解。
本代码求解：
\[ u_t+a u_x=0,\qquad x\in[0,1],\quad t\in[0,1] \]
周期边界：
\[ u(0,t)=u(1,t) \]
初始条件：
\[ u(x,0)=\sin(2\pi x) \]
对应精确解：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ u(x,t)=\sin[2\pi(x-at)] } $} \]
代码中的主要训练参数为：
三类训练点分别为：
\(N_f\) 用于 PDE 残差，\(N_i\) 用于初始条件，\(N_b\) 用于周期边界条件。
网络输入为：
\[ XT=\begin{bmatrix}x\\t\end{bmatrix} \]
代码中写成：
XT_f = dlarray(single([x_f; t_f ]),'CB');
XT_i = dlarray(single([x_i; t_i ]),'CB');
U_i = dlarray(single(u_i),'CB');
XT_b0 = dlarray(single([x_b0; t_b ]),'CB');
XT_b1 = dlarray(single([x_b1; t_b ]),'CB');其中 CB 表示： C：通道维度，这里是 \(x,t\)；
B：批量维度，即采样点编号
训练时必须通过：
来开启自动微分跟踪。
代码中建立网络：
含义是：
网络主要由全连接层和 tanhLayer 组成：
这里使用 tanh，是因为 PINN
需要网络输出关于输入变量可微。
在代码中，总损失为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ \mathcal{L}=\mathcal{L}_{PDE}+10\mathcal{L}_{IC}+\mathcal{L}_{BC} } $} \]
对应 MATLAB 语句：
三项含义为：
lossPDE：约束 \(u_t+c
u_x=0\)lossIC：约束 \(u_\theta(x,0)=\sin(2\pi x)\)lossBC：约束 \(u_\theta(0,t)=u_\theta(1,t)\)初始条件前乘以 10，是为了加强对初始波形的约束。
在 modelLoss 中，先计算配点处网络输出：
再计算输出对输入的导数：
方程残差为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ r=u_t+a u_x } $} \]
训练循环的核心代码为：
for epoch = 1:numEpochs
[loss,gradients,lossPDE,lossIC,lossBC] = ...
dlfeval(@modelLoss,net,XT_f,XT_i,U_i,XT_b0,XT_b1,c);
[net,trailingAvg,trailingAvgSq] = ...
adamupdate(net,gradients,trailingAvg,trailingAvgSq,...
epoch,learnRate);
end每一轮主要做三件事：
如果 loss 持续下降，说明网络逐渐同时满足方程、初始条件和边界条件。
训练完成后，在规则网格上预测：
[X,T] = meshgrid(x,t);
XT_test = dlarray(single([X(:)'; T(:)']),'CB');
U_pred = forward(net,XT_test);精确解为：
误差定义为：
\[ \colorbox{yellow}{$\displaystyle \color{black}{ e(x,t)=|u_\theta(x,t)-u(x,t)| } $} \]
传统 CFD 的思路通常是：
\[ \text{先离散方程，再求解离散系统} \]
典型方法包括：
PINN 的思路是：
\[ \text{先构造函数近似，再用方程残差训练网络} \]
PINN 不是传统 CFD 的简单替代品，更像是一种“数据 + 方程”的建模框架。
PINN 常见应用可以分为两类：
已知方程、参数和边界条件，求解物理场。
例如：已知 ()，求 (u(x,t))。
已知部分观测数据，反推未知参数或未知项。
例如：已知 (u(x,t))，反推 ()。
流体问题中，PINN 常用于：
一维热传导方程：
\[ T_t=\alpha T_{xx} \]
PINN 残差：
\[ r_\theta=T_{\theta,t}-\alpha T_{\theta,xx} \]
Burgers 方程：
\[ u_t+u u_x=\nu u_{xx} \]
PINN 残差：
\[ r_\theta=u_{\theta,t}+u_\theta u_{\theta,x}-\nu u_{\theta,xx} \]
不可压缩 Navier–Stokes 方程：
\[ \nabla\cdot\mathbf{u}=0 \]
\[ \mathbf{u}_t+\mathbf{u}\cdot\nabla\mathbf{u} =-\frac{1}{\rho}\nabla p+\nu\nabla^2\mathbf{u} \]
PINN 也存在明显挑战：
理解 PINN 时要避免一个误区：
它不是替代所有传统数值方法的通用工具，而是将物理方程约束引入神经网络训练的建模方法。
CFD模块四章围绕一个基本问题展开：
如何把流体力学方程可靠地变成计算机可求解的问题？
CFD 研究的是流场变量在空间和时间中的分布：
\[ \rho,\quad u,\quad v,\quad w,\quad p,\quad T \]
流体力学控制方程来自：
微分形式适合局部分析，积分形式更适合理解守恒量跨边界的变化。
\[ \frac{\partial}{\partial t}\int_\Omega U\,d\Omega + \int_{\partial\Omega} F\cdot n\,dS = \int_\Omega S\,d\Omega \]
可压缩流动常写成守恒型：
\[ \frac{\partial U}{\partial t} +\frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} +\frac{\partial H}{\partial z}=S \]
以二维 Euler 方程为例：
\[ U=\begin{bmatrix} \rho\\ \rho u\\ \rho v\\ \rho E \end{bmatrix} \]
实际观察和物理解释常使用：
\[ \rho,\quad u,\quad v,\quad p,\quad T,\quad M \]
数值推进常针对守恒变量，后处理和边界条件常需要恢复基本变量。
不同方程具有不同的信息传播方式：
方程性质决定：
一维 Euler 方程具有三族特征速度：
\[u-a,\qquad u,\qquad u+a \]
其中 \(a\) 为声速。
特征速度影响边界条件、信息传播方向和稳定性限制。
激波是薄层内发生剧烈变化的压缩波，在宏观计算中常表现为不连续。
对一维守恒律：
\[ \frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=0 \]
激波速度 \(s\) 满足：
\[ s[U]=[F] \]
激波可以是不连续的，但激波前后仍必须满足积分形式的守恒律。
边界条件决定计算区域如何与外部环境交换信息。
亚声速边界和超声速边界的信息传播方向不同，给定变量的个数和方式也不同（回忆拉瓦尔喷管问题）。
\[ U_i^n\approx U(x_i,t^n) \]
网格越细通常越准确，但计算量更大，稳定性限制也可能更严格。
用相邻网格点上的函数值近似导数。
\[ \frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1}-u_i}{\Delta x} \]
\[ \frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_i-u_{i-1}}{\Delta x} \]
\[ \frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1}-u_{i-1}}{2\Delta x} \]
不同差分格式具有不同的精度、耗散、色散和稳定性。
新时刻未知量可以直接由旧时刻计算得到：
\[ U^{n+1}=\mathcal{L}(U^n) \]
优点是实现简单，缺点是时间步长受稳定性限制较强。
新时刻未知量同时出现在方程两边：
\[ A U^{n+1}=b \]
优点是稳定性较好，缺点是每一步需要求解代数方程组。
将精确解代入离散格式后，剩余项称为截断误差。
若误差满足：
\[ \tau=O(\Delta x^p)+O(\Delta t^q) \]
则空间为 \(p\) 阶精度，时间为 \(q\) 阶精度。
精度阶数说明误差随网格加密和时间步减小而降低的速度。
当 \(\Delta x,\Delta t\to 0\) 时，离散格式应趋近原微分方程。
计算过程中的误差不能无限增长。
当网格加密时，数值解应趋近精确解。
把误差分解为 Fourier 模态：
\[ \varepsilon_i^n=G^n e^{\mathrm{i}k x_i} \]
其中 \(G\) 为放大因子。
若所有波数下均满足：
\[ |G|\le 1 \]
则误差不会随时间步无限放大。
该方法主要用于线性、常系数、均匀网格问题的稳定性分析。
对于一维对流方程：
\[ \frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0 \]
CFL 数为：
\[ \mathrm{CFL}=\frac{|a|\Delta t}{\Delta x} \]
一个时间步内，物理信息传播距离不能超过数值格式能够有效感知的范围。
CFL 条件是很多显式格式稳定计算的必要限制。
复杂物理区域可映射到规则计算区域：
\[ (x,y)\longleftrightarrow(\xi,\eta) \]
坐标变换会引入度量项，方程形式和数值离散都需要相应处理。
Lax–Wendroff 方法利用 Taylor 展开进行时间推进，再用方程本身把时间导数转化为空间导数。
\[ \frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0 \]
典型格式为：
\[u_i^{n+1}=u_i^n-\frac{\lambda}{2}(u_{i+1}^n-u_{i-1}^n) +\frac{\lambda^2}{2}(u_{i+1}^n-2u_i^n+u_{i-1}^n) \]
其中 \(\lambda=a\Delta t/\Delta x\)。
以前向差分为例：
\[ \bar{u}_i=u_i^n-\frac{a\Delta t}{\Delta x}(u_{i+1}^n-u_i^n) \]
再用反向差分校正：
\[u_i^{n+1}=\frac{1}{2}\left[u_i^n+\bar{u}_i-\frac{a\Delta t}{\Delta x}(\bar{u}_i-\bar{u}_{i-1})\right] \]
数值格式使波幅随时间衰减，常表现为间断被抹宽、峰值降低。
不同波数成分传播速度不一致，常表现为相位误差和间断附近振荡。
人工粘性用于在强梯度或激波附近增加耗散，抑制非物理振荡。
人工粘性可以稳定计算，但会降低局部分辨率，需要控制其作用范围和强度。
守恒形式的核心是保证离散意义下的质量、动量和能量平衡。对于含激波流动，守恒形式直接关系到激波位置和强度是否正确。
激波管问题包含膨胀波、接触间断和激波，是理解可压缩流动波系结构的重要例子。
\[ u_\theta(x,t)\approx u(x,t) \]
PINN 用神经网络表示未知解，并在训练中同时约束：数据、方程、初始/边界条件
导数由自动微分得到。