大型补给舰船液舱晃荡载荷计算方法

金 晶 薛鸿祥 唐文勇 张圣坤

上海交通大学海洋工程国家重点实验室，上海 200240

0 引 言

补给舰船用于在海上航行或停泊中为已方舰艇提供燃料、淡水、食品等消耗品，以及鱼雷、水雷、炮弹、导弹等武器，作战舰艇通过航行补给不仅能延长活动半径，提高在航率，增加海上活动时间，而且还避免了对固定基地的过分依赖。其中，大型综合补给舰船还可为远洋舰队以及航母编队提供强大的支援，受到各海军强国的重视。随着我国海军战略发展的需要，大型综合补给舰船的研究与建造也势在必行。

随着补给舰船吨位的增大，装载燃油等补给的液舱尺度也随之增大，在某些装载情况下，有可能会发生剧烈的液体晃荡谐振运动，引起冲击载荷，从而危及舱壁及舱室内部构件结构的安全。因此，在大型补给舰船的液舱结构设计中，应考虑可能产生的液体晃荡冲击载荷。液舱晃荡冲击载荷预报是当前的研究热点，并已取得一定的成果，采用Fluent等商用软件，可以计算液舱内流体自由液面的变化、流体质点的速度矢量和流场内的压强分布等。例如，丁仕风等［1］基于 Fluent，采用VOF法对速度突变情况下液化天然气船液舱内的晃荡进行了模拟和分析；张书谊等［2］采用Fluent对矩形液舱的横荡流体载荷进行了数值模拟并和实验结果进行了对比。但是，在初始设计阶段或审图阶段，熟练掌握商用CFD软件对于结构设计人员来说并不方便，因此，发展可靠且方便于工程应用的晃荡载荷计算方法，并结合试验数据进行修正仍是当前需重点研究的方向。

大型补给舰船的液舱截面可能为非矩形截面，边界复杂，也可能含有内部构件，因此，其晃荡载荷的计算方法应能处理非矩形边界，并能准确、稳定地预报液舱内部的速度和加速度时间历程，从而准确获得舱壁以及构件上的晃荡载荷。对于非矩形液舱，如棱形液舱的晃荡模拟，传统的VOF方法采用的是将矩形网格的折线近似处理为斜线，从而导致模拟不够准确，引起一定的误差［3］。沈猛等［4］通过引入部分单元参数（Partial Cell Parameter，PCP）的概念对VOF法进行了改进，用于处理棱形液舱，但自由表面重构采用的还是传统的Hirt-Nichols算法。由于Hirt-Nichols算法为零阶精度［5］，对晃荡自由表面的模拟精度较低，从而导致速度、加速度预报结果离散度较大，给液舱内部构件的晃荡载荷分析带来了困难。端木玉等［6］采用Youngs方法对矩形液舱大幅晃荡进行模拟，得出用Youngs法比用Hirt-Nichols法能更准确地处理液舱大幅晃荡的结论，但该文的算法适于处理矩形截面的液舱晃荡计算，对于棱形液舱等复杂边界的情况还有待进一步的完善。本文将以棱形液舱为对象，采用改进的VOF法，基于精度较高的Youngs法重构表面，并结合部分单元参数概念［7］（以下简称“PCP-Youngs法”）模拟晃荡运动，以有效改善自由表面的模拟效果，特别是提高舱内速度和加速度的预报精度以及数值稳定性，为在大型补给舰船液舱结构设计中加入对于晃荡载荷冲击影响的考虑奠定基础。

1 液体晃荡的计算模型

1.1 控制方程

对于如图1所示的棱形液舱边界处的部分不规则网格（网格内有部分障碍物），可采用网格的体积通度来描述网格中可以被流体通过的体积与整个网格的体积之比，用网格的面通度λ来描述网格面上能够被流体通过的面积与该面积之比。基于PCP概念，考虑体积通度λ以后，流体控制方程中的连续性方程、动量方程和体积分数传输方程可分别表示如下。

式（1）～式（3）中，u为流体速度矢量；p为流体压力；F为流体的体积分数；ρ为流体密度；μ为运动粘性系数，f为流体微元体所受到的体积力，该力是引起流体晃荡的主要激励，采用相对坐标系，表达式为：

式中，ω为相应时刻横摇角速度；rP为相对坐标中质点的位矢。

图1 部分单元参数的变量位置与定义Fig.1 The variable distribution and definition in PCP method

控制方程的离散采用有限差分法，差分网格采用交错形式。其中，网格的体积通度λ、压力P和流体体积分数F定义于网格中心，而网格的速度u，v和网格的上、右、下、左4个面的通度 λt，λr，λb，λl定义于网格的4条边上（图1）。

考虑体积通度和面通度后的连续方程具体可参见文献［8］。对于体积传输方程，本文采用了精度较高的Youngs法［9］。

Youngs法基于几何学原理，在单个单元网格内是用斜线段近似模拟界面，其可以分两种情况进行计算：一是计算单元为满网格，F=1；二是计算单元为自由表面网格，0＜F＜1。其中，对于第2种情况涉及界面重构的问题，需要先计算网格内界面的法向n=，其中

根据法向和网格内F函数的大小，可将界面分为20种情况，除=0 和0 所对应的 4种情况外，其余16种可通过对称和翻转运算后简化归并为4大类。确定了界面类型后，便可计算出单元各条边上的流体分数，进而根据n时刻的速度计算出通过上边、右边、底边和左边的流体体积分数的运输量（以输出为准），分别用 Ft，Fr，Fb和Fl表示。由4种情况的流体体积分数的运输量，可以确定n+1时刻的流体体积分数，从而将传统的Youngs的流体体积传输方程修改为：

1.2 自由表面条件

在自由表面上，动力学条件须满足表面应力条件，忽略表面张力，则表面应力条件为：

式中，un为垂直于自由表面的法向速度；ut为切向速度；p0为舱室内气体的压力。

对于自由表面网格速度的求解，采用常数外插与连续性条件和自由表面条件相结合的方法。对于自由表面边界速度，采用常数外插，但当判断出网格转换为不含自由表面边界时，则要求网格满足连续性条件，尽量避免压力和速度的突变。

2 算例及对比分析

为了验证基于PCP的Youngs法（PCP-Youngs法）对棱形液舱晃荡数值模拟的可行性和计算精度，与文献［10］中的棱形液舱模型晃荡试验数据进行了对比验证，同时，还将本文的结果与基于PCP的Hirt-Nichols法的计算结果进行了对比。

选取文献［10］中的30%和46%两种装载水平进行晃荡计算，两种装载工况的横摇激励周期分别为1.517 s和1.207 s，舱室模型的横摇幅值均为5.73°。

棱形液舱模型尺寸及各测试点的位置如图2所示。图中R1～R4分别为舱室边界的压力测试点，P1～P3，Q1～Q3分别为30%和 46%装载率下液面高度、2/3液面高度和1/3液面高度处的监测点，这些监测点位于液舱的中间纵向剖面处，用于监测速度和加速度。

图2 棱形液舱模型Fig.2 Prismatic liquid tank model

2.1 液面形态

装载率为30%和46%，网格为114×78时的液体波面形态分别如图3～图6所示（其中图3角隅处出现的毛刺是由于软件画图设置方面的因素，因软件在画图过程中采用的是简化的矩形的网格划分，并在整个网格内部按照网格的压力填充了颜色）。由图中可看出，采用Hirt-Nichols法会出现非物理性的飞溅，且自由表面也不连续，而采用Youngs法则未出现这些异常，自由液面相对光顺、连续。

图3 Hirt-Nichols法的液面形态（30%装载率）Fig.3 Liquid surface form of Hirt-Nichols method（30%filling ratio）

图4 Youngs法的液面形态（30%装载率）Fig.4 Liquid surface form of Youngs method（30%filling ratio）

图5 Hirt-Nichols法的液面形态（46%装载率）Fig.5 Liquid surface form of Hirt-Nichols method（46%filling ratio）

图6 Youngs法的液面形态（46%装载率）Fig.6 Liquid surface form of Youngs method（46%filling ratio）

2.2 舱室边界压力

图7和图8所示分别为30%装载率时测试点R1和R2的压力历程曲线对比，图9和图10所示分别为46%装载率时测试点R1和R2的压力历程曲线对比。30%和46%装载率下PCP-Youngs法和Hirt-Nichols法预报的R1，R2点冲击压力的最大值以及其与实验数据相比的相对误差如表1和表2所示，其中，30%装载率的最大值为图7和图8所示的3个周期内最大幅值的平均值，46%装载率与之类似。

图7 R1处压力（30%装载率）Fig.7 Pressure of R1（30%filling ratio）

图8 R2处压力（30%装载率）Fig.8 Pressure of R2（30%filling ratio）

图9 R1处压力（46%装载率）Fig.9 Pressure of R1（46%filling ratio）

图10 R2处压力（46%装载率）Fig.10 Pressure of R2（46%filling ratio）

表1 30%装载率下PCP-Youngs法和Hirt-Nichols法预报的R1，R2点处的冲击压力和相对误差Tab.1 The impact pressure and proportional error on R1and R2with PCP-Youngs method and H-N method（30%filling ratio）

表2 46%率下PCP-Youngs法和Hirt-Nichols法预报的R1，R2点处的冲击压力和相对误差Tab.2 The impact pressure and proportional error on R1and R2with PCP-Youngs method and H-N method（46%filling ratio）

由表中可见，PCP-Youngs法的计算结果要明显优于Hirt-Nichols法。从图中可看出，30%装载率下用两种方法计算得到的压力历程曲线都呈现出了双峰特征［11］，其中用PCP-Youngs法计算得到的压力曲线的冲击压力时刻略微提前，其主要原因是本文为直接加载的周期性的运动激励，而模型实验的激励则为由初始静止状态逐渐加至稳定周期的过程，因此，该现象符合实际情况。通过比较PCP-Youngs法和Hirt-Nichols法的结果发现，在30%装载率时，PCP-Youngs法在R2处的冲击压力峰值（图8）更接近于实验数据，相对误差也小一些（表1），而Hirt-Nichols法则仅在R1处所示的第2和第3个周期中更加接近于实验数据；在46%装载率时，PCP-Youngs法在R1和R2处的压力峰值均较接近于实验数据（图9、图10），预报的压力曲线也相对较光顺、连续。可见对于舱室边界压力的预报，从总体上而言，PCP-Youngs法要优于Hirt-Nichols法。

2.3 速度和加速度

图11和图12所示为30%装载率情况下，中间纵向剖面上Q2点处的速度和加速度的时间历程曲线对比。46%装载率下以及其他各测试点处的速度和加速度历程对比曲线均与之类似。由图中可看出，采用PCP-Youngs法得到的速度和加速度曲线更加光顺，扰动较小，而用Hirt-Nichols法得到的加速度曲线的稳定性则较差，存在较多震荡。这是由于Hirt-Nichols算法的Donor-Acceptor型零阶方法对晃荡自由表面模拟的精度太低，从而直接导致自由表面网格流体的传输计算存在一定的误差，且随着计算中误差的累积，内部流场速度场，特别是加速度场的数值模拟会存在较多的误差与扰动。

图11 Q2处水平速度（30%装载率）Fig.11 Velocity of Q2（filling ratio 30%）

图12 Q2处水平加速度（30%装载率）Fig.12 Acceleration of Q2（filling ratio 30%）

对速度和加速度的时间历程曲线进行两次中值、两次均值滤波，并将速度和加速度曲线与滤波后的曲线相减，得到速度和加速度曲线的数值波动曲线如图13和图14所示。表3所示为速度和加速度波动曲线的标准差。由表3可看出，PCP-Youngs法的数值波动约减小到了Hirt-Nichols法的1/3～1/4。

图13 Q2处水平速度波动（30%装载率）Fig.13 Velocity wave of Q2（filling ratio 30%）

图14 Q2处水平加速度波动（30%装载率）Fig.14 Acceleration wave of Q2（filling ratio 30%）

表3 速度和加速度曲线波动的标准差Tab.3 The standard deviation of velocity wave and acceleration wave

通过上述对比可以看出，与Hirt-Nichols法相比，PCP-Youngs法在自由液面、舱壁压力、舱室内部的速度和加速度预报等方面均有不同程度的改善，特别是有效改善了加速度预报的数值稳定性，为更加准确地预报液舱内部构件的晃荡载荷提供了有力保障。

3 网格敏感性分析

为了验证本文方法的可靠性和稳定性，选择30%装载率工况进行了网格敏感性分析，设定了3种网格大小情况：57×39，114×78和171×117。

选取底边舱斜板上R1处的压力进行对比分析，计算结果如图15和图16所示。由图可见，随着网格尺度的加密，冲击峰值时刻也略微提前，但相差不大。选取液舱中间纵向剖面上Q2处进行速度和加速度对比，其结果如图17和图18所示。由图可见，随着网格的细化，速度和加速度的历程曲线基本一致，只是在峰值处有细微变化。表4所示为不同网格密度下加速度和速度的峰值。由表可见，对于用不同网格密度计算出的速度和加速度，仅相差约1%～3%。由此可见，网格的粗细对于PCP-Youngs法的计算结果影响较小，特别是对速度和加速度的影响更小，其在实际工程计算中不需要采用很精细的网格便可得到较准确的数值模拟结果。

图15 R1处压力网格敏感性分析Fig.15 The pressure cell sensibility analysis of R1

图16 R1处压力网格敏感性分析局部放大图Fig.16 Zoomed picture in pressure peak of R1

图17 Q2处速度网格敏感性分析Fig.17 Cell sensibility analysis for velocity of Q2

图18 Q2处加速度网格敏感性分析Fig.18 Cell sensibility analysis for acceleration of Q2

表4 不同网格密度下速度和加速度的峰值Tab.4 The maximum and minimum value of velocity and acceleration in different cell numbers

本文采用 57×39，114×78和 171×117这 3种网格，用PCP-Youngs法对30%的装载工况进行了15个周期的计算，所耗费的时间分别约为8，40和150 min，可节省的计算资源非常可观。因此，在实际工程应用中，只需采用网格尺度为液舱边界长度2%左右的较粗网格，约10 min即可快速、准确地计算出液舱内部构件位置处的速度和加速度，从而快速确定液舱内部构件的晃荡载荷。

4 结 论

本文针对大型补给舰船液舱结构的晃荡冲击问题，考虑液舱为非矩形截面的复杂边界，以及可能含有内部构件的情况，基于改进的VOF法，采用Youngs法重构自由液面，并结合PCP概念，提出了一种精度更高且数值稳定的处理非矩形复杂舱室边界晃荡载荷的计算方法。分别对不同装载情况下的舱室边界压力、液舱中间纵向剖面处的速度、加速度以及网格敏感性等重要参数进行了对比分析，得到以下结论：

1）通过与模型试验对比，显示本文建立的PCP-Youngs法可有效改进自由表面的模拟效果；通过与实验数据进行对比，发现与Hirt-Nichols法相比，用PCP-Youngs法得到的压力预报相对误差可减小约2%～12%。

2）与Hirt-Nichols法相比，本文建立的PCPYoungs法不仅具有很好的数值稳定性，速度和加速度曲线有较强的周期性且精度更高，而且数值波动也减小到了Hirt-Nichols法的1/3～1/4左右。

3）网格敏感性分析显示，当网格尺度不大于液舱边界长度的2%左右时，数值计算的网格密度对计算结果影响较小，特别是对速度和加速度的时间历程曲线的影响更小，曲线几乎一致，不同网格之间的误差约为1%～3%，在实际工程中的误差许可范围之内。因此，在实际工程应用中不必采用极细的网格便可得到较满意的预报结果。

本文的方法可用于大型补给舰船液舱结构晃荡设计载荷的快速、准确预报。

［1］丁仕风，唐文勇，张圣坤.速度突变情况下液化天然气船液舱内晃荡问题的仿真［J］.上海交通大学学报，2008，42（6）：919-923.DING S F，TANG W Y，ZHANG S K.Simulation of liquid sloshing in the cabin caused by liquefied natural gas ship’s variable speed［J］.Journal of Shanghai Jiao Tong University，2008，42（6）：919-923.

［2］张书谊，段文洋.矩形液舱横荡流体载荷的Fluent数值模拟［J］.中国舰船研究，2011，6（5）：73-77.ZHANG S Y，DUAN W Y.Numerical simulation of sloshing loads on rectangular tank based on fluent［J］.Chinese Journal of Ship Research，2011，6 （5）：73-77.

［3］KKEEFSMAN K M T，FEKEN G，VELDMAN A E P，et al.A volume of fluid based simulation method for wave impact problems［J］.Journal of computational Physics，2005，206（1）：363-393.

［4］沈猛，王刚，唐文勇.基于改进VOF法的棱形液舱液体晃荡分析［J］.中国造船，2009，50（1）：1-9.SHEN M，WANG G，TANG W Y.Liguid sloshing analysis on prismatic tanks based on improved VOF method［J］.Shipbuilding of China，2009，50（1）：1-9.

［5］刘儒勋，王志峰.数值模拟方法和运动界面追踪［M］.合肥：中国科学技术大学出版社，2001.

［6］端木玉，朱仁庆.基于Youngs法的液舱晃荡大幅晃荡数值模拟［J］.船舶力学，2009，13（1）：9-18.DUAN M Y，ZHU R Q.Numerical simulation of violent liquid sloshing in a tank based on Youngs method［J］.Journal of Ship Mechanics，2009，13（1）：9-18.

［7］LIN P Z.A fixed-grid model for simulation of a moving body in free surface flows［J］.Computers and Fluids，2007，36（3）：549-561.

［8］朱仁庆.液体晃荡及其与结构的相互作用［D］.无锡：中国船舶科学研究中心，2002.ZHU R Q.Time domain simulation of liquid sloshing and its interaction with flexible structure［D］.Wuxi：China Ship Scientific Research Center，2002.

［9］端木玉，朱仁庆.流体体积方程的求解方法［J］.江苏科技大学学报：自然科学版，2007，21（2）：10-15.DUAN M Y，ZHU R Q.Method for solving fluid volume equation［J］.Journal of Jiangsu University of Science and Technology：Natural Science Edition，2007，21（2）：10-15.

［10］MIKELIS N E，MILLETR J K，TAYLOR K V.Sloshing in partially filled liquid tank and its effect on ship motions：numerical simulations and experimental verification［J］.The Royal Institution of Naval Architects Transaction，1984，126：267-281.

［11］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.