第46卷第3期2022年6月武汉理工大学学报(交通科学与工程版)Journalof Wuhan UniversityofTechnology(Transportation Science&Engineering)Vol46No3Jun2022非规则激励下浅水液舱晃荡数值模拟石留斌李廷秋"苏>(武汉理工大学船海与能源动力工程学院武汉430063)摘要：为探究非规则激励作用下浅水液舱晃荡运动规律，文中采用Boussinesq数值模型研究船舶横荡运动时液舱内部晃荡运动，使用JONSWAP谱构建外部非规则激励，完成非规则激励下二维液舱晃荡运动的数值模拟.数值计算结果表明，舱壁处浅水液舱晃荡运动能量大小与有义波高成正相关，且波能分布集中于奇数阶固有频率附近.关键词：非规则激励；浅水液舱晃荡;Boussinesq方程中图法分类号：U663.85doi：10.3963/j.issn.2095-3844.2022.03.0110引言随着液化石油气(LPG)船舶，液化天然气(LNG)船舶及超大型液化天然气(VLGC)船舶的需求不断增加，船舶液舱晃荡现象逐渐成为研究热点.FaltinsenE采用解析、数值，以及实验方法对液舱晃荡运动问题进行了系统研究.Ebrahimian 等閃采用边界元方法确定带有挡板的轴对称容器中对称和非对称晃荡固有频率和振型.Cho等〔3*运用有限元方法(FEM)，计算分析了在不同充液率和激励振幅时二维液舱晃荡运动，并将数值计算结果与线性理论计算结果进行了比较.Wu 等4基于二维NavierStokes方程运用有限差分法(FDM)，研究在纵荡和横荡耦合激励条件下，内部挡板结构对晃荡运动的影响并通过实验验证数值结果的可靠性.卫志军等5采用光滑粒子水动力学法(SPH)对二维矩形液舱晃荡运动进行研究，结果表明:SPH方法可以很好地模拟液体晃荡时水跃、破碎等自由液面的大变形运动.陈翔等6将移动粒子半隐式法(MPS)与图形处理器(GPU)并行加速技术相结合，对LNG船舶液舱晃荡数值模，将LNG与型液舱的晃荡运动进行对比.结果表明：在高充液率下LNG型液舱可以有效地减小晃荡运动幅值和壁面砰击压力，但在中低充液率下，LNG型液舱则会加剧晃荡运动.王庆丰等⑺采用去奇异边界元方法，在时域内建立液舱内不规则晃荡运动数值模型，模拟了规则激励下液舱晃荡运动并与解析解进行对比验证其准确性，在此基础上完成了非规则激励下液舱晃荡运动的数值模拟.宁德志等8采用完全非线性边界条件的时域数学模型对纵容中的晃荡，结:相对于数阶有，晃荡对奇数阶固有频率更为敏感.薛米安等9基于振动台对液舱晃荡运动进行实验研究，发现在波浪破碎等强非线性作用下，实际一阶共振频率大于理论推导的一阶固有频率.杨志勋等)10*采用缩尺比为1：1、1：2、1：3,载液率为20%的系列GTT液舱，进行长时间激励下二维液舱晃荡运动尺度效应研究，结果表明：气液密度比是导致模型和原型结果之间存在较大差异的重要原因.文中采用高精度全非线性Boussinesq型方程建立晃荡数值模型，其在非规则外部激励作用下液舱内部晃荡运动特征,对深入理解非线性浅水液舱晃荡运动机理具有重收稿日期：2022-03-16第一作者：石留斌(1997—),男，硕士生，主要研究领域为船船浅水液舱晃荡运动"通信作者：李廷秋(1962—)，男，博士，教授，主要研究领域为CFD与船船与海洋工程水动力学基金项目：国家自然科学基金重点国际合作研究项目(51720105011)第3期石留斌，等：非规则激励下浅水液舱晃荡数值模拟・439・要意义.1数值计算模型(1)(2)(3)选取长度为l 、静水水深为h 的二维矩形液舱，见图1.惯性坐标系oyz 原点位于静水面中点，固定于大地上，z 轴竖直向上.假定矩形水箱 内为无旋无黏不可压缩流体，流场中速度势/满 足下列条件：=2/ — 0■7— —v •兀>n0* /z + 0y/y — 0 z = 0◎ + ⑵▽/・▽ / + g 0 — 0 z = 0 (4)式(1)〜(4)分别表示拉普拉斯方程(=2 =22对+左)、物面不可穿透条件(定义在液舱内部 湿表面上)、自由表面运动学和动力学边界条件•v =!v b ,v c )为液舱速度，vf ,v °分别为水平方向和垂直方向速度分量;比为物面单位法向量，由液舱内指向液舱外；0为自由表面高度；t 为时间；g 为 速度图1二维液舱晃荡示意图随体坐标系固定于液舱上，随液舱一起运动.对于自由表面条件，将惯性坐标系下定义的时间 导数通过以下表达式转化为随体坐标系下的时间 导数.0 : 0— v ・=0(6)>t >t /故随体坐标系下的自由表面运动学和动力学边界条件为0t — v • = 0— /z + 0y /y — 0 z — 0 (7)/t — v ・▽ / + 1- ▽/・▽ / + g 0 — 0 z — 0(8 )本节后续方程均定义在随体坐标系下.总速度势/分为两个部分：/=,+ ,,其中速度势特解,满足拉普拉斯方程和物面不可穿透条 件•若仅考虑横荡运动和升沉运动，那么,—yv b + v(9)另一部分速度势，采用Boussinesq 方程进行求解.数值求解过程先将速度势①= /(y,0,t )y ，垂向速度w = (/z )=0和水平速度5 =(乌)z =0直接定义在自由表面上.基于链式法则定义导数为y — /(y,0,t )y — (/y )z —0 + (/z)z —0・ 0y —u + w 0y (10)— /(，0，t )t — C/t )z —0 + (/z )-0・ 0t—C/t )z —0 + @)v : =0 + /w — 0y /y * —(/t) z —0 + + (b + v) * (11)代入随体坐标系下自由表面边界条件可得：01 — v ・=0 + /z — 0y •/y —v b 0y ― ?y 0y + (1 + 0?)@ (12)1/t — ◎+ /z 0t — v ・ / —⑵▽/・▽ / — g 0 +/z )v b 0y — ?y 0y + (1 + 0?) @* — v b /y + vw —g 0 —1-/y+ 1-@2 (1 + 0?)(13)假设初始时刻液舱内部流场速度势/为0, 自由液面处于静止状态•基于式(12)〜式(13)获得时间导数0和？'，即可采用四阶龙格库塔方法迭代求解不同时刻0和？数值•式(12)和式(13) 中空间导数/y 和％采用五点中心差分法求解. 垂向速度@求解分为两部分@ = ,1z + ,2z ，其中 ,z 通过 对速度 ,2 获,2z 则通过Boussinesq 方程计算得到.根据文献[11*中的求解方法，速度势导数,l z 求解过程如下•首先速度势以及速度势垂向导数近似表达式 为式(14)和式(15)，式中算子丿01和丿02表达式为 式(16)和式(17) ,°表达式为式(18).,](y,z,t ) — ?01 : ,* +?02 : w "(14)w(y ，，) — —— J 02 =2 : , * +?01 : w "(15几-1+(—2+f ⑵)=2 (⑹02 — ° + (- 6 +=2 (7)° — z —z,z =—£(18)式中和w "为未知变量.式(14)满足自由表面速度势为?且式(15)满足底部速度势垂向导数为零：v — J JS : ,"+ J JS : w " z — 0(19)@ =_ 丿影=2 : ," + J b i : w " — 0 z —— h(20)・440・武汉理工大学学报(交通科学与工程版)2022年 第46卷式中J FS 和丿FS 为算子丿°】和 几定义在自由表面 上；？bi 和几为算子丿01和？02定义在液舱底部•将式(19)和式(20)联立可得矩阵如下.算子FS 和丿FS ,?bi 和？b2可以计算得到，卩初始 ，其余 递基于矩阵的解&和v "，代入式(22)可 处.,1z = — J 02 V 2 •7 ?FS ・ v " Z = 0(222非规则激励计算方法采用北海联合波浪谱(JONSWAP ),其波能谱密度函数S(&)表达「S(&) =(H 2&p &—5e —4盂)4 7"(23)式中：(&—&p ) 2"=e --2扎2(24)其中：&p 为JONSWAP 谱图值对应频率，当频率&5&p 时，参数2=0. 07；当频率&〉&p 时, 参数2=0.09.参数/取标准值3. 3,参数H s 为有义波高，系数(的取值 等式(25)成立.横标&,纵坐标为谱密度S .H 2 = 16j s(&)d & (25)不规则激励由规则激励叠加获得，其位移和 速度 达y () = (&()t 7 1() ) (26)v(t) =()&()cos(&()t 7 1 ( i ) ) (27)i式中)为幅值;1为随机相位.)与谱密度函数S(&)存在以下关系：)=/2S(&)・△&(28)式中：△&为频率间隔.考虑到开始阶段速度变化过大，添加渐进函数/(),对 谱峰周期 修正，以 二计算结果的稳定性.sin()-) t 0 t/() = , 2t、1 t 〉t 0t 0 = 10)/&p式中t 为前五个谱峰周期.则有最终结I 下.ytt) = /()(&()t 7 1()) (31)29)30)3() = /()#)()&()cos(&()t 7 1(i))(32)舱内液体一阶固有频率(&i=2. 04 rad/s ,根 据色)作为谱 ，则有谱 i 率&p =&i =2. 04 rad/s.频率&取值范围为0. 01〜5rad/s ,频率间隔△& = 0. 01 rad/s.参数/取标准值3. 3,有义波高H s = 0. 012 m.波能谱密度函数 S(&)化(见图2),不规则速度激励图像见图3,并选取其中t = 0〜30 s 表明渐进函数/()在前五个周期的缓冲效果(见图4).图2波能谱密度函数S(&)随频率变化图添加渐进函数＜(/)— 不添加渐进函数几/)|(_‘s ・ U n0 10 20 30 40 50 60 70 80 90 100t/s图3不规则速度激励随时间变化图(_‘s ・ U n添加渐进函数＜(/)……不添加渐进函霸同5 1015 20 25 30t/s图4不规则速度激励中添加与不添加渐进函数/()对比图3数值模拟结果3.1 规则激励验证该数值模型的有效性，在相同工况下与文献结果(采用固有 模块叠)进行对比.二维矩形液舱长度L = 1. 175 m,水深力=0.06 m,受到长度方向简谐激励，激励幅值)=0. 003 9 m.基于线性色 获 晃荡运动一阶固有频率&i=2. 04 rad/s ,激励频率按文 献选取 &p 1 = 0. 98&i ，&p 2 =1. 02&i ，&p 3 =1. 08&i .图5为激励别为 时二维液舱壁 高度变化对比图,横、纵坐标另I] 与激励周期 t/T 和高第3期石留斌，等：非规则激励下浅水液舱晃荡数值模拟・441・度与水深之比“/h.在激励&p1=0.98&1时，中可以观察到波波峰在舱壁I来回，舱壁波高时历曲线见图5a）.在&p2=1.02&i和&p3=1.08&i时，液舱分别观察到两波波峰和波波峰，波高时历曲线见图5b）〜c）.由图5，文中用数值模型能准确地有附近高度变化规律•3.2不规则激励数值模拟结果基于色散关系式计算获得液舱晃荡运动一阶至五阶固有频率分别为&i=2.04rad/s,&2= 4.03r a d/s,&3=5.93r a d/s,&4=7.71rad/s，&5 =9.34rad/s.下面对谱参数对液舱晃荡运动的影响，采用matlab数值模拟,程序运行时间t=600s,网格点数取N=39（长度C=1.175m）,时间步长=0.02s,matlab自带能量谱密度函数pwelch函数对液舱舱壁处的自高程历线谱.选取海浪谱谱峰频率&p=&i=2.04r a d/s,参数*=3.3#有波高H s=0.008m图6为舱壁处自由液面高程时历曲线图，图6中高程化呈规则趋.1.0r!=2.04/(rad•s1)_0'50100200300400500600t/s图6舱壁处自由液面高程时历曲线图（=&i）谱参数有义波高H s的取值其对液舱晃荡运动的影响，图7为&p=&i=2.04 r a d/s,参数*=3.3,三种有义波高H s舱壁处自高程历线谱图.图7：工况呈相规律.波谱密度图谱于一阶有附近#其阶有附近在部分波能，且波能逐阶递减的趋势，在阶固有波忽略不计•三种工况对比发现波能大小与有义波咼成正相关，有义波咼越大波大•…H5=0.012m—H5=0.008m—皆0004m图761X1「5-s4-S3-<2s1°1CD,CO3OJ A co512co/(rad•s_1)不同H s取值舱壁高程时历曲线谱分析图研究海浪谱谱&p与液舱晃荡运动的关系，图8谱谱别选取一阶:阶固有频率&p=&1，&p=&2，&p=&3p=&4，&p=&5，参数*=3.3,有义波咼H s=0.008m时壁高程时历曲线谱图•相对于数阶固有而言，液舱晃荡对奇数阶固有，液舱晃荡波集中于—?、三、五阶•10"rxio3,co,co3co A cd51214co/(rad•s_1)a)aj=aj=2.Q4rad/s2-OrxiO-31.51.00.5-21co1co,co3co A cd512〔co/(rad•s_1)b)aj=aj2=4.Q3rad/sco1co,co3co A co5121ml(rad•s_1)d)aj=aj=7.71rad/s(co,co3co A co51214co/(rad•s_1)c)aj=aj.=5.93rad/s,co,co3co A cd51214co/(rad•s_1)e)co=co=9.34rad/sIO3图8不同&取值舱壁处自由液面高程时历曲线谱分析图对谱等于共振，以&p=1.8rad/s和&p=2.3rad/s为例（见图9）,与一阶共振情况&p=&i=2.04rad/s对比,参数*=3.3,有义波咼H s=0.008m.结I明，工况下谱图像呈现出类似规律，谱1于一阶固有附近，其他阶固有附近在波能，且波能逐阶递减的趋势.a) 1.8rad/s「6■X104:4;1目-1a H舟0（■A人心」)24681012a>/(rad•s_1)「1.0's冷0.5024681012a>/(rad•s_1)b)a>=2.3rad/s图9壁高程历线谱图4结论1）当谱等于一阶固有，波能大・442・武汉理工大学学报(交通科学与工程版)2022年第46卷小与有义波高大小成正相关，有义波高越大波能越大.当谱峰频率等于或偏离一阶固有频率时，波能主要集中于各阶固有频率附近，波能在一阶固有频率附近最大且从一阶至五阶波能逐渐降低，五阶以后可忽略不计.2)当谱峰频率等于二阶及以上固有频率时#浅水液舱晃荡波能集中于三阶或五阶固有频率附近.当谱峰频率为二阶、三阶固有频率时，浅水液舱晃荡运动波能均集中于三阶固有频率附近.当谱峰频率为四阶、五阶固有频率时，浅水液舱晃荡运动波能均集中于五阶固有频率附近.参考文献FALTINSEN O M.Sloshing[M].Cambridge:Cam-bridgeUniversityPress#2009EBRAHIMIAN M#NOORIAN M A#HADDAD-POUR H A successiveboundaryelementmodelfor investigationofsloshingfrequenciesinaxisymmetric muti baffled containers[J*.Engineering Analysis with Boundary Elements#2013,37(2)：383-392.)*CHO J R#LEE H WJ.Non-linear finite element analysis of large amplitude sloshing flow in two-dimensional tank)*.International Journal for Numeri-calMethodsinEngineering#2004#614):514-531)*WU C#FALTINSEN O M#CHEN B.Numerical studyofsloshingliquidintankswithba f lesbytime-independentfinitedi f erenceandfictitiousce l method [J*.Computers and Fluids#2012,63:9-26.)*卫志军，张文首，王安良，等•基于SPH方法的二维矩形舱液体晃荡数值研究大连理工大学学报# 2014#54(6):597-603)*陈翔，万德成.MPS与GPU结合数值模拟LNG液舱晃荡力学学报，2019,51(3):714729.)*王庆丰，徐刚，王树齐，等•去奇异边界元方法在液舱晃荡模拟中的应用振动与冲击，2018,37(19)：69-7396)*宁德志，宋伟华，滕斌.纵摇容器中液体晃荡的非线性数值模拟船舶力学，2017,21(1)：15⑵2.)*薛米安，邢建建，陈奕超，等•基于振动台实验的液体晃荡激励参数敏感性研究大连理工大学学报# 2019,59(2):162-171.)0*杨志勋，徐潜岳，吴尚华，等.长时间规则激励下的二维液舱晃荡试验尺度效应研究中国造船#2017#58(2):156-164)11*SU Y#LIU Z Y Numerical modelofsloshingin rectangulartankbasedonBoussinesq-typeequations)J*OceanEngineering#2016(1):166-173NumericalSimulationofSha l ow-waterSloshingMotion UnderIrregularExcitationsSHI Liubin LI Tingqiu SU Yan(School of Naval Architecture#Ocean and Energy Power Engineering#Wuhan University of Technology#Wuhan430063,China)Abstract:In order to explore the sloshing motion law of shallow water tank under irregular excitation# Boussinesqnumericalmodelwasusedtostudythesloshing motioninsidethetank whentheshipis swaying JONSWAP spectrum was used to construct external irregular excitation#and the numerical simulationofsloshingmotionoftwo-dimensionalliquidtankunderirregularexcitationwascompleted The numerical results show that the sloshing energy of the shallow water tank at the bulkhead is posi-ivelycorrelated withthesignificant waveheight#andthe waveenergydistributionisconcentrated neartheodd-ordernaturalfrequencyKey words：irregular excitation；shallow-water sloshing；boussinesq-type equation。
