海底物体运动生成内波问题研究

李林杰, 赵彬彬, 段文洋

(哈尔滨工程大学 船舶工程学院,黑龙江 哈尔滨 150001)

海洋内波是指在海水稳定层化的海洋中产生的、最大振幅出现在海洋内部的波动。海洋内波对水下结构物的安全存在巨大威胁,如我国372号潜艇于2014年在远航期间因为遭遇内波而意外下沉,险些造成事故。因此,对内波的生成及其特性进行研究具有重要的意义。

对于物体运动激发内波问题,国内外学者都做了很多研究,可以分为物理实验和数值模拟2个方面。

实验室中的内波模拟以其成本低、可重复性强、容易实现以及容易获得内波相关特征参数的特性受到了更多研究人员的青睐。Schooley等[1]首次实现了关于内波的实验,该实验是针对弱线性分层流体中自驱动物体运动产生的紊流尾迹的研究。Arntsen[2]在分层实验水槽中系统全面地进行了分层流体中圆柱绕流实验,得到不同密度比、不同深度、不同Froude数下圆柱的阻力系数以及升力系数,以及产生的内界面波波高和波长的无量纲数值。马晖扬等[3]通过物理实验,研究了球体在密度突跃层附近运动时导致的自由面的特定变形,讨论了物体的位置和运动速度等因素对波形形态的影响。李万鹏等[4]对分层流体中螺旋桨效应激发内波进行了实验分析。陈科等[5]对1个球体和2个不同长径比细长体在拖曳运动下激发内波的时空特性进行了系列实验。高德宝等[6]在大型分层流水池中开展了拖曳体内波尾迹水下/水面特征试验研究,釆用电导率仪阵列测量了水下波高场及其随时间变化,获取了内波诱导的水面流场,并对内波尾迹水下/水面特征的相关性进行了分析。

相比于物理实验,基于内波的基本方程,利用计算机实现数值模拟内孤立波的生成、传播以及演化等相关特性更加容易。罗恒[7]进行了2层流体中水下运动潜艇激发表面尾迹及其与随机海面相互作用的数值模拟。段菲[8]研究了2层流体中SUBOFF潜艇在不同潜深位置处做非定常运动对内波波动过程的影响。黄璐[9]计算了SUBOFF潜艇在不同工况下近水面航行时产生的 kelvin 波,分析了不同海况以及姿态调整对自由表面波的影响。黄风来等[10-11]分析了连续密度分层流体中JoubertBB2潜艇在自由表面兴波过程。胡开业等[12]研究了在含跃层的分层流中,运动椭球的速度、潜深以及跃层强度对水面兴波的影响。何广华等[13-18]研究了SUBOFF潜艇在近水面航行时的兴波尾迹、兴波阻力、深水中航行的兴波尾迹、兴波阻力,以及附体对潜艇水动力特性的影响。除了定常运动,于祥等[19]研究了潜体的加减速运动对尾迹特性的影响。

然而,上述研究都没有考虑物体在海底运动时生成内波的情况。Grue等[20]对物体在海底运动激发内波问题进行了数值模拟,发现海底匀速运动的物体会在其前方激发内孤立波,随着运动时间增长,激发内孤立波的数目也变多。但是Grue等[20]没有开展物理实验,也没有经过其他数值模型的验证;此外,海底物体形状对内波产生的影响,尚未开展研究;Grue等[20]研究中上层液体顶部边界条件,用的刚盖近似来代替实际海洋中的这种自由面是否合理,有待进一步研究。针对上述的3个问题,本文给出了数值计算模型,介绍了控制方程、边界条件以及VOF方法监测内波界面和自由表面;给出了本文的数值计算结果以及分析。

1 数值计算模型

本文基于商业软件STAR-CCM+,对海底物体运动引发内波进行数值模拟研究,所使用的控制方程、流体域体积(volume of fluid,VOF)方法处理内波界面及边界条件将依次进行介绍。

本文研究的是分层流体中一个物体在海底运动激发内孤立波问题(如图1所示),上层流体厚度较大为h2,下层流体厚度较小为h1,一个物体(以半椭圆形物体为例)以一定速度在海底运动,物体的运动会对内波界面形成扰动,最终生成内孤立波。当上层液体表面边界条件设置为壁面时是刚盖假设,当在上层液体上方加入一层空气是自由液面,示意图如图1。

图1 海底物体运动生成内波示意

1.1 控制方程

STAR-CCM+可以选择层流、湍流、无粘的求解器,本文选择的是无粘求解器,对应的是求解无粘的欧拉方程,对应每一层流体其控制方程为:

(1)

(2)

(3)

式中:u、w分别代表水平、垂直方向的流体质点速度分量;p为压力;ρ为每一层流体的密度;g=9.81m/s2为重力加速度。

1.2 VOF方法处理内波界面

内孤立波形成于流体间的分层界面,可通过VOF法进行界面捕捉。VOF方法是通过计算每一个网格内每种流体所占的体积分数来处理界面问题,目前VOF在有限体积法中得到了广泛的应用[21]。

假设r(k)为第k相在某一单元网格所占的比例,根据质量守恒可得:

(4)

体积分数可以通过每一相的体积与总体积求出:

(5)

本文中当上层液体顶部边界条件是自由表面时,有三相,k=1、2、3,通过分别监测空气、下层液体体积分数为50%的等值面的垂向位置,即可知道自由表面、内波界面的起伏。当上层液体顶部边界条件是刚盖假设时,同理可得内界面的起伏变化。

1.3 边界条件

对于数值模拟物体的运动兴波问题,可以采用重叠网格,使物体本身运动,但本文要模拟的时间较长并且物体是匀速平移运动,若采用这种方式,则需要较长的计算域,所以本文采用的是相对运动的方法,让物体本身不动,在速度入口处给定液体一个反向速度,从而可以减少计算域长度。

入口设置为速度入口边界条件,出口设置为自由流出口边界条件,底部和顶部以及物体表面都设置为壁面边界条件,当考虑自由面效应时,在上层液面上方会增加一层空气,从而使上层液体表面成为自由液面,如图2所示。

图2 边界条件设置示意

2 数值计算及分析

首先使用Grue等[20]的算例进行对比验证,然后研究了海底运动物体的形状对生成内孤立波的影响,最后分析了自由表面近似为刚盖假设的合理性。

2.1 半椭圆形物体在海底运动产生内孤立波

Grue等[20]考虑的是一个二维情况下2层流体中半椭圆形物体在海底运动生成内孤立波问题(如图3所示):上层流体厚度h2=0.12 m,密度ρ2=787.3 kg/m3,下层流体厚度h1=0.03 m,密度ρ1=1 000 kg/m3,上下层厚度比为h2/h1=4,密度比是ρ2/ρ1=0.787 3,一个半椭圆形物体以一定速度U在海底运动,物体的运动会对内波界面形成扰动,最终生成内孤立波。

图3 海底物体运动生成内波示意

物体的运动速度是U=1.1c0,其中特征波速c0为:

(6)

为了验证本文数值计算结果的准确性,首先需要进行数值结果的收敛性分析。本文进行了3种不同数值加密策略的收敛性分析,在基本网格尺寸的基础上再在物体以及内波界面附件进行网格加密,每种策略的网格加密方式都一样,表1中Δx表示的是基本网格尺寸,Δt表示时间步长。

表1 收敛性分析

t(g/h1)1/2=1 080时刻下,3种数值加密策略的海底物体运动生成内孤立波的内波波面空间分布计算结果如图4所示。从图4可以看出,数值策略2和数值策略3得到的结果一样,而数值策略1得到的结果与其他数值策略的结果相差较大,说明本文采用的Δx=0.02 m,Δt=0.003 s可以得到收敛的数值结果。

图4 内波波面空间分布结果收敛性分析

本文将t(g/h1)1/2=1 080时刻下,STAR-CCM+收敛的海底物体运动生成内孤立波的内波波面空间分布计算结果与Grue等[20]的结果进行比较,结果如图5所示。

图5 数值收敛结果与Grue等[20]的结果对比

从图5可以看出,本文使用STAR-CCM+计算得到的海底物体运动生成内孤立波的内波波面空间分布计算结果与Grue等[20]的结果吻合得很好,说明本文的数值计算结果是准确的,同时也验证了Grue等[20]的结果。

图6展示了海底物体以一定速度运动,能够在内波界面激发内孤立波,随着时间的推移,产生的内孤立波的数量越来越多的过程。

图6 不同时刻半椭圆形物体运动激发的内孤立波

2.2 海底运动物体的形状对内波生成影响

为考虑海底运动物体的形状对内孤立波生成的影响,基于Grue等[20]的算例参数,本文使用矩形、半椭圆形、三角形物体来进行海底物体运动生成内孤立波的数值模拟,各个形状的物体都是左右对称的且长和高的尺寸都一样,分别是20h1和0.1h1(h1为下层液体厚度,算例中的具体值是0.03 m),如图7所示。为了更直观地比较矩形、半椭圆形、三角形物体产生的内孤立波的波幅大小以及内孤立波数量关系,可见表2。

表2 不同形状物体运动激发的内孤立波波幅和内孤立波数量

图7 运动物体的形状示意

图8为矩形、半椭圆形、三角形这3种形状的物体在海底运动到t(g/h1)1/2=1860时刻下生成的内孤立波波面空间分布的数值计算结果。从图8可以看出,物体海底在运动,可以在物体前方激发一系列内孤立波,在t(g/h1)1/2=1860时刻,矩形物体前方激发了5个内孤立波,内波波幅为1.01h1;半椭圆形物体前方正在激发第4个内孤立波,内波波幅为0.88h1;半椭圆形物体前方激发了3个内孤立波,内波波幅为0.83h1。

图8 t(g/h1)1/2=1860时刻下不同形状物体在海底运动激发的内孤立波

从图8和表2可以发现,矩形、半椭圆形、三角形这3种形状的物体在海底运动,随着时间的推移,都能在物体前方激发一系列内孤立波,当物体的长宽比确定时,在运动了相同的时间后,面积最大的矩形物体在海底运动激发的内孤立波波幅最大,且矩形物体运动激发的内孤立波数量最多;面积最小的三角形物体在海底运动激发的内孤立波幅最小,且三角形物体运动激发的内孤立波数量最少。

2.3 自由面效应的影响

Grue等[20]考虑的半椭圆形物体在海底运动产生内孤立波算例中,上层液体密度为ρ2=787.3 kg/m3,下层液体密度为ρ1=1 000 kg/m3,上下层液体的密度比是ρ2/ρ1=0.787 3,上层液体表面边界条件使用的是壁面边界条件,也就是采用了刚盖假设,但是Zhao等[22]发现对于2层流体间内孤立波生成问题,只有当上下层密度比ρ2/ρ1接近于1时,刚盖假设才是合理的,否则应该考虑自由面效应的影响。

因此,对于半椭圆形物体运动生成内孤立波的算例,本文做了一组保留自由面的对照组,即在上层液体上方增加层空气,从而使上层液体表面成为自由表面,在上下层液体密度比为0.787,t(g/h1)1/2=1 920时刻下的数值模拟结果如图9所示。

图9 当 ρ2/ρ1=0.787 3时,2种边界条件下半椭圆形物体运动激发的内孤立波

从图9可以看出,当上下层液体密度比为0.787时,海底半椭圆形物体运动到t(g/h1)1/2=1 920时刻,在刚盖假设条件下半椭圆形物体运动激发的内孤立波波幅比在自由表面边界条件下半椭圆形物体运动激发的内孤立波波幅更大,内孤立波传播的距离也更远,也就是自由面效应会减少海底物体运动激发的内孤立波的波幅和波速。

为了进一步研究上下层液体的密度比接近1时刚盖假设的适用情况,本文计算了在上下层液体密度比为0.95时,半椭圆形物体在海底运动同样运动到t(g/h1)1/2=1 920时刻,2种边界条件下物体运动激发的内孤立波波面空间分布结果,如图10所示。

图10 当 ρ2/ρ1=0.95时,2种边界条件下半椭圆形物体运动激发的内孤立波

从图10可以看出,在上下层液体密度比为0.95时,海底半椭圆形物体运动到t(g/h1)1/2=1 920时刻,刚盖假设和自由表面这2种上层液体顶部边界条件下,物体运动激发的内孤立波波形的空间分布结果几乎重合,所以当上下层液体的密度比接近1时,自由表面可以近似为刚盖。

此外,结合图9和图10可以发现,在t(g/h1)1/2=1 920时刻,当上下层液体的密度比为0.787时,半椭圆形物体前方正在激发第4个内孤立波,当上、下层液体的密度比为0.95时,半椭圆形物体前方正在激发第2个内孤立波,2种情况下的内孤立波数量相差1倍。因此,上下层液体的密度比越接近于1时,在同一时刻,海底物体运动激发的内孤立波数量越少。

3 结论

1)本文使用STAR-CCM+验证了在刚盖假设条件下,海底半椭圆形物体运动激发内孤立波,Grue等的数值模型计算结果的准确性。

2)矩形、半椭圆形、三角形物体在海底运动,都能够激发内孤立波。在运动了相同的时间后,矩形物体运动激发的内波波幅最大,内孤立波数量最多,三角形物体运动激发的内孤立波波幅最小,数量最少。

3)如果上、下层液体的密度比是0.78这样不接近1的情况,对于海底物体运动激发内孤立波问题,不可以将自由表面近似为刚盖,否则,刚盖假设条件下内孤立波的波幅和波速比自由表面条件下内孤立波的波幅和波速更大。