船舶在航行时，若发生大幅横摇运动，不仅会降低船员的舒适性, 影响船上设备的正常运转，严重时还会出现翻船的事故，从而造成巨大的生命财产损失。

通常认为运动阻尼主要由船体的辐射波所致。相较其他自由度运动，横摇受黏性效应的影响最大，因此仅考虑辐射阻尼的势流方法存在预测精度低的问题^[1]。横摇运动的准确预报依赖于横摇阻尼系数的获取。模型试验是确定横摇阻尼最直观且可靠的方法，但存在试验成本高的问题。除此之外，常用的数值评估方法主要有经验公式方法和计算流体动力学(computational fluid dynamics，CFD)方法。

横摇阻尼可主要分为：兴波阻尼、升力阻尼、摩擦阻尼、漩涡阻尼以及舭龙骨阻尼^[2]。具体而言，兴波阻尼和升力阻尼与横摇角速度呈正比，升力阻尼随着航速的增加而增加。摩擦阻尼、漩涡阻尼和舭龙骨阻尼与横摇角速度平方呈正比，摩擦阻尼的占比很小，舭龙骨阻尼的占比最大。针对这5个部分，Ikeda等^[3-5]根据船模试验结果总结出了详细的经验公式进行计算，后续进行了相应的简化^[6]。尽管此类方法具有快捷、实用的优势，但对船宽、吃水、重心高、附体尺寸等参数有限制，例如对某些重心位置较高、中横剖面系数较小的船型并不适用，且无法得到横摇时历曲线。而CFD方法可以准确评估流体黏性对横摇运动的影响，且不受船体几何参数的影响，具有良好的适用性。通过求解Navier-Stokes方程，获取丰富的运动、流场相关信息，对非线性问题的求解也有较高的准确度。近年来，CFD法在船舶横摇运动分析中，得到了越来越广泛的应用^[7-13]。

相较零航速，分析船舶在有航速工况下的横摇自由衰减对认清其在实际营运过程中横摇的产生机制，探索减摇方法，具有更重要的意义。现有研究中，Yang等^[7]模拟了DTMB5512船模在2个航速(Fn=0.138和Fn=0.280)下的横摇自由衰减过程，横摇周期和阻尼系数均与试验结果吻合较好。Irkal等^[10]针对某集装箱船模的数值模拟结果表明，随着航速增加到Fn=0.240，横摇阻尼系数单调上升，主要受升力的影响；而横摇周期和附加质量呈非单调变化。罗天等^[12]采用自主开发的naoe-FOAM-SJTU求解器模拟了DTMB5512船模的自由衰减和强迫横摇运动。结论表明兴波阻尼与漩涡阻尼对横摇运动的影响较大，且舭龙骨在一定的横摇角度范围内可以减小总阻尼力矩。马山等^[13]基于重叠网格技术对C11集装箱船在零航速和有航速下的横摇自由衰减进行了模拟，横摇衰减速度随航速的增加而增加。与零航速相比，Fn高于0.1后，网格密度不足以捕捉舭龙骨附近的流场变化，最终导致CFD模拟结果与试验结果符合不理想。

综上所述，相关研究集中在对中低航速工况横摇自由衰减运动和流场特性的分析，而对高航速下横摇阻尼系数的变化规律鲜有报道。另一方面，未关注初始横倾角度的影响，无法评估初始扰动对后续横摇运动的影响。本文针对DTMB5512船模，采用重叠网格方法模拟了其在3个航速下的横摇自由衰减运动，覆盖了低、中、高航速，每个航速分别对应3个初始横倾角，详细地分析了航速和初始横倾角对于横摇阻尼系数的影响。最后根据短期预报结果来评估CFD模拟结果用于横摇运动分析时的可靠性。

1 数值方法 1.1 数学方程

考虑不可压流动，流体密度ρ为常数，在笛卡儿坐标系下，采用Reynolds平均法对脉动项进行处理，控制方程可表示为^[14]：

$ \frac{\partial}{\partial x_i}\left(\rho u_i\right)=0 $

(1)

$ \frac{\partial}{\partial t}\left(\rho u_i\right)+\frac{\partial}{\partial x_j}\left(\rho u_i u_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial}{\partial x_j}\left(\mu \frac{\partial u_i}{\partial x_j}-\rho \overline{u_i^{\prime} u_j^{\prime}}\right)+F_i $

(2)

式中: u_i、u_j为速度分量，μ为动力黏度系数，$\rho \overline{u_i^{\prime} u_j^{\prime}}$为Reynolds应力，F_i为体积力。采用有限体积法离散控制方程，对对流项的离散采用二阶格式，对时间项的离散采用隐式不稳定方法，对压力-速度耦合方程的求解采用基于SIMPLE算法的分离流求解器，最后采用动态流体固体相互作用(dynamic fluid body interaction，DFBI)模型求解刚体六自由度运动。

湍流模型选用Realizable k-ε模型。与标准k-ε模型相比，其湍动黏度计算公式中引入了与旋转和曲率相关的内容，且ε方程中不具有奇异性，可有效用于各种不同类型的流动模拟。

采用流体体积(volume of fluid，VOF)方法来追踪自由液面的变化，其方程可表示为

$ \frac{\partial}{\partial t}\left(f_q\right)+\frac{\partial}{\partial x_i}\left(f_q u_i\right)=0 $

(3)

式中f_q为第q相流体在单元格中所占的体积分数。

1.2 几何模型

DTMB5512船模是ITTC(International towing tank conference)推荐的专门用于CFD仿真计算的标准船模之一，水池试验数据丰富，为DTMB5415的全相似几何模型，缩尺比为1∶ 46.588。如图 1所示，船体包含带有声呐罩的艏部、船两侧的舭龙骨和1个方艉，主要参数见表 1。

1.3 网格划分

本文建立的计算域如图 2所示，包括背景区域和重叠区域两部分。重叠网格技术允许背景域和重叠域相互重叠、嵌套，通过对整个流场计算域和具有相对运动的船体单独划分网格，可有效减少网格数量，节省计算时间，并在各子网格间进行插值来传递流场信息^[15]。

计算域坐标原点位于船艉垂线与水面相交处，x轴指向船艏为正，y轴指向左舷为正，z轴垂直水平向上为正。背景域尺寸x方向为-3.0L_pp~4.0L_pp，y方向为-2.0L_pp~2.0L_pp，z方向为-2.0L_pp~L_pp。重叠域尺寸为x方向为1.5L_pp，y方向为-2.0B~2.0B，z方向为-2.0D~2.0D。对于背景域，船艏正对的边界设置为速度入口，船艉正对的边界设置为压力出口，其余边界均设置为速度入口。为避免自由面兴波反射的影响，在压力出口设置阻尼区进行消波。

网格划分采用非结构化网格，包括了表面重构模型、切割体网格生成器模型和棱柱层网格生成器模型3类。为准确捕捉自由液面和船体周围的流场变化，对舭龙骨(见图 3(a))、自由液面(见图 3(b)) 及背景域重叠域交界处(见图 3(c))进行加密。船体外壳为不可滑移壁面，沿垂直船模壁面方向，共设置6层边界层网格(见图 3(d))，拉伸率为1.1，总厚度为0.01 m。

2 结果与讨论 2.1 数值模拟结果验证

为进行收敛性验证，设置了基础尺寸为0.04、0.06、0.08 m的3套网格，对应的网格数分别为12.90×10⁶、4.41×10⁶和2.08×10⁶。图 4对比了网格基础尺寸对横摇自由衰减模拟结果的影响，图中：φ为横摇角，t为横摇运动时间。相较0.08 m，网格基础尺寸减少到0.04、0.06 m，对应的横摇衰减曲线无明显差别。

另一方面，时间步长的选取也会影响计算结果。根据ITTC指南^[16]，一个运动周期内至少包含100个时间步长。图 5对比了时间步长对横摇自由衰减模拟结果的影响。随着时间步长的减小，横摇周期略有增加，且横摇衰减速度加快。当计算时间步长减少到0.005 0、0.002 5 s时，对应的横摇衰减曲线无明显差别。

对于这类有航速问题的数值模拟，通常是给船舶施加某初始横倾角后航行一段距离再释放，以保证船体周围扰动流场由瞬态逐渐过渡到稳态^{[8, 13]}。图 6对比了释放时间对横摇自由衰减模拟结果的影响。释放时间设定为2、5、10 s对应的横摇自由衰减曲线几乎无差异。

船舶运动中，横摇会与其他自由度运动存在一定的耦合作用，进而影响模拟结果。图 7对比了考虑不同自由度运动对数值模拟结果的影响，横摇运动的幅值和周期均会随之发生变化，具体的横摇幅值见表 2。综合来看，仅考虑横摇和纵摇两个自由度运动时，计算结果与试验结果间的偏差最大。而考虑横摇、横荡和垂荡三自由度运动时，计算结果与试验结果最为接近。

综上所述，兼顾计算精度和计算效率，后续数值模拟统一考虑横摇、横荡和垂荡三自由度运动，设置网格基础尺寸0.06 m，时间步长0.005 0 s，船模释放时间2 s。

2.2 横摇阻尼分析方法

典型横摇自由衰减曲线如图 8所示，可以得到相邻半个周期内横摇幅值的平均值φ_m与差值Δφ：

$ \Delta \varphi=\varphi_k-\varphi_{k+1} $

(4)

$ \varphi_{\mathrm{m}}=\left(\varphi_k+\varphi_{k+1}\right) / 2 $

(5)

以φ_m为横坐标，Δφ为纵坐标，即可得到给定初始条件下的横摇消灭曲线。由于本文仅考虑线性阻尼和二次阻尼，根据最小二乘法可拟合得到φ_m和Δφ之间的函数关系：

$ \Delta \varphi=a \varphi_{\mathrm{m}}+b \varphi_{\mathrm{m}}^2 $

(6)

横摇运动的非线性主要是由非线性横摇阻尼力矩造成的，对应的运动方程为

$ (I+A) \ddot{\varphi}+B_{\mathrm{L}} \dot{\varphi}+B_0 \dot{\varphi}^2+C \varphi=0 $

(7)

式中：I为横摇转动惯量，A为附加横摇转动惯量，B_L为线性横摇阻尼系数，B_Q为二次横摇阻尼系数，C为横摇恢复力矩系数。进一步地，利用能量法，船舶由φ_k摇至φ_k+1的过程中，近似认为其在一个周期里作简谐振荡运动，阻尼耗散的能量等于势能的减少，可以确定阻尼系数B_L、B_Q与式(6)中消灭系数a、b间的关系^[17]为

$ \left\{\begin{array}{l} B_{\mathrm{L}}=\frac{2 D h}{\omega_{\varphi} \pi} a \\ B_Q=\frac{3 D h}{4 \omega_{\varphi}^2} b \end{array}\right. $

(8)

式中：D为船舶排水量，N；h为稳性高，m；ω_φ为横摇固有频率，s^-1。

最后将二次阻尼线性化，得到任意横摇幅值φ_i对应的等效阻尼B_eq为

$ B_{\text {eq }}=B_{\mathrm{L}}+\frac{8}{3 \pi}\left(\omega_{\varphi} \varphi_i\right) B_Q $

(9)

2.3 航速对横摇自由衰减运动的影响

根据Irvine等^[18]的测量结果，试验结果(EXP)和CFD方法模拟结果对比见图 9~11。包括了3个航速(Fn=0.138，Fn=0.280，Fn=0.410)，每个航速分别对应3个初始横倾角(φ₀=10°，φ₀=15°，φ₀=20°)。低、中、高航速工况下，模拟得到的横摇幅值均与试验结果吻合良好。随着航速的增加，横摇衰减速度变快。以20°初始横倾角为例，Fn=0.138时，第4个周期对应的幅值为6.12°，而Fn=0.410航速下该幅值已经下降到2.67°，10 s后横摇运动趋于停止。

为更好地解释横摇自由衰减现象，图 12给出了不同航速下的兴波图。在此基础上，图 13直接对比了船艉(aft)和船艏(fore)附近自由液面升高情况，对应位置的坐标分别为(-0.5，0，0)和(2.65，0.25，0)。Fn=0.138时，航速较低，由船体运动产生的辐射兴波占据主导，船艏处自由波面升高呈现典型的衰减规律(见图 13(b))，与横摇运动趋势一致。随着航速增加，水流流经弯曲的船体时，沿船体表面的压力分布不同，导致船体周围的水面升高或下降，在重力和惯性的作用下形成了船行波^[19]。Fn= 0.410时，如图 12(c)所示，船艏兴波以散波为主，船艉兴波由横波和散波混合而成，伴随自由液面的快速上升，此时船艏处兴波时历曲线的衰减特性相较低航速已不明显(见图 13(b))，因此对于高航速下横摇自由衰减问题的求解，需要考虑兴波阻尼所引入的影响。

表 3为横摇周期T_φ的试验和模拟结果对比，航速越高，T_φ越小，所有航速下试验值与模拟值间的相对误差均在6%以内，吻合良好。尽管如此，由于误差会随着时间不断累积，导致两条横摇时例曲线在运动后期产生了相位差。

2.4 航速对横摇阻尼系数的影响

根据图 9~11中横摇自由衰减曲线，可进一步处理得到不同航速和不同初始横倾角对应的横摇阻尼系数，图 14和表 4对具体结果进行了汇总。考虑到横摇衰减运动在初始阶段的不稳定性，本文均选取第3个横摇幅值作为起始分析点。Fn=0.280航速下初始横倾角为10°时的试验数据质量较差。另一方面，如图 11所示，Fn=0.410航速下横摇衰减速度较快，试验仅记录了前6 s的衰减曲线，可供分析的有效数据点较少，影响后续的曲线拟合准确性，因此对于表 3中的上述4个工况，未对其试验结果开展横摇阻尼系数分析。通过分析可以得到如下结果：

1) 相同航速时，初始横倾角从10°增加到20°，对应的横摇消灭曲线趋势接近一致，因此由曲线拟合得到的线性阻尼系数差别不大。

2) 相同初始横倾角时，线性阻尼系数随着航速的升高而逐渐增加，高航速对应的横摇消灭曲线趋于线性。航速会引入兴波阻尼和升力阻尼的影响，二者均与横摇角速度呈正比，因此增加了线性阻尼在总阻尼中的占比。

3) 尽管CFD方法可以准确地预测横摇幅值，但得到的阻尼系数与试验结果存在差异，线性阻尼和二次阻尼呈现此消彼长的关系。横摇消灭曲线的拟合对横摇幅值的敏感性较高，会将两种方法得到的幅值间差异放大，从而对阻尼系数产生较大影响。

4) 根据式(9)将二次阻尼线性化后得到的等效阻尼系数如图 15所示，随着横摇幅值的增加而增加，且试验结果和模拟结果吻合良好。相同横摇幅值，高航速对应更大的等效阻尼系数，加快了横摇衰减速度。

2.5 横摇短期预报

除CFD方法，ITH(Ikeda，Tanaka，Himeno)方法在横摇阻尼分析中也得到了广泛应用，且具有求解快速的优势。本文将试验(EXP)、CFD方法和ITH方法得到的横摇阻尼系数作为输入，采用三维频域势流理论^[20]，计算船舶在不规则波中的横摇运动响应。根据短期预报结果来评估ITH方法的适用性，并进一步验证CFD方法得到的横摇阻尼系数用于不同海况和不同航速下运动分析时的可靠性。

表 5对横摇短期预报的海况进行了汇总。长峰不规则波浪谱选用JONSWAP谱。有义波高H_s覆盖了1.25~9.00 m的宽海况范围，对应的谱峰周期T_p和谱峰因子γ根据DNV(Det Norske Veritas)规范^[21]进行选取。

横摇短期预报结果如图 16所示，计算得到的单幅有义值为所有谱峰周期和浪向对应的最大值。相同海况，航速由Fn=0.138增加到Fn=0.280时，由于等效阻尼增大，横摇运动响应受到抑制。基于试验结果和CFD方法得到的单幅有义值十分接近，而ITH方法偏差较大，且差异性在高海况下表现得尤为明显。ITH方法通过一系列经验公式，对船体长、宽、吃水、重心、舭龙骨尺寸等几何参数进行计算来求解横摇阻尼。但对于有航速的工况，该方法无法准确捕捉流场的变化，以及评估升力所引入的影响。另一方面，Kawahara等^[6]的结果表明，C_m对漩涡阻尼影响较大，船体舭部的流体分离随着中横剖面系数C_m增加而加剧，导致漩涡阻尼快速升高。并指出当C_m高于0.90时，使用ITH方法评估横摇阻尼才更为合理。而本船的C_m为0.82，可能也会造成计算精度偏低。

3 结论

1) 由CFD方法得到的横摇幅值及横摇周期与试验结果吻合良好。航速由Fn=0.138增加到Fn= 0.410，横摇衰减速度明显变快，伴随着船体周围兴波的快速上升。

2) 相同航速时，初始横倾角从10°增加到20°，线性阻尼系数无明显变化；相同初始横倾角时，线性阻尼系数随着航速的升高而逐渐增加。

3) 横摇消灭曲线的拟合对横摇幅值的敏感性较高，导致由试验结果和CFD方法得到的阻尼系数间的差异会被放大。等效阻尼为线性阻尼和二次阻尼在对应横摇幅值下的组合，可以减弱横摇消灭曲线拟合带来的影响。相同横摇幅值，高航速对应更大的等效阻尼系数。

4) 横摇短期预报结果表明，Fn=0.138和Fn=0.280航速下，CFD方法得到的横摇阻尼系数精度均高于ITH方法，在高海况时表现得尤为明显。