0 引　言

横摇是船舶航行中最危险的运动形式之一，大幅横摇会直接危及航行安全。因此，准确预报横摇阻尼对评估船舶耐波性和保障航行安全具有重要意义。

目前，横摇阻尼预报方法主要包括经验公式法、模型试验法及数值模拟法。其中，经验公式法以Ikeda等^[1]提出的方法为基础，后经Himeno^[2]改进，该方法将横摇阻尼分解为旋涡阻尼、兴波阻尼、摩擦阻尼、舭龙骨阻尼和升力阻尼5个成分，显著提升了定量预报能力。值得注意的是，作为横摇阻尼的重要组成部分，舭龙骨阻尼的计算仍存在较大不确定性。Ikeda等^{[3 − 4]}通过研究平板在非定常运动中的阻力特性，发现阻力系数不仅随运动状态变化，还受到前次摆动幅值的影响，这一现象被称为“流动记忆效应”。总体而言，经验公式法虽计算便捷，但其预报精度有限，难以充分捕捉横摇阻尼的非线性特性及流动记忆效应等复杂物理机制。

在模型试验方面，张晓慧等^[5]基于KVLCC2船模研究了横摇幅值对阻尼的影响，揭示了幅值变化对旋涡特性的显著影响；Patrick等^[6]采用陀螺力矩发生器激励的自航模强制横摇试验，探讨了舵对横摇阻尼的影响及参数横摇运动；章东等^[7]基于横摇衰减数据建立了以角度和角速度为自变量的阻尼模型，并讨论了航速和流动记忆效应对横摇衰减的影响。

在数值模拟方面，Wilson等^[8]最早开发了非定常RANS方法用于预报大振幅横摇运动；BO等^{[9 − 10]}先后对DTMB5512模型开展了横摇衰减和强迫横摇模拟，系统分析了前进速度、横摇幅值及频率对阻尼的影响；Zhou等^[11]则通过多船型算例验证了CFD方法在零航速横摇阻尼计算中的有效性。此外，Gokce等^[12]和Isar等^[13]分别从粘性效应分解和CFD-势流耦合求解2个角度推进了横摇阻尼的精细化预报。

本研究以巴拿马型集装箱船（DTC）为对象，基于非定常雷诺平均纳维-斯托克斯（URANS）方法，模拟了不同激励力矩作用下的横摇过程。通过系统分析横摇运动响应、三维旋涡结构演化以及舭部流场的压力与速度分布，初步揭示了流动记忆效应等因素对横摇阻尼的影响机理。

1 数值计算方法 1.1 控制方程

本文采用非定常雷诺平均纳维-斯托克斯（URANS）方法对船体周围的粘性流场进行数值模拟。控制方程包括质量守恒方程（连续性方程）和动量守恒方程。

质量守恒方程为：

$ \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \vec{u})=0。$

(1)

动量守恒方程为：

$ \frac{\partial (\rho \vec{u})}{\partial t}+\nabla \cdot (\rho \vec{u}\vec{u})=-\nabla p+\nabla \cdot (\overline{\overline{\tau } })+\rho \vec{g}。$

(2)

式中：$ \rho $ 为流体密度；$ t $ 为时间；$ \vec{u} $ 为速度矢量；$ p $ 为静压；$ \rho \vec{g} $ 为重力体积力。$ \overline{\overline{\tau } } $ 为应力张量，其表达式为：

$ \overline{\overline{\tau } }={\mu }_{eff}\left[(\nabla \vec{u}+\nabla {\vec{u}}^{T})-\frac{2}{3}\nabla \cdot \vec{u}\overline{\overline{I} }\right]。$

(3)

式中：$ {\mu }_{eff}=\mu +{\mu }_{t} $ 为等效动力粘性系数，由流体分子粘性系数$ \mu $ 和湍流粘性系数$ {\mu }_{t} $ 共同构成；$ \overline{\overline{I} } $ 为单位张量。

为封闭上述URANS方程，湍流模型选用SST $ k\text{-}\omega $ 模型。该模型能有效地模拟近壁面流动，并精确捕捉在逆压梯度下的流动分离，适用于本研究中复杂的横摇运动与涡脱落过程。该模型通过求解湍动能$ k $ 和比耗散率$ \omega $ 的输运方程来确定湍流粘度$ {\mu }_{t} $ 。其中湍动能$ k $ 方程为：

${ \frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho {u}_{j}k)}{\partial {x}_{j}}={P}_{k}-{\beta }^{*}\rho \omega k+\frac{\partial }{\partial {x}_{j}}\left[(\mu +{\sigma }_{k}{\mu }_{t})\frac{\partial k}{\partial {x}_{j}}\right] 。}$

(4)

比耗散率$ \omega $ 方程为：

$ \begin{split}{\dfrac{\partial (\rho \omega )}{\partial t}+\dfrac{\partial (\rho {u}_{j}\omega )}{\partial {x}_{j}}=}&{\dfrac{\gamma }{{\nu }_{t}}{P}_{k}-\beta \rho {\omega }^{2}+\dfrac{\partial }{\partial {x}_{j}}\left[(\mu +{\sigma }_{\omega }{\mu }_{t})\dfrac{\partial \omega }{\partial {x}_{j}}\right]+}\\ & {2\rho (1-{F}_{1})\dfrac{{\sigma }_{\omega 2}}{\omega }\dfrac{\partial k}{\partial {x}_{j}}\dfrac{\partial \omega }{\partial {x}_{j}}。}\end{split} $

(5)

湍流粘度$ {\mu }_{t} $ 计算式为：

$ {\mu }_{t}=\frac{\rho {a}_{1}k}{\max ({a}_{1}\omega ,\Omega {F}_{2})}。$

(6)

式中：$ {P}_{k} $ 为湍动能产生项，$ {P}_{k}={\mu }_{t}\dfrac{\partial {u}_{i}}{\partial {x}_{j}}\left(\dfrac{\partial {u}_{i}}{\partial {x}_{j}}+\dfrac{\partial {u}_{j}}{\partial {x}_{i}}\right) $ ；$ \Omega $ 为涡量大小；$ {F}_{1} $ 和$ {F}_{2} $ 为混合函数；$ {\beta }^{*}、\beta 、\gamma 、{\sigma }_{k}、{\sigma }_{\omega }、{a}_{1} $ 均为模型常数。

1.2 计算模型和网格划分

本研究以巴拿马型集装箱船标模DTC（Duisburg Test Case）为研究对象。DTC为单桨集装箱船，具备球鼻艏、大外飘艏部及外伸艉部等特征。

该船模配备有方向舵和舭龙骨。方向舵采用NACA0018剖面，带有舵球，并绕轴线扭曲5°；舭龙骨则分五段布置于每舷，位于设计吃水及航速对应的流线上。本研究选用缩尺比为1∶80的船模，其主尺度参数见表1。

计算域设置如图1所示，速度入口边界设定于艏部前1.5倍垂线间长（$ {L}_{\mathrm{{pp}}} $）处，压力出口距艉部3$ {L}_{\mathrm{{pp}}} $ 处，两侧边界均设置为对称面，距离船体3$ {L}_{\mathrm{{pp}}} $ ，顶部和底部边界设置为速度入口，分别距离船体0.5$ {L}_{\mathrm{{pp}}} $ 和1$ {L}_{\mathrm{{pp}}} $ ，本文采用重叠网格技术实现船体六自由度运动，重叠网格区域（move）设计为可绕固定横摇轴旋转的长方体，通过动态插值界面与外部背景区域进行流场数据交换。

计算域网格采用切割体网格生成器进行划分，并引入了表面重构模型以进一步提升网格质量。为精确捕捉关键流动现象，对网格实施了局部加密策略：在舭龙骨、艏部及艉部等易发生流动分离的区域进行了加密，以解析复杂的旋涡结构生成与演化；在船体周围的自由液面区域实施了加密，旨在精确模拟横摇运动引发的兴波现象，网格划分如图2所示。

计算采用VOF（Volume of Fluid）方法追踪自由液面变化，使用分离求解器求解流体控制方程，压力-速度耦合则基于SIMPLE算法求解，空间离散采用二阶迎风格式，时间离散采用二阶隐式格式。

1.3 工况介绍

本文通过对船模施加激励力矩来实现船模的强制横摇运动，激励力矩为标准正弦，其表达式为：

$ {M}_{R}(\varphi )=A\sin (\omega t)。$

(7)

共7种不同激励力矩幅值，具体工况如表2所示。

1.4 数值计算模型验证

为验证本文计算模型的精度与可靠性，将不同工况的稳态横摇幅值与本课题组之前的试验结果^[14]进行对比，如图3所示。可以看出，试验和数值模拟结果吻合良好，整体误差小于3%，验证了本文计算模型的可靠性。

1.5 横摇阻尼计算方法

本文采用Froude能量法^[15]计算横摇阻尼。该方法基于能量守恒原理，通过直接计算一个完整运动周期内阻尼力矩所做的功来获得等效阻尼系数。Froude能量法主要用于分析自由横摇衰减试验数据，其核心思想是：在一个完整横摇周期内阻尼耗散的能量等于系统机械能的损失。在横摇运动峰值位置，系统动能为0，当船舶从峰值角度$ {\varphi }_{i} $ 运动至角度$ {\varphi }_{i+2} $ 时，系统的势能变化量等于其机械能变化量，即：

$ {E}_{P}=\Delta g\int_{{\varphi }_{i}}^{{\varphi }_{i+2}}\overline{GZ}(\varphi ){\mathrm{d}}\varphi。$

(8)

式中：$ \overline{GZ}(\varphi ) $ 为静稳性臂。

阻尼在一个周期内耗散的能量为：

$ {E}_{D}=\int_{{\varphi }_{i}}^{{\varphi }_{i+2}}{B}_{eq}{\mathrm{d}}\varphi。$

(9)

令两者相等，即可得到等效阻尼系数$ {B}_{eq} $ 。

本文中实现强制横摇时对船模施加了外部激励力矩$ M(\varphi ) $ ，根据能量守恒，此时阻尼耗散的能量等于系统机械能变化量与激励力矩做功之和，即：

$ {E}_{D}={E}_{P}+\int_{{\varphi }_{i}}^{{\varphi }_{i+2}}M(\varphi ){\mathrm{d}}\varphi。$

(10)

此时等效阻尼系数$ {B}_{eq} $ 为：

$ {B}_{eq}=\frac{\Delta g\displaystyle\int_{{\varphi }_{i}}^{{\varphi }_{i+2}}\overline{GZ}(\varphi ){\mathrm{d}}\varphi +\displaystyle\int_{{\varphi }_{i}}^{{\varphi }_{i+2}}M(\varphi ){\mathrm{d}}\varphi }{\displaystyle\int_{{t}_{1}}^{{t}_{2}}\overset{\cdot }{\varphi }{\mathrm{d}}t} 。$

(11)

当前周期的等效横摇角$ {\varphi }_{a} $ 定义为2个相邻峰值的平均值：

$ {\varphi }_{a}=\frac{|{\varphi }_{i}+{\varphi }_{i+2}|}{2} 。$

(12)

对于力矩激励的规则横摇运动，当运动达到稳态时，横摇幅值恒定，此时$ {\varphi }_{a} $ 取稳态横摇幅值。

2 数值模拟结果分析

通过对船模施加不同幅值的激励力矩，实现了各工况下的横摇运动。

2.1 横摇运动结果

在相同激励力矩形式下，各工况的横摇角变化趋势相似。此处以 R6 工况为例，展示横摇角的数值模拟结果。图4为 R6 工况的横摇运动时历曲线，力矩激励下的规则横摇运动曲线呈现2个典型阶段：初始的幅值增长阶段和随后的稳态横摇阶段。在初始阶段，正弦激励力矩持续对船模做功，船模吸收能量，导致横摇幅值逐渐增大。进入稳态阶段后，单个周期内激励力矩对船模所做的功与流体阻尼耗散的能量达到平衡，横摇运动幅值趋于稳定。

图5为不同激励力矩幅值对应的稳态横摇角。可以看出，稳态横摇幅值随激励力矩幅值的增加而增大，两者呈正相关关系。更大的激励力矩意味着每个周期内向系统输入更多的能量；为平衡这部分额外输入的能量，系统需要通过增大横摇幅值来增加阻尼耗能，从而建立新的能量平衡状态。

2.2 横摇阻尼变化规律

采用 1.5 节所述的、包含激励力矩项的 Froude 能量法计算各工况的横摇阻尼。此处将横摇运动曲线两个相邻正峰值之间的区间定义为一个完整周期。图6为R7工况等效阻尼系数随时间的变化。可以看出，等效阻尼系数的变化分为递增阶段和稳定阶段，分别对应横摇运动的幅值增长阶段和稳态横摇阶段。

稳态横摇阶段的等效阻尼系数与横摇角的关系表明，阻尼系数随横摇角的增大而增大。这是由于在横摇周期保持不变的情况下，更大的横摇幅值意味着更大的横摇角速度。这导致船体表面的相对速度增大，速度梯度增加，从而使摩擦耗散的能量增多，即摩擦阻尼增大。同时，更大的横摇角会引发更显著的自由表面兴波和更高的舭龙骨运动速度。更强的兴波导致兴波阻尼增加，而更高的舭龙骨速度则会激发更强的旋涡脱落，从而使旋涡阻尼增大。上述摩擦阻尼、兴波阻尼和旋涡阻尼的共同作用，最终导致横摇总阻尼随横摇角增大而增加。

2.3 流场分析 2.3.1 旋涡结构演化

船体所受水动力体现了流场对其的作用。对于横摇运动，船体通过向流场输入能量而损失机械能，这部分能量耗散即表现为横摇阻尼。在横摇流场中，能量主要储存于旋涡结构内。因此，旋涡结构的生成与演化过程，能够反映船体与流场之间能量交换的强弱，即横摇阻尼的大小。横摇运动具有周期性，其旋涡的生成与演化也呈现周期性特征。为清晰揭示此过程，对 R4 工况一个稳态横摇周期如图7所示，5个关键相位$ （t=0,T/8,T/4,3T/8,T/2） $ 的旋涡结构进行可视化分析。船舶横摇运动产生的旋涡主要集中于舭部区域，为获得更清晰的可视化效果，选取图8 所示的局部区域进行分析。

采用 Q 准则（令 $ Q=30 $ ）识别旋涡，得到如图9 所示的旋涡结构演化序列。

$ t=0 $ 时，舭部正自下向上运动，横摇角速度较大，舭龙骨运动速度较快，从而带动周围流体旋转，向流场输入能量。旋涡生成于舭龙骨下方并整体附着于其上，此时船体与流场能量交换强烈，船体横摇能量耗散较快，对应较大的横摇阻尼。

$ t=T/8 $ 时，横摇角速度减小，舭龙骨运动放缓。旋涡主体从舭龙骨表面脱离并向下游发展，附着区域收缩。这表明能量输入速率下降，船体与流场的能量交换减弱，横摇阻尼相应减小。同时，旋涡形态开始拉伸变形。

$ t=T/4 $ 时，横摇角达到极值，角速度为0。此时，起始涡已完全脱离舭龙骨，演化为一个清晰的、环状的涡结构，悬浮于舭部流场中。船体停止向该涡输入能量，旋涡进入自由的粘性耗散阶段，其携带的能量将逐步转化为流体内能。

$ t=3T/8 $ 时，横摇运动方向已反转，舭龙骨开始向下运动。一方面，上一周期残留的旧涡尚未完全耗散，其结构发生破碎与拉伸；另一方面，新的旋涡开始在舭龙骨上缘生成，但强度较弱。随着横摇角速度的增大，船体向流场输入能量的速率加快，新生旋涡强度迅速增长，并逐渐取代旧涡，至$ t=T/2 $ 时成为流场中主要的能量载体，其形态与$ t=0 $ 时刻的起始涡镜像相似，标志着一个完整演化周期的结束和下一个周期的开始。

值得注意的是，在横摇运动方向改变后的一段时间内，旧涡仍然存在于船体附近。旧涡的存在旧涡诱导的流场会改变新涡的生成环境和演化过程，这意味着当前周期的横摇阻尼不仅取决于当前的运动状态，也受到上一周期横摇运动的影响。

2.3.2 压力场分析

为探究旋涡结构与受力间的直接关联，在图8 所示舭龙骨中部（约 1/2 长度处）取一横截面进行分析，船首方向设定为垂直于纸面向外。采用随体坐标系，可视化该截面的压力场云图，此处的压力$ p $ 为$ p={p}_{{\mathrm{total}}}-{p}_{{\mathrm{static}}} $，$ {p}_{{\mathrm{total}}} $为总压力，$ {p}_{{\mathrm{static}}} $为静水压力，结果如图10 所示。

$ t=0 $ 时，舭龙骨后方出现一个范围集中、强度显著的极低压区，此低压区导致船体左右舷侧压力分布严重不对称，如图11所示，产生一个与横摇方向相反的压力力矩，即阻尼力矩的主要成分。低压区的强度和范围直接决定了瞬时阻尼力矩的大小。$ t=T/8 $ 时，随着旋涡脱离和发展，低压区也随之从船体表面向下游移动，强度有所减弱。压力场的变化与旋涡的演化同步，阻尼力矩相应减小。

低压区位置与2.3.1节分析的旋涡生成位置一致，表明低压区是由旋涡诱导产生的。从流体微团的受力平衡角度分析：作曲线运动的流体微团受离心力作用，有向外运动的趋势；为维持其旋转运动，流场内部需形成由外指向内的压力梯度力，以此作为向心力来平衡离心力。因此，压力沿旋转半径向外递增，旋转中心成为压力最低的区域。旋涡强度越大，则离心力越大，所需的向心力（即压力梯度）也越大，因此低压区的强度与旋涡强度呈正相关。

值得注意的是，当横摇运动方向改变后($ t=T/4 $ 时刻)，尽管横摇角速度为0，但由旧涡诱导的低压区依然清晰地存在于流场中。这意味着，阻尼力矩并不与角速度瞬时同步，而是存在相位滞后。船体在运动停止的瞬间，仍在承受来自上一周期旋涡的压力作用，这是流动记忆效应在受力层面的直接体现。$ t=3T/8 $ 时，流场压力分布复杂，旧涡的残余低压区与新涡诱发的低压区同时存在，导致低压区范围明显较大。这表明当前横摇周期内船体所受的水动力与上一周期的运动有关，上一周期的横摇幅值越大，则残留的旧涡强度就越强，对当前压力场的影响也越显著，进而对横摇阻尼产生更明显的影响。

2.3.3 速度场和涡量场分析

图12 为所述截面的速度云图与涡量云图。对比可知，$ t=0 $ 时刻的旋涡强度明显大于后续时刻，但此时船体附近流体速度的幅值却显著较小。理论上，旋涡强度越大，其诱导的速度场也越大。然而在$ t=0 $ 时刻，旋涡结构紧贴船体壁面，其诱导的速度场受到壁面无滑移边界条件的强烈抑制，导致近壁区域的速度幅值衰减。

此外，$ t=0 $ 时刻速度云图中显示有速度值的区域（速度大于0）范围较大。这是因为旋涡靠近壁面时，受船体表面阻挡，被其诱导的流体无法穿过壁面，只能沿船体表面方向扩散，从而扩大了受影响流体的区域。

$ t=T/4 $ 时刻，与$ t=0 $ 时刻形成鲜明对比。此时，尽管涡量强度因耗散而减弱，但在稍远离船体的区域（涡核与壁面之间），速度幅值却达到一个较高的水平。这是因为脱体的涡环其诱导的速度场不再受壁面直接约束，可以更完整的展现出来。$ t=3T/8 $ 时刻，虽然新生涡刚刚生成，但由旧涡诱导的背景速度场仍然清晰可见，舭部流体保持着较大的速度。这意味着，当船体开始新一轮运动时，它并非在静止流体中启动，而是在一个已有一定速度的流场中运动。这一预先存在的背景速度场，势必会改变船体与流体之间的相对速度，从而影响新旋涡的生成模式和强度，最终影响能量交换效率与横摇阻尼。这从流场演变的角度，进一步证实了上一次横摇运动对当前横摇运动的影响。

上述机理或许可以解释幅值增长阶段与稳态横摇阶段阻尼系数的差异。在稳态横摇时，前后周期的运动幅值相当，船体周围的流场结构及速度分布周期性重复，变化规律一致。而在幅值增长阶段，每一次横摇运动都起始于一个流体速度相对较低、流场相对平静的状态，然后过渡到速度更大、流场更紊乱的状态。这种从低速背景流场启动并建立强旋涡的过程，需要向流场输入更多的能量，意味着阻尼耗散更大，因此表现为更大的等效阻尼系数。

3 结　语

本文对正弦激励力矩作用下的船舶横摇运动进行了数值模拟研究。基于 STAR-CCM+ 软件建立了数值计算模型，通过与试验结果对比验证了模型的可靠性；分析了不同工况下的横摇阻尼与稳态横摇角，并对稳态横摇周期内不同时刻的三维旋涡结构演化、舭部截面压力场、速度场与涡量场进行了讨论。主要得到以下结论：

1）力矩激励下的强制横摇运动可分为幅值增长与稳态横摇2个阶段，分别对应等效阻尼系数的增加阶段与稳定阶段。稳态横摇阶段的结果显示，横摇角越大，对应的等效阻尼系数也越大。

2）三维旋涡结构演化分析表明，舭部旋涡的形态与横摇角速度密切相关：横摇角速度越大，旋涡在舭龙骨上的附着区域越大，船体与流场之间的能量交换越强，能量耗散速率越快。

3）船体周围压力分布与旋涡结构位置显著相关。旋涡会导致船体周围出现低压区，从而影响船体受力与横摇阻尼；当旋涡紧贴船体表面时，受壁面影响，其诱导速度场影响范围较大，但速度幅值相对较小。

4）流动记忆效应通过残留旋涡与背景流场对当前横摇运动产生持续影响。数值模拟揭示了明显的流动记忆效应，其影响主要体现在2个方面：一是压力影响，上一周期残留的旋涡通过其诱导的低压区持续作用于船体表面，改变压力分布，从而影响当前横摇运动的阻尼力矩；二是速度场影响，旧涡诱导的背景速度场为新周期舭龙骨运动提供了初始流动条件，影响新旋涡的生成强度与形态，进而改变船体与流体之间的能量交换效率。