稳定性是船舶设计和航行过程中是重要因素，而横摇运动影响船舶稳定性，所以减少横摇运动是船舶安全航行的关键。在船体上安装舭龙骨是常用的减摇方法，它可以增加船体的横摇阻尼。Ikeda等^[1]基于试验将横摇阻尼分成摩擦阻尼、兴波阻尼、旋涡阻尼、升力阻尼和舭龙骨阻尼。为了设计合理的舭龙骨，研究带有舭龙骨的船体横摇衰减模拟具有重要的意义。由于粘性对横摇运动具有重要的影响，横摇运动相比于其他自由度运动具有强烈的非线性，横摇运动的模拟比其他自由度运动更加困难。为了计算船舶的横摇运动，Ikeda等^[1]基于试验结果得到半经验公式，但是精度不高。传统的势流理论方法不考虑粘性，仅仅可以计算出兴波阻尼，但是无法计算关于粘性的阻尼。如果船体加上舭龙骨，它会产生更强的粘性效应，传统的势流理论方法计算会更加困难。随着计算机技术的发展，计算流体力学方法已经广泛地应用到船舶横摇的研究中。Wilson等^[2]利用RANS方法研究DTMB5512船型在不同航速下的横摇衰减运动，并且展示了横摇运动产生的流场信息。Yang等^[3]研究了船舶航行中的船舶横摇的问题，重点对强非线性进行了分析。罗天等^[4]利用自主开发的naoe-FOAM-SJTU求解器研究了S60与DTMB5415船模的自由与强迫横摇运动，发现旋涡阻尼所占的比例最大。Lavrov等^[5]研究了二维剖面强迫横摇运动，重点分析了横摇中泄涡的现象。对于带有舭龙骨的横摇模拟，Xie等^[6]基于FLUENT计算了带舭龙骨的柱体的横摇阻尼，并且得到了横摇过程中船体周围的压力分布图像，但没有对周围涡量场进行分析。蒋银等^[7]结合重叠网格方法对三维带有舭龙骨船体的横摇运动进行了数值模拟研究，重点研究了阻尼系数的频率相关性。Irkal等^[8]采用Flow3D软件研究了带有舭龙骨船体剖面的横摇阻尼运动，结合自己试验数据进行验证分析，但在数值模拟的后期出现了局部数值震荡。

本文基于OpenFOAM^[9]中InterDyMFoam求解器，建立粘性流体数值模型，考虑计算域边界辐射波反射的影响，对带有舭龙骨的二维船体横截面的横摇运动进行数值模拟。利用流体体积分数方法(volume of fluid, VOF)方法捕捉自由面，利用六自由度运动模型和动网格技术计算船舶运动，由于本文只研究横摇运动，所以只考虑其中一个自由度。采用湍流模型来模拟横摇中的湍流，对不同高度舭龙骨下的横摇展开模拟，分析舭龙骨附近涡量场。

1 两相流数值模型 1.1 控制方程

采用两相流模型求解水和空气运动，利用VOF方法捕捉自由液面，假设水和空气不可压缩。连续性方程，动量方程，相方程表示为^[10]:

连续性方程:

$ \nabla \cdot \mathit{\boldsymbol{U}} = 0 $

(1)

动量方程:

$ \begin{array}{*{20}{c}} {\frac{{\partial \rho \mathit{\boldsymbol{U}}}}{{\partial t}} + \nabla \cdot \left( {\rho \left( {\mathit{\boldsymbol{U}} - {\mathit{\boldsymbol{U}}_{\rm{g}}}} \right)\mathit{\boldsymbol{U}}} \right) - \nabla \cdot \left( {{\mu _{{\rm{eff}}}}\left( {\nabla \mathit{\boldsymbol{U}} + \nabla {\mathit{\boldsymbol{U}}^{\rm{T}}}} \right)} \right) = }\\ { - \nabla {p_{\rm{d}}} - \mathit{\boldsymbol{g}} \cdot \mathit{\boldsymbol{X}}\nabla \rho + \sigma \kappa \nabla \gamma + {f_{\rm{s}}}} \end{array} $

(2)

相方程:

$ \frac{{\partial \gamma }}{{\partial t}} + \nabla \cdot \left( {\mathit{\boldsymbol{U}}\gamma } \right) + \nabla \cdot \left( {{U_{\rm{r}}}\gamma \left( {1 - \gamma } \right)} \right) = 0 $

(3)

式中:ρ表示流体密度；U表示速度矢量；U_g表示网格速度；μ_eff是考虑湍流的有效动力粘度系数，表达式为μ_eff=ρ(ν+ν_t)，其中ν为流体运动粘度系数，ν_t为湍流粘度系数；p_d表示相对动压力，表达式为p_d=p-ρgX，X为相对于参考高度的位置矢量；σ为流体表面张力系数，本文取0.728 kg/s²；κ表示自由表面曲率，它可以计算表面张力，计算公式为$\kappa = - \nabla \cdot \left( {\frac{{\nabla \gamma }}{{\left| {\nabla \gamma } \right|}}} \right) $，其中“·”表示求散度；γ为流体体积分数，由VOF相方程计算得出，如果γ=0，则表示网格内是空气，如果0＜γ＜1，则表示网格内是流体交界面，如果γ=1，则表示网格内是水；f_s为阻尼项，用于消除横摇运动在远处产生的辐射波。

相方程中，U_r为用于界面压缩的速度场，只在自由面处起作用，压缩速度方向垂直于自由面，具体的表达式为$ {\mathit{\boldsymbol{U}}_r} = {C_\alpha }\left| \mathit{\boldsymbol{U}} \right|\frac{{\nabla \alpha }}{{\left| {\nabla \alpha } \right|}}$，其中$\frac{{\nabla \alpha }}{{\left| {\nabla \alpha } \right|}} $表示自由面的法向方向，|U |为压缩速度的大小，表达式为$\left| \mathit{\boldsymbol{U}} \right| = \min \left\{ {{C_a}\frac{{\left| \varphi \right|}}{{\left| {{\mathit{\boldsymbol{S}}_{\rm{f}}}} \right|}}, \max \left( {\frac{{\left| \varphi \right|}}{{\left| {{\mathit{\boldsymbol{S}}_{\rm{f}}}} \right|}}} \right)} \right\} $，表达式中φ表示速度通量；C_a为压缩因子，本文中取值1；S_f表示网格表面的面积矢量。

1.2 自由表面阻尼消波方法

本文采用二维剖面进行横摇数值模拟，计算模型示意图如图 1。

由于CFD模拟是在有限域中模拟，所以为了消除两边远处辐射波的反射所造成的数值误差，采用自由表面阻尼消波的方法，在动量方程里面添加阻尼项达到消波的目的:

$ {f_s} = \rho \theta \mathit{\boldsymbol{U}} + \rho \theta '\mathit{\boldsymbol{U}} $

(4)

式中:θ与θ'分别用于右边和左边的消波系数。

$ \begin{array}{*{20}{c}} {\theta = \left\{ \begin{array}{l} \frac{{x - {x_0}}}{{{x_1} - {x_0}}}{\theta _1},\;\;\;\;{x_0} \le x < {x_1}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le x < {x_0} \end{array} \right.}\\ {\theta ' = \left\{ \begin{array}{l} \frac{{x - {{x'}_0}}}{{{{x'}_1} - {{x'}_0}}}{\theta _1},\;\;\;\;{{x'}_1} \le x < {{x'}_0}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{x'}_0} \le x < 0 \end{array} \right.} \end{array} $

(5)

式中:θ₁为无因次阻尼系数，本文取θ₁为5；x₀、x₁、x'₀、x'₁用于确定消波长度，消波长度约为0.5倍波长，本文取1.0 m。

1.3 六自由度运动方程与动网格

为了计算横摇运动，通过求解六自由度运动方程来得到横摇的运动。

通过求解控制方程，可以得到船体周围流场的速度和压力，船体剖面受到的力F和力矩M可以通过边界表面积分得到:

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{F}} = m\mathit{\boldsymbol{g}} + \int_S {{\rm{d}}\mathit{\boldsymbol{\tau }}} \\ \mathit{\boldsymbol{M}} = r \times m\mathit{\boldsymbol{g}} + \int_S {\mathit{\boldsymbol{r}} \times {\rm{d}}\mathit{\boldsymbol{\tau }}} \end{array} \right. $

(6)

式中:m为船体质量，r为力臂，τ为应力张量。

得到船体表面作用力之后，代入船体运动方程可得船体运动:

$ \left\{ \begin{array}{l} m\mathit{\boldsymbol{\dot u}} = X{I_x}\mathit{\boldsymbol{\dot p}} = \mathit{\boldsymbol{K}}\\ m\mathit{\boldsymbol{\dot v}} = Y{I_y}\mathit{\boldsymbol{\dot q}} = \mathit{\boldsymbol{M}}\\ m\mathit{\boldsymbol{\dot w}} = Z{I_z}\dot r = \mathit{\boldsymbol{N}} \end{array} \right. $

(7)

式中: ${\mathit{\boldsymbol{\dot u}}} $、$ {\mathit{\boldsymbol{\dot v}}}$、${\mathit{\boldsymbol{\dot w}}} $分别为3个坐标轴方向的加速度，${\mathit{\boldsymbol{\dot p}}} $、$ {\mathit{\boldsymbol{\dot q}}}$、${\mathit{\boldsymbol{\dot r}}} $分别为3个坐标轴方向的角加速度。X、Y、Z和K、M、N分别为力和力矩的3个分量。I_x、I_y、I_z为船体绕质心的惯性矩。

求得船体加速度后，对于时间进行积分可以得到时间段内船体的位移，进而得到船体下一个时刻的位置。时间积分步进方法是通过Newmark算法实现。

Newmark算法^[11]的基本原理如下，下个时间步的加速度和速度为

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot u}}}_{t + 1}} = {{\mathit{\boldsymbol{\dot u}}}_i} + \left( {1 - \delta } \right){{\mathit{\boldsymbol{\ddot u}}}_i}\Delta t + \delta {{\mathit{\boldsymbol{\ddot u}}}_{i + 1}}\Delta t\\ {\mathit{\boldsymbol{u}}_{i + 1}} = {\mathit{\boldsymbol{u}}_i} + {{\mathit{\boldsymbol{\dot u}}}_i}\Delta t + \left[ {\left( {\frac{1}{2} - \beta } \right){{\mathit{\boldsymbol{\ddot u}}}_i} + \beta {{\mathit{\boldsymbol{\ddot u}}}_{i + 1}}} \right]\Delta {t^2} \end{array} $

(8)

式中δ和β是与精度和稳定性有关的参数，本文分别取0.5与0.25。

为了实现船体边界的运动，计算过程中采用动网格技术，本文采用的是类似于“随体网格”的方法，如图 2所示。变形网格层用来确定网格变形的区域，在船体附近的那层网格同船体一起运动，运动方向和大小和船体一致，而远离船体的最外层网格固定不动，之间的网格按照线性插值决定网格变形的大小。

值得注意的是要选取适当的距离避免求解过程网格过于形变造成计算发散，近距离值选用距离壁面相对近的位置，大概4~5层近壁面网格左右，而远距离值可以适当选用较大的距离。本文模拟计算选用近距离值为0.005 m，而远距离值选用0.4 m。

1.4 湍流模型

船体横摇过程中流场处于湍流状态，本文采用二方程的湍流模型进行处理，采用低雷诺数湍流模型k-ω模型进行计算。

k-ω双方程模型引入两个附加的传输方程，一个为湍流动能k，另一个为湍流耗散率ω。该模型由Wilcox等^[12]提出，其中k与ω分别为:

$ k = \frac{{\overline {{{u'}_i}{{u'}_i}} }}{2} $

(9)

$ \omega = \frac{\mu }{\rho }\overline {\left( {\frac{{\partial {{u'}_i}}}{{\partial {x_j}}}} \right)\left( {\frac{{\partial {{u'}_i}}}{{\partial {x_j}}}} \right)} /{C_\mu }k $

(10)

湍流动能和湍流耗散率的约束方程为:

$ \left\{ \begin{array}{l} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _{\rm{k}}}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{P_k} - {\beta ^ * }\rho k\omega \\ \frac{{\partial \left( {\rho \omega } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \omega {u_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _{\rm{ \mathsf{ ω} }}}}}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha \frac{\omega }{k}{P_{\rm{k}}} - \beta \rho {\omega ^2} \end{array} \right. $

(11)

P_k是湍动能生成项，定义为:

$ {P_{\rm{k}}} = \frac{{{\mu _{\rm{t}}}}}{\rho }\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)\frac{{\partial {u_i}}}{{\partial {x_j}}} $

(12)

式中μ_t为湍流粘度，计算公式为μ_t=ρk/ω。该模型的常数α=0.52, β=0.072, β^*=0.09, σ_k=0.5, σ_ω=0.5。

在运用湍流模型的时候要处理近壁面内的流动，近壁面处边界层内的流动可以分为粘性层，缓冲层和对数律层。网格高度的选取根据离壁面的无量纲距离y⁺确定。粘性层满足y⁺≤5。y⁺的定义为

$ {y^ + } = \frac{{{\mu _{\rm{ \mathsf{ τ} }}}y}}{\nu } $

(13)

式中:y为第一层网格节点高度; u_τ为摩擦速度；ν为流体运动粘度系数。

本文采用壁面函数法，由于模型尺度较小，将第一层网格放置在粘性层，靠近壁面最小的网格在1 mm，这样可以应用于低雷诺数湍流模型。

1.5 边界条件和初始条件

由于本文采用二维边界，所以前后两个边界定义为empty边界，表示模型是二维。除此之外，其他的边界有左右两边的出流边界outflow1与outflow2，上面的大气边界atmosphere，下面的底部边界bottom，还有二维船体剖面边界hull。

关于数值边界条件，在OpenFOAM中是通过设置各个物理场的边界来实现。

对于速度场，在两边出流边界和大气边界定义为自由出流边界条件，在底部速度为0，在船体表面为物面速度。对于压力场，在出流和底部边界定义为零梯度边界，在大气边界处压力为一个大气压，在船体表面为针对压力的第二边界条件。因为采用VOF方法捕捉自由面，所以要设置相场边界条件，在边界处设置为零梯度边界。另外关于湍流场的边界，除了在壁面处如船体和底部处采用壁面函数边界，其他边界设置为零梯度边界。

对于初始条件，本文计算自由横摇运动，初始时刻将船体绕质心旋转一定角度，船体初始速度为0，通过setField工具设置相场来定义自由面，水的体积分数为1，空气的体积分数为0。

1.6 数值方法

OpenFOAM将所有物理量定义在网格中心，采用有限体积法离散N-S方程。时间项采用一阶精度的隐式欧拉格式积分梯度项和拉普拉斯项采用二阶精度的中心差分格式。速度散度项采用二阶精度的线性迎风格式，体积分数散度项采用vanLeer格式，湍流动能和能量耗散的散度项采用一阶精度迎风格式。

对于瞬态问题，N-S方程的求解采用基于压力的隐式算子分裂算法(pressure correction method on unstructured grids, PISO)^[9]，在OpenFOAM中是通过pimple算法进行设置。在PISO算法中，每个时间步内首先通过求解动量方程进行速度预测，再通过求解压力方程进行速度矫正，称为内迭代，内迭代次数一般是2到3次，本文取2次，然后增加一个时间步，再进行相同的迭代，如此循环计算下去。

2 船体横摇计算与分析 2.1 算例设置

为了验证本文所建立的数值模型，将选取Irkal等^[8]采用的二维船体横摇试验展开数值模拟。取一艘10×10⁴ t巴拿马船型中的横剖面进行计算，缩尺比为1 :100。模型型宽0.3 m，型深为0.2 m，吃水为0.12 m，船体的重心距离基线0.08 m。本文采用舭龙骨高度分别为10、15和20 mm进行计算，宽度为3 mm。船体质量为20.88 kg，惯性矩为0.225 8 kg ·m²。算例的参数设置见表 1。

计算数值模型整体设置如图 3所示，x、y、z方向计算区域分别为12、0.58和1.01 m; 自由面位于距离底面0.6 m。尽管y方向有一定厚度，但这个方向只有一层网格，所以本文数值模型依然是二维模型，厚度取值为0.58 m是为了和试验厚度保持一致。为了更精确地计算船体附近的流动，对船体附近的网格采用局部加密。网格划分是利用OpenFOAM中snappyHexmesh工具完成，它可以较好地对复杂曲面进行网格划分。本文算例最大网格为16 mm，对船体附近和自由面进行加密，最小网格为1 mm。边界层最小厚度为0.35 mm，延展比为1.2。船体周围网格层次划分见图 4，舭龙骨附近网格划分见图 5。计算过程中，采用固定计算步，步长定为0.000 2 s。

2.2 网格收敛性分析

为了进行网格收敛性分析，采用BK1010算例，利用三套网格进行计算。最密网格船体周围网格尺寸在1 mm，网格数量约为36万，然后将尺寸按照$ \sqrt 2 $倍扩大两次形成粗网格，网格数依次约为18万和9万，计算结果见图 6。

从图中可以看出，不同网格下计算的横摇自由衰减曲线非常接近。说明本文采用的网格划分满足计算要求。为了减少计算时间，又考虑到计算精度，后面的算例将采用18万网格进行计算。

2.3 不同舭龙骨的计算结果

图 7为BK1020、BK1520和BK2020算例中数值模拟得到的横摇运动时间历程曲线和试验结果及数值计算结果的比较。

通过对比可以发现，本文数值模型结果和试验结果吻合良好，数值模型可以有效预报不同高度舭龙骨的减摇效果。此外，在某些细节处优于Mohsin等的模拟结果，尤其是横摇最后的衰减阶段。本文数值模型有效地模拟出了横摇运动10 s后的衰减趋势。其主要原因在于模型中自由表面阻尼层的建立，阻尼层可以有效吸收船体运动引发的辐射波，避免其在模型边界处反射，避免反射波造成的误差。因此，在横摇模拟中进行消波处理是十分重要的。

图 8将不同舭龙骨的横摇自由衰减曲线进行比较，可以看出舭龙骨高度越高横摇效果越好，所以工程中为了达到更好的减摇效果，可以增加舭龙骨高度，但受到舭龙骨根部应力的非线性增大，对结构安全性会造成影响，因此需要综合考虑。

2.4 阻尼分析

为了定量地分析舭龙骨减摇的效果，对阻尼系数进行求解。阻尼分析有线性和非线性模型两种，本文采用二次非线性模型。自由衰减模型公式为

$ \left( {I + A} \right)\ddot \phi + {B_1}\dot \phi + {B_2}\dot \phi \left| {\dot \phi } \right| + C\phi = 0 $

(14)

式中ϕ是横摇角，I是转动惯量，A是转动附加质量，B₁为线性阻尼系数，B₂为非线性阻尼系数，C为转动恢复力矩系数。

通过横摇衰减曲线可以得到两种阻尼系数，提取横摇曲线局部极值ϕ_k(k=1, 2, …)，可以得到:

$ {\phi _{{\rm{m}},{\rm{k}}}} = \frac{{\left( {{\phi _k} + {\phi _{k + 1}}} \right)}}{2} $

(15)

$ {\delta _{\rm{k}}} = \ln \left( {\frac{{{\phi _k}}}{{{\phi _{k + 1}}}}} \right) $

(16)

$ {\zeta _{{\rm{t}},{\rm{k}}}} = \frac{1}{{\sqrt {1 + {{\left( {2{\rm{ \mathsf{ π} }}/\delta } \right)}^2}} }} $

(17)

$ {B_{{\rm{t}},{\rm{k}}}} = 2\left( {I + A} \right){\omega _{\rm{n}}}{\zeta _{{\rm{t}},{\rm{k}}}} $

(18)

通过$ \frac{8}{{3{\rm{ \mathsf{ π} }}}}{\omega _n}{\phi _{m, {\rm{k}}}}$和B_{t, k}合并得到拟合曲线可以得到等效阻尼为:

$ {B_t} = {B_1} + {B_2}\frac{8}{{3{\rm{ \mathsf{ π} }}}}{\omega _{\rm{n}}}{\phi _{\rm{m}}} $

(19)

式中ω_n为横摇的自然频率。通过曲线的截距和斜率可以得到阻尼系数B₁和B₂。图 9为从衰减曲线得到的$\frac{8}{{3{\rm{ \mathsf{ π} }}}}{\omega _\rm{n}}{\phi _{\rm{m}, {\rm{k}}}}$与B_{t, k}的离散点和最小二乘法拟合的曲线，从曲线方程可以得到二次非线性阻尼系数B₂。

表 2为拟合得到的二次非线性阻尼系数与试验结果的误差比较，试验结果来自于文献[8]做的船体横剖面的实验。可以看出误差满足工程要求，说明了本文模型可以准确地预报横摇非线性阻尼。

2.5 舭龙骨附近流场分析

为了研究舭龙骨减摇的原理，将BK1025算例附近的流场与试验PIV(高速摄影技术)得到的流场进行比较。图 10为不同时间得到的涡量场比较图。涡量范围为-75~75 s^-1。

从图中可以看出数值模拟中捕捉到的漩涡运动和试验测量结果吻合良好，尤其是前三个时刻涡量大小和位置基本与试验均很好地吻合。在最后一个时刻，涡量位置预报基本准确，但数值预报涡量值小于试验实测值，这可能是数值扩散影响了涡量的持续存在。

此外，通过涡量场的变化可以发现，舭龙骨的存在很大程度上改变了水下的流场，在船体舭部造成了大量的泻涡，将船舶横摇运动的动能持续地转化为流体运动的动能，以一种阻尼的形式消耗了船舶能量，进而达到减摇效果。

3 结论

1) 本文提出的数值模型可以有效预报带舭龙骨船体的横摇自由衰减运动，数值预报结果与试验测试结果高度吻合，并准确预报了不同高度的舭龙骨的减摇效果。

2) 此类船舶水动力数值模型需要加入远场自由表面消波处理，阻尼层可以吸收辐射波，抑制辐射波的反射，最终显著改善自由衰减运动预报的精度。

3) 数值模拟发现，船舶在发生横摇运动过程中，舭龙骨尖端会不断产生泄涡，将船体的动能转化为流体的动能，有效耗散船体横摇能量，达到减摇效果。此现象与试验观测结果吻合良好，证实此模型可以有效预报舭龙骨涡旋的产生及运动。