简介：计算流体力学（CFD）二维对流方程在迎风C格式与A格式下的数值解，配套完整实习报告与Matlab编程实现。资源面向航空航天、土木建筑等领域需要掌握CFD差分格式与数值稳定性的学习者，尤其适合中山大学相关课程学生参考。内容覆盖迎风C格式、A格式、Lax-Wendroff格式以及时间分裂法四种离散方案，并针对Δt/Δx取0.5、1、2三种情形进行稳定性分析，还包含不同时间步下的u-x-y三维波形图与动画制作思路。压缩包共1个文件，为PDF格式，整体大小25.94MB，内含问题描述、初值条件、上下边界处理、误差传播理论分析及完整Matlab程序附录，可直接对照复现运行结果。已有319人浏览学习，适合正在攻读计算流体力学课程并需要Matlab编程验证差分格式的本科生及研究生使用。

1. 一道CFD上机题的完整解法：二维对流方程四种格式的Matlab实现

某学期，某高校航空航天学院的计算流体力学课布置了一次上机实习：在32×32的网格上求解二维对流方程，初值是一个高度为1的方形波，要求分别用迎风C格式、A格式以及自选的另外两种差分格式推进，改变γ=Δt/Δx为0.5、1、2，对比稳定性和波形演化，最后输出三维曲面图和动画。这份资源就是该上机题从差分公式到Matlab代码、从稳定性理论到边界处理的完整配套，包含四种格式的求解脚本、波形可视化以及动画制作模板，还给出了每种格式在三个γ取值下是否收敛的实测结论。对正在做CFD课设、或者想复现“二维对流方程在不同差分格式下稳定性条件”的开发者来说，可以直接套用并改参数。

我拿到资源之后，先把四个格式的脚本各膨胀了一遍，再从报错里把出流边界、循环范围、索引方向这些细节逐个摸清楚。下面按“先弄清楚格式，再逐段读代码，最后避开我踩过的坑”的顺序写。

2. 准备与格式框架：从初值方波到差分递推

2.1 计算域、网格与初值条件

资源把计算域设为 ( -16 \le x \le 16,\ -16 \le y \le 16 )，空间步长 (\Delta x=\Delta y=1)，所以得到 (33\times33) 的网格点。初值在 (-4\le x\le4,\ -4\le y\le4) 的方形区域内取 (u=1)，其余位置取 0，相当于一个孤立方波随对流方向平移。

dx = 1; dy = 1; x = -16:dx:16; % x方向网格 [X, Y] = meshgrid(x); % 生成二维网格坐标 U1 = zeros(size(X)); % 当前时间层的解 for i = 13:21 % 行索引范围 for j = 13:21 % 列索引范围 U1(i, j) = 1; % 方形波初值 end end

这里要注意，MATLAB的meshgrid生成矩阵时，行下标对应y方向、列下标对应x方向。循环里 i 从13到21、j 从13到21，正好对应物理坐标从 -4 到 4 的方形区间。把行列搞反的话，绘制surf时波形会旋转90度，后面的边界处理也会跟着错位。

2.2 迎风C格式：一阶差分如何体现“迎风”思想

迎风C格式的递推式是：

[ u_{i,j}^{n+1} = u_{i,j}^n - \frac{\Delta t}{\Delta x}(u_{i,j}^n - u_{i-1,j}^n) - \frac{\Delta t}{\Delta y}(u_{i,j}^n - u_{i,j-1}^n) ]

写成代码就是：

U2(i,j) = U1(i,j) ... - r*(U1(i,j) - U1(i-1,j)) ... % x方向用左偏迎风 - r*(U1(i,j) - U1(i,j-1)); % y方向用下偏迎风

由于对流项里速度为正，波沿x正方向和y正方向传播，离散时要从“信息来的方向”取差商，所以x方向取 (u_{i,j}-u_{i-1,j})，y方向取 (u_{i,j}-u_{i,j-1})。这种一阶迎风格式天然具备数值耗散，方波在推进过程中会被逐渐抹平，峰值减小、坡度变缓，这是格式特性，不是程序写错。

2.3 A格式、Lax-Wendroff格式与时间分裂法

A格式是二阶中心差分：

[ u_{i,j}^{n+1} = u_{i,j}^n - \frac{\Delta t}{2\Delta x}(u_{i+1,j}^n-u_{i-1,j}^n) - \frac{\Delta t}{2\Delta y}(u_{i,j+1}^n-u_{i,j-1}^n) ]

它不区分方向，直接取左右对称差分，精度比迎风格式高，但在二维条件下无条件不稳定，放大因子的模恒大于等于1，所以源程序里三种γ取值全部发散。运行结果里A格式的波很快就会变成高频振荡，继而溢出成NaN。

Lax-Wendroff格式在中心差分基础上补了二阶耗散项：

[ u_{i,j}^{n+1} = u_{i,j}^n - \frac{r}{2}(u_{i+1,j}^n-u_{i-1,j}^n) - \frac{r}{2}(u_{i,j+1}^n-u_{i,j-1}^n) + \frac{r^2}{2}(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n) + \frac{r^2}{2}(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n) ]

多出来的两项相当于人为扩散项，能压制一部分由中心差分激发的振荡，但二维Lax-Wendroff格式的稳定性条件比一维更苛刻，理论给出的稳定区间大约在 (\gamma \le 0.3536)。也就是说\gamma=0.5时它就已经超出稳定域了，只有把\gamma进一步压到0.3左右才勉强收敛，源程序里也附了γ=0.3的验证图。

时间分裂法把二维问题拆成两个连续的一维迎风更新。先在同一时间层内分别对y方向、x方向做半步推进，偶数时间层调换两个方向的顺序，消除人为方向偏差。它的递推关系是：

% 奇数时间层：先y后x U3(i,j) = U1(i,j) - r*(U1(i,j) - U1(i,j-1)); % 半步沿y方向 U2(i,j) = U3(i,j) - r*(U3(i,j) - U3(i-1,j)); % 半步沿x方向 % 偶数时间层：先x后y U3(i,j) = U1(i,j) - r*(U1(i,j) - U1(i-1,j)); % 半步沿x方向 U2(i,j) = U3(i,j) - r*(U3(i,j) - U3(i,j-1)); % 半步沿y方向

这里要提一下时间分裂法的边界处理。每个半步结束时，对出流边界做一次线性外插；两个半步用的是不同的差分方向，所以外插的边界也不同，不能整个时间步结束只做一次。资源里的代码在奇数层和偶数层分别写了外插语句，就是这个原因。

2.4 四种格式的稳定性条件与γ取值

γ本质上是 (\Delta t/\Delta x)，也叫CFL数。二维情况下要同时约束x、y两个方向的传播速度，所以稳定性条件往往比一维更紧。资源基于放大因子分析给出的结论如下表：

这个表对调试很有用。跑代码前先确认自己用的格式和γ取值落在哪个区间，再去看波形发没发散，能少走很多弯路。我最初跑Lax-Wendroff格式时，不知道二维稳定条件比一维的0.5更严，守着γ=0.5折腾好久，后来查了资源里的参考文献才知道临界值在0.3536附近。

3. 核心代码拆解：初始化、时间推进与边界处理的完整链路

3.1 从方形波到第一个时间层

源程序一上来做了几件事：清空工作区、加载cbrewer配色、定义网格、生成初值。这里把准备工作单独拿出来看：

clear, clc dx = 1; dy = 1; r = 0.5; % γ = Δt/Δx，要测试1或2时只改这里 x = -16:dx:16; [X, Y] = meshgrid(x); U1 = zeros(size(X)); % 当前时间层 U2 = zeros(size(U1)); % 下一时间层的工作数组 for i = 13:21 for j = 13:21 U1(i,j) = 1; end end

U2被提前初始化成和U1相同的尺寸，这是好习惯。时间推进时只往U2里填数，不重新分配内存，在网格变大时能省掉反复开辟矩阵的开销。(\gamma)的值单独用一个变量r控制，后面所有格式更新公式都引用它，切换0.5、1、2就只要改这一行。

3.2 主时间推进循环与格式切换

资源把C格式、A格式、Lax-Wendroff格式写在同一个脚本里，用注释切换。核心循环如下：

for n = 1:40 % 外层时间步 for i = 2:size(X,1) % 行循环，从2开始是为了避开i-1越界 for j = 2:size(Y,2) % 列循环，同样避开j-1 % 迎风C格式 U2(i,j) = U1(i,j) ... - r*(U1(i,j) - U1(i-1,j)) ... - r*(U1(i,j) - U1(i,j-1)); % 下面是A格式，和C格式互相注释切换 % U2(i,j) = U1(i,j) ... % - 0.5*r*(U1(i+1,j) - U1(i-1,j)) ... % - 0.5*r*(U1(i,j+1) - U1(i,j-1)); % 下面是Lax-Wendroff格式 % U2(i,j) = U1(i,j) ... % - 0.5*r*(U1(i+1,j) - U1(i-1,j)) ... % - 0.5*r*(U1(i,j+1) - U1(i,j-1)) ... % + 0.5*r^2*(U1(i+1,j) - 2*U1(i,j) + U1(i-1,j)) ... % + 0.5*r^2*(U1(i,j+1) - 2*U1(i,j) + U1(i,j-1)); end end % 后面是边界处理、绘图与帧保存 end

逐点双重循环是这里最直观的写法，每个网格点的更新都严格对应差分公式。但要留意，C格式的循环从2开始，因为公式里出现 (u_{i-1,j}) 和 (u_{i,j-1})；切到A格式或Lax-Wendroff格式时，出现的是 (u_{i+1,j})、(u_{i,j-1})，循环上界必须改成size(X,1)-1和size(Y,2)-1，否则最后一个索引越界。源程序注释里并不显眼，很容易忽略。

3.3 边界条件：入流置零、出流线性外插

边界处理是整个程序里最值得细读的部分。由于波向x正方向和y正方向传播，左边界和下边界是入流边界，直接置零：

U2(:,1) = 0; % 左边界 U2(1,:) = 0; % 下边界

右边界和上边界是出流边界。如果照抄C格式的“置零”逻辑，波传到边界后会被人为抹掉或者反弹回来。A格式和Lax-Wendroff格式采用线性外插法，假设边界外的值与内部最后两个点保持线性关系：

U2(:,end) = 2*U2(:,end-1) - U2(:,end-2); % 右边界外插 U2(end,:) = 2*U2(end-1,:) - U2(end-2,:); % 上边界外插

用这个式子，边界点不受入流条件的约束，波在边界处会“自然流出”计算域。C格式只有一步迎风，信息来自左下方向，对右边界和上边界没有依赖，所以不需要这套外插；而A格式和LW格式的差分方向包含左右、上下两侧，就必须把出流外插补上。时间分裂法的情况更特殊，每个半步结束都要根据当前半步的差分方向执行一次外插，顺序错一次，波就会在某条边界上堆积。

3.4 三维绘图与动画输出

可视化部分用surf绘制三维曲面，再通过getframe抓取每一帧，最后用movie或VideoWriter输出：

surf(X, Y, U2, 'EdgeColor', [0.0314 0.3176 0.6118]); colormap(RGB); % 如果cbrewer缺失，换成parula或jet axis([-20 20 -20 20 0 1]); % 固定z轴范围 m(n+1) = getframe; % 记录当前帧 % 需要保存成视频时，取消下面几行注释 % v = VideoWriter('C_upwind_0.5.mp4','MPEG-4'); % v.FrameRate = 15; % v.Quality = 90; % open(v); writeVideo(v,m); close(v);

axis的z轴固定在[0,1]很重要。如果不固定，波形一旦发散，曲面会瞬间冲出量程，y轴和z轴自动缩放后，看似“涨起来”的波反而掩盖了数值振荡的真实形态。固定坐标轴后，发散与否一眼就能看出来。

4. 避坑指南：从越界报错到“看起来正常的发散”

4.1 一运行就报“索引超出数组维度”

现象：把C格式注释掉、换成A格式之后，循环刚跑几步就报Index exceeds array bounds。

原因：A格式使用 (u_{i+1,j}) 和 (u_{i,j+1})，如果循环范围沿用C格式的2:size(X,1)，当i到达最后一行时 (i+1) 必定越界。

解决：切换格式时同步调整循环边界，A格式和Lax-Wendroff格式要从1:size(X,1)-1和1:size(Y,2)-1循环。代码层面看不出不对，但一跑就爆，这是最容易排查也最容易犯的错。

4.2 γ=1时C格式发散，老以为是程序写错了

现象：γ取0.5时波形正常，改成1之后，前20个时间步还算正常，后面波形开始出现锯齿状振荡，峰值迅速超过初值。

原因：这不是代码bug，而是格式的稳定性边界。迎风C格式二维稳定条件是 (\gamma\le0.5)，γ=1已经超出稳定域，误差被逐层放大。A格式更极端，无条件不稳定，γ=0.5也照样发散。

解决：先按资源里的结论把γ调回0.5，确认程序本身没有问题，再切换到γ=1、2去观察发散形态。这样既能看到数值现象，也不会怀疑自己的代码写错了。

4.3 右边界和上边界出现“波反弹”

现象：波推进到右边界或上边界后，没有流出计算域，反而在边界附近出现明显驻波或回波，影响内部区域的波形。

原因：出流边界被错误地置零了。置零相当于把边界外条件设成u=0，波到达边界时遇到一个“壁面”，数值上表现为反射。

解决：对A格式和Lax-Wendroff格式，要补上线性外插U2(:,end)=2*U2(:,end-1)-U2(:,end-2)和U2(end,:)=2*U2(end-1,:)-U2(end-2,:)。时间分裂法则在每个半步结束后分别执行对应方向的外插，不要等到整个时间步结束再一起处理。

4.4 cbrewer报错：未定义函数或变量

现象：脚本跑到RGB=cbrewer(...)时报错'cbrewer' is unrecognized。

原因：cbrewer是一套开源的ColorBrewer配色函数，并不在MATLAB内置工具箱里。资源作者用它来生成更柔和的曲面颜色，换到没有安装该函数的机器上自然报错。

解决：删掉或注释掉cbrewer相关两行，改回内置的colormap(parula)、colormap(jet)或colormap(jet)。它只影响配色，不影响数值结果。

4.5 方波越跑越矮、越跑越宽，是不是精度不够

现象：同样用C格式γ=0.5，方波推进到90个时间步以后，峰值从1降到0.6左右，波形边缘明显变缓。

原因：一阶迎风格式具有数值耗散性，相当于在方程里隐式引入了一个扩散项，耗散大小和 (\gamma)、网格尺寸有关。波形展宽、峰值降低是格式的固有行为，不是计算错误。

解决：想要更“锐利”的波面，一是把网格加密到 (\Delta x=0.5) 或更小，二是换用Lax-Wendroff这类二阶格式，但要注意后者稳定条件更苛刻，需要同时把γ压小、把时间步缩短。

5. 把代码改造成自己的CFD小实验台

5.1 差分格式函数化：一种格式一个接口

原脚本把格式切换写成注释，适合跑单一任务，但要做多格式对比时，反复改注释比较麻烦。我一般会把每种格式封装成独立函数，输入当前解和时间步比值，返回下一时间层的解：

function U2 = advanceUpwindC(U1, r) [nx, ny] = size(U1); U2 = zeros(nx, ny); U2(:,1) = 0; % 入流边界 U2(1,:) = 0; for i = 2:nx for j = 2:ny U2(i,j) = U1(i,j) ... - r*(U1(i,j)-U1(i-1,j)) ... - r*(U1(i,j)-U1(i,j-1)); end end end

主程序里只要写一行U2 = advanceUpwindC(U1, r)就够了，想换格式就换函数名。

5.2 加一个发散监视器，不让程序白跑

CFD脚本常遇到的尴尬是：循环跑了几百步，最后画图才发现早就在第几十步发散了。我通常会维护一个历史峰值数组：

history(n) = max(max(U2)); if history(n) > 1e4 fprintf('发散：第%d步峰值已超出阈值\n', n); break; end

这样一旦数值爆炸，程序立刻停下来并输出位置，不用傻等完整时间步。判断阈值可以按初值的量级设定，比如初值峰值是1，阈值取1e4比较充裕。

5.3 从逐点循环到向量化的性能优化

网格只有 (33\times33) 时，双重循环完全没问题；但把计算域扩展到 (129\times129) 或更大，逐点循环的速度就比较慢了。向量化写法可以直接替换内层循环：

U2(2:end, 2:end) = U1(2:end, 2:end) ... - r*(U1(2:end, 2:end) - U1(1:end-1, 2:end)) ... - r*(U1(2:end, 2:end) - U1(2:end, 1:end-1));

这段代码和迎风C格式的逐点循环数学上完全等价，但运行时速度能提升一个数量级。要特别小心行和列的切片方向，MATLAB矩阵第一维是行、对应y方向，切片错位会让波跑歪。我改造过一版，因为把行列切片写反，波形沿对角线移动，当时盯着surf图看了半天。后来每改完一次格式，都要先跑十个时间层，对比max(max(U2))和逐点循环的结果，一致了再继续。从那以后我每次做CFL数值实验都会强制走一遍这个流程，先函数化、再监视、再向量化，步步验证。希望帮到你。

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