别再死记硬背公式了！用Python手把手带你推导正激波方程（附完整代码）
用Python动态推导正激波方程：从理论到可视化的沉浸式学习
正激波方程是流体力学和空气动力学中的核心概念之一，但对于许多学习者来说，纯数学推导过程往往令人望而生畏。想象一下，当你面对一页页充满偏微分符号的推导过程时，是否曾感到困惑和挫败？这正是传统教学方法的一个痛点——它把生动的物理现象抽象成了冰冷的数学符号。本文将带你用Python重新探索这一经典理论，通过代码实现从基础方程到完整正激波关系的推导，并实时可视化参数变化，让抽象概念变得触手可及。
1. 准备工作：搭建Python科学计算环境
在开始推导之前，我们需要配置一个适合符号计算和科学可视化的Python环境。推荐使用Anaconda发行版，它集成了我们所需的所有关键库。
# 安装必要库（如果尚未安装） # conda install sympy matplotlib numpy -y核心工具库介绍：
- SymPy：Python的符号计算库，可以像在黑板上一样进行公式推导
- Matplotlib：经典的可视化工具，用于绘制激波前后的参数变化
- NumPy：提供数值计算支持，处理推导后的数值计算
import sympy as sp import matplotlib.pyplot as plt import numpy as np from IPython.display import display, Math # 用于美观显示公式提示：本文所有代码均在Jupyter Notebook环境下测试通过，这种交互式环境特别适合教学演示，可以分段执行并立即看到结果。
2. 从三大基本方程出发：建立符号体系
正激波推导的基础是流体力学三大守恒方程：连续性方程、动量方程和能量方程。我们首先用SymPy声明所有需要的符号变量：
# 定义符号变量 rho1, rho2 = sp.symbols('rho_1 rho_2') # 激波前后的密度 u1, u2 = sp.symbols('u_1 u_2') # 速度 p1, p2 = sp.symbols('p_1 p_2') # 压强 T1, T2 = sp.symbols('T_1 T_2') # 温度 gamma = sp.symbols('gamma', positive=True) # 比热比 a1, a2 = sp.symbols('a_1 a_2') # 声速 M1, M2 = sp.symbols('M_1 M_2') # 马赫数 R = sp.symbols('R') # 气体常数2.1 连续性方程的实现
连续性方程表达了质量守恒原理，对于正激波可以简化为：
continuity_eq = sp.Eq(rho1*u1, rho2*u2) display(continuity_eq)输出将显示美观的数学公式：ρ₁u₁ = ρ₂u₂
2.2 动量方程的代码表达
动量守恒方程描述了流体动量的变化与作用力之间的关系：
momentum_eq = sp.Eq(p1 + rho1*u1**2, p2 + rho2*u2**2) display(momentum_eq)2.3 能量方程的符号表示
能量方程考虑热力学第一定律，对于绝热流动可以表示为：
energy_eq = sp.Eq(h1 + u1**2/2, h2 + u2**2/2) # 其中h为比焓，对于理想气体h = cp*T = a^2/(gamma-1) energy_eq = energy_eq.subs({ h1: a1**2/(gamma-1), h2: a2**2/(gamma-1) }) display(energy_eq)3. 关键推导：声速公式与马赫数关系
声速在可压缩流动中扮演着关键角色。让我们从基础热力学关系出发，推导声速公式：
# 声速定义 a = sqrt(gamma*R*T) a = sp.sqrt(gamma*R*T1) display(Math(f"a_1 = {sp.latex(a)}")) # 马赫数定义 M = u/a M1_def = sp.Eq(M1, u1/a1) M2_def = sp.Eq(M2, u2/a2)通过这三个基础方程和定义，我们已经建立了完整的符号体系。接下来就可以进行正激波关系的推导了。
4. 正激波关系的完整推导
4.1 密度比关系推导
从连续性方程出发，结合声速和马赫数定义，我们可以推导出激波前后的密度比：
# 从连续性方程开始 density_ratio = sp.solve(continuity_eq, rho2)[0]/rho1 # 用马赫数表示速度 density_ratio = density_ratio.subs({ u1: M1*a1, u2: M2*a2 }) # 假设激波前后比热比相同，温度比可以表示为声速比的平方 density_ratio = density_ratio.subs(a2/a1, sp.sqrt(T2/T1)) # 最终密度比表达式 display(Math(r"\frac{\rho_2}{\rho_1} = " + sp.latex(density_ratio)))4.2 压强比关系推导
从动量方程出发，类似地可以推导压强比：
pressure_ratio = sp.solve(momentum_eq, p2)[0]/p1 pressure_ratio = pressure_ratio.subs({ u1: M1*a1, u2: M2*a2, rho2: rho1*density_ratio }).simplify() display(Math(r"\frac{p_2}{p_1} = " + sp.latex(pressure_ratio)))4.3 温度比关系推导
结合理想气体状态方程和前面推导的结果，可以得到温度比：
temperature_ratio = (p2/p1)/(rho2/rho1) temperature_ratio = temperature_ratio.subs({ p2/p1: pressure_ratio, rho2/rho1: density_ratio }).simplify() display(Math(r"\frac{T_2}{T_1} = " + sp.latex(temperature_ratio)))5. 可视化分析：激波参数变化规律
理论推导完成后，让我们用Matplotlib创建交互式可视化，直观展示激波前后参数如何随来流马赫数变化。
def plot_shock_relations(gamma=1.4): M1_values = np.linspace(1, 5, 100) # 马赫数范围1到5 # 计算各参数比 density_ratios = (gamma + 1)*M1_values**2 / ((gamma - 1)*M1_values**2 + 2) pressure_ratios = 1 + 2*gamma/(gamma + 1)*(M1_values**2 - 1) temperature_ratios = pressure_ratios / density_ratios plt.figure(figsize=(12, 8)) plt.subplot(3, 1, 1) plt.plot(M1_values, density_ratios, 'b-', linewidth=2) plt.title('正激波前后参数变化 (γ=%.1f)' % gamma) plt.ylabel('密度比 ρ₂/ρ₁') plt.grid(True) plt.subplot(3, 1, 2) plt.plot(M1_values, pressure_ratios, 'r-', linewidth=2) plt.ylabel('压强比 p₂/p₁') plt.grid(True) plt.subplot(3, 1, 3) plt.plot(M1_values, temperature_ratios, 'g-', linewidth=2) plt.ylabel('温度比 T₂/T₁') plt.xlabel('激波前马赫数 M₁') plt.grid(True) plt.tight_layout() plt.show() plot_shock_relations()这段代码会生成三个子图，分别展示密度比、压强比和温度比随来流马赫数的变化曲线。通过调整γ值（比热比），可以观察不同气体性质的差异。
6. 深入理解：为什么亚音速流动不形成激波？
从热力学第二定律的角度，激波过程必须是熵增的。我们可以通过代码计算熵变来验证这一点：
# 熵变计算 s2_s1 = sp.symbols('Delta_s') s2_s1_expr = sp.log((p2/p1)**(1/(gamma-1)) * (rho1/rho2)**(gamma/(gamma-1))) s2_s1_expr = s2_s1_expr.subs({ p2/p1: pressure_ratio, rho2/rho1: density_ratio }).simplify() display(Math(r"\Delta s = " + sp.latex(s2_s1_expr)))通过绘制熵变随马赫数的变化曲线，可以清楚地看到当M₁ < 1时，熵变为负值，这违反了热力学第二定律：
def plot_entropy_change(gamma=1.4): M1_values = np.linspace(0.5, 5, 500) entropy_changes = np.log( (1 + 2*gamma/(gamma + 1)*(M1_values**2 - 1))**(1/(gamma-1)) * ((gamma - 1)*M1_values**2 + 2)/((gamma + 1)*M1_values**2))**(gamma/(gamma-1)) ) plt.figure(figsize=(10, 6)) plt.plot(M1_values, entropy_changes, 'm-', linewidth=2) plt.axvline(x=1, color='k', linestyle='--') plt.axhline(y=0, color='k', linestyle='--') plt.title('正激波熵变随马赫数变化') plt.xlabel('激波前马赫数 M₁') plt.ylabel('熵变 Δs/R') plt.grid(True) plt.show() plot_entropy_change()这张图清晰地展示了为什么马赫数小于1时不满足激波形成的物理条件——因为会导致熵减，违反热力学第二定律。
7. 完整代码实现与交互式探索
为了便于读者实践，以下是完整的Python代码，包含了所有推导步骤和可视化功能：
import sympy as sp import matplotlib.pyplot as plt import numpy as np from IPython.display import display, Math class NormalShockCalculator: def __init__(self, gamma=1.4): self.gamma = gamma self._init_symbols() self._derive_relations() def _init_symbols(self): self.rho1, self.rho2 = sp.symbols('rho_1 rho_2') self.u1, self.u2 = sp.symbols('u_1 u_2') self.p1, self.p2 = sp.symbols('p_1 p_2') self.T1, self.T2 = sp.symbols('T_1 T_2') self.gamma_sym = sp.symbols('gamma', positive=True) self.M1, self.M2 = sp.symbols('M_1 M_2') self.R = sp.symbols('R') def _derive_relations(self): # 连续性方程 self.continuity_eq = sp.Eq(self.rho1*self.u1, self.rho2*self.u2) # 动量方程 self.momentum_eq = sp.Eq(self.p1 + self.rho1*self.u1**2, self.p2 + self.rho2*self.u2**2) # 能量方程 h1 = sp.symbols('h_1') h2 = sp.symbols('h_2') energy_eq = sp.Eq(h1 + self.u1**2/2, h2 + self.u2**2/2) # 用声速表示比焓 a1, a2 = sp.symbols('a_1 a_2') energy_eq = energy_eq.subs({ h1: a1**2/(self.gamma_sym-1), h2: a2**2/(self.gamma_sym-1) }) self.energy_eq = energy_eq # 推导密度比 self.density_ratio = (self.gamma_sym + 1)*self.M1**2 / ( (self.gamma_sym - 1)*self.M1**2 + 2) # 推导压强比 self.pressure_ratio = 1 + 2*self.gamma_sym/(self.gamma_sym + 1)*( self.M1**2 - 1) # 推导温度比 self.temperature_ratio = self.pressure_ratio / self.density_ratio # 推导激波后马赫数 self.M2_expr = sp.sqrt( (1 + (self.gamma_sym-1)/2*self.M1**2) / ( self.gamma_sym*self.M1**2 - (self.gamma_sym-1)/2)) def plot_shock_relations(self, M1_max=5): M1_values = np.linspace(1, M1_max, 100) # 计算各参数比 density_ratios = self._eval_expr(self.density_ratio, M1_values) pressure_ratios = self._eval_expr(self.pressure_ratio, M1_values) temperature_ratios = pressure_ratios / density_ratios M2_values = self._eval_expr(self.M2_expr, M1_values) plt.figure(figsize=(12, 10)) plt.subplot(2, 2, 1) plt.plot(M1_values, density_ratios, 'b-', linewidth=2) plt.title(f'正激波前后参数变化 (γ={self.gamma})') plt.ylabel('密度比 ρ₂/ρ₁') plt.grid(True) plt.subplot(2, 2, 2) plt.plot(M1_values, pressure_ratios, 'r-', linewidth=2) plt.ylabel('压强比 p₂/p₁') plt.grid(True) plt.subplot(2, 2, 3) plt.plot(M1_values, temperature_ratios, 'g-', linewidth=2) plt.ylabel('温度比 T₂/T₁') plt.xlabel('激波前马赫数 M₁') plt.grid(True) plt.subplot(2, 2, 4) plt.plot(M1_values, M2_values, 'k-', linewidth=2) plt.ylabel('激波后马赫数 M₂') plt.xlabel('激波前马赫数 M₁') plt.grid(True) plt.tight_layout() plt.show() def _eval_expr(self, expr, M1_values): # 将符号表达式转换为数值计算函数 f = sp.lambdify((self.M1, self.gamma_sym), expr, 'numpy') return f(M1_values, self.gamma) # 使用示例 shock_calc = NormalShockCalculator(gamma=1.4) shock_calc.plot_shock_relations()这个类封装了所有正激波关系的计算和可视化功能，读者可以轻松修改参数（如比热比γ）来探索不同条件下的激波特性。