简介这份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确认你的需求落在哪个入口上。按功能粗分发行包里的可编程入口大致覆盖这几类入口类别典型用途是否每步都调用边界条件类bc.f 所在类自定义入流波形、随时间变化的壁面条件、特殊出流是源项类在动量或连续方程里加体积力、质量源是物性类变密度、变粘度、非牛顿本构是初始化类初始水位、初始含沙浓度场否仅初始化输出/后处理类自定义监测点、额外写出物理量按输出步判断标准就三条这个量是不是每个时间步都要重算它是不是只作用在某一类边界上对应代码里的 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 kkprb,kprt这一行被注释掉只剩do 8000 jjprf,jprb循环头尾不配对。这种情况直接编译会报标签错误或者只跑一层网格。遇到就先按 Fortran 的老式标号循环规则补回c 循环头必须和结尾的标号语句成对出现 do 8000 kkprb,kprt do 8000 jjprf,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 是圆频率 phascdum2*xi(imax)-dum1*t wvcoscdum3*cos(phasc) ! 边界处的波面起伏分量 wvsincdum3*sin(phasc) ! 与垂向速度同相的分量 c 参考点 im1 处的相位用来算水平速度 wvcosrdum3*cos(dum2*x(im1)-dum1*t) iimax-1 ! 边界内侧一层 do 8000 kkprb,kprt do 8000 jjprf,jprb include ../comdeck/ijk.f ! ijk 索引换算 include ../comdeck/mijk.f ! imjk / im1 等偏移索引 c 自由面高程 静水位 波面起伏 fhtwflhtrwvcosc c 把高程差折算成单元内的体积分数并夹到 [0,1] fhtzmax(fhtw-z(k-1),zero)*rdz(k) fhtzmin(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 是波数 kdum3 是振幅 aflhtr 是静水深 dexp 项就是随水深衰减的指数因子。相位写作dum2*xi(imax)-dum1*t符号约定跟标准解一致说明它取的是沿 x 方向传播的波。一个容易被忽略的点u 同时赋给了 imjk 和 ijk 两个位置而 w 只赋给 ijk。这是因为水平速度存在交错网格上边界处的 u 需要同时覆盖边界面和内侧一层否则动量方程在边界上的通量会出现半格错位。这也是为什么必须引入 mijk.f——imjk 这类偏移索引不是手算的是从公共块里取的。4.2 自由面到体积分数fhtz 的两行夹逼fhtzmax(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 或者残留值。dum10 意味着圆频率为零相位不再随时间变化波面是一条静止的余弦曲线dum20 意味着指数衰减项恒等于 1速度沿水深不衰减相当于把一个深水波硬套到浅水dum30 则整个波面振幅为零什么都不会发生。三种情况表现出来的都是“算出来了但水面纹丝不动”很容易误判成代码没生效。所以每次改完参数先在前处理里确认三个 dum 值确实写进去了再跑。给值的物理含义要跟代码对上dum1 是 2π/Tdum2 是 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 的夹逼没做全或者干单元清零的逻辑漏了某个速度分量。边界上“进水”和“出水”不平衡体积就会漂。解决先把fhtzmin(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,, skiprows1) custom np.loadtxt(probe_custom.csv, delimiter,, skiprows1) 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(fmax|d_eta| {d.max():.4f} m) print(fmean|d_eta| {d.mean():.4f} m) print(fpeak ratio {eta_c[mask].max() / eta_b[mask].max():.3f})skiprows1跳过表头mask那行的 8.0 换成你自己的波浪周期目的是让统计窗口落在稳定段。peak ratio 是最有用的一个数只改振幅时它应该接近你设定的倍数偏离超过几个百分点基本可以断定边界上有单元没被正确赋 f或者干湿判断把波峰削掉了。从那以后我每次动 bc.f都强制走一遍“零改动编译 → 单参数扰动 → 三组对照”这条线哪怕只改了一个常数。吃过一次亏某次为了图快跳过了基线对照结果在一个 include 路径断掉的工程里查了两天最后发现跑的压根是旧二进制。希望帮到你。本文还有配套的精品资源点击获取
