流体力学核心贯通课 II

2025–2026 学年第二学期

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

任杰

综合算例准一维喷管流动

参考书

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

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

流动描述：收缩–扩张喷管中的等熵流动

流体从储气室流入喷管：

入口条件：
$\color{black}{p_0,T_0}$ （初速为零时=总压、总温）

流动沿喷管轴向发展

速度逐渐增加

喷管结构：

收缩段 → 流动加速
喉部（最小截面）：$\color{black}{M=1}$
扩张段 → 超声流

流动描述：收缩–扩张喷管中的等熵流动

关键假设：

准一维流动
- 尽管真实流动是二维
- 截面平均处理
- 忽略横向（$\color{black}{y}$ 向）变化

等熵（无粘、无热传递）

变量仅随 $\color{black}{x}$ 变化

流动描述：收缩–扩张喷管中的等熵流动

曲线 A 和 B：背压较高、喉部未达到声速、整个喷管为亚声速流动

曲线 C（临界状态）：喉部达到声速 $\color{black}{M=1}$、质量流量达到最大、流动刚好阻塞

曲线 D–F（激波在喷管内）：喉部已阻塞、扩张段超声、在扩张段内部形成正激波、激波后流动转为亚声速

曲线 G：激波位于出口、无法用单一正常激波匹配、出口外形成复杂斜激波结构

曲线 H（设计工况）：喉部 $\color{black}{M=1}$、扩张段完全超声、无激波

曲线 I：背压进一步降低、喷流完全超声、在喷管外形成膨胀波/压缩波结构

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❶ 从质量守恒出发

准一维定常流动中，质量流量保持不变：

\[ \dot{m}=\rho V A=\mathrm{const} \]

因此有：

\[ A=\frac{\dot{m}}{\rho V} \]

将速度写成马赫数形式：

\[ V=Ma \]

其中声速为：

\[ a=\sqrt{\gamma RT} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❷ 引入状态方程

由

\[ V=Ma,\qquad a=\sqrt{\gamma RT} \]

可得：

\[ V=M\sqrt{\gamma RT} \]

代入质量守恒关系：

\[ A=\frac{\dot{m}}{\rho M\sqrt{\gamma RT}} \]

利用状态方程：

\[ p=\rho RT \Rightarrow \rho=\frac{p}{RT} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❸ 整理为压力和温度形式

由

\[ A=\frac{\dot{m}}{\rho M\sqrt{\gamma RT}} \]

代入

\[ \rho=\frac{p}{RT} \]

得到：

\[ A=\frac{\dot{m}\sqrt{RT}}{pM\sqrt{\gamma}}=\frac{\dot{m}}{M\sqrt{\gamma}}\frac{\sqrt{RT}}{p} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❹ 代入等熵关系

对于等熵流动，有：

\[ \frac{T}{T_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-1} \]

\[ \frac{p}{p_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-\gamma/(\gamma-1)} \]

即：

\[ T=T_0 \left(1+\frac{\gamma-1}{2}M^2\right)^{-1} \]

\[ p=p_0 \left(1+\frac{\gamma-1}{2}M^2\right)^{-\gamma/(\gamma-1)} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❺ 得到面积关于马赫数的表达式

将等熵关系代入

\[ A=\frac{\dot{m}}{M\sqrt{\gamma}}\frac{\sqrt{RT}}{p}\Rightarrow A= \frac{\dot{m}}{p_0} \sqrt{\frac{RT_0}{\gamma}} \frac{1}{M} \left(1+\frac{\gamma-1}{2}M^2\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

也即：

\[ A=C\frac{1}{M} \left(1+\frac{\gamma-1}{2}M^2\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

其中：

\[ C=\frac{\dot{m}}{p_0}\sqrt{\frac{RT_0}{\gamma}} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❻ 引入临界面积

在喉部达到声速时：

\[ M=1 \]

对应面积记为：

\[ A=A^* \]

由

\[ A=C \frac{1}{M} \left(1+\frac{\gamma-1}{2}M^2\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

令 $\color{black}{M=1}$，得到：

\[ A^*=C \left(1+\frac{\gamma-1}{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}\Rightarrow A^*=C\left(\frac{\gamma+1}{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❼ 临界面积表达式

由

\[ A=C \frac{1}{M} \left(1+\frac{\gamma-1}{2}M^2\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

和

\[ A^*=C \left(\frac{\gamma+1}{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}} \]

取比值。

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❽ 最终关系式

取比值可得：

\[ \frac{A}{A^*}=\frac{1}{M} \left[ \frac{ 1+\frac{\gamma-1}{2}M^2 }{ \frac{\gamma+1}{2} } \right]^{\frac{\gamma+1}{2(\gamma-1)}} \]

整理为：

\[ \frac{A}{A^*}= \frac{1}{M} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M^2 \right) \right]^{\frac{\gamma+1}{2(\gamma-1)}} \]

两边平方，即得：

\[ \left(\frac{A}{A^*}\right)^2= \frac{1}{M^2} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M^2 \right) \right]^{\frac{\gamma+1}{\gamma-1}} \]

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❾ 物理含义

\[ \left(\frac{A}{A^*}\right)^2= \frac{1}{M^2} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M^2 \right) \right]^{\frac{\gamma+1}{\gamma-1}} \]

其中：

$\color{black}{A^*}$：喉部面积，对应 $\color{black}{M=1}$

👉 面积变化完全决定马赫数变化

👉 对于一个给定面积比 $\color{black}{A_e/A^*}$ 的喷管，要建立设计工况的亚声—超声等熵流动，所需的压强比必须是一个唯一确定的，等熵关系计算出的那个值。

这个压强比本质上是施加在流动上的一个边界条件；在实验中，它通常由入口处的高压气源和/或出口处的抽真空装置提供。

流动描述：收缩–扩张喷管中的等熵流动

等熵关系

等熵关系：压力比

\[ \frac{p}{p_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-\gamma/(\gamma-1)} \]

等熵关系：密度比

\[ \frac{\rho}{\rho_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-1/(\gamma-1)} \]

等熵关系：温度比

\[ \frac{T}{T_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-1} \]

流动描述：收缩–扩张喷管中的等熵流动

等熵关系的物理含义

这些关系表明，随着流速增加：

压力下降
密度下降
温度下降

一维等熵流动面积–速度关系： \[ \frac{dA}{A}=(M^2-1)\frac{dV}{V} \] 超声速区面积增大并不“减速”，因为此时膨胀导致的密度下降更强，为满足 $\rho V A=\mathrm{const}$，速度 $V$ 必须继续增大。

CFD求解：MacCormack方法

本节中，我们将使用 MacCormack 显式方法对准一维喷管流进行数值求解。

需要强调：

👉 该问题本身存在解析解（见上面的分析）

👉 这里进行数值求解的目的：

展示 CFD 方法
CFD 对于该问题属于 over-kill

理解数值流程

为复杂问题打基础

控制方程推导：从积分形式出发

连续方程 ❶ 积分形式与控制体

连续方程积分形式：

\[ \frac{\partial}{\partial t}\iiint \rho\,dV+\iint \rho \mathbf{V}\cdot d\mathbf{S}=0 \]

建立控制体，取喷管微元：

长度：$\color{black}{dx}$
面积：$\color{black}{A(x)}$

体积：

\[ dV=A\,dx \]

控制方程推导：从积分形式出发

连续方程 ❷ 体积分项

\[ \frac{\partial}{\partial t}\iiint \rho\,dV+\iint \rho \mathbf{V}\cdot d\mathbf{S}=0 \]

体积分项为：

\[ \frac{\partial}{\partial t}\iiint \rho\,dV \]

对准一维喷管微元，有：

\[ dV=A\,dx \]

因此：

\[ \frac{\partial}{\partial t}\iiint \rho\,dV = \frac{\partial(\rho A)}{\partial t}dx \]

控制方程推导：从积分形式出发

连续方程 ❸ 表面积分项

\[ \frac{\partial}{\partial t}\iiint \rho\,dV+\iint \rho \mathbf{V}\cdot d\mathbf{S}=0 \]

表面积分项（入口）：$-\rho V A$

表面积分项（出口）：$(\rho+d\rho)(V+dV)(A+dA)$

忽略高阶小量后，出口项展开为：

\[ \rho V A+\rho V\,dA+\rho A\,dV+A V\,d\rho \]

因此表面积分项为：

\[ \iint \rho \mathbf{V}\cdot d\mathbf{S}= \frac{\partial(\rho V A)}{\partial x}dx \]

控制方程推导：从积分形式出发

连续方程 ❹ 守恒形式

将体积分项和表面积分项代入：

\[ \frac{\partial}{\partial t}\iiint \rho\,dV+\iint \rho \mathbf{V}\cdot d\mathbf{S}=0 \]

得到：

\[ \frac{\partial(\rho A)}{\partial t}dx + \frac{\partial(\rho V A)}{\partial x}dx =0 \]

两边除以 $dx$，得到准一维连续方程的守恒形式：

\[ \frac{\partial(\rho A)}{\partial t} +\frac{\partial(\rho V A)}{\partial x}=0 \]

控制方程推导：从积分形式出发

连续方程 ❺ 非守恒形式

由守恒形式：

\[ \frac{\partial(\rho A)}{\partial t} +\frac{\partial(\rho V A)}{\partial x}=0 \]

由于：

\[ A=A(x),\qquad \frac{\partial A}{\partial t}=0\Rightarrow\frac{\partial(\rho A)}{\partial t} =A\frac{\partial\rho}{\partial t} \]

同时：

\[ \frac{\partial(\rho V A)}{\partial x} =A\frac{\partial(\rho V)}{\partial x} +\rho V\frac{dA}{dx} \]

控制方程推导：从积分形式出发

连续方程 ❻ 非守恒形式

代入连续方程：

\[ A\frac{\partial\rho}{\partial t} +A\frac{\partial(\rho V)}{\partial x} +\rho V\frac{dA}{dx}=0\Rightarrow \frac{\partial\rho}{\partial t} +\frac{\partial(\rho V)}{\partial x} +\rho V\frac{1}{A}\frac{dA}{dx}=0 \]

利用：

\[ \frac{\partial(\rho V)}{\partial x} =V\frac{\partial\rho}{\partial x} +\rho\frac{\partial V}{\partial x},\qquad \frac{1}{A}\frac{dA}{dx} =\frac{\partial(\ln A)}{\partial x} \]

得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial \rho}{\partial t} +V\frac{\partial \rho}{\partial x} +\rho\frac{\partial V}{\partial x} +\rho V\frac{\partial(\ln A)}{\partial x}=0 } $} \]

控制方程推导：从积分形式出发

动量方程

动量方程积分形式：

\[ \frac{\partial}{\partial t}\iiint(\rho u)\,dV +\iint(\rho u\mathbf{V})\cdot d\mathbf{S}= -\iint p\,dS_x \]

对于准一维喷管流，轴向速度记为：

\[ u=V \]

类似连续方程的处理，可得：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial(\rho V^2A)}{\partial x}= -\frac{\partial(pA)}{\partial x} +p\frac{\partial A}{\partial x} \]

控制方程推导：从积分形式出发

动量方程

由：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial(\rho V^2A)}{\partial x}= -\frac{\partial(pA)}{\partial x} +p\frac{\partial A}{\partial x} \]

展开并结合连续方程得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial V}{\partial t} +V\frac{\partial V}{\partial x}= -\frac{1}{\gamma} \left( \frac{\partial T}{\partial x} +\frac{T}{\rho}\frac{\partial\rho}{\partial x} \right) } $} \]

控制方程推导：从积分形式出发

能量方程 ❶ 积分形式

能量方程积分形式：

\[ \frac{\partial}{\partial t} \iiint \rho\left(e+\frac{V^2}{2}\right)dV +\iint \rho\left(e+\frac{V^2}{2}\right)\mathbf{V}\cdot dS= -\iint p\mathbf{V}\cdot dS \]

最终得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial T}{\partial t} +V\frac{\partial T}{\partial x}= -(\gamma-1)T \left[ \frac{\partial V}{\partial x} +V\frac{\partial(\ln A)}{\partial x} \right] } $} \]

控制方程推导：从积分形式出发

最终控制方程汇总

连续方程：

\[ \frac{\partial \rho}{\partial t} +V\frac{\partial \rho}{\partial x} +\rho\frac{\partial V}{\partial x} +\rho V\frac{\partial(\ln A)}{\partial x}=0 \]

动量方程：

\[ \frac{\partial V}{\partial t} +V\frac{\partial V}{\partial x}= -\frac{1}{\gamma} \left( \frac{\partial T}{\partial x} +\frac{T}{\rho}\frac{\partial \rho}{\partial x} \right) \]

能量方程：

\[ \frac{\partial T}{\partial t} +V\frac{\partial T}{\partial x}= -(\gamma-1)T \left[ \frac{\partial V}{\partial x} +V\frac{\partial(\ln A)}{\partial x} \right] \]

最终得到以$\rho$、$V$、$T$为基本变量的非守恒形式

MacCormack方法

预测步（forward）❶ 密度方程

\[ \frac{\partial \rho}{\partial t} +V\frac{\partial \rho}{\partial x} +\rho\frac{\partial V}{\partial x} +\rho V\frac{\partial(\ln A)}{\partial x}=0 \]

预测步采用前向差分计算当前时刻的时间导数。

密度方程：

\[ \left(\frac{\partial \rho}{\partial t}\right)_i^{t}= -V_i^{t}\frac{\rho_{i+1}^{t}-\rho_i^{t}}{\Delta x} -\rho_i^{t}\frac{V_{i+1}^{t}-V_i^{t}}{\Delta x} -\rho_i^{t}V_i^{t}\frac{\ln A_{i+1}-\ln A_i}{\Delta x} \]

MacCormack方法

预测步（forward）❷ 速度方程

\[ \frac{\partial V}{\partial t} +V\frac{\partial V}{\partial x}= -\frac{1}{\gamma} \left( \frac{\partial T}{\partial x} +\frac{T}{\rho}\frac{\partial \rho}{\partial x} \right) \]

速度方程：

\[ \left(\frac{\partial V}{\partial t}\right)_i^{t}= -V_i^{t}\frac{V_{i+1}^{t}-V_i^{t}}{\Delta x} -\frac{1}{\gamma} \left[ \frac{T_{i+1}^{t}-T_i^{t}}{\Delta x} +\frac{T_i^{t}}{\rho_i^{t}} \frac{\rho_{i+1}^{t}-\rho_i^{t}}{\Delta x} \right] \]

MacCormack方法

预测步（forward）❸ 温度方程

\[ \frac{\partial T}{\partial t} +V\frac{\partial T}{\partial x}= -(\gamma-1)T \left[ \frac{\partial V}{\partial x} +V\frac{\partial(\ln A)}{\partial x} \right] \]

温度方程：

\[ \left(\frac{\partial T}{\partial t}\right)_i^{t}= -V_i^{t}\frac{T_{i+1}^{t}-T_i^{t}}{\Delta x} -(\gamma-1)T_i^{t} \left[ \frac{V_{i+1}^{t}-V_i^{t}}{\Delta x} +V_i^{t}\frac{\ln A_{i+1}-\ln A_i}{\Delta x} \right] \]

MacCormack方法

预测步（forward）❹ 更新预测值

利用预测步得到的时间导数，更新预测值：

\[ \overline{\rho}_i^{t+\Delta t}= \rho_i^t + \left(\frac{\partial \rho}{\partial t}\right)_i^t\Delta t \]

\[ \overline{V}_i^{t+\Delta t}= V_i^t + \left(\frac{\partial V}{\partial t}\right)_i^t\Delta t \]

\[ \overline{T}_i^{t+\Delta t}= T_i^t + \left(\frac{\partial T}{\partial t}\right)_i^t\Delta t \]

MacCormack方法

校正步（backward）❶ 密度方程

\[ \frac{\partial \rho}{\partial t} +V\frac{\partial \rho}{\partial x} +\rho\frac{\partial V}{\partial x} +\rho V\frac{\partial(\ln A)}{\partial x}=0 \]

校正步采用预测值和后向差分重新计算时间导数。

\[ \left(\frac{\partial \overline{\rho}}{\partial t}\right)_i^{t+\Delta t}= -\overline{\rho}_i^{t+\Delta t} \frac{\overline{V}_i^{t+\Delta t}-\overline{V}_{i-1}^{t+\Delta t}}{\Delta x} -\overline{\rho}_i^{t+\Delta t}\overline{V}_i^{t+\Delta t} \frac{\ln A_i-\ln A_{i-1}}{\Delta x} -\overline{V}_i^{t+\Delta t} \frac{\overline{\rho}_i^{t+\Delta t}-\overline{\rho}_{i-1}^{t+\Delta t}}{\Delta x} \]

MacCormack方法

\[ \frac{\partial V}{\partial t} +V\frac{\partial V}{\partial x}= -\frac{1}{\gamma} \left( \frac{\partial T}{\partial x} +\frac{T}{\rho}\frac{\partial \rho}{\partial x} \right) \]

校正步（backward）❷ 速度方程

\[ \left(\frac{\partial \overline{V}}{\partial t}\right)_i^{t+\Delta t}= -\overline{V}_i^{t+\Delta t} \frac{\overline{V}_i^{t+\Delta t}-\overline{V}_{i-1}^{t+\Delta t}}{\Delta x} -\frac{1}{\gamma} \left[ \frac{\overline{T}_i^{t+\Delta t}-\overline{T}_{i-1}^{t+\Delta t}}{\Delta x} + \frac{\overline{T}_i^{t+\Delta t}}{\overline{\rho}_i^{t+\Delta t}} \frac{\overline{\rho}_i^{t+\Delta t}-\overline{\rho}_{i-1}^{t+\Delta t}}{\Delta x} \right] \]

MacCormack方法

校正步（backward）❸ 温度方程

\[ \frac{\partial T}{\partial t} +V\frac{\partial T}{\partial x}= -(\gamma-1)T \left[ \frac{\partial V}{\partial x} +V\frac{\partial(\ln A)}{\partial x} \right] \]

\[ \left(\frac{\partial \overline{T}}{\partial t}\right)_i^{t+\Delta t}= -\overline{V}_i^{t+\Delta t} \frac{\overline{T}_i^{t+\Delta t}-\overline{T}_{i-1}^{t+\Delta t}}{\Delta x} -(\gamma-1)\overline{T}_i^{t+\Delta t} \left[ \frac{\overline{V}_i^{t+\Delta t}-\overline{V}_{i-1}^{t+\Delta t}}{\Delta x} +\overline{V}_i^{t+\Delta t} \frac{\ln A_i-\ln A_{i-1}}{\Delta x} \right] \]

MacCormack方法

平均时间导数

预测步和校正步得到两个时间导数，对它们取平均：

\[ \left(\frac{\partial \rho}{\partial t}\right)_{\mathrm{av}}= \frac{1}{2} \left[ \left(\frac{\partial \rho}{\partial t}\right)_i^{t} + \left(\frac{\partial \overline{\rho}}{\partial t}\right)_i^{t+\Delta t} \right] \]

\[ \left(\frac{\partial V}{\partial t}\right)_{\mathrm{av}}= \frac{1}{2} \left[ \left(\frac{\partial V}{\partial t}\right)_i^{t} + \left(\frac{\partial \overline{V}}{\partial t}\right)_i^{t+\Delta t} \right] \]

\[ \left(\frac{\partial T}{\partial t}\right)_{\mathrm{av}}= \frac{1}{2} \left[ \left(\frac{\partial T}{\partial t}\right)_i^{t} + \left(\frac{\partial \overline{T}}{\partial t}\right)_i^{t+\Delta t} \right] \]

MacCormack方法

最终更新

用平均时间导数完成最终更新：

\[ \rho_i^{t+\Delta t}= \rho_i^t + \left(\frac{\partial \rho}{\partial t}\right)_{\mathrm{av}}\Delta t \]

\[ V_i^{t+\Delta t}= V_i^t + \left(\frac{\partial V}{\partial t}\right)_{\mathrm{av}}\Delta t \]

\[ T_i^{t+\Delta t}= T_i^t + \left(\frac{\partial T}{\partial t}\right)_{\mathrm{av}}\Delta t \]

MacCormack方法

时间步长：$\Delta t$ 应该取多大？❶

控制方程组在时间上是双曲型的，因此显式格式存在稳定性约束。

对准一维可压缩流，扰动信息的传播速度包括：

\[ V-a,\qquad V,\qquad V+a \]

其中：

$V$：流体质点随流动输运的速度
$a$：声速
$V+a$：沿下游传播的声学扰动速度
$V-a$：沿上游传播的声学扰动速度

因此，控制时间步长的不是单纯的 $\color{black}{V}$，而是最快的信息传播速度 $\color{black}{V+a}$。

MacCormack方法

时间步长：$\Delta t$ 应该取多大？❷

显式格式中，一个时间步内，数值信息主要在相邻网格点之间传递。

因此，数值格式能够分辨的信息传播速度量级为：

\[ \frac{\Delta x}{\Delta t} \]

物理中最快的信息传播速度为：$V+a$

为了使物理信息不跨过数值格式能够分辨的范围，需要满足：

\[ V+a\leq\frac{\Delta x}{\Delta t} \]

整理得到：

\[ \Delta t\leq\frac{\Delta x}{a+V} \]

MacCormack方法

时间步长：$\Delta t$ 应该取多大？❸

定义：

\[ C=\frac{(a+V)\Delta t}{\Delta x} \]

于是时间步长可写为：

\[ \Delta t= C\frac{\Delta x}{a+V} \]

其中： $C$ 为 CFL 数

对显式格式，要求：

\[ C\leq 1 \]

MacCormack方法

时间步长：$\Delta t$ 应该取多大？❹

一个时间步内，最快的物理信息传播距离不能超过一个网格。

也即：

\[ (a+V)\Delta t\leq \Delta x \]

本问题由非线性偏微分方程控制，因此：

严格的稳定性条件 $C\leq1$ 仅作为指导原则
但实践证明该条件仍然有效且可靠

显式格式要求 $\color{black}{\Delta t}$ 足够小，使得声波和流动输运共同造成的信息传播不越过一个网格。

MacCormack方法

局部时间步要求

由于：

$\Delta x$ 固定
$V,a$ 随空间变化

因此每个网格点对应的时间步长要求不同：

\[ (\Delta t)_i^t= C\frac{\Delta x}{a_i^t+V_i^t} \]

相邻点：

\[ (\Delta t)_{i+1}^t= C\frac{\Delta x}{a_{i+1}^t+V_{i+1}^t} \]

👉 各网格点的时间步长要求一般不同。

方案 1：局部时间步（Local Time Stepping）

每个网格点使用自己的时间步：

\[ (\Delta t)_i^t \]

每个点时间推进不同

解对应“非物理时间”

优点：收敛快（steady state）

缺点：不能描述真实瞬态过程、存在“时间扭曲”（time warp）

方案 2：全局时间步（推荐）

在所有网格点计算：

\[ (\Delta t)_1^t,\quad (\Delta t)_2^t,\quad \dots,\quad (\Delta t)_N^t \]

取最小值：

\[ \Delta t= \min\left[ (\Delta t)_1^t, (\Delta t)_2^t, \dots, (\Delta t)_N^t \right] \]

即：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \Delta t= \min_i \left( C\frac{\Delta x}{a_i+V_i} \right) } $} \]

所有点同步推进（对应同一物理时间）

优点：满足真实瞬态演化、数值解具有物理意义

缺点：收敛较慢

时间步长

思考

CFL 条件本质是信息传播约束：数值时间步不能超过物理信息传播所允许的尺度。

为什么分母是 $a+V$ 而不是 $V$？
因为可压缩流中最快信息速度是声波相对流体传播后的合速度，即 ${a+V}$。

在激波附近，$\Delta t$ 会如何变化？
激波附近速度、声速和梯度变化剧烈，通常需要更小的 $\Delta t$ 以保持稳定。

若采用隐式格式，CFL 条件是否仍然限制稳定性？
对线性稳定性的限制可放宽，但精度、非线性收敛和物理分辨率仍会限制 ${\Delta t}$。

边界条件（Boundary Conditions）

数值求解的另一个关键方面是边界条件。

如果没有在物理上正确设置边界条件，并在数值上正确实现它们，就无法得到正确的数值解。

网格点 1 和 N 是两个边界点

点 1：亚声入口边界（subsonic inflow boundary）

点 N：出口边界（outflow boundary）

边界条件（Boundary Conditions）

亚声入口边界

点 1 位于储气室附近，流体从储气室进入喷管。

入口流动为低速亚声流。

速度不能为 0，否则无质量流进入喷管。

因此，它不完全等同于理想储气室。

边界应给定什么？

特征线方法

对于一维非定常无粘流：方程为双曲型、存在三条特征线：

\[ \frac{dx}{dt}=V-a,\qquad \frac{dx}{dt}=V,\qquad \frac{dx}{dt}=V+a \]

物理意义：

$V-a$：向上游传播的声学扰动
$V$：随流体运动的熵波/接触波
$V+a$：向下游传播的声学扰动

边界条件的个数由进入计算域的特征线数目决定。

亚声入口边界（点 1）

入口满足：

\[ 0<V_1<a_1 \]

此时三条特征线：

$V_1-a_1$：离开计算域
$V_1$：进入计算域
$V_1+a_1$：进入计算域

👉 因此亚声入口有两条特征线进入计算域，需要指定两个边界变量： \[\rho_1=1,\qquad T_1=1\] 它们对应储气室总参数。速度 $V_1$ 由内部点决定。

\[ V_1=2V_2-V_3 \]

超声出口边界（点 N）

出口满足：

\[ V_N>a_N \]

此时三条特征速度为：

\[ V_N-a_N>0,\qquad V_N>0,\qquad V_N+a_N>0 \]

三条特征线均指向下游
在出口边界处全部离开计算域

👉 在超声出口，无需指定任何变量。

👉 所有变量都由内部解决定。

\[ V_N=2V_{N-1}-V_{N-2}, \qquad \rho_N=2\rho_{N-1}-\rho_{N-2}, \qquad T_N=2T_{N-1}-T_{N-2} \]

喷管几何与初始条件

喷管形状 $A=A(x)$ 为已知并固定，与时间无关。

本算例采用如下抛物线型面积分布：

\[ A=1+2.2(x-1.5)^2,\qquad 0\leq x\leq 3 \]

$x=1.5$：喉部

$x<1.5$：收缩段

$x>1.5$：扩张段

初始条件的设定

为了开始时间推进计算，需要给定：

\[ \rho(x),\qquad T(x),\qquad V(x) \quad \text{在 } t=0 \]

理论上：初始条件可以任意

实践中：必须合理选择

初始条件的设定

为什么要“聪明地”选择初始条件？

原因 1：初始条件越接近最终稳态，收敛越快，计算时间越短

原因 2：若初始条件偏离实际过远，初始时间导数会很大，早期计算可能产生巨大梯度，数值计算可能不稳定

可以将时间推进过程类比为：一根被拉伸的橡皮筋。

初始时：拉伸很大 → 强驱动力

随时间推进：逐渐放松 → 收敛减慢

👉 若拉得过度：数值解可能“断裂”（发散）

初始条件选取原则

应利用对问题的物理认知：

喷管中：

$\rho$、$T$ 递减
$V$ 递增

本例采用的初始条件：

\[ \rho=1-0.3146x \]

\[ T=1-0.2314x \]

\[ V=(0.1+1.09x)T^{1/2} \]

💻 代码讲解和演示 nozzle_nonconservative_transonic_supersonic.m

课堂互动：先预测，再计算

对于收缩—扩张喷管：

入口总压、总温固定
逐渐降低出口压力 $p_e$

请判断流动状态如何变化？

出口压力较高时，喷管内是否全场亚声速？
出口压力继续降低，喉部会发生什么？
再继续降低，扩张段是否可能出现激波？
如果出现激波，应优先使用守恒形式还是非守恒形式？

参考判断：

\[ \text{全亚声速}\rightarrow\text{喉部阻塞}\rightarrow\text{超声区/激波} \]

含激波问题中，守恒形式更可靠。

纯亚声速喷管流的 CFD 求解

PURELY SUBSONIC NOZZLE FLOW

纯亚声速喷管流

同一个收缩—扩张喷管，
由于两端压比不同，也可能形成全程亚声速的等熵流动。

Subsonic branch · Exit pressure · Boundary-condition sensitivity

纯亚声速喷管流的物理说明

在前面的算例中，我们研究了亚声速—超声速等熵喷管流。

本节讨论另一类不同的解：

纯亚声速（purely subsonic）等熵喷管流

其物理本质在于：虽然喷管仍然是收缩—扩张型，但由于喷管两端所施加的压比不同，流动在整个喷管内都保持亚声速。

因此，喉部也不会达到声速。

纯亚声速喷管流的物理说明 ❶

对于某一给定喷管，若出口压力足够高，则喷管中的流动虽然会在收缩段加速、在扩张段减速，但局部马赫数始终满足：

\[ M<1 \]

这就是纯亚声速喷管流。

在这种情况下，喉部不再是声速截面，因此需要区分两个面积：

$A_t$：喷管真实喉部面积，即几何最小面积
$A^*$：若该流动在某一截面达到 $M=1$ 时对应的临界面积

对于堵塞流动，喉部就是声速截面，因此 $A_t=A^*$。
但对于纯亚声速流动，喉部马赫数仍小于 1，因此二者不再相等。

纯亚声速喷管流的物理说明 ❷

等熵准一维喷管流满足面积—马赫数关系：

\[ \frac{A}{A^*} =\frac{1}{M} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M^2 \right) \right]^{\frac{\gamma+1}{2(\gamma-1)}} \]

对于亚声速流动，有：

\[ M<1 \]

代入面积—马赫数关系可知：

\[ \frac{A}{A^*}>1 \]

也就是说，任意一个亚声速截面的面积都大于其对应的临界面积。

纯亚声速喷管流的物理说明 ❸

特别地，在喷管喉部：$M_t<1$，因此：$\frac{A_t}{A^*}>1$，于是得到：

\[ \boxed{A^*<A_t} \]

物理含义是：
当前流量还不足以使喉部达到声速。若要把这股流动加速到 $M=1$，需要一个比真实喉部更小的“等效临界面积” $A^*$。

也就是说，$A^*$ 在这里仅仅是一个参考面积，它小于真实喉部面积。

一旦 $A^*$ 被确定，局部面积比 $A/A^*$ 仍可通过面积—马赫数关系确定局部马赫数：

\[ \left(\frac{A}{A^*}\right)^2 =\frac{1}{M^2} \left[ \frac{2}{\gamma+1} \left(1+\frac{\gamma-1}{2}M^2\right) \right]^{(\gamma+1)/(\gamma-1)} \]

由于这里要求的是纯亚声速解，因此求 $M$ 时应始终取亚声速分支。

纯亚声速喷管流的物理说明

一旦局部马赫数 $M$ 已知，局部压力比、密度比和温度比可继续由等熵关系得到。

压力比：

\[ \frac{p}{p_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-\gamma/(\gamma-1)} \]

密度比：

\[ \frac{\rho}{\rho_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-1/(\gamma-1)} \]

温度比：

\[ \frac{T}{T_0} =\left(1+\frac{\gamma-1}{2}M^2\right)^{-1} \]

纯亚声速喷管流的物理说明

因此，在纯亚声速喷管流中：

\[ \frac{A}{A^*}\rightarrow M\rightarrow \frac{p}{p_0},\frac{\rho}{\rho_0},\frac{T}{T_0} \]

几何关系仍然决定 $A/A^*$，而 $A/A^*$ 再决定 $M$，继而决定局部热力学状态。

与亚声速—超声速解的根本差别在于：

这里整个流场始终走的是亚声速分支。

数值求解设置

喷管几何

本节采用的喷管面积分布为：

\[ \frac{A}{A_t}=1+2.2(x-1.5)^2,\qquad 0\le x\le 1.5 \]

以及：

\[ \frac{A}{A_t}=1+0.2223(x-1.5)^2,\qquad 1.5\le x\le 3.0 \]

这里 $A_t$ 表示喷管喉部面积。

在纯亚声速情形下，喉部处并不满足 $M=1$，因此：

\[ A_t\ne A^*,\qquad A^*<A_t \]

控制方程

本节仍采用与前面亚声速—超声速喷管流相同的准一维非定常控制方程。

连续方程：

\[ \frac{\partial \rho}{\partial t} =-\rho\frac{\partial V}{\partial x} -\rho V\frac{\partial(\ln A)}{\partial x} -V\frac{\partial \rho}{\partial x} \]

动量方程：

\[ \frac{\partial V}{\partial t} =-V\frac{\partial V}{\partial x} -\frac{1}{\gamma} \left( \frac{\partial T}{\partial x} +\frac{T}{\rho}\frac{\partial \rho}{\partial x} \right) \]

能量方程：

\[ \frac{\partial T}{\partial t} =-V\frac{\partial T}{\partial x} -(\gamma-1)T \left[ \frac{\partial V}{\partial x} +V\frac{\partial(\ln A)}{\partial x} \right] \]

控制方程

从方程形式看，本节与前一算例并无本质区别。
真正改变纯亚声速解的是：

喷管几何 $A(x)$

初始流场分布

以及最关键的下游边界条件

边界条件的物理分析

入口边界：亚声速入口

入口仍然是亚声速流入边界，因此和前一算例一样，入口应允许一个变量“浮动”，其余变量由储气室条件给定。

这里取：

入口速度由内点外推
入口密度与温度固定为储气室值

即：

\[ \rho_1=1,\qquad T_1=1 \]

入口速度由线性外推给出：

\[ V_1=2V_2-V_3 \]

边界条件的物理分析

出口边界：亚声速出口

与前一算例不同，本节的出口改为亚声速出口。

因此在出口点 $N$，特征线分析表明：

有两个变量应允许“浮动”
但必须有一个变量被指定

从物理上说，这个被指定的量正是：

\[ p_N \]

即出口压力。

因为对于纯亚声速喷管流，要获得唯一解，必须给定跨喷管的压比，或者在固定 $p_0$ 的条件下给定出口静压 $p_e$。

出口压力边界的实现

数值方程中的未知量是：$\rho$、$V$、$T$，而不是 $p$。

因此需要借助无量纲状态方程：

\[ p=\rho T \]

于是，在出口点 $N$，压力边界条件可写为：

\[ p_N=\text{给定值} \]

由于控制变量是 $\rho_N$ 和 $T_N$，必须保证二者在每一时间步中始终满足状态方程：

\[ \rho_NT_N=p_N=\text{给定值} \]

出口边界条件的具体实现

方法一：先外推温度

先由内部点线性外推出口温度：

\[ T_N=2T_{N-1}-T_{N-2} \]

然后由状态方程求出口密度：

\[ \rho_N=\frac{p_N}{T_N} \]

这样便能保证出口压力满足给定值：

\[ \rho_NT_N=p_N \]

这里的关键不是单独指定 $\rho_N$ 或 $T_N$，而是通过 $p_N=\rho_NT_N$ 保证出口压力边界条件被严格满足。

出口边界条件的具体实现

方法二：先外推密度

也可以先由内部点线性外推出口密度：

\[ \rho_N=2\rho_{N-1}-\rho_{N-2} \]

再由状态方程求出口温度：

\[ T_N=\frac{p_N}{\rho_N} \]

两种做法在数值实现上都可以使用，本质上都是为了满足：

\[ p_N=\rho_NT_N \]

此外，出口速度通常由线性外推给出：

\[ V_N=2V_{N-1}-V_{N-2} \]

初始条件

初始流场采取如下较为简单的无量纲分布：

\[ \rho=1.0-0.023x \]

\[ T=1.0-0.009333x \]

\[ V=0.05+0.11x \]

这些初始条件给出了 $t=0$ 时刻的流场分布。

需要注意：这些初始条件只是计算的起点，并不要求严格满足最终的定常解。随着时间推进，流场会在控制方程和边界条件共同作用下逐渐收敛。

数值方法

对该纯亚声速喷管流，仍使用与前一算例完全相同的：

MacCormack 显式预测—校正格式

因此，在程序实现上，只需修改：

喷管形状

初始条件

下游边界条件

而不必改动差分算法本身。

采用压比为：$\frac{p_e}{p_0}=0.93$的纯亚声速算例。

💻 代码讲解和演示 nozzle_nonconservative_subsonic_pressure_exit.m

本节小结

纯亚声速喷管流与亚声速—超声速喷管流相比，最大的差异不在于控制方程或数值格式，而在于：

物理解属于纯亚声速分支

出口是亚声速出口，必须指定出口压力

出口压力边界的数值实现会深刻影响时间推进的稳定性

即使最终稳态边界条件是物理正确的，它在非定常时间推进过程中也不一定总是数值上合适。

这正是数值边界条件设计中最微妙、也最重要的问题之一。

准一维喷管流的守恒形式与非守恒形式比较

CONSERVATIVE VS. NONCONSERVATIVE FORM

守恒形式与非守恒形式

同一个物理问题，
采用不同的方程形式，可能得到完全不同的激波结果。

Conservation law · Shock capturing · Physical consistency

为什么要比较两种形式？

前面的准一维喷管流求解中，我们主要采用的是：

非守恒形式（nonconservative form）

其基本变量为：

$\rho$
$V$
$T$

这种形式物理直观，便于理解变量随时间和空间的演化。

另一种重要表达：守恒形式

在 CFD 中，还有另一种更重要的表达方式：

守恒形式（conservative form）

其基本变量为：

$\rho A$
$\rho V A$
$\rho E A$

它直接对应控制体中的质量、动量和能量守恒。

本节目标

本节的目的在于：

推导喷管流的守恒形式方程

使用 MacCormack 方法进行求解

与非守恒形式结果进行比较

揭示一个关键 CFD 结论

守恒形式在处理激波问题时具有决定性优势。

守恒形式控制方程推导

连续方程

准一维连续方程的守恒形式为：

\[ \frac{\partial(\rho A)}{\partial t} +\frac{\partial(\rho V A)}{\partial x}=0 \]

该方程表示：

控制体内质量的变化率 + 通过控制体边界的质量通量 = 0

其中 $\rho A$ 是单位长度内的质量，$\rho V A$ 是质量流率。

守恒形式控制方程推导

动量方程 ❶

从前面的准一维动量方程出发：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial(\rho V^2A)}{\partial x} =-A\frac{\partial p}{\partial x} \]

为了写成守恒形式，需要将右侧压力梯度项重新整理。

利用乘积求导：

\[ \frac{\partial(pA)}{\partial x} =A\frac{\partial p}{\partial x} +p\frac{\partial A}{\partial x} \]

因此：

\[ -A\frac{\partial p}{\partial x} =-\frac{\partial(pA)}{\partial x} +p\frac{\partial A}{\partial x} \]

守恒形式控制方程推导

动量方程 ❷

将压力项整理结果代入动量方程：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial(\rho V^2A)}{\partial x} =-\frac{\partial(pA)}{\partial x} +p\frac{\partial A}{\partial x} \]

将通量项合并到左侧：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial(\rho V^2A+pA)}{\partial x} =p\frac{\partial A}{\partial x} \]

也即：

\[ \frac{\partial(\rho V A)}{\partial t} +\frac{\partial\left[(\rho V^2+p)A\right]}{\partial x} =p\frac{dA}{dx} \]

右侧源项来自喷管面积变化。

守恒形式控制方程推导

能量方程

总能量定义为：

\[ E=e+\frac{V^2}{2} \]

准一维喷管流的能量守恒形式为：

\[ \frac{\partial(\rho E A)}{\partial t} +\frac{\partial(\rho V E A+pVA)}{\partial x}=0 \]

也可以写为：

\[ \frac{\partial(\rho E A)}{\partial t} +\frac{\partial\left[\rho V A\left(E+\frac{p}{\rho}\right)\right]}{\partial x}=0 \]

该方程表示总能量随流动输运，同时包含压力功的贡献。

守恒变量定义

定义守恒变量：

\[ U_1=\rho A \]

\[ U_2=\rho V A \]

\[ U_3=\rho E A \]

对应地，$\mathbf{U}$ 表示控制体内的守恒量。

通量与源项定义

定义通量：

\[ F_1=\rho V A \]

\[ F_2=\rho V^2A+pA \]

\[ F_3=\rho VEA+pVA \]

定义源项：

\[ \mathbf{J}= \begin{bmatrix} 0\\ p\dfrac{dA}{dx}\\ 0 \end{bmatrix} \]

统一向量形式

于是，守恒形式可统一写为：

\[ \frac{\partial\mathbf{U}}{\partial t} +\frac{\partial\mathbf{F}}{\partial x} =\mathbf{J} \]

这就是后续 MacCormack 方法直接离散的方程形式。

其中：

$\mathbf{U}$：守恒变量
$\mathbf{F}$：通量
$\mathbf{J}$：几何源项

MacCormack 方法：守恒形式

对守恒形式：

\[ \frac{\partial\mathbf{U}}{\partial t} +\frac{\partial\mathbf{F}}{\partial x} =\mathbf{J} \]

仍采用显式预测—校正思想。

不同之处在于：推进对象由原始变量变为守恒变量 $\mathbf{U}$。

MacCormack 方法：预测步

预测步采用前向差分：

\[ \overline{U}_i =U_i^n-\frac{\Delta t}{\Delta x}\left(F_{i+1}^n-F_i^n\right)+J_i^n\Delta t \]

其中：

$\overline{U}_i$：预测后的守恒变量
$F_{i+1}^n-F_i^n$：前向通量差
$J_i^n$：当前时刻源项

这是守恒形式下的预测更新。

MacCormack 方法：校正步

校正步采用后向差分：

\[ U_i^{n+1} =\frac{1}{2} \left[ U_i^n+\overline{U}_i -\frac{\Delta t}{\Delta x}\left(\overline{F}_i-\overline{F}_{i-1}\right) +\overline{J}_i\Delta t \right] \]

其中：

$\overline{F}_i-\overline{F}_{i-1}$：预测状态下的后向通量差
$\overline{J}_i$：预测状态下的源项

预测步与校正步平均后，得到二阶精度。

状态变量恢复

MacCormack 方法推进的是守恒变量：

\[ U_1=\rho A,\qquad U_2=\rho V A,\qquad U_3=\rho E A \]

但为了计算通量和物理量，需要恢复原始变量。

由守恒变量可得：

\[ \rho=\frac{U_1}{A}, \qquad V=\frac{U_2}{U_1}, \qquad E=\frac{U_3}{U_1} \]

\[ T=e=(\gamma-1)\left(\frac{U_3}{U_1}-\frac{\gamma}{2}V^2\right) \]

\[ p=\rho T \]

数值结果比较

现在比较两种方法：

1. 非守恒形式

基本变量：

\[ \rho,\quad V,\quad T \]

2. 守恒形式

基本变量：

\[ \rho A,\quad \rho V A,\quad \rho E A \]

无激波情况下

在无激波、解足够光滑的情况下：

两种形式的结果基本一致。

原因是：对于光滑解，守恒形式与非守恒形式在连续意义下可以相互转换。

存在激波情况下

在存在激波时，情况发生根本变化。

两种形式会出现本质差异。

非守恒形式的问题：

无法正确捕捉激波位置
可能出现错误解

守恒形式的优势：

满足守恒律
正确处理不连续
激波位置准确

物理本质

为什么守恒形式更重要？

激波本质上是守恒律的弱解。

经典解要求“每一点都满足微分方程”；弱解只要求“任意控制体内的守恒关系成立”。激波处物理量不连续，不能作为经典解，但只要满足跨激波的守恒跳跃条件，它就是守恒律允许的弱解。

激波处变量不连续，不能再依赖普通微分意义下的方程变形。

此时真正可靠的是跨越不连续面的积分守恒关系。

守恒形式的意义

守恒形式保证：

\[ \text{质量、动量、能量守恒} \]

非守恒形式在光滑区域通常没有问题。

但在激波处，非守恒形式可能破坏守恒关系。

因此，含激波问题中，方程是否写成守恒形式不是形式问题，而是物理正确性问题。

涉及不连续，特别是激波的问题，必须使用守恒形式。

方法	是否守恒	激波能力
非守恒形式	❌	差
守恒形式	✅	强

激波捕捉：人工黏性方法

ARTIFICIAL VISCOSITY

激波捕捉与人工黏性

守恒形式对于正确处理激波是必要的，
但并不充分。

激波附近仍可能出现非物理数值振荡

激波捕捉算例

准一维喷管内激波计算 ❶

本算例研究收缩—扩张喷管中的准一维可压缩流动：

\[ A=A(x),\qquad 0\le x\le3 \]

喷管面积分布为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ A(x)=1+2.2(x-1.5)^2 } $} \]

出口压力指定为：

\[ p_e=0.6784 \]

扩张段内部将出现一道正激波。

这类问题的关键不是预先指定激波位置，而是让数值格式在守恒方程、边界条件和适当数值耗散的共同作用下自动捕捉激波。

激波捕捉算例

内激波流动的基本图像 ❷

当喷管内部存在正激波时，典型流动过程为：

收缩段：流动加速，压力下降；
喉部：达到声速，$M=1$；
扩张段前部：继续加速为超声速；
正激波处：马赫数由超声速突降为亚声速，压力突增；
激波后：亚声速流动在扩张通道中减速，压力继续回升。

激波位置由喷管几何、入口总量、出口背压以及守恒方程共同决定（手动推导）。

激波捕捉算例

喷管几何与喉部位置 ❸

计算域和喷管面积为：

\[ x\in[0,3],\qquad A(x)=1+2.2(x-1.5)^2 \]

喉部位置由面积最小值确定：

\[ A_t=\min A(x),\qquad x_t=\operatorname*{arg\,min}_x A(x) \]

代码中对应：

A = 1 + 2.2*(x-1.5).^2;
dAdx = gradient(A,dx);

[At,ithroat] = min(A);
x_throat = x(ithroat);

激波捕捉算例

先算解析解 ❹

本算例中，出口压力被指定为：$p_e=0.6784$，这意味着喷管不能保持完整的亚声—超声等熵解。

当出口背压高于设计等熵出口压力、但又低于临界压力时，流动会出现：

喉部阻塞：$M_t=1$
扩张段前部：超声速等熵加速
扩张段内部：形成正激波
激波后：亚声速等熵减速并升压

解析解的目的不是替代 CFD，而是给出一个清楚的参照：
激波应该出现在什么位置？激波前后状态大约是多少？

激波捕捉算例

出口压力与激波位置的关系 ❺

对于完整的亚声—超声等熵解，出口压力不是任意指定的，而是由喷管面积比自动决定。

但对于含激波的喷管流，出口背压必须给定：

\[ \frac{p_e}{p_{01}}=0.6784 \]

这个背压决定了激波位置。

激波越靠下游，激波前马赫数越大，激波越强；
激波越靠上游，激波前马赫数越小，压力恢复过程也随之改变。

\[ p_e \Longrightarrow \text{激波强度} \Longrightarrow \text{激波位置} \]

\[ \text{几何}+p_e \Longrightarrow A_1 \]

激波捕捉算例

激波前后区域的记号 ❻

设正激波位于喷管扩张段某一截面，其面积为：

\[ A_1 \]

激波前状态记为 1，激波后状态记为 2：

\[ M_1,\ p_{01},\ A_1^* \qquad\Longrightarrow\qquad M_2,\ p_{02},\ A_2^* \]

激波前：从储气室到激波前为等熵流动，因此总压为 $p_{01}$。
激波后：从激波后到出口仍近似等熵，但总压降低为 $p_{02}$。

\[ p_{02}<p_{01} \]

\[ s_2>s_1 \]

\[ A_2^*>A_1^* \]

激波捕捉算例

质量流量约束 ❼

准一维等熵临界流量可写为：

\[ \dot{m}= \frac{p_0A^*}{\sqrt{T_0}} \sqrt{ \frac{\gamma}{R} \left( \frac{2}{\gamma+1} \right)^{(\gamma+1)/(\gamma-1)} } \]

对于给定入口总温条件，有：

\[ \dot{m}\propto p_0A^* \]

由于质量流量穿过激波前后保持不变：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ p_{01}A_1^*=p_{02}A_2^* } $} \]

激波会造成总压损失，因此 $p_{02}<p_{01}$。
为了保持相同质量流量，激波后的等效临界面积必须增大，即 $A_2^*>A_1^*$。

激波捕捉算例

出口面积比与出口马赫数 ❽

对于本喷管，出口面积与喉部面积之比为：

\[ \frac{A_e}{A_t}=5.95 \]

激波前流动已在喉部阻塞，因此：

\[ A_1^*=A_t \]

将出口压力和面积比组合起来：

\[ \frac{p_eA_e}{p_{01}A_t}= 0.6784\times5.95= 4.03648 \]

激波后到出口为亚声速等熵流动，于是：

\[ \frac{p_eA_e}{p_{02}A_2^*}= \frac{p_eA_e}{p_{01}A_t}= 4.03648 \]

激波捕捉算例

出口马赫数的求解 ❾

出口处等熵关系给出：

\[ \frac{A_e}{A_2^*}= \frac{1}{M_e} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M_e^2 \right) \right]^{\frac{\gamma+1}{2(\gamma-1)}} \]

以及：

\[ \frac{p_e}{p_{02}}= \left( 1+\frac{\gamma-1}{2}M_e^2 \right)^{-\gamma/(\gamma-1)} \]

联立后得到关于 $M_e$ 的方程：

\[ \frac{1}{M_e} \left( \frac{2}{\gamma+1} \right)^{\frac{\gamma+1}{2(\gamma-1)}} \left[ 1+\frac{\gamma-1}{2}M_e^2 \right]^{-1/2}= 4.03648\Longrightarrow \colorbox{yellow}{$\displaystyle \color{black}{ M_e=0.1431 } $} \]

激波捕捉算例

激波后的总压损失 ❿

由出口等熵关系：

\[ \frac{p_e}{p_{02}}= \left[ 1+\frac{\gamma-1}{2}(0.1431)^2 \right]^{-3.5}= 0.9858 \]

因此，激波前后总压比为：

\[ \frac{p_{02}}{p_{01}}= \frac{p_{02}}{p_e} \frac{p_e}{p_{01}}= \frac{0.6784}{0.9858}= 0.6882 \]

这说明正激波造成了明显的总压损失。
出口静压虽然由边界条件指定，但出口总压已经低于入口总压。

激波捕捉算例

由总压比确定激波前马赫数 ⓫

正激波总压比只由激波前马赫数 $M_1$ 决定：

\[ \frac{p_{02}}{p_{01}}= \left[ \frac{(\gamma+1)M_1^2}{(\gamma-1)M_1^2+2} \right]^{\gamma/(\gamma-1)} \left[ \frac{\gamma+1}{2\gamma M_1^2-(\gamma-1)} \right]^{1/(\gamma-1)} \]

令：

\[ \frac{p_{02}}{p_{01}}=0.6882 \]

求解得到激波前马赫数：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ M_1=2.07 } $} \]

也就是说，数值计算若正确捕捉激波，激波前的局部马赫数应接近 $2.07$。

激波捕捉算例

由激波前马赫数确定激波位置 ⓬

激波前为等熵超声速流动，因此局部面积比满足：

\[ \frac{A_1}{A_1^*}= \frac{1}{M_1} \left[ \frac{2}{\gamma+1} \left( 1+\frac{\gamma-1}{2}M_1^2 \right) \right]^{\frac{\gamma+1}{2(\gamma-1)}} \]

由于 $A_1^*=A_t$，代入 $M_1=2.07$ 得：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{A_1}{A_t}=1.790 } $} \]

因此，解析解预测：
正激波应位于喷管扩张段中面积比约为 $A/A_t=1.79$ 的位置。

激波捕捉算例

对应的几何位置 ⓭

本算例喷管面积为：

\[ A(x)=1+2.2(x-1.5)^2 \]

喉部面积为：

\[ A_t=1 \]

由解析解：

\[ \frac{A_1}{A_t}=1.790 \Longrightarrow 1+2.2(x_s-1.5)^2=1.790 \]

激波位于扩张段，故取 $x_s>1.5$：

\[ x_s=1.5+\sqrt{\frac{0.790}{2.2}}\Longrightarrow \colorbox{yellow}{$\displaystyle \color{black}{ x_s\simeq2.10 } $} \]

激波捕捉算例

激波前后状态 ⓮

由正激波关系，激波前后压力比为：

\[ \frac{p_2}{p_1}= 1+\frac{2\gamma}{\gamma+1}(M_1^2-1) \]

代入 $M_1=2.07$，得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{p_2}{p_1}=4.83 } $} \]

激波后的马赫数为：

\[ M_2= \left[ \frac{1+\frac{\gamma-1}{2}M_1^2} {\gamma M_1^2-\frac{\gamma-1}{2}} \right]^{1/2}= \colorbox{yellow}{$\displaystyle \color{black}{ 0.566 } $} \]

激波捕捉算例

解析解汇总 ⓯

对于指定出口压力：

\[ \frac{p_e}{p_{01}}=0.6784 \]

解析结果为：

出口马赫数：$M_e=0.1431$
激波前后总压比：$\dfrac{p_{02}}{p_{01}}=0.6882$；静压跃升：$\dfrac{p_2}{p_1}=4.83$
激波前后马赫数：$M_1=2.07$，$M_2=0.566$
激波位置面积比：$\dfrac{A_1}{A_t}=1.790$
激波位置：$x_s\simeq2.10$

后续数值结果应重点检查：
激波位置是否接近 $x\simeq2.10$，以及激波前后马赫数和压力跃升是否合理。

激波捕捉算例

守恒变量 ⓰

代码采用无量纲准一维守恒变量：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ U_1=\rho A } $} \]

\[ \colorbox{yellow}{$\displaystyle \color{black}{ U_2=\rho V A } $} \]

\[ \colorbox{yellow}{$\displaystyle \color{black}{ U_3=\rho A\left( \frac{T}{\gamma-1} +\frac{\gamma}{2}V^2 \right) } $} \]

其中：

\[ p=\rho T,\qquad M=\frac{V}{\sqrt{T}} \]

激波捕捉算例

守恒方程形式 ⓱

准一维喷管流守恒方程写为：

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

其中通量项与几何源项分别为：

\[ \mathbf{F}= \begin{bmatrix} F_1\\[3pt] F_2\\[3pt] F_3 \end{bmatrix}= \begin{bmatrix} U_2\\[6pt] \dfrac{U_2^2}{U_1}+\dfrac{\gamma-1}{\gamma} \left( U_3-\dfrac{\gamma}{2}\dfrac{U_2^2}{U_1} \right)\\[8pt] \gamma\dfrac{U_2U_3}{U_1} -\dfrac{\gamma(\gamma-1)}{2}\dfrac{U_2^3}{U_1^2} \end{bmatrix}, \qquad \mathbf{J}= \begin{bmatrix} 0\\[6pt] \dfrac{1}{\gamma}p\dfrac{dA}{dx}\\[6pt] 0 \end{bmatrix}. \]

激波捕捉算例

为什么必须使用守恒形式？ ⓲

守恒形式

离散的是通量差：

\[ F_{i+1/2}-F_{i-1/2} \]

可以保证质量、动量和能量在控制体意义下守恒。

非守恒形式

直接离散局部导数：

\[ A(q)q_x \]

在激波间断处没有唯一含义，可能得到错误激波位置和跳跃强度。

激波捕捉的核心要求是：数值解应收敛到满足 Rankine–Hugoniot 跳跃关系的弱解。

激波捕捉算例

MacCormack格式：预测步 ⓳

预测步采用前向差分：

\[ \bar{U}_i= U_i^n -\frac{\Delta t}{\Delta x} \left( F_{i+1}^n-F_i^n \right) +\Delta tJ_i^n +S_i^n \]

代码中对应：

for i = 2:Nx-1
    U_bar(:,i) = U(:,i) ...
        - dt/dx*(F(:,i+1)-F(:,i)) ...
        + dt*J(:,i);
end

S = artificial_viscosity(U,p,Cx);

for i = 2:Nx-1
    U_bar(:,i) = U_bar(:,i) + S(:,i);
end

预测步先利用当前时刻的流场向前推进，得到中间状态 $\bar{U}$。

激波捕捉算例

MacCormack格式：校正步 ⓴

校正步采用后向差分：

\[ U_i^{n+1}= \frac{1}{2} \left[ U_i^n+\bar{U}_i -\frac{\Delta t}{\Delta x} \left( \bar{F}_i-\bar{F}_{i-1} \right) +\Delta t\bar{J}_i \right] +\bar{S}_i \]

代码中对应：

for i = 2:Nx-1
    U_new(:,i) = 0.5*( ...
        U(:,i) + U_bar(:,i) ...
        - dt/dx*(F_bar(:,i)-F_bar(:,i-1)) ...
        + dt*J_bar(:,i) );
end

S_bar = artificial_viscosity(U_bar,p_b,Cx);

for i = 2:Nx-1
    U_new(:,i) = U_new(:,i) + S_bar(:,i);
end

MacCormack格式通过“预测+校正”获得二阶时间与空间精度，但在激波附近仍需要额外数值耗散。

激波捕捉算例

CFL时间步㉑

每一步时间步长由局部传播速度限制：

\[ \Delta t= \mathrm{CFL} \frac{\Delta x} {\max(|V|+a)} \]

其中无量纲声速为：

\[ a=\sqrt{T} \]

代码中对应：

a = sqrt(T);
dt = CFL*dx/max(abs(V)+a);

CFL条件反映了显式格式的时间步必须受物理信息传播速度限制。

激波捕捉算例

边界条件：入口㉒

入口采用简化的亚声速入口条件：

固定入口总量近似：

\[ \rho_{\mathrm{in}}=1,\qquad T_{\mathrm{in}}=1 \]

速度由内部外推：

\[ V_1=2V_2-V_3 \]

再由守恒变量重构入口状态：

\[ U_1=\rho A \]

\[ U_2=\rho VA \]

\[ U_3= \rho A \left( \frac{T}{\gamma-1} +\frac{\gamma}{2}V^2 \right) \]

激波捕捉算例

边界条件：出口㉓

出口为亚声速出口，指定静压：

\[ p_N=p_e \]

先外推质量变量和动量变量：

\[ (U_1)_N=2(U_1)_{N-1}-(U_1)_{N-2} \]

\[ (U_2)_N=2(U_2)_{N-1}-(U_2)_{N-2} \]

然后计算出口密度和速度：

\[ \rho_N=\frac{(U_1)_N}{A_N},\qquad V_N=\frac{(U_2)_N}{(U_1)_N} \]

由给定出口压力确定温度：

\[ T_N=\frac{p_e}{\rho_N} \]

最后更新能量变量：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ (U_3)_N= \rho_N A_N \left( \frac{T_N}{\gamma-1} +\frac{\gamma}{2}V_N^2 \right) } $} \]

亚声速出口不能简单外推所有变量。出口压力通过能量变量的重构进入计算。

激波捕捉算例

人工粘性：压力传感器㉔

对含激波流动，如果不加入人工粘性，数值解容易在激波附近出现非物理振荡，甚至导致残差不收敛。

压力二阶差分可以作为激波传感器：

\[ \phi_i= \frac{ \left|p_{i+1}-2p_i+p_{i-1}\right| }{ p_{i+1}+2p_i+p_{i-1} } \]

人工粘性项为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ S_i=C_x\phi_i \left( U_{i+1}-2U_i+U_{i-1} \right) } $} \]

代码中对应：

sensor = abs(p(i+1)-2*p(i)+p(i-1))/denom;

S(:,i) = Cx*sensor*(U(:,i+1)-2*U(:,i)+U(:,i-1));

压力二阶差分越大，说明局部压强变化越剧烈，人工粘性越强。

激波捕捉算例

人工粘性的作用㉕

好处

抑制激波附近振荡；
改善时间推进稳定性；
使压力和马赫数分布更平滑；
有助于获得定常解。

代价

激波会被抹宽；
梯度分辨率下降；
激波跳跃不再无限尖锐；
过强时会影响定量精度。

人工粘性的关键不是越大越好，而是在“稳定性”和“分辨率”之间取得平衡。

激波捕捉算例

数值结果应观察什么？㉖

在喷管激波捕捉算例中，建议重点观察：

压力分布 $p/p_0$：是否出现合理压力突升；
马赫数分布 $M$：是否从超声速跳到亚声速；
质量流量 $\dot{m}=\rho VA$：是否接近常数；
激波位置：是否与背压条件相匹配；
残差曲线：是否能够逐渐收敛；
人工粘性系数变化对解的影响。

激波捕捉算例

为什么数值解和解析解不会完全重合？㉗

解析解

激波是零厚度间断；
激波前后状态突跳；
满足精确正激波关系；
没有数值耗散影响。

数值解

激波被抹宽到若干网格；
人工粘性抑制振荡；
激波附近存在离散误差；
远离激波应与解析解接近。

因此，判断激波捕捉是否合理，不能只看激波点处是否完全重合，还要看上下游状态、激波位置和整体守恒性是否正确。

激波捕捉算例

本节小结㉘

喷管内正激波算例集中体现了 CFD 中几个核心思想：

时间推进求定常解；
守恒形式处理激波；
出口背压决定激波位置；
人工粘性稳定激波捕捉；
数值耗散与分辨率之间需要平衡。

一句话理解：

激波捕捉不是预先告诉程序激波在哪里，而是通过守恒方程、边界条件和适当数值耗散，让激波自然出现在正确位置附近。

激波捕捉算例

为什么非守恒形式不适合处理激波间断？ ❶

在光滑流动中，守恒形式和非守恒形式往往可以相互转化。

例如：

\[ \frac{\partial U}{\partial t}+\frac{\partial F(U)}{\partial x}=0 \]

可在光滑区域写成：

\[ \frac{\partial q}{\partial t}+A(q)\frac{\partial q}{\partial x}=0 \]

但激波不是光滑结构，而是变量发生跳跃的间断结构。
此时，非守恒形式与守恒形式不再等价。

激波捕捉算例

激波本质上是守恒律的弱解问题 ❷

激波前后，密度、速度、压力、温度都可能发生突变：

\[ \rho_1,u_1,p_1,T_1 \quad \longrightarrow \quad \rho_2,u_2,p_2,T_2 \]

但是，跨过激波的质量、动量和能量仍必须满足守恒关系。

也就是说，变量可以跳变，但控制体意义下的守恒不能被破坏。

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

考虑一维守恒律：

\[ \frac{\partial U}{\partial t}+\frac{\partial F(U)}{\partial x}=0 \]

设激波位置为：

\[ x=x_s(t) \]

激波速度定义为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ s=\frac{dx_s}{dt} } $} \]

激波两侧状态分别记为：

\[ U_1=U(x_s^-,t),\qquad U_2=U(x_s^+,t) \]

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

激波两侧状态分别记为：

\[ U_1=U(x_s^-,t),\qquad U_2=U(x_s^+,t) \]

对应通量为：

\[ F_1=F(U_1),\qquad F_2=F(U_2) \]

这里的 $s$ 不是声速，而是激波面在物理空间中的传播速度。
对于定常喷管内激波，激波位置不随时间移动，因此 $s=0$。

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

在包围激波的微小控制体上积分：

\[ \int_{x_s^-}^{x_s^+} \left( \frac{\partial U}{\partial t} +\frac{\partial F}{\partial x} \right)dx=0 \]

由于控制体边界随激波运动，需要考虑边界移动带来的通量贡献。
对激波两侧做积分平衡，可得：

\[ \frac{d}{dt} \int_{x_s^-}^{x_s^+}U\,dx+ F_2-F_1=0 \]

当控制体厚度趋于零时，积分项只保留激波移动扫过间断所造成的贡献（why）：

\[ \frac{d}{dt} \int_{x_s^-}^{x_s^+}U\,dx= -s\left(U_2-U_1\right) \]

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

设激波位置为 $x=x_s(t)$，速度为：

\[ s=\frac{dx_s}{dt} \]

在激波附近取一个随激波一起运动的微小控制体：

\[ [x_s(t)-\varepsilon,\;x_s(t)+\varepsilon] \]

积分量为：

\[ I(t)=\int_{x_s(t)-\varepsilon}^{x_s(t)+\varepsilon}U(x,t)\,dx \]

当 $\varepsilon$ 很小时，可近似认为：

\[ U(x,t)= \begin{cases} U_1, & x<x_s(t)\\ U_2, & x>x_s(t) \end{cases} \]

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

关键点：
积分区间内真正发生变化的，不是 $U_1$ 和 $U_2$ 本身，而是激波移动后，左状态区域和右状态区域的长度发生了变化。

假设在很短时间 $\Delta t$ 内，激波向右移动：

\[ \Delta x_s=s\Delta t \]

原来属于右状态 $U_2$ 的一小段区域，被左状态 $U_1$ 取代，长度为 $s\Delta t$。

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

因此积分量的变化为：

\[ \Delta I= U_1(s\Delta t)-U_2(s\Delta t)=-s\left(U_2-U_1\right)\Delta t \]

两边除以 $\Delta t$，并令 $\Delta t\to0$，得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{d}{dt} \int_{x_s^-}^{x_s^+}U\,dx= -s\left(U_2-U_1\right) } $} \]

也就是：

\[ \frac{dI}{dt}=-s[U] \]

激波捕捉算例

Rankine–Hugoniot 跳跃关系 ❸

因此：

\[ -s\left(U_2-U_1\right)+\left(F_2-F_1\right)=0 \]

即：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ s[U]=[F] } $} \]

其中：

\[ [U]=U_2-U_1,\qquad [F]=F_2-F_1 \]

Rankine–Hugoniot 关系说明：激波前后状态不能任意跳变，必须满足跨间断的质量、动量和能量守恒。

激波捕捉算例

守恒形式为什么可靠？ ❹

守恒形式计算的是通量差：

\[ \frac{\partial U}{\partial t} =-\frac{\partial F}{\partial x} \]

离散后对应：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ U_i^{n+1}= U_i^n -\frac{\Delta t}{\Delta x} \left( F_{i+1/2}-F_{i-1/2} \right) } $} \]

一个单元失去的质量、动量和能量，会作为通量进入相邻单元。
因此，守恒形式天然具有控制体意义上的守恒性。

激波捕捉算例

非守恒形式的问题在哪里？ ❺

非守恒形式通常包含如下项：

\[ A(q)\frac{\partial q}{\partial x} \]

在光滑流场中，这没有问题。

但在激波处，$q$ 发生跳跃，$\partial q/\partial x$ 不再是普通函数，而是类似 Dirac delta 的奇异分布。

此时 $A(q)$ 在间断处到底取左值、右值，还是某种平均值，并没有唯一规定。

激波捕捉算例

间断处的非守恒乘积没有唯一意义 ❻

在激波处：

\[ q= \begin{cases} q_L, & x<x_s\\ q_R, & x>x_s \end{cases} \]

于是：

\[ \frac{\partial q}{\partial x} \]

在普通函数意义下并不存在。

因此：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ A(q)\frac{\partial q}{\partial x} } $} \]

在间断处并没有天然唯一的数学含义。

这就是非守恒形式处理激波时最根本的问题。

激波捕捉算例

一个简单例子：Burgers 方程 ❼

守恒形式为：

\[ u_t+\left(\frac{u^2}{2}\right)_x=0 \]

在光滑区域，由链式法则可写成：

\[ u_t+uu_x=0 \]

在光滑区：

\[ \left(\frac{u^2}{2}\right)_x=uu_x \]

所以两种形式等价。

激波捕捉算例

但有间断时，二者不再等价 ❽

假设间断左、右状态分别为：

\[ u_L,\qquad u_R \]

守恒形式给出的激波速度由 Rankine–Hugoniot 条件决定：

\[ s= \frac{f(u_L)-f(u_R)}{u_L-u_R} \]

对于 Burgers 方程：

\[ f(u)=\frac{u^2}{2} \]

因此：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ s=\frac{u_L+u_R}{2} } $} \]

激波捕捉算例

非守恒形式可能给出错误的激波速度 ❾

如果直接使用非守恒形式：

\[ u_t+uu_x=0 \]

在间断附近，$u$ 应该取什么值并不明确。

取左值

\[ u=u_L \]

激波速度偏向左状态

取右值

\[ u=u_R \]

激波速度偏向右状态

取平均

可能依赖具体离散方法

因此，非守恒形式可能得到错误的激波传播速度和错误的激波位置。

激波捕捉算例

对喷管内激波意味着什么？ ❿

对于喷管内正激波，激波前后必须满足：

质量守恒；
动量守恒；
能量守恒；
正确的压力跳跃；
正确的马赫数跳跃；
正确的总压损失。

若采用非守恒形式，激波附近的跳跃关系可能被破坏。

激波捕捉算例

非守恒形式可能导致的问题 ⓫

激波位置错误

即使网格加密，数值激波也可能不收敛到正确位置。

激波强度错误

压力、密度、温度的跳跃幅值可能不满足 Rankine–Hugoniot 关系。

质量流量不守恒

喷管中 $\rho u A$ 可能在激波附近出现不合理变化。

收敛到错误弱解

计算也许稳定，但得到的不是物理正确的解。

激波捕捉算例

守恒形式 vs 非守恒形式 ⓬

对比项	守恒形式	非守恒形式
基本形式	$\dfrac{\partial U}{\partial t}+\dfrac{\partial F}{\partial x}=0$	$\dfrac{\partial q}{\partial t}+A(q)\dfrac{\partial q}{\partial x}=0$
光滑流动	适用	适用
激波间断	适用	不可靠
跳跃关系	自动满足 Rankine–Hugoniot 条件	不一定满足
激波位置	通常更可靠	可能偏移
物理守恒性	控制体意义下守恒	可能破坏守恒

激波捕捉算例

为什么守恒形式适合激波捕捉？ ⓭

守恒形式的优势在于：

基于控制体积分守恒；
能自然处理变量跳跃；
激波速度由通量跳跃决定；
可与人工粘性、通量分裂、Riemann 求解器结合；
更容易收敛到物理正确的弱解。

一句话理解：

激波不是普通的陡梯度，而是守恒律意义下的间断。

因此，激波捕捉必须尊重守恒形式。

激波捕捉算例

小结 ⓮

非守恒形式在光滑区域中可以使用，但在激波间断处会失去明确的数学意义。

真正决定激波前后状态的不是局部微分形式，而是跨间断的积分守恒关系。

因此，对于喷管内激波、激波管、高速外流、爆轰波等含间断问题，应优先采用守恒形式，并配合适当数值耗散或高分辨率激波捕捉格式。

小结与方法评价

SUMMARY AND PERSPECTIVE

准一维喷管流动算例总结

通过一个经典问题，
系统展示 CFD 方法的完整应用流程。

本章主要内容回顾

我们研究了两类喷管流动。

亚声速—超声速流

喉部达到 $M=1$
扩张段进入超声
不需要指定出口压力

纯亚声速流

全程 $M<1$
喉部不再是声速点
必须指定出口压力

控制方程形式

两种形式：

非守恒形式变量：

\[ \rho,\quad V,\quad T \]

特点：推导简单、物理直观。

守恒形式变量：

\[ \mathbf{U}= \begin{bmatrix} \rho A\\ \rho V A\\ \rho E A \end{bmatrix} \]

特点：满足守恒律，可处理不连续。

数值方法

采用：

MacCormack 显式预测—校正方法

主要特点：

二阶精度
显式格式
易实现
适合作为 CFD 入门算例

关键数值技术

时间推进求稳态：

\[ \frac{\partial}{\partial t}\rightarrow 0 \]

CFL 条件：

\[ \Delta t\leq\frac{\Delta x}{a+V} \]

边界条件：

亚声速入口：部分变量指定
亚声速出口：必须指定压力
超声出口：完全外推

初始条件对收敛速度影响很大。

激波处理

问题：

👉 激波附近出现振荡

解决：

👉 人工黏性

本质：

加入适当耗散项，抑制激波附近的非物理振荡。

核心 CFD 思想总结

定常问题不一定直接求解

可以转化为：非定常问题 + 时间推进

方程形式极其重要

平滑流：非守恒与守恒差别不大
含激波：非守恒与守恒差别巨大

核心 CFD 思想总结

边界条件决定解

尤其是亚声速出口压力。

数值稳定性是关键

需要：

CFL 约束
人工黏性

CFD 本质不是求“精确解”，而是求：

稳定且物理正确的解

CFD 实际流程

一个标准 CFD 流程如下：

建立物理模型：控制方程、状态方程

选择变量形式：守恒形式优先

离散化：空间差分、时间推进

设置边界条件：入口、出口、壁面

CFD 实际流程

给定初始条件

时间推进计算

收敛判断

结果验证：与解析解、实验或高可信计算结果对比

这也是从简单喷管算例走向复杂工程 CFD 的基本路径。

向更高级 CFD 发展的方向

高分辨率格式：

TVD
ENO / WENO

Riemann 求解器：

Roe
HLLC

隐式方法：

LU-SGS
ADI

多维 CFD：

本章总结

喷管 CFD = 守恒方程 + 时间推进 + 正确边界 + 必要耗散

CFD 的本质是：
在数值稳定性与物理真实性之间寻找平衡。

作业：守恒形式计算纯亚声速喷管流

参考课件中的准一维喷管流代码，采用守恒形式计算纯亚声速收缩—扩张喷管流，并与解析解对比。

要求：

使用守恒形式控制方程（无量纲，推导过程参考教材第七章）：

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

数值方法采用 MacCormack 预测—校正格式。
出口为亚声速出口，需要给定出口压力：

\[ p_N=p_e \]

计算结束后要求全场满足：

\[ M<1 \]

作业：守恒形式计算纯亚声速喷管流

喷管几何

计算区域：

\[ 0\le x\le 3 \]

面积分布：

\[ \frac{A}{A_t}=1+2.2(x-1.5)^2,\qquad 0\le x\le 1.5 \]

\[ \frac{A}{A_t}=1+0.2223(x-1.5)^2,\qquad 1.5\le x\le 3.0 \]

作业：守恒形式计算纯亚声速喷管流

基本算例

\[ \gamma=1.4,\qquad CFL=0.5,\qquad N=61 \]

出口压比可先取：

\[ p_e/p_0=0.95 \]

计算得到稳态后，画出：

$M(x)$；
$p(x)$；
$\rho(x)$；
$T(x)$；
残差随迭代步数变化曲线。

作业：守恒形式计算纯亚声速喷管流

解析解对比

利用等熵关系和面积—马赫数关系计算解析解。注意本题为纯亚声速流动，求解马赫数时始终取亚声速分支。

要求将数值解与解析解画在同一张图中，比较：

\[ M(x),\qquad p(x),\qquad \rho(x),\qquad T(x) \]

并计算一个误差指标，例如：

\[ E_M=\max_i|M_i^{num}-M_i^{exact}| \]

作业：守恒形式计算纯亚声速喷管流

参数影响分析

选择下面 3 类参数进行讨论。

1. 出口压比

建议取：

\[ p_e/p_0=0.93,\quad 0.95,\quad 0.97 \]

讨论出口压比变化对喉部马赫数、最大马赫数和流场分布的影响。

作业：守恒形式计算纯亚声速喷管流

参数影响分析

2. 出口边界处理方式

比较两种实现：

方式一：

\[ \rho_N=2\rho_{N-1}-\rho_{N-2},\qquad T_N=\frac{p_N}{\rho_N} \]

方式二：

\[ T_N=2T_{N-1}-T_{N-2},\qquad \rho_N=\frac{p_N}{T_N} \]

讨论哪种方式与解析解更接近。

3. 初始条件

尝试改变初始 $\rho,V,T$ 分布，讨论初始条件对最终稳态和收敛速度的影响。

作业要求

提交一份打印版作业，包含：

无量纲控制方程；
关键代码说明；
数值解与解析解对比图；
参数影响分析。

于 5月27日随堂提交打印版报告。

对比项	守恒形式	非守恒形式
基本形式	\(\dfrac{\partial U}{\partial t}+\dfrac{\partial F}{\partial x}=0\)	\(\dfrac{\partial q}{\partial t}+A(q)\dfrac{\partial q}{\partial x}=0\)
光滑流动	适用	适用
激波间断	适用	不可靠
跳跃关系	自动满足 Rankine–Hugoniot 条件	不一定满足
激波位置	通常更可靠	可能偏移
物理守恒性	控制体意义下守恒	可能破坏守恒

流体力学核心贯通课 II

综合算例 准一维喷管流动

参考书

流动描述：收缩–扩张喷管中的等熵流动

流动描述：收缩–扩张喷管中的等熵流动

流动描述：收缩–扩张喷管中的等熵流动

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❶ 从质量守恒出发

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❷ 引入状态方程

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❸ 整理为压力和温度形式

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❹ 代入等熵关系

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❺ 得到面积关于马赫数的表达式

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❻ 引入临界面积

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❼ 临界面积表达式

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❽ 最终关系式

流动描述：收缩–扩张喷管中的等熵流动

面积–马赫数关系 ❾ 物理含义

流动描述：收缩–扩张喷管中的等熵流动

等熵关系

流动描述：收缩–扩张喷管中的等熵流动

等熵关系的物理含义

CFD求解：MacCormack方法

控制方程推导：从积分形式出发

连续方程 ❶ 积分形式与控制体

控制方程推导：从积分形式出发

连续方程 ❷ 体积分项

控制方程推导：从积分形式出发

连续方程 ❸ 表面积分项

控制方程推导：从积分形式出发

连续方程 ❹ 守恒形式

控制方程推导：从积分形式出发

连续方程 ❺ 非守恒形式

控制方程推导：从积分形式出发

连续方程 ❻ 非守恒形式

控制方程推导：从积分形式出发

动量方程

控制方程推导：从积分形式出发

动量方程

控制方程推导：从积分形式出发

能量方程 ❶ 积分形式

控制方程推导：从积分形式出发

最终控制方程汇总

MacCormack方法

预测步（forward）❶ 密度方程

MacCormack方法

预测步（forward）❷ 速度方程

MacCormack方法

预测步（forward）❸ 温度方程

MacCormack方法

预测步（forward）❹ 更新预测值

MacCormack方法

校正步（backward）❶ 密度方程

MacCormack方法

校正步（backward）❷ 速度方程

MacCormack方法

校正步（backward）❸ 温度方程

MacCormack方法

平均时间导数

MacCormack方法

最终更新

MacCormack方法

时间步长：\(\Delta t\) 应该取多大？❶

MacCormack方法

时间步长：\(\Delta t\) 应该取多大？❷

MacCormack方法

时间步长：\(\Delta t\) 应该取多大？❸

MacCormack方法

时间步长：\(\Delta t\) 应该取多大？❹

MacCormack方法

局部时间步要求

方案 1：局部时间步（Local Time Stepping）

方案 2：全局时间步（推荐）

时间步长

综合算例准一维喷管流动