Flow-3D二次开发：bc.f自定义造波边界与编译调试避坑

简介：这份PPT围绕Flow-3D V9.4的二次开发与客制化操作展开，面向需要扩展三维流体模拟功能、修改边界条件或实现特定造波逻辑的工程技术人员与研究者，帮助其理解可编程子程序与执行流程的衔接方式。压缩包内仅含1个PPT文件，大小约1.11MB，内容以图文与代码示例相结合的方式组织，便于集中查阅关键步骤与注意事项。已有1345人学习，说明该主题在相关实践中具有一定关注度。材料从开发环境准备、Visual Studio工程打开、Source目录下Fortran子程序编辑，到编译生成可执行文件、前处理版本选择及运行验证均有覆盖，并以边界条件客制化为例展示波浪高度计算与流场速度更新的实现思路，同时提示备份原文件、注意dum赋值等细节。读者可借此建立Flow-3D客制化流程的整体框架，对照示例排查编译与运行中的常见问题，为后续按项目需求定制物理模型提供参考。

1. 一条对不上的波浪边界，逼出了 flow-3D 二次开发

接手一个港池波浪爬高的算例，卡住的地方往往不是网格，也不是湍流模型，而是入口那条波浪边界。软件自带的几个造波选项跟实测谱型对不上，参数改到怀疑人生，最后发现只能自己写。flow-3D 二次开发里最常被开刀的就是这块：客制化边界条件子程序，用 Fortran 把入射波的水位和水平、垂向速度按你的谱型算出来，再塞回求解器。Flow-3D V9.4 的发行包带了完整工程，路径很明确——Visual Studio 打开 Cumtom_double.sln，改 Source 目录下的 bc.f，编译出一个新求解器，前处理里把版本类型切成 Custom double precision。四步听着简单，每一步都有能耗掉一整天的细节，下面按我实际复现的顺序拆。

2. 挂载点选型：bc.f 在时间步里的位置与 9.1.3 子程序清单

2.1 为什么不改内核，只挂子程序

商业求解器的内核是编译好的二进制，源码不对外。Flow-3D 给出的折中方案是预留一批“空壳”入口：发行包里已经写好函数签名、公共块引用和调用位置，函数体留空或只做默认处理，用户把肉填进去，再整体重新编译链接。手册 9.1.3 节《Subroutines Included in Distribution》列的就是这份清单，它决定了你能改什么、不能改什么。

这种机制的好处很实际。第一，升级换代时改动面小——只要挂载点签名不变，你那套 Fortran 代码能平移到下一个版本，不用重写。第二，调试边界清晰，出问题基本能定位到“我这段代码”还是“我给的参数”，不会陷入改内核改出来的黑匣子。第三，多人协作时冲突可控，两个人分别改不同的子程序，合并时不会互相踩脚。

代价同样明确：你只能改它允许你改的地方。想让求解器换一种压力泊松解法、换一套湍流封闭，这条路走不通。所以选型的第一步不是打开编辑器，而是先翻 9.1.3，确认你的需求落在哪个入口上。

按功能粗分，发行包里的可编程入口大致覆盖这几类：

判断标准就三条：这个量是不是每个时间步都要重算；它是不是只作用在某一类边界上（对应代码里的 wr 标识）；计算它需不需要访问邻居网格（需要就得引入 comdeck 里的索引公共块）。三条都指向边界条件，那就落在 bc.f 上，别去动源项。

2.2 一个时间步里 bc.f 被调用的位置

理解调用时机比理解代码本身更重要。求解器每个时间步大致走这么一条链：推进自由面（更新 VOF 体积分数）→ 施加各类边界 → 组装并求解动量方程 → 解压力泊松方程做速度修正 → 进入下一个时间步。bc.f 挂在“施加边界”这一环，也就是说你写进去的 f、u、w 会作为该边界上的已知值参与离散，而不是算完之后再覆盖。

这带来两个必须记住的后果。其一，你在 bc.f 里给的速度是“边界条件值”，不是“边界附近的解”，两者在数值上会差一层，尤其是边界层网格很粗的时候，别拿它直接跟实测点流速做点对点比较。其二，f 的赋值会直接参与自由面重构，赋错一个量级，轻则水面抖动，重则质量不守恒、计算发散。所以造波这类问题，通常是一边改 bc.f，一边盯着总水体体积的时程曲线，看它有没有缓慢漂移。

还有一个容易忽略的细节：bc.f 在每个时间步会被调用多次，按边界块遍历。如果你写的循环没做范围裁剪，或者忘了用 wr 判断当前处理的是哪一类边界，代码会在不该生效的边界上反复写值，表现为“结果莫名奇妙多了个扰动”。

2.3 拿到示例先读结构，别急着改

发行包里带的那个造波示例，价值不在公式，而在结构。它示范了一件事：怎么在 Fortran 里取到边界处的坐标、怎么引用网格索引、怎么把物理量写回数组。这三件事的写法，所有边界条件客制化都通用。

我一般会先做一次“零改动编译”：一行代码不动，直接编一遍，确认工具链通、产物能生成、跑出来的结果跟自带求解器一致。这一步是后面所有排查的参照系，跳过它，后面任何异常你都分不清是环境问题还是代码问题。所谓血泪经验，就是曾经在一个连环境都没验证的工程上改了两天代码，最后发现是链接了一个旧的目标文件。

3. 从 Cumtom_double.sln 到新求解器：编译链的三段路

3.1 打开解决方案前先认清三件事

第一件，Cumtom_double.sln 这个名字里的拼写就是发行包里的原样，不是笔误，搜索时别自作主张改成 Custom。解决方案里通常挂着一个或多个项目，分别对应不同的构建配置（单精度/双精度、调试/发布），你要确认自己激活的是双精度那一套，因为前处理里最终选的版本类型必须跟它对齐。

第二件，这个工程依赖 Intel Fortran 编译器。Visual Studio 2005 只是外壳，真正把 .f 变成 .obj 的是 ifort。装 VS 的时候如果没勾 Fortran 集成，打开解决方案会看到项目类型不认识，或者编译时提示找不到 fortran 编译器。常见做法是先装 Fortran 编译器、再装 VS，顺序反了要修复安装。

第三件，解决方案文件里记录的是绝对或相对路径。整个源码树被移动过、或者从压缩包里只解出一部分，路径就会断。所以解压之后先别急着打开 .sln，先看目录结构对不对。

# 解压后的相对结构，缺一层都可能导致编译失败 # Custom_double.sln # Source/bc.f <- 用户可编辑的子程序 # comdeck/ijk.f <- 索引公共块 # comdeck/mijk.f # （其余项目文件与库目录） # 动手之前先留后悔药，纯 Windows 环境直接复制粘贴改名也一样 cp Source/bc.f Source/bc.f.bak cp Source/bc.f "Source/bc.f.orig.$(date +%Y%m%d)"

这段命令的重点不是 cp 本身，而是备份的粒度。bc.f 是你唯一要改的文件，也是唯一会在覆盖安装、清理工程、切换分支时被无声抹掉的文件。备份两份的意义在于：一份用于回滚到“改之前”，一份用于对照“我到底改了哪几行”。用 diff 对比 .bak 和当前文件，比凭记忆回想靠谱得多。

3.2 用 Visual Studio 打开并定位到 bc.f

打开 Cumtom_double.sln 之后，在解决方案资源管理器里展开 Source 文件夹，找到 bc.f 双击。如果文件在编辑器里显示为纯文本、没有语法高亮，说明 Fortran 的编辑器扩展没装上，不影响编译，但读代码会累。

定位到 bc.f 之后，第一件事不是改，是找到你要改的那一段。造波示例一般会在文件里留一段被注释掉的参考实现，或者一段带明显标记的代码块。发行包里的示例常被人裁剪过——比如do 8000 k=kprb,kprt这一行被注释掉，只剩do 8000 j=jprf,jprb，循环头尾不配对。这种情况直接编译会报标签错误或者只跑一层网格。遇到就先按 Fortran 的老式标号循环规则补回：

c 循环头必须和结尾的标号语句成对出现 do 8000 k=kprb,kprt do 8000 j=jprf,jprb include '../comdeck/ijk.f' include '../comdeck/mijk.f' c ... 单元内的赋值逻辑 ... 8000 continue

这里的 kprb/kprt 是 k 方向的起止索引，jprf/jprb 是 j 方向的起止索引，都来自公共块。8000 continue是这套老式写法的收尾标号，缺了它，编译器要么报错要么把后面的语句吞进循环，症状非常隐蔽。

3.3 编译、链接与产物确认

这一步最容易被“编译成功”四个字骗过去。解决方案编译成功只说明没有语法错误，不说明你链接进去的是新目标文件，更不说明前处理会用你这个产物。

编译时关注三件事。一是输出窗口里的源文件名，确认编译的确实是 Source/bc.f，而不是某个中间目录里的副本。二是链接阶段的库列表，客制化工程会链接求解器的静态库，库路径失效时报的是 LNK 系列错误，跟代码无关。三是产物时间戳，这是最硬的一条证据。

# 产物必须比源文件新，否则你跑的还是上一轮的二进制 find . -name "*.exe" -newer Source/bc.f -print # 顺手比对一下改动量，确认自己没有多改或漏改 diff Source/bc.f.bak Source/bc.f

第一条命令列出所有比 bc.f 更新的可执行文件。如果结果为空，说明编译产物没落地到你以为的位置，或者编译根本没跑成功只是 VS 报了“已跳过”。第二条命令是收尾检查，配合备份使用：改动量应该跟你预期的一致，多了说明误删，少了说明没保存。

产物就位之后，把新生成的求解器可执行文件放到求解器目录（或者按工程约定替换原文件），这一步做完，才轮到前处理出场。前处理文件里的版本类型必须选 Custom double precision，跟编译时的双精度配置严格对齐。选错了，跑的还是自带求解器，你会看到“代码改了但结果一模一样”的经典现象——这个坑后面单独说。

4. 造波边界逐行拆：dum 参数、fhtz 夹逼与指数衰减项

4.1 相位与两个速度分量：线性波理论怎么落到代码里

先把示例代码完整看一遍。下面是我整理过、补齐了循环结构的版本，注释是我加的：

c 只有垂向阻尼参数有效、且当前是第 6 类边界时才往下走 c 否则直接返回，避免污染其他边界 c if(dum3.lt.ztest .or. wr.ne.6) return c 相位：dum2 是波数，xi(imax) 是边界处的 x 坐标，dum1 是圆频率 phasc=dum2*xi(imax)-dum1*t wvcosc=dum3*cos(phasc) ! 边界处的波面起伏分量 wvsinc=dum3*sin(phasc) ! 与垂向速度同相的分量 c 参考点 im1 处的相位，用来算水平速度 wvcosr=dum3*cos(dum2*x(im1)-dum1*t) i=imax-1 ! 边界内侧一层 do 8000 k=kprb,kprt do 8000 j=jprf,jprb include '../comdeck/ijk.f' ! ijk 索引换算 include '../comdeck/mijk.f' ! imjk / im1 等偏移索引 c 自由面高程 = 静水位 + 波面起伏 fhtw=flhtr+wvcosc c 把高程差折算成单元内的体积分数，并夹到 [0,1] fhtz=max(fhtw-z(k-1),zero)*rdz(k) fhtz=min(fhtz,one) if(fhtz.gt.ztest) then f(ijk)=fhtz ! 部分填充单元 u(imjk)=dum1*wvcosr*exp(dum2*(zk(k)-flhtr)) ! 水平速度 u(ijk) =dum1*wvcosr*exp(dum2*(zk(k)-flhtr)) w(ijk) =dum1*wvsinc*exp(dum2*(z(k)-flhtr)) ! 垂向速度 else f(ijk)=zero ! 干单元清零 u(ijk)=zero w(ijk)=zero endif 8000 continue

对照线性波（Airy 波）理论的标准解：自由面 η = a·cos(kx−ωt)，水平速度 u = aω·cos(kx−ωt)·e^{k(z−d)}，垂向速度 w = aω·sin(kx−ωt)·e^{k(z−d)}。代码里的对应关系一目了然：dum1 是圆频率 ω，dum2 是波数 k，dum3 是振幅 a，flhtr 是静水深 d，exp 项就是随水深衰减的指数因子。相位写作dum2*xi(imax)-dum1*t，符号约定跟标准解一致，说明它取的是沿 +x 方向传播的波。

一个容易被忽略的点：u 同时赋给了 imjk 和 ijk 两个位置，而 w 只赋给 ijk。这是因为水平速度存在交错网格上，边界处的 u 需要同时覆盖边界面和内侧一层，否则动量方程在边界上的通量会出现半格错位。这也是为什么必须引入 mijk.f——imjk 这类偏移索引不是手算的，是从公共块里取的。

4.2 自由面到体积分数：fhtz 的两行夹逼

fhtz=max(fhtw-z(k-1),zero)*rdz(k)这一行是全段最需要理解的。fhtw 是算出来的自由面高程，z(k-1) 是当前网格单元底面的高程，两者相减就是这个单元被水淹没的高度。乘以 rdz(k)（即 1/dz，单元垂向尺寸的倒数），得到的是无量纲的填充比例，正好就是 VOF 里的体积分数 f。

问题在于这个比例可能超出 [0,1]：水位低于单元底面时差值为负，水位高过单元顶面时比值大于 1。所以先 max(..., zero) 把负数砍成 0，再 min(..., one) 把超出的截到 1。两行夹逼之后，f 就落回合法区间。

if(fhtz.gt.ztest)是干湿判断：体积分数大于阈值才认为这个单元是湿的，才赋速度。否则 f、u、w 全清零。ztest 这个阈值直接影响近底层的数值行为——给小了会留下一层几乎没水的“湿单元”，速度场里出现虚假的薄层流动；给大了会把刚被波浪爬升覆盖的单元判成干的，波峰推进时出现台阶。常见做法是取一个跟网格尺度挂钩的小量，改完之后一定要看自由面动画有没有锯齿状跳变。

4.3 comdeck 的两个公共块与索引边界

include '../comdeck/ijk.f'和include '../comdeck/mijk.f'是 Fortran 的老派写法，作用相当于把索引换算的公共块原地展开。ijk.f 把三维下标 (i,j,k) 压成一维数组下标，mijk.f 在此基础上给出 imjk、im1 这类带偏移的索引。

理解这两行有两层价值。一层是性能：手写(i-1) + (j-1)*imax + (k-1)*imax*jmax这种表达式既难读又容易写错，用公共块省事。另一层是正确性：交错网格里 u、v、w、f 的下标基准各不相同，u 在 x 方向半格偏移，w 在 z 方向半格偏移，直接用 ijk 去索引 u 会错位半格。示例里 u 用 imjk 和 ijk、w 只用 ijk，正是这个规则的体现。

还有一条硬约束：include 的路径是相对路径../comdeck/，它相对于当前源文件所在目录。这意味着你不能把 bc.f 单独拷到别处去改，也不能在项目属性里随意改工作目录，否则编译时报“找不到包含文件”。真想挪位置，就把 comdeck 目录一起挪，或者改 include 行。

4.4 dum 参数从哪来：前处理赋值与默认零陷阱

示例里所有物理量都由 dum1、dum2、dum3 驱动，但代码里没有一行给它们赋值。这是设计如此：这三个量由前处理传入。原文那句提醒写得很直白——若这里没有赋予 dum 的数值，就必须在前处理 prepin 中给定。

如果没给会怎样？Fortran 的公共块变量在未显式初始化时，行为依赖编译器和运行环境，通常是 0 或者残留值。dum1=0 意味着圆频率为零，相位不再随时间变化，波面是一条静止的余弦曲线；dum2=0 意味着指数衰减项恒等于 1，速度沿水深不衰减，相当于把一个深水波硬套到浅水；dum3=0 则整个波面振幅为零，什么都不会发生。三种情况表现出来的都是“算出来了，但水面纹丝不动”，很容易误判成代码没生效。

所以每次改完参数，先在前处理里确认三个 dum 值确实写进去了，再跑。给值的物理含义要跟代码对上：dum1 是 2π/T，dum2 是 2π/L（线性波色散关系下由水深和周期反算），dum3 是波高的一半。单位制跟前处理里其他量保持一致，别一个用米一个用毫米。

5. 避坑：五条从编译到发散的踩坑记录

5.1 编译期：路径断掉与产物没更新

现象一：打开解决方案后编译，报找不到../comdeck/ijk.f，但文件明明在那儿。

原因：include 是相对当前源文件目录解析的。压缩包只解出了一部分、或者有人把 Source 目录单独拷出来改，相对层级就断了。另一种情况是项目属性里配了额外的包含路径，把相对路径的解析基准改掉了。

解决：把整个源码树按原样解压，Source 和 comdeck 保持同级。确实要挪，就同步改 include 行，别依赖项目属性兜底——属性是跟着 .vcproj 走的，换台机器就丢。

现象二：改完 bc.f 编译，输出窗口刷了一堆“已是最新”，产物时间戳没变。

原因：VS 的增量编译判定看时间戳，系统时间异常、文件从网络盘拷过来带着旧时间戳，都会让它误判为不需要重编。

解决：先 clean 再 build，然后用find . -name "*.exe" -newer Source/bc.f -print确认产物确实比源文件新。这一步别省，它是唯一能证明“跑的是新代码”的证据。

5.2 运行期：代码生效了但结果没变

现象三：反复确认代码改对了、也编译通过了，跑出来的水面跟改之前一模一样。

原因：前处理文件里的版本类型没切成 Custom double precision，求解时调用的还是发行包自带的求解器。这是最高频的一条，因为整个流程里只有这一个下拉框决定了“用哪个二进制”。

解决：打开前处理，把版本类型改成 Custom double precision，跟编译配置严格对齐。改完存盘，重新生成输入文件再跑。顺手做个对照实验：故意在 bc.f 里把振幅乘 2，如果水面高度没有翻倍，就说明根本没走到你的代码。

现象四：波浪跑起来了，但总水体体积随时间缓慢单调变化。

原因：fhtz 的夹逼没做全，或者干单元清零的逻辑漏了某个速度分量。边界上“进水”和“出水”不平衡，体积就会漂。

解决：先把fhtz=min(fhtz,one)和 else 分支的三清零补全，再检查 u 的双点赋值有没有遗漏 imjk。跑的时候加一条总水体体积的时程输出，斜率不为零就是还有漏。

5.3 数值期：精度混用与阈值玄学

现象五：双精度工程编译成功，结果跟单精度版对比差了一大截，甚至算到中途发散。

原因：编译用双精度、前处理却选了单精度版；或者反过来的组合。两边的实数精度不一致，公共块里的数组长度和内存布局都会错位，读进去的 dum 值是乱的。

解决：编译配置、前处理版本类型、求解器可执行文件三者必须同精度。换机器复现时，把这三项写进算例说明里，别只留一个输入文件。

现象六：自由面动画出现锯齿状台阶，波峰推进时一格格跳。

原因：干湿阈值 ztest 取值不当。网格在垂向不均匀时，固定阈值会让不同层的判定尺度不一致。

解决：把 ztest 跟当地 dz 挂钩做无量纲化，或者干脆做一组敏感性试算，取结果对阈值不敏感的那个区间。这块有点玄学成分，别指望一次调对，凭经验加两次对照算例最稳。

6. 钉死结果：三组对照算例与一条回归命令

改完代码最怕的是“看着像对的”。我的做法是固定三组对照：第一组是零改动基线，用来确认工具链和输入文件没被动过；第二组只改 dum3（振幅），验证波高是否线性响应；第三组只改 dum1（周期），验证波长和相位推进速度。三组都跑通，才说明这段代码是被真正驱动起来的，而不是碰巧输出了个像波浪的东西。

比对不要看云图，看数字。把边界附近某个监测点的水位时程导成 CSV，用脚本算两组之间的最大绝对偏差和相对误差：

import numpy as np # 两列数据：时间, 水位；分别来自基线和客制化算例 base = np.loadtxt("probe_base.csv", delimiter=",", skiprows=1) custom = np.loadtxt("probe_custom.csv", delimiter=",", skiprows=1) t = base[:, 0] eta_b, eta_c = base[:, 1], custom[:, 1] # 截掉前 2 个波周期的启动段，避免初始化的瞬态污染统计 mask = t > 2 * 8.0 d = np.abs(eta_c[mask] - eta_b[mask]) print(f"max|d_eta| = {d.max():.4f} m") print(f"mean|d_eta| = {d.mean():.4f} m") print(f"peak ratio = {eta_c[mask].max() / eta_b[mask].max():.3f}")

skiprows=1跳过表头，mask那行的 8.0 换成你自己的波浪周期，目的是让统计窗口落在稳定段。peak ratio 是最有用的一个数：只改振幅时它应该接近你设定的倍数，偏离超过几个百分点，基本可以断定边界上有单元没被正确赋 f，或者干湿判断把波峰削掉了。

从那以后我每次动 bc.f，都强制走一遍“零改动编译 → 单参数扰动 → 三组对照”这条线，哪怕只改了一个常数。吃过一次亏：某次为了图快跳过了基线对照，结果在一个 include 路径断掉的工程里查了两天，最后发现跑的压根是旧二进制。希望帮到你。

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