简介计算流体力学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 y1)所以得到 (33\times33) 的网格点。初值在 (-4\le x\le4,\ -4\le y\le4) 的方形区域内取 (u1)其余位置取 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}^{n1} 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}^{n1} u_{i,j}^n - \frac{\Delta t}{2\Delta x}(u_{i1,j}^n-u_{i-1,j}^n) - \frac{\Delta t}{2\Delta y}(u_{i,j1}^n-u_{i,j-1}^n) ]它不区分方向直接取左右对称差分精度比迎风格式高但在二维条件下无条件不稳定放大因子的模恒大于等于1所以源程序里三种γ取值全部发散。运行结果里A格式的波很快就会变成高频振荡继而溢出成NaN。Lax-Wendroff格式在中心差分基础上补了二阶耗散项[ u_{i,j}^{n1} u_{i,j}^n - \frac{r}{2}(u_{i1,j}^n-u_{i-1,j}^n) - \frac{r}{2}(u_{i,j1}^n-u_{i,j-1}^n) \frac{r^2}{2}(u_{i1,j}^n-2u_{i,j}^nu_{i-1,j}^n) \frac{r^2}{2}(u_{i,j1}^n-2u_{i,j}^nu_{i,j-1}^n) ]多出来的两项相当于人为扩散项能压制一部分由中心差分激发的振荡但二维Lax-Wendroff格式的稳定性条件比一维更苛刻理论给出的稳定区间大约在 (\gamma \le 0.3536)。也就是说\gamma0.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两个方向的传播速度所以稳定性条件往往比一维更紧。资源基于放大因子分析给出的结论如下表差分格式稳定性分析关键结论理论稳定条件实测结果迎风C格式放大因子模在特定相位取得最大值(\gamma \le 0.5)γ0.5收敛γ1、2发散A格式(G\sqrt{1\gamma^2(\sin\beta\sin\theta)^2} )Lax-Wendroff二维条件下修正的稳定性判据(\gamma \le 0.3536)γ0.5、1、2发散γ0.3收敛时间分裂法C格式两个半步放大因子相乘后取最大模(\gamma \le 1)γ0.5、1收敛γ2发散这个表对调试很有用。跑代码前先确认自己用的格式和γ取值落在哪个区间再去看波形发没发散能少走很多弯路。我最初跑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 endU2被提前初始化成和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(i1,j) - U1(i-1,j)) ... % - 0.5*r*(U1(i,j1) - U1(i,j-1)); % 下面是Lax-Wendroff格式 % U2(i,j) U1(i,j) ... % - 0.5*r*(U1(i1,j) - U1(i-1,j)) ... % - 0.5*r*(U1(i,j1) - U1(i,j-1)) ... % 0.5*r^2*(U1(i1,j) - 2*U1(i,j) U1(i-1,j)) ... % 0.5*r^2*(U1(i,j1) - 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_{i1,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(n1) 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_{i1,j}) 和 (u_{i,j1})如果循环范围沿用C格式的2:size(X,1)当i到达最后一行时 (i1) 必定越界。解决切换格式时同步调整循环边界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 右边界和上边界出现“波反弹”现象波推进到右边界或上边界后没有流出计算域反而在边界附近出现明显驻波或回波影响内部区域的波形。原因出流边界被错误地置零了。置零相当于把边界外条件设成u0波到达边界时遇到一个“壁面”数值上表现为反射。解决对A格式和Lax-Wendroff格式要补上线性外插U2(:,end)2*U2(:,end-1)-U2(:,end-2)和U2(end,:)2*U2(end-1,:)-U2(end-2,:)。时间分裂法则在每个半步结束后分别执行对应方向的外插不要等到整个时间步结束再一起处理。4.4 cbrewer报错未定义函数或变量现象脚本跑到RGBcbrewer(...)时报错cbrewer is unrecognized。原因cbrewer是一套开源的ColorBrewer配色函数并不在MATLAB内置工具箱里。资源作者用它来生成更柔和的曲面颜色换到没有安装该函数的机器上自然报错。解决删掉或注释掉cbrewer相关两行改回内置的colormap(parula)、colormap(jet)或colormap(jet)。它只影响配色不影响数值结果。4.5 方波越跑越矮、越跑越宽是不是精度不够现象同样用C格式γ0.5方波推进到90个时间步以后峰值从1降到0.6左右波形边缘明显变缓。原因一阶迎风格式具有数值耗散性相当于在方程里隐式引入了一个扩散项耗散大小和 (\gamma)、网格尺寸有关。波形展宽、峰值降低是格式的固有行为不是计算错误。解决想要更“锐利”的波面一是把网格加密到 (\Delta x0.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数值实验都会强制走一遍这个流程先函数化、再监视、再向量化步步验证。希望帮到你。本文还有配套的精品资源点击获取