简介面向船舶CFD仿真工程师与Fluent中高级使用者这份资料演示了利用Ansys Fluent开展船体在逆海波浪中运动的瞬态仿真流程重点涵盖升沉与纵摇两自由度响应。资源共2个文件包含Fluent案例文件.gz压缩格式与配套的UDF源码.c文件整体约42.4MB可直接加载运行并查看运动限制与数值海滩等关键设置。仿真基于VOF模型求解明渠流动通过明渠波浪边界条件生成浅层波浪并采用动态网格技术模拟Wigley船体运动同时以UDF将船体约束为四自由度以聚焦升沉、纵摇等主要动态行为。案例还启用数值海滩选项以抑制出口波浪反射有助于获得更真实的瞬态结果。目前已有162人学习下载适合希望系统掌握船舶CFD瞬态建模、动态网格配置及UDF编写思路的工程师与研究者可结合案例文件快速上手省去从零搭建仿真模型的时间。1. 用 Fluent 做船体 CFD 瞬态仿真先想清楚波浪、运动和网格三件事用 Fluent 做船体 CFD 瞬态仿真多数人第一次翻车不是模型太复杂而是把静水阻力算得挺好一放进波浪里几十步就发散。船在波浪中的垂荡和纵摇本质上是自由液面、动网格、非定常湍流三件事耦合在一起其中任何一环设置不对残差就会在某个时间步突然爆表。这篇笔记按我实际拆船体瞬态仿真项目的顺序展开先讲计算域和网格怎么布局再讲 VOF 离散步与波浪入口如何设置然后给出时间步长和六自由度参数的组合方案最后把最容易踩的几个翻车点捋一遍顺便聊聊怎么验证结果才服众。适合准备接波浪中船体运动预报的 CFD 从业者也适合静态阻力会算、一碰动网格就没底的工程师。2. 计算域与网格剖分自由液面精度和动网格变形余量2.1 数值波浪水槽的尺寸以船长 L 为基准留余量船体瞬态仿真本质上是在 Fluent 里搭一个数值波浪水槽船放在槽中入口造波出口消波。计算域尺寸选大了浪费网格选小了边界反射波会绕回船身附近把阻力曲线和运动响应全部污染掉。这一节先定边界再把网格策略跟边界联系起来。深水波波长用 Airy 波理论估算公式为import math g 9.81 # 重力加速度单位 m/s^2 T 5.0 # 波浪周期单位 s lambda_wave g * T * T / (2 * math.pi) print(f深水波波长: {lambda_wave:.2f} m)逻辑上来自相速度 c λ/T 与 c sqrt(gλ/(2π)) 联立得到 λ gT²/(2π)。代码里 T 是输入的波周期lambda_wave 是输出波长。这个值后面所有网格加密和域深度都要引用不是算完看看就完。注意 g 取 9.81如果模拟的是缩尺船模需要保持模型尺度的物理参数一致。域尺寸我的经验取值如下边界位置基准量经验取值目的入口到船首船长 L2L~3L让入射波速度剖面充分发展出口到船尾船长 L4L~6L尾流波浪自然衰减避免反射侧边到船中纵剖面船长 L1.5L~2L侧反射影响降到可接受范围底部水深波长 λ0.5λ~1.0λ满足深水波假设防止底边界干扰举个例子船模长 3 m波周期 1.5 s对应的深水波长约 3.5 m那入口就要留 6~9 m出口留 12~18 m。如果实船长 100 m波周期 8 s波长接近 100 m出口要留 400~600 m。很多人在静水算例里用过的紧凑计算域到波浪瞬态算例里完全不行因为这时代价不再只是阻力而是波浪反射的物理问题。2.2 网格分级加密自由液面、近船体与背景区计算域大并不等于网格可以铺得均匀。瞬态波浪算例里的网格资源要集中在三个地方自由液面、船体边界层和动网格变形区。自由液面区域的单元尺寸要按波长锚定而不是按船长。自由液面的垂直方向单元尺寸取 λ/80~λ/100水平方向取 λ/100~λ/120。波高方向的经验是至少 10~20 个单元覆盖一个波高这样自由液面重构后不会出现明显锯齿。以 λ3.5 m 为例垂直单元尺寸约 0.035~0.044 m水平方向约 0.029~0.035 m。如果只做到 λ/40波浪传播两三个波长后波高可能只剩一半数值耗散大到没法看。各区域网格策略可以这样分配区域特征尺寸增长率关注点自由液面核心区垂直 λ/80~λ/100水平 λ/100~λ/1201.05波浪衰减控制在 5% 以内船体近场与边界层第一层由 y 决定层数 8~121.1~1.2阻力和升沉力数值稳定动网格变形区与背景网格平滑过渡避免突变的细长单元1.2防止负体积保证网格质量背景区可放到 λ/40 或更粗1.3~1.5减少网格总量加快迭代这里的增长率是指网格尺寸从加密区向外逐层放大的比例。增长率太大加密区外的网格过渡太急波浪通过界面时会发生反射和畸变增长率太小网格总量指数上升每步时间步长也被拖小。我用 1.05 做自由液面区其余区域用 1.2 到 1.5是在精确度和算量之间的折中。2.3 船体壁面 y 与第一层网格估算瞬态计算里边界层的处理决定了摩擦阻力和边界层分离是否可信。多数船体 CFD 算例我选 y≈30~300 这个范围配合壁面函数。若使用 k-omega SST 模型并开启增强壁面处理原则上可以做到 y≈1 附近但瞬态算例中边界层网格过密会显著拉低时间步长所以实际工程中以壁面函数路线为主。估算第一层网格高度可以用平板摩擦系数近似rho 998.2 # 水密度单位 kg/m^3 nu 1.0e-6 # 运动粘度单位 m^2/s U 5.0 # 航速单位 m/s L 10.0 # 船长单位 m Re U * L / nu # 雷诺数 Cf 0.0576 / Re**0.2 # 湍流平板摩擦系数估算 tau_w 0.5 * rho * Cf * U**2 u_tau (tau_w / rho)**0.5 y_plus_target 30 y_wall y_plus_target * nu / u_tau print(f第一层网格高度: {y_wall * 1000:.2f} mm)这个计算的关键是把 y 定义反解成实际第一层网格高度。代码中 Cf 是平板摩擦系数tau_w 是壁面切应力u_tau 是摩擦速度y_wall 是预计的第一层网格离壁面的距离。实际项目里网格生成时如果直接按经验拍一个尺寸很可能离目标 y 差好几倍。我的习惯是先做一次二维或粗网格试算拿到平均壁面切应力后再回代修正第一层高度。2.4 动网格还是重叠网格选型与纵摇 UDF 示例船体在波浪中要运动垂荡和纵摇是主导模态。Fluent 里实现船体运动有两条路线动网格加六自由度模型或者重叠网格。动网格路线是让船体边界随时间移动同时用扩散平滑和重构控制网格畸变重叠网格则是让船体及附近网格嵌入背景网格船体运动时背景网格完全不动。我做大幅运动时更倾向重叠网格但设置复杂度高而且部分 Fluent 版本中重叠网格与 Open Channel Wave 边界同时使用时会引出兼容性问题。比如有的版本里 Open Channel Wave 入口与重叠区域插值互为干扰波浪在入口附近就出现局部畸变。所以我的一般做法是中等运动幅度优先动网格算得稳、设置快运动幅度很大再切重叠网格前提是先用规则波空流场验证波能正常传播。以动网格强制纵摇为例可以加载一个 CG_MOTION UDF#include udf.h DEFINE_CG_MOTION(ship_pitch, dt, vel, omega, time, dtime) { real wave_freq 2.0 * M_PI * 0.5; /* 波浪圆频率弧度/秒 */ real pitch_amp 0.05; /* 纵摇幅值单位弧度 */ omega[1] pitch_amp * wave_freq * cos(wave_freq * time); }这段代码的逻辑是动网格求解器每个时间步调用一次 DEFINE_CG_MOTION将 omega[1]绕船体 Y 轴的角速度赋值为余弦形式。对角速度积分后纵摇角随时间按正弦变化。注意 UDF 里不需要写运动方程Fluent 会自行积分速度得到位移。若不是强制运动而是自由响应就改用内置 Six DOF 求解器在面板里填入质量、重心坐标和惯性矩Fluent 根据流场压力计算船体加速度和位置。这里的坑在于惯性矩至少要保证 Iyy绕横轴准确重心坐标必须基于动网格区域的局部坐标系填错的话算出来的纵摇固有周期会跟理论值差很远。动网格里的 Diffusion 参数也要单独调主要控制网格形变从船体向边界传播的速度。距离型扩散函数加上扩散系数 1.5 到 2.0是我常用的组合。扩散系数太小船体附近网格挤压过猛太大整个计算域网格都在漂移自由液面会跟着产生虚假波。3. VOF 模型与波浪生成波形入口、离散格式和湍流选型的配合3.1 VOF 体积分数输运与 Geo-Reconstruct 的底层逻辑船体瞬态仿真的自由液面只能用多相流模型描述。Fluent 里处理船体绕流的标准方案是 VOFVolume of Fluid其核心思路是给每个单元一个相体积分数 alpha。alpha1 表示该单元全为水alpha0 表示全为空气0 到 1 之间的单元就是自由液面界面单元。VOF 方法求解的是体积分数的输运方程数值格式天生会带着耗散如果不加控制界面会被抹成一条模糊过渡带波高和波形都会失真。Fluent 中解决这个问题的关键是离散格式选 Geo-Reconstruct。这个格式用分段线性方法重构界面把自由液面厚度压缩在 1~2 个单元范围内对波浪传播和船体附近的碎波捕捉都比默认的快速格式好得多。界面设置路径如下进入 Multiphase 模型面板选择 Volume of Fluid勾选 Open Channel Flow主相设为空气次相设为水Solution Methods 中把 VOF Scheme 从默认值改为 Geo-Reconstruct同时 Body Force Formulation 选 Implicit Body Force。前两项好理解第三项是数值稳定性的核心。第四项很多人忽略但瞬态计算中如果不勾选隐式体积力静水压力梯度和体力之间的平衡会在初始化阶段出现问题船体周围流场会出现虚假速度。3.2 Open Channel Wave 入口规则波的最省事方案波浪入口Fluent 提供了一个专门的 Open Channel Wave 边界条件这也是我推荐新手先试的路径。设置时在入口边界条件里选择 Multiphase 标签页下的 Open Channel Wave填入波高、波长、周期和水深再选波浪理论Airy 或 StokesFluent 会自动生成入口处的波面高度和速度分布。这套方案适合规则波特别是实验室船模试验中最常见的正弦规则波。需要注意几个隐含前提波陡不能太大Open Channel Wave 内置理论在中等波陡下还能用强非线性波就需要 UDF 或表达式自己铺入口水深必须满足波浪理论的适用范围对于深水波默认取水深大于 0.5λ 即可入口处湍流量也要单独设置湍流强度和粘度比默认值往往偏小不修正的话波面附近湍动能会被压低。3.3 表达式入口不规则波与多频率叠波规则波之外如果想做多频率叠加的不规则波Open Channel Wave 就不够用了一般用 Fluent Expression Language表达式在速度入口自定义波形。下面是一个 Airy 深水波的速度入口示例# Fluent Expression Language wave_amp 0.05 # 波幅单位 m omega 2 * PI / T # 圆频率单位 1/s k 2 * PI / lambda_wave # 波数单位 1/m u_inlet wave_amp * omega * exp(k * (z - water_level)) * cos(k * x - omega * t) w_inlet wave_amp * omega * exp(k * (z - water_level)) * sin(k * x - omega * t)这里的 u_inlet 是水平方向速度分量w_inlet 是垂向速度分量。指数项 exp(k * (z - water_level)) 保证速度随水深呈指数衰减符合深水波特性。z 减 water_level 一定是负数若写成正数速度会随高度发散。表达式的优势是可以把多个不同周期不同幅值的波简单相加模拟不规则海况的初步近似。缺点是调试比 Open Channel Wave 麻烦每次修改都要检查入口截面的速度分布不能只看波面云图。3.4 湍流模型配合k-omega SST 为什么排在首位瞬态波浪算例里湍流模型的选择对波浪衰减影响很大。我在船体瞬态算例里的默认选择是 k-omega SST它在逆压梯度区域和边界层分离的预测上比标准 k-epsilon 准确对升沉和纵摇的阻力响应更接近实验值。代价是计算量稍大而且对入口湍流参数更敏感。如果用 k-epsilon Realizable计算会便宜一些但对波浪的数值阻尼更明显。具体表现是波浪传播几个波长后波高明显降低船首附近的波面形状和实际不符。这个问题的根源是湍流模型的湍动能生成项在自由液面附近被高估把波能转化为湍动能。SST 在这一区域的耗散控制得更保守所以在船体瞬态仿真中优先选用。如果真的只能用 k-epsilon需要在入口把湍流强度调低到 1% 以下粘度比设置为 0.1 左右才能减小人为阻尼。4. 瞬态求解器参数时间步长、六自由度与松弛因子的组合策略4.1 时间步长的双约束CFL 条件与波浪周期瞬态船体仿真里时间步长受两个条件同时约束。一是局部 CFL 数二是每个波浪周期内至少要覆盖足够多的物理时间步。两者取小值才是可用时间步长。CFL 条件定义为delta_x 0.03 # 自由液面最小网格尺寸单位 m U_char 5.0 # 特征速度航速 波速单位 m/s CFL_max 1.0 delta_t_cfl CFL_max * delta_x / U_char print(fCFL 约束的时间步长: {delta_t_cfl:.5f} s)这里 delta_t_cfl 约等于 0.006 s。如果网格加密到 0.02 m时间步长就要降到 0.004 s 左右。第二个约束来自波浪周期。以波周期 T5 s 为例一个周期内建议至少覆盖 200 个时间步对应的步长约 0.025 s。如果按 CFL 算出来是 0.004 s那就取 0.004 s不要为了省步数硬放大步长。经验规律是当时间步长超过 T/150 时波浪入口产生的波面会和船体运动解耦船体响应出现数值振荡。4.2 六自由度运动设置质量、重心、惯性矩和浮力初始化自由运动的船体要在 Fluent 的 Six DOF 面板里设置物理属性。这块界面看起来简单填错一个量整个运动就错了。参数说明常见错误质量船体排水量对应的质量单位 kg直接用排水体积乘密度忘记去掉附连水重心位置基于动网格区域局部坐标系坐标原点选取不一致Iyy绕横轴的惯性矩决定纵摇固有周期用手算公式估错导致纵摇周期偏差初始波高初始化时自由液面高度与入口波面高度不一致产生瞬时冲击我一般会在初始化之前先做一次静水下的船体阻力计算把浮力验证一遍。具体做法是初始化后在 Report 里读出船体表面受到的总浮力和重力比较误差超过 1% 就检查压降初始化是否偏离静水压力分布。常见现象是船体重心偏高、浮力匹配不上船在计算前几步快速下沉之后整个流场就乱了。4.3 松弛因子与内部迭代次数不要迷信默认值瞬态船体算例的松弛因子和稳态算例不一样。Fluent 默认的动量松弛因子 0.7、压力 0.3在一开始运行时可能偏激进。我的常用区间是压力 0.2~0.3动量 0.5~0.6VOF 方程 0.2~0.25k 和 omega 各 0.6 左右。调低松弛因子会牺牲每个时间步内的收敛速度但可以让每个步内迭代更稳定避免残差在波浪作用期反复震荡。瞬态计算每个时间步内部建议跑 10~20 次迭代。如果 20 次迭代内残差仍不下降不要无限加迭代次数优先检查是否有网格负体积或波面附近过大压差。残差下降缓慢的一个常见被忽视原因是初始化没有采用体积力平衡设置流场中还残留静水压力梯度前几十步都在“消化”这个假压力。4.4 数据监测与自动保存算到一半不手忙脚乱瞬态计算长时间运行监测和保存策略要提前设计。我在开始算之前就会把三个监测点建好船体垂荡力、纵摇力矩、一个固定的波高探针。波高探针放在船首前一个船长距离处用于判断入口波浪是否按预期到达船体。把这三个量的时间曲线输出比看残差更直观。自动保存频率要结合网格量和时间步数确定。每 100 步保存一次通常太频繁写盘开销大每 500~1000 步保存一次比较稳妥万一计算中途崩了损失最多几百步数据。如果网格量在千万级保存文件体积大建议同时开自动导出摘要文件只记录监测量和时间文件很小用于崩溃后快速判断是在哪个阶段出问题。5. 高频翻车排查发散、波衰减、下沉和回流四个方向5.1 现象算到某个时间步压力残差突然爆表原因最常见的是局部网格变形过大船体附近的单元被拉成长条或压出负体积其次是时间步长偏大压力速度耦合迭代在波峰附近无法收敛。某一次我算垂荡幅值较大的工况在前 30 步残差都很平滑到第 31 步压力残差从 1e-3 跳到 1e3就是因为船体抬升瞬间局部网格被拉爆。解决第一步把时间步长减半看是否还在同一位置爆残差若依旧打开动网格设置里的 Remeshing把最小长度尺度和最大偏斜度重新设置让等于或超过阈值的单元被自动重构。这里的参数不是越大越好最小长度尺度设得太大重构频繁每步网格数变化剧烈收敛反而差。5.2 现象入口波浪沿传播方向衰减严重船前一个船长处的波高只剩六成原因自由液面网格太粗或者动量离散格式用的阶数偏低。瞬态计算里如果只开 First Order Upwind数值耗散会把波浪能量抹掉。另一个原因是湍流模型对自由液面附近涡粘过度波浪能量被湍流耗散。解决把自由液面区域网格细化到 λ/100 水平同时将动量和 VOF 的离散格式都改成 Second Order UpwindVOF 保持 Geo-Reconstruct。再把湍流模型的入口粘度比下调到 0.1 左右。做完这三处修改波高衰减通常能控制在 5%~10% 以内。5.3 现象船体在瞬态启动后持续下沉几分钟内垂向位移一直在增大原因这是初始化阶段浮力没有平衡。静水力压初始化用的是默认混合初始化没有按 Open Channel Flow 的静水压力分布赋值船体底部的压力场不等于实际应受到的浮力数值把船往下压。解决初始化方式改为 Hybrid Initialization并勾选 Open Channel Initialization让 Fluent 根据自由液面位置初始化静水压力分布。之后进入 Report 检查浮力确保与船体重力之差在 1% 以内。这个步骤我每次都在算瞬态之前强制做一遍相当于数据保底。5.4 现象出口附近回流明显尾流区域压力出现异常波动原因计算域出口距离船尾太近尾流尚未充分发展就被出口边界强行截断压力被反射回船体后侧。常见是把出口设置在 2L 附近以为足够了实际上 4L 都没办法完全消除反射。解决把出口延长到 5L 以上同时在这个延长区内逐步拉大网格尺寸人为制造一个数值阻尼区让波浪在到达出口前自然衰减。还可以在出口边界把湍流强度上调一点帮助波浪耗散但不要设得过高否则尾流数据失真。5.5 现象动网格区域某处始终出现负体积计算直接终止原因船体运动幅值超过了动网格变形区能承受的极限网格单元被压扁。常见是扩散系数设太小变形全集中在船体附近或者 Remeshing 的最小长度尺度设得比背景网格尺寸大重构时无法生成有效单元。解决把扩散系数调到 2.0 左右并选择基于距离的扩散函数。同时把 Remeshing 的最小长度尺度设置成背景网格尺寸的 0.3 倍左右最大偏斜度设为 0.9 以下。如果仍不稳定还需要增大动网格变形区范围让变形分摊到更大区域。6. 结果验证与批处理技巧RAO 对比和算量账本验证的第一步是在正式波况计算前先跑一个静水阻力算例。静水算例比瞬态算例简单如果静水阻力的曲线和实验值差超过 5%后面波浪中的运动响应没有参考价值。静水验证通过后再缓缓把波高从小幅值逐步抬高。不要一上来就上大波高因为船体的非线性响应在小波高下更容易和线性理论或实验数据对比定位问题更方便。后处理阶段的核心是把船体的垂荡和纵摇时间序列转成 RAOResponse Amplitude Operator。做法是取船体重心垂向位移的时间曲线做傅里叶变换在波浪频率处提取幅值除以入射波幅。纵摇也用同样方法。做这一步骤时有一个关键习惯一定要把时间序列截取到船体运动达到稳态之后的波段前几个周期的瞬态响应不能混入傅里叶变换。我一般会在达到稳态后至少保留 10 个波浪周期的数据。对比线性理论值时垂荡和纵摇 RAO 在长波区都应趋近于 1短波区衰减。如果算出的 RAO 在波浪频率附近有异常峰值基本可以判断是数值共振问题出在时间步长或网格分辨率上不是真实物理现象。批处理多工况时我习惯先在 Fluent 里手动算通一个波况然后导出 journal 文件替换波浪参数后批量提交。这样同一个几何网格下一套流程算多个波高、多个周期的工况效率提高明显。另一个常用的技巧是把监测数据和波高探针的位置写入表达式这样每次更换工况只需修改表达式里的波高和周期两个变量其余保持不变。从那以后我每次做船体瞬态仿真都会在算正式工况前强制走一遍静水验证和浮力检查确认没问题再切到波浪环境。这套流程帮我省掉了大量半夜翻车重算的时间。希望帮到你。本文还有配套的精品资源点击获取