简介本资源是一份面向计算流体力学CFD初学者与教学实践者的数值模拟代码包聚焦Maccormack格式求解二维前台阶扰流问题适用于高校流体力学课程实验、数值方法入门训练及非线性偏微分方程求解器验证场景。压缩包仅含1个核心文件——maccormack.cpp4KB完整实现了Maccormack显式时间推进算法涵盖网格初始化、预测-校正两步差分、固壁边界处理、以及纳维-斯托克斯方程的有限差分离散过程代码结构清晰、注释充分便于理解二阶精度时间积分与空间中心差分的协同实现机制。已有278人学习下载读者可直接编译运行获得前台阶下游典型涡脱落演化特征并导出数据用于Paraview等工具可视化分析是掌握经典显式格式编程实现与流动不稳定性机理分析的实用入门范例。 最近翻硬盘的时候翻出一个压在角落里很久的压缩包maccormack.zip。文件名朴实无华连个版本号都没标。解压一看里面是一套用MacCormack格式计算前台阶扰流问题的Fortran程序附带网格参数、初始场配置和一份极简的README。这个算例就是计算流体力学里大名鼎鼎的Woodward-Colella前台阶流动马赫3的超声速气流冲过一个台阶用来检验格式对激波、膨胀波、接触间断这些复杂结构的捕捉能力。我从这份老代码出发把MacCormack格式的原理、前台阶算例的实现细节、跑完之后的流场分析一次讲透。如果你正在学CFD数值方法或者手头刚好有一个类似的代码包不知道怎么下手这篇文章应该能帮你省不少时间。1. 压缩包里的东西一个经典算例的完整代码包1.1 先解压再说话maccormack.zip里有什么老式CFD代码包的结构通常很有规律。解压maccormack.zip之后你会看到这些文件maccormack.f主程序包含MacCormack格式的预测-校正推进循环mesh.f均匀网格生成子程序init.f初始场设置子程序全流场均匀来流bc.f边界条件子程序包括入流、出流、固壁反射artvis.f人工粘性模块这是格式能压住激波振荡的关键output.f结果输出通常输出密度、压力等量的文本或二进制数据README简单说明运行方式和参数设置解压这一步本身就可能踩坑。如果你在Linux下解压一个Windows传过来的包可能会遇到文件名乱码这通常是因为压缩包用了GBK编码而系统默认按UTF-8解压。用unzip -O GBK可以解决。如果提示file is not a zip file先别急着下结论用file maccormack.zip看一下真实文件类型很多所谓损坏的zip其实是下载不完整或者文件头被传输过程破坏了。分卷压缩的情况更麻烦比如有一个maccormack.z01和一个maccormack.zip必须把分卷文件放在同一目录下再从主zip文件解压缺一个分卷都会报错。如果压缩包确实损坏了Linux下可以用zip -FF尝试修复。说回代码本身。这套程序的语言是Fortran 77风格没有模块化、没有动态数组靠GOTO语句串起整个流程。现代人看这种代码可能觉得头疼但老代码有一个好处逻辑非常直白每个数组、每个循环都摆在明面上非常适合拿来做数值方法教学。1.2 前台阶算例为什么能成为格式试金石前台阶流动问题的出处是Woodward和Colella在1984年发表的一篇经典论文他们提出了几个测试激波捕捉格式的算例前台阶就是其中之一。计算域是一个长为3、高为1的矩形管道从x0.6处开始下壁面抬高了0.2形成一个向前突起的台阶。初场是整个管道里充满均匀的超声速气流无量纲参数取ρ1.4p1.0声速正好为1来流速度取3.0也就是马赫数3。为什么这个算例能火几十年因为它把超声速流动里能遇到的复杂结构几乎全凑齐了。气流撞上前台阶在台阶前缘产生一道强斜激波在台阶上角通道突然扩张又产生一簇膨胀波斜激波打到上壁面发生反射反射激波向下游传播与下壁面再次作用同时从台阶角点还会拖出一条接触间断也就是滑移线这条间断面两侧的气体经历了不同的压缩过程温度密度都不同随着流动推进会逐渐失稳卷起形成涡结构。一个合格的格式必须同时处理好强激波、膨胀波扇、激波反射和接触间断这几类特征。中心差分类的格式通常在激波附近产生振荡耗散太大的格式又会把接触间断抹成一团糊。前台阶算例把这些问题暴露得一清二楚所以它成了CFD领域公认的格式照妖镜。maccormack.zip里这个程序选择这个算例作为测试对象不是偶然的它就是想让你直观地看到MacCormack格式的优缺点。2. MacCormack格式预测-校正两步走的经典显式格式2.1 格式的运作逻辑MacCormack格式是Robert W. MacCormack在1969年提出的属于显式有限差分格式本质上是Lax-Wendroff格式的一种简化实现。它不直接构造复杂的时空耦合差分模板而是通过预测-校正两步在时间上达到二阶精度空间上也达到二阶精度。以守恒型方程组为例一维形式是$$\frac{\partial U}{\partial t} \frac{\partial F}{\partial x} 0$$MacCormack格式的预测步用前差$$U_i^{(p)} U_i^n - \frac{\Delta t}{\Delta x}(F_{i1}^n - F_i^n)$$校正步用后差$$U_i^{n1} \frac{1}{2}\left[U_i^n U_i^{(p)} - \frac{\Delta t}{\Delta x}\left(F_i^{(p)} - F_{i-1}^{(p)}\right)\right]$$两步合起来就等价于一个空间和时间都是二阶精度的显式格式。前差和后差的组合让格式具有了迎风偏向但又不像完全迎风格式那样带着方向性偏差这是MacCormack格式在当时受欢迎的重要原因。二维问题也很容易扩展把x和y两个方向的通量差同时算进去就行。maccormack.zip里的主程序就是二维版本我在第3节会具体拆解。有一个细节值得注意预测步先用前差还是先用后差会影响格式的色散特性。实际使用中很多人会在相邻时间步交替使用两种顺序也就是这一步先用前差下一步先用后差。这样做可以抵消单个方向的偏差让激波附近的振荡更对称。maccormack.zip里的程序有没有做这个交替看代码里的差分调用就能判断如果没有你可以自己加上试试效果差异在激波云图上能看出来。2.2 稳定性边界CFL条件不是数学游戏MacCormack格式是条件稳定的受CFL条件约束。对于一维线性对流方程要求$$CFL \frac{|a|\Delta t}{\Delta x} \leq 1$$二维情况下CFL条件更严格通常要求$$\Delta t \leq \frac{CFL}{\frac{|u|}{\Delta x} \frac{|v|}{\Delta y}}$$其中特征速度要取最大值。前台阶算例里来流马赫3声速为1那么气流方向的最大特征波速度约为uc4也就是4。取CFL0.5网格间距Δx0.0125时间步长就是$$\Delta t \frac{0.5 \times 0.0125}{4} \approx 0.00156$$如果要把流场推进到无量纲时间t4需要约2560步。如果CFL取到0.9步长约0.0028能减少到约1430步但稳定性余量小了很多。我在实际跑这个算例时CFL0.5最稳妥。取到0.8以上时台阶角点附近会慢慢累积微小振荡运行几百步之后可能突然发散。这个现象在对称性很好的网格上尤其明显因为数值误差不会互相抵消反而会在角点这种几何奇异处持续累积。2.3 人工粘性没有它格式就是没装保险的枪MacCormack格式本质上还是中心差分颜色对光滑流场没问题一旦遇到激波数值解会出现明显的过冲和欠冲严重时直接发散。前台阶算例里有强斜激波还有激波在壁面上的多次反射不加人工粘性根本跑不下去。经典做法是加二阶和四阶人工粘性项。二阶粘性在激波附近起平滑作用四阶粘性在光滑区域抑制高频振荡。控制逻辑通常是利用压力梯度作激波指示器压力梯度大的地方加大二阶粘性的权重同时关掉四阶粘性避免在激波附近引入不必要的耗散光滑区域则以四阶粘性为主保持背景的数值稳定性。人工粘性系数的取值经验我会在第5节详细讲。这里先记住一点人工粘性是MacCormack格式能不能用的分水岭调得好激波锐利且稳定调得太大所有流动细节都会被抹平等值线看过去像一块被揉过的面团。3. 从格式到程序前台阶算例的代码拆解3.1 网格与无量纲化的对应关系maccormack.zip里用的是均匀网格计算域x从0到3y从0到1。标准算例的常用分辨率是240×80加密版是480×160。台阶区域是x∈[0.6,3]y∈[0,0.2]这部分格点被标记为固壁不参与流场求解。无量纲化的好处是所有量都是点数值结果可以横向比较。取来流密度ρ1.4压力p1.0这样声速c√(γp/ρ)√(1.4×1/1.4)1马赫3的来流速度就是3.0。所有长度也做了归一化管道长3高1台阶高0.2。这套无量纲参数是Woodward-Colella原算例的标准设置几乎所有论文和代码包都沿用。网格生成本身没有技术含量但有一个容易忽略的点台阶角点(x0.6, y0.2)处于固壁标记的边界上它的处理方式直接影响计算能否稳定进行。老代码通常的做法是在角点附近几个网格内做特殊标记在每一时间步对角点周围施加额外的人工耗散或者对角点的通量做局部光滑处理。3.2 初始场与边界条件超声速流动的边界处理初始场很简单全流场设为均匀来流ρ1.4p1.0u3.0v0。这个初始状态不是定常解因为台阶突然挡在气流路径上从第一步开始就会产生压缩波。边界条件分四块左边界x0超声速入流所有物理量直接给定为来流值不需要任何外插右边界x3超声速出流用一阶外插让信息从内场传出去上边界y1固壁滑移边界法向速度为零下边界y0以及台阶表面y0.2x∈[0.6,3]固壁滑移边界同上下壁面MacCormack格式的边界处理最常用的是镜像网格法。在固壁外侧设一层虚拟格点法向速度取反号密度、压力、切向速度用内点值填充。比如下壁面y0处! 固壁镜像边界 rho(i,0) rho(i,1) p(i,0) p(i,1) u(i,0) u(i,1) v(i,0) -v(i,1)这样在预测步和校正步计算边界处通量时程序不需要特判直接当普通内点处理就行。镜像边界是显式格式最省事的固壁处理方式比通量计算法的逻辑简单得多。台阶内部格点也要处理这些格点不参与推进每个时间步结束后直接把它们的值设回初始状态或者干脆在循环里跳过。老程序更粗暴直接在循环范围里排除台阶区域台阶外的下表边界用镜像网格台阶内的格点冻结在原地。3.3 台阶角点整个程序最脆弱的地方如果你跑过前台阶算例就会发现流场最先出问题的位置几乎永远是台阶角点。MacCormack格式在角点附近会遇到一个尴尬的情况角点上方是流体左侧和下方是固壁几何结构的不连续性让数值通量在角点两侧不对称容易产生持续的数值扰动。我在第一次跑这个算例时没用任何特殊处理结果到数百步时角点附近密度出现负值程序直接崩了。后来加了一手在角点周围3×3的网格范围内把人工粘性系数临时调大2到3倍。这个操作虽然有点土但效果立竿见影。Woodward-Colella原论文里也提到了角点奇异性问题他们的处理更精细一些在角点附近对通量做局部的守恒修正。如果你用MacCormack格式最稳妥的组合是角点周围加强人工粘性同时在角点最邻近的格点上做一次保正修正比如强制密度不低于某个极小值。这样既不会破坏整个流场结构又能保证计算稳定。4. 实测结果分析马赫3下的流动结构4.1 流场演化从均匀流到复杂的激波反射图案程序跑起来之后第一个值得观察的时段是t1左右。这时候斜激波已经从台阶前缘向斜上方延伸打到上壁面后开始反射。台阶上角有一道清晰的膨胀波扇气流在这里加速转向等值线在这个区域会明显变稀这是膨胀波区的典型特征。到t2反射激波已经回到下壁面在下壁面再次反射。此时流场中能看到一个由斜激波和各级反射波构成的菱形反射图案。这个图案是前台阶算例的标志性结构几乎所有格式跑出来的结果都能看到它。继续推到t4接触间断开始卷起。这一条从台阶角点拖出来的滑移线在激波反射的反复作用下失去了稳定性卷成一串涡旋。这时候的密度云图会很好看也最能反映格式的耗散水平。4.2 MacCormack格式的性格缺陷在云图上暴露无遗MacCormack格式跑前台阶算例优点是代码简单、每一步的计算量小缺点也在云图上一眼可见。激波附近有可见的过冲和欠冲表现在密度云图上就是激波等值线附近出现细小的斑纹或锯齿这是中心差分格式的本质属性人工粘性只能压制不能根除。我在240×80网格下观测到激波位置的密度过冲大约有3%到5%如果不加人工粘性这个值会飙升到百分之十几然后直接崩溃。接触间断的耗散是另一个明显问题。MacCormack格式的数值耗散会把接触间断抹成一条弥散的过渡带而不是一条锐利的分界线。对比一些高分辨率TVD格式的结果MacCormack格式跑出的接触间断带明显更宽。如果你把t4的密度云图拉出来能清晰地看到接触间断卷起的涡旋结构发糊这就是耗散的代价。4.3 网格加密能改善什么、不能改善什么把网格从240×80加密到480×160激波和接触间断的锐利程度都会改善计算时间则按大约8倍增长因为网格数翻了4倍时间步数又翻了一倍。加密网格不能根治MacCormack格式的振荡倾向只能让它出现的尺度更小、更不易察觉。我在480×160网格上跑同样的算例激波过冲缩小到了2%左右接触间断的过渡带也明显变窄。这与理论预期是一致的MacCormack格式是二阶精度误差随网格尺度减小而减小。如果你需要的精度更高光靠加密网格会很吃力这时候就得考虑换格式或者至少换成带限制器的TVD类型格式。5. 让格式真正可用的参数调试经验5.1 CFL数的实际取值与精度权衡理论上CFL可以取到接近1我在实际调试中发现MacCormack格式加前台阶算例的组合CFL取0.5是甜点。这个值既保留了足够的稳定余量又不会因为步长太小让计算时间失控。CFL选得太保守比如取0.2除了增加计算时间外还会额外引入时间方向的耗散让激波位置产生微小的偏移。尤其在非定常流动中时间步长直接决定你能分辨的最高频率CFL取得太小等于自己把时间分辨率降下来了。所以我的建议是先按CFL0.5跑通确认流场结构正确后再试着把CFL往上调。如果调到0.7就出现振荡说明问题可能不在CFL本身而在人工粘性设置上。5.2 人工粘性系数怎么调二阶人工粘性系数ε2我通常的取值范围是0.01到0.1四阶系数ε4取0.001到0.01。这两个值是指数量级的概念不是精确值具体取值要根据网格分辨率和来流马赫数微调。调参路径有一个经验可循先用大粘性把计算跑稳比如ε20.05ε40.005确认流场结构合理。然后逐步减小系数每次跑完看激波区域的等值线有没有出现微小振荡。一旦出现振荡说明耗散不够把系数回调20%到30%就行。这个从大往小试的方法比一上来就对着小参数调要快得多。5.3 不同网格规模下的计算成本对比用现代单核CPU跑纯Fortran代码性能已经足够快。我实测了一下网格规模时间步长约推进到t4的步数单步耗时约总耗时约240×800.00156约256030ms约1.5分钟480×1600.00078约5120200ms约17分钟960×3200.00039约102401.5s约4.3小时注意单步耗时不是线性的它包含两层网格的循环同时边界处理、人工粘性、通量计算的开销也随网格规模增长。到了960×320这一档纯串行的Fortran 77代码已经显得吃力。如果你想跑更高分辨率就得考虑OpenMP并行或者干脆换用更高效的格式和更现代的代码结构。5.4 把老代码迁移到现代环境的一些建议maccormack.zip里的Fortran 77代码在现代编译器下通常可以直接编译但会有大量警告主要是数组长度不检查、隐式变量类型这类问题。我建议做几件事把固定大小数组改成可分配数组这样换网格规模时不用改代码重新编译。把GOTO循环改成do-enddo结构可读性提升一个档次。把通量计算、人工粘性、边界条件拆成独立子程序方便单独调试。如果你想用Python复现用NumPy做数组运算可以大幅简化代码。二维的预测和校正步直接对数组切片操作不需要写显式的双层循环代码量能少一半。代价是运行速度比Fortran慢一个数量级不过对于教学目的来说完全够用。我见过有人用Python版本的MacCormack格式配合matplotlib实时输出密度云图把整个流场演化做成了动画这比一堆文本数据直观得多。6. 常见问题排查解压、编译与运行时踩过的坑6.1 解压阶段的几个经典报错虽然标题是maccormack.zip但zip解压的问题在所有压缩包里都会遇到。最典型的几个file is not a zip file这通常不是解压工具的锅而是文件本身的问题。先用file命令看文件头确认它到底是zip还是rar还是压根就是个改名换后缀的文本文件。如果是下载中断导致的损坏重新下载一次远比修复省事。分卷压缩包maccormack.z01加maccormack.zip必须把所有分卷放在同一个目录下目录名和文件名都不能改然后对主zip文件解压。有的解压工具需要手动指定分卷的命名规则遇到报错就检查一下文件名是否连续。解压后文件名乱码在Windows下打包的zip文件名默认是GBK编码Linux下解压会显示乱码。用unzip -O GBK指定编码即可。macOS下打包的zip则相反有时会出现__MACOSX目录和一堆.DS_Store文件这些都可以直接删掉。压缩包损坏但确实需要有价值数据Linux下可以用zip -FF damaged.zip --out fixed.zip尝试修复。注意这是暴力修复结果不一定完整但总比没有强。6.2 编译与运行阶段的数值崩溃排查Fortran老代码编译时最容易踩的是隐式类型坑。比如以i、j、k开头的变量默认是整数其他字母开头默认是浮点数你要是写过index 1.5这种赋值编译器不会报错运行结果就是错的。所有现代工程代码都应该加implicit none老代码没加的话建议自己主动加上。运行时遇到Floating point exception先检查三个地方初始场有没有未初始化的变量、密度有没有变成负数、人工粘性有没有把值推到发散。一个常见操作是在推进循环里加一个密度下界保护if (rho(i,j) 1.0e-6) rho(i,j) 1.0e-6这个操作不改变物理结果但能防止数值误差导致密度变成负值从而引发一连串的NaN。输出NaN时第一反应应该是检查CFL数和人工粘性系数而不是怀疑物理模型。我把CFL从0.5调到0.9测试过一次第800步左右开始出现NaN整个流场瞬间崩溃。老代码的输出文件很大建议每个时间步只输出几个关键量比如全场密度最大值和最小值这样能快速判断问题出在哪个时间段。6.3 结果后处理从二进制数据到密度云图老代码的输出格式五花八门有的直接写ASCII文本几十万行看着都头疼有的用无格式的二进制需要知道记录结构才能读。maccormack.zip里用的是文本格式每行是一个格点的坐标和密度、压力值这种格式最笨但也最好解析。我建议写一个Python脚本读取然后用matplotlib的contourf函数画密度云图。密度云图对激波和接触间断都有很好的分辨力是前台阶算例的首选输出量。云图的色标推荐用coolwarm或者viridis前者能直观看出高低密度区域后者色觉更均匀。如果想把激波位置定量化可以沿某条水平线提取密度剖面观察激波处的密度跳跃位置和宽度。不同格式、不同网格的结果放在同一张图上对比数值耗散的大小一目了然。写在最后跑通maccormack.zip里的程序之后回头再看这套MacCormack格式你会发现它就像一把结构简单的瑞士军刀不算精致但非常实用。它能处理前台阶这类复杂算例能给你一套完整的流场图像代价是你必须理解它的性格——什么时候加人工粘性CFL取多少角点怎么保护。这些经验不是看公式能看出来的必须亲手跑一遍、让程序崩溃几次才能记住。我的体会是老代码包是最好的学习材料但使用方式不是复制粘贴然后直接跑而是先花半天读懂每个数组、每个循环的作用然后故意改坏几个参数用让程序崩掉的方式来理解它的边界。等你亲手把密度云图从一团乱麻调到漂亮的激波菱形结构数值格式在你眼里就不再是公式而是一个有脾气、能交流的东西了。本文还有配套的精品资源点击获取