用Python动态推导正激波方程从理论到可视化的沉浸式学习正激波方程是流体力学和空气动力学中的核心概念之一但对于许多学习者来说纯数学推导过程往往令人望而生畏。想象一下当你面对一页页充满偏微分符号的推导过程时是否曾感到困惑和挫败这正是传统教学方法的一个痛点——它把生动的物理现象抽象成了冰冷的数学符号。本文将带你用Python重新探索这一经典理论通过代码实现从基础方程到完整正激波关系的推导并实时可视化参数变化让抽象概念变得触手可及。1. 准备工作搭建Python科学计算环境在开始推导之前我们需要配置一个适合符号计算和科学可视化的Python环境。推荐使用Anaconda发行版它集成了我们所需的所有关键库。# 安装必要库如果尚未安装 # conda install sympy matplotlib numpy -y核心工具库介绍SymPyPython的符号计算库可以像在黑板上一样进行公式推导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, positiveTrue) # 比热比 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(fa_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(gamma1.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-, linewidth2) plt.title(正激波前后参数变化 (γ%.1f) % gamma) plt.ylabel(密度比 ρ₂/ρ₁) plt.grid(True) plt.subplot(3, 1, 2) plt.plot(M1_values, pressure_ratios, r-, linewidth2) plt.ylabel(压强比 p₂/p₁) plt.grid(True) plt.subplot(3, 1, 3) plt.plot(M1_values, temperature_ratios, g-, linewidth2) 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(gamma1.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-, linewidth2) plt.axvline(x1, colork, linestyle--) plt.axhline(y0, colork, 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, gamma1.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, positiveTrue) 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_max5): 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-, linewidth2) plt.title(f正激波前后参数变化 (γ{self.gamma})) plt.ylabel(密度比 ρ₂/ρ₁) plt.grid(True) plt.subplot(2, 2, 2) plt.plot(M1_values, pressure_ratios, r-, linewidth2) plt.ylabel(压强比 p₂/p₁) plt.grid(True) plt.subplot(2, 2, 3) plt.plot(M1_values, temperature_ratios, g-, linewidth2) plt.ylabel(温度比 T₂/T₁) plt.xlabel(激波前马赫数 M₁) plt.grid(True) plt.subplot(2, 2, 4) plt.plot(M1_values, M2_values, k-, linewidth2) 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(gamma1.4) shock_calc.plot_shock_relations()这个类封装了所有正激波关系的计算和可视化功能读者可以轻松修改参数如比热比γ来探索不同条件下的激波特性。