简介面向计算流体力学课程学习与数值方法入门这套MATLAB代码围绕一维激波管问题给出可运行的求解示例覆盖Lax-Friedrichs、Lax-Wendroff、Roe与ENO等常用格式适合需要理解双曲守恒律、通量分裂、近似黎曼求解与高分辨率激波捕捉的读者对照学习。压缩包共46个m文件约22KB内容以MATLAB函数脚本为主按不同格式组织主程序、通量差分、空间离散、边界条件、物理量与守恒量转换、精确解计算、精度阶验证及结果绘图等模块便于比较各格式在Sod问题中的表现。已有3327人学习说明其在教学与练习场景中具有一定参考价值。读者可借助这些脚本复现激波管算例观察激波、接触间断与稀疏波的数值模拟结果并在此基础上修改初始条件、网格与格式参数理解数值通量、边界处理及误差阶验证的基本流程。各脚本之间接口相对清晰便于按格式拆分阅读与二次修改。1. 激波管求解不是“跑个曲线”这份一维 ENO Roe/Lax 代码能落地在哪做可压缩流课程设计时最怕拿到一份“能出图”的激波管代码却说不清密度平台为什么被抹平。这份“一维激波管问题的数值求解代码”把一维 Euler 方程、ENO 重构和 Roe/Lax 两类通量放在同一个可运行框架里适合复现 Sod 激波管、观察激波与接触间断也适合作为 CFD 入门到网格收敛验证的底稿。它不解决工程湍流问题但能把重构阶数、CFL、熵修正、边界条件这些参数落到一张密度曲线上。适合正在做课程设计、准备 CFD 基础题或需要一个可改的 Riemann 求解器对照实验的人。2. 先把物理和离散讲透Sod 激波管、ENO 重构与 Roe/Lax 的选型边界2.1 从 Euler 方程到激波管初值四个守恒量和状态方程一维激波管不是“画一条密度线”那么简单它验证的是可压缩 Euler 方程在间断附近的离散行为。守恒量通常写成 (U[\rho,\rho u,E]^T)通量 (F[\rho u,\rho u^2p,u(Ep)]^T)状态方程用理想气体 (p(\gamma-1)(E-\frac12\rho u^2))。这套代码的关键价值是把“物理初值、空间重构、数值通量、时间推进”拆成可替换的模块而不是把所有东西塞进一个主循环。我一般会先确认包内是否把原始变量和守恒量分开存放因为后面画图、算误差、调 Roe 平均都依赖这个边界。Sod 激波管的标准初值是计算域 (x\in[0,1])间断在 (x0.5)左态 (\rho1,u0,p1)右态 (\rho0.125,u0,p0.1)(\gamma1.4)输出时刻 (t0.2)。这个初值同时包含膨胀扇、接触间断和激波适合检查 ENO 重构是否在激波附近产生振荡也能看出 Roe 通量对接触间断的分辨能力。下面这段初始化代码可以直接对照包内入口参数确认它用的是原始变量还是守恒量。import numpy as np def sod_initial(x, gamma1.4): # x: 单元中心坐标长度 N rho np.where(x 0.5, 1.0, 0.125) u np.zeros_like(x) p np.where(x 0.5, 1.0, 0.1) # 总能量内能 动能 E p / (gamma - 1.0) 0.5 * rho * u * u # 守恒量按列存放便于后续通量计算 U np.stack([rho, rho * u, E], axis1) return U代码里的x 0.5是硬编码间断位置实际使用时可以改成参数gamma不要随手改成 1.0 或 2.0除非你明确知道状态方程和声速公式都同步改了。U的第二列是动量而不是速度第三列是总能量而不是压力很多“结果对不上”的翻车都出在画图时把守恒量当原始变量用。常见做法是再写一个cons_to_prim把 (\rho,u,p) 还原出来再做曲线。2.2 为什么用 ENO 而不是固定模板重构阶数、限制器与边界ENO 的核心不是“高阶一定更好”而是“在间断附近自动换模板”。固定模板在激波处容易把左右状态插到同一条跨越间断的多项式上密度曲线会出现非物理振荡ENO 通过比较候选模板的光滑程度选一个不跨间断的 stencil。二维/三维代码里还要做特征变量重构一维激波管虽然简单但如果包内直接在守恒量上做 ENO接触间断附近仍然可能被抹平。我一般会先确认重构变量原始变量重构实现简单特征变量重构对接触间断更友好但代码量更大。重构阶数常见取 ENO2 或 ENO3。ENO2 模板短、边界处理简单适合第一次跑通ENO3 在光滑区精度更高但对边界 ghost cell 和模板越界更敏感。下面这段只是 ENO2 模板选择的示意目的是让读者看懂“比较两侧差分、选更光滑的一侧”这个逻辑不要直接把它当生产代码。实际包内还要处理界面左右状态、边界外推和特征变量变换。import numpy as np def eno2_demo(v): # v: 单元平均量长度 N返回界面左右状态示意 N len(v) left np.zeros(N 1) right np.zeros(N 1) for i in range(1, N): # 两个候选斜率左偏和右偏 s_left v[i] - v[i - 1] s_right v[i 1] - v[i] if i 1 N else s_left # 选差分绝对值更小的一侧避免跨间断 if abs(s_right) abs(s_left): left[i] v[i] 0.5 * s_right right[i] v[i] - 0.5 * s_right else: left[i] v[i] 0.5 * s_left right[i] v[i] - 0.5 * s_left return left, right这段代码的参数说明很短v是单元平均left/right是界面两侧状态。真正要注意的是边界i0和iN附近不能直接套中心模板常见做法是降阶到一阶重构或者用零梯度 ghost cell。ENO 的“高阶”只在模板合法时成立边界处理错了激波后照样会起振。另一个坑是重构后没有做正性检查密度或压力出现负值后面 Roe 平均开根号直接出 NaN。2.3 时间推进与 CFL从半离散到 RK 的落地参数空间离散完之后时间推进常见用 RK2 或 RK3。显式格式的 CFL 条件大致是 (dt \mathrm{CFL}\cdot dx/\max(|u|a))其中 (a\sqrt{\gamma p/\rho}) 是声速。CFL 取 0.30.5 比较稳ENO3 或 Roe 熵修正较激进时可以再降。很多人第一次跑就把 CFL 拉到 0.9结果残差先震荡后发散回头怀疑 Roe 通量写错其实只是时间步越界。下面给一个 RK2 的骨架参数只有dt和rhs但rhs内部必须已经包含 ENO 重构、数值通量和边界条件。def rk2_step(U, dt, rhs): # U: 当前守恒量rhs(U) 返回 dU/dt k1 rhs(U) k2 rhs(U dt * k1) return U 0.5 * dt * (k1 k2)rhs的返回值形状要和U一致否则广播会把错误藏得很深。RK2 对一维激波管通常够用但如果包内用 ENO3 且时间步接近 CFL 上限建议换 RK3 或把 CFL 降到 0.3。参数表可以按下面这张表核对参数常用值作用调大后的典型后果CFL0.30.5控制时间步残差震荡、负压、NaNENO 阶数2 或 3空间重构精度边界处理不足时激波后起振网格数 N200/400/800分辨率接触间断被抹平误差表看不出阶数gamma1.4空气状态方程声速和激波速度同时偏掉输出时刻0.2Sod 标准对比激波位置对不上文献提示先固定gamma1.4、t0.2、x∈[0,1]这三个量再调网格和 CFL。三个基准量一变后面所有对比都失去意义。3. 把代码跑起来解压、编译、设参数到第一张密度曲线3.1 环境与依赖检查编译器、Python 绘图、输出格式拿到.rar后别急着双击运行先看目录结构和依赖。常见的一维激波管代码包会带 C/Fortran 源码、Makefile、输入参数文件和 README也有可能是 MATLAB 或 Python 脚本。第一步是解压并列出文件确认入口、编译方式、输出格式。下面命令只做环境检查不写死具体文件名实际以包内为准。# 解压到独立目录避免污染原目录 unzip 一维激波管问题的数值求解代码.rar -d shocktube cd shocktube # 先看文件类型和目录层级别急着编译 find . -maxdepth 2 -type f | sort file $(find . -maxdepth 2 -type f | head -n 5)find用来找 README、Makefile、输入文件和源码file用来判断是文本、压缩包还是可执行文件。如果包内是 C通常需要g或clang如果是 Fortran需要gfortran如果是 MATLAB确认.m文件的主脚本。Python 绘图依赖numpy和matplotlib没有的话先装好否则后面读结果只能看数字。输出格式常见是空白分隔的文本列顺序多为x rho u p或x rho u p e先head看一眼再写读取脚本。# 如果已有输出文件先确认列数和注释行 head -n 5 result.dat awk NR1{print NF} result.datNF是列数列数不对就不要继续画图。常见翻车是把x rho u p当成x rho u e压力曲线看起来“像那么回事”但激波后能量对不上。还有的包输出的是守恒量rho rho*u E需要先还原原始变量。3.2 入口参数怎么设网格数、CFL、重构阶数、Riemann 求解器入口参数一般有四类网格数、CFL、重构阶数、Riemann 求解器类型。如果程序支持命令行参数可以直接传如果只读输入文件就改同名变量。下面命令是模板shocktube要换成实际可执行文件名参数名也要按包内 README 调整。# 常见命令行入口模板没有命令行参数就去输入文件里改 ./shocktube --nx 400 --cfl 0.4 --tend 0.2 --recon ENO2 --flux Roe --out result.dat--nx 400表示 400 个单元第一次跑不建议直接上 1600先看 200 和 400 的差异--cfl 0.4是显式时间步系数Roe 通量加 ENO2 通常能稳住--tend 0.2是 Sod 标准输出时刻--recon ENO2和--flux Roe是两组可替换模块。如果包内参数名是N、CFL、T、scheme按实际改不要硬套。常见做法是先跑ENO2 Lax-Friedrichs确认框架能出图再换ENO2 Roe看接触间断是否更锐利最后才上 ENO3。3.3 一次完整运行编译、执行、读结果、画图编译时优先用包内 Makefile没有 Makefile 再用通配符编译所有源文件。下面命令不写死源文件名避免包内文件结构变化后命令失效。# 优先使用包内构建脚本 make -j4 # 没有 Makefile 时编译当前目录所有 C 源文件 g -O2 -stdc17 ./*.cpp -I./include -o shocktube-O2是常用优化级别调试阶段可以换-O0 -g-stdc17按包内代码要求调整-I./include只在有头文件目录时需要。编译报错先看第一个错误不要从最后一行倒着修。运行后得到结果文件用下面脚本读入并画密度和压力。import numpy as np import matplotlib.pyplot as plt # 把文件名改成包内实际输出常见列顺序x rho u p data np.loadtxt(result.dat, comments#) x, rho, u, p data[:, 0], data[:, 1], data[:, 2], data[:, 3] fig, ax plt.subplots(1, 2, figsize(10, 4)) ax[0].plot(x, rho, lw1.5, labelrho) ax[0].set_xlabel(x) ax[0].set_ylabel(density) ax[0].legend() ax[1].plot(x, p, lw1.5, labelp) ax[1].set_xlabel(x) ax[1].set_ylabel(pressure) ax[1].legend() fig.tight_layout() fig.savefig(sod_profile.png, dpi150)comments#用来跳过输出文件里的注释行data[:, 1]是密度列data[:, 3]是压力列如果包内列顺序不同改索引即可。画图前先确认x范围是不是 0 到 1间断是不是在 0.5 附近rho左平台是不是接近 1、右平台是不是接近 0.125。第一张图不要求漂亮先看有没有负值、NaN、明显对称破缺。如果密度曲线在激波后出现高频毛刺先降 CFL 或换 ENO2不要一上来就怀疑 Roe 通量。注意输出文件里的u如果是动量rho*u画速度曲线会得到一条量级完全不对的线。先看 README 或源码里的输出函数别靠猜。4. Roe 与 Lax-Friedrichs 怎么选通量函数、熵修正与交叉验证4.1 Roe 通量熵修正、特征值下限与激波稳定性Roe 通量的优点是接触间断分辨率好缺点是膨胀激波附近可能违反熵条件低速或近声速区域特征值接近零时会出现非物理解。常见做法是加 Harten 熵修正把绝对值小于delta的特征值抬起来。delta取 1e-6 到 1e-3 都有人用太小压不住膨胀激波太大又把接触间断抹平。下面是一个一维 Euler 的 Roe 通量骨架包含 Roe 平均、波强度和熵修正读者可以对照包内实现看差异。import numpy as np def cons_to_prim(U, gamma1.4): rho U[0] u U[1] / rho p (gamma - 1.0) * (U[2] - 0.5 * rho * u * u) return rho, u, p def euler_flux(U, gamma1.4): rho, u, p cons_to_prim(U, gamma) return np.array([rho * u, rho * u * u p, u * (U[2] p)]) def roe_flux(UL, UR, gamma1.4, delta1e-6): rhoL, uL, pL cons_to_prim(UL, gamma) rhoR, uR, pR cons_to_prim(UR, gamma) HL (UL[2] pL) / rhoL HR (UR[2] pR) / rhoR sqL, sqR np.sqrt(rhoL), np.sqrt(rhoR) # Roe 平均 u_tilde (sqL * uL sqR * uR) / (sqL sqR) H_tilde (sqL * HL sqR * HR) / (sqL sqR) rho_tilde sqL * sqR a_tilde np.sqrt(max((gamma - 1.0) * (H_tilde - 0.5 * u_tilde**2), delta)) # 特征值 lam np.array([u_tilde - a_tilde, u_tilde, u_tilde a_tilde]) abs_lam np.abs(lam) # Harten 熵修正近零特征值抬起来 for k in range(3): if abs_lam[k] delta: abs_lam[k] 0.5 * (lam[k] * lam[k] / delta delta) # 波强度 drho rhoR - rhoL du uR - uL dp pR - pL alpha np.array([ (dp - rho_tilde * a_tilde * du) / (2.0 * a_tilde * a_tilde), drho - dp / (a_tilde * a_tilde), (dp rho_tilde * a_tilde * du) / (2.0 * a_tilde * a_tilde) ]) # 右特征向量 R np.array([ [1.0, 1.0, 1.0], [u_tilde - a_tilde, u_tilde, u_tilde a_tilde], [H_tilde - u_tilde * a_tilde, 0.5 * u_tilde * u_tilde, H_tilde u_tilde * a_tilde] ]) diss R (abs_lam * alpha) FL euler_flux(UL, gamma) FR euler_flux(UR, gamma) return 0.5 * (FL FR) - 0.5 * dissdelta是熵修正阈值不是网格间距rho_tilde用左右密度几何平均不是算术平均alpha的三个分量对应三道波。这段代码只处理一个界面实际包内会对所有界面循环。Roe 通量跑 Sod 时接触间断通常比 Lax-Friedrichs 锐利但如果delta取太大接触间断会变宽看起来像“网格不够”。如果出现负压先把 CFL 降到 0.3再检查熵修正是否生效最后才怀疑重构越界。4.2 Lax-Friedrichs 通量数值耗散、CFL 与鲁棒性Lax-Friedrichs 通量更稳代价是数值耗散大接触间断容易被抹成一条斜坡。它适合第一次跑通框架、检查边界条件、对照 Roe 结果。局部 Lax-Friedrichs 的耗散系数取左右最大特征速度下面函数可以直接接到同一套界面循环里。def lax_friedrichs_flux(UL, UR, gamma1.4, alphaNone): FL euler_flux(UL, gamma) FR euler_flux(UR, gamma) if alpha is None: rhoL, uL, pL cons_to_prim(UL, gamma) rhoR, uR, pR cons_to_prim(UR, gamma) aL np.sqrt(gamma * pL / rhoL) aR np.sqrt(gamma * pR / rhoR) alpha max(abs(uL) aL, abs(uR) aR) # alpha 越大越稳也越耗散 return 0.5 * (FL FR) - 0.5 * alpha * (UR - UL)alpha是局部最大特征速度取全局最大值也能跑但接触间断会更糊。Lax-Friedrichs 的 CFL 仍然要小于 1常用 0.4 左右不要因为“它稳”就把 CFL 拉到 0.8。它和 Roe 的差别不在时间推进而在界面耗散同一套 ENO 重构下Lax 的密度曲线在接触间断处更宽激波位置可能略后移这是数值耗散的正常表现不是代码错。4.3 用同一组初值做交叉验证残差、激波位置、接触间断选 Roe 还是 Lax不要靠感觉用同一组初值、同一套网格、同一输出时刻做交叉验证。先跑ENO2 Lax确认程序能稳定出图再跑ENO2 Roe比较密度曲线在接触间断处的宽度、激波位置和压力平台。下面脚本假设两个结果文件网格一致、列顺序都是x rho u p算密度 L1 差异和激波附近的最大梯度位置。import numpy as np roe np.loadtxt(roe.dat) lax np.loadtxt(lax.dat) # 网格不一致时先插值到同一坐标 assert np.allclose(roe[:, 0], lax[:, 0]), 两套结果网格不一致先插值 l1_rho np.mean(np.abs(roe[:, 1] - lax[:, 1])) # 用最大密度梯度粗略定位激波 i_roe np.argmax(np.abs(np.gradient(roe[:, 1], roe[:, 0]))) i_lax np.argmax(np.abs(np.gradient(lax[:, 1], lax[:, 0]))) print(L1 density difference:, l1_rho) print(shock x Roe/Lax:, roe[i_roe, 0], lax[i_lax, 0])np.gradient用坐标和密度算梯度最大梯度位置可以粗略当激波位置L1差异不是越小越好而是要看差异集中在哪个区域。如果差异主要在接触间断说明 Roe 和 Lax 的耗散差异符合预期如果差异在膨胀扇内部检查熵修正和重构变量。交叉验证的目不是让两条曲线完全重合而是确认每个差异都能用通量耗散解释。5. 避坑与排查一维激波管代码最容易翻车的 5 个地方5.1 现象密度或压力出现负值、NaN现象很直接跑几十步后rho出现负数或者压力开根号得到 NaN。原因通常有三个CFL 太大导致时间步越界Roe 通量没有熵修正膨胀扇附近特征值过零ENO 重构在激波附近选到跨间断模板重构出负密度。解决顺序是先降 CFL 到 0.3再检查 Roe 的delta是否大于 0最后看重构后有没有做正性检查。常见做法是重构后如果rho 0或p 0把该界面降阶到一阶重构或者对原始变量做限幅。不要用“把负密度改成小正数”这种后悔药它会把质量问题藏到后面。5.2 现象激波后出现高频振荡激波后密度曲线出现细密毛刺网格越密越明显换 Lax 通量会好一些。原因多半是重构变量选错直接在守恒量上做高阶 ENO激波处能量和动量耦合容易产生非物理振荡另一个原因是边界 ghost cell 用了反射边界激波到边界后反射回来干扰内部。解决方法是先把 ENO 降到 2 阶确认毛刺是否消失再把重构变量换成特征变量或原始变量最后检查左右边界是否为零梯度/出流。如果毛刺只出现在边界附近优先改边界条件不要改通量。5.3 现象接触间断被抹得很平密度曲线上接触间断变成一条宽斜坡压力曲线却还比较锐利。原因通常是用了 Lax-Friedrichs 通量数值耗散大或者网格太粗200 个单元在接触间断处只有几个点还有可能是输出时刻不是 0.2接触间断已经移动到边界附近。解决方法是换 Roe 通量、把网格加到 400 或 800、确认t0.2并检查画图时用的是密度而不是动量。接触间断本身对耗散很敏感ENO2 Roe 在 400 网格下通常能看到明显平台但不要期待它像激波一样锐。5.4 现象CFL 稍微调大就发散CFL 从 0.4 调到 0.6 后残差迅速上升最后溢出。原因可能是时间推进用 RK2但空间重构是 ENO3稳定域不匹配也可能是边界条件把反射波放回来局部特征速度变大还有可能是 Roe 熵修正阈值太小近声速区域特征值接近零。解决方法是先固定 CFL0.3 跑通再以 0.05 为步长往上试ENO3 配 RK3或者 ENO2 配 RK2边界改成零梯度Roe 的delta从 1e-6 调到 1e-4 看是否改善。不要只看最后一步的 CFL显式格式的稳定域是全局的。5.5 现象结果和文献对不上密度曲线形状“差不多”但激波位置偏后、接触间断偏前或者左右平台值不对。原因通常是基准量错了x范围不是 0 到 1间断不在 0.5gamma不是 1.4输出时刻不是 0.2或者画图时把rho*u当成了u。解决方法是先打印入口参数和输出列名再核对初始条件用 200 网格先看左右平台左平台密度应接近 1右平台接近 0.125压力左平台接近 1右平台接近 0.1。如果平台值对但位置偏检查网格坐标是单元中心还是节点很多代码输出节点值画图时却按单元中心解释。6. 进阶验证与复用用网格收敛表把 Sod 结果变成可对比的底稿跑出一张漂亮的密度曲线只是开始真正让这份代码有复用价值的是网格收敛表。我一般固定ENO2 Roe跑 200、400、800 三套网格再用 1600 网格结果作为参考把粗网格插值到参考网格坐标上算密度 L1 误差。下面脚本可以直接改文件名使用interp只对密度列做插值避免把激波位置抹平。import numpy as np def l1_against_reference(ref_file, cur_file): ref np.loadtxt(ref_file) cur np.loadtxt(cur_file) # 把参考解插值到当前网格坐标只比较密度 rho_ref np.interp(cur[:, 0], ref[:, 0], ref[:, 1]) return np.mean(np.abs(cur[:, 1] - rho_ref)) ref_file sod_1600.dat for name in [sod_200.dat, sod_400.dat, sod_800.dat]: print(name, l1_against_reference(ref_file, name))把结果填成下面这张表比单看曲线更有说服力。表格里的“待填”不是偷懒而是让你用自己的输出替换避免拿别人的数字当结论。网格 N密度 L1 误差激波位置接触间断中心备注200待填待填待填先看趋势400待填待填待填误差应下降800待填待填待填接触间断变窄如果 400 到 800 的 L1 误差没有明显下降先别怀疑 ENO 阶数检查输出时刻和参考解插值方式如果激波位置随网格加密反而漂移检查边界条件和 CFL 是否随网格同步调整。验证完 Sod 后可以把初值换成左右密度比更大、压力比更小的激波管变体观察 Roe 熵修正和 ENO 模板选择是否还稳。从那以后我每次拿到这类一维激波管代码都先跑 200 网格的 Sod再补 400 和 800 的误差表最后才改初值和通量这套顺序能少走很多“曲线像但数不对”的弯路。希望帮到你。本文还有配套的精品资源点击获取 一维激波管ENO重构与Roe/Lax通量求解代码实践
简介面向计算流体力学课程学习与数值方法入门这套MATLAB代码围绕一维激波管问题给出可运行的求解示例覆盖Lax-Friedrichs、Lax-Wendroff、Roe与ENO等常用格式适合需要理解双曲守恒律、通量分裂、近似黎曼求解与高分辨率激波捕捉的读者对照学习。压缩包共46个m文件约22KB内容以MATLAB函数脚本为主按不同格式组织主程序、通量差分、空间离散、边界条件、物理量与守恒量转换、精确解计算、精度阶验证及结果绘图等模块便于比较各格式在Sod问题中的表现。已有3327人学习说明其在教学与练习场景中具有一定参考价值。读者可借助这些脚本复现激波管算例观察激波、接触间断与稀疏波的数值模拟结果并在此基础上修改初始条件、网格与格式参数理解数值通量、边界处理及误差阶验证的基本流程。各脚本之间接口相对清晰便于按格式拆分阅读与二次修改。1. 激波管求解不是“跑个曲线”这份一维 ENO Roe/Lax 代码能落地在哪做可压缩流课程设计时最怕拿到一份“能出图”的激波管代码却说不清密度平台为什么被抹平。这份“一维激波管问题的数值求解代码”把一维 Euler 方程、ENO 重构和 Roe/Lax 两类通量放在同一个可运行框架里适合复现 Sod 激波管、观察激波与接触间断也适合作为 CFD 入门到网格收敛验证的底稿。它不解决工程湍流问题但能把重构阶数、CFL、熵修正、边界条件这些参数落到一张密度曲线上。适合正在做课程设计、准备 CFD 基础题或需要一个可改的 Riemann 求解器对照实验的人。2. 先把物理和离散讲透Sod 激波管、ENO 重构与 Roe/Lax 的选型边界2.1 从 Euler 方程到激波管初值四个守恒量和状态方程一维激波管不是“画一条密度线”那么简单它验证的是可压缩 Euler 方程在间断附近的离散行为。守恒量通常写成 (U[\rho,\rho u,E]^T)通量 (F[\rho u,\rho u^2p,u(Ep)]^T)状态方程用理想气体 (p(\gamma-1)(E-\frac12\rho u^2))。这套代码的关键价值是把“物理初值、空间重构、数值通量、时间推进”拆成可替换的模块而不是把所有东西塞进一个主循环。我一般会先确认包内是否把原始变量和守恒量分开存放因为后面画图、算误差、调 Roe 平均都依赖这个边界。Sod 激波管的标准初值是计算域 (x\in[0,1])间断在 (x0.5)左态 (\rho1,u0,p1)右态 (\rho0.125,u0,p0.1)(\gamma1.4)输出时刻 (t0.2)。这个初值同时包含膨胀扇、接触间断和激波适合检查 ENO 重构是否在激波附近产生振荡也能看出 Roe 通量对接触间断的分辨能力。下面这段初始化代码可以直接对照包内入口参数确认它用的是原始变量还是守恒量。import numpy as np def sod_initial(x, gamma1.4): # x: 单元中心坐标长度 N rho np.where(x 0.5, 1.0, 0.125) u np.zeros_like(x) p np.where(x 0.5, 1.0, 0.1) # 总能量内能 动能 E p / (gamma - 1.0) 0.5 * rho * u * u # 守恒量按列存放便于后续通量计算 U np.stack([rho, rho * u, E], axis1) return U代码里的x 0.5是硬编码间断位置实际使用时可以改成参数gamma不要随手改成 1.0 或 2.0除非你明确知道状态方程和声速公式都同步改了。U的第二列是动量而不是速度第三列是总能量而不是压力很多“结果对不上”的翻车都出在画图时把守恒量当原始变量用。常见做法是再写一个cons_to_prim把 (\rho,u,p) 还原出来再做曲线。2.2 为什么用 ENO 而不是固定模板重构阶数、限制器与边界ENO 的核心不是“高阶一定更好”而是“在间断附近自动换模板”。固定模板在激波处容易把左右状态插到同一条跨越间断的多项式上密度曲线会出现非物理振荡ENO 通过比较候选模板的光滑程度选一个不跨间断的 stencil。二维/三维代码里还要做特征变量重构一维激波管虽然简单但如果包内直接在守恒量上做 ENO接触间断附近仍然可能被抹平。我一般会先确认重构变量原始变量重构实现简单特征变量重构对接触间断更友好但代码量更大。重构阶数常见取 ENO2 或 ENO3。ENO2 模板短、边界处理简单适合第一次跑通ENO3 在光滑区精度更高但对边界 ghost cell 和模板越界更敏感。下面这段只是 ENO2 模板选择的示意目的是让读者看懂“比较两侧差分、选更光滑的一侧”这个逻辑不要直接把它当生产代码。实际包内还要处理界面左右状态、边界外推和特征变量变换。import numpy as np def eno2_demo(v): # v: 单元平均量长度 N返回界面左右状态示意 N len(v) left np.zeros(N 1) right np.zeros(N 1) for i in range(1, N): # 两个候选斜率左偏和右偏 s_left v[i] - v[i - 1] s_right v[i 1] - v[i] if i 1 N else s_left # 选差分绝对值更小的一侧避免跨间断 if abs(s_right) abs(s_left): left[i] v[i] 0.5 * s_right right[i] v[i] - 0.5 * s_right else: left[i] v[i] 0.5 * s_left right[i] v[i] - 0.5 * s_left return left, right这段代码的参数说明很短v是单元平均left/right是界面两侧状态。真正要注意的是边界i0和iN附近不能直接套中心模板常见做法是降阶到一阶重构或者用零梯度 ghost cell。ENO 的“高阶”只在模板合法时成立边界处理错了激波后照样会起振。另一个坑是重构后没有做正性检查密度或压力出现负值后面 Roe 平均开根号直接出 NaN。2.3 时间推进与 CFL从半离散到 RK 的落地参数空间离散完之后时间推进常见用 RK2 或 RK3。显式格式的 CFL 条件大致是 (dt \mathrm{CFL}\cdot dx/\max(|u|a))其中 (a\sqrt{\gamma p/\rho}) 是声速。CFL 取 0.30.5 比较稳ENO3 或 Roe 熵修正较激进时可以再降。很多人第一次跑就把 CFL 拉到 0.9结果残差先震荡后发散回头怀疑 Roe 通量写错其实只是时间步越界。下面给一个 RK2 的骨架参数只有dt和rhs但rhs内部必须已经包含 ENO 重构、数值通量和边界条件。def rk2_step(U, dt, rhs): # U: 当前守恒量rhs(U) 返回 dU/dt k1 rhs(U) k2 rhs(U dt * k1) return U 0.5 * dt * (k1 k2)rhs的返回值形状要和U一致否则广播会把错误藏得很深。RK2 对一维激波管通常够用但如果包内用 ENO3 且时间步接近 CFL 上限建议换 RK3 或把 CFL 降到 0.3。参数表可以按下面这张表核对参数常用值作用调大后的典型后果CFL0.30.5控制时间步残差震荡、负压、NaNENO 阶数2 或 3空间重构精度边界处理不足时激波后起振网格数 N200/400/800分辨率接触间断被抹平误差表看不出阶数gamma1.4空气状态方程声速和激波速度同时偏掉输出时刻0.2Sod 标准对比激波位置对不上文献提示先固定gamma1.4、t0.2、x∈[0,1]这三个量再调网格和 CFL。三个基准量一变后面所有对比都失去意义。3. 把代码跑起来解压、编译、设参数到第一张密度曲线3.1 环境与依赖检查编译器、Python 绘图、输出格式拿到.rar后别急着双击运行先看目录结构和依赖。常见的一维激波管代码包会带 C/Fortran 源码、Makefile、输入参数文件和 README也有可能是 MATLAB 或 Python 脚本。第一步是解压并列出文件确认入口、编译方式、输出格式。下面命令只做环境检查不写死具体文件名实际以包内为准。# 解压到独立目录避免污染原目录 unzip 一维激波管问题的数值求解代码.rar -d shocktube cd shocktube # 先看文件类型和目录层级别急着编译 find . -maxdepth 2 -type f | sort file $(find . -maxdepth 2 -type f | head -n 5)find用来找 README、Makefile、输入文件和源码file用来判断是文本、压缩包还是可执行文件。如果包内是 C通常需要g或clang如果是 Fortran需要gfortran如果是 MATLAB确认.m文件的主脚本。Python 绘图依赖numpy和matplotlib没有的话先装好否则后面读结果只能看数字。输出格式常见是空白分隔的文本列顺序多为x rho u p或x rho u p e先head看一眼再写读取脚本。# 如果已有输出文件先确认列数和注释行 head -n 5 result.dat awk NR1{print NF} result.datNF是列数列数不对就不要继续画图。常见翻车是把x rho u p当成x rho u e压力曲线看起来“像那么回事”但激波后能量对不上。还有的包输出的是守恒量rho rho*u E需要先还原原始变量。3.2 入口参数怎么设网格数、CFL、重构阶数、Riemann 求解器入口参数一般有四类网格数、CFL、重构阶数、Riemann 求解器类型。如果程序支持命令行参数可以直接传如果只读输入文件就改同名变量。下面命令是模板shocktube要换成实际可执行文件名参数名也要按包内 README 调整。# 常见命令行入口模板没有命令行参数就去输入文件里改 ./shocktube --nx 400 --cfl 0.4 --tend 0.2 --recon ENO2 --flux Roe --out result.dat--nx 400表示 400 个单元第一次跑不建议直接上 1600先看 200 和 400 的差异--cfl 0.4是显式时间步系数Roe 通量加 ENO2 通常能稳住--tend 0.2是 Sod 标准输出时刻--recon ENO2和--flux Roe是两组可替换模块。如果包内参数名是N、CFL、T、scheme按实际改不要硬套。常见做法是先跑ENO2 Lax-Friedrichs确认框架能出图再换ENO2 Roe看接触间断是否更锐利最后才上 ENO3。3.3 一次完整运行编译、执行、读结果、画图编译时优先用包内 Makefile没有 Makefile 再用通配符编译所有源文件。下面命令不写死源文件名避免包内文件结构变化后命令失效。# 优先使用包内构建脚本 make -j4 # 没有 Makefile 时编译当前目录所有 C 源文件 g -O2 -stdc17 ./*.cpp -I./include -o shocktube-O2是常用优化级别调试阶段可以换-O0 -g-stdc17按包内代码要求调整-I./include只在有头文件目录时需要。编译报错先看第一个错误不要从最后一行倒着修。运行后得到结果文件用下面脚本读入并画密度和压力。import numpy as np import matplotlib.pyplot as plt # 把文件名改成包内实际输出常见列顺序x rho u p data np.loadtxt(result.dat, comments#) x, rho, u, p data[:, 0], data[:, 1], data[:, 2], data[:, 3] fig, ax plt.subplots(1, 2, figsize(10, 4)) ax[0].plot(x, rho, lw1.5, labelrho) ax[0].set_xlabel(x) ax[0].set_ylabel(density) ax[0].legend() ax[1].plot(x, p, lw1.5, labelp) ax[1].set_xlabel(x) ax[1].set_ylabel(pressure) ax[1].legend() fig.tight_layout() fig.savefig(sod_profile.png, dpi150)comments#用来跳过输出文件里的注释行data[:, 1]是密度列data[:, 3]是压力列如果包内列顺序不同改索引即可。画图前先确认x范围是不是 0 到 1间断是不是在 0.5 附近rho左平台是不是接近 1、右平台是不是接近 0.125。第一张图不要求漂亮先看有没有负值、NaN、明显对称破缺。如果密度曲线在激波后出现高频毛刺先降 CFL 或换 ENO2不要一上来就怀疑 Roe 通量。注意输出文件里的u如果是动量rho*u画速度曲线会得到一条量级完全不对的线。先看 README 或源码里的输出函数别靠猜。4. Roe 与 Lax-Friedrichs 怎么选通量函数、熵修正与交叉验证4.1 Roe 通量熵修正、特征值下限与激波稳定性Roe 通量的优点是接触间断分辨率好缺点是膨胀激波附近可能违反熵条件低速或近声速区域特征值接近零时会出现非物理解。常见做法是加 Harten 熵修正把绝对值小于delta的特征值抬起来。delta取 1e-6 到 1e-3 都有人用太小压不住膨胀激波太大又把接触间断抹平。下面是一个一维 Euler 的 Roe 通量骨架包含 Roe 平均、波强度和熵修正读者可以对照包内实现看差异。import numpy as np def cons_to_prim(U, gamma1.4): rho U[0] u U[1] / rho p (gamma - 1.0) * (U[2] - 0.5 * rho * u * u) return rho, u, p def euler_flux(U, gamma1.4): rho, u, p cons_to_prim(U, gamma) return np.array([rho * u, rho * u * u p, u * (U[2] p)]) def roe_flux(UL, UR, gamma1.4, delta1e-6): rhoL, uL, pL cons_to_prim(UL, gamma) rhoR, uR, pR cons_to_prim(UR, gamma) HL (UL[2] pL) / rhoL HR (UR[2] pR) / rhoR sqL, sqR np.sqrt(rhoL), np.sqrt(rhoR) # Roe 平均 u_tilde (sqL * uL sqR * uR) / (sqL sqR) H_tilde (sqL * HL sqR * HR) / (sqL sqR) rho_tilde sqL * sqR a_tilde np.sqrt(max((gamma - 1.0) * (H_tilde - 0.5 * u_tilde**2), delta)) # 特征值 lam np.array([u_tilde - a_tilde, u_tilde, u_tilde a_tilde]) abs_lam np.abs(lam) # Harten 熵修正近零特征值抬起来 for k in range(3): if abs_lam[k] delta: abs_lam[k] 0.5 * (lam[k] * lam[k] / delta delta) # 波强度 drho rhoR - rhoL du uR - uL dp pR - pL alpha np.array([ (dp - rho_tilde * a_tilde * du) / (2.0 * a_tilde * a_tilde), drho - dp / (a_tilde * a_tilde), (dp rho_tilde * a_tilde * du) / (2.0 * a_tilde * a_tilde) ]) # 右特征向量 R np.array([ [1.0, 1.0, 1.0], [u_tilde - a_tilde, u_tilde, u_tilde a_tilde], [H_tilde - u_tilde * a_tilde, 0.5 * u_tilde * u_tilde, H_tilde u_tilde * a_tilde] ]) diss R (abs_lam * alpha) FL euler_flux(UL, gamma) FR euler_flux(UR, gamma) return 0.5 * (FL FR) - 0.5 * dissdelta是熵修正阈值不是网格间距rho_tilde用左右密度几何平均不是算术平均alpha的三个分量对应三道波。这段代码只处理一个界面实际包内会对所有界面循环。Roe 通量跑 Sod 时接触间断通常比 Lax-Friedrichs 锐利但如果delta取太大接触间断会变宽看起来像“网格不够”。如果出现负压先把 CFL 降到 0.3再检查熵修正是否生效最后才怀疑重构越界。4.2 Lax-Friedrichs 通量数值耗散、CFL 与鲁棒性Lax-Friedrichs 通量更稳代价是数值耗散大接触间断容易被抹成一条斜坡。它适合第一次跑通框架、检查边界条件、对照 Roe 结果。局部 Lax-Friedrichs 的耗散系数取左右最大特征速度下面函数可以直接接到同一套界面循环里。def lax_friedrichs_flux(UL, UR, gamma1.4, alphaNone): FL euler_flux(UL, gamma) FR euler_flux(UR, gamma) if alpha is None: rhoL, uL, pL cons_to_prim(UL, gamma) rhoR, uR, pR cons_to_prim(UR, gamma) aL np.sqrt(gamma * pL / rhoL) aR np.sqrt(gamma * pR / rhoR) alpha max(abs(uL) aL, abs(uR) aR) # alpha 越大越稳也越耗散 return 0.5 * (FL FR) - 0.5 * alpha * (UR - UL)alpha是局部最大特征速度取全局最大值也能跑但接触间断会更糊。Lax-Friedrichs 的 CFL 仍然要小于 1常用 0.4 左右不要因为“它稳”就把 CFL 拉到 0.8。它和 Roe 的差别不在时间推进而在界面耗散同一套 ENO 重构下Lax 的密度曲线在接触间断处更宽激波位置可能略后移这是数值耗散的正常表现不是代码错。4.3 用同一组初值做交叉验证残差、激波位置、接触间断选 Roe 还是 Lax不要靠感觉用同一组初值、同一套网格、同一输出时刻做交叉验证。先跑ENO2 Lax确认程序能稳定出图再跑ENO2 Roe比较密度曲线在接触间断处的宽度、激波位置和压力平台。下面脚本假设两个结果文件网格一致、列顺序都是x rho u p算密度 L1 差异和激波附近的最大梯度位置。import numpy as np roe np.loadtxt(roe.dat) lax np.loadtxt(lax.dat) # 网格不一致时先插值到同一坐标 assert np.allclose(roe[:, 0], lax[:, 0]), 两套结果网格不一致先插值 l1_rho np.mean(np.abs(roe[:, 1] - lax[:, 1])) # 用最大密度梯度粗略定位激波 i_roe np.argmax(np.abs(np.gradient(roe[:, 1], roe[:, 0]))) i_lax np.argmax(np.abs(np.gradient(lax[:, 1], lax[:, 0]))) print(L1 density difference:, l1_rho) print(shock x Roe/Lax:, roe[i_roe, 0], lax[i_lax, 0])np.gradient用坐标和密度算梯度最大梯度位置可以粗略当激波位置L1差异不是越小越好而是要看差异集中在哪个区域。如果差异主要在接触间断说明 Roe 和 Lax 的耗散差异符合预期如果差异在膨胀扇内部检查熵修正和重构变量。交叉验证的目不是让两条曲线完全重合而是确认每个差异都能用通量耗散解释。5. 避坑与排查一维激波管代码最容易翻车的 5 个地方5.1 现象密度或压力出现负值、NaN现象很直接跑几十步后rho出现负数或者压力开根号得到 NaN。原因通常有三个CFL 太大导致时间步越界Roe 通量没有熵修正膨胀扇附近特征值过零ENO 重构在激波附近选到跨间断模板重构出负密度。解决顺序是先降 CFL 到 0.3再检查 Roe 的delta是否大于 0最后看重构后有没有做正性检查。常见做法是重构后如果rho 0或p 0把该界面降阶到一阶重构或者对原始变量做限幅。不要用“把负密度改成小正数”这种后悔药它会把质量问题藏到后面。5.2 现象激波后出现高频振荡激波后密度曲线出现细密毛刺网格越密越明显换 Lax 通量会好一些。原因多半是重构变量选错直接在守恒量上做高阶 ENO激波处能量和动量耦合容易产生非物理振荡另一个原因是边界 ghost cell 用了反射边界激波到边界后反射回来干扰内部。解决方法是先把 ENO 降到 2 阶确认毛刺是否消失再把重构变量换成特征变量或原始变量最后检查左右边界是否为零梯度/出流。如果毛刺只出现在边界附近优先改边界条件不要改通量。5.3 现象接触间断被抹得很平密度曲线上接触间断变成一条宽斜坡压力曲线却还比较锐利。原因通常是用了 Lax-Friedrichs 通量数值耗散大或者网格太粗200 个单元在接触间断处只有几个点还有可能是输出时刻不是 0.2接触间断已经移动到边界附近。解决方法是换 Roe 通量、把网格加到 400 或 800、确认t0.2并检查画图时用的是密度而不是动量。接触间断本身对耗散很敏感ENO2 Roe 在 400 网格下通常能看到明显平台但不要期待它像激波一样锐。5.4 现象CFL 稍微调大就发散CFL 从 0.4 调到 0.6 后残差迅速上升最后溢出。原因可能是时间推进用 RK2但空间重构是 ENO3稳定域不匹配也可能是边界条件把反射波放回来局部特征速度变大还有可能是 Roe 熵修正阈值太小近声速区域特征值接近零。解决方法是先固定 CFL0.3 跑通再以 0.05 为步长往上试ENO3 配 RK3或者 ENO2 配 RK2边界改成零梯度Roe 的delta从 1e-6 调到 1e-4 看是否改善。不要只看最后一步的 CFL显式格式的稳定域是全局的。5.5 现象结果和文献对不上密度曲线形状“差不多”但激波位置偏后、接触间断偏前或者左右平台值不对。原因通常是基准量错了x范围不是 0 到 1间断不在 0.5gamma不是 1.4输出时刻不是 0.2或者画图时把rho*u当成了u。解决方法是先打印入口参数和输出列名再核对初始条件用 200 网格先看左右平台左平台密度应接近 1右平台接近 0.125压力左平台接近 1右平台接近 0.1。如果平台值对但位置偏检查网格坐标是单元中心还是节点很多代码输出节点值画图时却按单元中心解释。6. 进阶验证与复用用网格收敛表把 Sod 结果变成可对比的底稿跑出一张漂亮的密度曲线只是开始真正让这份代码有复用价值的是网格收敛表。我一般固定ENO2 Roe跑 200、400、800 三套网格再用 1600 网格结果作为参考把粗网格插值到参考网格坐标上算密度 L1 误差。下面脚本可以直接改文件名使用interp只对密度列做插值避免把激波位置抹平。import numpy as np def l1_against_reference(ref_file, cur_file): ref np.loadtxt(ref_file) cur np.loadtxt(cur_file) # 把参考解插值到当前网格坐标只比较密度 rho_ref np.interp(cur[:, 0], ref[:, 0], ref[:, 1]) return np.mean(np.abs(cur[:, 1] - rho_ref)) ref_file sod_1600.dat for name in [sod_200.dat, sod_400.dat, sod_800.dat]: print(name, l1_against_reference(ref_file, name))把结果填成下面这张表比单看曲线更有说服力。表格里的“待填”不是偷懒而是让你用自己的输出替换避免拿别人的数字当结论。网格 N密度 L1 误差激波位置接触间断中心备注200待填待填待填先看趋势400待填待填待填误差应下降800待填待填待填接触间断变窄如果 400 到 800 的 L1 误差没有明显下降先别怀疑 ENO 阶数检查输出时刻和参考解插值方式如果激波位置随网格加密反而漂移检查边界条件和 CFL 是否随网格同步调整。验证完 Sod 后可以把初值换成左右密度比更大、压力比更小的激波管变体观察 Roe 熵修正和 ENO 模板选择是否还稳。从那以后我每次拿到这类一维激波管代码都先跑 200 网格的 Sod再补 400 和 800 的误差表最后才改初值和通量这套顺序能少走很多“曲线像但数不对”的弯路。希望帮到你。本文还有配套的精品资源点击获取