波浪上船体与液舱晃荡耦合的非线性时域模拟

李超凡， 朱仁传， 周文俊

(1.上海交通大学 海洋工程国家重点实验室， 上海 200240； 2.中国船舶及海洋工程设计研究院， 上海 200011)

船舶行业的发展和研究与工程需要密切相关。如今LNG、LPG和一些特种液货船舶等的使用比例日益增大，但该类船舶在航行状态时剧烈的舱内流体晃荡却对营运安全性提出考验，尤其当船舶没有装满时，船舶稳性发生改变，舱内液体受船体运动影响产生剧烈晃荡，由此产生的晃荡力与船体外部所受的波浪力相互耦合，给运动预报带来困难。此外，对某些特种载液船型来说，人们有时需要知道其在波浪中工作时某一时刻以后若干秒内船体的运动幅值、位置、姿态等，对此展开更加准确有效的运动时历计算研究十分必要。Vassalos[1]、Zaraphonitiset[2]等较早地基于势流理论模拟破损舱进水问题，但其忽略内部液舱晃荡的非线性影响。Santos等[3-4]研究了随机波浪中破损船舶时域运动的数值方法，Newman[5]开发了WAMIT用于计算线性化的波浪中液舱晃荡与船体耦合运动。Kim等[6]采用脉冲响应函数法求解线性的船舶时域运动，对非线性液舱晃荡采用有限差分法进行模拟。Tang等[7]计算了时延函数对运动的影响。Spanos和Papanikolaou[8]，Liu和Papanikolaou[9]应用频域零航速自由面格林函数考虑航速修正，在瞬时湿表面上进行压力积分获得非线性FK力和回复力进行模拟。Ahmed、Hudson等[10]也采用了类似的方法。Zhao[11-12]和Hu等[13]基于势流理论处理内部晃荡问题，证实结果的可行性。Zou[14]等用时域势流求解船体运动，商用CFD软件处理内部晃荡。Huang[15]提出了一种基于线性理论的能量耗散条件，在非线性模型中考虑了晃荡流的流动粘性效应。洪亮等[16]对船体内外流场均采用时域势流理论求解，将船体运动与液舱流体晃荡构建耦合运动方程。唐恺等[17]结合辐射和绕射问题的线性解建立船舶非线性时域预报方法。李裕龙等[18-19]基于全时域势流理论研究有航速时的液舱晃荡耦合问题。

本文基于势流理论提出了一种计算在波浪中船体与液舱晃荡耦合运动的非线性时域混杂法，重点考虑外流场非线性对液舱晃荡耦合运动影响。文中内部液舱晃荡采用时域Rankine元法求解，将舱内液体流动引起的力矩进行修正；对时域运动方程，辐射力与绕射力采用完整船舶计算的结果，对入射波力和静恢复力考虑瞬时湿表面的非线性影响；将液舱晃荡力矩和非线性入射波力、静恢复力与绕射波浪激励力放在时域运动方程的右端建立耦合运动方程，实现对规则波或不同工况下载液船舶运动时历的求解。本文中有航速的工况采用移动脉动源方法计算频域结果，计算分析对象为一载液S175船，对其进行有无液舱的线性时域模拟与非线性时域模拟，给出运动幅值响应算子RAO曲线，并与实验值进行对比验证。

1 载液船舶运动的基本理论

1.1 参考坐标系的建立

大地坐标系为o0x0y0z0，随船舶一起以相同速度运动而构成表征船舶摇荡位移和姿态的基准的坐标系为参考坐标系oxyz，该坐标系的oxy平面与静水面重合。与船体固结，随船体一起摇荡的坐标系为动坐标系o′x′y′z′。各坐标系如图1所示。

图1 船体运动坐标系Fig.1 The coordinate system for ship motion

在液舱内部由于考虑流体晃荡问题，对其采用三维时域势流理论求解，为此需提前确定时域势流计算中的边界条件。若在大地坐标系或参考坐标系中计算，则自由面条件的表达较为简单，但舱壁面边界条件十分复杂。因此，本文中液舱先在oxy面与舱内平均液面重合的动坐标系中求解，此时舱壁上速度势的法向偏导数均为0，使外界激励的影响体现在网格数更少的自由面的条件中，将液舱晃荡在动坐标系中产生的力矩模拟计算后，再通过运动和力在坐标系中的转换实现船体运动与液舱晃荡的耦合计算。液舱计算所在坐标系如图2所示。

图2 液舱内液体晃荡坐标系示意Fig.2 The three ordinates describing sloshing in the liquid tank

1.2 基于脉冲响应函数法的频-时域转换

本文辐射力部分采用频域转时域方式求得[20]。在时域计算时采用卡明斯的脉冲响应法思想，把其时间历程看成一系列瞬时的脉冲运动组成，由此建立的船体运动微分方程可表达为：

(1)

频域的水动力系数中附加质量μij和阻尼系数λij在时域的转换关系为：

(2)

(3)

则时延函数的表达式为：

(4)

1.3 非线性Froude-Krylov力和静恢复力

在时域理论角度，其相对频域的优势还体现在能够反映瞬变或非线性结果的影响，由于船舶存在内部液舱晃荡问题，在计算时考虑一定的非线性影响十分必要。本文在计算时通过瞬时湿表面上的流体压力积分得到非线性入射波力和静回复力，波幅为ζa的简谐规则波的入射势为：

i(ωt-k(xcosβ+ysinβ))]

(5)

非线性回复力根据其物理含义直接计算，船舶静水回复力等于实际湿表面的流体静水力与船舶正浮平均湿表面流体静力之差：

(6)

(7)

以上计算需求取瞬时湿表面，为此需要先对船体网格进行坐标转换。由于运动是动坐标系相对参考坐标的，而船体网格在动坐标系下表示，所以先将网格点坐标位置转移到参考坐标系下便于用波面来截船体表面。图3为瞬时湿表面船体网格示意图。

图3 瞬时湿表面船体网格Fig.3 Panels for transient wet surface

坐标转换的方式为：

(8)

(9)

1.4 边界元法求解时域液舱晃荡

考虑液舱内部非定常的流体晃荡问题，本文采用基于三维势流理论的边界元法对此进行时域数值分析。为减少计算的复杂度，液舱内部在动坐标系下求解，此时舱壁上的速度势法向导数均为0，之后通过坐标转换实现与船体运动的耦合。场内速度势满足的定解条件：

(10)

根据液舱尺寸和计算时装载深度进行网格划分实现模型的空间离散，场内任意一点的速度势可根据格林定理得到，选取三维简单格林函数，则速度势可表达为：

(11)

时间离散方式采用中心差分格式，计算表达，

(12)

(13)

设已知上一时刻的自由面速度势与波面升高，根据此刻的速度势φ和波面升高ζ，再经自由面条件计算此刻每个离散点的∂φ/∂t和∂ζ/∂t，实现时间步的递进。

1.5 基于时域混杂法的船舶时域运动方程

考虑弱散射的时域混杂法船舶时域运动的方程为：

(14)

图4 非线性时域耦合运动计算流程Fig.4 Flowchart of non-linear prediction for coupled motion in time domain

2 载液船舶与计算模型

本文选取一加载方形液舱的S175高速集装箱船作为算例模型[14]，船体横剖线图如图5所示。该模型试验由邹康等在中国船舶科学研究中心耐波性水池中完成，实船与模型的主尺度对比如表1所示，缩尺比为1∶55。

图5 S175船型线图Fig.5 Body plan of S175

表1 S175实船与模型主尺度Table 1 Principal particulars of S175 and its model

液舱位于第9站至第13站之间，关于船舶中纵剖面对称，重心位置与船模重心重合，如图6所示。该液舱大小为0.6 m×0.3 m×0.25 m，此时舱内装载液体深度为0.125 m。

图6 液舱大小尺寸与位置Fig.6 The size and location of tank

由于考虑了非线性入射波力和非线性静恢复力的影响，船体网格需划分到水线以上一段距离，在gambit中建立四节点船体网格，液舱网格采用三角形面元划分，如图7所示。

图7 船体瞬时湿表面网格与液舱网格示意Fig.7 Panels for ship and tank

3 耦合计算结果与分析

3.1 迎浪工况

图8为零航速迎浪不同频率下有无液舱的纵摇和垂荡运动时历曲线，从考虑非线性因素的计算结果来看，液舱晃荡力矩对纵摇运动的影响较小，随着频率的减小这种差异也逐渐变小。有无液舱的垂荡运动时历峰值存在一个小的差值，但整体零航速迎浪的运动所受影响不大。

图8 零航速迎浪有无液舱纵摇和垂荡运动时历结果Fig.8 The pitch and heave motions in head sea when Fn=0

对于数值计算得到的船体运动时历，需要进行频谱分析结合傅里叶级数展开法对各时历曲线进行拟合，由此将时域计算结果转换到频域进行试验值对比。

已知任意变量可表示为：

n=1,2,…,n

(15)

以迎浪工况波长1.5 m，波幅0.1 m有液舱时船舶纵摇运动时历为例。待运动稳定后，舍去前段选取之后的若干个周期结果进行傅里叶级数拟合，X5可展开为：

X5=X5t0+X5t1cos(ωet+γ1)+X5t2cos(ωet+γ2)…

(16)

(17)

式中：X5为纵摇运动结果；ωe为遭遇频率；X5t0为傅里叶级数展开的定常项即平均值；X5t1为一阶运动响应幅值；γ1为一阶运动响应相位，以此类推。

将时域运动通过傅里叶变换取到3阶精度，相关参数如表2所示。从拟合结果可以分析认为除主频率外其余噪声受非线性计算影响，其结果相对线性得到修正。拟合曲线如图9所示。

表2 傅里叶级数拟合参数Table 2 Fourier series fitting parameters

图9 傅里叶级数展开示例Fig.9 Example of Fourier series expansion method

图10将考虑非线性影响的有无液舱运动与线性时域程序计算结果和船模试验值进行比较，可以发现线性时域计算结果在中高频率不及本文所用方法的结果，再结合试验值可以看出考虑非线性影响后对纵摇的准确预报有了较为明显的改善。

图10 零航速迎浪下非线性和线性数值计算结果与实验值的纵摇RAO对比Fig.10 Comparison of pitch RAO in head sea by experiment and linear and nonlinear calculation when Fn=0

图11～12中给出傅汝德数为0.275时迎浪不同频率下经傅里叶变换得到的有无液舱的纵摇和垂荡运动RAO曲线，试验值为ITTC在1978年给出的S175在Fn=0.275时无液舱的RAO，线性数值曲线为不加载液舱的计算结果。

图11 Fn=0.275迎浪下非线性和线性数值计算结果与实验值的纵摇RAO对比Fig.11 Comparison of pitch RAO in head sea by experiment and linear and nonlinear calculation when Fn=0.275

从图11、12中可以看出，在考虑回复力和入射波力的非线性情况下，与线性处理结果对比均能得到比较满意的纵摇运动预报，并且由于液舱内流体晃荡和波浪诱导力的相互作用，在某些频率范围内相互抵消减小了船体的运动幅值。垂荡情况下线性时域方法所得的峰值结果较大，低频时装载后的非线性方法计算所得结果较未加载状态时小，当频率增大到一定范围时，装载液舱后的运动幅值又较无载液时大，其原因在于液舱诱导力矩与波浪诱导力矩叠加造成。

图12 Fn=0.275迎浪下非线性和线性数值计算结果与实验值的垂荡RAO对比Fig.12 Comparison of heave RAO in head sea by experiment and linear and nonlinear calculation when Fn=0.275

3.2 横浪工况

对零航速横浪时有无液舱的S175在线性与非线性的情况下进行模拟，图13中给出零航速横浪入射波长为9.1、20、40.1 m时有无液舱的横摇和垂荡运动时历曲线。从横摇运动可以看出，在高频时，运动幅值变化不大，非线性的模拟结果较线性有明显改善；在接近共振频率附近时运动幅值明显减小，原因与液舱晃荡力与波浪力在相位和量级上的变化有关；当频率继续逐渐降低时，有液舱的横摇幅值开始增大，船舶的共振频率发生改变。对垂荡运动来说其运动幅值的变化不大，可见液舱晃荡力矩对垂直方向的影响较小。

图13 零航速横浪有无液舱横摇和垂荡时历Fig.13 The roll and heave motions in beam sea when Fn=0

图14分别给出了不同波长的入射波激励下，在液舱晃荡与船体耦合运动过程中，液舱诱导力矩与波浪力的时历对比。其对应的无因次化频率分别为1.5、1、0.7。

在短波时液舱晃荡力与波浪激励力量级相差较大，此时液舱晃荡对船体运动造成的影响很小，非线性与线性结果差异不明显。当频率位于不加载液舱船舶的共振频率附近时，液舱晃荡力约在波浪力幅值的1/10～1/3内，此时对船体运动影响较大，当波浪为20 m时，可以看出液舱晃荡力的相位与波浪激励力存在约180°的相位差，此时二者叠加船舶所受合外力的幅值减小，而船舶横摇运动的幅值也有了明显的减小。当波长逐渐增加，到达载液船舶的共振频率附近时，如图14(c)，可以看出此时波浪液舱晃荡力与波浪激励力仍处于一个量级，且相位发生改变，线性与非线性外流场运动传递后造成的液舱晃荡力基本一致，平稳后与波浪激励力可见峰峰叠加的现象，船体所受总的力增加，此时船舶横摇运动更为剧烈，对应频率下船舶横摇的 RAO 也有所增加，船舶容易处于危险状态。

图14 横浪工况横摇波浪力与液舱晃荡力比较Fig.14 Comparison between sloshing force and wave force in beam sea

图15给出计算多个入射波频率后经傅里叶变换得到的有无液舱船舶在线性和非线性下的横摇运动RAO曲线。液舱流体晃荡使得S175 船模的横摇共振频率区间发生了偏移，船舶加载液舱后横摇惯性半径变大，重心位置和初稳性高也发生改变，横摇固有周期增大，相应的共振频率减小，这与模拟结果一致。线性计算与非线性差距不大，其横摇运动峰值主要与横摇的一次和二次阻尼系数相关。船舶进水后，非线性计算结果相较线性结果有一定的改善，与圆形散点的实验值更为接近，但在无液舱的共振频率附近的模拟仍不是很理想。

图15 横浪下非线性和线性有无液舱数值计算结果与实验值的横摇RAO对比Fig.15 Comparison of roll RAO in beam sea by experiment and linear and nonlinear calculation

4 结论

1)模型虑及瞬时湿表面影响下的入射力与回复力的非线性，内部流动在液舱坐标系下求解，提高了液舱晃荡耦合运动模拟的精度，证实考虑非线性后的结果能得到更为准确的模拟。

2)迎浪有无液舱的纵摇和垂荡运动随航速的变换趋势大致相同，随频率的增大运动幅值差距较小。零航速横浪时的横摇运动在液舱作用下，船体横摇幅值减小，峰值位置前移。

3)耦合运动RAO的变化与液舱晃荡力和波浪诱导力的相位、量级改变有关，当两者相差不大，相位相差180°左右时，横摇运动的剧烈程度得到改善。

4)波浪上船舶液舱晃荡的时域运动预报在线性条件下有一定的局限性，但完全非线性的模拟计算代价巨大，在工程应用上还有一定的距离，而本文所用方法计算快速高效，在帮助处理载液船舶的检验或前期设计中能提供有力帮助。