基于势流理论的内孤立波与立管作用数值研究1)

胡英杰 邹 丽,†,2) 孙 哲 金国庆 马鑫宇

* (大连理工大学船舶工程学院,辽宁大连 116024)

† (工业装备结构分析国家重点实验室,辽宁大连 116024)

引言

海洋内波是一种发生在海面以下的波浪,具有波幅大、周期长的特征,内波在海洋中经常以大波幅的内孤立波及波列的形式传播,在世界范围的海域内广泛存在[1-3].在传播过程中,内孤立波的流场会在海洋表面产生明显的特征,因而使用合成孔径雷达(SAR)对海面特征的捕捉成为对其观测的重要依据[4-6].同时,使用声学多普勒流速仪(ADCP)等仪器的海洋实测也在多个区域观测到了内孤立波的存在[7-8].通过遥感和实地观测表明我国南海也是内孤立波频发的海域[9].在海洋资源开发过程中,海洋石油和深海矿产输送立管长期在海上作业,海洋立管在内孤立波流场的作用下会产生剧烈而大幅度的动力响应,对管道的构形、结构强度、作业效率与安全都产生严重影响[10].海洋内孤立波不仅对水下潜器和结构物影响巨大,还会对海洋平台等水面结构物产生影响.内孤立波在远距离传播过程中,会在海洋表面引导出明显的海水流动,因而会造成海洋平台的长距离漂移,对锚泊结构是一个重大考验,影响正常生产作业,造成严重的经济损失[11].在安德曼海作业的一个石油钻井机曾被海洋内孤立波的流场推移了30.48 m,严重影响作业安全[12].

近年来,内孤立波对海洋立管的作用引起了更多学者的关注.文献[13-15]针对内孤立波对海洋立管的作用问题开展了系统深入的研究,通过KdV 和mKdV 方程计算了内孤立波的流场特征,采用动力有限元模型得到了顶张力立管的动力响应,并讨论了内波参数对管道响应的影响.文献[16-18]采用试验和数值模拟方法研究了内孤立波与海流共同作用下管道与海洋平台桩柱的动力响应问题,获得了不同入射角度及不同波浪参数条件下的管道变形特征.Wang等[19]采用非线性有限元模型计算了弱约束深海矿产输送管道在内孤立波作用下的变形与运动特性,讨论了内孤立波不同入射角度及水面平台运动对管道的影响.Wang等[20]采用试验方法研究了内孤立波对小尺度桩柱的作用力分布特征,讨论了不同水深条件下桩柱的受力情况,分析了水动力系数选取与KC数的关系.Fan等[21]基于欧拉-伯努利梁理论,考虑了海流、表面波和内孤立波的共同作用,发现内孤立波的存在明显增大了立管的变形包络,并且分层流体的密度差也会对立管变形产生影响.李朝玮等[22]采用准静态分析法研究了非均匀海流、表面波和内孤立波作用下隔水管的变形,发现内孤立波对立管的横向位移和应力都有不可忽视的影响且表面波和内波相位角分别在12°和0°时对立管的影响最大.Zhang等[23]采用数值模拟和试验方法研究了平台运动和内孤立波共同作用下顶张力立管的涡激振动特性,发现内孤立波的波幅越大、平台升沉运动的频率越高,立管的涡激振动越明显.Wang等[24]在分层流水槽中开展模型试验分析了斜向入射的内孤立波的波流结构特征及其对水下细长结构的载荷作用特性.文献[25-28]采用莫里森经验公式、模态分解和回归分析等方法对内孤立波对柱体及管道的动力载荷进行了深入分析.

内孤立波流场速度和加速度信息的准确获取是立管载荷确定的关键,进而影响整个计算模型的准确性.对于大尺度的海洋内孤立波,无法开展实验获得流场信息,同时求解NS 方程的CFD 方法计算量较大.目前对大尺度内孤立波与立管作用的研究中,内孤立波流场信息多是基于解析模型近似给定的.本文工作基于势流理论采用多域边界单元法建立数值模型求解了内孤立波的流场特征,较为高效又能实时准确地计算流场信息.依据流场信息采用莫里森方程给定立管载荷,采用动力学有限元方法计算了顶张力立管和弱约束立管的动力响应,以期为海洋立管设计在考虑内孤立波作用时提供参考.

1 计算模型

1.1 分层流体运动势流计算模型

内孤立波与立管作用如图1 所示,两层流体的深度分别为h1,h2;密度分别为ρ1,ρ2.下凹型内孤立波位于密度界面分层处.

图1 内孤立波与立管作用示意图Fig.1 Schematic diagram of the interaction between internal solitary wave and the riser

两层流体都假定为无旋,无黏性,不可压缩的,因此两层流体中存在流动势函数φk满足拉普拉斯方程

在分层流体界面处的运动学边界条件为上下层流体沿界面法线方向的速度连续,因此有

其中η2为内孤立波波面,n为计算区域Ω1,Ω2在界面处的外法线方向矢量.

在分层流体界面处的动力学边界条件为上下层流体区域在界面处的压力连续,因此有

其中g 为重力加速度.海底满足不可穿透条件

由于内孤立波传播时在海洋表面产生的波动较小,为了提高计算效率,本文对上表面采用刚盖假定,因此上表面的边界条件为

通过格林第二定理,可将上下两层流体区域Ω1,Ω2的控制方程拉普拉斯方程转化为两个边界上的曲线积分方程

其中r为源点,q为场点,G(r,q)=-1/(2π)ln|r-q|为格林函数,c(r)=1-α(r)/(2π),α(r)为场点r处的单元夹角,Γk为上下层区域边界.

由于内孤立波位于密度分层的流场中,所以计算过程需要涉及两个流体区域,进而要将单独区域的边界单元法扩展到多个区域[29-30].对于计算区域Ω1,Ω2的边界分别采用n个线单元进行离散,从而积分方程 (6) 可以写成如下的形式

根据两个区域的关系,可以将边界点分为公共边界点和自身边界点,从而将积分方程进一步处理为

其中,上标c和s分别指公共边界点和自身边界点.

波面位置和速度势的更新格式分别根据运动学和动力学边界条件给定

内孤立波初始波形的给定依据强非线性的MCC 模型,实现对大波幅内孤立波的准确模拟.

1.2 管道动力响应模型

在内孤立波由远方向立管传播过程中,立管处的流场速度及加速度是实时变化的,因此不能用简单的均匀流或剪切流静力假定求解,需要采用动力学进行求解.

海洋立管的动力学控制方程

其中,M,C,K分别为质量阵、刚度阵和阻尼阵,x(t),分别为管道的变形位移、速度和加速度.Q(t)为载荷向量,根据内孤立波流场的速度、加速度以由莫里森方程给定.

在计算中需要考虑管道轴向张力的影响,因此管道的刚度阵由两个部分组成,分别是弹性刚度阵KE和几何刚度阵KG,即

其中L为管道单元的长度.

阻尼阵采用瑞利阻尼,质量阵与刚度阵的线性组合形式

其中 α 和 β 为比例系数,根据模态频率给定,分别为

ω1和ω2分别为管道的一阶和二阶频率,ζ为管道的阻尼比.

采用Newmark-β法求解动力学控制方程 (11),获得管道的动力响应.

2 内孤立波载荷

在通过对边界积分方程的求解获得区域边界上的速度势及其法向导数值后,流场内部速度势可以用方程 (15) 表达.对积分方程 (15) 进行微分,即可将速度势的积分方程转化为关于流场速度的积分方程 (16),对该方程进行离散求解即可获得内孤立波流场内部的速度分布特征,即

在进行计算求解时,当内部的待求点与边界节点接近时,格林函数会表现出近奇异特征,导致计算精度降低甚至使计算失效.本文的研究中采用了单元子分法[31-32]来解决这一问题,得到了良好的计算结果.为了验证计算模型的准确性,与文献中的内孤立波水平速度垂向分布试验结果[33]进行了对比.试验在分层流体水槽中开展,上层流体深度h1=15 cm,下层深度h2=62 cm,密度分别为ρ1=0.999 kg/m3,ρ2=1022 kg/m3,内孤立波无量纲波幅为a/h1=0.22.计算结果与分层流体内孤立波传播试验流场的测量结果吻合良好如图2 所示.计算得到的内孤立波流场速度分布如图3 所示,计算结果表明内孤立波流场速度具有明显的剪切特征,上层流体的质点速度与内孤立波的传播方向相同,而下层流体则与之相反,且上层流体质点速度的绝对值要明显大于下层流体.图4 为内孤立波流场的水平速度在上层流体区域z/h1=0.2 处沿水平方向的分布情况.可以发现水平速度的最大值位于波谷位置(x/h1=25),且随着与波谷距离的增大而迅速衰减.图5为z/h1=0.2 处流体质点加速度在波长方向的分布情况,加速度的最大值位于波谷两侧,波谷位置处的加速度为零.

图2 水平速度垂向分布与文献[33]结果对比Fig.2 Comparison between the calculation result about the vertical distribution of horizontal velocity with the result of Ref.[33]

图3 内孤立波流场速度分布Fig.3 Velocity distribution of internal solitary wave flow field

图4 内孤立波水平速度沿水平方向的分布Fig.4 Horizontal velocity distribution of internal solitary wave

图5 内孤立波水平加速度分布Fig.5 Horizontal acceleration distribution of internal solitary wave

由于深海矿产输送管道具有大长细比特征,因此采用经典的莫里森公式作为桥梁将非线性内孤立波数值计算模型和管道响应的计算模型连接起来,从而实现内孤立波对管道作用的准确求解.

由于内孤立波对管道的作用过程中管道变形较大,所以要考虑管道的位移以及管道倾角的影响,因此内孤立波对管道作用力的莫里森公式写成如下形式

其中,Cm为惯性力系数,本文计算取为2.0;Cd为拖曳力系数,取为1.2;Un为垂直于轴向的内孤立波流场速度,为对应的加速度;ρ为海水密度,x˙ 和x¨ 分别为管道运动的速度和加速度.根据前面内孤立波流场速度和加速度的计算结果及莫里森方程可知,当管道的直径相对较小时,对管道响应起主要作用的载荷为内孤立波对管道的拖曳力,惯性力的影响相对较小.

3 管道动力响应

3.1 顶张力立管动力响应

在海洋石油开采等领域,顶张力立管得到了广泛的应用,本文计算了内孤立波对顶张力立管的作用特性.为了验证本文计算模型的准确性,将计算结果与文献[34]结果进行了对比.上下层流体深度分分别为60 m 和412 m,密度分别为ρ1=1025 kg/m3,ρ2=1028 kg/m3,管道内外径分别为0.26 m 和0.21 m,内孤立波波幅为75 m.对比结果如图6 所示,可以发现两者吻合较好.由于本文内孤立波流场的给定方法与文献有一定区别,所以动力响应结果稍有偏差,管道顺流向最大位移处的偏差约为3.34%.

图6 管道流向最大位移与文献[34]计算结果对比Fig.6 Comparison of the maximum displacement of the pipe in flow direction with the calculation results in the Ref.[34]

本文计算的顶张力立管的几何参数和物理参数见表1.为了讨论顶张力大小对管道动力响应的影响,顶张力T分别取为400 kN,500 kN,600 kN,700 kN,800 kN.

表1 顶张力立管参数Table 1 Parameters of the top tension riser

计算中采用的分层流体参数见表2.为了研究内孤立波参数对管道动力响应的影响,内孤立波的波幅分别给定为40 m,50 m,60 m,70 m,80 m.

表2 流体分层参数Table 2 Parameters of stratified fluids

不同内孤立波波幅时管道的变形特征如图7 所示,可以发现立管的动力响应对波高敏感,随着波高的增加,管道的流向位移明显增大.同时,由于上层流体速度明显大于下层,且在所研究问题中拖曳力远大于惯性力,因此管道顺流向的最大位移发生在上层区域.

图7 不同内孤立波入射波幅时顶张力立管的流向变形特征Fig.7 Characteristics deformation of top tension riser in flow direction with different incident amplitudes of internal solitary waves

当内孤立波波幅为60 m 时不同顶张力条件下的立管变形特征如图8 所示,可以发现顶张力对管道的流向变形具有重要影响.在相同内孤立波参数的条件下,施加的顶张力越大,管道的变形越小.图9展示了不同张力时管道流向最大变形的变化趋势.这是由于随着顶张力T的增加,作用在管道单元上的有效张力Te也逐渐增大,通过方程 (13) 可知,管道的几何刚度阵KG值将随之增大,进而提高管道的整体刚度.因此,在相同参数的内孤立波作用下,顶张力越大,管道的顺流向变形越小.

图8 不同顶张力时管道的变形特征Fig.8 Deformation characteristics of the riser under different top tensions

图9 管道最大流向位移随顶张力大小的变化关系Fig.9 The relationship between the maximum displacement of the riser in flow direction and the value of the top tension

图10 展示了z=-40 m 位置处管道的应力在不同内孤立波波幅条件下的变化特征.管道应力随着波幅的增大而显著增大,同时波幅较小的内孤立波作用时管道的最大应力时刻出现滞后.根据分层流体K d V 理论,内孤立波的波速可以表示为,其中c0和c1是与分层流体厚度和密度有关的量,因此当分层参数确定之后,内孤立波的波速与波幅密切相关,波幅越大的内孤立波具有更大的波速.因此,当波幅较小时,内孤立波传播到管道位置处的时间显著增加,管道的最大应力出现的时刻明显滞后.

图10 不同波幅内孤立波作用下管道的应力变化Fig.10 Stress changes of the riser under the action of solitary waves with different amplitudes

3.2 弱约束立管动力响应

本文讨论的弱约束立管是指上端与水面设备刚固连接,下端为自由端的海洋立管,一般在深海采矿等领域得到应用.计算所采用的管道几何和物理参数如表3 所示.

表3 弱约束立管参数Table 3 Parameters of the weakly constrained riser

为了讨论内孤立波参数对弱约束立管的作用规律,计算了不同波幅内孤立波作用下的立管动力响应如图11 所示.从图中可以发现内孤立波波幅的变化对立管变形作用明显,立管的流向位移随初始内孤立波波幅的增大而显著增大.

图11 不同波幅内孤立波作用下弱约束立管的变形特征Fig.11 Characteristics deformation of weakly constrained riser in flow direction with different incident amplitudes of internal solitary waves

对于深海矿产输送管道,输送矿物浓度的变化会导致内部流体密度的变化,因此本文讨论了内部流体不同密度对管道响应特征的影响.不同内部流体浓度条件下管道变形特征如图12 所示.从图中发现,管道内部流体密度对弱约束管道的流向变形影响较小.

图12 内部流体密度对弱约束管道位移的影响Fig.12 Influence of inner fluid densities on displacement of weakly constrained riser

4 结论

本文基于势流理论采用多域边界单元法建立了内孤立波流场的计算模型,内孤立波波面采用全非线性动力学和运动学边界条件,准确获得了内孤立波流场的速度分布特征.上下层的流体质点速度方向相反,对于下凹型内孤立波,上层流体速度与波的传播方向相同,下层速度与之相反,且上层流体质点速度明显大于下层.

根据计算获得的内孤立波流场信息,采用莫里森方程确定了内孤立波对海洋立管的载荷输入,计算了顶张力立管和弱约束立管在内孤立波作用下的动力响应特征.管道的顺流向位移和管道应力随内孤立波波幅的增大而显著,同时顶张力通过改变几何刚度阵而对管道的动力响应产生明显影响.由于上层流体速度明显大于下层,且在所研究问题中拖曳力远大于惯性力,因此管道顺流向的最大位移发生在上层区域.同时计算了不同波幅和不同内部流体密度条件下弱约束管道的动力响应,发现内孤立波波幅对弱约束管道的动力响应影响显著,而管道内部流体密度的影响则相对较小.

本文中内孤立波对立管的作用载荷根据流场信息采用莫里森方程给定,没有考虑管道的涡激振动响应特征.为了更加全面分析管道在内孤立波作用下的动力响应规律,需要在今后的研究中进一步改进,为实际工程中的海洋立管设计和安全性评估提供更加充分的依据.