企业官网建设流程全解析

简介：本资源是一份面向计算流体力学（CFD）初学者与教学实践者的数值模拟代码包，聚焦Maccormack格式求解二维前台阶绕流问题——该经典基准案例广泛用于验证算法对激波、分离区及涡结构的捕捉能力。压缩包仅含1个核心文件（4KB），即完整可编译的C++源码maccormack.cpp，涵盖网格初始化、预测-校正型时间推进、中心差分空间离散、固壁边界处理及结果输出逻辑，代码结构清晰、注释充分，便于理解Maccormack格式二阶精度实现机制与非线性流动建模流程。已有278人学习下载，适合作为CFD入门课程实验、数值方法课程设计或自主复现经典流动案例的可靠起点，读者可直接编译运行、修改参数观察流场演化，并结合Paraview等工具开展后处理分析。

1. 问题从哪来：为什么一个 90 年代的算例到今天还有人翻出来

“maccormack.zip”这个名字，老搞CFD的一看就知道是什么东西——大概率是某位前辈早年用MacCormack格式算前台阶超声速绕流（forward-facing step）留下的程序包。这个算例在上世纪80年代末到90年代初非常流行，几乎是每个做可压缩流数值模拟的人都绕不过去的“入门关卡”。哪怕放到今天，我依然觉得这个算例是检验一个格式好坏、一套代码是否靠谱的试金石。

先说清楚这个算例解决什么问题。它模拟的是超声速来流（通常取马赫数3.0）冲过一个二维矩形管道中凸起的台阶。台阶高度0.2，计算域长3.0、高1.0，台阶位于x=0.6处。来流从左侧均匀进入，右侧出流。流动经过台阶角点时，会产生膨胀扇；在台阶顶部以及上壁面之间，会形成一道弯曲的弓形激波，激波在上壁面反射，反射激波又与台阶顶部产生的斜激波相互作用，最终在下游形成非常复杂的波系结构，包括λ形激波交汇、滑移线、接触间断等。

这个算例之所以经典，是因为它的流场结构极为丰富，但几何又极其简单。一套代码能不能稳定捕捉激波、能不能抑制驻点附近的振荡、能不能正确处理角点奇异性，跑一遍前台阶就全暴露了。MacCormack格式作为二阶精度的显式时间推进格式，在当年是绝对的主力，用它跑这个算例，既有历史意义，也有现实的教学价值。

我自己后来在给学生讲CFD时，也经常把这个zip解压出来当作教学素材。它不依赖任何商用软件，不涉及复杂的网格生成，只需要一份朴素的Fortran或C代码，就能在笔记本上复现教科书里那张经典的密度云图。这种“手搓求解器”的体验，是黑箱软件永远给不了的。

如果你刚接触CFD，或者正在学有限差分法，这个算例非常适合作为第一个完整项目来复现。它麻雀虽小，五脏俱全：有双曲型方程组、有非线性通量、有激波、有边界处理，跑通它，你对数值格式的理解会提升一个台阶。

2. 控制方程与格式原理：MacCormack为什么能处理可压缩流

2.1 二维欧拉方程与无量纲化

前台阶问题用无粘可压缩流模型，也就是二维欧拉方程。守恒形式如下：

[ \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} = 0 ]

其中：

[ \mathbf{U} = \begin{bmatrix} \rho \ \rho u \ \rho v \ E \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \rho u \ \rho u^2 + p \ \rho u v \ u(E + p) \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} \rho v \ \rho u v \ \rho v^2 + p \ v(E + p) \end{bmatrix} ]

这里 (\rho) 是密度，(u)、(v) 是x、y方向速度分量，(p) 是压力，(E) 是单位体积总能。对于理想气体，有：

[ E = \frac{p}{\gamma - 1} + \frac{1}{2}\rho (u^2 + v^2) ]

比热比 (\gamma) 取1.4。整个计算域用均匀来流条件做无量纲化，来流马赫数3.0，密度、速度、压力都归一到量纲为一的量级。这样的好处是数值量级接近1，不容易出现精度丢失，也方便比较不同代码的结果。

2.2 预测-校正：MacCormack格式的半步推进思想

MacCormack格式本质上是Lax-Wendroff格式的一种变体，1969年由MacCormack提出。核心思想是用“预测-校正”两步完成一个时间层的推进：

1. 预测步（Predictor）：先做半个时间步，空间导数用前向差分，得到预测值 (\mathbf{U}^{n+1/2})。
2. 校正步（Corrector）：在预测值基础上，再做半个时间步，空间导数改用后向差分，得到最终修正值 (\mathbf{U}^{n+1})。

以一维标量方程为例，预测步写为：

[ \mathbf{U}_i^{n+1} = \mathbf{U}i^n - \frac{\Delta t}{\Delta x} \left( \mathbf{F}{i+1}^n - \mathbf{F}_i^n \right) ]

校正步写为：

[ \mathbf{U}_i^{n+1} = \frac{1}{2} \left[ \mathbf{U}i^n + \mathbf{U}i^{n+1} - \frac{\Delta t}{\Delta x} \left( \mathbf{F}{i}^{n+1} - \mathbf{F}{i-1}^{n+1} \right) \right] ]

注意校正步中通量 (\mathbf{F}^{n+1}) 是用预测值算出来的。这样做的好处是：单步的截断误差在时间上是二阶精度，空间上也是二阶精度，而且实现起来非常直观，不像隐式格式那样需要解大型线性方程组。

2.3 为什么选MacCormack而不是别的高分辨率格式

现在TVD、WENO、Roe、AUSM这些格式满地走，为什么还要回头看MacCormack？有一个很现实的原因：它在光滑区域表现出色，计算效率极高，代码极其简洁。对于教学、格式对比、快速验证算法想法，MacCormack几乎是零成本的起点。

但它也有致命缺点：在激波附近会产生明显的数值振荡（Gibbs现象）。这是因为二阶中心型格式天然缺少足够的数值耗散。所以实际计算时，必须加入人工粘性项来抑制振荡。这一步在maccormack.zip里应该有体现——没有人工粘性，前台阶算例根本跑不出能看的图，流场会早早点崩溃。

这里补一个关键认知：MacCormack格式不是激波捕捉格式，它是“激波容忍格式”。要让它稳定运行，人工粘性的调参是核心中的核心。

3. 离散实现与边界处理：从网格到时间步进

3.1 计算域与网格生成

前台阶算例的几何很简单，网格一般用均匀笛卡尔网格。计算域取 (0 \le x \le 3)，(0 \le y \le 1)，台阶位于 (0.6 \le x \le 3)、(0 \le y \le 0.2) 的区域。也就是说，台阶占了计算域的底部右侧。网格数量通常用300×100，或者加密到600×200。

网格生成不需要任何复杂工具，直接循环赋值就行：

nx = 300 ny = 100 dx = 3.0d0 / dble(nx) dy = 1.0d0 / dble(ny) do j = 1, ny y(j) = (dble(j) - 0.5d0) * dy end do

注意台阶面内部的网格点要被标记为“固体”，在流场推进时跳过更新。

3.2 初始条件与来流设置

整个流场初始化为均匀来流比自由流参数。无量纲化后的来流条件为：

• (\rho = 1.0)
• (u = 1.0)（对应马赫数3.0，声速约0.333）
• (v = 0.0)
• (p = 1/(\gamma M^2) \approx 0.0794)

这个初场足够简单，时间推进后会自然发展出激波结构。要说明的是，前台阶算例不依赖小扰动启动，本身就是非定常演化问题，直接让流动自己“撞”上台阶就行。

要注意时间步长受CFL条件限制：

[ \Delta t = \text{CFL} \cdot \frac{\min(\Delta x, \Delta y)}{\max(|u| + c, |v| + c)} ]

其中 (c) 是当地声速，CFL数一般取0.5~0.8。用300×100网格时，时间步长大约在 (10^{-4}) 量级，推进到无量纲时间t=4.0，大约需要数万步。我的经验是每100步输出一次密度场，整个计算在普通笔记本上跑几分钟就能完成。

3.3 边界条件的正确打开方式

前台阶算例的边界条件看似简单，实际细节很多。我展开说一下：

左边界（入流）：固定为来流值，每个时间步直接赋值。

do j = 1, ny rho(1, j) = rho_inf u(1, j) = u_inf v(1, j) = 0.0d0 p(1, j) = p_inf end do

右边界（出流）：超声速出流时，所有物理量用外插，最简单是用零阶外插，即边界值等于相邻内点值。对超声速问题这完全足够，因为扰动不会逆流传播。

上壁面和台阶壁面：壁面采用反射边界条件。具体做法是：在壁面外设一层“虚拟网格”，令虚拟网格的法向速度取相反数，切向速度、密度、压力与内点相同。反射边界条件能天然满足无穿透条件，实现起来也非常简单。

这里有个容易踩的坑：壁面并不是整个计算域的下边界，而是分成两段——(0 \le x < 0.6) 的入流段下边界是对称边界，(0.6 \le x \le 3) 的下边界是台阶壁面。前者在y方向上设置为对称条件（相当于虚拟网格的v分量取反号），后者才是真正的固体壁面。如果不区分这两段，流场在台阶上游就会出错。

3.4 人工粘性的引入方式

MacCormack格式在激波附近振荡，需要加入人工粘性来“抹平”振荡。程序中常见的是在更新后的守恒变量上加上一个四阶人工耗散项：

[ \mathbf{U}i^{new} = \mathbf{U}i^{new} + \epsilon \left( \mathbf{U}{i+1} - 4\mathbf{U}i + 6\mathbf{U}{i-1} - 4\mathbf{U}{i-2} + \mathbf{U}_{i-3} \right) ]

这个表达式其实是对压力场二阶导数的离散，系数 (\epsilon) 一般取0.01~0.05。注意人工粘性不能加太大，否则会过度抹平激波，让流场细节失真。这个参数在maccormack.zip里通常已经调好，但如果自己复现，需要仔细测试：

• (\epsilon) 取值太小，激波附近出现明显振荡；
• (\epsilon) 取值太大，激波被抹平得很宽，滑移线细节也看不到了；
• 最佳值通常在0.02左右，且需要同时在x和y两个方向分别施加。

实战建议：调人工粘性系数时，可以每100步输出一帧密度云图，用视频的方式观察振荡情况。只看最终结果很难判断是哪个环节出了问题。

4. 实操复现：解压maccormack.zip后应该怎么跑起来

4.1 解压后的文件结构

拿到maccormack.zip之后，先看一下里面的文件结构。典型的90年代风格代码包会包含这几个文件：

• maccormack.f或maccormack.c：主程序
• init.f：初始化子程序
• boundary.f：边界条件子程序
• output.f：结果输出子程序
• README或readme.txt：简短说明
• run.sh或makefile：编译运行脚本

我当年解压后第一件事是看README，结果里面只写了三行字，大意是“用gfortran编译，运行，得到density.dat”。这对新手其实不够友好，但对有经验的人来说正好——代码本身才是王道。

4.2 编译与运行环境

这个代码年代久远，但用的都是标准Fortran语法，几乎不需要改动就能在现代编译器上编译。推荐用gfortran：

gfortran -O2 -o maccormack maccormack.f init.f boundary.f output.f ./maccormack

如果遇到未定义的变量或函数，大概率是把common块写法弄混了。90年代的Fortran代码经常用common /block/ var1, var2这种老语法，现代gfortran默认仍支持，但如果开启了对标准语法更严格的编译选项，可能会报错。建议不用-std=f2008之类选项，用默认的gnu混合模式编译即可。

运行结束后，会生成类似density_0001.dat这样的文件。这些就是每个输出时刻的密度场数据，可以用Python的matplotlib或Tecplot后处理。

4.3 后处理：把数据变成可以放进论文的图

因为输出的是二维数组，用Python读取和绘图最简单：

import numpy as np import matplotlib.pyplot as plt data = np.loadtxt('density_0001.dat') plt.imshow(data, origin='lower', cmap='jet', extent=[0, 3, 0, 1]) plt.colorbar(label='Density') plt.xlabel('x') plt.ylabel('y') plt.title('Density field - Forward-facing step') plt.savefig('density_step.png', dpi=150)

注意这里要设置origin='lower'，否则图像会上下颠倒，y轴方向反了会让人误以为激波结构位置不对。还有一个常见问题是读取顺序——Fortran按列存储数组，如果代码里输出的是按(j, i)顺序的数组，Python读入后可能需要转置，即data.T。

4.4 典型流场结构与物理对照

时间推进到t=4.0左右，密度云图应该呈现以下经典结构：

1. 台阶角点（x=0.6, y=0.2）处产生膨胀扇，流动在这里加速、密度降低；
2. 台阶顶部上游不远处，一道弯曲的弓形激波从台阶前缘向上游延伸；
3. 弓形激波在上壁面反射，形成反射激波向下游传播；
4. 反射激波与台阶顶部下方延伸的滑移线交汇，形成复杂的λ形结构；
5. 在台阶壁面附近，存在一个低速回流区，这在无粘计算中虽然不如有粘计算明显，但也能从密度分布中看到痕迹。

如果跑出来的图里，激波位置明显偏后或偏前，首先检查边界条件——尤其是台阶角点附近是否被错误地当作壁面处理。角点本身是一个奇点，很多代码在这里会报NaN，常见的解决办法是把它当作普通壁面点处理，不特殊对待，只要人工粘性足够，数值上能稳定推进。

5. 常见问题速查表与调试经验

5.1 高频bug与排查方向

5.2 三大经典陷阱

陷阱一：角点奇点处理

台阶角点（x=0.6, y=0.2）是几何上的奇异点，流动在这里发生拐折，物理上压力是间断的。数值上，角点处如果处理不当，会产生持续的扰动并向全流场扩散。比较稳妥的做法是：角点周围几个网格不做特殊处理，依靠人工粘性压制。不要尝试在角点处做光滑处理，那反而会引入额外的数值误差。

陷阱二：人工粘性方向搞反

人工粘性需要在x和y方向分别施加。如果只在x方向加，y方向上的振荡会积累，表现为竖直方向上的条带状伪影。如果只在y方向加，激波会被拉成水平的宽条。正确做法是各向同性施加，或者至少保证两个方向的系数一致。

陷阱三：边界条件更新顺序

正确的顺序是：先更新内部点，再更新边界点，最后施加人工粘性。如果先把边界点赋好值，再做内部点推进，边界处的数值会被内部计算覆盖掉，导致边界条件失效。这种bug很难直接看出来，因为计算结果“看起来差不多”，但细节会出差错。

5.3 代码里值得留意的细节

翻看maccormack.zip的源码时，我建议大家重点看这几个地方：

• 预测步和校正步的通量计算是否使用同一组状态量；
• 人工粘性是否基于压力梯度做缩放，还是直接加在密度上；
• 台阶内部网格点是否被“跳过”更新，还是被当作壁面点处理；
• 输出频率与总推进步数是否匹配，避免最后一步才输出导致磁盘里只有一套数据。

这些细节决定了代码的稳健性和结果质量。从教学角度看，故意在代码里留几个bug让你去调试，也是老派程序员带新人的常用方法。

6. 格式对比与扩展思考：跑通之后还能做什么

6.1 MacCormack格式与其他主流格式的差异

先用一张表对比MacCormack格式与当前主流格式的核心差异：

MacCormack格式在光滑流动区域有不错的精度，但一到激波附近就需要人工粘性，而人工粘性的大小直接影响结果。相比之下，Roe格式通过Riemann问题的精确解构造通量，天然带有迎风耗散，不需要人工调参。WENO则通过模板自适应选择，在高阶精度和激波捕捉之间取得了很好的平衡。

如果你用MacCormack格式跑通前台阶后，再用Roe或WENO格式跑同一个算例，对比激波厚度和流场细节，会非常直观地感受到格式之间的分辨率差距。这也是我推荐大家做的一个标准练习。

6.2 从二维到三维、从无粘到有粘的扩展路径

跑通前台阶算例后，可以沿着三条路径进行扩展：

1. 加粘性：在欧拉方程基础上加入Navier-Stokes粘性项，并给壁面设置无滑移边界条件。这时台阶后方的回流区会变得更加明显，能够观察到更多真实流动特征。

2. 改几何：把前台阶改成后台阶、楔形、膨胀管，或者加一个对称面模拟管道半截面。几何变了，波系结构完全不同，对格式的考验也各不相同。

3. 换格式：将MacCormack格式替换为Roe或HLLC格式，对比同一算例的激波分辨率、收敛速度、振荡水平。这个对比实验是理解各类格式特性的最佳方式。

6.3 一些进阶建议：如何把这份代码用出更大价值

我个人在实际使用中体会最深的一点是：不要只把maccormack.zip当成一个“跑出云图”的工具，而应该把它当成一个实验平台。我自己做过的几个扩展包括：

• 在代码中加入量纲恢复模块，把无量纲结果还原成真实流动参数；
• 加入滑移线跟踪功能，自动识别流场中的接触间断；
• 把人工粘性系数做成动态可调的，在激波附近自动增大、在光滑区自动减小；
• 把二维代码改写成三维版本，用于简单的超声速进气道初步设计。

每一次修改都会踩新的坑，但正是这些坑让人对数值格式的理解不断深化。MacCormack格式虽然“老”，但它的简洁性决定了它不会被淘汰——它始终是学习差分格式最好的“解剖样本”。

最后再分享一个小技巧：跑前台阶算例时，不要把输出间隔设得太稀。时间推进到t=2.0之前，流场结构还在快速演化，每隔0.1个无量纲时间单位输出一帧，后期可以用这些帧做成动图，能看到激波从台阶前缘逐渐形成、反射、相互作用的完整过程。这个动图对学生理解非定常波系演化的帮助，远比一张最终状态的云图要大。

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