基于SPH方法的二维矩形舱液体晃荡数值研究

卫志军，张文首，王安良，董玉山，胡方源，岳前进

（大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连 116024）

0 引 言

随着超大型浮式液化天然气（floating liquid natural gas，FLNG）船的研发，超大型船载储液舱内的液体晃荡问题成为FLNG 船设计中的关键问题.储液舱内液体晃荡是当外激激励频率与容器内部液体自由液面的固有频率接近时，液体产生的剧烈共振运动.

研究液体晃荡问题主要有理论、实验和数值3种方法.晃荡的高度随机性和非线性，使得理论方法具有一定局限性［1］.实验方法在晃荡问题的研究中仍然是最可靠的方法.不过物理实验设备的成本非常昂贵，而且实验的周期都相对很长，所以实验方法的经济性较差.随着计算机技术的快速发展和计算机容量的迅速扩大，计算机的计算能力得到很大提升且对物理现象的模拟也更加精确，可以采用数值方法开展长时间大量的仿真研究.因此，数值方法可以很好地弥补前二者的缺点，成为研究晃荡问题的主要方法之一.基于网格的数值模拟方法已经在计算力学领域被广泛应用，例如有限元法、有限体积法和有限差分法等.液体晃荡具有强非线性自由液面（波浪破碎和翻卷等），基于网格法的数值方法在模拟自由液面的强非线性时存在相当大的困难且对网格质量有很强的依赖性［2］.而光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）方法是一种无网格Lagrangian算法.该方法于1977年分别由Lucy和Gingold等首次同时提出，起初这种方法主要用来解决星体引力天体物理方面所涉及的相关力学问题［3－4］.当时，这种方法在精度和稳定性方面还存在不足，所以它并未得到广泛的应用.后来Monaghan对SPH 方法进行了改进，首次将其运用到了模拟溃坝、钻削，以及造波问题［5］.SPH 方法不仅省去了网格划分的烦琐步骤，还能很好地模拟自由液面大变形问题.因此，SPH 方法被广泛应用于爆炸、弹道和晃荡等问题的研究.

影响储液舱内液体晃荡的因素很多，如载液率、激励方式、激励振幅和激励频率等［1，6］.很多学者采用SPH 方法对储液舱内液体晃荡开展了研究.陈正云等采用SPH 方法研究二维液舱大振幅非谐振频率激励下液体晃荡的自由液面的波形变化［7］.崔岩等采用SPH 方法研究二维矩形液舱在共振频率附近纵荡、纵摇，以及二者综合效应下波高的变化［8］.李大鸣等采用SPH 方法研究不同横荡振幅对二维矩形液舱非共振频率激励下自由液面的影响［9］.此外，李大鸣等还采用SPH 方法研究远离共振频率的激励频率下，液舱内中高载液率对自由液面的影响［10］.

上述研究工作仅针对整体自由液面的波形或者固定点波高的数值结果与实验结果的对比.而液体晃荡最重要的危害是会产生瞬间增大的砰击荷载，该荷载作用于结构，可能会导致结构失效，因此分析液体晃荡引起的砰击荷载更具有工程意义.Delorme等采用SPH 方法模拟了纵摇激励作用下共振频率附近低载液率流体晃荡［11］，分别与实验测量的最大冲击压力和冲击压力时程对比，但没有考虑水平激励下液体晃荡荷载的特性.基于SPH 方法，Zhu等研究了在纵荡激励、非共振频率和低载液率工况下，液体砰击荷载特性［12］.

本文采用SPH 方法在二维矩形液舱内分别针对大振幅横荡和横摇激励，研究低载液率和自由液面的一阶固有频率f1附近（共振频率和非共振频率）工况下液体晃荡特性，并开展相应的实验研究以验证SPH 方法模拟自由液面和砰击荷载的合理性.

1 SPH 方法基本理论

1.1 控制方程

SPH 理论中将N－S的3个基本方程（连续方程、动量方程和状态方程）进行如下离散：

连续方程

式中：ρi是流体密度；Wij是与粒子i和j的位置相关的核函数；m是质量；v是速度，vij＝vi－vj；xi表示粒子i的坐标.

动量方程

状态方程

式（3）是动量方程的黏性形式.其中pi和pj分别是粒子i和粒子j处的压强，可以通过状态方程（5）得到；g是重力加速度.式（4）是Monaghan型人工黏性hj）；α和β是标准常数，本文中α的范围是0.02～0.10，并且β＝0.式（5）中ρ0 为流体的初始密度；cs是声速；γ是常数，一般取7.

1.2 边界条件

对于固壁边界效应可以用多种方法来引入，本文采用动力边界条件，式（2）、（3）和（5）这3个控制流体粒子的方程同样可以用来描述边界粒子，但这只是计算得到了边界粒子合适的密度，而边界粒子的运动仅依赖于其被定义的运动状态.为了防止边界粒子的非物理穿透，密度低于初始密度的边界粒子，在每一时间步中，其密度都会被自动重置成初始密度.

1.3 密度初始化

为了抑制压力场数值计算的噪声，引入Shepard filter算法，在每计算30 步后，通过式（6）和（7）做一次滤波处理.该算法会在当前时间步对每个粒子做密度重分配，从而光滑压力场.这种处理方法已经在Colagrossi等的工作中被验证［13］.

2 数值算例

2.1 数值模型的建立

采用SPH 方法模拟二维矩形液舱内液体在自由液面共振频率附近的晃荡.数值实验所用二维模型液舱长（L）×高（H）为97.0cm×92.7 cm.冲击压力监测点P1设置在静水面处，如图1所示.设定流体粒子的初始间距为15 mm.计算时间步长和数据采集频率分别为1.0×10－5s和50Hz.

2.2 实验模型的建立

为了验证SPH 方法模拟液体晃荡的合理性，本文在大连理工大学工业装备与结构分析国家重点实验室开展了一系列与数值实验对应的实验研究.实验采用的大尺度六自由度运动平台设计承重荷载为12t，可模拟单个自由度、多个自由度耦合的规则和不规则运动.实验系统如图2所示.

图1 模型几何尺寸（单位：cm）Fig.1 The dimension of model（unit：cm）

图2 晃荡模型实验系统Fig.2 The system for sloshing model experiments

为了尽可能降低数值实验和物理实验的差别，应保证物理模型液舱内流场二维性，物理模型液舱的宽度设计主要考虑以下几个方面［6］：

（1）液舱宽度应足够小，以保证全局流场尽量表现出二维特性，避免三维流场效应（漩涡等）的影响.

（2）液舱宽度要足够大，以保证可以忽略沿着液舱宽度方向上产生的边界层厚度对测试结果的影响［14］.

为了保证流场的二维性，欧洲多家研究机构在全球晃荡模型实验基准（sloshing model test benchmark，SMT）下开展了大量的数值计算，并给出SMT 实验模型液舱合理尺寸范围［15］.

基于上述因素，本文采用的矩形液舱宽（B）为15.8cm，物理模型的舱长和高与数值模型一致.该模型液舱采用透明的有机玻璃制成，厚度为10mm，可视为刚性模型.

2.3 模拟工况

设计的横荡和横摇简谐激励信号如下式：

式中：ηa 为激励振幅；f为外激激励频率；t为计算时间；φ为初相位；Y为激励位移.本文中横荡激励振幅η2a ＝30mm；横摇激励振幅η4a ＝3°，横摇中心在舱底中心处；初相位φ＝π／2；计算时间t＝30s.

通过边界条件、线性自由液面运动和动力条件求解二维Laplace控制方程，得到二维矩形液舱自由液面的固有频率计算公式：

式中：g为重力加速度；h为载液高度；L为液舱的特征尺寸；n为阶数；ωn为自由液面固有角频率.

式中：fn是n阶自由液面固有频率.

静水面高度（h）为18.6cm，即模型液舱载液率（h／H）为20%.详细的实验工况见表1.所有物理实验工况均与SPH 数值实验相同.

表1 实验工况Tab.1 Test condition

2.4 砰击压力监测

采用压力监测半径内所有粒子的压力平均值代表监测点的压力值.压力监测半径R为3倍的初始粒子间距（如图3所示）.

图3 压力监测点及其计算半径Fig.3 Pressure monitoring point and its calculated radius

物理实验中采用孔径为5mm 的薄膜型压阻式压力传感器监测P1点处的冲击压力［6］.

3 结果分析与讨论

3.1 自由液面波形对比

图4给出了横荡激励和1.0f1下，冲击侧壁之前的SPH 模拟的全局自由液面波形和相应的实验波形.在共振频率下，液体在冲击侧壁之前，波面发生了明显的卷波现象.因此，波面水头砰击侧壁的瞬间会夹杂一些气团.通过数值模拟发现气团仅发生在1.0f1附近.

横摇激励和1.1f1下，4个连续时刻中SPH模拟的全局自由液面波形和实验结果具有较好的一致性（如图5所示）.数值模拟的冲击侧壁前后的波面均可与实验波形吻合.

图5 横摇激励和1.1f1 下4 个连续时刻（18.0、18.2、18.4、18.6s）SPH 和 实验的自由液面对比Fig.5 The comparison of free surface from SPH and experiments for four continuous moments （18.0，18.2，18.4，18.6 s）under 1.1f1and roll forced excitation

3.2 冲击压力对比

对20%的低载液率工况，液体对垂直侧壁的砰击较严重.分别对横荡和横摇激励下，数值计算和实验监测的压力监测点P1处的压力时程进行对比分析和讨论.

横荡激励下，4个激励频率在P1处监测到的压力时程曲线对比如图6所示.除了共振频率1.0f1以外，SPH 在其他激励频率下模拟所得压力时程曲线与实验结果吻合良好.在共振频率下，数值结果与实验结果出现较大差异的原因主要是因为共振频率下，液体砰击壁面的物理过程非常复杂，波面水头砰击侧壁时会夹杂一些气团（如图4所示）.物理实验中，由于气体和液体的可压缩性，气团会对液体冲击壁面起到缓冲作用，从而有效降低了冲击压力.而本文所采用的是单相SPH数值模拟，没有考虑气体相，所以液体在冲击壁面时，就不可能有气团来起到缓冲作用.因此建议采用两相流来模拟共振频率下液体剧烈的晃荡运动.

当激励频率为1.2f1时，数值和实验均体现出“拍”现象（beating effect）.此外，在非共振频率1.1f1处，放大一个冲击过程，如图7所示.SPH可以模拟砰击荷载的双峰冲击特性.第一个峰值代表液体冲击舱壁产生的砰击荷载（动压力）.第二个峰值是由于液体砰击壁面后，沿着舱壁爬升，当运动速度为零时，动能最小，势能最大导致的.因此第二个峰值表征静水力.此外，从图7中可以看出，在同一个周期内，SPH 数值模拟得到的压力第二个峰值要高于实验结果，而且这种现象普遍存在于相应激励频率下的整个冲击过程.出现这种现象的主要原因，是由于实验所用的模型液舱具有一定宽度，因此会有三维效应.通过实验观测发现液体在沿垂直壁面爬升运动过程中会出现沿宽度方向的发散运动，这种发散会使得回落到P1位置处的液体水团质量减小，从而减低了液体水团回落时对测点P1处的冲击压力.

图7 横荡激励和1.1f1 下压力监测点P1处的压力时程放大图Fig.7 Enlarged view of the pressure time history at P1under 1.1f1and sway forced excitation

横摇激励下，压力监测点P1 处的压力时程对比如图8所示.横摇激励的冲击荷载较横荡激励下更显著.SPH 模拟的数值结果较实验结果要大，分析认为，主要原因是液舱三维效应及液体剧烈砰击舱壁时瞬间产生的气体和液体压缩作用.

4 结 论

（1）SPH 方法可以很好地模拟低载液率横荡和横摇激励下全局自由液面的波形，如水跃、破碎波等自由液面的大变形现象.

（2）SPH 方法模拟共振频率下的砰击荷载比实验结果要大，主要原因是本文忽略气体相影响.因此建议采用两相流来模拟共振频率下液体剧烈的晃荡运动.

（3）SPH 方法可以很好地模拟横荡和横摇激励下非共振频率处晃荡荷载的主要特点.因此，该方法可用于模拟和评估非共振区域晃荡荷载的特性.

（4）采用SPH 方法对共振频率激励下储液舱内液体晃荡荷载的特性进行研究发现，液体砰击舱壁前会发生水跃现象，且晃荡冲击荷载具有双峰冲击特性.

［1］ 卫志军，岳前进，阮诗伦，等.矩形液舱晃荡冲击载荷的试验机理研究［J］.船舶力学，2012，16（8）：885－892.WEI Zhi－jun，YUE Qian－jin，RUAN Shi－lun，etal.An experimental investigation of liquid sloshing impact load on a rectangular tank［J］.Journal of Ship Mechanics，2012，16 （8）：885－892.（in Chinese）

［2］ ZHU Ren－qing，WU You－sheng，Incecik A.Numerical simulation of liquid sloshing—a review［J］.Shipbuilding of China，2004，45（2）：14－27.

［3］ Lucy L B.A numerical approach to the testing of the fission hypothesis ［J］.The Astronomical Journal，1977，82：1013－1024.

［4］ Gingold R A，Monaghan J J.Smoothed particle hydrodynamics— Theory and application to nonspherical stars［J］.Monthly Notices of the Royal Astronomical Society，1977，181：375－389.

［5］ Monaghan J J.Simulating free surface flows with SPH ［J］.Journal of Computational Physics，1994，110（2）：399－406.

［6］ 卫志军，岳前进，张文首，等.大尺度储舱液体晃荡砰击压力测量方法研究［J］.中国科学：物理学 力学 天文学，2014，44（7）：746－758.WEI Zhi－jun，YUE Qian－jin，ZHANG Wen－shou，etal.Experimental investigation of violent slamming pressure in large－scaled tank［J］.Scientia Sinica：Physica，Mechanica ＆ Astronomica，2014，44（7）：746－758.（in Chinese）

［7］ 陈正云，朱仁庆，祁江涛.基于SPH 法的二维液体大幅晃荡数值模拟［J］.船海工程，2008，37（2）：44－47.CHEN Zheng－yun，ZHU Ren－qing，QI Jiang－tao.Numerical simulation of sloshing in two dimensional liquid tank based on SPH method ［J］.Ship ＆Ocean Engineering，2008，37（2）：44－47.（in Chinese）

［8］ 崔 岩，吴 卫，龚 凯，等.二维矩形水槽晃荡过程的SPH 方法模拟［J］.水动力学研究与进展，2008，23（6）：618－624.CUI Yan，WU Wei，GONG Kai，etal.Numerical simulation of sloshing in two dimensional rectangular tanks with SPH ［J］.Chinese Journal of Hydrodynamics，2008，23 （6）：618－624.（in Chinese）

［9］ 李大鸣，陈海舟，张建伟，等.基于SPH 法的二维矩形液舱晃荡研究［J］.计算力学学报，2010，27（2）：369－374.LI Da－ming，CHEN Hai－zhou，ZHANG Jian－wei，etal.A study of 2Dliquid sloshing in rectangle tanks based on SPH method［J］.Chinese Journal of Computational Mechanics，2010，27（2）：369－374.（in Chinese）

［10］ 李大鸣，李玲玲，高 伟，等.液舱内液体晃荡的SPH 模型［J］.哈尔滨工程大学学报，2012，33（1）：37－41.LI Da－ming，LI Ling－ling，GAO Wei，etal.A model of liquid sloshing in tanks based on the smoothed particle hydrodynamics （SPH）method［J］.Journal of Harbin Engineering University，2012，33（1）：37－41.（in Chinese）

［11］ Delorme L，Colagrossi A，Souto－Iglesias A，etal.A set of canonical problems in sloshing，Part I：Pressure field in forced roll—comparison between experimental results and SPH ［J］.Ocean Engineering，2009，36（2）：168－178.

［12］ ZHU Ren－qing，CHEN Zheng－yun，WANG Quan.Numerical simulation of 2Dsloshing in liquid tanks based on SPH method ［C］／／ASME 2010 29th International Conference on Ocean，Offshore and Arctic Engineering，OMAE2010.New York：American Society of Mechanical Engineers，2010：961－966.

［13］ Colagrossi A，Landrini M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics ［J］.Journal of Computational Physics，2003，191（2）：448－475.

［14］ Faltinsen O M，Firoozkoohi R，Timokha A N.Steady－state liquid sloshing in a rectangular tank with a slat－type screen in the middle：Quasilinear modal analysis and experiments ［J］.Physics of Fluids，2011，23（4）：042101.

［15］ Loysel T，Chollet S，Gervaise E，etal.Results of the first sloshing model test benchmark ［C］／／Proceedings of the 22nd （2012）International Offshore and Polar Engineering Conference.Mountain View：International Society of Offshore and Polar Engineers，2012：398－408.