简介面向空气动力学与计算流体力学学习者SSE.zip聚焦激波间断问题用于理解气动力计算中强间断的数值处理方法。压缩包共22个文件大小仅770KB包含Visual Studio工程文件.sln/.vcxproj、C源码.cpp、可执行程序.exe、数据文件.dat、输入说明.txt等既有工程配置也有运行模块便于直接查看整体结构、编译调试或运行验证。已有110人学习下载适合具备一定流体力学基础、希望学习激波捕获或间断求解算法的学生与工程师。内容以数值模拟源码为主附带必要的输入输出文件能帮助读者从代码层面掌握SSE格式在建模仿真中的应用并依据示例数据对照分析激波位置与气动力参数的变化。由于压缩包体积小也适合快速上手试验自行修改算例参数观察激波间断的捕捉效果是一份轻量但完整的入门实践资料。1. 气动力计算遇到激波间断时SSE 是最后那层兜底标题里 SSE.zip 这份代码核心不是某个新格式而是可压缩气动力求解器里最常被低估的兜底组合Shock sensor激波传感器、Slope limiter限制器和 Entropy fix熵修正三个模块合起来才叫 SSE。做气动力计算的人大概都有过这种经历网格够密、格式够高阶翼型绕流在跨声速段却出现非物理膨胀激波或者高超声速钝体头部弓形激波前面长出一截红玉状凸起。这两种现象和网格关系不大瓶颈都在对激波间断的数值处理。这篇文章按我做 CFD 求解器时的习惯把 SSE 拆开讲哪里失效、怎么补、参数怎么设、用什么验证给一套能直接抄到自研代码或 OpenFOAM 自定义通量里的做法。适合正在给层流/湍流气动力程序加激波捕捉的工程师或想搞懂 Roe 通量为什么总在声速点闯祸的研究生。2. SSE 要处理的三个激波间断问题膨胀激波、红玉失稳与过冲2.1 激波在稀疏网格里是“被允许的间断”不是被算出来的波用有限体积法解可压缩欧拉方程时激波从来不是被“画”出来的而是由守恒律的弱解自然允许的速度、密度、压力跳跃。只要通量公式是守恒的激波就会以几个网格宽度自动形成。这个薄层本质上是一个数值间断左右两边状态满足 Rankine-Hugoniot 跳跃条件。问题也出在这里因为激波是间断任何基于 Taylor 展开的离散都不天然适应它。高阶重构在激波两旁会产生 Gibbs 震荡低阶格式又会把激波抹得太宽直接影响气动力。对激波间断的处理几十年来的共识是通量格式提供基础耗散限制器负责在激波附近降回一阶激波传感器负责告诉限制器“哪里降阶”熵修正负责保证这个间断是物理上允许的激波而不是非物理膨胀。这套组合正好对应 SSE 的三个字母。下面把三个失效场景分别说清楚因为它们各自触发条件不同混在一起调参会很痛苦。2.2 声速点膨胀激波Roe 通量的熵缺陷Roe 格式是对两个相邻状态做线性化用 Roe 平均矩阵近似通量 Jacobian。它最大的问题在于当 Roe 平均状态的特征值中有某一个接近零时数值耗散会退化。最典型的是跨声速流动的声速点附近通量计算允许熵减弱的间断通过于是在本该光滑的等熵加速区里出现一串阶梯状的非物理“膨胀激波”。翼型跨声速吸力峰、喷管喉道这两类气动力场景最容易中招。熵修正不是加物理模型而是保证数值耗散在特征值过零时不下掠某个下限。最常见的是 Harten 开关函数[ h(\lambda)\begin{cases} |\lambda|, |\lambda|\ge\delta\ \frac{\lambda^2\delta^2}{2\delta}, |\lambda|\delta \end{cases} ]其中 (\delta) 通常取 (c(|u|a))比例系数 (c) 在 0.02 到 0.2 之间。这个修正把 Roe 格式在声速点的耗散拉回有限量代价是稍微多抹一点激波。3.3 节的代码就是按这个公式实现的跑一遍 Sod 激波管就能直观看到区别。2.3 强激波前的红玉失稳和激波后的过冲分工不同高超声速钝体头部弓形激波在对称轴附近会出一个著名的二维失稳计算得到的激波在轴线处凸起像一颗红宝石所以叫红玉carbuncle现象。这是 Riemann 解算器在强激波背后驻点区域的奇异性问题和边界条件无关用 Roe 这类低耗散通量就会出现。温和版本是激波后方压力过冲沿壁面产生假的高热流气动热计算对压力/温度梯度的误差是按平方放大的。我的处理习惯不对通量格式做外科手术而是让激波传感器在强压缩区把线性特征值的耗散阈值提高到平时的 2 到 5 倍相当于给驻点区域的低速特征加一道阻尼。激波传感器可以用 Ducros 涡量/散度传感器也可以用 Jameson 压力第二差分传感器输出都是 0 到 1 之间的连续指数标明一个网格是否处于间断区域。限制器是另一层防线MUSCL 重构里的 van Leer 或 minmod 限制器控制激波附近斜率不产生新的极值。注意限制器只在局部斜率上工作它区分不了膨胀激波和物理激波所以替代不了熵修正。三个模块的分工如下表。失效现象常见位置数值机理SSE 应对模块参考阈值膨胀激波声速点附近Roe 特征值过零耗散消失EEntropy fix(\delta0.02\sim0.2\times(红玉失稳carbuncle钝体对称轴弓形激波驻点横向扩散不足驻点数值奇性SShock sensor 增强耗散传感器输出大于 0.05 时激活激波后过冲/振铃强激波下游壁面高阶重构产生新极值SSlope limitervan Leer 或 minmod表格里的阈值只是起点。传感器阈值要配合具体算例的来流马赫数调参考文献里常见做法是先跑一阶格式拿到无振荡解作为参考再逐步放大二阶重构并同时提高传感器放大系数。3. 在 1D 欧拉方程上复现最小 SSESod 激波管全程3.1 守恒方程与半离散有限体积一维欧拉方程是验证 SSE 的最小场景一个运动的激波、一个接触间断、一个稀疏波三种间断特征一次集齐。守恒向量和通量向量写在一起就是[ \frac{\partial}{\partial t}\begin{bmatrix}\rho\\rho u\E\end{bmatrix} \frac{\partial}{\partial x}\begin{bmatrix}\rho u\\rho u^2p\u(Ep)\end{bmatrix}0 ]有限体积离散后每个网格单元的守恒量更新为[ U_i^{n1}U_i^n-\frac{\Delta t}{\Delta x}\left(F_{i1/2}-F_{i-1/2}\right) ](F_{i1/2}) 是界面通量由 Roe 格式给出。时间步长不写死每一步用当前最大波速算(dt \mathrm{CFL}\cdot dx / \max(|u|a))其中 (a\sqrt{\gamma p/\rho})。CFL 取 0.4 对 Roe 加熵修正足够稳。3.2 MUSCL 重构、minmod 限制器与激波传感器接在一起二阶实现里每个界面 (i1/2) 需要左右两个状态。MUSCL 做法是把单元中心值向界面外推再套限制器防止出现新极值。对均匀网格中心差值项为[ U_{i1/2}^L U_i \frac{1}{2}\mathrm{minmod}(U_i-U_{i-1},\ U_{i1}-U_i) ]minmod 的规则是左右差同号时取绝对值小者异号时取零。van Leer 限制器更平滑实际气动力程序里我更常用 van Leer实验代码里 minmod 更简单且对 Sod 激波管已经足够。激波传感器用压力二阶差分做这是 Jameson 型传感器最常见形式。对单元 (i) 计算[ s_i \frac{|p_{i1}-2p_ip_{i-1}|}{p_{i1}2p_ip_{i-1}\varepsilon} ]光滑区这个量是 (O(\Delta x^2))激波处直接到 (O(1))。界面 (i1/2) 上取 (s\max(s_i,s_{i1}))这样它在激波两侧都能及时激活又不会在光滑区误报。3.3 带熵修正与传感器增强耗散的 Roe 通量完整代码下面这段代码可以在普通笔记本上 10 秒内跑完。它把熵修正放在 Roe 特征值上把激波传感器接在线性场的耗散上三者都齐了属于 SSE 的最小闭环。import numpy as np GAMMA 1.4 CFL 0.4 N 200 # 内部网格数 DX 1.0 / (N-1) T_END 0.14 # Sod 激波管经典终止时间 # 左右各 2 个 ghost cell总数组长度 M NG 2 M N 4 def prim_to_cons(rho, u, p): E p/(GAMMA-1.0) 0.5*rho*u*u return np.array([rho, rho*u, E]) def cons_to_prim(U): rho U[0] u U[1]/rho E U[2] p (GAMMA-1.0)*(E - 0.5*rho*u*u) return rho, u, p def physical_flux(U): rho, u, p cons_to_prim(U) E U[2] return np.array([rho*u, rho*u*up, u*(Ep)]) def harten_abs(lam, delta): # Harten 熵修正特征值过零时耗散不落到 0 if abs(lam) delta: return abs(lam) return (lam*lam delta*delta) / (2.0*delta) def roe_flux(UL, UR, sensor, delta_base0.05, sensor_amp0.3): FL physical_flux(UL) FR physical_flux(UR) dU UR - UL # 守恒量跳跃 rhoL, uL, pL cons_to_prim(UL) # 原始变量左状态 rhoR, uR, pR cons_to_prim(UR) # 原始变量右状态 # Roe 平均密度、速度、总焓、声速 sL, sR np.sqrt(rhoL), np.sqrt(rhoR) u (sL*uL sR*uR) / (sL sR) HL 0.5*uL*uL GAMMA/(GAMMA-1.0)*pL/rhoL HR 0.5*uR*uR GAMMA/(GAMMA-1.0)*pR/rhoR H (sL*HL sR*HR) / (sL sR) a np.sqrt(max((GAMMA-1.0)*(H - 0.5*u*u), 1e-12)) # 转向特征变量跳跃 rho_avg sL * sR # sqrt(rhoL*rhoR) drho dU[0] dM dU[1] dE dU[2] du (dM - u*drho) / rho_avg dp (GAMMA-1.0)*(dE - u*dM 0.5*u*u*drho) a2 a*a dW1 (dp - rho_avg*a*du) / (2.0*a2) dW2 drho - dp / a2 dW3 (dp rho_avg*a*du) / (2.0*a2) # 三个特征值u-a, u, ua lam1, lam2, lam3 u-a, u, ua # 熵修正下限 delta 随激波传感器放大 delta (delta_base sensor_amp*sensor) * (abs(u) a) abs1 harten_abs(lam1, delta) abs2 harten_abs(lam2, delta) abs3 harten_abs(lam3, delta) # 右特征向量r1, r2, r3 对应三个特征场 r1 np.array([1.0, u-a, H-u*a]) r2 np.array([1.0, u, 0.5*u*u]) r3 np.array([1.0, ua, Hu*a]) # 耗散项 sum |lambda_k| * dW_k * r_k dissip abs1*dW1*r1 abs2*dW2*r2 abs3*dW3*r3 return 0.5*(FL FR) - 0.5*dissip def apply_bc(U): # 零梯度外推内部单元 2..N1 是有效计算域 U[:, 0] U[:, 2] U[:, 1] U[:, 3] U[:, N2] U[:, N1] U[:, N3] U[:, N1] return U def build_smoothness(p): # Jameson 型激波传感器压力二阶差分归一化 # 返回 smooth[k]表示数组第 k1 个单元的压力光滑度 num np.abs(p[:-2] p[2:] - 2.0*p[1:-1]) den p[:-2] 2.0*p[1:-1] p[2:] 1e-12 return 2.0 * num / den # 初始条件Sod 激波管间断在 x0.5 x np.linspace(0.0, 1.0, N) rho0 np.where(x 0.5, 1.0, 0.125) u0 np.zeros(N) p0 np.where(x 0.5, 1.0, 0.1) U np.zeros((3, M)) U[:, NG:NGN] prim_to_cons(rho0, u0, p0) t 0.0 while t T_END: U apply_bc(U) rho, u, p cons_to_prim(U) dt CFL * DX / np.max(np.abs(u) np.sqrt(GAMMA*p/rho)) if t dt T_END: dt T_END - t smooth build_smoothness(p) sensor_face np.zeros(M-1) # 界面 j 位于单元 j 和 j1 之间 for j in range(1, M-2): sensor_face[j] max(smooth[j-1], smooth[j]) Fface np.zeros((3, M-1)) for j in range(M-1): Fface[:, j] roe_flux(U[:, j], U[:, j1], sensor_face[j]) for i in range(NG, NGN): # 内部索引 2..N1 U[:, i] - dt / DX * (Fface[:, i] - Fface[:, i-1]) t dt # 画密度分布验证 import matplotlib.pyplot as plt rho_out, u_out, p_out cons_to_prim(U[:, NG:NGN]) plt.plot(x, rho_out, .-) plt.xlabel(x); plt.ylabel(density) plt.title(Sod shock tube with Roe SSE) plt.savefig(sod_roe_sse.png, dpi150)几点必须说明。代码中smooth[k]表示以数组第 k1 个单元为中心的压力二阶差界面传感器取左右两个中心传感器的最大值这样激波界面两侧都会增强耗散。sensor_amp只作用于线性场特征值 (u) 的熵修正下限不会改变接触间断和稀疏波的耗散主体。ghost cell 用零梯度外推对 Sod 这种纯波传播问题足够干净。运行后密度曲线应当是左边稀疏波光滑过渡中间接触间断约 2 到 3 个网格完成跳跃右边激波约 2 个网格完成跳跃。如果把delta_base调成 0.0稀疏波区尾部会出现可见的小台阶那就是膨胀激波熵修正的作用一眼就能看出来。3.4 三个必调参数delta、sensor_amp 与 CFL参数作用推荐范围调参信号delta_base熵修正下限系数0.02 ~ 0.15太小出现膨胀激波太大抹平接触间断sensor_amp激波处线性场耗散放大0.1 ~ 0.5太大降低激波分辨率太小红玉抑制不住CFL显式时间步稳定因子0.3 ~ 0.6超过格式稳定上限直接震荡发散delta_base的判断标准是接触间断过渡点数。Sod 激波管 200 个网格下过渡点超过 4 个就说明熵修正系数偏大小于 2 个且稀疏波里有锯齿说明系数偏小。sensor_amp在 1D 算例中影响不大主要留给 2D 弓形激波时用。CFL 是最后调的我不建议用大于 0.6 的显式步长去“补”格式的稳定性问题。4. 把 SSE 落到二维气动力网格弓形激波与气动力系数4.1 结构网格上的方向分裂与激波传感器二维化一维公式在二维里按方向分裂扩展对每个坐标方向分别计算 Roe 通量特征值里的速度换成对应方向的法向速度分量声速不变。关键是激波传感器不能按方向独立读取压力否则斜激波在某一方向上的压力变化会被漏掉导致传感器在该方向上不激活。我一般用二维压力拉普拉斯形式[ s_{i,j} \frac{|p_{i-1,j}p_{i1,j}p_{i,j-1}p_{i,j1}-4p_{i,j}|} {p_{i-1,j}p_{i1,j}p_{i,j-1}p_{i,j1}4p_{i,j}\varepsilon} ]这个量和一维版本行为一致光滑区 (O(\Delta x^2))跨过弓形激波时 (O(1))。计算界面通量时取左右网格传感器较大值方向和维度无关直接复用 3.3 节roe_flux里的sensor参数。注意斜激波在结构网格上可能穿过多个界面二维传感器激活区的宽度会比一维宽 1 到 2 个网格这是正常现象。4.2 高超声速钝体弓形激波的参数调节顺序二维算例最容易踩的是同时调多个参数出问题不知道是谁引起的。我的顺序固定三步。第一步先只开熵修正delta_base0.1sensor_amp0跑 200 步看对称轴压力分布。如果轴线驻点压力出现隆起或凹陷说明红玉已经在酝酿。第二步保持delta_base不变把sensor_amp从 0.1 提到 0.3看驻点压力剖面的凸起是否被压平。第三步如果来流马赫数大于 8sensor_amp0.5仍然压不住我会在驻点附近把 Roe 通量切换为 HLLE 格式而不是继续加大耗散。HLLE 的耗散结构对驻点更友好代价是边界层内分辨率变差所以只限驻点附近几层网格使用。这套顺序对应的直观逻辑是膨胀激波在收敛段声速点附近红玉在驻点位置不同、触发条件不同。先修熵条件再修驻点奇性可以避免把两个机制的参数混在一起。飞片、升力体、进气道唇口这类高超声速气动力算例我都按这个顺序处理。4.3 气动力系数对间断处理敏感时的监控指标激波位置和驻点压力直接决定压力系数分布进而影响升力、阻力和压心位置。红玉不消除时驻点 Cp 会偏低整体阻力系数在迭代过程中出现周期波动激波后过冲不消除时壁面热流会偏高而壁面热流误差在校验算例里通常比压力误差放大一个量级。我一般会在每 100 步输出一个激波位置监测值沿壁面找压力梯度最大点import numpy as np def shock_index(p_wall): # p_wall: 沿壁面的一维压力序列 grad np.abs(np.diff(p_wall)) idx int(np.argmax(grad)) return idx, grad[idx]这个量有两个用途判断激波是否在迭代中来回走动以及做网格收敛性检查。粗网格上激波位置偏差 1 到 2 个网格是正常的如果细化网格后激波位置单调漂移超过 5 个网格问题多半不在 SSE 而在远场边界或网格分布此时继续调sensor_amp没有意义。气动力系数对间断处理敏感的另一个典型表现是升力系数残差在某一马赫数区间突然出现低频振荡这时候先看激波位置监测值是否同频摆动再决定是否调整传感器放大系数。5. 验证 SSE 是否生效熵残差、激波位置与三个常踩的坑5.1 用总熵单调性验证熵修正生效不依赖解析解最通用的验证是跟踪总熵。对 Sod 激波管这类含激波的问题总熵在数值解中可能因为人工耗散而缓慢增加但绝不允许在声速点出现明显下降。总熵按如下方式计算def total_entropy(rho, p, dx): s rho * np.log(p / np.power(rho, GAMMA)) return np.sum(s) * dx每 100 步打一次total_entropy。如果曲线在某个区间出现台阶状下降优先怀疑熵修正没生效检查harten_abs里的delta是否在特征值过零时真的被替换而不是只算没被使用。另一个常见错误是把熵修正用到了重构后的界面状态上而没有在 Roe 平均特征值上生效那样修正等于没做。5.2 三个常踩的坑第一个坑是delta_base设到 0.3 以上。接触间断会被抹平成 8 个网格以上激波位置仍然正确但壁面热流系统性偏高。判断标准Sod 激波管密度曲线上接触间断过渡点超过 4 个就回落系数。第二个坑是sensor_amp没有乘上传感器信号而是全局加在abs(u)上。这会污染边界层和尾迹区摩阻偏低。我一般会在纯层流平板算例上把sensor_amp置 0 和置 0.5 各跑一次摩阻差异超过 5% 就说明传感器的选择性有问题优先检查二维传感器分母里是否漏掉了压力偏移项。第三个坑是只在声速点开熵修正不在驻点增强传感器耗散。高超声速钝体的红玉不是声速点熵条件问题而是驻点驻止流的奇异性所以驻点附近必须有传感器增强耗散兜底。判断方法很简单把sensor_amp从 0 改成 0.3若驻点压力剖面无变化说明传感器信号在驻点处为零需要检查压力二阶差分在驻点网格上的激活条件。这三个坑对应三种表现接触间断过抹、边界层摩阻偏低、驻点红玉不消。跑任何气动力算例前先用这一节的总熵监测和激波位置监测过一遍确认 SSE 三个模块各自都在工作再去相信后面的气动力系数结果。本文还有配套的精品资源点击获取