CFD格式选型指南：上风、中心差分与QUICK的二维稳态对流扩散对比

搞CFD的人迟早会撞上同一个问题：同一个二维稳态对流扩散方程，迎面风、中心差分、高阶格式分别去解，画出来的场图好像描述的是两个物理问题。一开始我以为是自己代码写错了，后来才发现问题出在离散格式对对流项的处理方式上。

这篇文章围绕的就是这么个算例：用上风方案（Upwind）、一阶中心差分、二阶中心差分（工程上我用QUICK实现）三种格式，在矩形计算域上求解二维稳态对流扩散方程，并用Matlab做对照实验。我尽量把离散原理、代码坑、报告编排都讲透，适合正在学数值传热学、计算流体力学的朋友，也适合做课程设计或项目验证时想快速出图的同学。

1. 这个算例真正要回答的问题：为什么同一方程换种离散方式结果会差这么多

1.1 控制方程与数值任务边界

二维稳态对流扩散方程的标准形式是：

[ \frac{\partial (u\phi)}{\partial x}+\frac{\partial (v\phi)}{\partial y}

\Gamma\left(\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}\right)+S ]

这里 (\phi) 可以代表温度、浓度、焓等等被动标量；(u)、(v) 是已知速度场，本算例取常数；(\Gamma) 是扩散系数；(S) 是源项。方程左边的对流项是“运输”，右边扩散项是“抹平”，两个效应互相竞争。

这个方程在CFD里的地位，相当于“最小测试平台”。几乎所有有限体积法的离散内核，都是从它开始验证的。你别看它形式简单，它对数值格式提出的要求其实相当苛刻：对流项是双曲型的，扩散项是抛物型的，两种特性叠加后，格式的稳定性、精度、有界性就成了一组互相拉扯的指标。

我们最常见的困惑是：用中心差分在高Peclet数下解出来的曲线会出现“锯齿状”摆动，用迎风格式又总感觉峰被“磨平”了。这恰恰说明，问题不在求解器，而在离散方案本身对待对流项的数学处理。这篇文章的实验，就是把这三种格式放在完全相同的网格、边界条件和源项下，让差距直接暴露出来。

1.2 网格Peclet数：一票否决的稳定性指标

在动手写代码前，必须先聊清楚一个无量纲数——网格Peclet数。定义如下：

[ Pe_x = \frac{|u|\Delta x}{\Gamma},\quad Pe_y = \frac{|v|\Delta y}{\Gamma} ]

它衡量的是“在一个网格步长内，对流输运相对于扩散输运的强度”，跟工程上常见的雷诺数、全局Peclet数不是一回事。网格Peclet数小，意味着扩散在每个网格里都能充分把信息抹匀，离散格式很难惹出麻烦；网格Peclet数大，意味着对流压过扩散，格式的“方向性”就变得至关重要。

中心差分格式之所以出名地“脆”，就是因为它对网格Peclet数有一个硬门槛：

[ Pe_x \le 2,\quad Pe_y \le 2 ]

超过这个门槛，离散方程的某个邻居系数会变成负数，数值解就可能出现非物理振荡。而迎风格式没有这个门槛，无条件保持系数符号正确，代价是精度降一阶。这套算例真正想让你直观看到的，就是这个门槛前后的剧烈变化。

我建议你在报告里明确写出 (Pe_x)、(Pe_y) 的定义，而不是只写“网格很密”“扩散很小”。很多人把网格Peclet数和全局Peclet数混着用，评审一眼就能看出概念不扎实。

1.3 关于“一阶、二阶中心差分”的命名，先对齐一下

标题里的“上风及一阶、二阶中心差分方案”，不同教材的叫法其实并不完全统一。我在这篇里的处理方式是：把最简单的中心差分叫“基础中心差分”（CDS），把引入节点曲率修正的高阶中心型格式叫“二阶中心差分”（本文用QUICK实现）。三种格式的界面插值公式如下表，公式比名字可靠：

如果你在教材里看到别的叫法，不用纠结，下面给的公式、代码、数值结论才是真正要对齐的东西。做数值实验最忌讳的就是只记名字、不记公式。

2. 三种离散格式的推导与系数矩阵变化

2.1 控制体积视角下的统一框架

Matlab里组装二维离散矩阵，我最推荐从控制体积通量平衡的角度写代码：对每个节点P的控制体，看东、西、南、北四个界面上对流和扩散通量如何流入流出，然后把所有对未知量的贡献累加到系数矩阵里。这个思路的最大好处是，换格式时只需要改“界面值怎么重构”，其他结构完全不动。

以东界面e为例，总通量拆成两块：

[ F_e^{conv} = u\phi_e,\quad F_e^{diff} = -\Gamma\frac{\phi_E-\phi_P}{\Delta x} ]

扩散通量在均匀网格上永远长一个样，所以三种格式只在对流项的界面重构上有区别。这也解释了为什么要用控制体积方式组织代码——你是在“拼装”同一个几何框架下的不同格式，而不是每次推一套新公式。

2.2 UDS：为什么它永不振荡，却总被嫌弃精度低

上风格式的核心思想直白到近乎粗暴：界面上的值，直接取上游节点的值。假设 (u>0)，东界面 (\phi_e=\phi_P)，西界面 (\phi_w=\phi_W)。代入控制体方程后，得到的五对角系数是：

[ \begin{aligned} a_E &= \frac{\Gamma}{\Delta x^2} \ a_W &= \frac{\Gamma}{\Delta x^2}+\frac{u}{\Delta x} \ a_N &= \frac{\Gamma}{\Delta y^2} \ a_S &= \frac{\Gamma}{\Delta y^2}+\frac{v}{\Delta y} \ a_P &= a_E+a_W+a_N+a_S \end{aligned} ]

所有邻居系数都是正数，没有任何“负系数反噬”的可能，所以UDS非常抗造。但代价藏在截断误差里。把 (\phi_e=\phi_P) 在界面位置做泰勒展开：

[ \phi_P \approx \phi(x_e)-\frac{\Delta x}{2}\frac{\partial \phi}{\partial x}+O(\Delta x^2) ]

也就是说，它等效于在真实对流项旁边塞进了一个额外的扩散项，数值扩散系数大约是：

[ \Gamma_{num} \approx \frac{u\Delta x}{2} ]

一旦真实扩散系数 (\Gamma) 比 (\frac{u\Delta x}{2}) 小得多，数值扩散就喧宾夺主，把真实物理抹平了。二维斜风时还有交叉扩散项，实际数值扩散会更严重，这是UDS最让人头疼的地方。

2.3 基础中心差分（CDS）：二阶精度背后的正系数门槛

基础中心差分的界面值取左右节点的算术平均：

[ \phi_e = \frac{\phi_P+\phi_E}{2} ]

代入控制体方程后，系数变为：

[ \begin{aligned} a_E &= \frac{\Gamma}{\Delta x^2}-\frac{u}{2\Delta x} \ a_W &= \frac{\Gamma}{\Delta x^2}+\frac{u}{2\Delta x} \ a_N &= \frac{\Gamma}{\Delta y^2}-\frac{v}{2\Delta y} \ a_S &= \frac{\Gamma}{\Delta y^2}+\frac{v}{2\Delta y} \end{aligned} ]

对比UDS就会发现，CDS在 (a_E) 项里多了一个负向的 (\frac{u}{2\Delta x})。当 (\Gamma) 不够大，也就是网格Peclet数超过2的时候，(a_E<0) 出现了。负系数意味着数值上“下游节点反向影响上游”，破坏了标量输运的有界性和单调性，于是解开始在空间上振荡。

这个门槛正是中心差分精度再高也压不住的东西。它不是你代码写得不够好，而是离散方程本身失去了“正系数保护”。报告里写这部分时，我强烈建议画一个“(a_E) 随 (Pe_x) 变化”的曲线图，比写十句文字都管用。

2.4 二阶中心差分（QUICK实现）：界面重构加入曲率修正

QUICK格式的思路是：界面值不再用直线连接左右节点，而是用上游两个点、下游一个点做二次插值。假设 (u>0)，东界面取：

[ \phi_e = \frac{3}{8}\phi_E+\frac{6}{8}\phi_P-\frac{1}{8}\phi_W ]

注意它本质上是在基础中心差分 ((\phi_P+\phi_E)/2) 上加了一个曲率修正项 (-(\phi_E-2\phi_P+\phi_W)/8)，相当于用二阶导数信息修正线性插值，所以标题里叫它“二阶中心差分”是合理的，工程实现就是QUICK。

但QUICK有个新手很容易忽略的麻烦：它不再只涉及左右两个邻居。把东、西两个界面合并后，忽略y方向，x方向对流项贡献为：

[ \frac{u}{\Delta x}\left[ \frac{3}{8}\phi_E+\frac{3}{8}\phi_P-\frac{7}{8}\phi_W+\frac{1}{8}\phi_{WW} \right] ]

(\phi_{WW}) 冒出来了，二维情况下矩阵从五点格式膨胀成九点格式。此外，最靠近边界的第二层内部节点没有完整的 (\phi_{WW}) 可用，必须降阶处理——我在算例里对该处退化为CDS，否则程序直接越界。这个细节是初学者最容易踩的坑，后面第3章会重点讲。

3. Matlab核心实现：制造解、稀疏组装与迭代收敛

3.1 为什么用“制造解”（MMS）做基准

既然要做误差分析，就得有一个精确解当裁判。工程上最标准的做法是Method of Manufactured Solutions（MMS），中文一般叫“制造解”：你先自己挑一个光滑函数当精确解，比如：

[ \phi_{ex}(x,y) = \sin(2\pi x)\sin(\pi y) ]

然后把它代回原方程，反推出必须添加的源项 (S(x,y))：

[ S(x,y) = u\frac{\partial \phi_{ex}}{\partial x} +v\frac{\partial \phi_{ex}}{\partial y} -\Gamma\left( \frac{\partial^2 \phi_{ex}}{\partial x^2} +\frac{\partial^2 \phi_{ex}}{\partial y^2} \right) ]

边界值同样由精确解直接给定。这样做的好处很明显：无论格式怎么换，误差都能定量计算；收敛阶也能通过网格加密实测出来。做课程设计时，评审最吃这套——误差表一摆，格式优劣一目了然。

我建议初学时就拿这个正弦解练手，因为源项可以手推出来，不用引入符号计算；等代码跑通了，再换一个带陡峭梯度的双曲正切型解析解，去观察虚假扩散对峰值的“砍头”效应，那才是所有报告里最有说服力的一张图。

3.2 系数矩阵组装：一次循环搞定三种格式

说句实话，二维稳态问题的矩阵组装并不复杂，只要控制好“哪一项移到了右边”。我采用 (A\phi=b) 的约定：对流项散度直接累加正系数到A中，扩散项因为从方程左端移项，累加负系数。核心循环如下：

% 参数与网格 Lx = 1; Ly = 1; Nx = 40; Ny = 40; h = Lx / Nx; u = 1; v = 0.5; Gamma = 0.01; nx = Nx + 1; ny = Ny + 1; nnode = nx * ny; alpha = Gamma / h^2; idx = @(i,j) (j-1) * nx + i; scheme = 'QUICK'; % 'UDS' / 'CDS' / 'QUICK' % MMS解析解与源项 phi_ex = @(x,y) sin(2*pi*x) .* sin(pi*y); S_mms = @(x,y) u*2*pi.*cos(2*pi*x).*sin(pi*y) ... + v*pi.*sin(2*pi*x).*cos(pi*y) ... + Gamma * ((2*pi)^2 + pi^2) .* sin(2*pi*x).*sin(pi*y); % 稀疏矩阵预分配 A = sparse(nnode, nnode); b = zeros(nnode, 1); for i = 2:Nx for j = 2:Ny x = (i-1)*h; y = (j-1)*h; P = idx(i,j); E = idx(i+1,j); W = idx(i-1,j); N = idx(i,j+1); S = idx(i,j-1); % 扩散项（标准五点） A(P,P) = A(P,P) + 4*alpha; A(P,E) = A(P,E) - alpha; A(P,W) = A(P,W) - alpha; A(P,N) = A(P,N) - alpha; A(P,S) = A(P,S) - alpha; % 对流项按格式区分 switch scheme case 'UDS' A(P,P) = A(P,P) + u/h + v/h; A(P,W) = A(P,W) - u/h; A(P,S) = A(P,S) - v/h; case 'CDS' A(P,E) = A(P,E) + u/(2*h); A(P,W) = A(P,W) - u/(2*h); A(P,N) = A(P,N) + v/(2*h); A(P,S) = A(P,S) - v/(2*h); case 'QUICK' if i > 2 && i < Nx && j > 2 && j < Ny WW = idx(i-2,j); % e/w面合并后的x方向对流贡献 A(P,P) = A(P,P) + u/h*3/8; A(P,E) = A(P,E) + u/h*3/8; A(P,W) = A(P,W) - u/h*7/8; A(P,WW)= A(P,WW)+ u/h*1/8; % n/s面合并后的y方向对流贡献（v>0） A(P,P) = A(P,P) + v/h*3/8; A(P,N) = A(P,N) + v/h*3/8; A(P,S) = A(P,S) - v/h*7/8; A(P, idx(i,j-2)) = A(P, idx(i,j-2)) + v/h*1/8; else % 近边界退化为CDS A(P,E) = A(P,E) + u/(2*h); A(P,W) = A(P,W) - u/(2*h); A(P,N) = A(P,N) + v/(2*h); A(P,S) = A(P,S) - v/(2*h); end end b(P) = S_mms(x, y); end end % 边界节点直接强制为精确解 for i = 1:nx A(idx(i,1),:) = 0; A(idx(i,1),idx(i,1)) = 1; b(idx(i,1)) = phi_ex((i-1)*h, 0); A(idx(i,ny),:) = 0; A(idx(i,ny),idx(i,ny)) = 1; b(idx(i,ny)) = phi_ex((i-1)*h, Ly); end for j = 1:ny A(idx(1,j),:) = 0; A(idx(1,j),idx(1,j)) = 1; b(idx(1,j)) = phi_ex(0, (j-1)*h); A(idx(nx,j),:) = 0; A(idx(nx,j),idx(nx,j)) = 1; b(idx(nx,j)) = phi_ex(Lx, (j-1)*h); end % 求解（也可以改成SOR迭代） phi = A \ b;

运行这个脚本后，你应该画出 (\phi) 的等值线图、中心剖线图，并计算与精确解的误差。如果一切正常，低Peclet数下三种格式几乎重合，高Peclet数下CDS立刻出现波浪状振荡。这个“临界点出现”本身就是报告的核心结论。

3.3 SOR迭代怎么写，以及它在高Pe时的表现

虽然上面示例用了Matlab的反斜杠直接求解，但如果你想在报告里展示“迭代收敛性”，写一个逐点Gauss-Seidel加SOR是最直观的做法：

w = 1.5; tol = 1e-8; maxIt = 20000; phiIter = zeros(nnode, 1); for it = 1:maxIt phiOld = phiIter; for i = 2:Nx for j = 2:Ny P = idx(i,j); s = b(P) - A(P,:) * phiIter; phiIter(P) = phiIter(P) + w * (s / A(P,P)); end end if norm(A*phiIter - b, inf) < tol, break; end end

有一点值得强调：中心差分在高Peclet数下，矩阵非对角元素出现正值，不再满足M矩阵性质，此时Gauss-Seidel的收敛速度会急剧下降，甚至彻底发散。你会在迭代日志里看到残差像过山车一样忽大忽小。这是“格式稳定性”在代数层面的直接投影，不是求解器抽风。报告中如果能放一组“不同格式的残差历史曲线”，会非常有说服力。

4. 数值实验：四个典型Peclet数下，三种格式怎么分化

4.1 实验矩阵怎么设计

控制变量法：固定网格步长 (h=1/40)，固定速度和域大小，只改扩散系数 (\Gamma)，就能把网格Peclet数压到不同区间。

这套实验矩阵足够支撑报告里“格式适用性”的讨论。我还建议额外做一组“固定物理参数加密网格”的实验，把 (N) 从40加密到80、160，观察同一格式的误差是否按预期阶数下降。

4.2 从结果图里能读到什么

以工况B（(Pe_x=2.5)）为例，最典型的画面是：沿流动方向切开一条中心剖线，CDS的解在锋面位置出现了明显的锯齿状超射和欠射，数值甚至越过了解析解的上下界；UDS的曲线光滑，但锋面明显被拉宽，峰值高度下降；QUICK在同样的剖线上几乎贴着解析解走。

到了工况C（(Pe_x=12.5)），对比就更戏剧化了：CDS的等值线已经出现棋盘格状的虚假波纹，整张图看过去像被某种高频噪声污染了；UDS依然稳定闭合，却把尖锐的峰形几乎削成了平顶；QUICK也开始出现局部振荡，但幅度比CDS小得多。

这说明一个非常本质的问题：对流占主导时，格式的“单调性”和“精度”很难兼得。UDS牺牲精度换单调性，CDS牺牲单调性换精度，QUICK试图在曲率修正上找回平衡，但依然有天花板。这张对比图放报告里，评审基本不用看文字就能抓住结论。

4.3 误差表与实测收敛阶

下面这组数据是我跑代表性参数得到的趋势，不同解析解或边界处理会有数值差异，但规律稳定：UDS大致一阶，CDS稳定二阶，QUICK在边界退化处理下介于二阶和三阶之间。

收敛阶的计算方法很简单：

[ p \approx \log_2\left(\frac{E_{40}}{E_{80}}\right) ]

把两张网格的 (L_2) 误差代进去即可。值得留意的是，QUICK的实测阶数偏低，主要原因就是我在近边界退化为CDS，边界附近的二阶误差拖了整体后腿。如果你想让QUICK名正言顺地跑出三阶，边界必须用单侧三次插值或解析通量来配合，否则别在报告里宣称“严格三阶”，说实话反而更专业。

4.4 一个容易被忽略的坑：残差收敛不等于解正确

我见过不少人看到残差降到1e-10，就认定算完了。但高Peclet数下的CDS可以做到代数残差很低，场里却满是振荡。原因是振荡由负系数结构本身维持，残差未必能反映网格尺度的虚假波纹。

正确的做法是始终留一个基准解，每步迭代后顺手算一下 (L_\infty) 误差，而不是只看残差。MMS在这里的价值立刻体现出来了——你不光知道“方程解没解完”，还知道“离真解到底差多少”。这个习惯对日后接触复杂工程软件非常有用。

5. 报告写作：图表怎么组织，结论怎么说才站得住

5.1 一图一结论，别把云图当装饰品

审报告的人最怕看到连续三张没头没尾的等值线云图。我的经验是：每张图放出来之前，先问自己“读者从这张图能得出什么唯一结论”。如果是“CDS在高Pe下出现锯齿振荡”，那就只放中心剖线的放大图并标出振荡区域，比整场云图更聚焦；如果是“UDS峰值被抹平”，就放峰值附近的三格式对比线，同时画一条解析解虚线。

配色也值得注意。所有格式的等值线必须用同一个色标范围，比如都固定在 ([-1,1])。否则CDS的振荡会因为色标太宽而被“压掉”，UDS的抹平又可能因为色标太窄而被夸大。这个小细节最容易被忽视，但直接影响结论可信度。

5.2 误差表和收敛阶表才是报告的地基

数值实验报告如果只贴图不贴表，颗粒度远远不够。建议至少两张表：

一是“固定网格、扫描Peclet数”的误差表，行是工况，列是三种格式的 (L_2) 误差和最大误差；二是“网格加密”的收敛阶表，列是不同网格密度下的误差和实测阶数。

写表格时把“无法计算”四个字直接写进CDS高分Pe的行里，反而加分。因为这说明你确实关注到了格式失效的边界，而不是把振荡结果也当误差填进去。

5.3 结论部分怎么写得专业而不空

报告最后的结论，不要用“通过实验证明了……”这种套话。更自然的写法是直接陈述适用边界：

“在本算例参数范围内，UDS在全部工况下保持有界，但峰值误差随网格Peclet数增大而显著上升；CDS在 (Pe_x\le2) 时保持二阶精度，超过后出现非物理振荡；QUICK在 (Pe_x=2.5) 时优势明显，在 (Pe_x=12.5) 时也出现局部振荡，且边界降阶使其名义精度未能完全兑现。总体而言，中等Peclet数范围内QUICK是本算例最均衡的格式。”

这种结论每一句都有实验作为依据，没有空话，评审挑不出毛病。

6. 几个容易翻车的地方和扩展思路

6.1 迎风方向与速度符号

UDS代码里，我默认 (u>0,v>0)，所以对流项只修改了西侧和南侧系数。一旦速度反向，迎风节点就变了，(u<0) 时东界面应取 (\phi_E)，代码里对应的系数方向要整体镜像。很多同学只改速度符号不改离散系数，结果算出完全背离物理的场图。建议写一个内部函数处理“界面上游节点索引”，不要让方向逻辑散落在外层循环里。

6.2 QUICK边界退化问题

我已经提到，QUICK在边界附近没有足够节点构造完整插值，直接截断会失去二阶以上的精度。这里再补充一句：退化的位置和方式不同，实测收敛阶会漂移。如果你想复现“接近三阶”的漂亮数据，可以考虑在边界层使用高阶单侧公式，或者在MMS框架下把边界通量也解析化。对于课程设计来说，明确注明“边界处退化为CDS，整体阶数因此受限”就足够了。

6.3 向非均匀网格、非稳态、非线性方向扩展

这套代码的框架稍加修改就能扩展：非均匀网格需要重写扩散通量中的 (\Delta x)；非稳态问题加上时间推进项后，可以做简单的显式或Crank-Nicolson对比；引入非线性源项后，可以体验Newton迭代与格式稳定性的交互。还有一条路线是用TVD限流器替代固定格式，它能同时保住单调性和二阶精度，是工程上更实用的选择。

在扩展过程中，我个人的体会是：做数值实验最容易被骗的，恰恰是“解算出来了”这个错觉。每次跑完都问自己一句——“我拿什么证明这个解接近真解？”如果回答不上来，说明测试基准还没搭好。这篇算例里MMS和Peclet数扫描，就是回答这个问题的最低成本方案。把这个习惯保持下去，后面做再复杂的流动问题，你都不会在格式选型上吃亏。