简介：面向流体力学研究者、CFD工程师及高年级本科生，提供一维Navier-Stokes方程粘性激波结构的精确解与数值实现。Navier-Stokes方程本身多为非线性偏微分方程组，解析解稀少，而粘性激波恰能体现黏性耗散下的流动突变过渡，这份资源正好围绕该问题给出可运行的求解程序与后处理脚本，能够作为检验有限差分、有限元及谱方法精度的对照基准，也可用于课程设计或论文中的CFD算例验证。压缩包共11个文件、199KB，以f90、m、dat、exe为主：f90为求解主程序，m负责绘制激波密度、速度剖面，dat保存计算数据，exe便于直接运行复现；另有pdb调试符号和dsp、dsw、opt等工程配置文件，方便二次修改与断点调试。已有296人学习下载。借助该包可获得完整的一维粘性激波数值实验流程，理解激波内部从上游到下游状态量的连续过渡，并学会用解析解评估各类数值格式的分辨率与耗散特性，是流体力学与CFD学习者很实用的参考资料。
1. 为什么 NS 方程里会存在一个能写出解析解的“粘性激波结构”
如果把一维定常激波放进 Navier-Stokes(NavierStokes) 方程的坐标系里，速度剖面并不是流体力学教科书里画的那条竖直线，而是一条由粘性主导的光滑过渡曲线。激波不再是一个数学间断，而是一段厚度只有几十纳米到几微米的连续过渡层，这就是 Viscous_shock_structure，也是 NS 方程里极少数能写出闭式解的流动之一。粘性激波能被解析求解，核心原因是动量通量积分和 Rankine-Hugoniot 条件把压力、密度和速度锁在同一个一阶微分关系里；马赫数接近 1 时，速度剖面直接退化成双曲正切函数。本文从方程化简开始，给出可复现代码，并说明参数量级和 CFD 验证里常见的坑。
2. 从一维 Navier-Stokes 方程到粘性激波的结构解
2.1 激波参考系下的守恒积分与粘性应力
一维定常、无外力、无化学反应的压缩性 NS 方程，在激波静止的参考系里可以压缩成三个守恒积分。连续方程给出质量通量常数：
[ \rho u = m ]
动量方程把压力、动量通量和粘性应力合并成一个守恒量：
[ m u + p - \mu_{\mathrm{eff}} \frac{du}{dx} = C_M ]
其中 (\mu_{\mathrm{eff}}=4\mu/3) 来自 Stokes 假设下的纵向粘性应力系数。能量方程同样可以写成通量形式：
[ m\left(c_p T+\frac{u^2}{2}\right) - \kappa \frac{dT}{dx} = C_E ]
这三个积分式的价值在于，它们把“激波内部发生了什么”从偏微分方程问题降成了常微分方程问题。对理想气体 (p=\rho R T)，再利用 (m=\rho u)，压力可以改写成 (p=mRT/u)，所以未知量只剩下 (u(x)) 和 (T(x))。
从实现的稳健性看，我一般不会直接去积分二阶动量方程，而是先利用上游条件定出 (C_M) 和 (C_E)，再处理一阶系统。这里需要记住一个容易忽略的点：跨过激波时总焓 (c_pT+u^2/2) 在上游和下游相等，但在激波内部并不恒等于常数，因为热传导项 (-\kappa T') 在波内不为零。若先把总焓当成常数代入动量方程，求出来的激波厚度会偏小，强激波时误差更明显。
2.2 为什么激波结构能写成“精确解”
可压缩 NS 方程一般没有解析解，但粘性激波结构是例外。原因在于激波两侧的边界条件不是随便给的，而是由 Rankine-Hugoniot 关系强制固定的。动量守恒积分在 (x\to-\infty) 时给出 (C_M=m u_1+p_1)，在 (x\to+\infty) 时又要求 (C_M=m u_2+p_2)。把这个差写成余量函数：
[ S(u)=p(u)+m(u-u_1)-p_1 ]
它同时满足 (S(u_1)=0) 和 (S(u_2)=0)。也就是说，激波结构解本质上是寻找一条从 (u_1) 单调过渡到 (u_2) 的轨道，而这一条轨道被两个守恒常数卡死了。
如果把温度和压力之间的关系再代入，动量方程可以整理成更直观的形式：
[ \mu_{\mathrm{eff}} \frac{du}{dx} = S(u) ]
当速度跳跃比较小，也就是上游马赫数接近 1 时，(S(u)) 可以近似展开成两个根的乘积：
[ S(u)\approx -C(u-u_1)(u-u_2), \quad C>0 ]
这样动量方程就变成了标准的一维 Burgers 型方程。这个模型虽然牺牲了强激波下的精确性，但抓住了粘性激波最本质的对流-扩散平衡。工程上做算法验证时，先在这个模型里把代码跑通，再上完整 NS 求解器，是常见且稳妥的做法。
2.3 弱激波极限下的 tanh 解与粘性长度
对上面的 Burgers 型方程积分，得到闭式解：
[ u(x)=\frac{u_1+u_2}{2}-\frac{u_1-u_2}{2}\tanh\left(\frac{x-x_0}{2\delta}\right) ]
其中 (x_0) 是激波中心位置，由坐标平移自由性决定；(\delta) 是激波结构的特征长度。从量纲上看，(\delta) 与 (\mu_{\mathrm{eff}}/m) 同阶，也就是粘性系数与质量通量的比值。空气中常温常压下这个值大约在几十纳米量级。
下面这串代码可以直接估算不同条件下 (\delta) 的量级：
import numpy as np def shock_delta(M1, gamma=1.4, mu=1.8e-5, T1=300.0, R=287.0, p1=101325.0): """ 估算粘性激波的特征长度 delta。 M1: 上游马赫数 返回 delta，单位米 """ rho1 = p1 / (R * T1) a1 = np.sqrt(gamma * R * T1) m = rho1 * M1 * a1 # 质量通量 mu_eff = 4.0 / 3.0 * mu # Stokes 假设 return mu_eff / m print(shock_delta(2.0))代码里最关键的参数是质量通量 (m)，它由上游马赫数直接决定。(\mu_{\mathrm{eff}}) 用的是 Stokes 假设下的纵向粘性系数，如果你做二维或三维模拟，需要用完整的应力张量重新写等效系数，不能照抄这个一维表达式。
下表列出了常温常压空气在不同马赫数下的特征长度量级：
| (M_1) | (u_2/u_1) | (\delta) 量级 | (4\delta) 量级 |
|---|---|---|---|
| 1.5 | 0.537 | 38 nm | 152 nm |
| 2.0 | 0.375 | 29 nm | 116 nm |
| 3.0 | 0.259 | 19 nm | 77 nm |
注意 (\delta) 是“粘性长度”，不是激波的视觉宽度。对于 tanh 解，最大斜率出现在 (x_0) 处，斜率为 ((u_1-u_2)/(4\delta))，所以如果把激波宽度定义为速度差除以最大斜率，宽度正好是 (4\delta)。这个关系后面做数值验证时会反复用到。
2.4 Prandtl 数和完全可压情形的边界
上面给出的 tanh 解来自弱激波近似，严格说它没有完整包含能量方程。要做完整可压缩 NS 方程的粘性激波精确解，普朗特数会闯进来。最经典的完全解称为 Becker 解，它假设 (\mathrm{Pr}=3/4)，并允许粘性系数随温度变化。在这个条件下，能量方程和动量方程可以再次解耦，最终仍然得到一个与温度梯度耦合的一阶系统，只是形式上更复杂。
很多人在自编 CFD 程序时直接把弱激波 tanh 解当成完整 NS 解来验证二阶格式，这在 (M_1) 小于 1.3 左右时误差不明显，但 (M_1) 超过 2 后，激波内部温度和密度的非对称性会显著偏离 tanh 曲线。我的建议是：先确认模拟对象的马赫数范围，再决定使用弱激波近似解还是数值求解完整的一阶守恒方程系统。不要把 tanh 解这个“标尺”本身校歪了。
3. 用 Python 把粘性激波剖面复现出来：最小可运行代码
3.1 参数初始化：从 M1、T1、p1 到 Rankine-Hugoniot
生成解析剖面的第一步，是用正激波关系求出下游速度。对理想气体，上游马赫数 (M_1) 和速度比满足：
[ \frac{u_2}{u_1}=\frac{2+(\gamma-1)M_1^2}{(\gamma+1)M_1^2} ]
这里不要用等熵关系代替 Rankine-Hugoniot，激波是强耗散过程，等熵假设会在马赫数稍高时引入明显误差。代码里先设置气体常数和上游状态，再计算 (u_2)，最后用 (\mu_{\mathrm{eff}}/m) 得到 (\delta)。
3.2 生成解析剖面并做宽度自检
下面这段代码完整生成粘性激波的速度剖面，并自动验证最大斜率与厚度的关系：
import numpy as np import matplotlib.pyplot as plt # 气体与上游状态 gamma = 1.4 M1 = 2.0 T1 = 300.0 # K p1 = 101325.0 # Pa R = 287.0 # J/(kg K) cp = gamma * R / (gamma - 1) # 上游声速与速度 a1 = np.sqrt(gamma * R * T1) u1 = M1 * a1 rho1 = p1 / (R * T1) m = rho1 * u1 # 质量通量 # 正激波关系 u2 = u1 * (2.0 + (gamma - 1.0) * M1 * M1) / ((gamma + 1.0) * M1 * M1) # 粘性长度 mu_eff = 4.0 / 3.0 * 1.8e-5 delta = mu_eff / m # 生成剖面 x0 = 0.0 x = np.linspace(-15.0 * delta, 15.0 * delta, 3000) u_profile = 0.5 * (u1 + u2) - 0.5 * (u1 - u2) * np.tanh((x - x0) / (2.0 * delta)) # 用最大斜率自检宽度 du = np.gradient(u_profile, x) slope_max = np.max(np.abs(du)) width_num = (u1 - u2) / slope_max width_theory = 4.0 * delta print(f"u1={u1:.2f} m/s, u2={u2:.2f} m/s") print(f"delta={delta*1e9:.2f} nm") print(f"数值宽度={width_num*1e9:.2f} nm, 理论宽度={width_theory*1e9:.2f} nm") plt.plot(x * 1e6, u_profile, lw=2) plt.xlabel("x (μm)") plt.ylabel("u (m/s)") plt.title(f"Viscous shock structure, M1={M1}") plt.grid(True) plt.show()先说参数逻辑：delta = mu_eff / m是粘性尺度的定义，直接把粘性应力和对流质量通量放在一起比较，量纲正好是长度。x = np.linspace(-15*delta, 15*delta, 3000)取 15 倍特征长度，意味着两端已经接近均匀流，边界截断误差可以忽略。
自检部分用np.gradient计算最大斜率，再按 tanh 剖面的理论关系反推宽度。如果输出中两个宽度接近，说明生成的剖面确实是解析解的正确采样；如果偏差超过 1%，多半是网格没有覆盖到足够远处，或者剖面两端过于接近，边界处的导数没有衰减到零。对于强激波，速度剖面会略微不对称，宽度自检误差会变大，这属于预期行为，不是代码 bug。
3.3 剖面参数怎么看
| 参数 | 含义 | 对剖面的影响 |
|---|---|---|
| (x_0) | 激波中心位置 | 平移整条曲线，不影响形状和宽度 |
| (\delta) | 粘性特征长度 | 控制剖面拉伸，也决定最大斜率 |
| (u_1-u_2) | 激波强度 | 决定速度跳变和剖面整体高度差 |
| (M_1) | 上游马赫数 | 通过正激波关系决定 (u_2) 和 (\delta) |
实际使用中，(x_0) 经常被设置为网格中点。如果模拟区间远大于激波厚度，(x_0) 的选取不影响验证结论，因为解析解本身具有平移不变性。但要注意：CFD 计算的激波位置取决于初场和边界条件，比较解析解时，需要先对数值剖面做对齐，再计算误差。最稳妥的方式是让激波中心固定在一个已知网格点上，比如把初场的激波中心放在计算域正中间。
4. 激波厚度、马赫数和网格分辨率的三角关系：坑都在这里
4.1 一张表看不同马赫数下的激波厚度量级
用前面代码里的公式，可以快速生成不同马赫数下的厚度表，直接指导网格设计。典型空气状态，(T_1=300\text{K})，(p_1=101325\text{Pa})，分子粘性取 (\mu=1.8\times10^{-5}\text{Pa·s})。
for M in [1.2, 1.5, 2.0, 3.0]: delta_m = shock_delta(M) print(f"M={M:4.1f} delta={delta_m*1e9:7.2f} nm " f"4delta={4*delta_m*1e9:7.2f} nm")结果如下：
| (M_1) | (\delta) | 理论激波宽度 (4\delta) | 建议网格尺度 |
|---|---|---|---|
| 1.2 | 48 nm | 192 nm | ≤ 20 nm |
| 1.5 | 38 nm | 152 nm | ≤ 15 nm |
| 2.0 | 29 nm | 116 nm | ≤ 12 nm |
| 3.0 | 19 nm | 77 nm | ≤ 8 nm |
这个量级常常让做可压缩流的人意外：激波厚度不是自由参数，而是由分子输运性质决定的。用求解 Euler 方程的激波捕捉格式时，数值激波宽度通常被强制压到 2 到 3 个网格；但如果求解的是真实 NS 方程，激波宽度必须从物理尺度中来，网格太粗会把物理激波直接抹掉，误差会以伪熵增的形式表现出来。
4.2 为什么网格雷诺数是第一个要检查的量
CFD 里判断一个网格能否分辨粘性激波，最直接的指标不是网格数，而是网格雷诺数：
[ Re_{\Delta x}=\frac{u_1\Delta x}{\nu} ]
其中 (\nu=\mu/\rho_1)。对比激波厚度 (\delta)，可以发现：
[ Re_{\Delta x} = \frac{\Delta x}{\delta} ]
这里取 (\mu_{\mathrm{eff}}) 和 (\mu) 的差异忽略不计。所以网格雷诺数本质上就是“一个网格占几个激波特征长度”。经验上，要在激波内部放 5 个以上网格点，需要 (\Delta x < \delta)；如果只要求厚度相对误差小于 10%，一般需要 (\Delta x \approx 0.2\delta\sim0.3\delta)。
很多求解器用了高阶格式，比如五阶 WENO，看起来激波很“锐利”，但锐利来自格式的人工耗散，并不代表真实粘性厚度。用解析解对照时，不要只看激波位置，要用最大斜率或者积分宽度去对厚度，否则很容易把数值效应误判成物理结果。
4.3 人工粘性把“精确解”伪装成误差
一阶迎风格式相当于在 NS 方程里额外加入一个人工扩散项，等价粘性系数大概为：
[ \mu_{\text{num}} \sim \rho u \Delta x/2 ]
以 (M_1=2)、(\delta=29\text{nm}) 为例，若网格取 (290\text{nm})，人工粘性带来的宽度会比物理宽度大一个数量级。更隐蔽的是，某些高阶格式的人工耗散依赖马赫数，导致不同马赫数下同一套网格的有效激波厚度不一致。这时如果拿弱激波解析解去校准格式，会得到一个马赫数相关的“伪验证结论”。
我的经验是：先在同一网格上分别跑 (M_1=1.5) 和 (M_1=3.0) 两组算例，看厚度误差是否随马赫数漂移。真正的粘性激波解析解对马赫数的依赖是确定的，如果两组算例的厚度误差方向相反，几乎都可以归因于格式人工粘性或边界条件处理不当。
5. 用解析激波厚度给 CFD 结果做“内嵌标尺”，反推数值扩散
5.1 从最大斜率提取数值激波宽度
做完 CFD 模拟后，从输出数据里提取激波剖面并计算厚度，是最快的验证动作。已知 tanh 解中，速度差除以最大斜率等于 (4\delta)，于是可以从数值剖面反推等效粘性长度：
from numpy import gradient, abs, argmax, diff def numeric_width(x_cfd, u_cfd): """ 从CFD剖面计算激波宽度。 返回宽度、最大斜率位置、等效粘性长度。 """ du = gradient(u_cfd, x_cfd) idx = argmax(abs(du)) width = (u_cfd.max() - u_cfd.min()) / abs(du[idx]) delta_num = width / 4.0 return width, idx, delta_num这个函数为什么可靠：tanh 剖面的最大斜率固定出现在中心点，而且中心点的速度正好是上下游的平均值。如果数值剖面中心速度偏离平均值，说明激波结构不对称，此时再用这个公式算等效宽度，只能当作粗估；更严格的思路是对 (\partial u/\partial x) 做积分包络，但工程验证中最大斜率法足够直观。
5.2 反推数值粘性并换算成数值雷诺数
得到 (\delta_{\text{num}}) 后，可以反推求解器实际“感知”的粘性系数：
[ \mu_{\text{num}} = m \delta_{\text{num}} ]
再用它和物理 (\mu_{\mathrm{eff}}) 的比值判断格式的人工扩散水平：
[ R_{\mu} = \frac{\mu_{\text{num}}}{\mu_{\mathrm{eff}}} ]
当 (R_\mu) 接近 1 时，说明网格已经足以分辨物理粘性激波；当 (R_\mu) 大于 5 时，激波宽度主要由格式人工扩散控制。我常用的一个验收线是：在高阶格式中，(R_\mu) 应小于 2；在稳定的二阶格式中，小于 3 即可认为激波处的 NS 方程已基本被分辨。若超出这个范围，优先加密网格，而不是调大格式的耗散参数。
这个方法比直接看激波位置误差更有效，因为激波位置可以通过相同的通量函数轻易保持一致，而厚度则暴露了求解器在激波内部的真实耗散水平。把 (R_\mu) 随网格加密的收敛曲线打出来，如果斜率接近格式的理论阶数，说明粘性激波结构已经进入网格收敛区；斜率掉下来时，通常就是解析解本身的适用边界——激波变强、剖面不再对称，需要换用完整 Becker 解或直接求解耦合的一阶方程组。
本文还有配套的精品资源，点击获取
