简介面向船舶流体力学学习者与工程人员的 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 做基准取值区间如下表缩尺船模和实船尺度都适用只是实船尺度下要额外注意边界层网格量会爆。方向建议范围理由上游入口到艏柱1.02.0 Lpp让来流有距离发展避免入口面直接压在兴波区艉柱到下游出口3.05.0 Lpp艉部伴流和兴波需要空间耗散短了出口回流明显中纵剖面到侧边界1.52.0 Lpp侧向阻塞效应控制在 1% 以内水线面到顶部边界0.51.0 Lpp兴波高度有限太高压根没网格白占资源水线面到底部边界1.01.5 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 的区域整体设为水的体积分数 1z 0 设为 0界面附近留一层过渡。这一步漏了整个域都是空气船体等于在天上飞阻力曲线会从第一个时间步就开始发散。界面重构格式有三个常用选项选型理由如下格式适用场景代价Geo-Reconstruct自由面清晰、网格较规整船体静水阻力首选每步计算量大界面附近网格质量差时容易出毛刺Compressive网格偏多面体、界面容易拉花时的稳健选择界面略厚兴波细节会被抹一点CICSAM强对流、自由面破碎明显如艏部翻卷收敛慢时间步要压得更小表面张力在船体阻力算例里通常可以关掉船模尺度的韦伯数在 10⁴ 量级表面张力对总阻力的贡献远小于 1%但如果你的算例要复现艏部翻卷、甲板上浪这类现象就得开 CSF 模型并把表面张力系数设成 0.072 N/m。开与不开的差别在阻力上体现不出来在自由面形态上很明显别拿阻力曲线去判断该不该开。时间格式先用一阶隐式跑通等流场稳定、自由面形态合理之后再切二阶隐式重跑。上来就二阶前期残差震荡会让你分不清是格式问题还是网格问题。2.3 湍流模型与壁面处理k-ω SST 还是 Realizable k-ε船体绕流的湍流模型选型核心看两件事艉部逆压梯度区的分离预测以及壁面处理方式和你的 y 是否匹配。k-ω SST 在逆压梯度下表现稳艉部分离点和伴流宽度跟试验对得比较好是船体阻力算例的默认选择。它对壁面距离敏感低雷诺数处理要求 y 在 1 附近壁面函数处理要求 y 落在 30100中间那段过渡区5 y 30是它的软肋尽量别落进去。Realizable k-ε 更稳健收敛快对网格质量的容忍度高适合前期试算和批量工况筛选但它对分离区的预测偏保守艉部伴流容易算厚阻力偏高个 3%5% 是常有的事。如果算例涉及模型尺度下的层流到湍流转捩艏部前缘一段可以考虑 γ-Re θ 转捩模型代价是网格得按低雷诺数要求做第一层高度要压到 y 1网格量至少翻一倍。多数工程算例不做这一步直接用全湍流然后在结果里注明这个假设。壁面处理方式在 Near-Wall Treatment 里选选 Standard Wall Functions 就得保证第一层网格落在对数律区选 Enhanced Wall Treatment 就得保证 y 足够低。两者混用是最常见的错误之一后面避坑章节会展开。求解器类型选压力基Pressure-Based瞬态压力-速度耦合用 PISOPISO 比 SIMPLE 多一次邻接修正时间步推进时收敛更稳代价是每步多花约 30% 的时间。PISO 的 skewness correction 和 neighbor correction 各设 1 就够设大了收敛反而慢。3. 网格不是越细越好船体表面网格、边界层与 y 的账网格这一章的核心不是“划多密”而是“密在哪、第一层多高、质量指标守住哪几条线”。船体瞬态算例的网格量通常在 300 万到 800 万之间超过一千万而阻力还在漂八成不是网格量的问题。3.1 全局网格与局部加密区的划分船体这类细长体全局用切割体网格cutcell或者 trimmer 比较省事船体表面贴一层棱柱层自由面附近做各向异性加密。加密区我一般分四块第一块是船体表面按船长给基础尺寸Lpp 4 m 的船模表面网格尺寸给 0.020.04 m曲率大的地方艏部、艉部、舵前缘用 curvature 或者 proximity 再细化一档。第二块是自由面水线面上下各 0.30.5 倍波高做体加密网格尺寸取船体表面尺寸的 0.5 倍。自由面加密不够兴波波形会明显偏平波高被数值耗散吃掉。第三块是艉部尾流区从艉柱往后 23 Lpp、横向 0.5 Lpp 的盒子里给一档中等加密网格尺寸取表面尺寸的 1.52 倍。伴流场的涡结构靠这块网格托住。第四块是艏部兴波区艏柱前方 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/μ反推目标 y第一层高度 y壁面处理方式11.6×10⁻⁵ m0.016 mmEnhanced Wall Treatment低雷诺数304.9×10⁻⁴ m0.49 mmStandard Wall Functions 下限508.2×10⁻⁴ m0.82 mmStandard Wall Functions 常用1001.6×10⁻³ m1.6 mmStandard Wall Functions 上限再高就丢分离细节棱柱层增长比给 1.151.2层数 1015 层总厚度要盖住边界层厚度。边界层厚度按 δ ≈ 0.37·Lpp·Re^(-1/5) 估约 0.037 m棱柱层总厚度做到 0.040.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 时间步长怎么定库朗数与通过时间双约束船体瞬态算例的时间步有两个约束取小的那个。第一个是船长通过时间约束一个船长被船体走过的过程至少要分 100200 个时间步否则船体运动的历史效应被跳步跳过。Δt ≤ Lpp/(U·N)N 取 150Lpp 4 m、U 1.5 m/s 时 Δt ≤ 0.0178 s。第二个是库朗数约束自由面附近的网格最密Δx 最小局部库朗数 Co U·Δt/Δx 建议控制在 15。自由面加密区网格尺寸 0.01 mU 1.5 m/s 时Co 5 对应 Δt 0.033 sCo 1 对应 Δt 0.0067 s。两个约束一夹Δt 落在 0.0050.015 s 是比较合理的区间。带强自由面破碎或者 6DOF 大幅运动的算例往 0.005 s 压只算静水阻力、船体固定不动的算例0.010.015 s 够用。全局库朗数在求解器里设 25 作为自适应上限让它自己在局部自动缩步。自适应缩步不是免费的缩得太频繁会让总计算时间翻倍所以自定义时间步的上限还是按上面的估算给自适应只当兜底。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 支持一阶到五阶斯托克斯波和椭圆余弦波参数表如下参数取值示例说明波高 H0.08 m按缩尺比换算对应实船有效波高周期 T1.1 s与波长一起决定色散关系水深 d2.5 m深水波 d/L 0.5波长 L1.89 m由色散关系反算深水 L gT²/(2π)波向沿 -X与船体航向相反入口面设成 Open Channel Wave顶部设压力出口并把 open channel 选项打开这样气相能自由进出。造波边界下面要留够 12 个波长的过渡区波形在这个距离内从理论解过渡到数值解太短了波还没长起来就撞到船体上。消波放在出口段用数值海滩numerical beach或者多孔介质阻尼区。数值海滩的长度取 1.52 个波长阻尼系数从 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 以下只是必要条件真正的判据是阻力曲线的滑动平均在最近 35 个船长通过时间里变化小于 1%升沉和纵摇的振幅包络稳定自由面监测点的波高不再漂移。瞬态算例的残差曲线本来就是震荡的拿定常那套残差标准去卡永远等不到“收敛”。统计平均的窗口也要说清楚阻力曲线前 23 个船长通过时间是启动段不能进平均取后面 58 个通过时间做时间平均同时记录脉动均方根。只报一个平均值不报脉动等于把瞬态算成了定常。5. 避坑与排查五个把我按在椅子上重算的现场这一章的五条都是实际算例里撞出来的按“现象 → 原因 → 解决”写顺序大致对应算例推进的先后前面几条偏设置后面几条偏结果。5.1 现象初始化后第一步就发散残差冲到 1e3原因基本锁定在自由面初始化。整个域被默认填成水或者填成空气船体初始湿表面跟实际差了一个船宽第一步的压力场修正量巨大求解器直接崩。解决初始化之后手动 patch 区域。用 Mark/Adapt 按 z 坐标切出 z 0 的体区域patch 成水的体积分数 1z 0 patch 成 0界面上下各 0.5 倍波高的薄层 patch 成 0.5 作为过渡。patch 完之后在初始化面板里检查水的总质量跟理论排水量对一下差太多说明切区域的时候把船体内部也算进去了得先把船体实体区域从 patch 范围里减掉。5.2 现象阻力比试验值系统性高 10% 以上加密网格只降 1%原因多半是壁面处理方式和 y 不匹配。第一层高度给到 0.5 mmy 落在 10 附近正好在 k-ω SST 的过渡区壁面函数用不上、低雷诺数解析又不够壁面切应力被高估阻力整体抬高。解决先在后处理里画 y 沿船长的分布确认实际落点。要么把第一层压到 0.02 mm 走低雷诺数路线网格量涨 30%要么把第一层放到 1 mm 以上走壁面函数路线阻力精度略降但稳定。两条路选一条不要停在中间。改完之后阻力通常会往试验值靠 58 个百分点。5.3 现象自由面算着算着变成一坨模糊的过渡带兴波波形消失原因是界面格式和时间步不匹配。用 Geo-Reconstruct 配 0.03 s 的时间步自由面处库朗数到了 15 以上界面重构来不及数值扩散把界面抹厚了。解决先算自由面加密区的最小网格尺寸和局部流速反推满足 Co 5 的时间步把 Δt 压到 0.008 s 以内如果算例本身对流很强把界面格式换成 Compressive它对大库朗数的容忍度更好。两者一起改界面厚度能回到 23 个网格以内。5.4 现象6DOF 放开后船体原地打转或者直接沉底原因是 UDF 里的约束宏逻辑搞反了或者转动惯量量纲错了。约束宏 TRUE 表示锁死很多人按“TRUE 表示生效”的直觉填结果把想放开的垂荡锁了把想锁死的艏摇放开了船体在流场里乱转。解决在 UDF 里加一行打印第一个时间步把六个自由度的约束状态和惯量值输出到控制台逐个核对。惯量的核对办法是跟理论估算比把船模当成均匀长方体Izz ≈ m·(L² B²)/12L 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把流场跑出一个合理的初场再切到瞬态细步长。这样能省掉瞬态启动段那 23 个船长通过时间的迭代。批量工况用 journal 脚本串起来每个工况算完自动导出阻力、升沉、纵摇的时间序列不要等全部算完再一个个点开后处理那时候你已经记不清哪个工况对应哪组参数了。还有一个小技巧把自由面监测点的波高时间序列导出来做 FFT跟造波边界的理论频率对一下。如果主频对不上说明入口波形在传播过程中被数值耗散改了这时候别急着调求解器先回去查自由面加密和消波区设置。这个检查花五分钟比在阻力曲线上猜半天有用。从那以后我每次开新算例都强制先跑一个 1 个船长通过时间的短算例只看三件事自由面有没有保持住、y 分布落点对不对、阻力曲线有没有明显的启动震荡。这三件事过了再放长算能挡掉大部分后面要花几天才能定位的问题。希望帮到你。本文还有配套的精品资源点击获取