简介NASA VOFVolume of Fluid方法源自 NASA 开发的二维两相流界面捕捉程序用于追踪气液等互不相溶流体间的自由表面可在不人为定义界面的前提下模拟喷射、波浪、气泡与液滴破碎等复杂流动现象。这份公开资料为“最终版 NASA VOF V13.0”完整源码压缩包仅含1个C文件整体大小约12KB核心模块集中在一个文件中便于逐行阅读、断点调试或移植到自研求解器。尽管文件名带CPP描述指出其原始计算核心采用Fortran编写推测已被封装或转换为C接口以方便集成。代码基于有限体积法和Navier-Stokes方程覆盖体积分数离散、自由表面重构、时间推进及边界条件处理等关键步骤对研究VOF算法、调试两相流求解器或改写并行版本均有直接参考价值。目前已有97人学习下载是一份轻量但典型的算法示例源码适合具备流体力学和数值计算基础的研究者、工程师及CFD初学者快速理解界面捕捉核心思想。1. NASA-VOF 这套老 Free Fortran 代码为什么做界面捕捉绕不开它打开标题里这个压缩包的人多半第一反应是编译报错。文件名写着 NASA-VOF、VOF2D、界面捕捉解压后却是一堆没有扩展名或 .f 结尾的 Fortran 77 源码读行用标号续行看第 6 列甚至还有 H 格式这种现代编译器早就移除的语法。这套东西没有教程、没有 Makefile 现场演示有的只是上世纪流传下来的二维 VOFVolume of Fluid两相流求解器。可它恰恰是很多做溃坝、液滴碰撞、气泡上升模拟的人翻箱底也要找出来的参考实现F 函数输运、界面几何重构、压力修正投影法界面捕捉的骨架在这里都有而且只有几千行读起来比现在动辄几万行的开源框架轻松得多。这篇笔记要解决的问题很直接这个包是什么、怎么编译、参数怎么填、为什么跑出来界面像糊了一样。适合两类人——新手想用老代码快速验证一个算例熟手想拿它当对照实现来检查自己的求解器。先说结论它能跑但别指望开箱即用坑全在编译选项、单位制和输入读取顺序上下面逐条拆。2. VOF 界面捕捉的数学内核F 函数、PLIC 重构与自由面推进2.1 F 函数与界面单元体积分数怎么映射自由面VOF 方法的核心是一个标量场 F叫体积分数或相体积比。每个网格单元里F1 表示该单元完全被主相占据F0 表示完全被次相占据而 0F1 的单元就是界面单元。所谓界面捕捉就是不直接追踪一条连续的自由面曲线而是用每个单元里的 F 值把自由面“藏”在网格里。这套老代码用的正是这个思路它的开源实现通常就叫 vof2d二维直角网格上的自由面演进全靠 F 输运方程[ \frac{\partial F}{\partial t} \nabla \cdot (\mathbf{u} F) 0 ]这里 F 是一个被动标量它跟着流场走但 F 的梯度又反过来决定界面法向进而影响表面张力与压力场所以 F 输运与 Navier-Stokes 方程是耦合推进的。老代码的常见做法是在每个时间步先解速度-压力耦合再用修正后的速度场显式推进 F最后做几何重构。这套流程里我最想强调的一点是F 只是单元内的一个数值自由面本身的几何细节全部靠重构算法补出来因此重构质量决定了整个界面捕捉的成败。用体积分数而不是拉格朗日粒子或 Level Set 距离函数最大优势在于拓扑变化几乎不需要额外处理。液滴合并、液桥断裂、气泡破碎这类强非线性行为在 VOF 框架下只是 F 的数值在相邻单元之间重新分配不需要重划网格也不存在两个界面段交叉的尴尬。这也是 vof2d 代码虽然老但依然是很多人验证自研两相求解器的基准之一。2.2 从体积分数到界面段PLIC 几何重构是这套代码的灵魂F 值本身说明不了界面朝向。两个单元同样都是 F0.5一个可能是水平自由面另一个可能是竖直自由面光看数值无法区分。所以 VOF 家族里最关键的步骤是界面重构。最粗糙的施主-受主Donor-Acceptor法只用 F 值大小粗判界面方向而这套老代码沿用了更接近 PLICPiecewise Linear Interface Calculation分段线性界面重构的思路先算界面法向再在单元内放置一条直线段使该线段与单元围成的面积恰好等于 F 乘以单元体积。法向由 F 的梯度给出[ \mathbf{n} \frac{\nabla F}{|\nabla F|} ]实际离散时用周围单元的 F 值做中心差分示意代码如下它展示的只是核心斜率计算的骨架c 示意二维 PLIC 法向分量差分 NX (F(I1,J) - F(I-1,J)) / (2.0D0 * DX) NY (F(I,J1) - F(I,J-1)) / (2.0D0 * DY) SLEN DSQRT(NX*NX NY*NY) IF (SLEN .GT. 1.0D-12) THEN UNX NX / SLEN UNY NY / SLEN ELSE UNX 0.0D0 UNY 1.0D0 ENDIF逻辑说明这里先对界面单元0F1做梯度中心差分得到界面法向分量 NX 和 NY归一化后得到单位法向 UNX 和 UNY。SLEN 接近零时说明该单元周围 F 场几乎没有梯度此时强制指定一个默认法向避免除以零。参数说明DX 和 DY 是当前单元的网格宽度老代码里这两个值通常在输入数据里统一给定而不是逐单元存储。差分系数 2.0D0 对应二阶中心差分网格不均匀时这套代码不适用因为它默认均匀网格。单位法向随后会交给界面几何模块用来计算该单元内直线段的截距这一步在源码里通常对应一个名为斜率计算或几何流量的子程序。2.3 表面张力与压力修正界面捕捉如何参与动量方程光有界面几何还不够自由面两侧的流体性质不同表面张力必须折算成体积力进入动量方程。这套老代码的常见处理属于连续表面力CSF模型的雏形不要在界面单元上用 Delta 函数硬点力而是把表面张力按体积分数梯度散布到界面附近几个单元里。曲率计算依靠 F 场二阶导数也就是说界面捕捉的数值噪声会直接传导到曲率上进而扰动压力场。很多人说老代码“玄学”多半是卡在这里界面看起来还行压力场却振荡。动量方程里的压力修正用投影法先预测速度、再解压力 Poisson 方程、最后修正速度使其满足连续方程。F 函数在这条链里的位置很微妙——它不直接参与速度修正但它的重构结果决定下一时刻速度场的边界条件与表面张力源项。所以调参数时不要孤立地看某个系数时间和空间上的耦合误差会在几十步后放大成界面的非物理波纹。理解了这层耦合关系后面编译和调参时才知道为什么某些参数组合必炸。3. 编译 vof2d 老代码最小命令、报错处理和运行目录布局3.1 先别急着编译从压缩包到固定格式源码的识别解压前先确认两件事源码是固定格式还是自由格式以及有没有独立的输入文件。固定格式 Fortran 的行首有严格含义第 1 列是注释位第 6 列是续行标记语句标签放在第 1~5 列。用现代编辑器打开时如果默认把制表符当成空白顶格写入编译大概率报错。我的习惯是先跑一遍 file 命令和 grep而不是直接盲编译unrar x nasa-vof.rar cd nasa-vof file *.f *.for 2/dev/null grep -n READ\|FORMAT\|DATA vof2d.f | head -40逻辑说明unrar 解压后file 命令的作用是确认源码是 ASCII 文本而不是某种编码防止后续编辑器与编译器读到乱码。grep 的作用是快速定位输入读取语句READ和格式化输出语句FORMAT以及 DATA 初始化块。参数说明head -40 限制只输出前 40 行目的是先粗略了解输入变量的数量不至于刷屏。看到 READ 语句越多说明输入卡片行序越敏感改参数时越要以源码里的实际读取顺序为准而不是靠文件名猜测。这一步还有一个隐藏作用确认有没有用 Hollerith 格式形如 14H 开头的字符串常量有的话现代编译器必须加 legacy 选项否则会直接报错。3.2 最小编译命令与三个高频报错老代码最常见的编译武器是 gfortran关键在于标准选项。用默认的 GNU Fortran 标准去编 Fortran 77 老代码基本是在给自己找不痛快。推荐的最小命令我一般写成这样mkdir -p ~/vof_case/run cd ~/vof_case/run cp ../nasa-vof/vof2d.f . gfortran -stdlegacy -ffixed-line-length-none -O2 -o vof2d vof2d.f逻辑说明-stdlegacy 让编译器接受 Fortran 77 时代的过时语法包括 Hollerith 常量、老式 DO 循环和带标号的 DO 终止语句。ffixed-line-length-none 取消固定格式下的行长度限制老代码里经常有超过 72 列的长行不加这个参数会在你最意料不到的地方报“行太长”。O2 是优化级别VOF 这类显式时间步进的求解器开 O2 是安全且明显的性能提升。参数说明vof2d.f 是源文件名如果你手上是 split 成多个子程序文件的版本就直接列出所有 .f 文件或者用 *.f 通配。输出名 vof2d 可自定但运行时要和后面输入文件名对上。如果遇到“Label at … is not in the same block”这类报错常见原因是源码里用了老式 DO 嵌套加 -stdlegacy 就能解决。另一个高频报错是整数溢出或者数值常量越界症状是编译在某个 DATA 行上突然中止提示“Integer too big for its kind”。这不是逻辑错误而是老代码里写了大整数常量当作位标记加一个 -fno-range-check 就能绕过。第三个报错是隐式类型警告老代码通常不开 IMPLICIT NONE变量名首字母 i~n 会被当成整型如果源码里混入了以 i 开头的实数变量运行时会出诡异的截断错误但编译期不会提示。遇到这种问题只能靠读源码上下文判断类型这也是为什么我不建议新手把老代码直接当黑匣子跑。3.3 把输入输出拆出源码目录留着干净的比较基线编译通过只是第一步。老代码运行时的读写路径很别扭它通常硬编码一个输入文件名比如 INPUT 或 fort.5输出也直接写成 fort.7 或 OUT 之类的名称。建议不要在原目录里跑因为一旦输出文件覆盖了源码分布后面改了参数找不到上一次的结果后悔药都没有。我一般把输入数据和输出目录全部剥离mkdir -p ~/vof_case/input ~/vof_case/output cd ~/vof_case/run ln -sf ../input/PROB.IN INPUT ./vof2d ../output/run.log 21 tail -f ../output/run.log逻辑说明ln -sf 建立软链接把输入文件映射成老代码期望的固定文件名这样不用改源码就能指定算例。后台执行并把标准输出重定向到 run.log用 tail -f 持续观察是否有 NaN 或压力发散提示。参数说明PROB.IN 是自定义的输入文件名字可以随算例改INPUT 是老代码默认读取的文件名具体以你手上源码里 OPEN 语句为准有些版本读 fort.10 或直接 READ(5,*)只能靠源码确认。重定向符号 21 把错误流也并入日志处理崩溃时这条日志是唯一的现场记录。4. 输入文件与参数卡片读懂网格、时间步和物性跑通溃坝基准算例4.1 数据卡片的行序与 READ 顺序眼睛和手要对齐这套代码的输入文件没有统一的标准化格式不同流传版本的行序差异很大。最可靠的办法不是找文档而是直接看源码里的 READ 语句。常见结构是前三行读标题与网格控制接下来读边界条件、流体物性、时间步控制最后是输出频率。下面是我按常见顺序整理出的带注释版本它适合用来建立直觉但不保证和你手上的版本逐行一致以源码为准# 用 READ 顺序做卡片模板不是通用标准 Dam-break benchmark 120 60 1.0 0.5 # nx ny length height 0.001 0.01 5000 # dt 输出的时间间隔 总步数 2 2 1 1 # 左右边界类型 上下边界类型 1.0 1000.0 # 密度比 次相密度 0.001 0.5 # 运动粘度 表面张力系数 1 # 初始液面高度占计算域比例逻辑说明第一行是算例名字仅作输出标志第二行给出网格数与计算域尺寸决定 DX 和 DY第三行是时间步长与输出频率边界类型那行通常用整数编号对应自由滑移、无滑移、出流等边界再往下是两相密度和粘度以及表面张力系数。参数说明DT 取多大直接受 CFL 条件约束显式格式下一般要求 u*DT/DX 小于 0.3 左右。输出时间间隔如果太大界面演化过程会丢细节总步数则要按实际物理时间估算。边界类型编号在不同版本间含义不一一定要对照源码里的边界设置模块。4.2 六个必调参数和量纲坑表格比文字直观把反复折磨人的参数集中成一张表照着调能省很多时间。参数作用常见参考值量纲/单位坑NX / NY界面捕捉的空间分辨率100~300网格宽高比接近 1 才稳DT临时推进步长1e-4~1e-3无量纲时别直接填秒We或表面张力系数控制表面张力与惯性比值0.1~10 按特征尺度换算无量纲数物性直接填必炸密度比两相惯性差异10~1000老代码常读的是比值不是绝对值输出步数结果保存频率每隔 10~50 步太频繁磁盘爆满初始液面高度溃坝算例初始条件0.5~0.9 倍高度与网格数对应别小于 3 个网格这里面的量纲坑最值得展开。vof2d 这类老代码默认用特征尺度和特征速度做了无量纲化直接把水的 998kg/m³、1e-6m²/s 填进去压力 Poisson 方程的迭代必然发散。正确做法是先确定特征长度 L 和特征速度 U把实际物性换成一维的 We 数和 Re 数再填卡。没有说明文档时看源码里压力方程和表面张力项有没有除以特征项一眼就能判断代码是否无量纲化。4.3 批量生成算例的小工具用 Python 改卡片找参数区间做参数扫描时手工改卡片效率太低。我一般写个小脚本读模板替换参数并输出到不同目录顺便把 DT 和网格数的 CFL 校验做掉# 批量生成 vof2d 输入卡片并做 CFL 粗校验 import shutil pmax 300.0 # 期望的最大入口流速 for nx in [120, 180, 240]: dx 1.0 / nx # 按 CFL 条件反推 DT 上限取安全系数 0.2 dt 0.2 * dx / pmax text fcase nx{nx}\n{nx} {nx//2} 1.0 0.5\n{dt:.6f} 0.005 20000\n with open(fcase_{nx}/PROB.IN, w) as fp: fp.write(text)逻辑说明这段脚本把空间分辨率作为一个扫描变量每轮按网格宽度重算 DT保证不同网格算例的 CFL 数一致这样结果对比才有说服力。计算得到的 dt 用格式化字符串写入卡片避免浮点数尾部多余位数影响老代码的固定列读取假设。参数说明nx//2 是粗略的纵横比控制仅适用于高度 0.5 的算例pmax 需要根据上游边界条件估计取值偏大时 DT 会偏小、总步数增加取值偏小则可能超出稳定边界。0.2 的安全系数是为了给老代码的表面张力显式项留余量CFL 接近理论值 0.5 时界面会出现齿状抖动。5. 避坑编译过了却跑出鬼界面——五个常见翻车现场与排查5.1 现象输出文件能记事本打开用绘图工具却全乱码原因老代码的输出写的是 FORTRAN 格式化记录每条记录前面有长度头且行内字段宽度由 FORMAT 语句里的 A 格式和空格控制以字节为单位解析时读到的全是偏移。解决按字符逐行读取并用正则拆字段而不是用二进制的二进制读取函数。读取脚本里要严格按列索引切分比如固定从第 11 列到第 20 列取 F 值。我见过某开发者在一个“模拟项目X”里反复怀疑网格不对最后发现只是读文件时跳过了记录头。5.2 现象压力场迭代到几千步后飘出 NaN界面直接消失原因物性填了有量纲的实际值无量纲代码里的压力 Poisson 方程把量级差异放大到无法收敛。密度比 1000 配合有量纲粘度显式表面张力项的稳定性条件被彻底破坏。解决回到无量纲化把 We 数和 Re 数作为输入。换算公式很简单We ρ U² L / σRe ρ U L / μ。如果代码注释里写的是 DIN 之类的变量名先在源码搜索变量在压力方程中的用法确认它到底是绝对值还是比值。5.3 现象编译报错停在 FORMAT 行的 H 格式上原因Hollerith 常量比如 5HHELLO 表示字符串 HELLO在 Fortran 90 标准中已被移除现代 gfortran 默认标准下直接拒绝。解决编译命令加 -stdlegacy这是最省事的路。如果你不想动编译选项可以把源码里的nH...改成对应的带引号字符串但要注意行宽和列对齐改错会让所有后续格式错位。我的原则是老代码优先改编译环境不要轻易改源码。5.4 现象液滴始终不合并两个界面之间始终隔着一个单元的次相原因F 重构时把表面张力按 CSF 模型散布到 3×3 邻域而界面间隙只有 1~2 个网格受拉普拉斯压力影响中间那层薄流体被表面张力撑住形成“虚拟膜”。解决增加界面单元附近的网格分辨率或者适当降低表面张力系数让合并发生的拉普拉斯压力超过膜阻力。这个翻车现场最像玄学因为界面看起来完全正常物理过程却不对只能靠网格收敛性判断。5.5 现象自由面碎裂成雾状小液滴质量守恒残差在 1%原因F 输运没有用差分格式限制器老代码为了速度用了简单的迎风差分在界面切向速度较大的位置产生数值色散。解决检查时间步是否卡在 CFL 边沿必要时把 DT 再降一半。另外看网格宽高比超过 1:2 的单元会让法向梯度计算产生方向性偏差F 的小液滴状撕裂就是典型症状。把网格改成接近正方形能明显改善。6. 用质量守恒和液面高度验证界面捕捉顺带一个输出转图技巧验证老代码算得对不对我习惯只盯两个硬指标。第一个是全局质量守恒残差把每个时间步的 F 场全度量相加理论上总质量应守恒残差公式为[ M_{err} \frac{ \sum_{i,j} |F^{n1}{i,j} - F^{n}{i,j}| }{ \epsilon \sum_{i,j} F^{n}_{i,j} } ]稳定运行的代码这个值应该随时间波动但不超过 0.5%如果超过 2%要么 F 输运有质量问题要么边界条件把流体推出了计算域。第二个是液面高度时间序列溃坝类算例可以直接对比无量纲高度 h/H 与无量纲时间 t√(g/H) 的实验关联线这是最直观的界面捕捉验证。老代码输出的是文本我可以写个小脚本把 F 场转成网格化视图# 把老代码输出的 F 场文本转成二维数组并做云图 import numpy as np import matplotlib.pyplot as plt rows np.loadtxt(F_000100.txt, skiprows1) nx, ny int(rows[:, 0].max()), int(rows[:, 1].max()) field np.zeros((ny, nx)) field[rows[:, 1].astype(int)-1, rows[:, 0].astype(int)-1] rows[:, 2] plt.contourf(field, levels[0.05, 0.5, 0.95], cmapjet) plt.gca().set_aspect(equal)逻辑说明老代码输出的每行对应一个单元的 (i, j, F) 三元组脚本先把三个字段拆出来再按坐标填回二维网格。contourf 只画三个区间级别因为 VOF 界面捕捉我们关心的是 0.5 等值线的大致位置和破碎形态而非逐格纹理。参数说明skiprows1 跳过输出文件的标题行levels 三个刻度的含义分别对应次相、界面、主相。如果输出文件行数不对多半是前面的记录头没跳过回到 5.1 排查。我用这套路径跑过不止一次溃坝算例说实话第一次看到 F 场云图里那条干净的自由面曲线时还是会意外——几千行老代码简单几个参数居然能把界面捕捉做到这个程度。踩过最狠的坑就是物性单位白干了两天从那以后每次跑之前先花十分钟看源码里的无量纲项不省这一步。希望帮到你。本文还有配套的精品资源点击获取
