简介这份PDF面向学习计算流体力学与数值格式的本科生、研究生及Matlab实践者围绕二维对流方程在迎风C格式、A格式、Lax-Wendroff格式和时间分裂法下的数值求解展开重点处理∆t/∆x与∆t/∆y取0.5、1、2时的稳定性讨论并给出u-x-y曲面图与动画制作思路。压缩包仅含1个PDF文件大小约25.94MB正文按问题描述、实验原理、差分格式、边界条件、运行结果、结果分析和附录程序组织便于顺着报告脉络复现计算流程。已有319人学习适合作为CFD上机实习或课程设计的参考。读者可从中获得四种格式的离散表达式、入流出流边界处理方式、不同Courant数下的结果对比与稳定性分析以及可直接借鉴的Matlab程序结构和绘图方法用于核对自身代码与实验现象。1. 同一个方程、四种写法、三种结局稳定性不是玄学二维对流方程 ∂u/∂t ∂u/∂x ∂u/∂y 0初值是一个中心 8×8 的方形平台网格 Δx Δy 1计算域 [-16,16]²推进 90 步。这是计算流体力学里最简单的一个 CFD 算例Matlab 三十行就能写完。可一旦把 γ Δt/Δx 分别取 0.5、1、2四种差分格式给出的结果完全不是一回事迎风 C 格式在 0.5 稳住、在 1 和 2 处爆掉A 格式三档全爆Lax-Wendroff 也全爆要把 γ 压到 0.3 才看得见一个像样的波形时间分裂法下的 C 格式却在 0.5 和 1 都收敛直到 2 才散。四个格式写的是同一个方程代码骨架几乎一样差别只在模板取哪几个邻居、系数乘多少。所以「稳定性」从来不是调参调出来的玄学而是格式本身写死了的一条硬边界。搞懂这条边界怎么来的才算真正把 CFD 入门而不是把有限差分当黑匣子调着玩。2. 放大因子 G冯·诺依曼分析如何把稳定性压成一道复平面几何题2.1 误差和真解走同一个差分方程这是整套推导的地基差分格式可以写成 u^{n1} L(u^n) 的形式L 是一个线性算子。设数值解为 U精确解为 u误差 ε U - u。因为方程和格式都是线性的且两者满足同样的边界条件把两个式子相减ε 满足的差分方程和 U 满足的完全一致——只是初值从平台函数换成了初始误差分布。这一步是整个稳定性分析的基础误差不受解答本身影响它会独立地按格式规定的规律演化。既然误差自己会走那么判断一个格式行不行就等价于判断「任意初始误差在反复施加 L 之后会不会被放大到无穷」。如果存在某个误差模态每次迭代都被放大一点点90 步、900 步之后就会盖过真解我们在图上看到的就是数值爆炸而不是波形传播。常见的实现思路是只在常系数、周期边界下讨论这件事因为这时 L 是一个循环矩阵特征函数是离散傅里叶模态误差可以干净地分解成一堆平面波的叠加。这也是冯·诺依曼分析的适用前提绝大多数 CFD 格式的第一轮筛选都靠它。2.2 模态代入把逐点递推换成逐波数缩放取单个误差模态 ε^n_{i,j} A_n · exp(i(k1·x_i k2·y_j))其中 k1、k2 是两个方向的波数A_n 是第 n 步的复振幅。把它代进误差传播方程左右两边 exp 因子会整体约掉剩下的只有振幅之间的比例关系定义这个比例为放大因子 G A_{n1} / A_n。以迎风 C 格式为例误差方程是ε^{n1}{i,j} ε^n{i,j} - γ(ε^n_{i,j} - ε^n_{i-1,j}) - γ(ε^n_{i,j} - ε^n_{i,j-1})代入后ε^n_{i-1,j} ε^n_{i,j}·exp(-i·k1·Δx)ε^n_{i,j-1} ε^n_{i,j}·exp(-i·k2·Δy)。记 β k1·Δx、θ k2·Δy整理得到G 1 - 2γ γ·exp(-iβ) γ·exp(-iθ)β 和 θ 各自扫过 [-π, π]就覆盖了网格能分辨的全部波数。稳定性要求对每一个 (β, θ)|G| ≤ 1 都成立只要存在一组波数让 |G| 1对应模态就会被指数放大格式不稳定。这个结论在常系数周期边界下是充要的作为工程判据足够用。用几行 Matlab 可以把某个格式的 |G| 极值直接扫出来比手推快得多也方便验证推导有没有算错% 扫描全部波数求放大因子的最大模迎风 C 格式 / A 格式 / LW 格式通用 gamma 0.5; % gamma dt/dx dt/dy beta linspace(-pi, pi, 601); theta linspace(-pi, pi, 601); [B, T] meshgrid(beta, theta); G_C 1 - 2*gamma gamma*exp(-1i*B) gamma*exp(-1i*T); % 迎风 C G_A 1 - 1i*gamma*(sin(B) sin(T)); % A 格式 G_LW G_A - 2*gamma^2*(sin(B/2).^2 sin(T/2).^2); % Lax-Wendroff fprintf(C : %.6f\n, max(abs(G_C (:)))); fprintf(A : %.6f\n, max(abs(G_A (:)))); fprintf(LW : %.6f\n, max(abs(G_LW(:))));gamma 就是网格雷诺数意义上的库朗数B、T 是两个方向的无量纲波数。把 gamma 从 0.1 扫到 2打印三条曲线你会看到 C 格式的 |G|max 在 γ 0.5 处从 1 拐头向上A 格式从 γ 稍微大于 0 开始就永远大于 1LW 格式的拐点在 0.3536 附近。三条曲线三种命运后面几章就是把它们各自的原因拆开。2.3 |G|≤1 的代数与几何双重读法把迎风 C 格式的 G 写成矢量叠加的形式最直观G (1 - 2γ, 0) (γ·cosβ, -γ·sinβ) (γ·cosθ, -γ·sinθ)第一个矢量是实轴上长度 1 - 2γ 的线段后两个是模长都为 γ、角度分别由 β 和 θ 决定的旋转矢量。当 β θ 时两个旋转矢量同向它们的水平分量叠加、垂直分量也叠加这是最容易冲出去的方向两个 γ 长的箭头并列甩向左边正好抵消掉实轴段的收缩也最容易越过单位圆。也可以老老实实算一遍。令 β θ、c cosβ、s sinβ|G|² (1 - 2γ 2γc)² (2γs)² 1 - 4γ 8γ² 4γ(1 - 2γ)·cγ ≤ 0.5 时系数 1 - 2γ ≥ 0最大值在 c 1 处取到代入得 |G|²max 1恰好压在单位圆上γ 0.5 时系数变负最大值挪到 c -1|G|²max (1 - 4γ)²模长 |G|max |1 - 4γ|。γ 1 时是 3γ 2 时是 7。这意味着每推进一个时间步最快那个误差模态就放大 7 倍90 步之后是 7⁹⁰ 量级surf 图上只剩几根戳破坐标轴的尖刺其余部分被压成一张平板。注意发散在图上的表现形式常常不是「数值变成 NaN」而是「极少数点数值巨大、其余点看起来还是零」。看到这种形态先别怀疑初值优先怀疑放大因子越界。3. 迎风 C 格式γ≤0.5 这条线究竟卡在哪个方向3.1 依赖域迎风模板把信息流向写进了网格迎风 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})它只取 (i,j)、(i-1,j)、(i,j-1) 三个点。这三个点正好位于 (i,j) 点的上游特征线从西北方向来信息沿 x、y 正方向传播模板就往负方向偏一格。这种「偏」不是随意选择而是把 PDE 的特征方向编码进离散模板让数值依赖域尽量贴近真实依赖域。一维情况下CFL 条件要求特征线在 Δt 内跨过的距离不超过一个网格即 γ ≤ 1。二维时两个方向同时推进Δt 内波在 x 方向走 γ 格、y 方向走 γ 格加起来不能超过一个网格的对角尺度于是条件变成 γ_x γ_y ≤ 1。取 γ_x γ_y γ就得到 γ ≤ 0.5。这跟上一章代数推导的结果完全吻合——两条路同一个答案。迎风 C 格式是一阶精度的代价是数值耗散明显课件里那个方形平台在 γ 0.5 下推进 90 步后边缘被抹圆、峰值被削低看起来像被风吹散的一团雾。但耗散的另一面是稳它老老实实把高频误差压下去不会反过来放大。3.2 入流给值、出流外插边界处理决定了波形能不能走出去边界条件是这类算例里最容易翻车的地方比格式本身还容易。波向 x、y 正方向传播那么左边界i 1和下边界j 1是入流边界物理上信息从外面流进来直接给定值 u 0右边界i end和上边界j end是出流边界信息从域内流出去不能强行给定值否则波撞上墙会反弹回来形成一圈假的回波。常见的处理是线性外插把内部最后两点的斜率延长一格u_{i1,j} 2u_{i,j} - u_{i-1,j} u_{i,j1} 2u_{i,j} - u_{i,j-1}这个式子的含义是「边界外那一格的值按线性趋势推」——等价于二阶的单侧差分允许波平滑地离开计算域而不产生明显反射。它的误差量级比格式本身低一阶不会污染内部解是出流边界最省事也最常用的做法。实现顺序要特别注意内点循环跑完之后再统一补边界而且近边界那一圈内点的循环范围要相应缩小。如果一边循环一边改边界后面的点会用到已经更新过的邻居等价于在格式里混进了半隐式成分稳定性分析就不再适用。3.3 Matlab 实现与三档 γ 的对照下面是迎风 C 格式的主体循环。网格用 meshgrid(x, y) 生成注意此时 X(i,j) x(j)、Y(i,j) y(i)第一维对应 y、第二维对应 x方台的索引区间要和这个约定对齐dx 1; dy 1; gamma 0.5; % gamma dt/dx dt/dy x -16:dx:16; y -16:dy:16; [X, Y] meshgrid(x, y); % X(i,j)x(j), Y(i,j)y(i) U1 zeros(size(X)); U1(13:21, 13:21) 1; % 初值x,y 均落在 [-4,4] 的方形平台 U2 zeros(size(U1)); for n 1:90 for i 2:size(X,1) for j 2:size(X,2) % 迎风 C 格式只用上游的 (i-1,j) 与 (i,j-1) U2(i,j) U1(i,j) - gamma*(U1(i,j) - U1(i-1,j)) ... - gamma*(U1(i,j) - U1(i,j-1)); end end U2(1,:) 0; U2(:,1) 0; % 入流边界 U2(end,:) 2*U2(end-1,:) - U2(end-2,:); % 出流x 向外插 U2(:,end) 2*U2(:,end-1) - U2(:,end-2); % 出流y 向外插 U1 U2; endgamma 是库朗数取 0.5 时正好踩在稳定边界上dx、dy 固定为 1 时它就等于时间步长。索引 13:21 对应 x 从 -16 起算第 13 到第 21 个点即 [-4, 4] 共 9 个网格点。入流边界的赋值放在内点循环之后是为了保证边界值覆盖掉循环里可能算出的无效值出流外插也在同一位置顺序是先 x 后 y两者互不依赖。把 gamma 改成 1 或 2同一个脚本跑出来的结果就是上一章算过的 3ⁿ 和 7ⁿ。这不是步长没调好是格式在 γ 0.5 时对某些波数的 |G| 已经大于 1任何初值、任何边界处理都救不回来。4. A 格式、Lax-Wendroff 与时间分裂γ 从 0 到 1 的稳定性谱系4.1 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})代入模态后两个 exp 相减变成正弦G 1 - i·γ·(sinβ sinθ)实部恒为 1虚部是一个纯虚数于是 |G|² 1 γ²(sinβ sinθ)²。除了 sinβ sinθ 0 的少数模态其余模态的模长全部大于 1γ 越大、偏离越多。γ 0.5、β θ π/2 时|G|² 1 4×0.25 2单步放大 1.414 倍90 步后放大约 10¹¹ 量级。这个结果可以直接在复平面上读懂。把格式写成 G 1 z 的形式z -iγ(sinβ sinθ) 落在虚轴上。显式欧拉的绝对稳定域是以 -1 为圆心、半径为 1 的圆盘它只在原点处与虚轴相切。任何纯虚的 z 都会把 1 z 推出单位圆只有一个点除外。换言之中心差分 显式时间推进这个组合天生就把特征值放在了稳定域之外跟 γ 取多少无关。A 格式的不稳定是无条件的不是「参数没调好」。实践里看到 A 格式在 γ 很小时前几步还算正常那是初值本身光滑、高频误差还没被激发出来的假象。跑够步数误差一定会长出来。4.2 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 格式一模一样后两项是新增的阻尼项形式上就是两个方向的离散拉普拉斯算子。放大因子相应变成G 1 - i·γ(sinβ sinθ) - 2γ²(sin²(β/2) sin²(θ/2))虚部由 -iγ(sinβ sinθ) 贡献随 γ 线性增长阻尼项由 -2γ²(...) 贡献随 γ² 增长。两者在 γ 较小时后半部分压得住前半部分γ 一大会反过来。平衡点在哪就决定了稳定域的上限。二维的情况比一维严格得多。一维 Lax-Wendroff 的稳定条件是 γ ≤ 1但二维方形网格上对角波 β θ 会让两个方向的虚部同相叠加达到单个方向的两倍而阻尼项的叠加跟不上这个速度。按已知的判据二维稳定性条件写成 K1^{2/3} K2^{2/3} ≤ 1K1 K2 γ 时解出 γ ≤ 2^{-3/2} ≈ 0.3536。这个数比 0.5 还小比一维的 1 小了近三倍。数值实验能验证这一点γ 取 0.5、1、2 三档全炸把 γ 降到 0.3波形立刻恢复正常。0.3 0.3536刚好在安全区内。代价也很清楚——同样推进 90 个时间步γ 0.3 意味着物理时间只走了 27而迎风 C 格式在 γ 0.5 时走了 45。想推进到同样的物理时刻Lax-Wendroff 要多花将近一倍的时间步这是二阶精度在二维买的单。4.3 时间分裂法两个一维迎风串联稳定域反而放宽到 γ≤1时间分裂法换了个思路不直接离散二维方程而是把它拆成两个一维方程先沿 y 推进一次再用中间结果沿 x 推进一次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})每一步都是标准的一维迎风格式。因为两步是算子的串联放大因子就是两个一维放大因子的乘积G [1 - γ(1 - e^{-iθ})] · [1 - γ(1 - e^{-iβ})]分别算每个因子的模。以 y 方向那个为例展开得 |·|² 1 - 2γ 2γ² 2γ(1 - γ)cosθ。cosθ 1 时值恒为 1cosθ -1 时值为 (1 - 2γ)²只要 γ ≤ 1 就不超过 1。中间取值更小。所以单个因子在 0 ≤ γ ≤ 1 内满足 |G| ≤ 1乘积自然也不超过 1。稳定域就这样从 0.5 放宽到了 1。原因不复杂二维迎风 C 格式里存在交叉项x、y 的推进同时发生误差在两个方向上的贡献会叠加。分裂法把它们错开了任何时刻只在一个方向上推进交叉项根本不存在只剩两个独立的一维问题。一维迎风的条件是 γ ≤ 1乘积的条件也是 γ ≤ 1。这就是「拆开反而更稳」的全部秘密。实现上有一个细节。既然每步只沿一个方向推进另一步的推进方向在时间上是分开的那么第一步求出的中间量在第二步里会被完整消费掉边界处理也要分两次做。代码如下r 1.0; % 分步推进时的系数取 1 仍稳定 for n 1:90 % 第一步沿 y 方向推进对应 i 方向索引 for i 2:size(X,1) for j 1:size(X,2) U3(i,j) U1(i,j) - r*(U1(i,j) - U1(i-1,j)); end end U3(end,:) 2*U3(end-1,:) - U3(end-2,:); % 中间步出流外插 % 第二步沿 x 方向推进对应 j 方向索引 for i 1:size(X,1) for j 2:size(X,2) U2(i,j) U3(i,j) - r*(U3(i,j) - U3(i,j-1)); end end U2(:,end) 2*U2(:,end-1) - U2(:,end-2); % 整步出流外插 U2(1,:) 0; U2(:,1) 0; % 入流边界 U1 U2; endr 是分步推进的系数理论上每一步用半时间步更严谨取全步时稳定性判据会相应收紧算例里直接取 r 0.5 或 1 都能出结果这正是分裂法宽容度的体现。U3 是中间层注意它只在两步之间存活不要复用成下一轮的初值。两次出流外插必须分开做第一步的外插作用于 y 方向的边界第二步的外插作用于 x 方向混在一起会让两个方向的外插互相污染。分裂法的代价是时间精度。两个一维算子串联交换性误差是一个 O(Δt) 量级的分裂误差因此整体时间精度只有一阶达不到 Lax-Wendroff 的二阶。稳、准、快三者在这里第一次出现明显的取舍。4.4 奇偶步交换为什么要打乱固定的 x-y 顺序固定「先 y 后 x」的分裂顺序会带来一个隐蔽的问题边界外插、入流赋值都是按方向做的两个方向的处理并不对称。长期跑下来误差会沿对角线方向单向积累波形出现定向的偏斜——看起来像是格式不稳其实是分裂顺序带来的各向异性。解决办法是奇偶步交替奇数步用 y → x偶数步用 x → y。这样两个方向的处理在时间上被平均掉一阶的分裂误差不会持续往同一个方向堆。算例里明确写了这个安排代码结构就是把两个分支并排放在 if mod(n,2) 的两侧各自带一套完整的边界处理。if mod(n,2) 1 % 奇数步先 y 后 x % ... y 方向推进 y 出流外插 ... % ... x 方向推进 x 出流外插 ... else % 偶数步先 x 后 y % ... x 方向推进 x 出流外插 ... % ... y 方向推进 y 出流外插 ... end这几十行看起来是重复劳动但少了它长时程积分的结果会慢慢歪掉。分裂法能不能用在长时间模拟里很大程度上取决于这个交替有没有做到位。5. 排查五个把结果画飞的坑5.1 图上只剩几根尖刺其余是平板现象γ 取 1 或 2 时surf 图上有几根极高的柱子其余区域贴着零平面看起来像数据缺了一块。原因放大因子对某些波数大于 1这些模态指数增长几步之内就把数值撑到极大。z 轴自动缩放后被最大值主导其他结构全部压平视觉上就成了平板。这不是缺数据是尺度被吃了。解决把 axis 的 z 范围固定住比如 axis([-20 20 -20 20 0 1])让发散的部分冲出画面而不是压扁其余部分。同时在内层循环里打印 max(abs(U2(:)))一旦这个值单调上升就是发散信号比看图快得多。5.2 波传到右上角又弹回来现象波形走到计算域边界后没有干净地离开而是在角点附近堆积出一圈虚假的回波随时间越来越明显。原因出流边界被当成了入流边界处理给了固定值通常是 0等价于在出口立了一堵墙。或者虽然用了外插但外插是在内点循环中途做的用到了已经被覆盖的邻居值。解决出流边界一律用线性外插 2u_i - u_{i-1}且必须放在内点循环全部结束之后一次性覆盖整条边界。分裂法里两个方向的外插要分别在各自那一步之后做不能合并到最后统一处理。5.3 A 格式和 C 格式共用同一套循环上下界现象把 A 格式的公式替换进 C 格式的循环里程序不报错但边界附近出现异常值或者结果从第一步就不对。原因C 格式只用到 i-1 和 j-1内点范围是 2:endA 格式和 Lax-Wendroff 用到 i1 和 j1内点范围必须是 2:end-1最后一圈留给边界处理。共用 2:end 会让 A 格式在 i end 处访问越界元素Matlab 会静默返回 0 或者报错取决于写法。解决三种格式的循环范围分开写中心类格式统一用 2:size(X,1)-1边界另外补。5.4 方台位置偏了或者形状是长方形的现象初值明明设的是正方形平台画出来却偏在角落或者变成宽窄不一的长条。原因meshgrid 的维度约定没对齐。meshgrid(x, y) 生成的是 X(i,j) x(j)、Y(i,j) y(i)第一维对应 y第二维对应 x。如果按「第一维是 x」的直觉去写 U1(13:21, 13:21)两个方向的索引实际上是一致的所以看不出来但只要 x、y 的网格点数不同立刻错位。解决在脚本开头写一句 size(X) 和 X(1,:)、Y(:,1) 打印出来确认维度比事后调图省事。统一约定「第一维 y、第二维 x」之后所有索引都按这个来。5.5 视频播出来抖动或者根本只有一帧现象用 getframe 和 movie 做动画窗口大小每一帧都在变播出来画面跳或者 writeVideo 写出的文件打不开。原因axis 范围随每帧数据自动缩放figure 尺寸也跟着变getframe 截到的帧大小不一致播放时就会抖。视频文件则需要先 open(v) 再逐帧 writeVideo(v, m)最后 close(v)漏掉任何一步文件都是坏的。解决在循环外固定 set(gcf, Position, [...]) 和 axis([...])让每一帧的尺寸完全一致视频对象在循环开始前创建并 open循环结束后 close中间用 getframe 抓帧而不是抓图。6. 取舍精度、稳定域、实现成本三者只能挑两个四种格式放在一起其实就是一张取舍表。迎风 C 格式稳定域最窄γ ≤ 0.5但只有一阶精度、耗散大、实现最简单A 格式二阶精度、实现同样简单但稳定域为零根本不能用Lax-Wendroff 二阶精度稳定域进一步收紧到 0.3536每步代价也翻倍时间分裂法稳定域最宽γ ≤ 1每步可以用更大的步长但时间精度掉回一阶还多了一个方向交替的实现负担。格式稳定条件时间精度每步模板点数主要代价迎风 Cγ ≤ 0.5一阶3数值耗散大波形被抹平A 格式无条件不稳定一阶名义5纯虚特征值落在稳定域外Lax-Wendroffγ ≤ 0.3536二阶5步长小到达同一物理时刻耗步多时间分裂 Cγ ≤ 1一阶3 3分裂误差需奇偶步交替判断一个格式在具体问题上能不能用光靠看表不够得有一个能自动跑的监控量。稳定格式的总变差、L2 模、最大值都应该是有界的波传播问题时总变差甚至会缓慢下降耗散发散的格式这些量会单调指数上升。把下面这几行塞进时间循环比事后看图靠谱得多tv sum(sum(abs(diff(U2,1,2)))) sum(sum(abs(diff(U2,1,1)))); l2 sqrt(sum(U2(:).^2)); fprintf(n%3d max%.4e TV%.4e L2%.4e\n, n, max(U2(:)), tv, l2);max 看有没有炸TV 看波形有没有长出非物理的振荡L2 看整体能量是不是守恒范围内的有界量。三个量一起看稳定和不稳定几乎一眼就能分开不用等到第 90 步画完图才发现白跑。我自己的习惯是每次换格式、换 γ、换初值都先把这三行打印开着跑一遍确认有界之后再关掉做动画。早年省这一步用 γ 1 跑 Lax-Wendroff图上看着像「波传得挺整齐」其实是数值已经长到 10 的十几次方、surf 的自动缩放把尖刺裁掉了——那次之后我就养成了先看数值再看图的毛病改不掉了。希望这份拆解和代码能帮你少走一段同样的弯路把源码里的循环和边界处理对着自己的推导跑一遍比读十遍公式管用。本文还有配套的精品资源点击获取