简介面向流体力学研究者、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)压力可以改写成 (pmRT/u)所以未知量只剩下 (u(x)) 和 (T(x))。从实现的稳健性看我一般不会直接去积分二阶动量方程而是先利用上游条件定出 (C_M) 和 (C_E)再处理一阶系统。这里需要记住一个容易忽略的点跨过激波时总焓 (c_pTu^2/2) 在上游和下游相等但在激波内部并不恒等于常数因为热传导项 (-\kappa T) 在波内不为零。若先把总焓当成常数代入动量方程求出来的激波厚度会偏小强激波时误差更明显。2.2 为什么激波结构能写成“精确解”可压缩 NS 方程一般没有解析解但粘性激波结构是例外。原因在于激波两侧的边界条件不是随便给的而是由 Rankine-Hugoniot 关系强制固定的。动量守恒积分在 (x\to-\infty) 时给出 (C_Mm u_1p_1)在 (x\to\infty) 时又要求 (C_Mm u_2p_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 C0 ]这样动量方程就变成了标准的一维 Burgers 型方程。这个模型虽然牺牲了强激波下的精确性但抓住了粘性激波最本质的对流-扩散平衡。工程上做算法验证时先在这个模型里把代码跑通再上完整 NS 求解器是常见且稳妥的做法。2.3 弱激波极限下的 tanh 解与粘性长度对上面的 Burgers 型方程积分得到闭式解[ u(x)\frac{u_1u_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, gamma1.4, mu1.8e-5, T1300.0, R287.0, p1101325.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.50.53738 nm152 nm2.00.37529 nm116 nm3.00.25919 nm77 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}{(\gamma1)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(fu1{u1:.2f} m/s, u2{u2:.2f} m/s) print(fdelta{delta*1e9:.2f} nm) print(f数值宽度{width_num*1e9:.2f} nm, 理论宽度{width_theory*1e9:.2f} nm) plt.plot(x * 1e6, u_profile, lw2) plt.xlabel(x (μm)) plt.ylabel(u (m/s)) plt.title(fViscous 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_1300\text{K})(p_1101325\text{Pa})分子粘性取 (\mu1.8\times10^{-5}\text{Pa·s})。for M in [1.2, 1.5, 2.0, 3.0]: delta_m shock_delta(M) print(fM{M:4.1f} delta{delta_m*1e9:7.2f} nm f4delta{4*delta_m*1e9:7.2f} nm)结果如下(M_1)(\delta)理论激波宽度 (4\delta)建议网格尺度1.248 nm192 nm≤ 20 nm1.538 nm152 nm≤ 15 nm2.029 nm116 nm≤ 12 nm3.019 nm77 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_12)、(\delta29\text{nm}) 为例若网格取 (290\text{nm})人工粘性带来的宽度会比物理宽度大一个数量级。更隐蔽的是某些高阶格式的人工耗散依赖马赫数导致不同马赫数下同一套网格的有效激波厚度不一致。这时如果拿弱激波解析解去校准格式会得到一个马赫数相关的“伪验证结论”。我的经验是先在同一网格上分别跑 (M_11.5) 和 (M_13.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 解或直接求解耦合的一阶方程组。本文还有配套的精品资源点击获取
粘性激波结构解析解：从NS方程到CFD网格验证的标尺
简介面向流体力学研究者、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)压力可以改写成 (pmRT/u)所以未知量只剩下 (u(x)) 和 (T(x))。从实现的稳健性看我一般不会直接去积分二阶动量方程而是先利用上游条件定出 (C_M) 和 (C_E)再处理一阶系统。这里需要记住一个容易忽略的点跨过激波时总焓 (c_pTu^2/2) 在上游和下游相等但在激波内部并不恒等于常数因为热传导项 (-\kappa T) 在波内不为零。若先把总焓当成常数代入动量方程求出来的激波厚度会偏小强激波时误差更明显。2.2 为什么激波结构能写成“精确解”可压缩 NS 方程一般没有解析解但粘性激波结构是例外。原因在于激波两侧的边界条件不是随便给的而是由 Rankine-Hugoniot 关系强制固定的。动量守恒积分在 (x\to-\infty) 时给出 (C_Mm u_1p_1)在 (x\to\infty) 时又要求 (C_Mm u_2p_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 C0 ]这样动量方程就变成了标准的一维 Burgers 型方程。这个模型虽然牺牲了强激波下的精确性但抓住了粘性激波最本质的对流-扩散平衡。工程上做算法验证时先在这个模型里把代码跑通再上完整 NS 求解器是常见且稳妥的做法。2.3 弱激波极限下的 tanh 解与粘性长度对上面的 Burgers 型方程积分得到闭式解[ u(x)\frac{u_1u_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, gamma1.4, mu1.8e-5, T1300.0, R287.0, p1101325.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.50.53738 nm152 nm2.00.37529 nm116 nm3.00.25919 nm77 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}{(\gamma1)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(fu1{u1:.2f} m/s, u2{u2:.2f} m/s) print(fdelta{delta*1e9:.2f} nm) print(f数值宽度{width_num*1e9:.2f} nm, 理论宽度{width_theory*1e9:.2f} nm) plt.plot(x * 1e6, u_profile, lw2) plt.xlabel(x (μm)) plt.ylabel(u (m/s)) plt.title(fViscous 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_1300\text{K})(p_1101325\text{Pa})分子粘性取 (\mu1.8\times10^{-5}\text{Pa·s})。for M in [1.2, 1.5, 2.0, 3.0]: delta_m shock_delta(M) print(fM{M:4.1f} delta{delta_m*1e9:7.2f} nm f4delta{4*delta_m*1e9:7.2f} nm)结果如下(M_1)(\delta)理论激波宽度 (4\delta)建议网格尺度1.248 nm192 nm≤ 20 nm1.538 nm152 nm≤ 15 nm2.029 nm116 nm≤ 12 nm3.019 nm77 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_12)、(\delta29\text{nm}) 为例若网格取 (290\text{nm})人工粘性带来的宽度会比物理宽度大一个数量级。更隐蔽的是某些高阶格式的人工耗散依赖马赫数导致不同马赫数下同一套网格的有效激波厚度不一致。这时如果拿弱激波解析解去校准格式会得到一个马赫数相关的“伪验证结论”。我的经验是先在同一网格上分别跑 (M_11.5) 和 (M_13.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 解或直接求解耦合的一阶方程组。本文还有配套的精品资源点击获取