简介：面向计算流体力学初学者与高校课程实训，这份PDF资料围绕二维对流方程的数值解法，完整呈现一阶迎风C格式、A格式、Lax-Wendroff格式及时间分裂法下的Matlab实现。内容按实习报告体例编排：先给出初值条件（方形波在-16到16区域传播）与边界条件（左、下边界入流，右、上边界出流），再推导各差分格式和误差传播方程，重点分析时间步长与空间步长比值取0.5、1、2时的稳定性，并对迎风C格式在高γ下发散的原因给出理论解释。文档绘制了四种格式在不同时刻的u-x-y三维图，部分结果制作为动画，直观展示波的传播、格式耗散和相位误差；附录附有完整实验程序代码，便于读者复现、验证并开展参数拓展。资源共1个PDF文件，大小约25.94MB，排版规范，层次清晰，既适合航空航天、力学专业上机实习参考，也是CFD入门与Matlab编程自学的实用材料。目前已有318人学习下载，该实习报告结构完整、理论分析与程序结合紧密，参考价值较强。
1. 为什么二维对流方程值得用 MacCormack 格式手写一遍
做 CFD 的人迟早会撞上这道坎：一维对流方程看着简单，换成二维之后，数值格式的稳定性、边界处理和网格效应全都变了味道。标题里这个「AC 格式」在多数教材和开源代码里指的就是 MacCormack 格式，一种基于预测-校正两步的时间显式格式。它不靠隐式迭代，也不用解大型稀疏矩阵，用 Matlab 写出来往往只有几十行核心循环，却能清楚展示二维对流问题里最要命的数值耗散与色散行为。对正在学有限差分、准备课程作业或者想验证自己求解器正确性的人来说，这是一个性价比极高的实验对象。
二维对流方程和解一维最大的不同在于：流动方向不再和网格线对齐，斜向传播的波会被数值格式各向异性地扭曲。用 MacCormack 格式跑一个高斯峰的对流算例，你会直观看到峰被拉伸、产生轻微拖尾振荡，而这些现象在教科书的一维算例里很难暴露。这篇文章就沿着「方程离散 → Matlab 实现 → 边界与网格处理 → 参数调试 → 验证技巧」这条路走一遍，中间会穿插我实际调试时踩过的坑。适合已经会写一维对流代码、想往二维推进的读者，也适合被作业逼着必须从零搭一个能跑的求解器的人。
2. 先把控制方程和 AC 格式的离散逻辑理顺
2.1 二维对流方程的守恒与非守恒形式
二维对流方程的标准形式是：
[ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} + b \frac{\partial u}{\partial y} = 0 ]
其中 (a) 和 (b) 分别是 (x)、(y) 方向的对流速度，(u) 是被输运的标量场，比如温度、浓度或者某个被动标量。这个方程在物理上描述的是一个形状保持不变的信号以速度 ((a, b)) 沿直线运动。做数值验证时，大家最喜欢用高斯波包或者正弦波做初值，因为理论解就是初值的刚体平移，可以直接拿误差做定量评估。
如果流动速度本身不是常数，而是依赖 (u) 或者空间位置，那就必须用守恒形式：
[ \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} + \frac{\partial g(u)}{\partial y} = 0 ]
(f(u)=au)、(g(u)=bu) 时的守恒形式和非守恒形式在连续情形下等价，但离散后会有细微差别。我们这里只处理常系数线性问题，所以直接用非守恒形式就够了。不过写代码的时候我建议还是把通量函数拆开，方便以后扩展成非线性问题。
MacCormack 格式的思路来源于 Lax-Wendroff 格式：先用向前差分做一次预测，再用向后差分做一次校正，时间方向保留二阶精度。它属于显式格式，存在 CFL（Courant-Friedrichs-Lewy）稳定性限制，二维问题的稳定条件比一维更严格，后面会专门说。
2.2 预测步、校正步和最终更新式的推导
以一维为起点更容易建立直觉。MacCormack 格式的预测步用向前空间差分：
[ \bar{u}i^{n+1} = u_i^n - a \frac{\Delta t}{\Delta x} \left( u{i+1}^n - u_i^n \right) ]
校正步改用向后空间差分：
[ u_i^{n+1} = \frac{1}{2} \left[ u_i^n + \bar{u}_i^{n+1} - a \frac{\Delta t}{\Delta x} \left( \bar{u}i^{n+1} - \bar{u}{i-1}^{n+1} \right) \right] ]
把两步合起来的效果等价于：时间二阶精度、空间二阶精度，并带有一个与 Lax-Wendroff 形式相近的数值耗散项。换到二维，关键是预测步和校正步在 (x)、(y) 两个方向要错开差分方向，否则格式的稳定性会退化。最常见的做法是预测步对 (x) 用向前差分、对 (y) 用向前差分，校正步对 (x) 用向后差分、对 (y) 用向后差分，也就是全向前预测、全向后校正。
也可以让两个方向相互错开，比如预测步在 (x) 方向向前、(y) 方向向后，校正步反过来。这种做法能改善斜向传播时的相位误差，但写代码时必须一致。我自己一般用全向前-全向后，因为边界逻辑最简单，数值表现对大多数人做验证已经足够。
二维情形下，如果记预测步得到的中间量为 (u^*)，那么完整更新是：
[ u^*{i,j} = u^n{i,j} - \frac{a\Delta t}{\Delta x} \left( u^n_{i+1,j} - u^n_{i,j} \right) - \frac{b\Delta t}{\Delta y} \left( u^n_{i,j+1} - u^n_{i,j} \right) ]
[ u^{n+1}{i,j} = \frac{1}{2} \left[ u^n{i,j} + u^_{i,j} - \frac{a\Delta t}{\Delta x} \left( u^{i,j} - u^*{i-1,j} \right) - \frac{b\Delta t}{\Delta y} \left( u^_{i,j} - u^_{i,j-1} \right) \right] ]
注意 (u^) 需要存储在完整的网格空间里，因为校正步会用到 (u^{i-1,j}) 和 (u^*{i,j-1})，这意味着不能只开一个二维数组就原地更新完。
2.3 CFL 条件在二维显式格式里的变化
一维显式格式的 CFL 条件是 (\frac{a\Delta t}{\Delta x} \le 1)，也就是 CFL 数小于等于 1。二维情况下，MacCormack 格式的稳定条件从理论上分析会变成：
[ \frac{|a|\Delta t}{\Delta x} + \frac{|b|\Delta t}{\Delta y} \le 1 ]
注意这里是「加和」而不是「各自小于 1」。如果两个方向的 CFL 数都取 0.5，加起来已经是 1.0，刚好在边界上。实际调试时我不会贴着理论值跑，通常让两个方向各自的 CFL 数都在 0.4 左右，这样总和约 0.8，留出安全余量。原因是边界处理、非光滑初值都会引入额外的高频分量，这些分量在理论上并不被标准稳定性分析覆盖，真跑到理论边界时很容易在边界附近先爆掉。
如果你用的是不均匀网格，(\Delta x) 和 (\Delta y) 各点不同，CFL 条件要逐点检查。最省事的做法是找整个网格上最小的 (\Delta x) 和 (\Delta y) 来计算全局时间步长，这样保证所有点都安全，但代价是计算量变大。对教学算例来说，均匀网格就够用了。
2.4 为什么选择 MacCormack 而不是 Lax-Wendroff 或迎风格式
这里值得花一段说清楚选型理由。Lax-Wendroff 格式在二维情况下需要同时计算四个交叉导数项，推导和编码都比 MacCormack 繁琐；迎风格式一阶精度耗散太大，高斯峰跑不了几十步就变成一座平头山，用来做格式验证看不出真实精度。MacCormack 格式正好在编码复杂度、精度和稳定性之间取了平衡：没有交叉导数项，只需要两次简单差分，却依然有时间和空间二阶精度。
还有一个容易被忽略的点：MacCormack 格式适合后面扩展成非线性问题。你把常数 (a)、(b) 换成依赖 (u) 的通量函数，预测-校正的结构完全不用动，只需要把差分里的 (a u) 换成 (f(u)) 即可。这也是很多人在学习从线性对流到 Burgers 方程时，会拿 MacCormack 当跳板的原因。
3. 用 Matlab 从零搭一个二维 AC 格式求解器
3.1 参数设置与网格初始化的最小框架
不管用什么语言写 CFD，第一步永远是设置物理参数、网格参数和时间步进参数。下面这段代码是一个可以直接跑通的最小框架：
% 二维对流方程 MacCormack 格式求解 % 参数设置 Lx = 1.0; % x 方向计算域长度 Ly = 1.0; % y 方向计算域长度 nx = 101; % x 方向网格数 ny = 101; % y 方向网格数 dx = Lx / (nx - 1); % x 方向网格间距 dy = Ly / (ny - 1); % y 方向网格间距 a = 0.5; % x 方向对流速度 b = 0.3; % y 方向对流速度 cflx = 0.4; % x 方向 CFL 数 cfl_y = 0.4; % y 方向 CFL 数 dt = min(cflx * dx / abs(a), cfl_y * dy / abs(b)); t_end = 0.5; % 总计算时间 nt = ceil(t_end / dt); % 时间步数，向上取整 dt = t_end / nt; % 重新微调 dt，保证恰好落在 t_end x = linspace(0, Lx, nx); y = linspace(0, Ly, ny); [X, Y] = meshgrid(x, y); % 初值：高斯波包，中心位于 (0.3, 0.4) sigma = 0.05; u0 = exp(-((X - 0.3).^2 + (Y - 0.4).^2) / (2 * sigma^2)); u = u0;这段代码有几个关键点需要说明。首先，初值的空间尺度 (\sigma) 要和网格间距匹配，一般至少要有 5 个网格点覆盖波包的半宽。如果 (\sigma = 0.05) 而网格间距是 0.01，那么波包半宽大约覆盖 5 个点，跑起来勉强能接受；再小，数值色散会把波形撕碎。其次，dt先用两个方向的 CFL 约束取最小值，再重新调整让总时间恰好是整数步。这样后处理时对比理论解的位置会非常方便，不需要做时间插值。
3.2 核心时间步进循环：预测-校正两步的完整代码
下面这段是核心，我建议直接敲一遍而不是复制，敲的过程中你会注意到数组索引方向的问题：
% 分配预测步的临时数组 ustar = zeros(ny, nx); % 主时间循环 for n = 1:nt % 预测步：x 和 y 方向都用向前差分 ustar(2:end-1, 2:end-1) = u(2:end-1, 2:end-1) ... - a * dt / dx * (u(3:end, 2:end-1) - u(2:end-1, 2:end-1)) ... - b * dt / dy * (u(2:end-1, 3:end) - u(2:end-1, 2:end-1)); % 校正步：x 和 y 方向都用向后差分 u(2:end-1, 2:end-1) = 0.5 * (u(2:end-1, 2:end-1) + ustar(2:end-1, 2:end-1) ... - a * dt / dx * (ustar(2:end-1, 2:end-1) - ustar(1:end-2, 2:end-1)) ... - b * dt / dy * (ustar(2:end-1, 2:end-1) - ustar(2:end-1, 1:end-2))); % 边界更新 u(1, :) = u(2, :); % y 方向下边界 u(end, :) = u(end-1, :); % y 方向上边界 u(:, 1) = u(:, 2); % x 方向左边界 u(:, end) = u(:, end-1); % x 方向右边界 end这段代码的索引逻辑值得仔细讲。Matlab 的数组索引是(行, 列)，对应空间是(y, x)。u(2:end-1, 2:end-1)表示去掉四边后的内部网格点；u(3:end, 2:end-1)是同一批点的右邻居，对应 (x) 方向向前差分；u(2:end-1, 3:end)是上邻居，对应 (y) 方向向前差分。校正步里u(1:end-2, 2:end-1)是左邻居，u(2:end-1, 1:end-2)是下邻居。
边界条件这块，我上面用的是零梯度外推，也就是 Neumann 边界。这样做的好处是波包从边界出去时不会产生明显的反射，适合短时间验证。如果做长时间模拟，零梯度边界会轻微改变波包的形状，但你用一维算例对照会发现问题不大。另一种常见选择是周期边界，适用于信号永远不出域的测试。
3.3 结果可视化与理论解对比
数值算完不对比理论解等于白算。由于是常系数线性方程，二维对流方程的理论解就是初值整体平移：
% 理论解：把初值中心平移到新位置 x_shift = mod(0.3 + a * t_end, Lx); y_shift = mod(0.4 + b * t_end, Ly); u_exact = exp(-((X - x_shift).^2 + (Y - y_shift).^2) / (2 * sigma^2)); % 误差分析 err = u - u_exact; l2_err = sqrt(sum(err(:).^2) / sum(u_exact(:).^2)); max_err = max(abs(err(:))); fprintf('L2 相对误差: %e\n', l2_err); fprintf('最大绝对误差: %e\n', max_err); % 绘图对比 subplot(1, 3, 1); contourf(X, Y, u, 20); title('数值解'); axis equal; colorbar; subplot(1, 3, 2); contourf(X, Y, u_exact, 20); title('理论解'); axis equal; colorbar; subplot(1, 3, 3); contourf(X, Y, err, 20); title('误差分布'); axis equal; colorbar;注意mod函数处理的是平移出界后再从另一边进来的情形，配合零梯度边界并不完全自洽。严格来说，零梯度边界下理论解应该只保留还在域内的部分，越界的波包部分不算误差。一个更严谨的做法是把理论解的初始中心直接平移，不取模，然后把域外部分的解析值截断到 0，如下：
x_shift = 0.3 + a * t_end; y_shift = 0.4 + b * t_end; u_exact = exp(-((X - x_shift).^2 + (Y - y_shift).^2) / (2 * sigma^2)); u_exact(X < 0 | X > Lx | Y < 0 | Y > Ly) = 0;这样对比时，误差集中在波包靠近边界那侧，你能清楚看到边界处理引入的误差量级。对二阶格式，你的 L2 相对误差一般应该在 (10^{-3}) 到 (10^{-2}) 之间；如果误差到了 (10^{-1})，说明初值分辨率不够或者时间步长过大。
4. 稳定性边界、精度验证与网格收敛性测试
4.1 CFL 数变化对解的宏观影响
CFL 数不是越小越好，也不是靠近 1 越好。把cflx和cfl_y从 0.2 调到 0.8，你会看到三种典型现象：CFL 很小时误差主要由空间离散产生，而时间离散误差很小，所以波形保持得不错但计算步数多；CFL 在 0.5 附近时总误差最小，因为时间离散误差和空间离散误差部分抵消了；CFL 接近 1 时，如果初值不够光滑，波包尾部会出现小锯齿，再往前一点就直接发散。
我调试时常做的实验是固定网格、固定初值、只改 CFL 数，然后把 L2 误差画成曲线。你会发现误差不是单调变化的，而是一条U形曲线。这在 CFD 里很常见，不光是 MacCormack 格式，很多显式格式都有误差相消的现象。所以别盲目把 CFL 往小了压，否则不仅计算变慢，精度还可能反而变差。
4.2 网格收敛性测试：如何证明你的代码没有低级错误
判断一个格式写没写对，最可靠的方法是网格收敛性测试。如果你把网格间距减半，时间和空间都是二阶格式，那么误差应该缩小为原来的约四分之一。具体做法是固定 CFL 数不变，让网格从 (41 \times 41) 加密到 (81 \times 81)、(161 \times 161)，同时把时间步长按比例缩小，保持 CFL 数不变：
nx_list = [41, 81, 161]; ny_list = [41, 81, 161]; cfl_target = 0.4; for k = 1:length(nx_list) nx = nx_list(k); ny = ny_list(k); dx = Lx / (nx - 1); dy = Ly / (ny - 1); dt = min(cfl_target * dx / abs(a), cfl_target * dy / abs(b)); nt = ceil(t_end / dt); dt = t_end / nt; % 重新初始化高斯波包，运行求解器 % 计算 L2 误差，记录到数组 end计算收敛阶的公式为：
[ p = \frac{\log(E_{coarse} / E_{fine})}{\log(2)} ]
理想情况下 (p) 应该接近 2。如果算出来是 1 或者 0.5，常见原因有三个：一是边界条件实现成了一阶精度，二是初值的高斯波包在粗网格上没有被充分分辨，三是dt没有随网格加密而缩小，导致时间误差成为主导。
这里要特别提醒：网格收敛性测试要求初值在每一层网格上都能被解析。高斯波包的 (\sigma) 必须和网格无关地用同一个值，这样加密网格实际上是在提高波包的分辨率。如果 (\sigma) 随网格变化，收敛性曲线就会一团糟。
4.3 非均匀网格下的 CFL 陷阱
遇到需要局部加密的工况时，CFL 条件必须逐点判断。假设某个区域网格间距是全局最小值的五分之一，那么局部 CFL 数会放大五倍。一个容易踩的坑是只按全局平均间距选dt，结果在加密区直接爆炸。排查这类问题时，有一个好习惯是在代码里主动计算全局最大 CFL 数并打印出来：
cfl_actual = max(abs(a) * dt / dx, abs(b) * dt / dy); fprintf('实际最大 CFL 数: %f\n', cfl_actual);如果打印出来的值超过 0.9，即使没炸也说明你在拿稳定性余量换效率，风险很高。非均匀网格下更稳妥的做法是在计算前先扫描整个网格，找出最小间距，然后基于它选dt。
5. 边界处理与典型翻车现场：四个高频踩坑记录
5.1 数组索引方向写反导致的结果镜像
现象：计算出的波包移动方向完全反了，明明 (a > 0) 应该向右移动，结果波包往左跑。原因：在预测步或校正步里，向前差分和向后差分的索引写反了。Matlab 里u(3:end) - u(2:end-1)是向右看，对应正向流动；如果写成u(2:end-1) - u(1:end-2)，就变成向左看。解决：在写好循环后，先用一个只有 3 个格点的小规模数组手动算一遍，把预测步和校正步的每个差分拆开打印，对照公式逐项检查。
5.2 边界点没有参与预测-校正，导致边界附近波形扭曲
现象：波包远离边界时波形正常，一旦靠近边界就开始出现非对称变形，甚至局部凹陷。原因：边界点只做了简单赋值，没有参与格式本身的时间更新。零梯度外推假设边界点等于内点，但这个假设只在波包垂直入射边界时才近似成立。解决：如果要做更长时间的模拟，改用迎风单边差分更新边界点；或者直接扩大计算域，把边界推远，让波包在关心的区域内不与边界发生交互。我通常会把计算域扩大 30% 到 50%，只统计内部区域的结果。
5.3 用surf查看结果时波形「撕裂」，误以为格式不稳定
现象：等值线图平滑美观，但用surf(u)看三维曲面时，波包表面出现密密麻麻的细小锯齿，看起来像是数值振荡。原因：不是格式的问题，而是surf的渲染方式夸大了网格尺度上的微小色散。解决：先用contourf或pcolor确认等值线是否光滑，如果等值线正常，那这个锯齿就是视觉假象；也可以用smooth类的后处理把可视化用的数据做滑平均，但注意不要拿处理后的数据去算误差。
5.4 计算时间步数过多导致误差累积异常
现象：把t_end设得很大，比如从 0.5 改到 5.0，结果误差不是线性增长，而是加速恶化。原因：零梯度边界条件允许波包部分能量离开计算域，但理论解对比时如果你还按「平移不衰减」算，误差自然越来越大。这不是格式不稳定，而是边界条件与理论解模型不匹配。解决：要么限制t_end使波包始终在域内，换周期边界配合周期性理论解，要么把误差计算限定在波包仍保持高幅值的区域。最省心的是让波包跑完整个域后仍然留在内部，用一组固定参考时间点做对比。
6. 进阶验证技巧：用一维参考解和模态分解定位误差来源
当你确认二维代码基本能跑，但误差总觉得偏大时，有一个很实用的排查手段：把所有网格设为一行，比如ny = 3，并把b设为 0，此时二维求解器应该退化成一维求解器。拿它和一维单独写的 MacCormack 结果的误差对比。如果两者对不上，说明二维数组索引或者y方向的残差处理有潜在问题，这和初值无关，纯粹是编码层面的 bug。
第二个技巧是观察误差的空间分布模式。MacCormack 格式是二阶精度，它的主导误差项是色散误差，表现为高频振荡尾巴在波包后方拖出一条长须。如果你的误差分布不是这种模式，而是整体抬高、形状对称变形，那更可能是耗散过大，常见原因是网格太粗或波包太小。误差模式的好坏能帮你快速判断该加密网格还是换更高分辨率的初值。
我个人的一个习惯是，每次改完代码都先跑一组一维退化测试和网格收敛性测试，数据全部存档，再开始下一个参数的调试。这样就能在后续调整边界条件或转非均匀网格时，快速判断是不是引入了新的 bug。格式本身是成熟工具，真正出问题的几乎都是索引、边界和时间步选择。如果你照着上面的流程走一遍，应该能在一小时内跑通高斯波包算例并拿到接近二阶的收敛阶数据。MacCormack 格式虽然老，但作为理解显式格式、预测-校正思想、二阶精度特性的载体，它依然是最顺手的工具之一。希望帮到你。
本文还有配套的精品资源，点击获取
