Fluent船体CFD瞬态仿真：VOF、6DOF与y+参数避坑

简介：面向船舶流体力学学习者与工程人员的 Fluent 船体瞬态仿真案例，聚焦逆海波浪中 Wigley 船体的升沉与纵摇响应分析。案例采用 VOF 模型求解明渠流动，配合明渠波浪边界生成浅层波浪，并以动态网格刻画船体运动；UDF 将运动约束为四自由度，出口设置数值海滩以抑制非物理波浪反射，适合用来理解船舶耐波性仿真的建模思路与关键设置。压缩包共 2 个文件，包含 1 个 gz 算例文件和 1 个 c 格式 UDF 源码，整体约 42.4MB，可直接在 Fluent 中载入复现，或按需调整边界、网格与运动参数。已有 163 人学习，对希望熟悉 VOF 明渠波浪、动态网格及 UDF 二次开发的中级用户，可作为一个较完整的对照样例，便于排查设置细节、验证波浪传播与船体运动结果，并观察四自由度约束与数值海滩对结果的影响。

1. 从一条跑偏的阻力曲线说起：Fluent 船体 CFD 瞬态仿真的适用边界

某开发者拿一条 4 米缩尺船模做 Fluent 船体 CFD 瞬态仿真，定常算出来的总阻力比水池试验高 18%，改成瞬态后曲线先剧烈震荡二十几个周期才慢慢收住，最后落在高出 6% 的位置。他一开始以为是网格不够密，加密到八百万又重算一遍，结果只挪动了 1.5%。问题不在网格，而在他把瞬态当成了“精度更高的定常”来用。船体绕流本身就是强非定常的：艏部兴波、艉部伴流、自由面破碎、船体升沉纵摇，这些东西在定常里被时间平均抹平了，瞬态的价值是把它们还原出来，代价是你必须把时间步、动网格、自由面初始化这些环节一个个摆对。这篇笔记面向已经会用 Fluent 基本操作、准备把船体算例从定常推到瞬态的从业者，也适合刚接手船模阻力或耐波性课题、需要一份可复现参数表的人。

2. 开算之前先把水槽搭对：计算域、边界条件与两相设置

船体瞬态算例翻车，十次里有六次不是求解器的问题，是域和边界一开始就搭歪了。这一章把计算域尺寸、四类边界面的摆法、VOF 自由面的初始化与界面格式、湍流模型选型四件事按落地顺序过一遍，每一步都给可抄的参数区间和背后的理由。

2.1 计算域尺寸：几个船长才够用

船体绕流的计算域不是越大越好，太大了网格量摊薄，自由面附近反而分不到足够的分辨率；太小了出口回流会把艉部压力场污染掉，阻力曲线出现那种“算着算着往上爬”的漂移。我一般按垂线间长 Lpp 做基准，取值区间如下表，缩尺船模和实船尺度都适用，只是实船尺度下要额外注意边界层网格量会爆。

速度入口放在上游面，压力出口放下游面，侧边界和顶部用对称面（symmetry），底部和船体表面用无滑移壁面。这里有个细节值得多说一句：侧边界用对称面而不是壁面，是为了避免侧壁面长出的边界层反过来挤压船体附近的流场；但如果你的算例带横向来流或者斜航，侧边界就得改成速度入口，否则物理上就错了。压力出口的静压参考点建议设在自由面附近、远离船体的位置，出口面不要正对着艉部涡脱落区。

阻塞比要算一下，船体中横剖面浸湿面积除以计算域横截面积，控制在 3% 以内比较稳妥，超过 5% 就得把侧边界往外推。这个数字很多人不算，等到阻力和试验对不上才回头查，那时候网格都划完了，返工成本很高。

2.2 四类边界面怎么摆，VOF 初始化怎么设

两相流设置里，主相选空气，次相选水，两相的物性按试验水温给。水的密度 998.2 kg/m³、黏度 1.003e-3 Pa·s 是 20 摄氏度的常用值，缩尺船模的试验水池通常控温，按实际水温改一下更靠谱。空气按不可压理想气体处理没必要，船体绕流的马赫数极低，直接给常数密度 1.225 kg/m³ 就行，省一层迭代。

自由面的初始化是瞬态算例第一个容易翻车的地方。初始化之后必须用 region adaption 或者 patch 把 z < 0 的区域整体设为水的体积分数 1，z > 0 设为 0，界面附近留一层过渡。这一步漏了，整个域都是空气，船体等于在天上飞，阻力曲线会从第一个时间步就开始发散。

界面重构格式有三个常用选项，选型理由如下：

表面张力在船体阻力算例里通常可以关掉，船模尺度的韦伯数在 10⁴ 量级，表面张力对总阻力的贡献远小于 1%；但如果你的算例要复现艏部翻卷、甲板上浪这类现象，就得开 CSF 模型并把表面张力系数设成 0.072 N/m。开与不开的差别在阻力上体现不出来，在自由面形态上很明显，别拿阻力曲线去判断该不该开。

时间格式先用一阶隐式跑通，等流场稳定、自由面形态合理之后再切二阶隐式重跑。上来就二阶，前期残差震荡会让你分不清是格式问题还是网格问题。

2.3 湍流模型与壁面处理：k-ω SST 还是 Realizable k-ε

船体绕流的湍流模型选型，核心看两件事：艉部逆压梯度区的分离预测，以及壁面处理方式和你的 y+ 是否匹配。

k-ω SST 在逆压梯度下表现稳，艉部分离点和伴流宽度跟试验对得比较好，是船体阻力算例的默认选择。它对壁面距离敏感，低雷诺数处理要求 y+ 在 1 附近，壁面函数处理要求 y+ 落在 30～100，中间那段过渡区（5 < y+ < 30）是它的软肋，尽量别落进去。

Realizable k-ε 更稳健，收敛快，对网格质量的容忍度高，适合前期试算和批量工况筛选；但它对分离区的预测偏保守，艉部伴流容易算厚，阻力偏高个 3%～5% 是常有的事。

如果算例涉及模型尺度下的层流到湍流转捩（艏部前缘一段），可以考虑 γ-Re θ 转捩模型，代价是网格得按低雷诺数要求做，第一层高度要压到 y+ < 1，网格量至少翻一倍。多数工程算例不做这一步，直接用全湍流，然后在结果里注明这个假设。

壁面处理方式在 Near-Wall Treatment 里选，选 Standard Wall Functions 就得保证第一层网格落在对数律区，选 Enhanced Wall Treatment 就得保证 y+ 足够低。两者混用是最常见的错误之一，后面避坑章节会展开。

求解器类型选压力基（Pressure-Based），瞬态，压力-速度耦合用 PISO，PISO 比 SIMPLE 多一次邻接修正，时间步推进时收敛更稳，代价是每步多花约 30% 的时间。PISO 的 skewness correction 和 neighbor correction 各设 1 就够，设大了收敛反而慢。

3. 网格不是越细越好：船体表面网格、边界层与 y+ 的账

网格这一章的核心不是“划多密”，而是“密在哪、第一层多高、质量指标守住哪几条线”。船体瞬态算例的网格量通常在 300 万到 800 万之间，超过一千万而阻力还在漂，八成不是网格量的问题。

3.1 全局网格与局部加密区的划分

船体这类细长体，全局用切割体网格（cutcell）或者 trimmer 比较省事，船体表面贴一层棱柱层，自由面附近做各向异性加密。加密区我一般分四块：

第一块是船体表面，按船长给基础尺寸，Lpp = 4 m 的船模，表面网格尺寸给 0.02～0.04 m，曲率大的地方（艏部、艉部、舵前缘）用 curvature 或者 proximity 再细化一档。

第二块是自由面，水线面上下各 0.3～0.5 倍波高做体加密，网格尺寸取船体表面尺寸的 0.5 倍。自由面加密不够，兴波波形会明显偏平，波高被数值耗散吃掉。

第三块是艉部尾流区，从艉柱往后 2～3 Lpp、横向 0.5 Lpp 的盒子里给一档中等加密，网格尺寸取表面尺寸的 1.5～2 倍。伴流场的涡结构靠这块网格托住。

第四块是艏部兴波区，艏柱前方 1 Lpp、水面上下 0.5 Lpp 的区域加密一档，艏波陡度大，网格粗了直接抹平。

这四块加密叠起来，网格量增长是可控的；反过来，如果全域统一加密，网格量翻倍而阻力只动 1%，性价比很低。

3.2 第一层高度反算：从 y+ 目标推棱柱层参数

这一步是船体网格里最需要算账的地方，凭感觉给第一层高度，十次有八次落在过渡区。反算流程是：先估壁面摩擦速度，再按 y+ 定义反推第一层高度。

以 Lpp = 4.0 m、航速 U = 1.5 m/s、水温 20 摄氏度为例，雷诺数 Re = U·Lpp/ν = 1.5 × 4.0 / 1.003e-6 ≈ 6.0×10⁶。用 ITTC 1957 摩擦系数公式估一下：

Cf = 0.075 / (log10(Re) - 2)^2 = 0.075 / (6.778 - 2)^2 ≈ 0.00329

壁面切应力 τw = 0.5·ρ·U²·Cf ≈ 0.5 × 998.2 × 2.25 × 0.00329 ≈ 3.69 Pa，摩擦速度 uτ = sqrt(τw/ρ) ≈ 0.0608 m/s。第一层高度 y 和 y+ 的关系是 y+ = ρ·uτ·y/μ，反推：

棱柱层增长比给 1.15～1.2，层数 10～15 层，总厚度要盖住边界层厚度。边界层厚度按 δ ≈ 0.37·Lpp·Re^(-1/5) 估，约 0.037 m，棱柱层总厚度做到 0.04～0.05 m 比较稳。增长比设 1.3 以上，相邻层尺寸跳变太大，棱柱层容易在船体曲率大的地方塌陷或者外翻。

这里有个容易被忽略的点：y+ 是沿船体表面分布的，不是常数。艏部前缘流速高、边界层薄，y+ 偏大；艉部边界层厚，y+ 偏小。给第一层高度时按艏部最不利位置算，然后在后处理里画一条 y+ 沿船长的分布曲线，确认 90% 以上的表面落在目标区间内。

3.3 网格质量的三条硬指标

切割体网格的质量检查里，有三条线我每次都看：

最低正交质量（Minimum Orthogonal Quality）不低于 0.15，低于 0.1 的区域基本都会在瞬态推进中出负体积，尤其是船体曲率突变处和棱柱层到四面体的过渡区。

最大长宽比（Aspect Ratio）在边界层里可以放宽到 1000 以上，那是棱柱层的正常形态；但在自由面加密区和尾流区，长宽比超过 200 就得查一下，往往是加密盒的尺寸过渡太急。

最大偏斜度（Skewness）控制在 0.85 以下，超过 0.9 的单元集中在艉部和舵叶后缘，瞬态算到涡脱落阶段最容易在这些位置炸掉。

质量不达标时的处理顺序是：先调加密盒的尺寸过渡比（相邻加密区尺寸比控制在 2 以内），再调棱柱层的增长比和层数，最后才考虑局部重构。直接全域重划，多半是把问题挪了个位置。

网格无关性验证放到最后一章讲，因为那一步的流程和提速技巧绑在一起，单独拎出来说更清楚。

4. 瞬态推进：时间步、6DOF 动网格与造波的落地设置

前面三章的参数定完，进入瞬态推进阶段。这一章的四件事按依赖顺序排：时间步长决定了你后面所有设置的成本，6DOF 决定了船体动不动、怎么动，造波决定了入口边界给什么，监测决定了你什么时候能停。

4.1 时间步长怎么定：库朗数与通过时间双约束

船体瞬态算例的时间步有两个约束，取小的那个。

第一个是船长通过时间约束：一个船长被船体走过的过程至少要分 100～200 个时间步，否则船体运动的历史效应被跳步跳过。Δt ≤ Lpp/(U·N)，N 取 150，Lpp = 4 m、U = 1.5 m/s 时 Δt ≤ 0.0178 s。

第二个是库朗数约束：自由面附近的网格最密，Δx 最小，局部库朗数 Co = U·Δt/Δx 建议控制在 1～5。自由面加密区网格尺寸 0.01 m，U = 1.5 m/s 时，Co = 5 对应 Δt = 0.033 s，Co = 1 对应 Δt = 0.0067 s。

两个约束一夹，Δt 落在 0.005～0.015 s 是比较合理的区间。带强自由面破碎或者 6DOF 大幅运动的算例，往 0.005 s 压；只算静水阻力、船体固定不动的算例，0.01～0.015 s 够用。

全局库朗数在求解器里设 2～5 作为自适应上限，让它自己在局部自动缩步。自适应缩步不是免费的，缩得太频繁会让总计算时间翻倍，所以自定义时间步的上限还是按上面的估算给，自适应只当兜底。

4.2 6DOF 与重叠网格：船体开始动了之后

船体固定不动的阻力算例不需要动网格，自由面用 VOF 就够了。一旦要放开升沉和纵摇（耐波性、波浪中的阻力增加），就得引入 6DOF 求解器，动网格的处理方式有两种：网格重构（smoothing + remeshing）和重叠网格（overset）。

网格重构的优点是设置简单，缺点是船体大幅运动时重构频繁，守恒误差累积，长航时算例会出现总质量漂移；重叠网格的优点是组件网格刚性运动、不重构，守恒性好，缺点是背景网格和组件网格之间的插值面要单独处理，初始化阶段容易出孤儿单元。

船模的 6DOF 属性用 UDF 给，质量、转动惯量、约束方向都在里面定义：

#include "udf.h" /* 船模六自由度属性：只放开垂荡和纵摇，其余锁死 */ DEFINE_SDOF_PROPERTIES(ship_sdof, prop, dt, time, dtime) { /* 排水质量，按缩尺船模实测值填，单位 kg */ prop[SDOF_MASS] = 12.5; /* 绕三个轴的转动惯量，单位 kg·m^2，务必核对量纲 */ prop[SDOF_IXX] = 1.85; prop[SDOF_IYY] = 1.90; prop[SDOF_IZZ] = 0.42; /* 约束：横荡、横摇、艏摇锁死，只留垂荡和纵摇 */ prop[SDOF_ZERO_TRANS_X] = TRUE; prop[SDOF_ZERO_ROT_X] = TRUE; prop[SDOF_ZERO_ROT_Z] = TRUE; /* 重心位置，相对船中，单位 m */ prop[SDOF_LOAD_M_X] = 0.0; prop[SDOF_LOAD_M_Y] = 0.0; prop[SDOF_LOAD_M_Z] = 0.02; }

代码的逻辑很直白：把船模当成一个刚体，告诉求解器它有多重、转起来有多难、哪些方向不许动。三个容易出错的地方分别说。质量单位是 kg，不是排水量吨，填 12.5 对应 12.5 kg 的缩尺船模，填成 12500 就是量纲错。转动惯量的量纲是 kg·m²，缩尺船模的惯量通常在零点几到几之间，如果算出来是几千，八成是把 mm 当成了 m。约束宏在求解器里是取反逻辑，TRUE 表示该方向被锁死，这一点跟直觉相反，第一次写很容易把想放开的方向锁掉，船体一动不动还以为是流场的问题。

重叠网格的设置顺序是：先划背景网格（整个计算域），再划组件网格（包住船体的一层），组件网格的边界设成 overset 类型，然后在 Fluent 里创建 overset interface，把组件网格边界和背景网格关联起来。初始化之后检查孤儿单元（orphan cells）数量，超过总网格量的 0.1% 就得回头调组件网格的包络范围。

4.3 造波与消波：入口边界与数值海滩

静水阻力算例的入口给均匀来流就行，波浪中的算例得在入口造波。Fluent 内置的 Open Channel Wave BC 支持一阶到五阶斯托克斯波和椭圆余弦波，参数表如下：

入口面设成 Open Channel Wave，顶部设压力出口并把 open channel 选项打开，这样气相能自由进出。造波边界下面要留够 1～2 个波长的过渡区，波形在这个距离内从理论解过渡到数值解，太短了波还没长起来就撞到船体上。

消波放在出口段，用数值海滩（numerical beach）或者多孔介质阻尼区。数值海滩的长度取 1.5～2 个波长，阻尼系数从 0 线性增加到最大值，避免波在出口面反射回来污染整个域。不做消波，出口反射波会和入射波叠加，形成驻波，阻力曲线出现周期性的大幅震荡，看起来像物理现象，其实是边界条件造的假。

4.4 监测量与“算完了”的判断

监测项至少设四个：船体总阻力（X 方向力）、升力（Z 方向力）、纵倾力矩（绕 Y 轴）、自由面在某几个监测点的高度。

阻力系数按 Cd = Fx/(0.5·ρ·U²·S) 归一化，S 取静水湿表面积。缩尺船模的 S 大约 1.2 m²，U = 1.5 m/s、Fx = 5.4 N 时，Cd ≈ 0.004，量级对得上。

“算完了”的判断不能只看残差。连续性残差降到 1e-4 以下只是必要条件，真正的判据是：阻力曲线的滑动平均在最近 3～5 个船长通过时间里变化小于 1%，升沉和纵摇的振幅包络稳定，自由面监测点的波高不再漂移。瞬态算例的残差曲线本来就是震荡的，拿定常那套残差标准去卡，永远等不到“收敛”。

统计平均的窗口也要说清楚：阻力曲线前 2～3 个船长通过时间是启动段，不能进平均；取后面 5～8 个通过时间做时间平均，同时记录脉动均方根。只报一个平均值不报脉动，等于把瞬态算成了定常。

5. 避坑与排查：五个把我按在椅子上重算的现场

这一章的五条都是实际算例里撞出来的，按“现象 → 原因 → 解决”写，顺序大致对应算例推进的先后，前面几条偏设置，后面几条偏结果。

5.1 现象：初始化后第一步就发散，残差冲到 1e3

原因基本锁定在自由面初始化。整个域被默认填成水或者填成空气，船体初始湿表面跟实际差了一个船宽，第一步的压力场修正量巨大，求解器直接崩。

解决：初始化之后手动 patch 区域。用 Mark/Adapt 按 z 坐标切出 z < 0 的体区域，patch 成水的体积分数 1；z > 0 patch 成 0；界面上下各 0.5 倍波高的薄层 patch 成 0.5 作为过渡。patch 完之后在初始化面板里检查水的总质量，跟理论排水量对一下，差太多说明切区域的时候把船体内部也算进去了，得先把船体实体区域从 patch 范围里减掉。

5.2 现象：阻力比试验值系统性高 10% 以上，加密网格只降 1%

原因多半是壁面处理方式和 y+ 不匹配。第一层高度给到 0.5 mm，y+ 落在 10 附近，正好在 k-ω SST 的过渡区，壁面函数用不上、低雷诺数解析又不够，壁面切应力被高估，阻力整体抬高。

解决：先在后处理里画 y+ 沿船长的分布，确认实际落点。要么把第一层压到 0.02 mm 走低雷诺数路线（网格量涨 30%），要么把第一层放到 1 mm 以上走壁面函数路线（阻力精度略降但稳定）。两条路选一条，不要停在中间。改完之后阻力通常会往试验值靠 5～8 个百分点。

5.3 现象：自由面算着算着变成一坨模糊的过渡带，兴波波形消失

原因是界面格式和时间步不匹配。用 Geo-Reconstruct 配 0.03 s 的时间步，自由面处库朗数到了 15 以上，界面重构来不及，数值扩散把界面抹厚了。

解决：先算自由面加密区的最小网格尺寸和局部流速，反推满足 Co < 5 的时间步，把 Δt 压到 0.008 s 以内；如果算例本身对流很强，把界面格式换成 Compressive，它对大库朗数的容忍度更好。两者一起改，界面厚度能回到 2～3 个网格以内。

5.4 现象：6DOF 放开后船体原地打转或者直接沉底

原因是 UDF 里的约束宏逻辑搞反了，或者转动惯量量纲错了。约束宏 TRUE 表示锁死，很多人按“TRUE 表示生效”的直觉填，结果把想放开的垂荡锁了，把想锁死的艏摇放开了，船体在流场里乱转。

解决：在 UDF 里加一行打印，第一个时间步把六个自由度的约束状态和惯量值输出到控制台，逐个核对。惯量的核对办法是跟理论估算比：把船模当成均匀长方体，Izz ≈ m·(L² + B²)/12，L = 4 m、B = 0.6 m、m = 12.5 kg，算出来约 17 kg·m²，跟 UDF 里填的值差一个量级以上就得查单位。

5.5 现象：出口附近出现回流，阻力曲线缓慢向上漂移

原因是计算域下游太短，或者压力出口离艉部涡脱落区太近。回流把已经耗散的涡又卷回计算域，相当于给船体加了一个额外的推力或者阻力。

解决：把下游边界从 3 Lpp 推到 5 Lpp，或者在出口前加一段多孔介质阻尼区把回流压住。改完之后看一眼出口面的速度矢量，回流区面积占出口面 5% 以内算合格，超过 10% 就得继续往外推。这条改起来网格量涨得不多，因为下游区域的网格可以给得很粗。

6. 验证与提速：三套网格、一次外推和一个批量提速的习惯

算例跑通不等于结果可信。船体瞬态仿真的验证我固定走三步：网格无关性、时间步无关性、试验对比外推。三步都过一遍，一个工况大概多花两天，但能省掉后面反复返工的时间。

第一步网格无关性，划三套网格，细中粗的网格量比大约 4:2:1，基础尺寸按 2^(1/3) 的比例缩放。三套都跑同一个工况，取最后 5 个船长通过时间的阻力平均值，用网格收敛指数（GCI）评估。三套网格的阻力差如果落在 2% 以内，中网格就可以作为后续批量工况的基准；如果粗网格和细网格差 8% 以上，说明加密方向不对，得回去查是自由面没加密够还是边界层没托住。

第二步时间步无关性，用基准网格跑三个时间步，比如 0.02 s、0.01 s、0.005 s，看阻力的时间平均值和脉动均方根怎么变。平均值变化小于 1%、脉动均方根变化小于 5%，就可以把中间那个时间步定为标准步长。这一步很多人跳过，直接用经验值，结果批量工况里有的算例震荡、有的算例平稳，还以为是工况差异，其实是时间步和当地网格不匹配。

第三步试验对比，缩尺船模的试验阻力要外推到实船尺度再比。用 ITTC 1978 方法做三因次换算：

Ct,s = Cf,s + (Ct,m - Cf,m) + ΔCf

其中 Cf 用 ITTC 1957 公式，ΔCf 取粗糙度补偿系数，缩尺船模的试验值先扣掉摩擦阻力分量，剩下的兴波阻力和形状阻力按尺度因子换算。数值结果也走同一套流程，两边在同一尺度上比，总阻力系数差控制在 3% 以内算合格。直接把模型尺度的数值结果和实船试验值比，差 15% 都算正常，但这种比较没有意义。

提速方面，我一般这么做：先用定常或者准定常（大时间步、冻结 6DOF）把流场跑出一个合理的初场，再切到瞬态细步长。这样能省掉瞬态启动段那 2～3 个船长通过时间的迭代。批量工况用 journal 脚本串起来，每个工况算完自动导出阻力、升沉、纵摇的时间序列，不要等全部算完再一个个点开后处理，那时候你已经记不清哪个工况对应哪组参数了。

还有一个小技巧：把自由面监测点的波高时间序列导出来做 FFT，跟造波边界的理论频率对一下。如果主频对不上，说明入口波形在传播过程中被数值耗散改了，这时候别急着调求解器，先回去查自由面加密和消波区设置。这个检查花五分钟，比在阻力曲线上猜半天有用。

从那以后我每次开新算例，都强制先跑一个 1 个船长通过时间的短算例，只看三件事：自由面有没有保持住、y+ 分布落点对不对、阻力曲线有没有明显的启动震荡。这三件事过了再放长算，能挡掉大部分后面要花几天才能定位的问题。希望帮到你。

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