简介面向计算流体力学初学者与高校课程实训这份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 格式先用向前差分做一次预测再用向后差分做一次校正时间方向保留二阶精度。它属于显式格式存在 CFLCourant-Friedrichs-Lewy稳定性限制二维问题的稳定条件比一维更严格后面会专门说。2.2 预测步、校正步和最终更新式的推导以一维为起点更容易建立直觉。MacCormack 格式的预测步用向前空间差分[ \bar{u}i^{n1} u_i^n - a \frac{\Delta t}{\Delta x} \left( u{i1}^n - u_i^n \right) ]校正步改用向后空间差分[ u_i^{n1} \frac{1}{2} \left[ u_i^n \bar{u}_i^{n1} - a \frac{\Delta t}{\Delta x} \left( \bar{u}i^{n1} - \bar{u}{i-1}^{n1} \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_{i1,j} - u^n_{i,j} \right) - \frac{b\Delta t}{\Delta y} \left( u^n_{i,j1} - u^n_{i,j} \right) ][ u^{n1}{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 格式虽然老但作为理解显式格式、预测-校正思想、二阶精度特性的载体它依然是最顺手的工具之一。希望帮到你。本文还有配套的精品资源点击获取