张 珍,陈明辉,邵 菲,谢永和∗

(1.浙江海洋大学 船舶与海运学院,浙江 舟山316022;2.中国舰船研究设计中心,湖北 武汉430064)

载液货船在航行过程中引起的液舱晃荡问题一直以来都备受关注。强非线性以及随机性使得液舱晃荡一直是流体力学领域极具挑战的研究课题之一。对液舱晃荡引起的晃荡载荷、结构响应以及两者耦合效应的研究是保证船舶航行的重要依据,也为载液船液舱设计提供了参考价值。

目前液舱晃荡的研究主要集中在理论分析,试验研究以及数值模拟3 种途径。最早源于ABRAMSON[1]对线性势流理论方面的总结,随后FALTINSEN[2]采用非线性势流理论对水平激励作用下的矩形液舱晃荡进行了探究。相比理论分析,试验探究可以为液舱晃荡提供直观借鉴,FABRIZIO PISTANI[3]开展晃荡模型试验,通过对比不同箱体晃荡冲击压力,发现棱形箱体斜面设计可有效减缓冲击压力。BRAEUNIG 等[4]对横荡激励下LNG 液舱内的晃荡现象进行了试验和分析。虽然模型试验可以直观展现晃荡剧烈程度,但在部分条件限制下很难完成,且花费的成本较高。随着社会的进步和计算机的广泛运用,数值模拟方法备受青睐。其中有基于自由表面跟踪技术的探究,如:BATTAGLIA[5]通过水平(Level Set)法追踪自由液面对三维矩形液舱晃荡进行探究。SHAMSODDINI[6]基于ISPH 法分析了挡板对于减轻液舱晃荡的影响。还有根据计算网格的研究,如:JIANG[7]等基于Open FOAM 模拟棱形液舱晃荡,ZHAO[8]等基于重叠网格法对部分装载的LNG 液舱进行了研究。

近年来,液舱晃荡问题研究越来越广泛和趋于成熟,然而,大部分学者是将液舱刚性化,但在实际问题中,液舱其实是弹性结构,导致晃荡引起的冲击压力会受到舱壁弹性的影响。随着水弹性力学的发展,陆续有学者开展了弹性液舱晃荡的研究。朱仁庆和刘艳敏[9]将有限体积法、显示积分法、ALE 法相结合对弹性液舱晃荡进行数值模拟探究,不仅考虑了外部激励,还考虑了装载率变化对弹性液舱内晃荡载荷与响应的影响。ZHANG[10]采用MPSFEM 耦合方法对晃荡冲击载荷引起的弹性结构响应进行数值分析,改变材料的杨氏模量探究弹性舱内晃荡现象,分析冲击压力变化和结构在时域频域的动态响应特征。蒋梅荣等[11]基于振动台试验对弹性壁液舱的共振情况做了探究,将液舱置于简谐激励振动下,研究液舱内液体晃动现象。

本工作基于数值模拟方法对刚性和弹性液舱晃荡问题进行充分探究。通过改变外部激励频率参数,模拟刚性和弹性液舱晃荡过程,对比分析因晃荡而产生的对舱壁冲击压力以及现象图,探究外部激励对刚性与弹性液舱晃荡的影响。

1 理论方法

1.1 控制方程

本文主要是求解N-S方程,计算隐式不定常以及不可压缩的粘性流体在流场中的一系列问题。N-S方程如下所示:

1.2 自由液面捕捉方法

本工作应用流体体积(VOF)方法[12],通过追踪每一个网格中水和体积的体积分数来捕捉两相流的相界面。该方法假设网格分辨率足以求解相之间的交界面位置和形状。

交界面的相分布和位置由相体积分数的场来描述,根据相体积分数的值,可以区分网格单元中不同相或流体是否存在,相φ的体积分数定义如下:

其中,w和a分别表示水和空气,本文采用φ＝0.5作为水和空气间的相界面。应用VOF可以实现良好的质量守恒以及提高计算效率。

1.3 FVM 方 法

本工作采用的数值模拟核心方法是有限体积法(FVM)[13],FVM 法需要离散化求解区域,并将其划分为有限大小的离散网格。每个网格节点以特定方法形成围绕该节点的控制体积。因此,其核心是将控制微分方程集成到控制体积中:

其中,等式左边第一项表示控制变量V中特征变量φ的总量随时间变化。n·(ρφu)是沿n方向的流率,n是由于对流流动而在控制体积表面外部的法线方向。n·(Γ·gradφ)表示φ向外部控制量的扩散率。

1.4 流固耦合(FSI)

流固耦合[14]用来描述移动或变形的固体与流体内流场或外流场之间相互作用、相互影响的问题。当流动的流体与固体结构接触时,固体会受到应力与应变的作用,这些力会使结构产生变形。应力使作用于表面的单位面积力的度量。

Von Mises应力[15]是很重要的韧性失效度量,表示偏应力的大小,与平均应力无关:

如果材料处于单轴应力状态,则Von Mises应力等于应力的非零分量。通常,通过加载大致处于单轴应力状态的试样来研究材料屈服。式(9)中,S1,S2,S3分别为第一,二,三主应力。其大概的含义是在一定的变形条件下,当材料的单位体积形状改变的弹性位能(又称弹性形变能)达到某一常数时,材料就屈服;即当单元体的形状改变比能达到一定程度,材料开始屈服。

2 模型验证

为验证数值模拟方法的可靠性,采用西班牙造船厂Navantia建造的138 000 m2的LNG 液舱纵剖面做横摇模拟,该剖面按照1∶50的比例模型化(如图1),液舱的长度为L＝0.90 m,液舱的高度为H＝0.58 m,液舱内部的液体高度h＝0.093 m。在图(1)中,位于液舱左侧固壁上位于自由液面处的一点P点为压力监测点。该液舱模型的运动方程表示为

图1 液仓物理模型Fig.1 Physical model of tank

其中,θmax＝4°,为角位移振幅,周期T＝2.112 s,此为液舱的自然共振周期。

液舱横摇运动的中心选在液舱底部的中心位置,检测位于液舱左侧固壁上自由液面处一点P的压力随液舱横摇运动的变化。为了便于同实验值进行比较,对时间和压力进行了无量纲化处理,压力无量纲化处理公式为P＝P/ρgh,其中,P为P点处的压力值,ρ为液舱内液体的密度,h为液舱内的液体高度;g为重力加速度。

由于本工作采用的数值方法需要进行网格划分,因此在计算前进行了网格独立性验证,从而来确保精度。图2列出了15 s内最大压力随网格数变化的情况。采用了4种网格数来进行验证,从图中可以看出,最大误差在3%左右,因而可以验证网格的无关性。

图2 网格无关性验证Fig.2 Grid independence verification

为验证时间步长的独立性,对不同时间步长的压力历时曲线进行了对比探究。图3分别显示了压力在步长为0.001,0.001 5和0.002 s的历时曲线。从图中可见,最大误差在0.28%左右,因此,验证步长的独立性。而本文采取0.001 s的时间步长来更精准探究液舱晃荡的影响。

图3 P 点压力历时曲线Fig.3 Time history of pressure at the point P

图4展现了液舱纵剖面横摇晃荡前十五个周期的角度和液舱左侧壁面上P点处的压力时历值,包括实验结果与数值模拟结果。从图中可以明显看出液舱左侧固壁上P点的压力在自由表面波拍击壁面时刻发生徒增现象。压力最大值约为1 kPa,压力值误差较小,且数值模拟曲线波动基本与试验曲线一致,结果吻合的很好。从而验证了本工作采用的数值模拟方法在解决液舱晃荡问题方面是相当有效的。

图4 P 点的压力历时曲线Fig.4 Time history of pressure at the point P

3 数值模拟

本工作采用液舱几何尺寸如图5所示,液舱长为1.2 m,高为0.6 m,舱内水深0.12 m,对应充水率为20%,右侧和底部舱壁共设置3个测压点,位置如图1所示。液舱沿x方向做如下运动:其中,A＝0.06 m,为液舱运动幅值,由解析理论得到液舱晃荡一阶共振频率ω＝2.975 s－1,为验证外部激励频率对液舱晃荡的影响,外部激励频率设置如表1所示的3种工况。f为此装载下液舱晃荡的共振频率。

图5 液舱几何模型Fig.5 Geometric model of the tank

表1 计算工况Table 1 Cases for simulation

3.1 刚性液舱晃荡

为了计算方便,本节将液舱简化为二维结构。将舱壁视为刚性。网格采用四面体棱柱层网格。图6给出了不同激励频率下晃荡流体的运动情况。从图6可以看出,激励频率对晃荡流体晃荡特征具有显著的影响。随着激励频率变大,液体冲击左侧舱壁的高度越高,出现碎波现象越来越严重。从图6(h)可见,液面卷入空气比另两种工况严重,而到图6(i)所示时刻时,出现了明显不同,且伴有液体飞溅现象,这种飞溅现象因激励频率变大,而更严重。

图6 晃荡过程现象图Fig.6 Impact process of the tank sloshing

图7～9为3个压力点所监测的压力变化情况。图7可见,P1压力点所监测压力差异较小。在7 s前,3种工况曲线基本重合,7 s后,呈现的双峰特性的第二个峰值附近有些许差异。在12.5 s左右,0.4f状态下,曲线出现明显波动,而波动未在另两种工况下发生。在14 s时,2倍的激励频率作用下,产生的压力比另两种工况产生的压力值要大,约为1 k Pa,另两种工况值则低于1 kPa。0.4f工况下,压力最大值约为2 000 Pa,而2f工况下,最大压力达到近2 250 Pa。

图7 P 1 压力点压力变化历史曲线Fig.7 Time history of pressure at the point P 1

图8为P2压力点所监测压力的对比图。0.4f工况下,最大压力值约1 750 Pa,而2f工况下,最大值约为1 950 Pa。14 s时,曲线有明显差异,2f工况下压力大于其他两种工况。压力曲线仍然呈双峰特性。

图8 P 2 压力点压力变化历史曲线Fig.8 Time history of pressure at the point P 2

图9为P3压力所监测压力对比图,因为本研究涉及两相流,所以压力值出现了负值的情况。从图中可见,2f工况下,液面卷入空气情况较严重。作用P3压力点的最大压力为f工况下,达到1 220 Pa,而2f工况下最大压力低于1 200 Pa。在11.3 s左右,f工况下压力比另两种工况下压力要大,但在13.6 s左右,2f工况下最大压力大于其余两种工况。此外,2f工况下,曲线波动没有另两种工况平滑。

图9 P 3 压力点压力变化历史曲线Fig.9 Time history of pressure at the point P 3

综合图7～9可见,外部激励频率对液舱晃荡有一定的影响,频率越大,晃荡越明显,自由液面出现碎波,液体飞溅现象越严重。在7 s前,3种工况压力历时曲线基本吻合,但在7 s后,曲线开始发生变化,2f工况下,曲线的波动次数比其他两种工况要多,导致曲线不如另两种工况平滑。在14 s附近,P1,P2点压力差异较明显。13.5 s左右,P3点的峰值差异较明显。由于两相流的作用,P3值出现了负值情况,这是由于空气对其影响而产生的。

3.2 弹性液舱晃荡

为探究外部激励频率对弹性液舱晃荡的影响,壳体单元采用铝材质,因为铝的特性相对较柔,可以达到模拟的目的。液舱舱壁厚度为10 mm,具体材料参数如表2所示。图10为20%充液率状态下三维液舱模型图。因为流固耦合是非常复杂的问题,所以在对模型进行网格划分时,固体域和流体域采用不同网格:固体域部分采用四面体网格,流体域部分采用多面体网格。图11为模型的网格划分图。

图10 液舱计算模型Fig.10 Tanks calculation model

图11 模型网格划分图Fig.11 Mesh used in the simulation

表2 材料参数表Table 2 Material parameters

图12为考虑舱壁的弹性影响,在不同激励频率作用下,P1点和P2的压力历时曲线的比较。从图12中可以看出,不同频率作用下,压力时历曲线复杂程度均有所不同。对于弹性液舱,当外部激励频率恰为共振频率时,液舱晃荡更剧烈,产生的抨击压力也比另两种工况要大。当所加的外部激励频率远离共振频率时,不论是增大还是减小,对P1,P2点作用的抨击压力会减小,且前后两个脉冲数值相差不明显。当外部激励频率为0.4f时,压强时历曲线比2f作用下产生的时历曲线更平缓。

图12 P 1 和P 2 点的压力时历曲线对比Fig.12 Comparison for time history of pressure at the point P 1 and P 2

图13是当外部激励频率等于共振频率时,刚性液舱与弹性液舱上的3个压力监测点所监测的压力曲线对比图。图13(a)、(b)分别是P1点和P2的压力历时曲线对比图,P1点在刚性液舱上的最大压力约为2 180 Pa,而在弹性液舱上最大压力约为1 980 Pa。最大压力值减小了9.17%。P2点在刚性液舱上的最大压力约为1 900 Pa,而在弹性液舱上最大压力约为1 750 Pa,压力减小了7.9%。可见舱壁的弹性对晃荡冲击载荷起到了一定的缓和作用,且舱底弹性的缓和作用比舱两壁的缓和作用更明显。液舱晃荡引起的冲击载荷会产生双峰值,弹性液舱中,第二次产生的峰值相较于第一次峰值减小的更明显。

图13(c)是P3压力点的对比图,在前9 s,弹性液舱中的冲击载荷比刚性液舱中的值要小,但在9 s后,弹性液舱中P3压力点比刚性液舱的峰值大,在11.3 s左右,晃荡压力增大约27.5%。15 s内,弹性液舱中最大压力峰值约为1 380 Pa,刚性液舱中最大峰值约为1 230 Pa,增大约10.9%。但在弹性液舱中,晃荡引起的空气卷入现象要改善很多,从图中可见,弹性舱中,P3压力点上无负值出现。且在弹性液舱中,压力历时曲线表现的更光滑,历时曲线形式比起刚性液舱更简单。

图13 当外部激励频率等于共振频率时,刚弹性液舱上不同压力点时历曲线对比Fig.13 Comparison for different pressure points in elastic and rigid tanks when the external excitation frequency is equal to resonance frequency

综合分析图13的3幅对比图可见,液舱的弹性对冲击载荷有一定的缓和作用,且舱底的弹性缓和作用比舱两壁缓和作用要更大。通过对比平均最大压力峰值可见,弹性液舱中冲击载荷比刚性液舱中冲击载荷平均减小了约8.5%。此外,在弹性液舱中,峰值产生的相位相较于刚性液舱中有一定的延时作用。

图14(a)是当外部激励频率等于共振频率作用下,自由液面达到左舱壁最高处时的现象图。对比图6(b)中的第一幅图可见,在弹性液舱中,晃荡引起的自由液面爬壁情况比刚性液舱中液面爬壁高度要低。图6(b)表示在自由液面达到左舱壁最高处时,作用在液舱上的范式等效应力。从图中可见,最大应力为0.106 12 Pa,存在于左舱壁上。因此,在设计液舱结构的过程中,需要对此部分进行加强,防止疲劳引起舱壁结构破坏。

图14 自由液面达到左舱壁最高处现象图以及瞬时Stress Von Mises图Fig.14 Impact process of the free surface arriving to the left wall highest place and the instantaneous Stress Von Mises

4 结 论

建立了刚性与弹性液舱的数值模型,模拟了不同外部激励频率作用下液舱晃荡过程,主要结论如下:

1)采取的数值方法是可靠且准确的,并且该方法可以很好的模拟液舱晃荡过程中的碎波,液体粒子飞溅现象。

2)外部激励频率的变化对液舱晃荡有一定的影响。在刚性液舱中,激励频率越大,产生的碎波,液体飞溅现象越严重,但因晃荡引起的冲击载荷却相差不大。

3)在弹性液舱中,当外部激励频率等于液舱晃荡的共振频率时,作用在舱壁上的平均压力最大。当外部激励频率远离共振频率时,无论是增大还是减小,冲击载荷会减小。

4)液舱的弹性作用对冲击载荷有一定的减缓作用,且舱底的缓和作用比舱两壁更明显。在设计液舱结构时,因为舱两壁的应力比其他地方的应力要大,所以需要注意舱两壁的加厚或材料的选取。