双浮板液舱晃荡特性的数值研究

于曰旻

(1. 海南大学 土木建筑工程学院，海南 海口 570228； 2. 海南大学 南海海洋资源利用国家重点实验室，海南 海口 570228)

液舱晃荡是载液舱在外界激励下，在部分载液率的条件下，舱内自由液面产生的剧烈波动现象。当外界激励频率接近于载液舱内自由液面的固有频率时，特别是接近于一阶固有频率时，晃荡变得尤为剧烈，巨大的冲击载荷作用于舱壁和舱顶，往往会引起局部结构的破坏，甚至会引起整个舱体的失稳。近些年，液舱晃荡一直是工程界和学术界的研究热点，如何有效地抑制液舱晃荡也成为近几年的研究难点。

液舱晃荡的研究方法主要包括：理论分析、数值模拟和模型试验。除了带有少量解析解[1]的情况之外，液舱晃荡的研究大都采用模型试验和数值模拟的方法。薛米安等[2]采用模型试验的方法研究了浅水条件下液舱晃荡冲击压力的特性，结果显示一阶固有频率略大于解析解，三阶固有频率略小于解析解，而且液舱晃荡产生的冲击压力随着激励幅值的增加而增大。杨志勋等[3]对于3个不同比尺模型，开展了高载液率液舱晃荡的模型试验研究，冲击压力和对应波形结果表明，原型(比尺1∶1)液舱中的晃荡波携带大量气泡，尤其在静水面附近，因此冲击效果更加显著。此外，唐洁等[4]通过物理模型试验，对弹性侧壁液舱和刚性液舱内液体晃荡问题进行了研究，分析两种相对液深比条件下弹性侧壁液舱与刚性液舱内不同测点的波面、振幅谱和晃动波高。基于OpenFOAM开源代码，李金龙等[5]对比了捕捉自由液面的代数流体体积(VOF)方法和几何VOF方法，以及在液舱晃荡波面数值模拟方面的准确性。刘东喜等[6]通过耦合水平集与流体体积(CLSVOF)方法追踪自由液面，对矩形液舱内三层液体晃荡进行了数值研究，当外界激励频率接近于内部交界面的最低阶固有频率时，其波面会大于自由液面的波动幅度，而且流体界面会产生开尔文—亥姆霍兹不稳定性。基于移动粒子半隐式(MPS)方法，张雨新等[7]对二维矩形液舱晃荡问题进行了数值验证，结果表明MPS方法能够很好地计算晃荡产生的拍击压力。

目前，液舱晃荡的制荡结构主要分为固定式制荡结构和浮动式制荡结构。基于Ansys Fluent，徐毅等[8]研究了带有底部固定挡板圆柱形卧式储罐内的液舱晃荡问题，结果显示，随着挡板高度的增加，舱壁的冲击压力逐渐减小，舱内波面也随之变得平缓。基于MPS方法，Tsukamoto等[9]分析了浅水和中等水深条件下液舱晃荡载荷与底板或者侧壁肋条的相互作用。不同水深条件下，与光滑液舱相比，底板肋条和侧板肋条都起到了制荡作用，而侧板肋条制荡效果更加显著。通过使用有限单元—光滑例子耦合方法，Zhang等[10]数值研究了底板上固定的多个弹性板与液舱晃荡的流固耦合问题。相对于双垂直板、双垂直T型板和双水平—垂直板，三垂直板制荡效果更加显著。Kolaei和Rakheja[11]使用有限单元法来研究在自由液面处带有柔性薄膜的液舱晃荡问题，在水平圆柱形液舱中，随着薄膜张力的增加，固有频率向高频区移动。通过模型试验和解析方法，Ning等[12]和Zhang等[13]分析了在自由液面处布置多层浮球的液舱晃荡响应，随着浮球层数的增加，晃荡波面爬高和冲击压力都随之减小。固定式制荡结构与舱壁形成一体，安装和维修难度较大，为了起到有效的制荡效果，结构较为复杂，这在一定程度上减少了载液容积。浮动式制荡结构可以随着自由液面的变化而变化，但是当激励幅值较大时，往往难以控制。

采用数值模拟的方法研究双浮板液舱晃荡特性。首先详细介绍了数值模型，然后进行了网格收敛性验证，对光滑矩形液舱和双浮板液舱晃荡波面响应进行了分析和讨论，最后给出了相关的结论。

1 数值模型

为了研究双浮板液舱晃荡特性，采用大型商用软件Ansys Fluent[14-15]进行二维数值模拟。不可压缩流体的连续性方程为：

(1)

式中：ui为i方向流速。

动量方程为：

(2)

式中：ρ为流体密度；p为流体压强；ν为流体运动黏性系数；fi为外界作用力。

Ansys Fluent使用有限体积法离散和求解纳维—斯托克斯方程，流体体积模型(VOF)[16]用于捕捉自由表面，边界条件均为固壁。对于强迫激励下的液舱晃荡，可以采用动量源方法[17]。6自由度求解器[18]使用力和力矩来计算物体重心的平动或者转动，通过6自由度求解器可以限制浮板相对于液舱只发生垂荡运动，结合动网格技术可以模拟双浮板与晃荡波面的耦合运动。动网格技术[19]可以用于移动或者旋转物体周围网格并实时更新，双浮板周围的网格均采用三角形网格，网格光顺和网格重划分方法得以配合使用，实时更新内部网格来阻止负网格的产生。

2 网格收敛性分析

采用的计算液舱为矩形液舱，液舱和浮板的尺寸和布置见图1。静止时，浮板位于自由液面处，它的密度为液体密度的一半，通过6自由度求解器，使得两个浮板只能相对于液舱发生垂荡运动，每块浮板到最近侧壁的距离均为L/6(L为液舱的长度)。

图1 计算液舱尺寸和布置Fig. 1 Dimensions of the clean tank and layout

光滑矩形液舱的运动方程为：

U=U0sinωt

(3)

式中：U为外界激励速度；U0为外界激励速度幅值；ω为外界激励圆频率。

分别采用3种不同空间步长的网格离散计算区域，进行网格收敛性验证[20]，分别是M1、M2和M3，离散后的节点总数和单元总数见表1。

表1 计算网格Tab. 1 Computational meshes

在载液率为50%，激励幅值为2 mm，激励频率为0.84 Hz的条件下，对3个不同步长的计算网格分别开展了双浮板液舱晃荡的数值模拟，并提取了距离舱壁17 mm处的数值浪高仪得到的波面时历进行了网格收敛性验证，见图2。结果显示，3种不同网格的最大爬高较为一致，它们之间的相对误差(相对于M1)均小于4%。以上对比表明，该数值模型对于网格离散不敏感，因此，为了确保数值模拟的准确性和节省计算时间，将采用网格M2开展数值模拟。

图2 不同网格波面时历对比(载液率50%，激励幅值2 mm，激励频率0.84 Hz)Fig. 2 Wave elevation against time for three different computational meshes (50% filling rate,2 mm excitation amplitude,0.84 Hz excitation frequency)

3 结果和讨论

光滑矩形液舱n阶固有频率的解析解为[21]：

(4)

式中：fn为第n阶固有频率；kn为第n阶波数；g为重力加速度；d为水深。

数值模拟的计算工况见表2，载液率为50%，激励幅值为2 mm，激励频率范围为0.6～1.1 Hz(等间距为0.05 Hz)，并在一阶固有频率附近进行加密，每组运行总时间为120 s。

表2 计算工况Tab. 2 Computational cases

3.1 光滑液舱最大爬高的频率响应

基于摄动理论，单向平动激励下，光滑矩形液舱晃荡波面的线性解为[22]：

(5)

(6)

当t→∞时，Bn为无穷大，此时共振现象发生。

此外，通过MTS双向振动台开展了光滑液舱晃荡模型试验，台面尺寸为3 m×3 m，激励频率范围为0～50 Hz。光滑液舱的内部尺寸为0.9 m×0.3 m×0.91 m(长×宽×高)，由2 cm厚的透明亚克力板制成，固定于双向振动台台面上。浪高仪放置于液舱内，到液舱侧壁的距离为1.7 cm，用来实时记录波面。

图3给出了载液率为50%，激励幅值为2 mm，光滑液舱晃荡沿舱壁爬高的最大值对频率的响应曲线。图3中横坐标为激励频率，纵坐标为沿舱壁爬高的最大值(数值浪高仪距离舱壁17 mm)。从图3中可以看出，解析解在一阶固有频率0.89 Hz处的爬高最大值趋于无穷大，而数值模拟结果与试验结果吻合较好，爬高最大值先是随着外界激励频率的增加而增大，在一阶固有频率0.87 Hz附近达到最大值后，开始迅速减小。数值模拟结果和试验结果得到的一阶固有频率均略小于解析解。这里的解析解是基于势流理论计算得到的固有频率，忽略了波浪破碎、黏性等因素的影响。

图3 光滑液舱爬高最大值的频率响应曲线(载液率50%，激励幅值2 mm)Fig. 3 Maximal wave runup against forcing frequency in the clean tank (50% filling rate,2 mm excitation amplitude)

3.2 双浮板液舱晃荡的制荡效果

固定式隔板需要整体固定于舱壁上。水平固定式隔板往往需要布置多层，来抑制不同载液率条件下的液舱晃荡，因此加大了安装和维修难度。垂直固定式隔板将液舱分为两个或者多个子液舱，这在一定程度上减少了液舱容积。双浮板位于自由液面处，相对于液舱只产生垂荡运动，可以有效地抑制不同载液率下的液舱晃荡，其安装和维修难度相对较小，液舱的容积率也相对较大。

图4为载液率为50%，激励幅值为2 mm，双浮板液舱晃荡沿舱壁爬高的最大值对频率的响应曲线。当外界激励频率小于自由液面的一阶固有频率前，最大爬高随着激励频率的增加而增大；当外激励频率接近于0.84 Hz时，爬高达到最大值；之后，当外激励频率继续增大，最大爬高逐渐减小。通过数值模拟结果显示，光滑液舱晃荡自由液面的一阶固有频率为0.87 Hz，与此相比，在光滑液舱中添加双浮板之后，液舱晃荡现象得到了明显的抑制，最大爬高减小了83%。此外，最低阶共振频率从光滑液舱的0.87 Hz减小为0.84 Hz，说明添加双浮板后，与光滑液舱相比，一阶共振频率会向低频区移动[23]。

图4 双浮板液舱爬高最大值的频率响应曲线(载液率50%，激励幅值2 mm)Fig. 4 Maximal wave runup against forcing frequency in the tank with dual floating plates (50% filling rate,2 mm excitation amplitude)

图5为数值模拟得到的双浮板液舱晃荡在一个激励周期内不同时刻(T/4、2T/4、3T/4和T)的波面和双浮板位置。如图5所示，在t=T/4时刻，左侧浮板达到垂荡响应的最大值，此时右侧浮板位于其垂荡较小值；当t=2T/4时刻，舱内液体自由液面接近于静水位，双浮板均接近于垂荡的平均位置；随着时间的增加，在t=3T/4时刻，左侧浮板达到垂荡响应的较小值，而右侧浮板位于其垂荡最大值；当t=T时刻，舱内液体自由液面又接近于静水位，双浮板也均接近于垂荡的平均位置，完成一个循环过程。此外，通过4个不同时刻的波面，液舱晃荡模式发生了改变，由光滑液舱的驻波模式[24]变为双浮板的U管模式[25]，这也在一定程度上起到了抑制液舱晃荡的作用。

图5 双浮板液舱晃荡在不同时刻的波面(载液率50%，激励幅值2 mm，激励频率0.84 Hz)Fig. 5 Wave elevation in the tank with dual floating plates during different instants (50% filling rate,2 mm excitation amplitude,0.84 Hz excitation frequency)

4 结 语

为了研究双浮板液舱晃荡特性，基于Ansys Fluent建立了数值模型。通过动网格技术、6自由度求解器和动量源方法开展了双浮板液舱晃荡数值模拟，并进行了网格收敛性分析，通过数值浪高仪得到了载液率为50%，激励幅值为2 mm，较宽范围激励频率下液舱晃荡的波面，并对其进行了分析，得到以下主要结论：

1) 对于光滑液舱而言，液舱晃荡最大爬高和一阶共振频率的数值模拟和模型试验结果吻合较好，而最低阶固有频率比解析解略小。

2) 与光滑液舱相比，添加双浮板后，液舱晃荡的最大爬高得到明显抑制，对于最小固有频率处，最大爬高减小了83%，此外，一阶共振频率向低频区移动。

3) 不同时刻的波面显示，双浮板液舱晃荡模式从光滑液舱的驻波模式变为U管模式，这在一定程度上起到了抑制液舱晃荡的作用。