一维激波管CFD数值求解：Roe通量、ENO重构与验证

简介这份MATLAB程序包面向学习计算流体力学与双曲守恒律数值方法的本科生、研究生及科研人员围绕一维激波管问题给出可运行的求解框架帮助理解Lax-Friedrichs、Lax-Wendroff、Roe与ENO等格式的构造差异及激波捕捉效果适合课程教学与数值方法入门。压缩包共46个文件均为.m脚本约22KB按不同数值格式组织主程序、通量差分、空间离散、边界条件、物理量守恒量转换、精确解与误差精度分析、动态绘图等模块整体结构紧凑便于逐项对照阅读与二次修改。已有3327人学习说明其在教学与入门实践中具有一定参考价值。读者可借此完成Sod激波管算例的编程实现观察不同格式对间断的分辨能力并利用内置绘图与精确解模块进行结果对比、收敛性检查与排错适合作为课程作业、数值实验或算法复现的起点。1. 一行初值、四道波这个激波管小项目到底在算什么一维激波管问题说白了就是给一根管子左右两段气体各设一套密度、压力、速度中间用一张膜隔开膜一破气体就开始互相挤压然后你会看到向右跑的激波、跟在后面的接触间断、向左退的稀疏波和一道声波。这套解析解是少有的能写出闭式表达式、又能把数值格式的问题全都暴露出来的算例所以做 ENO、Lax、Roe 这些格式验证时几乎绕不开它。这份一维激波管问题的数值求解代码包干的正是这件事把一维欧拉方程离散到有限体积框架上用 Roe 近似 Riemann 求解器和 Lax-Friedrichs 通量做对照再用高阶 ENO 重构把界面状态抬到三阶精度最后把密度、压力、速度曲线和解析解画在一起。适合刚接触 CFD、想验证自己写的通量函数对不对的人也适合想拿一份干净骨架改成二维、加限制器、换时间格式的熟手。2. 项目骨架与数值选型一维欧拉方程怎么落到有限体积框架上2.1 目录结构与数据流这类激波管代码的工程结构通常不会太复杂但如果组织得干净二次开发会省很多事。我拿到的这份代码大体上分成四层参数与常量层、网格与初值层、数值格式层、输出层。各层之间只通过结构体传数据没有全局变量乱飞这一点在课程设计级别的代码里算是比较讲究的。打开目录常见的组织方式是这样的shocktube/ ├── include/ │ ├── params.h # 全局参数gamma、CFL、网格数、格式开关 │ ├── euler.h # 原始变量与守恒量互转、通量计算 │ ├── flux.h # Roe / Lax-Friedrichs 通量函数声明 │ ├── recon.h # ENO 重构接口 │ └── solver.h # 时间推进主循环 ├── src/ │ ├── main.cpp # 读参数、建网格、跑循环、写文件 │ ├── euler.cpp # 状态方程与特征值 │ ├── flux_roe.cpp # Roe 平均与通量差分分裂 │ ├── flux_lf.cpp # Lax-Friedrichs 通量 │ ├── recon_eno.cpp # ENO 模板选择与多项式重构 │ └── solver.cpp # RK3 时间推进 ├── config/ │ └── sod.cfg # Sod 激波管初值、格式选择、结束时间 └── output/ └── (运行时生成 .dat 文件)数据流是单向的main.cpp读sod.cfg填充Params结构体按网格数生成单元中心坐标调用初值函数把左右两段状态铺上去然后进入solver.cpp的时间循环。循环内部每一步做四件事ENO 重构得到左右界面状态、调用通量函数得到数值通量、用通量差更新守恒量、按 CFL 条件重算时间步。循环结束后把所有单元的中心值按列写到.dat画图脚本再去读。提示如果目录里没有单独的config/参数大概率硬编码在main.cpp顶部或一个params.h里第一次跑之前先翻一遍这两个地方把网格数、CFL、结束时间、格式开关找出来。这套结构的好处是格式可替换。flux.h里声明的是同一套函数签名Roe 和 Lax-Friedrichs 各有各的实现主循环里用一个枚举或者字符串开关切换。想加一个新格式只要实现同样签名的函数在开关里补一个分支就行不用动时间推进的代码。我后来在别的项目里沿用这个思路把 HLLC、AUSM 也挂上去改动量确实只有几十行。2.2 控制方程、状态方程与无量纲化一维欧拉方程写成守恒形式是∂U/∂t ∂F(U)/∂x 0守恒量 U 和通量 F 分别是U [ρ, ρu, E]ᵀF [ρu, ρu² p, u(E p)]ᵀ其中 E 是单位体积总能量理想气体状态下 E p/(γ−1) ½ρu²。γ 取 1.4 是空气的常规值代码里一般写成const double gamma 1.4;或者从配置文件读。这里有个新手容易忽略的点绝大多数激波管代码会做无量纲化处理。管长取 1初始密度和压力取 1 和 1 的量级声速 a √(γp/ρ) 也就是 1 左右特征时间大约就是 1。这样一来结束时间设成 0.2 左右就能看到完整的波系结构CFL 数取 0.40.8 都算合理范围。如果不做无量纲化把管长设成真实的 1 米、压力设成 101325 帕量级差出好几个数量级浮点误差和 CFL 判断都会变得很难看。我一般会在params.h里显式写清楚量纲约定哪怕只是一行注释。血泪经验是过几个月再回来看忘了当时用的什么单位一个dt算出来小得离谱排查半天才发现是声速的量级不对。2.3 三种通量格式的分工Roe、Lax-Friedrichs 与 ENO 重构这份代码里三个名词不是并列关系而是各管一段。Roe 格式负责的是界面通量。它把左右两个状态的差分解到三个特征方向上每个方向按特征速度的符号决定信息往哪边传。优点是接触间断和激波的分辨率高缺点是熵修正不做的话在声速点附近会出现非物理的膨胀激波这也是后面避坑章节要重点说的。Lax-Friedrichs 格式负责的是对照。它的形式简单粗暴取左右通量的平均再减去一个正比于状态差的数值耗散项耗散系数取界面上最大的特征速度。稳定、好实现、几乎不会崩代价是把接触间断抹得很厚。放在代码里主要是给 Roe 当参照物跑一遍对比图就能直观看出低耗散格式和高耗散格式的差别。ENO 重构负责的是空间精度。它不直接改通量而是在算通量之前用几个相邻单元的守恒量拟合一个多项式取多项式在界面处的值作为左右状态。三阶 ENO 会在两个候选模板里比较差商大小选那个更光滑的避免跨过激波做插值产生振荡。三者串起来的顺序是ENO 重构 → 得到界面左右状态 → Roe 或 LF 算通量 → RK3 推进。3. 核心模块拆解一个时间步里从 Roe 平均到 RK3 更新的完整链路3.1 网格、初值与边界条件网格用的是单元中心有限体积第 i 个单元的中心坐标 x_i (i 0.5)Δx界面在 iΔx 处。Sod 问题的初值很标准膜放在 x 0.5左边是 (ρ, u, p) (1.0, 0.0, 1.0)右边是 (0.125, 0.0, 0.1)。// 按单元中心坐标铺初值N 为单元总数x 为单元中心数组 void initSod(int N, const double* x, double gamma, double* rho, double* u, double* p) { for (int i 0; i N; i) { if (x[i] 0.5) { // 高压段 rho[i] 1.0; u[i] 0.0; p[i] 1.0; } else { // 低压段 rho[i] 0.125; u[i] 0.0; p[i] 0.1; } } }这段代码里没有做任何平滑过渡界面处密度从 1 直接跳到 0.125。有人会担心这么大的跳变会不会一上来就崩其实不会因为数值通量本身带耗散第一、二步会把跳变抹成几格宽的过渡带之后才是格式分辨率的真正考验。参数 gamma 在这里没用到但保留在签名里方便以后换气体。边界条件一般用两种固定远场值的 Dirichlet 边界或者零梯度外推。激波管两端在计算时间内波还没跑到用外推就够写法是把边界单元的守恒量复制到虚拟单元。注意外推要作用在守恒量上而不是原始变量上先外推再转原始变量顺序反了在低压区可能出现负密度。3.2 Roe 平均与通量差分分裂的实现细节Roe 格式的核心是那个加权平均状态。用 √ρ 做权重把左右状态的 u 和总焓 H 平均出来再由此算出平均声速 ã。// Roe 平均输入左右原始变量输出平均速度和平均声速 void roeAverage(double rhoL, double uL, double HL, double rhoR, double uR, double HR, double gamma, double uTilde, double aTilde) { double sL std::sqrt(rhoL); double sR std::sqrt(rhoR); double den sL sR; uTilde (sL * uL sR * uR) / den; double HTilde (sL * HL sR * HR) / den; aTilde std::sqrt((gamma - 1.0) * (HTilde - 0.5 * uTilde * uTilde)); }这里 HL、HR 是左右状态的总焓 (E p)/ρ。注意 ã² 的表达式不是直接对声速做加权平均而是从平均焓里减掉动能的一半再乘 (γ−1)这个形式能保证 Roe 格式满足 U 和 F 的跳跃条件是 Roe 格式具有精确捕捉接触间断能力的根源。拿到平均量和特征值 λ₁ ũ−ãλ₂ ũλ₃ ũã 之后要把左右守恒量之差 ΔU 投影到三个特征向量上得到波强度 α₁、α₂、α₃最后组装数值通量// Roe 数值通量FL/FR 为左右物理通量UL/UR 为左右守恒量 void roeFlux(const double* UL, const double* UR, const double* FL, const double* FR, double gamma, double* Fnum) { double rhoL UL[0], rhoR UR[0]; double uL UL[1] / rhoL, uR UR[1] / rhoR; double EL UL[2], ER UR[2]; double pL (gamma - 1.0) * (EL - 0.5 * rhoL * uL * uL); double pR (gamma - 1.0) * (ER - 0.5 * rhoR * uR * uR); double HL (EL pL) / rhoL, HR (ER pR) / rhoR; double uT, aT; roeAverage(rhoL, uL, HL, rhoR, uR, HR, gamma, uT, aT); // 三个特征值这里省掉了熵修正实际使用要补 double lam1 uT - aT, lam2 uT, lam3 uT aT; double drho rhoR - rhoL; double drhou UR[1] - UL[1]; double dE UR[2] - UL[2]; // 波强度把 ΔU 投影到特征向量上 double alpha2 drho - dp_over_a2(drho, drhou, dE, uT, aT, gamma); double alpha1 (dp_over_a2(drho, drhou, dE, uT, aT, gamma) - rhoT * (uR - uL) / aT) * 0.5; double alpha3 (dp_over_a2(drho, drhou, dE, uT, aT, gamma) rhoT * (uR - uL) / aT) * 0.5; // 耗散项按 |λ| α r 求和逐分量组装 for (int k 0; k 3; k) { double diss std::fabs(lam1) * alpha1 * rvec1(k, uT, aT) std::fabs(lam2) * alpha2 * rvec2(k, uT, aT) std::fabs(lam3) * alpha3 * rvec3(k, uT, aT); Fnum[k] 0.5 * (FL[k] FR[k]) - 0.5 * diss; } }这段代码里dp_over_a2是从 ΔU 反解压力跳变的辅助函数rhoT是 Roe 平均密度 √(ρL·ρR)rvec1/2/3是三个右特征向量分别是 [1, ũ−ã, H̃−ũã]ᵀ、[1, ũ, ½ũ²]ᵀ、[1, ũã, H̃ũã]ᵀ。逻辑上就是把状态差分解成三束波每束波按自身传播方向和速度大小决定耗散多少然后从中心通量里扣掉。参数 gamma 影响压力反解和特征向量改气体模型时这里是第一处要动的地方。注意上面代码里特征值直接取了绝对值没做熵修正。这是教学代码的常见写法跑 Sod 问题一般也能出图但在声速点附近可能出现小的非物理振荡后面避坑章节会给出修法。3.3 ENO 重构与三阶 Runge-Kutta 时间推进ENO 重构做的是对每个界面在左右两侧各用三个候选模板拟合二次多项式比较二阶差商的大小挑最光滑的那个取多项式在界面处的值。三阶 ENO 的候选模板是 {i−2, i−1, i}、{i−1, i, i1}、{i, i1, i2} 三组。// 三阶 ENO 重构返回界面 i1/2 左侧的重构值 double eno3Left(const double* q, int i, double dx) { // 三组模板对应的二次多项式在界面处的值 double q0 ( 2.0*q[i-2] - 7.0*q[i-1] 11.0*q[i]) / 6.0; double q1 (-1.0*q[i-1] 5.0*q[i] 2.0*q[i1]) / 6.0; double q2 ( 2.0*q[i] 5.0*q[i1] - 1.0*q[i2]) / 6.0; // 用二阶差商判断光滑程度 double d0 std::fabs(q[i-2] - 2.0*q[i-1] q[i]); double d1 std::fabs(q[i-1] - 2.0*q[i] q[i1]); double d2 std::fabs(q[i] - 2.0*q[i1] q[i2]); if (d0 d1 d0 d2) return q0; if (d1 d2) return q1; return q2; }逐分量对 ρ、ρu、E 做重构就能得到界面左侧的守恒量向量把索引反过来就是右侧。这里有个细节值得说重构作用在守恒量还是原始变量上社区里两种做法都有。守恒量重构实现简单、守恒性好但在接触间断附近容易出现密度和压力的小幅负值原始变量重构特征变量重构更讲究能缓解这个问题代价是要多做几次变量转换。这份代码用的是守恒量重构配合 CFL 取 0.4 左右基本稳。时间推进用三阶 SSP-RK三步的系数是 1、1/43/4、1/32/3// 三阶 SSP-RKL 为空间算子返回 dU/dt void rk3Step(std::vectordouble U, double dt, SpatialOp L) { std::vectordouble U1(U.size()), U2(U.size()), k1(U.size()), k2(U.size()); L(U, k1); for (size_t i 0; i U.size(); i) U1[i] U[i] dt * k1[i]; L(U1, k2); for (size_t i 0; i U.size(); i) U2[i] 0.75*U[i] 0.25*(U1[i] dt*k2[i]); L(U2, k1); for (size_t i 0; i U.size(); i) U[i] (1.0/3.0)*U[i] (2.0/3.0)*(U2[i] dt*k1[i]); }三个阶段的系数不是随便凑的SSP 性质保证只要前向欧拉在 CFL 约束下稳定这个三阶组合也稳定。参数 dt 由 CFL 条件给出dt CFL · Δx / max(|u| a)最大值在整场上取。每次 RK 子步之后都要重算一次因为中间状态的速度可能变化。4. 跑通与验证编译参数、运行选项和结果曲线怎么对4.1 编译与依赖这类代码的依赖通常只有标准库加一个绘图脚本。C 部分用 CMake 或者一个手写 Makefile 都能编我一般先用最笨的方式确认能过g -O2 -stdc17 -Iinclude src/*.cpp -o shocktube ./shocktube config/sod.cfg-O2别省激波管虽然网格不多但如果网格数调到几千、时间步上万不开优化跑起来会明显变慢。-stdc17看代码里有没有用std::filesystem或者结构化绑定没有的话 C11 也能过。编译报错最常见的是头文件路径。#include params.h这种写法要求-Iinclude如果代码里写的是#include ../include/params.h就得从src/目录里编译或者把路径改掉。另外有些版本用了omp.h做并行忘了加-fopenmp会报一堆未定义符号这个错误信息很长但根因就一个。画图部分一般是 Python 脚本matplotlibnumpy两件套就够import numpy as np import matplotlib.pyplot as plt # 读数值解和解析解两列x 和变量值 num np.loadtxt(output/density.dat) exa np.loadtxt(output/density_exact.dat) plt.plot(num[:, 0], num[:, 1], o, ms3, labelnumerical) plt.plot(exa[:, 0], exa[:, 1], -, lw1.5, labelexact) plt.xlabel(x); plt.ylabel(density) plt.legend(); plt.tight_layout() plt.savefig(output/density.png, dpi150)ms3让数值解的点不至于糊成一片dpi150保证截图发出去看得清。4.2 运行参数与输出格式配置文件一般是键值对或者简单的文本行重点盯这几个参数参数典型值作用改动的后果N2001000单元总数加倍后激波更锐但 dt 减半耗时约四倍CFL0.40.8时间步系数超过 1 大概率发散低于 0.2 精度提升有限但很慢T_end0.2结束时间太大波会撞到边界受边界条件影响flux_schemeroe/lf通量格式换格式直接影响间断分辨率recon_order1 / 3重构阶数一阶很耗散三阶激波更窄但更容易振荡gamma1.4比热比改气体时用空气别动输出一般是每个输出时刻一个.dat列是 x、ρ、u、p、E有的版本还会多输出一列马赫数。文件名里带时刻方便画动画。如果只输出最后一个时刻配置文件里会有一个output_interval之类的参数或者干脆就是循环结束后写一次。提示跑之前先把output/目录建好很多代码不会自动创建目录写文件失败时不一定报错结果就是一个空文件夹白等半天。4.3 用解析解验证密度、压力、速度三条曲线Sod 问题的解析解是自相似的按 ξ (x − 0.5)/t 分段给出。代码包里如果带了exact_solution.cpp或者一个 Python 脚本直接拿来对比如果没带用误差范数也能做粗略验证。判断跑得对不对看三个特征密度曲线上激波位置大约在 x ≈ 0.85接触间断在 x ≈ 0.68稀疏波头部在 x ≈ 0.27。数值解的激波会被抹开两三格这是正常的数值耗散网格加密后会变窄。如果激波位置明显偏移先查 dt 是不是算错了CFL 条件里用的是全场最大特征速度不是局部值。压力曲线上接触间断两侧压力应该相等密度不等。如果压力在接触间断处出现明显跳变说明格式的接触间断捕捉能力有问题Roe 格式正常不会这样出问题多半是特征向量或者波强度公式写错了。速度曲线上稀疏波区域应该是光滑过渡接触间断处速度连续。如果速度在接触间断处出现台阶同样指向特征分解的错误。我一般会算一个 L1 误差随网格数收敛的表一阶格式误差按一阶速率降三阶重构理论上能到三阶实际因为激波处降阶整体速率会落在 1 到 2 之间这是正常现象不用怀疑代码。5. 避坑与排查激波管代码最容易翻车的五个地方5.1 现象跑几十步后密度出现负值程序输出 nan原因Roe 格式在声速点附近没有做熵修正膨胀激波处特征值趋近于零耗散项消失数值解出现非物理跳跃接着压力被算成负数声速开根号就完蛋。解决对特征值加一个小的熵修正常用的 Harten 修正形式是 |λ| 替换成 (λ² δ²)/(2δ)当 |λ| δ 时δ 取 0.1 倍当地声速左右。代码里找到std::fabs(lam1)那几行包一个函数替换掉。别把 δ 取太大太大会把接触间断抹平0.050.1 倍声速是比较稳的区间。5.2 现象换成 Lax-Friedrichs 后能跑换回 Roe 立刻崩原因两种格式的稳定域不一样。LF 的数值耗散大能压住很多问题Roe 耗散小对 CFL 更敏感。如果 CFL 按 LF 调到 0.9换 Roe 就超了。解决换格式时把 CFL 降到 0.4 重新试稳定之后再往上加。另外检查一下时间步是不是每步重算的有的代码图省事在初始化时算一次 dt 就不动了格式一换最大特征速度变了dt 就不对了。5.3 现象接触间断处密度曲线出现一堆小振荡原因重构用的是守恒量而不是特征变量跨过接触间断时多项式拟合把密度和压力的跳变混在一起产生过冲。三阶 ENO 的模板选择在接触间断附近也可能选错。解决两个方向。一是把重构改到特征变量上先做局部特征分解对每个特征分量单独重构再组装回去效果好但代码量大二是加限制器在 ENO 外面套一层 WENO 权重或者 TVD 限制把过冲压下去。快速验证的话先把重构阶数降到一阶看振荡是否消失确认问题出在重构环节。5.4 现象结果整体还行但激波位置比解析解偏左或偏右原因边界条件或者初始界面位置搞错了。初始界面必须严格落在网格界面上如果网格数是奇数x 0.5 正好落在某个单元中心跳变就被分到了两个单元里激波位置会偏半格。解决网格数取偶数保证 0.5 落在界面上。另外检查左右状态的赋值条件用的是还是这个符号在界面处决定跳变归属虽然影响只有一格但对比解析解时会看得出来。5.5 现象程序跑完不报错但输出文件是空的或者只有一行原因输出逻辑写在循环里但没加刷新或者输出间隔参数设成了比总步数还大的值循环一次都没触发写文件。也有可能是文件路径用了相对路径运行目录不对文件写到别处去了。解决先确认输出目录存在再把输出间隔改成每 10 步写一次试跑。C 的ofstream在析构时才保证 flush如果程序中途异常退出缓冲区里的内容会丢重要输出后面手动加一句flush()。6. 二次开发换格式、加限制器、推到二维验证怎么做想把这份代码改成自己的东西我建议按「先验证、再替换、后扩展」的顺序来别一上来就大改。第一步是建立自己的验证基线。把 Sod 问题跑通把密度、压力、速度三条曲线和解析解存成图片记下网格数、CFL、误差范数。以后任何改动先跟这组基线比误差没变差才算通过。这一步看着笨但能省掉大量「改了之后不知道是不是改坏了」的纠结。第二步是替换通量格式。flux_roe.cpp和flux_lf.cpp的函数签名是一致的加 HLLC 只要新建一个flux_hllc.cpp实现同样的接口在solver.cpp的 switch 里加一个分支。HLLC 的关键是估计三个波速 SL、SM、SR中间那道接触波的速度 SM 用压力平衡条件解出来。改完拿基线一比接触间断应该比 Roe 更锐激波处略差一点这是正常的。第三步是加限制器。如果目标是高分辨率捕捉激波又不想振荡把 ENO 换成 WENO 是常见路线。三阶 WENO 的模板选择和光滑指示子计算量比 ENO 大但并行性好、边界处理简单。实现的时候注意光滑指示子里的 ε 不要取太小1e-6 是常用值太小会让权重在均匀区也剧烈变化。第四步是推到二维。一维的roeFlux是逐面调用的二维里把它改成按面法向旋转坐标先把左右状态投影到法向和切向用一维公式算法向通量再旋转回来。这一步最容易错的是切向速度的处理切向速度在 Roe 平均里要跟着一起平均但特征向量里对应的分量不一样抄公式的时候容易漏。最后说个我自己的习惯每次改完格式我都会用一阶重构 很小的时间步跑一遍确认空间算子本身没问题再往上加阶数。血泪经验是有一回直接在三阶框架里调通量振荡来源根本分不清是重构还是通量来回折腾了一整天才定位到是特征向量里少了一项。从那以后我每次动数值格式都强制走一遍「一阶验证 → 三阶对比 → 加密网格」的流程虽然慢一点但不会在错误的路上越走越远。希望这份拆解能帮你把这份代码跑起来也帮你少踩几个我踩过的坑。本文还有配套的精品资源点击获取