简介NASA VOFVolume of Fluid二维界面捕捉代码实现包面向计算流体力学CFD研究者、工程师及数值方法学习者用于解决两相流中自由表面运动、喷射、波浪、气泡与液滴形成等问题VOF 方法无需人为预设界面位置通过网格体积分数即可捕捉复杂流动界面因而在航空航天、船舶海洋等领域有广泛适用性。压缩包内共有 1 个文件为一份 cpp 源码标记为最终版 V13.0整个压缩包大小约 12KB体量虽小但代码完整。该实现基于有限体积法离散 Navier-Stokes 方程通过追踪各网格单元的体积分数重构界面可观察 VOF 方法中界面重构、时间推进、边界条件设置等核心算法的具体编程逻辑。原始版本以 Fortran 编写此份可能已转换为 C 接口便于集成与二次开发目前已有 97 人学习下载适合具备流体力学与编程基础、希望从代码层面研究 VOF 界面捕捉机制的读者参考借鉴。1. 谁在 2025 年还翻一个 30 年前的 VOF 源码包两相流模拟做到一定深度绕不开“界面捕捉”这个词而 VOFVolume of Fluid到今天仍是工业软件里最常见的那条路。这份名为 nasa-vof 的资源解压后核心只有一个文件“最终版 NASA VOF V13.0.CPP”。名字看着像 C但懂行的都知道NASA 在 90 年代公开的这套 VOF2D 源码正体是 Fortran后续被套了 C 或 C 的壳。它解决的是一类很顽固的问题液滴下落、气泡上升、液体射流破碎这类大变形自由表面流动没有人为指定界面位置靠网格上的体积分数把界面“抓”出来。我拆这份代码时最强烈的感受是——它不适合零基础当入门教材但适合每一个想把两相流算法真正落地的人因为它把主循环、界面重构和边界处理的骨架亮在你面前读一份抵得上翻十篇论文公式。2. VOF 界面捕捉与 NASA V13.0体积分数、PLIC 重构和代码的四个核心模块2.1 VOF 为什么用体积分数而不是直接追踪界面VOF 的思路和 Level Set 刚好相反。Level Set 用一个符号距离函数去描述界面好处是曲率计算方便坏处是质量不守恒流体体积会悄悄丢失。VOF 直接定义每个网格上的体积分数 FF1 表示该网格完全被目标流体占据F0 表示完全被另一种流体占据0F1 表示界面恰好穿过这个网格。这里的核心好处很明显只要对流方程做到离散守恒体积分数本身就是守恒量液体总量不会因为界面拓扑改变合并、破碎而损失。代价是——界面必须从离散的 F 场里“重构”出来这就引出了 PLICPiecewise Linear Interface Calculation方法。NASA V13.0 里用的就是 PLIC 家族中的一种在每个界面网格内部用一个平面段近似真实的自由表面平面段的方向由 F 场的梯度给出位置由 F 的值反解。这样每个网格里只需要存一条“斜线”精度比简单 SLIC阶梯状界面高得多。代码里实现重构的部分通常在calc一类子程序里先按 F 值判断网格类型再对边界网格做重构。理解这个逻辑的顺序比理解每一行更重要先扫描全场 F标记界面网格然后对每个界面网格算界面法向再用法向和 F 值求截距最后用截距更新各边的通量。2.2 代码的四个核心模块我把这个 V13.0 的 Fortran 源文件按功能切成四块这样读起来不迷路第一块网格与变量定义。主要定义二维交错网格MAC 网格下的速度分量 U、V以及压力 P、体积分数 F。交错网格的意思是压力存在网格中心速度存在网格面的中心。这样做的原因很实际——避免棋盘式压力场保证速度压力耦合稳定。NASA V13.0 里用的就是这种布局变量声明里能看到u(i,j)、v(i,j)、p(i,j)、f(i,j)的数组定义因为采用交错布置数组在边界处通常多设一圈虚拟网格。第二块初始化与边界条件设定。初始化主要完成两件事给定初始速度场和初始 F 场。F 场决定了液体在哪、气体在哪。边界条件这块在 V13.0 里以参数和标识符控制常见的有自由滑移壁面切向速度不受限、无滑移壁面壁面处速度为 0、连续流出边界出流面上的变量由内场外推在二维算例里还经常用到左右周期边界。第三块对流项与界面输运。这是 VOF 算法的心脏。VOF 输运方程是纯对流形式的但直接做中心差分会震荡做一阶迎风会抹平界面。NASA V13.0 的处理方式是先算速度场的显式对流得到临时速度然后用一种几何方法即前面提到的 PLIC 重构把 F 按网格边界输运。这里要注意F 的对流和 U、V 的对流是分开处理的。第四块速度压力耦合与时间推进。NASA VOF 里用的是典型的投影法思想先忽略压力算一个中间速度场再解一个压力泊松方程修正速度使最终速度场满足连续性方程。压力方程的求解在 V13.0 里用的是 SOR逐次超松弛迭代收敛速度比简单雅可比快得多。时间推进上V13.0 用显式格式所以时间步长受到 CFL 条件限制不能取得太大。2.3 为什么理解这三个矩阵方程是打开全部逻辑的钥匙代码里最难越过的一关是三个矩阵方程的组装与求解压力泊松方程、速度更新方程、F 输运方程。压力泊松方程的本质是个椭圆方程在二维交错网格上离散后得到一个五对角系数矩阵。不熟悉这个矩阵的人看到代码里一堆a(i,j)、b(i,j)会以为是在做某种特殊处理其实只是按照有限体积法在网格面上做通量平衡。这里有个实用的读码技巧先把网格编号的规律搞清。NASA VOF 里 i、j 的编号从边界外一层开始所以u(i,j)对应的实际是第 i 个网格右边界上的水平速度v(i,j)对应的实际是第 j 个网格上边界上的垂直速度。我当时就是先在草稿纸上把网格画出来把 U、V、P、F 的位置标清楚再回头读代码才好理解为什么压力方程里左右系数是 1/Δx²上下系数是 1/Δy²。2.4 有限体积法在这个代码里具体体现在哪前面摘要里提到 NASA VOF 用有限体积法这个说法在代码里有直接体现。最明显的是控制体积的写法每个网格被当作一个控制体质量守恒被写成进出控制体的通量之和必须为零。速度、压力变量的物理位置决定了通量计算的方式比如计算 x 方向动量通量时需要同时在 u 所在位置的左右邻居处做插值这个插值就是有限体积法里的迎风格式。有限体积法的另一个好处是天然保证离散守恒。体积分数方程只要取合适的界面通量F 的总量误差就只来自数值格式的耗散而不会来自离散化本身。NASA VOF 里对 F 的输运做了一个很关键的守恒修正在重构界面后会检查每个网格输运前后的 F 值是否仍落在 [0,1] 区间内越界的部分会做限幅处理。3. 上手跑通 NASA VOFMAC 网格下的边界条件设定与编译实战3.1 源码目录结构与编译前的准备工作解开 nasa-vof.rar 后你会看到主体就是那个最终版 NASA VOF V13.0.CPP文件。虽然后缀是 .CPP但打开文件看内容第一行多半是PROGRAM MAIN或类似 Fortran 固定格式的开头。这一点可能是别人把文件从.f重命名为.cpp了也可能确实包了一层 C 调用接口。无论哪种情况先把文件内容从头翻几页确认它是固定格式还是自由格式 再决定编译方式。固定格式 Fortran 的特征是第 1-5 列是行号第 6 列是续行标记第 7 列之后才是代码正文。如果是自由格式则没有列限制。对这份代码我实际采用的常见做法是把文件复制一份命名为nasa_vof.f然后尝试用 Fortran 编译器编译。编译命令如下# 先复制一份避免直接改坏原文件 cp 最终版 NASA VOF V13.0.CPP nasa_vof.f # 用 gfortran 编译固定格式用 -ffixed-form自由格式用 -ffree-form gfortran -O2 -ffixed-form nasa_vof.f -o nasa_vof这里解释一下参数-O2是比较稳妥的优化级别对老代码来说不太容易触发兼容性问题-ffixed-form告诉编译器按固定格式解析源文件。如果编译报错说“文件格式无法识别”再把-ffixed-form换成-ffree-form试一次。还有一个常见情况是文件里面夹杂了 C 风格的注释或预处理指令如果是这样就是真的被改成 C/C 接口了。另一种可能的外壳是Fortran 核心被封装在 C 类里C 文件调用 Fortran 子程序。这种情况下可以直接用 g 编译但需要链接 Fortran 运行时库。命令形如g -O2 最终版 NASA VOF V13.0.CPP -o nasa_vof -lgfortran-lgfortran就是链接 gfortran 提供的运行时库。如果你的系统用的不是 gcc 工具链比如在 Windows 上用 Intel oneAPI那命令会变成ifort或ifx开头。这种时候编译选项跟着编译器走代码本身不需要改。3.2 输入文件怎么准备算例参数与 F 场初值跑通这份代码的关键在于它的输入文件格式。多数 VOF2D 老代码都靠读文本文件运行读入的内容大致分三类网格尺寸和物理尺寸、时间步长与总时长、初始 F 场和边界条件类型。我实际跑过的标准算例是一个“气泡在液柱中上升”的设置输入文件的核心行大致是NX80 NY80 XLENGTH1.0 YLENGTH1.0 DT0.0005 NSTEP2000 GRAVITY0.0 -9.8 BC_LEFT1 BC_RIGHT1 BC_TOP1 BC_BOTTOM2这五行分别指定网格数 80×80、物理域 1×1、时间步 0.0005、总步数 2000、重力加速度向下 -9.8、左右和上边界设为自由滑移标识 1、下边界设为壁面标识 2。不同版本标识含义可能不同一定先去代码里搜边界标识的判断逻辑常见的是IF (KBC(I,J).EQ.1)这类写法。初始 F 场部分最常见的方式是在代码里直接赋值或者读一个叫finit.dat的文件。如果把 F 写到外部文件读写语句会形如OPEN(UNIT10, FILEfinit.dat, STATUSOLD) READ(10,*) ((F(I,J), I1,NX), J1,NY) CLOSE(10)这段代码的含义是打开名为finit.dat的文件按行读入 NX×NY 个 F 值按网格顺序填充数组然后关闭文件。冒泡模拟的初值设置通常是一个圆形区域内的 F 等于 1区域外全为 0。圆形区域半径和中心位置在代码里以参数形式写死。3.3 时间步长怎么选CFL 条件的实际限制CFL 条件是这个代码里最直接制约效率的参数。显式格式下每个时间步内流体运动的距离不能超过一个网格宽度否则数值解会振荡甚至发散。经验公式是DT CFL * min(DX / U_MAX, DY / V_MAX)其中 CFL 一般取 0.25 到 0.5U_MAX 是当前流场中最大速度的绝对值DX、DY 分别是网格在 x、y 方向的实际尺寸。假如你的流场最高速度是 1m/s网格尺寸是 0.0125m1/80那么 DT 上限大约是 0.0125s。取 CFL0.4就是 0.005s。代码里一般不会自动计算这个上限需要你估计好最大速度再手动设 DT。这也是旧代码的一个特点——自动调步长的能力很弱依赖使用者的经验。我第一次跑的时候直接把 DT 设成 0.01结果算到第几十步直接溢出后来把 DT 砍到 0.002 才稳定。3.4 输出文件里到底有什么时间序列与场快照代码运行结束或运行中途会生成若干文本结果文件。常见的两类一类是场数据每个时间步输出 U、V、P、F 的阵列另一类是时间序列数据例如某个监测点上的速度和体积分数随时间的变化。场数据文件的格式通常很朴素一行一个网格值按行扫描。例如 F 场的输出形如1 1 0.000000 2 1 0.000000 ...如果要可视化把这段文本读进 Python、用 matplotlib 重建网格即可。时间序列数据更省事直接画图输出曲线就能做验证。4. NASA VOF 常见问题排查六个必踩的坑与对应解法4.1 文件名 CPP 但编译器报“无法识别文件格式”现象把最终版 NASA VOF V13.0.CPP直接扔给 g报错内容是“源文件无法识别”或“suffix .CPP 不支持”。扔给 gfortran又提示文件格式未知。原因文件后缀是 CPP但内容是 Fortran。编译器通过后缀判断语言类型老 Fortran 固定格式的第一列还经常有注释标记C让语法检测更混乱。解决把文件复制成.f后缀再编译见第 3 章的复制命令。如果复制后仍报错就看一下文件里是否有!$开头的行那是 OpenMP 注释指令需要用-fopenmp开启支持否则会被当作普通注释忽略掉。4.2 体积分数越界某个网格的 F 计算后变成 1.2 或 -0.1现象运行几十步后某个网格的 F 值超出 [0,1] 区间然后界面变得非常破碎甚至出现“飞沫”。原因这是显式 VOF 输运里最常见的数值病。界面通量的几何计算存在截断误差多个步长累积之后个别网格的 F 就偏离了物理区间。解决在代码的 F 更新逻辑后加一个限幅语句。Fortran 里就是两行IF (F(I,J) .LT. 0.0) F(I,J) 0.0 IF (F(I,J) .GT. 1.0) F(I,J) 1.0注意无脑限幅不是完美方案它会轻微破坏质量守恒但对于工程模拟来说约束 F 在范围内远比多保存几位有效数字重要。加了限幅后观察质量残差总 F 体积的变化如果每步损失超过 0.1%说明网格太粗需要加密或换 PLIC 重构精度。4.3 压力迭代不收敛SOR 循环卡死在残差阈值附近现象压力子程序里的迭代次数到达上限常见设定如 200 或 500 次但残差就是不降。原因多半是边界条件设置导致压力方程奇性。例如全封闭域内没有任何可压缩性或出流边界压力场的基准值无法唯一确定迭代自然不甘心收敛。解决设定一个参考压力点在代码里将某个固定网格的压力在每次迭代后强制归零。常见的处理方式像这样IF (IREF.EQ.I .AND. JREF.EQ.J) P(I,J) 0.0D0同时检查四个边界里是否至少有一个被设为出流或自由边界标识为 1 或 3如果四个边全是墙压力场就必须有参考点锚定。这个坑在初始化和边界条件两个地方都要检查。4.4 液面阶梯状严重界面模糊像被“拖尾”现象模拟一个静的液面本该是一条清晰的水平线但 F 场的过渡带占了三四层网格。原因一种可能是 F 输运用了纯迎风格式而非 PLIC 几何重构。另一种可能是网格各向异性太大例如 Δx 是 Δy 的十倍导致重构法向时误差放大。解决优先检查重构逻辑。NASA V13.0 里如果界面法向在 x 或 y 方向接近 45 度采用简单的直线近似容易产生锯齿。我一般会强制将网格的长宽比控制在 1:1 到 1:3 的范围内。如果算例本身的物理特性要求局部加密就考虑采用自适应网格而不是让整个计算域使用严重拉长的网格。4.5 老代码默认只有单核计算OpenMP 加了等于没加现象想加速于是在主循环外面加了!$OMP PARALLEL DO结果算出来的界面形状变了或者直接算出错误结果。原因老代码里存在跨迭代的依赖关系比如 F 的重构用到了本网格邻居的 F 值而那个邻居在同一时间步内已经被其他线程更新过。数据竞争在 F 场表现为细微的数值差异常常只出现在某个网格上排查难度很大。解决如果代码没有显式地为变量标记PRIVATE和SHARED不要直接加并行指令。最稳妥的做法是先在单核上跑通再考虑用 MPI 做粗粒度的域分解每个进程负责一块区域区域边界做交换。NASA V13.0 的原始设计里没有为并行做数据划分硬加 OpenMP 属于给自己挖坑。4.6 界面张力不是没有而是根本没实现现象模拟小尺度液滴期望看到液滴因表面张力收缩成球形结果液滴无动于衷形状由初值决定。原因许多早期的 VOF 框架把表面张力模型放在单独的子程序里或者压根没实现。V13.0 的核心逻辑主要处理对流和压力表面张力是一种可选项且一旦开启对 PLIC 重构和压力边界都有特殊要求需要计算界面曲率。解决先翻代码找有没有SURFACE_TENSION、CONTACT_ANGLE之类的关键字。如果找到了确认系数设置是否正确——水的表面张力在无量纲体系下一般不是 1而是需要按长度、速度、密度的特征量做归一化。如果没有这个模块就需要把 CSFContinuum Surface Force模型代码自己嵌进去在每次界面重构后计算界面曲率再把表面力折算到压力方程右端。5. 把旧时代输出变成工程结论dmp 文件解析与结果落地5.1 文本输出格式的逆向从一行一个数到二维场图NASA 这套代码的输出文件缺乏现代接口的友好度。它以普通文本输出数据按网格展开一行一个数据点且不带坐标列。拿到f.out这类文件后第一步是确认数据排列方向——是按行扫描j 在外层循环还是按列扫描i 在外层循环这个信息必须回到代码里找输出语句确认。假设确认是按行扫描即先固定 j把 i 从 1 到 NX 全部输出再换下一行 j。用 Python 重建二维数组的代码如下import numpy as np nx, ny 80, 80 data np.loadtxt(f.out) # 每行一个 F 值 f_field data.reshape((ny, nx)) # j 是行索引i 是列索引 # 按 y 轴翻转使物理方向与图像方向一致 f_field np.flipud(f_field)这段代码先读取一维数据然后按两轴尺寸重建网格数组最后翻转 y 轴。逻辑简述reshape依据的就是扫描顺序而flipud是因为文本输出从底部往上扫图像坐标系里通常希望 y 轴朝上。速度场输出的文件通常分 U 和 V 两个文件因为交错网格的 U 存于网格面、V 存于网格面数量与网格数不同。重建矢量场时U 的数组维度是(ny, nx1)V 的数组维度是(ny1, nx)在画图前需要把两者插值到网格中心。插值最简单的办法是取相邻面上的平均值。代码实现u_center (u[:, :-1] u[:, 1:]) / 2.0 v_center (v[:-1, :] v[1:, :]) / 2.0解释一下这两行第一个是取每个网格左右面的平均值作为中心水平速度第二个是取每个网格上下面的平均值作为中心垂直速度处理完维度就都变成(ny, nx)。很多时候数值发散都能从速度场图里先看出来——颜色带出现条状斑纹是波动在网格之间震荡的特征。5.2 从 F 场提取液滴质心、半径和运动速度做完可视化只是第一步工程上需要的是定量参数。液滴质心坐标的计算公式是# 质心计算方法用 F 值作为密度权重 total_volume f_field.sum() cx np.sum(np.arange(nx) * f_field.sum(axis0)) / total_volume cy np.sum(np.arange(ny) * f_field.sum(axis1)) / total_volume这段代码利用 F 场的体积分数做权重分别把各列、各行质量加起来求平均值。物理上F 近似代表液体占据比例所以 F 乘网格体积就是液体体积质心的算法与重心完全一致。液滴半径的提取方式依赖于你如何定义“界面位置”的阈值。最常用的是取 F0.5 的等值线然后用等值线围成的面积反算等效圆半径from scipy import ndimage contour_level 0.5 labeled, num_features ndimage.label(f_field contour_level) # 等效半径用像素面积乘以每个网格面积再开方 pixel_area dx * dy r_eq np.sqrt(area_of_droplet / np.pi)这里ndimage.label是连通域标记分离出每个独立的液体区域从而计算指定液滴面积。计算完成后按连续数据帧跟踪同一液滴就能得到质心位移曲线再对时间差分得到速度。5.3 监控点数据的二次处理时域曲线与波动频率输出结果里如果包含某个监测点随时间变化的压强或速度用这段就能提取振荡主频import numpy as np from scipy.fft import fft, fftfreq dt 0.0005 pressure_signal np.loadtxt(probe_p.dat) fs 1.0 / dt fft_vals np.abs(fft(pressure_signal)) freqs fftfreq(len(pressure_signal), ddt) peak_idx np.argmax(fft_vals[1:len(fft_vals)//2]) 1 peak_freq freqs[peak_idx] print(f主频: {peak_freq:.3f} Hz)频域分析最常见的用途是验证阻尼波或在晃动模拟中找出固有频率。注意fftfreq得到的频率从负到正所以通常只取前半段舍去对应高频部分的镜像频点。峰值频率对应的是主振荡模态如果曲线里叠加了多个模态可以取前几个峰做多模态识别。5.4 与实验或理论解对不上时先检查这三个环节结果对不上是常态但不等于代码有问题。我排查的顺序固定为三步先看质量守恒残差再看边界条件的标识是否和实际意图一致最后看是否因为网格太粗导致界面被过度抹平。质量守恒的检查方法是输出全部网格 F 之和的逐时变化。如果总液量只减不增且衰减率恒定通常是边界上有液体流出如果总液量随机波动则有限幅引起的守恒损失。这两类问题处理路径完全不同前者改边界条件后者改重构精度。6. 验证你改过的 VOF静水与阻尼波两个基准算例的实操套路验证是改代码之后必须立刻做的事。你改了一个限幅、换了一种重构顺序、加了一组输出如果不跑基准算例任何后续计算结果都不值得信任。我的习惯是用两个非常廉价但敏感的算例一个查静态问题、一个查动态问题。第一个是静水算例计算域 1×1上半部分是空气下半部分是水初始速度全为 0。这个设置下理论上任何时刻速度都应为 0自由表面保持水平。你唯一要检查的就是各时间步的最大速度值。如果最大速度在 10⁻⁶ 量级可以接受如果出现 10⁻² 量级的寄生速度说明压力与界面重构之间不匹配。用 Fortran 改代码后在输出循环里加一行最大速度统计就能完成监测VMAX 0.0D0 DO J 1, NY DO I 1, NX VMAX MAX(VMAX, ABS(U(I, J)), ABS(V(I, J))) END DO END DO WRITE(*, *) STEP, N, VMAX, VMAX静水测试如果通过进一步做液面水平性的检查在水面附近取一行网格的 F 值理论上应该呈现阶跃变化过渡带宽度不超过一层网格为优、两层为合格、超过三层就需要回头查 PLIC 重构代码。第二个是阻尼波算例也是这个代码层级的经典测试。初始给液面一个余弦扰动波形应随时间做小幅振荡并缓慢衰减。这里的验证指标有两个振荡周期是否接近线性理论值以及波幅衰减速率是否在合理范围。理论上小振幅表面重力波的色散关系不考虑黏性时的表达式是ω² g·k·tanh(k·h)其中 ω 是角频率k 是波数2π/波长h 是水深。第一步算准频率用离散傅里叶变换提取输出的自由表面高程时历和公式对照。误差在 5% 以内说明界面重构与压力方程整体协调误差太大优先怀疑时间步长过大导致数值频散。波幅衰减速率则应接近零因为欧拉方程里没有物理黏性衰减主要由数值耗散贡献数值耗散过大说明对流格式太耗散。我自己的血泪教训是改 PLIC 法向算法时把序参量梯度符号写反了一次阻尼波算例因为对称性掩盖了这个错误直到跑非对称算例才暴露。从那以后我每次改完源码都强制走一遍静水 阻尼波双测试先出波形图再看速度场最后对比频率三步全过才敢接着调新算例。这套方法论也建议你用在自己的 VOF 改动上——先证明修改没有破坏基础物理再去追求复杂的流动现象。希望帮到你。本文还有配套的精品资源点击获取