随机荷载下陡波形立管的非线性动力分析

顾洪禄， 李效民， 郭海燕， 崔 鹏， 刘 震， 李福恒

(中国海洋大学 工程学院,山东 青岛 266100)

近年来，随着海洋油气开采逐渐向深水、超深水过渡，简单钢悬链线立管(steel catenary riser，SCR)顶部张力水平显著升高，此外，剧烈的顶部浮体运动导致悬挂点、触地点处极易产生疲劳破坏[1]。因此，SCR可能已经不再适用于深海油气开采[2]。为改善SCR的疲劳和强度问题，波形立管(缓波形立管、陡波形立管)应运而生，通过在SCR中间管段安装一定数量的浮力块，形成拱弯段，波形立管大幅度降低了顶张力的水平，解耦了顶部浮体与触地点之间的运动。其中，陡波形立管构型陡、跨距小，近地段与海床接近垂直，能有效避免海床稳定性问题。并且，陡波形立管对管内流体变密度的情况适应性较高，在深海油气开采中具有显著的优越性[3]。

然而，陡波形立管分析设计需要面临整体构型复杂、边界条件处理难度高、计算模型难收敛、内部流体非线性耦合等众多突出难点。此外，陡波形立管在环境荷载下极易产生大位移、大转角，具有显著的几何非线性，进一步增加了模型建立以及结构动力分析的难度。因此，对陡波形立管进行非线性动力分析具有重要的学术价值和工程实际意义。

国内外学者对顶张力立管、SCR的静动力响应、管土相互作用等进行了广泛的研究[4-12]。但是对于波形立管的研究相对较少，且主要集中于缓波形立管。

Li等[13]忽略抗弯刚度的影响，对缓波形立管动力响应进行了初步分析，该研究将各管段简化为悬链线，不能够考虑海流以及内流的影响。Santillan等[14-15]将陡波形立管简化为细长弹性杆，利用有限差分法对立管进行了静动力分析。基于相同的模型，通过在缓S形和陡S形立管拱弯点处施加集中力模拟浮子段，对两种立管进行了参数敏感性分析，这与陡波形立管强几何非线性以及浮子段均匀受力的实际情况存在差异。Sun等[16]利用OrcaFlex对陡波形立管的动力响应进行了参数敏感性分析。Yang等[17]对不同参数下的缓波形立管疲劳寿命进行了一系列参数敏感性分析。Wang等[18-19]考虑内部流体流动和管土相互作用，基于小变形梁、非线性大变形梁理论建立了缓波形立管二维非线性动力模型，利用有限差分法求解立管动力响应，然后对缓波形立管静力性能进行了一系列参数敏感性分析。此外，对缓波形立管在简谐顶部运动下的动力响应进行了研究[20]，但是以上研究仅仅局限在二维平面内。李艳等[21]应用三维集中质量模型，对缓波形立管在顶部浮体作简谐运动时的动力响应进行了研究。宋磊建等[22]计算了规则波作用下深海脐带缆在悬链线布局和缓波形布局下张力、曲率的动态响应。

上述研究均基于传统杆、梁理论建立波形立管力学模型。众所周知，在传统杆、梁理论中，需要将物理参数由局部坐标转换到整体坐标，这会显著增加计算时间和对计算机CPU的要求。近年来，一种基于绝对坐标的特殊弹性杆理论被广泛应用于波形立管的数值模拟中。

基于大变形杆理论和有限单元法，Ruan等[23]对弹性海床上的缓波形脐带缆的静力力学性能进行了一系列参数分析。基于相同的理论，Ruan等[24]和Cheng等[25]对缓波形立管在顶部浮体运动、波流荷载、内部流体和海床摩擦共同作用下的静动力响应进行了一系列参数敏感性分析。Kim等[26]对简单悬链线立管与缓波形立管在相同环境荷载下的整体静动力响应进行了比较分析，但并未对缓波形立管进行深入研究。Qiao等[27]基于大变形梁理论和力平衡原则建立了陡波形立管运动方程，结合有限差分法和打靶法对陡波形立管在海流作用下的静力性能进行了分析。Liu等[28]基于柔性杆理论建立了陡波形输流立管的运动方程，并对陡波形立管的静动力响应进行了一系列参数敏感性分析。

上述研究中，直接在全局坐标系中得到包含所有几何非线性的立管控制方程，避免了繁琐的坐标转换，极大节省了计算资源。本文正是基于此类模型——柔性杆模型[29]建立陡波形立管运动方程。

此外，从上述关于波形立管研究的讨论可以看出，目前大多数研究忽略了内部流体的作用，将顶部浮体运动、波浪简化为简谐运动，这与实际海况不符。并且目前对于陡波形立管的动力分析十分缺乏，关于陡波形立管在随机波浪、顶部浮体随机运动和内部流体共同激励下的非线性动力分析尚未见报道。

本文基于柔性杆模型，充分考虑陡波形立管几何非线性，建立了深水陡波形立管运动方程，采用有限元法离散立管运动方程，并利用Newton-Raphson法和Newmark-β法分别计算立管的静、动力响应，进一步计算得到了陡波形立管在内部流体、随机波浪、顶部浮体静态偏移及动态响应等多种因素共同作用下非线性动力响应，并对立管应力水平进行一系列参数敏感性分析。

1 基本理论

1.1 柔性杆模型

如图1所示，杆的位形由杆轴线位置表示，在三维笛卡尔坐标系中，空间曲线r(s,t)是弧长与时间的函数，表示柔性杆变形后轴线位置状态。

图1 柔性杆坐标系示意图

假设变形前后杆的弧长不发生改变；杆上任意一点的内力状态可以完全由作用在轴线上的合力和力偶表示；忽略转动惯量、剪切变形、均布扭矩和分布外力矩的作用；假定杆可伸长且伸长量为小量。根据动量守恒与动量矩守恒可得[29-30]

(1)

M′+r′×F+m=0

(2)

式中:ρ、mf为杆、内部流体单位质量；q为杆单位长度所受外力；m为杆单位长度所受外力矩。

由Euler-Bernoulli梁理论，内力矩M可以表示为

M=r′×EIr″+Hr′

(3)

式中:EI为截面弯曲刚度；H为扭矩，在实际海况下，柔性立管不受均布扭矩作用，其自身扭矩也可忽略。因此其运动方程可以写为

(4)

当考虑杆的伸长量为小量，可以认为杆件变形条件如下

(5)

式中:EA为杆件轴向刚度；λ为拉格朗日乘子。

综上，式(4)与式(5)构成柔性杆基本控制方程。

1.2 随机波浪与顶部浮体运动模拟

由于JONSWAP谱适用范围更广，适用于不同成长阶段的波浪，本文选取DNV RP-C205 (2007)[31]中推荐的JONSWAP谱

(6)

(7)

式中:Aγ=1-0.287ln(γ)为正则化因子；HS为特征波高；ωp=2π/TP为峰值角频率,TP为峰值周期；γ为无量纲峰形参数，当γ=1时，JONSWAP谱简化为P-M谱，本文按照如下公式进行取值

然后基于Longuet-Higgins模型[32]对随机波浪进行模拟。

顶部浮体运动与海面运动状况直接相关，顶部浮体对于特征波高的频率响应一般用运动幅值响应算子进行估算。顶部浮体响应采用Sexton等[33]推荐的运动模型表示

(8)

式中:等式右端第一项表示顶部浮体平均静偏移;第二项表示一阶低频响应;第三项表示二阶高频响应；SL、TL分别为顶部浮体慢漂单边幅值、周期；慢漂运动与波浪之间相位差αL一般取0；Sn为顶部浮体对特定波浪(周期Tn=2π/ωn，幅值An)的幅值响应；kn、φn分别为波数、波浪相位角；αn为顶部浮体与波浪周期之间的相位差。

随机波浪荷载通过莫里森方程求解；由于顶部浮体的尺寸与水深相比可以忽略不计，因此将顶部浮体运动转化为立管动边界问题。

1.3 有限元模型

将柔性杆运动方程和变形条件分别写成张量形式，利用三次Hermit插值函数进行离散，运用Galerkin方法得到矩阵形式的微分方程如下

(9)

(10)

通过Newton-Raphson法和Newmark-β法求解式(9)、式(10)便可以得到立管的静动力响应。

立管截面应力由下式进行计算

(11)

式中:A为立管横截面积；E为立管弹性模量。

2 模型验证

基于第1章所述理论，利用MATLAB编写陡波形立管非线性动力响应分析程序DRSWR。为验证本文程序的有效性，选取文献[34]陡波形立管模型，具体参数如图2、表1所示。在本文中下降段、浮子段、悬挂段的划分单元个数分别为10、30、40，进一步细化单元数并不能显著提高计算精度，并且会极大地消耗计算时间。数值模拟持续时间为500 s，时间步为0.25 s，以上参数在后续计算中不发生变化。利用本文程序对陡波形立管进行分析计算，将本文结果与相同参数下OrcaFlex计算结果进行对比验证。

表1 陡波形立管参数

图2 陡波形立管构型(m)

2.1 随机波浪动力分析验证

本文选取文献[35]中波浪数据如下：特征波高Hs=6.5 m，峰值周期TP=12.82 s。利用DRSWR对陡波形立管在随机波浪作用下的应力包络图进行计算，将本文计算结果与同条件下OrcaFlex计算结果进行对比。如图3所示，可以看出，两者计算结果吻合程度较好。

(a) 最大应力包络图

2.2 顶部浮体动力分析验证

利用DRSWR对立管在顶部浮体运动下的顶张力时程曲线进行计算。顶部浮体振幅为10 m，周期为27 s，沿X方向运动。本文计算结果与OrcaFlex计算结果进行对比，如图4所示，两者计算结果基本一致。

图4 顶部张力时程变化图

经过上述对比分析，一定程度上验证了DRSWR的有效性。

3 陡波形立管非线性动力分析

本文对陡波形立管在随机波浪、顶部浮体运动共同激励下的动力响应进行了分析计算。

选用Li等研究中的环境荷载数据，具体参数如下：随机波浪特征波高HS=6.5 m，峰值周期TP=12.82 s；顶部浮体平均静偏移S0=5 m，慢漂单边幅值SL=5 m，慢漂周期TL=200 s，内部流体密度ρf=0，内部流体速度v=0。以上数据为本文计算的基本环境荷载数据，后续计算分析中，无特殊说明，均采用以上环境参数。图5绘制了DRSWR模拟的随机波浪高度、顶部浮体运动时程曲线，可以看出两者运动均具有明显的随机性。

(a)

图6给出了陡波形立管在随机波浪、顶部浮体作用下的应力包络图，可以看出悬垂段、浮子段应力变化幅度较大，而下降段变化幅度较小。这主要是由于波浪入射速度随水深呈指数衰减，海面处波浪最为剧烈，此外，顶部浮体运动直接作用于悬挂点，因此悬垂段应力变化较为剧烈；由于添加了浮力块，导致浮子段受到更高水平的水动力荷载，因此该段应力变化幅值较大；浮子段解耦了顶部浮体运动与立管底端之间的动态响应，导致下降段应力幅值变化并不明显。

由以上分析可以看出陡波构型降低悬挂点处张力水平，减轻顶部浮体负荷的同时，还可以减小井口附近立管应力变化幅值，提高井口处立管疲劳寿命。除此之外，沿立管长度存在两个应力极大值点，分别位于拱弯点(760 m)、垂弯点(1 192.5 m)处，这主要是两点初始曲率较大所致。

应力水平是立管强度校核的重要依据，并对立管疲劳损伤具有重要的参考价值。为进一步增进对陡波形立管在随机波浪、顶部浮体共同作用下动力响应的理解，本文对随机波浪(特征波高、峰值周期)、顶部浮体运动(平均静偏移、慢漂运动)、内部流体(内流密度、内流速度)等因素对立管应力的影响进行了参数敏感性分析。

3.1 随机波浪参数敏感性分析

实际工程中，随机波浪参数由长期观测结果整理成的波浪散布图确定。本文为便于研究波浪参数对动力响应的影响，对参数的选取进行了简化。

3.1.1 特征波高参数敏感性分析

保持其他参数不变，选取特征波高HS=6.5 m、16.5 m、26.5 m，分别计算陡波形立管最大应力包络图。

如图7所示，分别给出了沿立管长度的最大应力包络图以及浮子段、悬垂段最大应力包络图的放大图，由于下降段应力水平较低并且变化幅值较小，不再针对该管段进行相关分析。由图7可以看出立管整体应力水平随特征波高的增加而升高，波高越大该趋势越为明显；各管段中，悬垂段对该参数变化敏感度最高，并且越靠近悬挂点敏感度越高；此外，随特征波高的增加，悬垂段应力变化的非线性特性逐渐增强。

(a)

3.1.2 峰值周期参数敏感性分析

保持其他参数不变，分别计算峰值周期TP=2.82 s、7.82 s、12.82 s时立管的最大应力包络图，得到峰值周期对立管应力的影响，如图8所示。

(a)

可以看出立管应力随峰值周期增加逐渐降低，在低周期下，该趋势更为显著；在低周期波浪作用下，立管应力水平急剧升高，且最大应力包络图曲线不再圆滑并出现震荡，应力的非线性明显增强；悬垂段对峰值周期变化的敏感度仍高于其他管段，且悬挂点附近敏感度最高。

由以上分析可以看出，陡波形立管对随机波浪的变化敏感度较高，大波高、低周期的随机波会导致立管整体应力水平急剧升高且表现出显著的非线性特性；悬垂段对随机波浪的变化敏感度最高，并且越靠近悬挂点，敏感度越高，主要是由于越接近海面，波浪作用越显著。因此在对立管进行强度校核时，应充分考虑立管所处海域波浪状况，并采取必要的消波措施。

3.2 顶部浮体运动参数敏感性分析

3.2.1 平均静偏移参数敏感性分析

保持其他参数不变，当顶部浮体静偏移S0=5 m、55 m、105 m时，分别绘制立管最大应力包络图，如图9所示。

(a)

由图9可以看出，关键点(拱弯点、垂弯点、悬挂点)附近应力对该参数变化较为敏感；随静偏移的增加，拱弯点、垂弯点附近应力逐渐降低，主要是由于这两处具有较高的初始曲率，应力水平主要由曲率控制，静偏移增大会导致这两点处构型更为平缓，曲率降低；随静偏移的增加，悬挂点附近应力升高，是由于悬挂点应力水平主要由张力控制，静偏移增大，立管整体张紧，导致悬挂点产生更高水平的张力。

3.2.2 慢漂运动参数敏感性分析

图10、图11分别绘制了不同慢漂幅值、不同慢漂周期下立管最大应力包络图。保持其他参数不变，分别选取慢漂幅值SL=5 m、20 m、35 m，慢漂周期TL=100 s、150 s、200 s。

(a)

(a)

由图10、图11可以看出顶部浮体慢漂运动只对悬挂点附近应力水平影响较大，应力水平随慢漂幅值SL的增加而升高，随周期TL的增加逐渐降低，主要原因是：顶部浮体直接作用于悬挂点，剧烈的顶部浮体慢漂运动会导致悬挂点附近应力水平升高。

由以上分析可以看出剧烈的随机波浪、顶部浮体运动都会造成悬挂点处高水平的应力，因此实际工程中应在悬挂点处设置弯曲加强构件或喇叭口等辅助装置，增强该处的抗弯刚度；此外，顶部浮体慢漂运动对立管整体应力的影响程度远低于平均静偏移，这是因为陡波构型提高了对顶部浮体瞬时运动的适应性。

3.3 内部流体参数敏感性分析

3.3.1 内流密度参数敏感性分析

立管在工作期间内部通常会充满多种流体介质，如石油、天然气、海水、泥沙等，不同密度的内部流体会对立管产生不同程度影响。本文选取3组不同内流密度：1.29 kg/m3、800.00 kg/m3、998.00 kg/m3，得到了不同内流密度下立管应力响应，如图12所示。

(a)

由图12可以看出，内流密度的变化对立管整体应力水平影响显著；随内流密度的增加，立管各段均出现明显的极值点；拱弯点附近区域应力水平随内流密度的增加逐渐减小；触地点、垂弯点附近区域变化趋势与拱弯点变化趋势正好相反。主要是因为内流密度增加使得立管湿重增加，导致浮子段拱弯部分构型更为平缓，拱弯点曲率降低；悬垂段整体构型弯曲程度变高，垂弯点曲率增大；由于下降段整体湿重的增加，出现一定幅度屈曲，曲率升高。

3.3.2 内流速度参数敏感性分析

立管在运行期间，内部流体以一定速度流动时与立管会产生耦合振动现象，对立管动力响应产生一定程度的影响。本文选取内流密度为800 kg/m3，内流速度为0、5 m/s、10 m/s时，分别计算立管最大应力包络图，如图13所示。

(a)

可以看出，立管整体应力水平随内流速度的增加而升高，尤其是当内流在高速流动时，趋势更为明显。悬垂段对于内流速度变化的敏感度明显高于其他管段，并且内流速度对立管弯曲段应力存在较为明显的影响，主要是由于内部流体流动会对弯曲管段内壁产生一定的流动压力，导致管段产生小幅度的弯曲振动耦合。

由以上分析可以看出，内部流体对各关键点处的应力影响显著，因此在对陡波形立管进行强度校核时，应对内流作用给予足够的重视。

3.4 多因素变化下陡波形立管动力响应分析

为进一步探究随机波浪和顶部浮体运动对陡波形立管动力响应的影响，本文选取3组不同的随机波浪工况(如表2所示)和顶部浮体慢漂幅值(SL=5 m、20 m、35 m)，对多因素共同变化下的关键点的最大应力进行了计算分析，计算结果如表3～表5所示。

表2 随机波浪工况

表3 SL=5 m时陡波形管关键点最大应力

表4 SL=20 m时陡波形管关键点最大应力

表5 SL=35 m时陡波形管关键点最大应力

在不同随机波浪条件下，关键点最大应力都随慢漂幅值增加逐渐增加；同样在不同慢漂幅值下，关键点最大应力随特征波高增加逐渐增加。这与3.1.1节单因素变化的结论是一致的。因此，可以看出多因素变化下，特征波高、慢漂幅值对陡波形立管动力响应的影响并没有发生改变。但是，关键点最大应力随峰值周期变化趋势与3.1.2节不同，是因为各因素对动力响应的影响作用相反时，各因素之间存在竞争机制，该算例中峰值周期对动力响应的影响明显小于特征波高和慢漂幅值。

此外，通过对比可以看出，特征波高、慢漂幅值对各个关键点影响的程度均为：悬挂点>拱弯点>垂弯点，这与3.1节～3.3节单因素变化下的结论也是一致的。

通过以上分析，我们可以看出，在多因素变化下，如果各因素变化同时增强动力响应，因素之间是相互独立的，变化趋势与单因素分析结果是一致的。反之，如果各因素变化趋势对动力响应的作用相反，单因素变化分析结果在多因素变化下则需要进一步分析确定。

4 结 论

本文基于柔性杆模型，充分考虑陡波形立管的几何非线性、内部流体非线性耦合，利用MATLAB编写计算程序DRSWR，对陡波形立管在随机波浪、顶部浮体共同作用下的非线性动力响应进行了计算分析，并得到了随机波浪(特征波高、峰值周期)、顶部浮体运动(平均静偏移、慢漂运动周期、慢漂运动幅值)、内部流体(内流密度、内流速度)等对陡波形立管应力的影响规律，主要结论如下：

(1) 在随机波浪、顶部浮体共同作用下，陡波形立管浮子段整体应力水平、应力变化幅值最高；应力极值点位于拱弯点、垂弯点。

(2) 陡波形立管对于随机波浪变化的敏感度较高，悬挂点处应力对随机波浪变化最为敏感，大波高、低周期的波浪会导致立管整体应力水平急剧升高。

(3) 陡波形立管对于顶部浮体运动敏感度明显低于随机波浪，浮体平均静偏移对立管的影响大于顶部浮体慢漂运动；随平均静偏移增大，立管应力极值减小；大幅值、低周期的慢漂运动会引起悬挂点应力显著增大。

(4) 陡波形立管对内流密度、内流速度表现出较高的敏感度，尤其是各关键点处敏感度最高；高密度、高流速的内部流体会使陡波形立管处于高应力状态。