★ FEATURED ARTICLE 简介这是一份面向计算流体力学初学者与高校课程复习者的Matlab上机实习报告围绕二维对流方程在迎风C格式、A格式、Lax-Wendroff格式及时间分裂法下的数值求解展开帮助读者理解差分格式构造、边界条件处理与稳定性判定。整份资料仅含1个PDF文件压缩包约25.94MB内含问题描述、初值条件、实验原理、运行结果、结果分析与附录程序等完整章节可直接对照复现。作者按ΔxΔy1、γ取0.5、1、2分别计算90个时间步长给出四种格式的u-x-y波形图与动画素材并借助误差传播方程推导放大因子解释C格式仅在γ0.5时稳定、其余发散的结论。附带的Matlab程序设计与实现细节适合作为CFD课程作业参考、格式对比笔记与编程排错依据已有319人学习。1. 计算流体力学里的一个反直觉同一个 γ两种命运计算流体力学CFD里最像玄学的一件事是同一个参数 γ Δt/Δx 在不同格式下判若两人。做二维对流方程 ∂u/∂t ∂u/∂x ∂u/∂y 0 的数值解时迎风 C 格式在 γ1 已经炸成一片时间分裂法却在 γ1 给出几乎不衰减的方波平移而 A 格式无论 γ 取 0.5 还是 2理论上一开始就不该指望它收敛。同一道题、同一张 Matlab 网格结论完全相反原因不在代码写得好不好而在格式本身的放大因子长什么样。这张地图大致是这样先用差分把偏导换成网格上的加减法再用 von Neumann 放大因子判断这个换法靠不靠谱最后用边界条件和可视化把结论坐实。它适合两类人刚做完 CFD 上机、只会照着公式敲代码却说不清为什么发散的以及复习计算流体力学、想把四种格式串成一条线的。2. 方程与离散骨架33×33 网格上的一块方波二维对流方程的四种差分格式看起来公式一大堆其实骨架只有三层控制方程给物理网格和初值给舞台边界条件给“从哪进、从哪出”。这三层任何一层含糊后面的稳定性分析就是在解释一个错误实验。梳理下来会发现真正决定你能不能复现出报告里那几张曲面的不是公式背得准不准而是索引、步长换算和边界位置有没有对齐。2.1 控制方程、初值与网格把偏导搬到格点上要算的方程是标量二维对流方程两个方向的波速都取 1也就是信息沿 x 正方向和 y 正方向同时传播∂u/∂t ∂u/∂x ∂u/∂y 0−16 ≤ x ≤ 16−16 ≤ y ≤ 16t 0初值是一个位于中心的方波在 −4 ≤ x ≤ 4、−4 ≤ y ≤ 4 内 u1其余为 0。空间步长固定为 Δx Δy 1于是每条边上有 33 个节点网格是 33×33方波占据的索引区间是 13:21因为从 −16 数起第 13 个节点正好落在 −4 上。这个换算关系一定要自己数一遍索引错一位后面所有图都是错的。要特别注意 γ 的定义。报告里把 γ 记成 Δt/Δx Δt/Δy由于 Δx Δy 1γ 在数值上就等于时间步长 Δt。所以“γ 取 0.5、1、2”实际是说时间步长取 0.5、1、2考试里如果 Δx 不是 1就绝不能把 γ 和时间步长混为一谈这是最常见的换算翻车点。计算总长按要求跑到 90 个时间步取帧和绘图频率则另算。提示Δx 固定为 1 只是教学网格的方便选择稳定性判据里的 γ 永远是 Δt/Δx 这个无量纲比值换网格后必须重新算。2.2 从偏导到差分把更新式看成一张账单四种格式的差别本质上是对 ∂u/∂x 和 ∂u/∂y 用几个点去近似。迎风 C 格式只取上游一点A 格式取左右两点求中心差Lax-Wendroff 在中心差基础上再加一个二阶修正项时间分裂法则是把二维算子拆成两个一维算子分半步走。写成代码循环骨架完全一样只有中间那行更新式不同dx 1; dy 1; % 空间步长网格共 33×33 个节点 r 0.5; % r 即 gamma dt/dx dt/dy x -16:dx:16; [X,Y] meshgrid(x); % 生成二维坐标网格 U1 zeros(size(X)); % n 时刻的解 U2 zeros(size(X)); % n1 时刻的解 U1(13:21,13:21) 1; % 初值方波x,y 均落在 [-4,4] for n 1:90 for i 2:size(U1,1) % 从 2 开始第 1 行/列是入流边界 for j 2:size(U1,2) U2(i,j) U1(i,j) - r*(U1(i,j)-U1(i-1,j)) ... - r*(U1(i,j)-U1(i,j-1)); % 迎风 C 格式 end end U2(:,1) 0; U2(1,:) 0; % 入流边界左列、下行给定为 0 U1 U2; % 层间推进U2 整体覆盖 U1 end这段代码里三个参数值得逐个确认size(U1,1)是行数对应 ysize(U1,2)是列数对应 xU1(13:21,13:21)的先行后列顺序不能反r一旦改成 1 或 2同一段代码就跑出完全不同的结果。U1 U2是整块替换而不是累加写错成U1 U1 U2会得到一条指数增长的曲线看起来像“格式发散”其实是内存复用写错了。2.3 边界条件入流给值出流外插两个方向的波速都大于 0意味着信息只从左边界和下边界进入计算域从右边界和上边界离开。所以左边界第 1 列和下边界第 1 行是入流边界直接给 u 0右边界和上边界是出流边界不能硬给固定值否则波打上去会反射回来把内部的解污染掉。迎风 C 格式只用到 i−1、j−1出流边界天然不需要额外处理但 A 格式和 Lax-Wendroff 格式用到了 i1、j1必须在出流侧补一个条件报告里用的是线性外插u(i1,j) 2u(i,j) − u(i−1,j)u(i,j1) 2u(i,j) − u(i,j−1)它的含义是“边界外面那一层按内部最后两点的斜率线性延伸出去”等于告诉格式波继续往外走不要回头。这段逻辑写成代码只有两行但它是 A 格式和 Lax-Wendroff 能不能跑出“边界干净”的曲面的关键。U2(:,end) 2*U2(:,end-1) - U2(:,end-2); % 右边界x 方向出流线性外插 U2(end,:) 2*U2(end-1,:) - U2(end-2,:); % 上边界y 方向出流线性外插 U2(:,1) 0; % 左边界入流 U2(1,:) 0; % 下边界入流顺序上建议先算外插、再置入流否则外插会把刚置零的入流点当成内部点参与计算。这两组语句谁先谁后在 γ 较小时差别肉眼看不出来等 γ 接近稳定上限时就会变成边界上一条明显翘起来的边。3. 四种格式对比迎风 C、A、Lax-Wendroff 与时间分裂同一道题摆四种格式比较的维度其实只有三个精度阶数、数值耗散和稳定域。精度决定波峰会不会被削平耗散决定削得多快稳定域决定你能不能活着算到第 90 步。中间这两章的重点不在于把公式抄全而在于看到一条公式就能预判它会画出什么形状的曲面——这件事在考场上比会推导更值钱在机房调试时更是救命。3.1 迎风 C 格式用一阶精度换来的安全感迎风 C 格式的更新式是u(n1,i,j) u(n,i,j) − γ[u(n,i,j) − u(n,i−1,j)] − γ[u(n,i,j) − u(n,i,j−1)]它做的事情非常朴素把 x 方向的增量和 y 方向的增量分别扣掉。因为只取上游一点它天然“顺着波走”稳定条件是四个格式里最宽松的之一γ ≤ 0.5。代价是一阶精度带来的强数值耗散。方波的四个角会被磨圆波峰高度随时间下降看起来像被热传导抹过一遍。用修改方程的话说一维迎风等价于在原方程右边加了一项 (1−γ)aΔx/2 · u_xx 的数值粘性γ 越接近 1这个粘性越小γ 1 时一维迎风是精确平移一点不糊。但二维迎风 C 是两个方向耗散的叠加稳定上限被压到 0.5所以二维里根本到不了“零耗散”那个甜点这也是为什么报告里 γ0.5 的曲面已经明显掉高度、边界也不再锐利。复现时最该盯住的现象是γ0.5 时方波整体向右上平移、峰高缓慢下降、四个角变钝γ1 和 γ2 时数值迅速抬高并出现正负振荡几十步之内就冲出坐标范围。3.2 A 格式中心差分漂亮但无条件不稳定A 格式中心差分的更新式把上游差值换成了左右两点的中心差u(n1,i,j) u(n,i,j) − (γ/2)[u(n,i1,j) − u(n,i−1,j)] − (γ/2)[u(n,i,j1) − u(n,i,j−1)]从截断误差看它是二阶精度比迎风 C 漂亮得多对称、无偏、不引入方向性耗散。但实验里 γ 取 0.5、1、2 三个值全部发散这不是调参能救的而是格式本身的缺陷。原因在放大因子上G 1 − iγ(sinβ sinθ)取模得 |G| sqrt(1 γ²(sinβsinθ)²)只要 γ0 且波数不为零|G| 就严格大于 1。也就是说中心差分对任何非零波数都在放大误差网格越细、初值越光滑它照样放大只是放大得慢一点。这类格式在数值上毫无“后悔药”可吃。一个很实用的观察点A 格式跑出来的图不是简单地“变成 NaN”而是先出现高频锯齿状振荡再整体抬高。如果只看前几步的图觉得“还行”很容易误判成稳定等跑到第 40 步才发现已经飞了。3.3 Lax-Wendroff二阶精度的代价是 0.3536Lax-Wendroff 格式在中心差的基础上补了一个二阶项u(n1,i,j) u(n,i,j) − (γ/2)[u(n,i1,j) − u(n,i−1,j)] − (γ/2)[u(n,i,j1) − u(n,i,j−1)] (γ²/2)[u(n,i1,j) − 2u(n,i,j) u(n,i−1,j)] (γ²/2)[u(n,i,j1) − 2u(n,i,j) u(n,i,j−1)]这一项可以理解为在时间方向做了二阶展开所以它既是二阶精度又比 A 格式稳。但二维的稳定域比一维窄得多。按 K1 Δt/Δx、K2 Δt/Δy 的定义二维 Lax-Wendroff 的判据是K1^(2/3) K2^(2/3) ≤ 1当 K1 K2 γ 时条件变成 2γ^(2/3) ≤ 1解出 γ ≤ 2^(−3/2) ≈ 0.3536。这就解释了为什么 γ0.5、1、2 三个值全部发散——它们都在这条线之上而补做的 γ0.3 实验就收敛了。注意一维 Lax-Wendroff 的稳定上限是 γ ≤ 1把同一套代码直接搬到二维稳定域会缩到 0.3536这是最容易踩的坑代码没改错是判据换了。收敛之后也不能高兴太早。Lax-Wendroff 没有迎风那种耗散方波边缘会出现明显的色散波纹峰前后各拖出一串小幅振荡这是二阶格式的典型痕迹不是程序写错了。3.4 时间分裂法二维拆成两个一维γ1 反而最准时间分裂法把二维问题拆成两步一维问题每个半步只处理一个方向u(n1/2,i,j) u(n,i,j) − γ[u(n,i,j) − u(n,i,j−1)] u(n1,i,j) u(n1/2,i,j) − γ[u(n1/2,i,j) − u(n1/2,i−1,j)]代码上它多了一个中间变量 U3 存放半时间步的解每步先算一个方向用半时间步的结果再算另一个方向。实验里 γ0.5 和 γ1 都收敛γ2 发散稳定条件是 γ ≤ 1。最反直觉的地方在这里γ1 时一维迎风格式恰好退化成精确平移 u(n1,i) u(n,i−1)耗散为零两个方向各平移一格整体就是把方波沿对角线挪一步形状几乎不变。所以时间分裂法在 γ1 反而给出了四种格式里最干净的曲面而不是像迎风 C 那样在 γ1 直接炸掉。这条结论几乎每年都会有人在复习时记反。还有一个实现细节值得抄进笔记为了避免固定按“先 y 后 x”的顺序造成方向上的系统偏差程序在奇数步用 y→x 顺序、偶数步用 x→y 顺序。这一步不影响稳定性判据放大因子是两个一维因子相乘交换顺序不改变模但会影响图形的对称性尤其是方波的对角线走向。4. 稳定性分析von Neumann 放大因子怎么一步步算稳定性这件事实验图像只能告诉你“炸了”或“没炸”说不清临界值在哪。von Neumann 方法的价值在于它把“会不会炸”变成一个可以手算的不等式让你在写代码之前就知道 γ 该取多少而不是跑完 90 步再回头猜。的推导套路只有三步四种格式走的是同一条流水线区别只在代入后整理出的 G 长什么样。4.1 误差传播方程与三步法设数值解和真解之间的误差为 ε由于方程是线性的误差本身满足和原差分格式完全相同的方程。以迎风 C 格式为例误差传播式就是把 u 换成 ε 的那条式子ε(n1,i,j) ε(n,i,j) − γ[ε(n,i,j) − ε(n,i−1,j)] − γ[ε(n,i,j) − ε(n,i,j−1)]第一步设误差按单一波数传播ε(n,i,j) A^n · exp[i(k1·x k2·y)]。第二步把相邻点的误差用同一个 A^n 表示例如 ε(n,i−1,j) A^n·exp[i(k1(x−Δx) k2·y)]这样所有项都能约掉空间因子。第三步剩下的是层间放大因子 G A^(n1)/A^n它是一个复数稳定性判据就是 |G| ≤ 1 对所有波数成立。记 β k1Δx、θ k2Δy把指数用欧拉公式拆成实部和虚部迎风 C 的放大因子整理成G 1 − 2γ γ(e^(−iβ) e^(−iθ))这一步是整个推导里最容易写错的地方符号漏一个负号后面所有结论都会反。建议整理完先代两个极端波数自检——β θ 0零波数应该 |G| 1和 β θ π最细锯齿波通常给出最坏情况。想验证自己的推导对不对不用肉眼看公式直接拿几行 Matlab 把 |G| 的最大值扫出来就行beta linspace(-pi, pi, 201); theta linspace(-pi, pi, 201); [B,T] meshgrid(beta, theta); % 覆盖所有波数组合 glist 0:0.05:1.2; for k 1:numel(glist) r glist(k); G 1 - r*(1-exp(-1i*B)) - r*(1-exp(-1i*T)); % 迎风 C 的放大因子 amax(k) max(abs(G(:))); % 全波数上的最大模 end plot(glist, amax, glist, ones(size(glist))); % 与 |G|1 的水平线比较图上曲线第一次触到 1 的那个 γ就是稳定上限迎风 C 会给出 0.5。把G ...那一行换成 A 格式的1 - 1i*r*(sin(B)sin(T))曲线会整体在 1 之上直观看到“无条件不稳定”。4.2 四个格式的 |G| 与稳定阈值把四种格式的放大因子并排放稳定性的差别一眼就出来了格式放大因子 G|G| 情况稳定条件迎风 C1 − 2γ γ(e^(−iβ) e^(−iθ))锯齿波 βθπ 时 G 1−4γγ ≤ 0.5A中心差1 − iγ(sinβ sinθ)sqrt(1 γ²(sinβsinθ)²) ≥ 1无条件不稳定Lax-Wendroff二阶展开模较复杂二维稳定域收窄K1^(2/3)K2^(2/3) ≤ 1等步长下 γ ≤ 0.3536时间分裂 C[1−γ(1−e^(−iβ))]·[1−γ(1−e^(−iθ))]两个一维因子模都不超过 1γ ≤ 1迎风 C 的临界值这样得到取最坏的锯齿波 β θ πG 1 − 4γ要求 |1 − 4γ| ≤ 1得 0 ≤ γ ≤ 0.5。γ 0.5 时进一步可验证 |G|² 0.5[1 cos(β−θ)] ≤ 1 对所有波数成立所以 0.5 是既不充分条件里的紧界理论值与实验完全对上。A 格式的结论更干脆|G|² 1 γ²(sinβsinθ)²只要 γ0总能取到 sinβsinθ ≠ 0模严格大于 1无条件不稳定。这也解释了为什么三个 γ 值全都发散——不是取值不够小是根本没有安全区。时间分裂法要分开看两个因子一维迎风因子的模平方是 1 − 2γ(1−γ)(1−cosβ)只要 0 ≤ γ ≤ 1 就不超过 1两个因子相乘模也不会超过 1。所以 γ ≤ 1 是它的稳定界γ2 时单个因子的模已经大于 1乘起来发散与实验一致。顺带解释了 γ1 的“零耗散”此时 1−cosβ 的系数 2γ(1−γ) 恰好为零|G| 恒等于 1。4.3 γ0.5/1/2 的实验现象与理论对照把报告里的图像现象和理论判据并排放是复习时最省力的一张表格式γ0.5γ1γ2理论临界值迎风 C稳定峰被抹平发散发散0.5A中心差发散发散发散无无条件不稳定Lax-Wendroff发散发散发散0.3536补测 γ0.3 收敛时间分裂 C稳定稳定且几乎无耗散发散1这张表里有两处最容易记反Lax-Wendroff 在 γ0.5 是发散的虽然 0.5 听上去“不大”时间分裂法在 γ1 是稳定的而且是最准的一档。判断依据永远回到 |G| ≤ 1而不是“γ 越小越安全”这种直觉——γ 太小反而会放大迎风格式的耗散让方波糊得更快。5. 易漏点实验报告里最容易翻车的五个地方同一份代码有人跑出干净的曲面有人跑出一团噪声差别往往不在格式而在这五个细节。它们有个共同特点程序不报错、Matlab 不警告、图还能画出来所以极难自查。下面按“错在哪 → 为什么容易漏 → 怎么补”逐条过每条都是可以直接照着改的。5.1 入流边界到底该置哪一行易漏点把U2(end,:) 0当成入流边界条件写上而end对应的是 y 最大的那一行也就是出流侧。为什么漏入流和出流在公式里只是“左/右”“上/下”的区别写代码时1和end长得同样顺手而且置零之后程序照跑图也照出只是边界附近多了一条奇怪的平台容易被当成数值耗散。怎么补波速为正时入流永远是第 1 列和第 1 行写成U2(:,1)0; U2(1,:)0;出流永远是end那一侧用线性外插。改完后对着 γ0.5 的曲面看四个边入流侧应该是一条干干净净的零线出流侧应该看到波“走出去”而不是堆在边上。5.2 出流外插的 end-2 写成 end-1易漏点外插式写成2*U2(end-1,:) - U2(end-1,:)两个点写成同一个点。为什么漏这样写等价于U2(end-1,:)也就是“边界值等于里面一格”零梯度边界看上去也合理甚至还更稳只有在波传到边界时才会发现反射比预期明显但那时候已经很难判断是格式问题还是边界问题。怎么补严格按 u(i1) 2u(i) − u(i−1) 写第二项必须是end-2两个方向的边界各检查一遍。一个笨但有效的自查方法临时把外插换成U2(end,:) U2(end-1,:)跑一次如果两张图几乎一样说明你的波还没走到边界这时候不能算“边界条件验证通过”。5.3 循环从 2 开始却漏掉第 1 行的更新易漏点内层循环写for i 2:size(U1,1)于是 i1 和 j1 的整行整列永远保持初值或上一层的旧值忘了每步都要重新赋 0。为什么漏初值里这两行本来就是 0前几十步看不出异常但入流边界旁边那一圈如果不每步刷新等于边界条件只在第 0 步生效后面全靠旧数据撑着。怎么补把边界赋值放在时间循环内部、内层循环之后每步都执行一次顺序上先外插后置零。想确认有没有漏把U2(:,1)单独打印出来看是否恒为 0比盯着曲面图可靠得多。5.4 时间分裂法的半时间步变量与奇偶顺序易漏点两个半步共用一个数组或者忘了在奇数步和偶数步交换 x/y 的先后顺序。为什么漏半步的结果如果直接覆盖原数组第二步就会拿已经被修改过的邻居去算等于把两个方向耦合在一起格式悄悄变成了另一种东西稳定性判据也就不再适用。怎么补用三个数组分工——U1 是整步初值U3 存半时间步结果U2 存整步结果每步结束再U1 U2。顺序上用mod(n,2)判断奇偶奇数步先 y 后 x、偶数步先 x 后 y并在两个半步各自补上出流外插。这样跑出来的 γ1 结果应该是一个几乎不失真的对角平移方波如果峰高掉得厉害先查这一步。5.5 把 0.3536 当成拍脑袋的数易漏点Lax-Wendroff 直接用一维的 γ ≤ 1 去设参数跑到 γ0.5 发现发散怀疑代码有 bug。为什么漏一维和二维的稳定域不是同一个判据二维多了一个方向的波数耦合条件从线性收紧成 2/3 次幂形式数值上从 1 掉到 0.3536差了三倍。怎么补把判据 K1^(2/3) K2^(2/3) ≤ 1 写在代码注释里等步长时直接记 γ ≤ 2^(−3/2) ≈ 0.3536调参时先取 0.3 确认格式本身没问题再往临界值靠。凡是“理论说稳、实验说炸”的情况第一反应都应该是回去核对判据的维度而不是重写循环。6. 复习清单从上机报告到考场推导的自查路径把这份材料当复习提纲用最有效的顺序不是从头读公式而是拿一张白纸按“方程 → 网格 → 格式 → 放大因子 → 临界值 → 现象”六步默写一遍每一步都问自己能不能说出下一层为什么这么定。下面这张表是我自己复盘时用的检查项左边是动作中间是判据右边是出错时的典型症状对着走一遍基本能定位到问题在哪一层。检查项判据出错时的典型表现网格与初值对齐33×33方波在 13:21方波偏离中心或尺寸不对γ 与 Δt 换算γ Δt/ΔxΔx1 时数值相等换网格后稳定性结论全部失效入流/出流位置第 1 行列为入流end 侧外插边界堆波、反射明显迎风 C 临界值γ ≤ 0.5锯齿波给出γ1 时误以为“应该稳”分裂法临界值γ ≤ 1γ1 最精确记成“γ 越小越准”LW 二维判据K1^(2/3)K2^(2/3) ≤ 1用一维的 γ ≤ 1 设参数动画取帧固定 z 轴范围帧率一致发散时图像自动缩放看不出爆炸默写完公式再用代码验一遍比纯看公式记得牢。最小验证脚本不用长把放大因子扫一遍、把边界打印出来看两件事就够% 1) 扫出迎风 C 的稳定上限应停在 0.5 附近 r 0.5; [B,T] meshgrid(linspace(-pi,pi,101)); G 1 - r*(1-exp(-1i*B)) - r*(1-exp(-1i*T)); fprintf(gamma%.2f, max|G|%.4f\n, r, max(abs(G(:)))); % 2) 确认入流边界每步都被刷新 assert(all(U2(:,1) 0) all(U2(1,:) 0), 入流边界没有每步置零);第一段把 γ 换成 0.6、0.8、1.0 各跑一次max|G| 会依次越过 1临界点就“长”在手上了第二段的断言放在时间循环里能第一时间抓住 5.3 里那个“只在第 0 步生效”的边界。这两件事加起来不超过二十行却比反复调参省掉一晚上。还有一个查漏补缺的笨办法把四种格式的曲面按 γ0.5、1、2 摆成 3×4 的九宫格Lax-Wendroff 补一张 γ0.3不看标题只看形状先自己判断哪个稳、哪个炸再翻回第 4 章对答案。能全部判对说明放大因子那套推导是真的懂了而不是背下了四行结论。我自己的教训是第一次做这个上机时把出流外插的end-2写成了end-1图看着挺正常直到答辩前被问“边界为什么有反射”才翻出来白跑了两天的动画。从那以后我每次写二维循环都强制先把四个边界各打印一行、确认每步都在刷新再开始调 γ。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
