ABSTRACT

Based on the Hamilton principle, a three-dimensional finite element dynamic model is established in this paper for investigating the vortex-induced vibration (VIV) of a flexible riser in deep water in kilometers. The model includes the fluid flowing within the riser, and the wake outside the riser is modeled by the van der Pol wake oscillators which are distributed along the riser length and coupled with the lateral accelerations of the riser. In addition to setting an entering speed, the upper end of the riser can be in movements with the platform, while the lower end can either move and swing with the additional weight (such as blowout preventer, BOP), or be connected to the underwater wellhead. Using the proposed models, the VIVs of therisersbared or partially wearing with buoyant blocks were simulated in the currents caused by tide, wind and waves. The effects of the entry speed, BOP weight and water depth on the deviation of the lower end of the riser, the maximum tensile stress at the upper end and the maximum bending stress throughout the riser are discussed especially. This study deals with a typical internal and external fluid-solid coupled vibration problem, involving in the complex initial-boundary values, and aims to providing an effective tool in calculation and analysis for the subsequent related researches.

Keywords:Deep Water, Flexible Riser, Vortex-Induced Vibration, 3D-Model, Van der Pol Wake Model, Fluid-Structure Interaction

深水立管入水过程中涡激耦合振动特性分析

张颖，朱卫平^*

上海大学，力学与工程科学学院，上海市应用数学和力学研究所，上海市能源工程力学重点实验室， 上海

收稿日期：2020年3月20日；录用日期：2020年4月23日；发布日期：2020年6月15日

摘 要

本文基于哈密顿原理，建立了一种柔性立管涡激振动(VIV)的三维有限元动力学模型。该模型考虑了管内流体的流动，而管外尾迹(wake)由一组沿立管变化的并与立管侧向加速度相耦合的范德波尔振子模型(van der Pol oscillator wake model)加以模拟，水深以千米计。立管上端可以设置下入速度，可以随平台移动和摆动，下端则可以悬挂重物(如防喷器)和自由运动，也可以与水下井口连接。基于该模型，本文模拟了在一定的潮流、风流和波浪条件下，带或不带浮力块的隔水管下入过程中的动力学响应，重点讨论了悬挂重量、进入速度、入水深度对隔水管下端落点偏离、上端最大拉应力及全管最大弯曲应力的影响。本文处理的是一种典型的内外流固耦合振问题和复杂初边值问题，为后续有关研究提供了有效的计算分析工具。

关键词 :深水，柔性隔水管，涡激振动，三维模型，范德波尔尾迹模型，流固耦合

Copyright © 2020 by author(s), Yangtze University and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

石油作为当今社会最常用的能源之一，也是经济和社会稳定发展的重要战略资源。陆地及浅海石油资源的勘探已经日渐成熟，然而对于更深水域石油资源的争夺与勘探日趋变为各国的战略需求。而当石油勘探转移到海上领域，就会利用到动态定位平台。动态定位平台可以工作的最大水深在理论上是不受限制的，但实际上存在两个限制：船用立管的长度限制，张紧器和升沉补偿器的性能限制。目前，地球上拥有丰富的海上石油和天然气资源，并且可能大多数都在深水之下，这驱使勘探工作进入到深海甚至超深水领域，一些研究已拓展到3000 m的范围 [1]。深水钻探立管是连接在海底井口和浮动钻探平台之间的重要设备。因此，迫切需要对海洋深水和超深水钻井设备进行开发 [2]。而在开发进程中，如果载荷计算不准确，强度设计不合理，将给深水钻井作业带来安全隐患。此前在墨西哥湾和巴西海的钻井过程中，就因为进料塔的故障造成了巨大的经济损失和环境破坏 [3] [4]。

作为深海石油勘探必不可少的设备(如图1)，钻井立管是从海床井口到钻探平台的主要连接通道，与泥浆的控制出口相连，在深水中承受多种载荷，包括风，浪，潮和尾迹的流体动力，泥浆和管道的重力和浮力，以及平台张紧器施加的张力，立管内外压差和上下温差等，这些状况都与水深有关。在深水中，比如水深达3000 m，立管的纵横比(L/D)将高达6000，是典型的三维柔性液固相互作用的结构。当水流流过垂直立管时，通常有涡从立管表面频繁脱落，从而形成涡流。涡流激励立管在横流(cross-flow, CF)方向振动，而来流则引起立管在流线(in-line, IL)方向的运动。涡激振动(vortex-induced vibration, VIV)是一种典型的自激振动，如果涡脱落频率接近立管的固有频率之一，则会发生频率锁定现象 [5]，这将加速结构的疲劳损伤，可能导致无法估量的损失。

计算流体动力学(CFD)直接数值模拟(DNS) [6] [7] [8] [9] [10] 被认为是预测细长立管VIV现象的理想方法。但是，由于三维DNS对计算资源的需求太高，从而限制了对在实际雷诺数和实际纵横比的情况下油气深水勘探立管的模拟，导致这种贴合实际工况条件下的模拟至今仍然很少。为了满足当前的要求，本文使用van der Pol振子方程(范德波尔方程)来模拟VIV在结构上的涡激效应，这可以大大节省计算需求并能捕获涡激效应的主要特征。利用范德波尔振子来模拟涡激振动对结构的影响已经被前人广泛地讨论 [5] [11] [12]，如Wang等人 [13] 为研究内部输运流体的柔性管道涡激振动引起的屈曲，发展了一种三维计算模型和伽辽金解，其中与涡动力相关的非定常水动力用两个范德波尔振子进行模拟，立管为两端铰支，L/D = 200。然而，当轴立管直径为非均匀，部分安装浮力块，以及更复杂的载荷和边界条件时，例如：在受到剪切流和非均匀压力，以及立管因海上平台漂浮而振动，下端悬重摇摆和向下运动等情况下，伽辽金方法已难以解决。因此在本文中，我们首先基于哈密顿原理建立相应的动力学方程，然后采用有限元法(FEM)离散方程，将涡激振子也离散到各结点上，从而让复杂的载荷和边界条件变得易于处理，便于模拟和分析局部带有浮力块的立管在下入过程中含涡激耦合效应的力学行为。本文提出的模型和解法具有模拟纵横比大于6000的立管对所有上述载荷和边界条件的VIV响应的能力，为以后的海洋钻探系统研究和应用提供了强大的计算平台。

图1. 作业立管的结构、载荷及坐标系意图

2. 模型与求解

2.1. 立管振动模型

本文分析模型中，定义水平面上为y和z轴，垂直方向上为x轴(如图1)，根据哈密顿原理，考虑管内流体压力、密度、流速和温差的动力学方程可表示为

$\begin{array}{l}{\displaystyle {\int }_0^L[\left(m_r+m_f+m_a\right)\left({\ddot{u}}_x\delta u_x+{\ddot{u}}_y\delta u_y+{\ddot{u}}_z\delta u_z\right)}+J_x{\ddot{\theta }}_x\delta {\theta }_x+\cdots \\ +2u_i{\rho }_f{A}_i\left({\stackrel{\dot{}}{u^{\prime }}}_y\delta u_y+{\stackrel{\dot{}}{u^{\prime }}}_z\delta u_z\right)+c\left({\stackrel{\dot{}}{u}}_x\delta u_x+{\stackrel{\dot{}}{u}}_y\delta u_y+{\stackrel{\dot{}}{u}}_z\delta u_z+R_r{\stackrel{\dot{}}{\theta }}_x\delta {\theta }_x\right)+\cdots \\ \text{+}\left[\left(u_i^2{\rho }_f+p_i\right)A_i-T_0\right]\left({u^{″}}_y\delta u_y+{u^{″}}_z\delta u_z\right)+E{A}_r{\left({\epsilon }_x^0+1\right)}^2{u^{\prime }}_x\delta {u^{\prime }}_x+\cdots \\ +G{I}_x{{\theta }^{\prime }}_x\delta {{\theta }^{\prime }}_{x}EI\left({u^{″}}_y\delta {u^{″}}_y+{u^{″}}_z\delta {u^{″}}_z\right)-\left(f_x\delta u_x+f_y\delta u_y+f_z\delta u_z\right)]\text{d}x=0\end{array}$ (1)

式(1)中L为立管长度， $m_r={\rho }_r{A}_r$ 是立管的单位长度上的质量， ${\rho }_r$ 是立管的材料密度， $A_r$ 是立管的截面积； $m_a={\rho }_w\pi D^2/4$ 是立管单位长度上的附加质量，D为立管外径， $R_r=D/2$ (此处暂不考虑浮力块，如果有浮力块，可由FEM处理)； ${\rho }_w$ 为立管周围水的密度(此处假设水密度是恒定的，如果不是恒定的，则由FEM处理)； $m_f={\rho }_f{A}_i$ 是立管内单位长度上的流体的质量， ${\rho }_f$ 为立管内的流体密度， $A_i$ 为立管内截面积， $u_i$ 是管内流速， $p_i$ 是管内压力； $u_x$ 、 $u_y$ 和 $u_z$ 和分别是立管在x轴，y轴和z轴方向上的平移位移； ${\stackrel{\dot{}}{u}}_y=\partial u_y/\partial t$，${\stackrel{\dot{}}{u}}_z=\partial u_z/\partial t$，${\stackrel{\dot{}}{u}}_x=\partial u_x/\partial t$ ； ${u^{″}}_y={\partial }^2{u}_y/\partial x^2$，${u^{″}}_z={\partial }^2{u}_z/\partial x^2$ ； ${\theta }_x$ 和 ${{\theta }^{\prime }}_x=\partial {\theta }_x/\partial x$ 是立管横截面相对于柱体轴线的角位移和角速度， ${{\theta }^{\prime }}_x=\partial {\theta }_x/\partial x$ 是扭转应变；E和G分别是立管的杨氏模量和剪切模量，I和I_x分别是立管横截面相对于直径和x轴的惯性矩，J_x是立管单位长度质量对x轴的惯性矩； ${\epsilon }_x^0$ 是由组装和加工引起的立管的初始轴向应变(包括管道的温度应变 ${\alpha }_t\Delta t_r$ 其中 ${\alpha }_t$ 为线性膨胀系数， $\Delta t_r$ 为管子长度方向的温差)， ${u^{\prime }}_x=\partial u_x/\partial x={\epsilon }_x$ 除 ${\epsilon }_x^0$ 以外的因素引起的轴向应变；f_x，f_y和f_z分别是重力，浮力和管道与水动力的相互作用力；T₀是立管两端的附加张力，不包括立管的重力和浮力或其他力(例如摩擦和阻尼)，因此其大小不会随长度而变化，例如，如果立管下端悬挂BOP，则T₀等于BOP的重量； $c=c_s+c_f$ 是阻尼系数，其中c_s是结构阻尼 (结构中的粘性耗散)，c_f是流体阻尼。

采用空间细长梁单元将立管离散为有限元模型。设其形函数为N_x(x)，对式(1)在一个单元上进行积分，可得

$\left[m\right]\left\{{\ddot{f}}_{i;j}\right\}+\left[c\right]\left\{{\stackrel{\dot{}}{f}}_{i;j}\right\}+\left[k\right]\left\{f_{i;j}\right\}=\left\{F_{i;j}\right\}$, (2)

其中

$\left\{f_{i;j}\right\}={\left[\begin{array}{cccccccccccc}{u}_{xi}& u_{yi}& u_{zi}& {\theta }_{xi}& {\theta }_{yi}& {\theta }_{zi}& u_{xj}& u_{yj}& u_{zj}& {\theta }_{xj}& {\theta }_{yj}& {\theta }_{zj}\end{array}\right]}^{\text{T}}$, (3)

$\left[m\right]={\displaystyle {\int }_0^l{\left[N_x\right]}^{\text{T}}\left[\begin{array}{cccc}{m}_r+m_f+m_a& 0& 0& 0\\ 0& m_r+m_f+m_a& 0& 0\\ 0& 0& m_r+m_f+m_a& 0\\ 0& 0& 0& J_x\end{array}\right]\left[N_x\right]\text{d}x}$, (4)

$\left[c\right]=\left[c_i\right]+\left[c_f\right]+\left[c_s\right]$, (5)

$\left[c_i\right]=2u_i{\rho }_f{A}_i{\displaystyle {\int }_0^l{\left[N_x\right]}^{\text{T}}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right]\left[{N^{\prime }}_x\right]\text{d}x}$, (6)

$\left[c_f\right]=c{\displaystyle {\int }_0^l{\left[N_x\right]}^{\text{T}}\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& R_r\end{array}\right]\left[N_x\right]\text{d}x}$, (7)

$\left[c_s\right]=2{\zeta }_1\left[m\right]{\Omega }_{r1}-2\frac{\left({\Omega }_{r1}^2\left[m\right]-\left[k_b\right]\right)\left({\zeta }_1{\Omega }_{r1}-{\zeta }_2{\Omega }_{r2}\right)}{{\Omega }_{r1}^2-{\Omega }_{r2}^2}$ (8)

$\left[k\right]=\left[k_a\right]+\left[k_b\right]+\left[k_t\right]$, (9)

$\left[k_a\right]={\displaystyle {\int }_0^l{\left[N_x\right]}^{\text{T}}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& \left(u_i^2{\rho }_f+p_r\right)A_i-T_0& 0& 0\\ 0& 0& \left(u_i^2{\rho }_f+p\right)A_i-T_0& 0\\ 0& 0& 0& 0\end{array}\right]\left[{N^{″}}_x\right]\text{d}x}$, (10)

$\left[k_b\right]={\displaystyle {\int }_0^l{\left[{N^{″}}_x\right]}^{\text{T}}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& EI& 0& 0\\ 0& 0& EI& 0\\ 0& 0& 0& 0\end{array}\right]\left[{N^{″}}_x\right]\text{d}x}$, (11)

$\left[k_t\right]={\displaystyle {\int }_0^l{\left[{N^{\prime }}_x\right]}^{\text{T}}\left[\begin{array}{cccc}E{A}_r{\left({\epsilon }_x^0+1\right)}^2& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& G{I}_x\end{array}\right]\left[{N^{\prime }}_x\right]\text{d}x}$, and (12)

$\left\{F_{i;j}\right\}={\displaystyle {\int }_0^l{\left[N_x\right]}^{\text{T}}{\left[\begin{array}{cccc}{f}_x& f_y& f_z& 0\end{array}\right]}^{\text{T}}}\text{d}x$. (13)

注意式(5)和式(6)中的 $\left[c_i\right]$ 表示由于管道内流体的流动和管道横截面转动而引起的科氏惯性效应，是一种表观阻尼；式(7)中的c是流体阻尼；式(8)中的 ${\Omega }_{r1}$ 和 ${\Omega }_{r2}$ 分别是立管的最低第一和第二角频率， ${\zeta }_1$ 和 ${\zeta }_2$ 分别是该结构的一阶和二阶阻尼比，一般 ${\zeta }_1=0.001\text{-}0.03$，${\zeta }_2=\min \left({\Omega }_{r2}/{\Omega }_{r1},5\right){\zeta }_1$ ；式(13)是单元的结点载荷列向量，对于垂直放置的立管和垂直向上的x轴，我们有

$f_x=-f_b\left(m_r+m_a\right)g+f_{bb}{m}_{b}g$, (14)

其中，将浮力系数 $f_b$ 定义为 $f_b=1-\pi D^2{\rho }_w/4\left(m_r+m_a\right)$，将浮力块的浮力系数 $f_{bb}$ 定义为 $f_{bb}=1-{\rho }_b/{\rho }_w$，

${\rho }_b$ 和 $m_b$ 分别为浮力块密度和浮力块质量。而 $f_y$ 和 $f_z$ 是由于涡旋脱落和涡激效应而由外部流体动力学产生的，可以表示为 [14]

$f_y=f_D+f_{D0}-\frac{f_L}{u}{\stackrel{\dot{}}{u}}_y$,$f_z=f_L+\frac{f_D}{u}{\stackrel{\dot{}}{u}}_y$.(15)

在这里， $f_{D0}$ 是平均阻力， $f_D$ 是涡流引起的波动阻力， $f_L$ 是管道上的升力。 这些可以表示为

$f_D=\frac{1}{2}{C}_{Di}{\rho }_{w}D{u}^2$,$f_{D0}=\frac{1}{2}{C}_{D0}{\rho }_{w}D{u}^2$,$f_L=\frac{1}{2}{C}_L{\rho }_{w}D{u}^2$, (16)

其中， $C_{Di}$ 是随时间变化的涡流诱导阻力项， $C_{D0}$ 是两端固支管道的平均阻力系数， $C_L$ 是升力系数。 $C_{Di}$ 和 $C_L$ 可以表示为

$C_{Di}=C_{Di0}\frac{p\left(x,t\right)}{2}$,$C_L=C_{L0}\frac{q\left(x,t\right)}{2}$, (17)

其中 $C_{Di0}$ 和 $C_{L0}$ 是两端固支立管相关的发生涡旋脱落的非稳态阻力和升力系数。 $p\left(x,t\right)$ 和 $q\left(x,t\right)$ 分别代表在y和z方向上的涡激变量，其模型将在下一节详述。

顺便指出，在式(16)中，当立管被外径为D_b的浮力块所包裹时， $D=D_b$ ；如果立管或立管的一部分中没有浮力块，则D即为立管的外径。

此外，当立管的挠度变大时，挠度斜率将改变载荷的分布，因此，应将式(14)和(15)所表示的结点力分量可修正为

$f_{xm}=f_x\cos {\alpha }_i$,$f_{ym}=f_y-f_x\sin {\alpha }_i\cos {\beta }_i$,$f_{zm}=f_z-f_x\sin {\alpha }_i\sin {\beta }_i$ (18)

其中， $\sin {\alpha }_i=\Delta r_{yz}/\Delta r_{xyz}$，$\cos {\beta }_i=\Delta r_y/\Delta r_{yz}$，$\sin {\beta }_i=\Delta r_z/\Delta r_{yz}$，$\Delta u_x=u_{xj}-u_{xi}$，$\Delta u_y=u_{yj}-u_{yi}$，$\Delta u_z=u_{zj}-u_{zi}$，$\Delta r_{yz}=\left|u_{yj}-u_{yi},u_{zj}-u_{zi}\right|$，$\Delta r_{xyz}=\left|u_{xj}-u_{xi},u_{yj}-u_{yi},u_{zj}-u_{zi}\right|$。在模拟过程中，每个时间步对式(18)的值进行

更新。修正后更能反映实际情况。

2.2. 范德波尔尾迹振子模型

式(17)中两个尾迹变量 $p\left(x,t\right)$ 和 $q\left(x,t\right)$ 为无量纲量，可由以下两个van der Pol方程控制 [14]：

$\ddot{q}+{\epsilon }_z{\Omega }_s\left(q^2-1\right)\stackrel{\dot{}}{q}+{\Omega }_s^{2}q=\frac{A_z}{D}{\ddot{u}}_z$, (19)

$\ddot{p}+2{\epsilon }_y{\Omega }_s\left(p^2-1\right)\stackrel{\dot{}}{p}+4{\Omega }_s^{2}p=\frac{A_y}{D}{\ddot{u}}_y$, (20)

其中， ${\ddot{u}}_y$ 和 ${\ddot{u}}_z$ 为立管横向振动加速度，由方程(2)确定； ${\epsilon }_y$ 、 ${\epsilon }_z$ 、 $A_y$ 和 $A_z$ 是通过实验估算的无量纲参数，是涡流脱落角频率，定义如下 [14]：

(21)

称为斯托哈尔(Strouhal)常数，与雷诺数Re有关，在范围内(本文所属范围)，可近似取0.2 [15] [16]。

2.3. 方程求解

系统包括两组基本方程。一组为管道在内部流体流动和外部流体尾迹激励共同作用下的振动方程，如式(1)~(2)所示；另一组为受管道横向振动加速度激励的尾迹振子方程，如式(19)和(20)所示。式(2)中的非齐次项经过式(15)~(18)的非线性转换与式(19)和(20)相耦合。因此，要分析立管或尾迹的动力学特征，必须联合求解方程(2)~(20)。具体求解步骤如下：

i) 求解方程(2)齐次式的特征值，获取立管的最低第一和第二角频率和，然后通过式(8)计算；

ii) 将初始条件()设置为，，，以及；

iii) 计算式(4)~(7)和(9)~(18)，并将其代入公式(2)中；

iv) 用Newmark方法求解方程式(2)，提取和，将它们和式(21)代入式(19)和(20)中；

v) 用Runge–Kutta方法求解方程(19)和(20)，提取和，然后将它们代入方程(13)~(18)和(2)中，完成一个循环。让(可设)。转回到步骤iv)继续循环计算，直到。

下标，N表示有限元结构中的节点数；，是从0到60s的时间步长总数。通常，在求解式(2)之前，必须将局部单元坐标表示的所有相关矩阵和向量转换为全局结构坐标，并集成到结构矩阵和向量中。

2.4. 模型验证和参数设定

在模拟隔水管下入过程的动力学响应之前，我们通过将本模型的计算结果与已有的实验结果和数值模拟结果进行比对，来验证本模型的正确性和调试参数值。为此，我们选择的模拟对象其结构参数和材料参数与一个经典实验里的完全相同 [17]：无浮力块，立管长度，管道外径，壁厚，立管长度与直径之比，质量比，外部流体的质量密度，内部流体的质量密度，管道的杨氏模量。此外，来流平行于全局坐标下的y轴，并且附加张力施加到立管的每个端部。立管的两端均由铰链支撑，可以自由旋转。垂直立管的固有频率由埃克森美孚公司在MARINTEK测试得到 [18]，并由L. Wang等人的解析解 [13] 和E. Wang等人 [14] 在ANSYS软件中使用分析方法和FEM进行了计算.这些结果和本模型获得的结果汇总在表1。本模型的计算结果和有限元计算的结果 [14]，二者的相对误差Error₅₄ < 1%，与其他结果的误差也在工程许可的范围内。

表1. 固有频率对比

^*Error_ij (i= 3, 5; j = 2, 4)是第i列到第j列时数据的相对误差。

表2. 最大振幅均方根对比

^*Error_ij(i= 3, 5; j = 2)是第i列到第j列时数据的相对误差。

接下来，我们检查立管的振动幅度。与Lehn [18] 实验中的实验组#1105进行比较。除了来流速度设定为均匀且大小为0.42 m/s，实验组＃1105中的其他参数与参考文献 [18] 中的参数相同。表2总结了以无量纲形式表示的IL(y)方向和CF(z)方向振幅的最大均方根值(rms)，和，这些结果分别取自埃克森美孚公司在MARINTEK的实验 [18]，Huang等人 [18] 和Wang等人 [14] 的有限元计算，以及本文方法的计算结果。由表2可见，本方法的结果与实验吻合度最好。我们在对比计算中采用的参数值与有关文献 [18] 中大致相同，但是在对比计算中进行了微调，现将最佳值列于表3，作为今后计算时本模型的参数值。

表3. 模型主要参数值

以上与实验，计算和解析解方法得到的结果比较证明，本模型包括算法和编程是可信且有效的。在下一节中，我们将使用本方法来计算和讨论钻探立管在深水区下落时的动力学响应，以及带浮力块的钻井立管的力学行为分析。

3. 钻探立管在深水区下入时的动力学响应

3.1. 计算方案

假设水深分别为500 m、1000 m、1500 m、2000 m、2500 m和3000 m，立管在水下的全长L与对应水深相等，立管外径 [19] ，内径，固体立管的杨氏模量，泊松比，管道的质量密度；管内注满海水，其质量密度，管内流体输运速度；立管进入深海速度、0.15 m/s和0.3 m/s；防喷器在空气中的重量分别为Pex = 1000、2000、3000 kN；管外流体(海水)质量密度为常数。浮力块的质量密度，浮力块的外径；浮力块从立管底端向上连续布置，其总长度为L_b (由于外部流体的速度通常会随着深度的增加而降低，降低模块的安装位置能够减小流体动力对管柱的力学影响)。我们以浮重比R_BW度量浮力块的作用，本文仅考虑，即安装浮力块后能使立管上端张紧器的静态张力减小一半，并由此R_BW的值推算L_b。

另外，假设在深水条件下，由波浪引起的水质点运动的速度遵循平面波的线性理论，例如，在xy平面内，在水深x和时刻t，速度u_y为 [20]

, (,), (22)

其中H_w是波高，T_w是波周期，波数为，g为重力加速度，角频率为。波长定义为。此外，由风和潮流引起的水质点速度分别为 [20]。

, , (), (23)

其中v_m和v_t是水面上风流速和潮流速。在本文中，设，，，，并且假定潮流速的方向与y轴平行，风引起的流速和波浪引起的流速方向一致且与z轴平行。其他主要参数值见表3。

3.2. 代表性结果

为了对本模型模拟情况有一个大致了解，图2(a)~(h)给出了带有浮力块的立管进入到1000米水深时几个典型的动力响应图形和曲线，此时上端送入速度。浮重比，对应的浮力块总长(图2(a)和图2(b)绿色原点组成的长度)。图2的图题中简要说明各子图(a)~(h)的含义。值得注意的是，增加浮力块虽然可以减轻张紧器和平台的负担(如图2(i))，但会增加海流对立管的尾迹效应(如图2(f)~(h))，增大下端横向位移(如图2(a)~(d))和局部弯曲应力(如图2(e)，图2(j))，从而增加下入立管对水下井口的对准难度，还可能加速立管的疲劳损伤。

图2. 带有浮力块的立管进入到1000米水深时典型的动力响应结果。(a) 一个时间点上立管在xy平面内的变形；(b) 一个时间点上立管在xz平面内的变形；(c) 各时间点立管在xy平面内的变形；(d) 各时间点立管在xz平面内的变形；(e) 各时间点立管最大弯曲应力值和所在位置；(f) van der Pol变量p和q在不同的时间点围绕立管的分布情况；(g) p沿立管轴线的分布随时间的演化；(h) q沿立管轴线的分布随时间的演化；(i) 立管轴向拉、压应力时程曲线(顶部应力s_Nxt、中部应力s_Nxm和底部应力s_Nx1)；(j) 立管中部弯曲应力时程曲线(弯矩M_y引起的应力s_Mym，弯矩M_z引起的应力s_Mzm，由合弯矩M_b引起的应力s_Mbm)

3.3. 立管上端送入速度和下端悬挂重量的影响

本节假设水深(管长) L = 3000 m，管下端悬挂重量分别为Pex = 1000 kN、2000 kN和3000 kN，不加浮力块，流速等其他条件如在3.1所述。当立管下入到接近3000 m时，设立管上端的送入速度分别为u₀ = 0，u₀ = 0.15和u₀ = 0.3 m/s，以此为初始条件模拟60 s。考察这后20 s~60 s立管的下端最大横向位移、上端的最大拉伸应力和全管最大弯曲应力，比较结果如图3和图4所示。可见，在悬吊重量相等的情况下，u₀的增加会使最大y方向的位移u_y略微降低(图3(a))，而最大u_z几乎没有受到影响(图3(b))。随着u₀的增加，最大弯曲应力(图4(a))和最大拉伸应力(图4(b))略有下降。此外，Pex的增加使得最大u_y (图3(a))和最大拉伸应力(图4(a))都稍微下降，并使得最大弯曲应力也显著降低，比如，与Pex = 1000 kN相比，Pex = 2000 kN时下降27.6%，Pex = 3000 kN时下降53.8%。这是因为IL流向是在y方向上，导致垂直立管向该方向弯曲，而较大的Pex大大降低了IL流弯曲立管的弯曲曲率并有使立管向铅垂线靠拢的趋势，因此，u_y和弯曲应力减小。对于较大的u₀和Pex下的拉应力的降低，前者是由于水的轴向阻尼作用，而后者是由于悬重可以抑制轴向振动的幅值而发生的。从图4中可明显得出Pex越小，振幅越大，最大拉应力越大；Pex越大，振幅越小，最大拉应力越小。此外，较大的Pex增加了负z向的最大位移，如图3(b)所示。立管z向位移可以归结为流体动力f_z，即等式(15)~(17)中定义的尾流动力。由于较大的Pex减小了立管的弯曲曲率，且垂直直管的尾迹效应大于倾斜直管 [21]，因此会产生u_z。然而，最大的u_z与最大的u_y数值量级相比很小，可以忽略不计。

因此，适当增加下入速度和悬重可以降低最大弯曲应力和最大拉力应力，只要上部平台和张紧器足够强，几乎不会对立管进入产生负面影响。

图3. (a) 立管下端在来流(y)方向上的最大位移u_y和(b) 在横流(z)方向上的最大位移u_z随底端悬重(Pex)和上端送入速度(u₀)的变化

图4. (a) 立管最大拉应力和(b) 最大弯曲应力随底端悬重(Pex)和上端送入速度(u₀)的变化

3.4. 不同水深和悬重的影响

设水深分别为500 m，1000 m，1500 m，2000 m，2500 m，3000 m。当立管下入到对应的水深时，停止送入(u₀ = 0)，以此为初始条件再模拟60 s，考察这后20 s~60 s立管的下端最大横向位移、上端的最大拉伸应力和全管最大弯曲应力。立管上不加浮力块，下端悬重不同，分别为Pex = 1000 kN，2000 kN，3000 kN，海况等其他条件如第3.1小节所述。比较结果如图5和图6所示。

从图5(a)可以看出，当Pex为常数时，最大u_y随着L的增大而增大。在L = 3000 m时的最大u_y是在L= 500 m时的2.5倍。随着L的增大，最大u_y的增大率明显减小。此外，当L保持恒定时，最大u_y随着Pex的增大而减小，因为较大的Pex减小了底端位移(例如，当L = 3000 m时，Pex = 3000 kN时的最大u_y是Pex = 1000 kN时的0.9倍)。

图5. (a) 立管底端的y方向(u_y)位移和(b) z方向(u_z)位移随管长度(L)和悬重(Pex)的变化

图6. (a) 立管最大拉应力和(b) 最大弯曲应力随管长(L)和悬重(Pex)的变化

与图5(a)中的u_y相比，图5(b)中最大u_z的变化较为复杂。u_z在CF方向上，主要受复杂的尾迹控制，而u_y在IL方向上，主要受来流影响。比较(a) (b)中u_z和u_y的数量级，此时尾流效应可忽略不计。

图6(a)为最大拉伸应力随着L的增大而增大，符合预期结果；但随着吊重的增加，最大拉应力略有下降，此现象可归因于吊重抑制了立管的轴向振动。

从图6(b)可以看出，随着L的增大，最大弯曲应力基本保持不变，但随着悬挂重量的增加而显著减少。其原因可以归结为最大弯曲应力取决于立管的曲率，前者，增加L可能会增加横向位移，但曲率本身可以保持不变；后者，是因为吊重可以将弯曲的立管拉直，从而降低立管的曲率和弯曲应力。

4. 结论

本文针对内部输运流体的柔性立管发展了一种新的VIV (涡激振动)三维动力学模型，以便能模拟立管在实际环境条件下进入深度为3000 m甚至更深的水中时的动力学响应。重点讨论了悬重、下入速度以及水深对立管底端横向位移(事关防喷器的落点与下井口的对准和控制)以及立管上端最大拉应力(涉及张紧器和平台的提升能力)和最大弯曲应力(涉及立管的疲劳损伤)的影响，并得出如下结论：

1) 在一般海况下，将钻探立管下入到深水中时，降低下入速度对落点偏移量没有明显的影响，但可以少许减小局部最大弯曲应力；

2) 尽管深水需要较长的钻探立管，但长度大的立管自重也大，增加了下垂能力，因此不会增大落点的偏移量；

3) 当水深一定时，增加悬重可以减小落点偏移量，但立管的上端拉力会随悬挂重量的增加而增加，因此需要更强劲的张紧器和排水量更大的平台；

4) 增加浮力块虽然可以减轻张紧器和平台的负担，但会增加海流对立管的尾迹效应，从而加大落点偏离和最大弯曲应力，因此，浮力块加与不加以及加多少必须综合考虑和优化；

5) 本文的仿真结果与现有实验数据具有很好的一致性，但对更复杂情况，本文设置的载荷偏理想化，有待在实践进一步细化和验证。

基金项目

国家自然科学基金(U1663205)。

文章引用

张 颖,朱卫平. 深水立管入水过程中涡激耦合振动特性分析
Dynamic Characteristics of Drilling Risers as Being Entered into Deep Water[J]. 石油天然气学报, 2020, 42(02): 124-137. https://doi.org/10.12677/JOGT.2020.422024

参考文献