企业官网建设流程全解析

简介：这份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.4～0.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, i+1}、{i, i+1, i+2} 三组。

// 三阶 ENO 重构，返回界面 i+1/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[i+1]) / 6.0; double q2 = ( 2.0*q[i] + 5.0*q[i+1] - 1.0*q[i+2]) / 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[i+1]); double d2 = std::fabs(q[i] - 2.0*q[i+1] + q[i+2]); if (d0 <= d1 && d0 <= d2) return q0; if (d1 <= d2) return q1; return q2; }

逐分量对 ρ、ρu、E 做重构，就能得到界面左侧的守恒量向量；把索引反过来就是右侧。这里有个细节值得说：重构作用在守恒量还是原始变量上，社区里两种做法都有。守恒量重构实现简单、守恒性好，但在接触间断附近容易出现密度和压力的小幅负值；原始变量重构（特征变量重构更讲究）能缓解这个问题，代价是要多做几次变量转换。这份代码用的是守恒量重构，配合 CFL 取 0.4 左右基本稳。

时间推进用三阶 SSP-RK，三步的系数是 1、1/4+3/4、1/3+2/3：

// 三阶 SSP-RK，L 为空间算子（返回 dU/dt） void rk3Step(std::vector<double>& U, double dt, SpatialOp L) { std::vector<double> 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 -std=c++17 -Iinclude src/*.cpp -o shocktube ./shocktube config/sod.cfg

-O2别省，激波管虽然网格不多，但如果网格数调到几千、时间步上万，不开优化跑起来会明显变慢。-std=c++17看代码里有没有用std::filesystem或者结构化绑定，没有的话 C++11 也能过。

编译报错最常见的是头文件路径。#include "params.h"这种写法要求-Iinclude，如果代码里写的是#include "../include/params.h"，就得从src/目录里编译或者把路径改掉。另外有些版本用了<omp.h>做并行，忘了加-fopenmp会报一堆未定义符号，这个错误信息很长但根因就一个。

画图部分一般是 Python 脚本，matplotlib+numpy两件套就够：

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", ms=3, label="numerical") plt.plot(exa[:, 0], exa[:, 1], "-", lw=1.5, label="exact") plt.xlabel("x"); plt.ylabel("density") plt.legend(); plt.tight_layout() plt.savefig("output/density.png", dpi=150)

ms=3让数值解的点不至于糊成一片，dpi=150保证截图发出去看得清。

4.2 运行参数与输出格式

配置文件一般是键值对或者简单的文本行，重点盯这几个参数：

输出一般是每个输出时刻一个.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.05～0.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 平均里要跟着一起平均，但特征向量里对应的分量不一样，抄公式的时候容易漏。

最后说个我自己的习惯：每次改完格式，我都会用一阶重构 + 很小的时间步跑一遍，确认空间算子本身没问题，再往上加阶数。血泪经验是有一回直接在三阶框架里调通量，振荡来源根本分不清是重构还是通量，来回折腾了一整天才定位到是特征向量里少了一项。从那以后我每次动数值格式都强制走一遍「一阶验证 → 三阶对比 → 加密网格」的流程，虽然慢一点，但不会在错误的路上越走越远。希望这份拆解能帮你把这份代码跑起来，也帮你少踩几个我踩过的坑。

本文还有配套的精品资源，点击获取