波浪中船体与液舱晃荡耦合运动的时域数值计算

洪亮，朱仁传，缪国平，范菊，李裕龙

(上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240)

载液船舶航行时，液舱随船体在波浪中运动使其内部的液体产生晃荡，与此同时液舱晃荡力作用于船体，使得船舶的运动姿态发生变化.近年来，随着海洋油气资源开发力度的加大，LNG和LPG等特种液货船型的研制和广泛应用，液舱与船舶之间的作用越来越受人们重视，剧烈的舱内液体晃荡可能会给船体带来巨大的危害，特别是液舱部分装载时，大幅晃荡引起的砰击压力对船体结构可能造成严重的破坏.因此，在载液船舶设计阶段液舱晃荡与船舶运动之间的耦合作用也是需要考虑的重要因素之一.对于船体运动与液舱晃荡耦合作用的问题，国内外已有许多学者进行了分析和研究，如 Kim［1］，Rognebakke和Faltinsen［2-3］，Newman［4］等，对液舱流体晃荡的模拟方法大体分为线性化的频域法(如Newman［5］)和非线性时域法(如Kim［6-7］).从现有的研究结果来看，线性船体运动理论方法基本能够满足船舶运动与液舱晃荡之间耦合效应分析研究的需要.

本文对船舶在波浪上运动和液舱液体晃荡，即对船体内外流场问题均采用了势流理论方法求解，其中波浪中船体水动力和时延函数采用切片法和脉冲响应函数法计算获得，舱内液体非线性晃荡采用时域边界元法计算，并最终建立了在波浪中船体与液舱流体晃荡耦合的运动方程.文中就S175加载方形液舱在迎浪、横浪等不同工况下液舱流体晃荡及其与船体运动耦合分别进行了计算模拟与验证研究.研究中发现，当不出现液面破碎等强非线性现象时，非线性时域边界元法能够给出较好的液舱流体晃荡波形和压力;S175的液舱加载50%的液体时，船体与液舱晃荡耦合运动时历结果能清晰地反映液舱晃荡对船体运动的影响，运动RAO能反映出不同频率液舱对船体运动的影响程度;载液S175船横摇运动RAO能准确给出船体有无加载液舱时共振频率的偏移现象.

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

1.1 载液船舶在波浪上的运动

航行船舶的运动计算所采用的参考坐标系为oxyz，如图1所示，oxy平面与静水面重合，oz轴垂直向上，船舶以定速U0航行，在波浪激励下作六自由度运动.

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

直接求解满足定解条件的速度势不是一件很容易的事，基于切片理论的STF法可以不必求解三维有航速相应的边值问题，只要求解船舶各个横剖面的二维零航速水动力问题即可，STF法的详细介绍可见文献［8-9］.

若将船体横剖面的速度势φ沿船长方向积分即可获得整船的附加质量和阻尼系数:

那么船体五自由度运动方程可表示如下:

式中:ω为频率;ξj为船体j模态的运动位移;mij为包含了液舱内液体质量的船体广义质量矩阵;μij(ω)、λij(ω)为频域附加质量和阻尼系数;Cij为船体回复力系数(ω)为作用在船体上的广义波浪力;j= 2，3，…，6分别对应的运动模态为横荡、垂荡、横摇、纵摇和艏摇.切片法不能计及船体纵荡运动.

1.2 基于势流理论的液舱晃荡问题

采用三维势流理论对液舱内非定常流体晃荡问题进行求解，液舱流体晃荡是采用时域边界元法进行模拟计算的.

液舱内流体运动计算采用的坐标系如图2所示，是固定在液舱内的局部坐标系，其坐标原点固定在静止时自由液面的几何中心，x轴正方向与液舱长度方向平行指向右端，y轴正方向与液舱宽度方向平行指向前端，z轴正方向为竖直向上，在计算过程中该坐标系始终与液舱长宽高保持平行，其原点保持不变.假定流体不可压缩，无粘，流动无旋，则流体运动速度势满足拉普拉斯方程，即控制方程满足:

图2 液舱坐标系示意Fig.2 The coordinate system of the tank

三维流体晃荡运动的流场的自由表面上动力学边界条件为

其运动学边界条件:

壁面条件:

式中:φ为舱内流体运动速度势，Ω为液舱角速度矢量，V为液舱线速度矢量，r为位置矢量，n为边界外法向，ζ为自由面升高.

本文采用简单格林函数法来求解Laplace方程，取1/r为格林函数，则流场中任意一点的速度势可以表述如下:

式中:α为源点处的欧拉角，通过对空间的离散求得每个时间步中的速度势，再通过对自由表面动力学条件进行中心差分获得下一个时间步的自由面条件，从而使得液舱晃荡模拟可以在时域中反复进行下去.本文的空间离散采用了线性面元，积分计算中源点与场点重合时积分为假性奇点，采用解析法进行了计算，源点与场点不重合时采用了七点高斯积分公式计算.计算中根据每个时间步获得的速度势，通过拉格朗日积分计算获得静压力与动压力.

1.3 脉冲响应函数方法的船舶运动时域计算

在线性系统中，任意激励可以写为脉冲响应函数和激励的卷积积分形式［10］:

式中:x(t)为在输入h(t)下的系统响应，F(t)为单位脉冲输入下的脉冲响应函数.将以上概念推广到船舶运动问题，时域船体运动方程有着如下形式:

式中:Mij、μij(∞)和Cij分别代表载液船体质量，无穷大遭遇频率的附加质量以及船舶回复力系数;Kij(τ)为时延函数;(t)为作用在船体上的外部时域波浪力;b41和b42为考虑粘性作用的非线性船舶横摇阻尼系数，可以用经验公式、物理实验等方法获得.

时延函数与无穷大遭遇频率附加质量可由以下频域转时域方法获得:

式中:μij(∞)指频率无穷大时船体的附加质量.

求解运动方程时，横荡与首摇模态的运动可采用数值弹簧技术来抑制其慢漂现象.所加入的横荡和首摇两种运动模态弹簧刚度表达如下:

式中，周期Ti远大于波浪激励周期.

2 载液船舶在波浪上的时域耦合算法及计算步骤

船舶与液舱晃荡耦合运动方程可表示为

在计算模拟过程中，船体当前时刻的运动规律传递到液舱，然后液舱晃荡模拟程序计算出该时刻液舱内部流场的速度势，再通过对舱壁上的压力积分得到液舱的晃荡力和力矩FSlosh，并添加至式(13)的右端，采用四阶龙格库塔法进行耦合求解，此时获得的整船运动位移又作为下一时刻液舱流体晃荡的条件进行计算.由于船体运动预报与液舱晃荡是在不同坐标系中进行的，所以二者之间力与运动需要在这2个坐标系中进行转换，转换方法可参考文献［9］.这样反复循环便可以得到船体与液舱耦合运动的时历，这就是时域下耦合液舱晃荡的船舶全局运动的求解方法.具体的运算步骤如下:

1)读取船型液舱参数及频域水动力数据;

2)根据来波工况计算波浪力时历;

3)计算时延函数，确定计算时长;

4)在时域下求解船体液舱耦合运动方程;

5)转换整船运动规律到液舱，实时给出液舱晃荡的边值条件;

6)进行液舱流体晃荡数值计算;

7)传递液舱激励力(矩)至船体液舱耦合运动方程，转到步骤4)或结束.

3 液舱流体晃荡及其与船体运动耦合的计算与分析

3.1 液舱晃荡数值模拟与验证

液舱流体晃荡力对船体运动直接产生影响，因此非线性液舱流体晃荡问题的准确计算是船舶与液舱晃荡耦合运动系统模拟的前提.基于2.2节中所述时域边界元法，对2种不同尺寸的方形液舱进行了编程计算模拟，并将计算结果与实验、计入粘性的CFD结果进行了比较.

3.1.1 液舱晃荡数值模拟

方形液舱模型尺寸为1 m×1 m×1 m，所装的液体深度为0.5 m，本文对3种工况分别进行了模拟，其中模型液舱受到的激励分别为2个单自由度振荡和一个组合振荡，具体见表1.

表1 模型液舱的受到激励Table 1 Working condition of tank

图3(a)中所示的是工况1中液舱自由面最右端中心点的波面升高时历及其与Faltinsen的试验结果［10］的对比，该工况所选择的振荡频率接近液舱的固有频率.可以看到，在自由液面出现破碎之前数值计算结果与试验结果吻合良好.

图3 液舱右端中心点自由面波高时历Fig.3 Time history of the free surface elevation of tank

图3(b)、(c)分别是工况2、3中液舱自由液面最右端中心点处的波面升高时历，以及考虑了粘性影响的CFD计算结果.可以看出2工况中势流计算结果和考虑了粘性的真实流体的模拟结果吻合良好，两工况模拟计算过程中皆未出现自由液面破碎现象.

3.1.2 晃荡液舱舱壁压力计算验证

用来进行压力验证的方形液舱模型尺寸为1.2 m×1.2 m×0.6 m，其液体装载深度为0.36 m，受到的激励为横荡运动:x=0.015sin(2.475t).

记液舱右端壁面中线上2点为P1、P2，其中P1距离底面0.3 m，P2距离底面0.426 m.模拟计算所得P1、P2压力时历及其与实验值的比较，见图4所示.由图可以看出，液舱晃荡的数值程序计算结果与实验值吻合良好.

图4 P1、P2点处压强时历Fig.4 Time history of the pressure at P1and P2

上述计算结果说明基于时域势流理论的边界元法在不出现自由液面破碎等强非线性现象时可以较准确的进行液舱流体晃荡时域模拟，从而为船体与液舱流体晃荡耦合运动的准确计算奠定了基础.

3.2 船体运动与液舱晃荡的耦合时域数值计算

3.2.1 模型试验与时延函数

模型试验［11］是在中国船舶科学研究中心耐波性水池中进行的，试验模型是耐波性标准船模S175，模型与实船的缩尺比为1∶55，主尺度参数见表2，型线图船模加载的方形液舱尺寸和安装位置见图5.液舱长600 mm，宽300 mm，高250 mm，舱内液体的深度为125 mm液舱位于第9站到13站，重心位置与船模重心重合.

表2 S175实船和船模主尺度Table 2 Principal dimensions of S175 and its modelm

如上文所述，求解船体与液舱流体晃荡耦合时域运动，需预先得到加载液舱的S175船模的频域计算结果，进而利用脉冲响应函数方法获得时延函数.鉴于篇幅，本文只给出部分阻尼系数和时延函数计算结果，图6为横摇频率阻尼系数以及横摇和纵摇的时延函数.耦合运动方程中的非线性横摇阻尼系数b41和b42是根据实验［11］所得到的模型自摇衰减曲线通过最小二乘法得到的.

图5 S175船型线及所加载液舱的尺寸与位置Fig.5 Body plan of S175 and size and location of tank

图6 线性横摇阻尼系数以及横摇，纵摇模态的时延函数Fig.6 Linear damping coefficient of roll and retardation functions of roll and pitch

3.2.2 耦合计算结果与分析

图7给出了迎浪工况下载液S175船模纵摇和垂荡的运动时历.可以看出无论是高频还是低频的入射波激励，液舱流体晃荡对船舶纵向运动的影响不大.图8给出了迎浪工况下，数值和实验所得到的S175船模纵摇RAO，从图中可以看出数值计算与模型实验的结果吻合较好，同时也可以看出迎浪工况下液舱晃荡对船舶纵摇运动的影响不大.

图7 迎浪工况下船体纵摇与垂荡运动时历Fig.7 Time history of pitch and heave motion of ship in head sea

图9分别给出了横浪工况下载液S175船模横摇和垂荡的运动时历.后者处于S175船模的横摇共振频率附近.可以看出对于远离自振频率的高频入射波激励，液舱流体晃荡对船体横向运动的影响不大，而对于自振频率附近的低频入射波激励，液舱流体晃荡明显减小了船体横摇运动的幅值.出现这种现象的原因在于不同入射波激励下，液舱流体晃荡力和波浪诱导力相比有数量级上的差别和相位差.在高频入射波激励时，晃荡力远小于波浪诱导力，对船体运动的影响有限;而在低频入射波激励时，液晃荡力与波浪诱导力处于同一数量级，而且存在180°左右的相位差，在耦合运动的计算中二者相互抵消，明显减小了S175船模的横摇幅值.

图8 迎浪工况下数值与实验所得的S175船模纵摇RAOFig.8 Comparison of pitch RAO of ship in head sea by experiment and calculation

图9 横浪工况下船体横摇与垂荡运动时历Fig.9 Time history of roll and heave motion of ship in beam sea

图10 横浪时S175船模横摇RAOFig.10 Roll RAO of S175 model in beam sea

图11 横浪时船舶横摇时历(周期T=4.735 s，波幅=0.1 m)Fig.11 Time history of roll of S175 model in beam sea (T=4.735s，ζa=0.1m)

图10给出了横浪工况下，数值和实验所得到的S175船模横摇RAO，从图中可以看出数值计算与模型实验的结果吻合较好，同时液舱流体晃荡使得S175船模的横浪共振频率区间发生了偏移.由于对液舱流体晃荡问题采用了基于时域势流理论的边界元法求解，无法很好处理液舱大幅度晃荡时所出现的液面波随等强非线性现象，耦合运动的计算无法达到稳态，图11给出了加载液舱后S175船模共振频率附近的某个频率下横摇运动的时历，可以看出S175船模横摇运动幅值在不断增加，当幅值达到一定程度之后计算便无法进行下去.此外较低的共振频率对应的波浪周期较长，S175船模达到稳态运动所需时间也随之增加，计算误差的积累也是导致计算发散的原因之一.如果能够计算到稳态，其实际运动幅值预计会比图107中所示的幅值要大.所以在图10中，加载液舱的S175船模的横摇RAO的数值结果在共振频率附近的值比实验值要小.

4 结论

本文对船舶在波浪上运动和液舱液体晃荡，即对船体内外流场问题均采用了势流理论方法求解，其中波浪中船体水动力和时延函数采用切片法和脉冲响应函数法计算获得，舱内液体非线性晃荡采用时域边界元法计算，并最终建立了在波浪中船体与液舱流体晃荡耦合的运动方程.文中就S175加载方形液舱在迎浪、横浪等不同工况下液舱流体晃荡及其与船体运动耦合分别进行了计算模拟与验证研究.得出如下结论:

1)当不出现液面破碎等强非线性现象时，非线性时域边界元法能够给出较好的液舱流体晃荡波形和压力.

2)S175的液舱加载50%的液体时，船体与液舱晃荡耦合运动时历结果能清晰地反映液舱晃荡对船体运动的影响.在迎浪工况下，船体运动以纵摇和垂荡为主，由于纵向波浪诱导力矩较大，液舱晃荡力矩与之相比很小，此工况下液舱流体晃荡对于船舶纵向运动的影响不大;在横浪工况下，船舶运动以横摇和垂荡为主，当液舱晃荡力矩与波浪诱导力矩数量级相同时，液舱流体晃荡船体横摇运动产生较大的影响，二者的相位差则决定了液舱流体晃荡的影响效果(增大或减小).从运动RAO中能反映出不同频率液舱对船体运动的影响程度.其中载液S175船横摇运动RAO能准确给出船体有无加载液舱时共振频率的偏移现象.

3)本文只是针对单个液舱且其重心与船舶重心重合的情况进行了计算分析，就已经可以看出液舱流体晃荡对于船体横摇运动有着不小的影响.可以预计，当液舱重心与船舶重心不重合甚至加载多个液舱时，液舱流体晃荡力的变化更大，液舱流体晃荡对船舶横摇运动也会有着更大的影响.

4)本文方法具有较高的计算效率，数值计算与试验结果吻合良好，为设计载液船舶或减摇水舱等的前期设计提供了快速有效的分析方法和技术手段.这种方法的局限性在于无法处理自由液面破碎等强非线性现象，虽然可以较好的预报船体与液舱晃荡耦合运动的共振频率区间，但是无法准确预报共振频率附近的运动幅值.

［1］KIM Y.A numerical study on sloshing flows coupled with ship motion——the anti-rolling tank problem［J］.Journal of Ship Research，2002，46(1):52-62.

［2］FALTINSEN O M，ROGNEBAKKE O F.Sloshing［C］// International Conference on Ship and Shipping Research，Venice，Italy，2000.

［3］ROGNEBAKKE O F，FALTINSEN O M.Coupling of sloshing and ship motions［J］.Journal of Ship Research，2003，47(3):208-221.

［4］NEWMAN J N.Wave effects on vessels with internal tanks［C］//The 20th Workshop on Water Waves and Floating Bodies.Longyearbyen，Norway，2005.

［5］LEE C H，NEWMAN J N.Computation of wave effects using the panel method［M］.Cambridge:MIT Press，2005: 211-251.

［6］KIM Y，SHIN Y S，LIN W M，et al.Study on sloshing problem coupled with ship motion in waves［C］//The 8th International Conference on Numerical Ship Hydrodynamics.Busan，Korea，2003.

［7］KIM Y，SHIN Y S，LEE K H.Numerical study on slosh-induced impact pressures on three-dimensional prismatic tanks［J］.Applied Ocean Research，2004，26(5):213-226.

［8］SALVESEN N，TUCK E O，FALTINSEN O M.Ship motions and sea loads［J］.Trans.of Society of Naval Architects and Marine Engineers，1970，78:250-287.

［9］刘应中，缪国平.船舶在波浪上的运动理论［M］.上海:上海交通大学出版社，1986:151-171，28-31.

［10］FALTINSEN O M.A numerical nonlinear method of sloshing in tanks with two-dimensional flow［J］.Journal of Ship Research，1978，22(3):193-202.

［11］邹康，缪泉明，朱仁庆.耦合液舱晃荡的船舶运动性能研究［D］.江苏:江苏科技大学，2009:37-55.

ZOU Kang，MIAO Quanming，ZHU Renqing.Motion Performance Research of Ship with Sloshing Tanks［D］.Jiangsu:Jiangsu University of Science And Technology，2009: 37-55.