一维无粘 Burgers 方程的激波形成问题：MacCormack 格式求解

一维无粘 Burgers 方程是研究可压缩流激波捕捉格式时最经典的一维原型问题。本文完整介绍初值间断如何演化成一条向右传播的激波这一经典算例出处J.D. Anderson《Computational Fluid Dynamics: The Basics with Applications》第 7 章内容包括问题描述、CFD 理论背景、MATLAB 与 C 源码、运行结果数据与可视化图像。1. 问题描述1.1 控制方程在一维均匀网格上求解非线性双曲型方程无黏 Burgers 方程∂u/∂t ∂/∂x (u²/2) 0, 0 ≤ x ≤ 2其中f(u) u²/2为通量未知量u可理解为速度无黏 Burgers 方程与可压缩 Euler 方程的非线性结构同源是后者的标量原型。1.2 初始条件Riemann 问题初值在x 1处存在一个间断u(x, 0) 1.0, 0 ≤ x ≤ 1 u(x, 0) 0.5, 1 x ≤ 21.3 边界条件Dirichletu(0, t) 1.0, u(2, t) 0.51.4 理论解Rankine–Hugoniot 激波位置左侧高速流体u 1不断追赶右侧低速流体u 0.5初值间断将演化成一条向右传播的激波。由 Rankine–Hugoniot 关系激波速度为两侧状态的平均特征速度s (a(u_L) a(u_R)) / 2 (1 0.5) / 2 0.75激波位置随时间线性右移x_s(t) 1 s·t 1 0.75·t1.5 数值参数计算域 L 2 网格节点 nx 81Δx L/(nx-1) 0.025 时间步数 nt 120Δt 0.0125模拟至 t 1.5 s CFL 数 ν max|u|·Δt/Δx 1×0.0125/0.025 0.5≤ 1稳定2. CFD 理论背景2.1 无粘 Burgers 方程最简的非线性双曲守恒律一维守恒律的标准形式为∂u/∂t ∂f(u)/∂x 0当通量取f(u) u²/2时即为一维无粘 Burgers 方程。它是非线性的即使初值光滑特征线也会在有限时间内相交形成梯度趋于无穷的间断激波即弱解。因此它是理解可压缩 Euler 方程中激波捕捉问题的最小原型也是 Anderson 教材里从线性对流方程过渡到非线性方程组Euler 方程的桥梁算例。本题初值是典型的 Riemann 问题左右两个常状态由一个间断隔开。两侧特征速度a(u) u分别为 1 与 0.5左侧更快因此该间断必然演化成激波而非稀疏波若左侧慢于右侧则退化为稀疏波。激波传播速度取两侧状态的算术平均s 0.752.2 MacCormack 预估–校正格式MacCormack 格式是一种显式、二阶精度的时间/空间离散实现非常简单先做一次后差预估再做一次前差校正并取平均。记F u²/2。预估步后差通量F_i^n 0.5·(u_i^n)² ū_i u_i^n − (Δt/Δx)·(F_i^n − F_{i−1}^n) F̄_i 0.5·(ū_i)²校正步前差通量u_i^{n1} 0.5·[ u_i^n ū_i − (Δt/Δx)·(F̄_{i1} − F̄_i) ]在内点区间i 2, …, nx−1做预估i 1, …, nx−1做校正之外两端节点每步被 Dirichlet 边界值1.0 / 0.5重新赋值“钉扎”。整个时间步推进可紧凑地写成整向量形式ubar[1.0,u(2:80)-(dt/dx)*(F(2:80)-F(1:79)),0.5];% 预估内点 2..80边界钉扎Fbar0.5*(ubar.*ubar);unew0.5*(u(1:80)ubar(1:80)-(dt/dx)*(Fbar(2:81)-Fbar(1:80)));% 校正内点 1..80u[unew,0.5];2.3 CFL 稳定性条件对显式格式时间步必须满足CFL 条件ν max_i |u_i| · Δt / Δx ≤ 1本题取Δt 0.0125、Δx 0.025得ν 0.5格式稳定程序运行全程不发散。2.4 数值耗散与色散Gibbs 型振荡/过冲MacCormack 格式二阶精度、不含人工耗散在强间断附近会产生典型的小幅Gibbs 型振荡/过冲Anderson 称之为色散误差。尤其当激波在t ≈ 1.33 s后离开右边界而右边界仍被钉扎在u 0.5时末段节点会积累较明显的过冲见第 5 节数据。这正是后来引入人工黏性、通量限制器TVD以及 Riemann 求解器Godunov 族等技术的动机。3. MATLAB 源码以下给出该算例的完整 MATLAB 实现网格与参数 → MacCormack 时间推进并保存 t 0 / 0.25 / 0.5 / 1.0 / 1.5 五个快照→ 多时刻剖面绘图 → 激波位置数值检测并与 Rankine–Hugoniot 理论值对比。可直接在 MATLAB / Octave 中运行。%% % 一维无粘 Burgers 方程 —— 激波形成问题% MacCormack 预估-校正格式% 参考: J.D. Anderson, CFD: The Basics with Applications, Ch.7%% clear;clc;close all;%% 1. 网格与时间参数L2;% 计算域长度nx81;% 网格节点数nt120;% 时间步数dxL/(nx-1);% 空间步长dt0.0125;% 时间步长xlinspace(0,L,nx);% 网格坐标CFLmax(1)*dt/dx;% 初始 CFL 数fprintf(CFL %.3f (1 即稳定)\n,CFL);%% 2. 初始条件Riemann 问题u0.5*ones(1,nx);u(x1)1.0;% 保存若干时刻用于画图snap_t[00.250.51.01.5];% 想要观察的物理时刻snap_ucell(numel(snap_t),1);snap_u{1}u;% t 0isnap1;t_now0;%% 3. 时间推进MacCormackforn1:nt F0.5*u.^2;% 当前时刻通量% ---- 预估步空间后差 ----ubaru;ubar(2:end)u(2:end)-dt/dx*(F(2:end)-F(1:end-1));ubar(1)1.0;ubar(end)0.5;% 边界保持Fbar0.5*ubar.^2;% 预估通量% ---- 校正步空间前差 ----unewu;unew(1:end-1)0.5*(u(1:end-1)ubar(1:end-1)...-dt/dx*(Fbar(2:end)-Fbar(1:end-1)));unew(end)0.5;% 更新uunew;t_nowt_nowdt;% 保存特定时刻的解ifisnapnumel(snap_t)abs(t_now-snap_t(isnap1))dt/2isnapisnap1;snap_u{isnap}u;endend%% 4. 结果可视化figure(Color,w,Position,[100100700450]);hold on;colorslines(numel(snap_t));fork1:numel(snap_t)plot(x,snap_u{k},LineWidth,1.8,Color,colors(k,:),...DisplayName,sprintf(t %.2f,snap_t(k)));endgrid on;box on;xlabel(x);ylabel(u(x,t));title(一维 Burgers 方程激波形成);legend(Location,northeast);ylim([0.31.1]);%% 5. 激波位置与理论值对比Rankine–Hugoniotxs_numx(u(1.00.5)/2[diff(u)0,0]);% 数值激波位置xs_theory10.75*max(snap_t);% 理论位置fprintf(\n t%.2f 时数值激波位置 x_s%.3f, 理论值 x_s%.3f\n,...snap_t(end),xs_num(1),xs_theory);4. C 源码本节给出与上述 MATLAB 数值逻辑完全等价的 C 实现便于工程集成与批量运行。工程由 5 个文件组成文件作用burgers_vec.h入口函数声明burgers_vec.cpp数值求解主体时间推进 激波检测main.cpp程序入口main()rtwtypes.h跨平台基础类型定义CMakeLists.txt可选的 CMake 构建脚本以上 5 个文件构成一个完整工程入口函数为main()其输出见第 5 节。4.1burgers_vec.h/* * burgers_vec.h * 一维无粘 Burgers 方程激波形成问题MacCormack 预估-校正——头文件 */#ifndefBURGERS_VEC_H#defineBURGERS_VEC_H#includecmath#includecfloat#includecstddef#includecstdlib#includecstdio#includertwtypes.h#includevector#includestring/* 入口函数声明 */externvoidburgers_vec();#endif4.2burgers_vec.cpp/* * burgers_vec.cpp * 一维无粘 Burgers 方程激波形成问题MacCormack 预估-校正——数值实现 * 与第 3 节 MATLAB 实现的数值逻辑一致 * 预估步后差、校正步前差、端点每步钉扎边界值 1.0 / 0.5 */#includeburgers_vec.h#includecmathvoidburgers_vec(){doubleL;/* 计算域长度 */intnx;/* 网格节点数 */intnt;/* 时间步数 */doubledx;/* 空间步长 */doubledt;/* 时间步长 */doublex[81];/* 网格坐标 */doubleCFL;/* CFL 数 */doubleu[81];/* 解向量 */intn;/* 时间步循环 */doubleF[81];/* 当前通量 */doubleubar[81];/* 预估解 */doubleFbar[81];/* 预估通量 */doubleunew[80];/* 校正后的内点 */doublexs_num;/* 数值激波位置 */intfound;/* 是否已找到 */inti;doublexs_theory;/* 理论激波位置 */L2.0;nx81;nt120;dxL/(nx-1);/* 0.025 */dt0.0125;doubletv1(L-0.0)/(double)(nx-1);/* linspace(0,L,nx) 步长 */for(intk0;k81;k){x[k]0.0(double)k*tv1;}CFLdt/dx;printf(CFL %.3f (1 stable)\n,CFL);/* 初始条件u 0.5 全域x 1 处再加 0.5 - 1.0 */for(intk0;k81;k){u[k]0.5;}for(intk0;k81;k){u[k]u[k]0.5*(x[k]1.0);}/* MacCormack 时间推进 */for(n1;n120;n){/* 当前通量 F u.^2 / 2 */for(intk0;k81;k){F[k]0.5*(u[k]*u[k]);}/* 预估步后差内点 k1..79两端钉扎 */ubar[0]1.0;for(intk1;k80;k){ubar[k]u[k]-(dt/dx)*(F[k]-F[k-1]);}ubar[80]0.5;/* 预估通量 */for(intk0;k81;k){Fbar[k]0.5*(ubar[k]*ubar[k]);}/* 校正步前差内点 k0..79最右端钉扎 */for(intk0;k80;k){unew[k]0.5*((u[k]ubar[k])-(dt/dx)*(Fbar[k1]-Fbar[k]));}for(intk0;k80;k){u[k]unew[k];}u[80]0.5;}/* 激波位置数值检测第一个满足 u(i)0.75 且 u(i1)u(i) 的网格点 */xs_num0.0;found0;for(i1;i81;i){/* 内点 1..80 */if(u[i]0.75){if(u[i1]u[i]){if(found0){xs_numx[i];found1;}}}}xs_theory1.00.75*1.5;/* x_s(1.5) 1 0.75×1.5 */printf(t %.2f: numerical x_s %.3f, theory x_s %.3f\n,1.5,xs_num,xs_theory);}说明C 中数组下标为 0-based对应 MATLAB 的 1-based 节点偏移 1预估/校正均只在 79~80 个内点上计算两端节点每步由边界值重写因此与 MATLAB 切片写法逐点等价。x 1的逻辑判断在 C 中返回 0/1与 MATLAB 中(x1)用作数值一样。4.3main.cpp/* main.cpp —— 程序入口 */#includeburgers_vec.hintmain(){burgers_vec();return0;}4.4rtwtypes.h/* rtwtypes.h —— 跨平台基础类型定义 */#ifndefRTWTYPES_H#defineRTWTYPES_H#includecstddef#includecstdinttypedefdoublereal_T;typedeffloatreal32_T;typedefintint32_T;typedefint16_tint16_T;typedefint8_tint8_T;typedefuint32_tuint32_T;typedefuint16_tuint16_T;typedefuint8_tuint8_T;typedefuint8_tboolean_T;typedefcharchar_T;#endif/* RTWTYPES_H */4.5CMakeLists.txtcmake_minimum_required(VERSION 3.10) project(burgers_solver CXX) set(CMAKE_CXX_STANDARD 14) add_executable(burgers_app main.cpp burgers_vec.cpp ) target_include_directories(burgers_app PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})5. 运行结果与数据5.1 程序输出运行第 4 节的 C 实现数值格式与第 3 节 MATLAB 实现一致两者关键数值输出相同得到CFL 0.500 (1 stable) t 1.50: numerical x_s 0.450, theory x_s 2.125即CFL dt/dx 0.0125 / 0.025 0.500 显式格式稳定ν ≤ 1 理论激波 x_s(t1.5) 1 0.75×1.5 2.125 Rankine–Hugoniot 关系 数值检测 x_s 0.450 判据第一个命中的网格点见 5.3 说明5.2 各快照时刻剖面数据各时刻剖面的最小/最大值max u已包含 MacCormack 格式的过冲时刻 t (s)min umax u含过冲形态0.000.51.000000阶跃初值x 1 处间断1.0/0.50.250.51.157157间断处出现数值抹平与小幅振荡0.500.51.130594波形左陡右缓、整体向右推进1.000.51.131811波形逼近右边界1.500.51.331206激波已离开右边界末段节点出现明显过冲5.3 结果说明1CFL 与稳定性。ν 0.5 ≤ 1满足显式格式的稳定条件程序运行全程不发散解的单调推进行为与理论一致。2激波位置检测为何输出 0.450而不是理论值 2.125理论激波位置x_s(t) 1 0.75·t在t 1.5时为2.125已大于右边界x 2。这意味着激波其实在t_exit ≈ (2 − 1) / 0.75 ≈ 1.333 s就已经离开计算域t 1.5 时理论上激波位于域外。而程序中的检测判据只能返回网格坐标最大 2.0同时t 1.5 时u 1的高位平台上布满约1e-12量级的浮点噪声扫描会先在这些噪声处命中u(i1) u(i)于是输出了第一个命中点x 0.450第 19 号节点。也就是说激波位置判据本身受网格坐标上限 平台噪声影响输出的 0.450 不能直接当作物理激波位置物理上激波以s 0.75匀速右移并在约 1.33 s 时离开计算域。若需在激波出域后继续跟踪应把判据限定在波形前缘靠近右端附近或改为记录每个时间步的传播前沿。3MacCormack 格式的过冲。t 1.5 时末段节点最大过冲u_max − 1 ≈ 0.331。这是二阶无耗散格式在激波推出边界 Dirichlet 钉扎条件下的典型色散振荡与 2.4 节的理论背景吻合在更长时间积分或更粗网格下会变得更明显工程中常配合人工黏性 / 通量限制器使用。6. 可视化结果以下图像由运行数据绘制左图为 t 0 / 0.25 / 0.5 / 1.0 / 1.5 五个时刻的整场剖面右图为 t 1.5 时刻最终剖面与激波检测结果。图 1不同时刻的速度剖面激波的形成与传播图 2t 1.5 s 最终剖面与激波检测图中可见完整物理过程初始阶跃 → 间断处轻微的 Gibbs 型振荡 → 波形整体向右推进并逐渐逼近右边界t 1.5 时高位平台几乎贯穿全域理论激波位置2.125已超出右边界图中以竖线标出离散判据先命中了平台噪声点x 0.450。7. 小结该算例完整展示了非线性双曲守恒律从初值间断 → 激波形成 → 传播/离开计算域的物理过程是理解可压流激波捕捉格式TVD、WENO、Godunov 等最直接的入门原型。MacCormack 格式实现简洁、二阶精度但本身不含耗散在强间断处会产生振荡/过冲——这正是后续引入人工黏性、通量限制器或 Riemann 求解器的动机。通过 MATLAB 与 C 两套等价实现可以清晰对照向量化切片写法与逐点循环写法的对应关系两套代码的关键数值输出完全一致结果可复现、可审计。8. 参考资料J.D. Anderson,Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill, 1995第 7 章激波捕捉与 Burgers 方程算例。R.J. LeVeque,Numerical Methods for Conservation Laws, Birkhäuser, 1992守恒律数值方法与 Riemann 问题。