0.   引言

置于海床表面的海底管道是深海油气工程的重要组成部分，其在位稳定性至关重要。管道的在位稳定性受到波流荷载及管道与周围土体相互作用的影响，因此如何准确评估水-土-管耦合作用下的管道在位稳定性是深海油气工程需要关注的关键问题之一，其中存在的难点在于准确评估管-水相互作用、管-土相互作用以及两者间的耦合作用^[1]。

关于管-水相互作用方面的研究，其关键在于如何准确计算作用在管道上的波流荷载。最初的研究主要针对局部流场不受扰动下的波流荷载，Morison等^[2]发现波和流作用在圆柱形物体上时，会产生水平方向的拖曳力、惯性力，及竖直方向上的升力，并提出了相应的计算方法，即Morison方程。然而管道的存在会改变其周围局部流场，流场的变化会进一步改变水动力荷载大小，而Morison方程并未考虑该变化带来的影响。为进一步分析管道周围流场变化对管道所受水动力荷载的影响，Jacobsen等^[3]与Lambrakos等^[4]进行了一系列管-水相互作用试验，基于试验结果将Morison方程中恒定系数修正为关于时间的函数，提出了考虑局部流场影响下的修正Morison方程。Fyfe等^[5]进一步对修正Morison方程中的力系数进行傅里叶分解，并提出了相应的计算方法，通过与试验以及其他计算方法的结果进行对比，证明了该方法更具准确性。

关于管-土相互作用方面的研究，其重点在于如何准确计算管道侧向运动过程中的土体抗力。早期工程中利用库仑摩擦理论计算土抗力^[6]，但计算算结果往往与真实情况间存在较大差异，Lyons等^[7]和Wagner等^[8]通过模型试验对管道在土体中的受力情况进行分析，发现使用库仑摩擦理论得到的土抗力与试验结果存在较大差异。Verley等^[9-10]通过回归分析法对大量试验数据进行分析，提出了工业界广泛使用的Verley管土相互作用模型。Bruton等^[11]及Cheuk等^[12]研究了管道在黏土中发生较大侧向位移时的土抗力大小，得出了相应的计算模型。在理论研究方面，Tian等^[13]基于塑性理论，在二维荷载空间内针对钙质砂建立了土抗力-位移“宏单元”模型。Wang等^[14]在研究管道在砂土中大位移条件下的受力特性时，发现了侧向土抗力随位移增大而出现的软化现象。同样针对大位移条件，Wang等^[15-16]开展了管土模型试验，研究了土阻力与管道初始埋深及管重的关系。Wang等^[17]对护堤的形成机制及其对横向土阻力的影响进行了研究，提出了一种计算护堤土阻力的方法。其中Verley等^[9-10]提出的计算模型所需参数少，适用土质多，工业界广泛采用该模型对海底管道的稳定性进行分析。

准确评估海底管道的在位稳定性需综合考虑水-土-管三者的耦合相互作用。Youssef等^[18]开展了离心机试验，对现有侧向管土相互作用模型进行了完善修正。高福平^[19]针对流-固-土多因素耦合作用下海底管道的失稳机理进行了分析阐述，包括管道悬空触发的渗透侵蚀破坏机制与涡激振动响应。在此基础上，徐万海等^[20]和赵瑞等^[21]分别进行了海底管道涡激振动响应试验装置开发及预报模型建立。在评估管道的在位稳定性时，现有研究通常采用对管道二维断面进行力学分析，然而对于多因素耦合作用下三维海底管道的稳定性分析方法的研究还十分匮乏。为此，本文通过有限元软件ABAQUS中的DLOAD与UEL子程序，搭建了由管-水相互作用计算程序和管-土相互作用计算程序组成的计算分析平台，并分别基于南海一年一遇和百年一遇的海况条件，对海底管道的三维在位稳定性进行了分析。

1.   水-土-管耦合计算分析平台

本文搭建的水-土-管耦合分析平台主要包含管-水作用及管-土作用两部分程序，其中管-水作用程序用来模拟海洋环境并施加水动力荷载，管-土作用程序用以及计算管道运动过程中收到的水平向土体抗力。本文通过有限元分析软件ABAQUS中的DLOAD与UEL子程序将上述两部分程序进行有机结合，形成了水-管-土耦合分析计算平台，如图 1所示。该计算平台由输入文件提供初始运行参数，ABAQUS作为程序接口，将数据参数分别提交给管-水作用与管-土作用程序进行计算，计算结果储存在相关输出文件中。

图  1  计算平台工作逻辑

Figure  1.  Working logic of computing platform

下载: 全尺寸图片 幻灯片

1.1   管-水相互作用计算程序

置于海床土表面的管道，会受到由海浪与洋流联合诱发的水平荷载${F_{\text{H}}}$与竖直荷载${F_{\text{V}}}$，如图 2（a）所示^[22]。海浪复杂且无规律性的波，为对其进行特征分析，使用傅里叶分解法将其分为数个不同振幅、频率、方向和相位的规则波浪^[22-23]，如图 2（b）所示，并使用有效波高与有效周期对海浪进行表征。假设一个海浪由N个波浪组成，将其中波高最大的N/3个波的波高平均值作为该海浪的有效波高，与有效波高类似，有效周期为N个波浪中波高最大的N/3个波的周期的平均值：

图  2  作用于海底管道的水动力荷载

Figure  2.  Hydrodynamic loads acting on submarine pipelines

下载: 全尺寸图片 幻灯片

式中：${H_{\text{s}}}$为有效波高；$N$为组成海浪的波数；${H_i}$为单个波的波高；$ {T_{\text{s}}} $为有效周期；${T_i}$为单个波的周期。

通过下式计算三维海面的高度变化情况：

式中：$\eta $为自由面距平均水面的瞬时高度；${a_{i,j}}$为波浪振幅；${k_i}$，${k_j}$为波浪系数；${\omega _i}$为角频率；$\theta $表示波浪方向；${\varepsilon _{i,j}}$为波浪相位值，代表波浪的能量谱，本文使用Pierson-Moskowitz谱（PM谱）以及基于PM谱提出的JONSWAP谱对能$S({\omega _i})$量谱进行描述；$ D({\omega _i},{\theta _j}) $为定向扩展函数，如果不考虑方向扩散函数，即假设总能量沿一个方向运动，此时的计算结果较为保守。为更准确地计算水动力荷载，在模拟海浪时应采用定向扩展函数。

在波浪运动时，水质点会进行圆周运动，运动的半径随着水深的增加而减小，直到达到波基面，由波浪带来的影响消失。波基面的深度通常等于1/2波长，其范围一般在20~100 m。

为了估算水面以下深度为z处来自波浪的水动力载荷，需计算波浪引起的波浪速度与加速度，计算表达式为

式中：d为总水深；$\alpha $为波浪方向与管道沿线法线的夹角。

作用于管道上的水平水动力荷载包括拖曳力${F_{\text{D}}}$和惯性力${F_{\text{I}}}$，竖直方向的荷载为上升力${F_{\text{L}}}$，使用傅里叶分析方法计算荷载的一般计算公式为

式中：${U_{\text{e}}}$为管道附近的水流速；${U_{\text{P}}}$为管道在水平方向的速度；$ {\dot U_{\text{P}}} $为管道水平运动的加速度；${C_D}$和${C_{\text{L}}}$分别为拖曳力系数和上升力系数；${C_{\text{a}}}$为附加质量系数；$\rho $为水的密度。

本文将上述荷载计算方法通过FORTRAN语言代码化，建立了管-水相互作用计算程序。图 3为使用该程序计算的一组水动力荷载，从图中可以看出惯性力、拖曳力和上升力随时间和管长的变化。

图  3  水动力荷载计算结果

Figure  3.  Calculated results of hydrodynamic loads

下载: 全尺寸图片 幻灯片

本文将程序的模拟结果与Jacobsen等^[3]进行的水动力试验及本程序的计算结果进行了比较，如图 4所示，可以看出两者的变化情况与数值大小均有较强的一致性（图 4中Jacobsen等^[3]通过开展试验获得的水动力数据是通过传感器测量获得的，而由于传感器的存在，其水动力结果无法避免会受到一定干扰，与真实值之间存在一定差异。因此，本文模拟结果与其试验结果会存在一定差异性），从而证明了该程序结果的准确性。

图  4  水动力计算结果比较

Figure  4.  Comparison of calculated hydrodynamic results

下载: 全尺寸图片 幻灯片

1.2   管-土相互作用计算程序

管-土相互作用计算程序基于Verley等^[9-10]提出的模型开发，该模型因其计算结果的准确性，且适用于不同性质的海床土，已被海底管道设计行业广泛使用^[24-25]。Verley等^[9-10]采用了PIPESTAB，AGA和DHI 3个来源的侧向管土作用试验数据，以此为基础进行了量纲分析和回归分析，提出了土抗力与位移的变化模型，如图 5所示。

图  5  Verley管土相互作用模型

Figure  5.  Verley's pipe-soil interaction model

下载: 全尺寸图片 幻灯片

模型假定管道所受来自土体的水平抗力$F$分为被动土抗力$ {F_{\text{r}}} $和摩擦力$ {F_{\text{f}}} $两部分，并认为摩擦力部分不对该模型计算产生影响^[26]。

图 5模型描述了管道位移与土阻力之间的对应关系。该模型主要分为4个区域（关于原点对称）：

（1）$y \leqslant {y_1}$时，处于弹性阶段，在此阶段没有能量变化，管道埋深不发生变化，土体阻力随着管道水平位移线性增加。

（2）${y_1} < y \leqslant {y_2}$时，随着位移的增加，管道埋深开始逐渐增加，土阻力随着管道前方的土堆大小的增加而增大，y₂按经验值取0.5倍管径。

（3）${y_2} < y \leqslant {y_3}$时，土阻力下降的同时，埋深也将在管道位移达到y₃时线性减少至$y = {y_2}$对应埋深值的一半。

（4）$y > {y_3}$时，土阻力和管道埋深均保持恒定。

相关计算公式如表 1，2所示，其中${\gamma '_{\text{s}}}$为土的浮重度，D为管道直径，W_s为单位管长重量，${F_{\text{L}}}$为波浪在单位管长上产生的上浮力，s_u为土的不排水抗剪强度。

表  1  海床为砂土时对应的计算公式

Table  1.  Corresponding formulae for seabed of sandy soil

参数	公式	编号
峰值土阻力	$ {F_{{\text{r2}}}} = \left\{ \begin{array}{l} {{\gamma '}_{\text{s}}}{D^2}(5.0 - 0.15{\kappa _i}){\left( {\frac{{{z_2}}}{D}} \right)^{1.25}}{\text{ (}}{\kappa _i} \leqslant 20){\text{ }} \hfill \\ 2.0{{\gamma '}_{\text{s}}}{D^2}{\left( {\frac{{{z_2}}}{D}} \right)^{1.25}}~~~~~~~~~~~~~~~~~~~({\kappa _i} > 20) \hfill \\ \end{array}{} \right. $	(9)
$ {y_1} $对应土阻力	${F_{{\text{r1}}}} = 0.3{F_{{\text{r2}}}}$	(10)
$ {y_3} $对应土阻力	$ {F_{{\text{r3}}}} = \left\{ \begin{array}{l} {{\gamma '}_{\text{s}}}{D^2}(5.0 - 0.15{\kappa _i}){\left( {\frac{{{z_3}}}{D}} \right)^{1.25}}{\text{ (}}{\kappa _i} \leqslant 20){\text{ }} \hfill \\ 2.0{{\gamma '}_{\text{s}}}{D^2}{\left( {\frac{{{z_3}}}{D}} \right)^{1.25}}~~~~~~~~~~~~~~~~~~~({\kappa _i} > 20) \hfill \\ \end{array}{} \right. $	(11)
初始埋深	$ {z_i} = 0.037D{\left( {\frac{{{{\gamma '}_{\text{s}}}{D^2}}}{{{W_{\text{s}}}}}} \right)^{ - 0.67}} $	(12)
$ {y_2} $对应埋深	$ {z_2} = 0.23D{\left[ {\frac{E}{{{{\gamma '}_s}{D^3}}}{\kappa _i}^{ - 1}{{\left( {\frac{y}{D}} \right)}^{ - 0.5}}} \right]^{0.31}} + {z_i} $	(13)
$ {y_3} $对应埋深	$z_3= \begin{cases}{\left[0.82-3.2\left(\frac{z_2}{D}\right)\right] z_2} & \left(\frac{z_2}{D} \leqslant 0.1\right) \\ 0.5 z_2 & \left(\frac{z_2}{D}>0.1\right)\end{cases}$	(14)
土体刚度参数	$ {\kappa _i} = \frac{{{{\gamma '}_{\text{s}}}{D^2}}}{{{W_{\text{s}}} - {F_{\text{L}}}}} $	(15)
能量	$ E(t) = \int_0^t {{F_r}{\text{d}}s} $	(16)

下载: 导出CSV 

| 显示表格

表  2  海床为黏土对应的计算公式

Table  2.  Corresponding formulae for seabed clay soil

参数	公式	编号
峰值土阻力	$ {F_{{\text{r2}}}} = 4.13{s_u}D{G^{ - 0.392}}{\left( {\frac{{{z_2}}}{D}} \right)^{1.31}} $	(17)
$ {y_1} $对应土阻力	$ {F_{{\text{r2}}}} = 4.13{s_{\text{u}}}D{G^{ - 0.392}}{\left( {\frac{{{z_{\text{i}}}}}{D}} \right)^{1.31}} $	(18)
$ {y_3} $对应土阻力	$ {F_{{\text{r3}}}} = 4.13{s_{\text{u}}}D{G^{ - 0.392}}{\left( {\frac{{{z_3}}}{D}} \right)^{1.31}} $	(19)
初始埋深	$ {z_i} = 0.0071D{(S{G^{0.3}})^{3.2}} + $ $ 0.062D{(S{G^{0.3}})^{0.7}} $	(20)
$ {y_2} $对应埋深	$ {z_2} = 0.296{\xi ^{0.32}}{S^{0.637}}D $	(21)
$ {y_3} $对应埋深	$ {z_3} = 0.489S{G^{0.54}} $	(22)
能量	$ E(t) = \int_0^t {{F_{\text{r}}}{\text{d}}s} $	(23)
其他参数	$ S = \frac{{{W_{\text{s}}} - {F_{L} }}}{{D{s_{\text{u}}}}} $	(24)
	$ G = \frac{{{s_{\text{u}}}}}{{D{\gamma _{\text{s}}}}} $	(25)
	$ \xi = \frac{E}{{{s_{\text{u}}}{D^2}}} $	(26)

下载: 导出CSV 

| 显示表格

随着管道位移的变化，管道模型也会根据位移的变化实时更新。当在0到y₁间的范围内移动时，即使管道在外界荷载作用下运动方向发生改变时，模型也不会发生改变。当管道位移大小超过y₁时，模型原点会向管道移动的方向平移大小$y - {y_2}$的距离，与此同时，F_r2值会降低为当前位移对应的F_r值，其余关键点的值也会随之调整。且若在此时管道发生转向，更新后的模型如图 5所示，该部分内容均由Verley等^[9-10]设定并验证。

本文除了将上述内容代码化以外，为还原管道实际运动状况，还加入了管道上浮判定以及下落位置判定。当水流上浮力大于管重时，程序判定管道进入自由漂浮状态，此时模型将不参与计算，直到上浮力小于管重后再次落入土层；当管道由漂浮状态再次落入土层，若管道落入之前运动形成的沟槽范围内时，会调用管道上浮之前存储的节点状态数值，否则，程序将会以初始模型来重新计算。

图 6将Verley等^[9]通过在砂土海床中模拟得到的管道埋深随时间变化的计算结果与本文程序计算结果进行了对比。由图 6可知，两者在趋势性与数值大小方面均拥有较强的一致性，说明该程序准确反映了原模型。但因在提取位移曲线中的数据作为输入参数时存在误差，故图 6中两曲线无法完全对应。

图  6  计算平台与Verley模型模拟结果对比

Figure  6.  Comparison of simulated results between computing platform and Verley's model

下载: 全尺寸图片 幻灯片

1.3   计算分析平台功能实现

本文主要涉及DLOAD和UEL两个子程序，其中DLOAD是ABAQUS中定义随时间和空间位置或其他变量变化的荷载的子程序，用于生成水动力荷载并施加在管道上；UEL子程序负责自定义ABAQUS中有特殊性质的单元，并可以根据ABAQUS主程序传入的信息更新单元对模型不平衡力以及刚度矩阵的贡献，用于将上述管土作用模型耦合入ABAQUS自定义单元。在该计算平台中，管-水相互作用与管-土相互作用计算程序分别通过DLOAD和UEL子程序整合入ABAQUS。

如图 7所示，在开始分析之前，采用梁单元模拟管道，水动力分布荷载通过DLOAD作用于管道，并通过附加在管单元节点上的UEL单元实现土抗力与管道埋深的计算。

图  7  管道模型的组成

Figure  7.  Composition of pipeline model

下载: 全尺寸图片 幻灯片

计算平台的调动逻辑如下：

（1）分析开始时，ABAQUS通过DLOAD调用管-水相互作用计算程序，根据设定的时间长度与管道模型生成水动力荷载数据，并将这些数据输出保存在数据库中，便于后续调用。

（2）在随后的每个时间增量步中，ABAQUS调用管-水作用程序，根据新的管道位置更新数据。

（3）通过调用水动力荷载计算库将水动力荷载施加在管道上后，在ABAQUS中通过刚度矩阵计算增量位移。

（4）利用当前节点的增量位移，在UEL中更新管土刚度矩阵并传递给ABAQUS。

（5）ABAQUS使用牛顿迭代法，重复步骤（3）和步骤（4），直到收敛。

除此之外，计算平台中还使用了ABAQUS内置接口用户子程序UEXTERNALDB管理外部的数据库，以实现读取，输出与更新等功能，如图 1所示。

为观察模块耦合前后计算结果的差异，在相同参数设置情况下，本文使用耦合后的计算分析平台与管-土计算程序分别计算土抗力变化情况，如图 8所示，虽然二者的计算结果在曲线变化趋势上有较强的一致性，但耦合管-水程序后，在上浮力的影响下，土阻力较耦合前数值大小及正负皆有差异，且存在因管道漂浮而导致土阻力为0的阶段，表明即使位移及其他参数相同，上浮力对计算结果的影响依然比较明显，即管-水程序除了通过水平力影响管道位移外，上浮力对于研究管土作用也十分关键。

图  8  程序耦合前后的结果比较

Figure  8.  Comparison of results before and after program coupling

下载: 全尺寸图片 幻灯片

2.   算例分析

2.1   参数设置

通常海浪的周期为0.5~25 s，为完整反映管道在水动力作用下的运动，本节使用计算平台对海底管道在时长为180 s的波浪作用下的稳定性进行了分析。管道全长1000 m，将其划分为长度5 m的管道单元，共200个，在201个节点上设置UEL单元。假设管道埋置与水深50 m的黏土海床上，管道两端设置为完全固定。具体参数如表 3，4所示。

表  3  波流参数

Table  3.  Parameters of waves and currents

波浪类型	波谱	波速/ (m·s-1)	波浪周期/ s	水流流速/ (m·s-1)
不规则波	Jonswap	2.7	4.3	0

下载: 导出CSV 

| 显示表格

表  4  管土参数

Table  4.  Parameters of pipelines and soil

管道直径/m	管长/m	管道浮重度/kN·m-3)	土的类型	土不排水强度/kPa	管土摩擦系数
1.0	1250	10.0	黏土	3	0.6

下载: 导出CSV 

| 显示表格

2.2   计算结果展示与分析

本节选取了南海在一般海况（一年一遇）与极端海况条件（百年一遇），如表 5所示，分析了海底管道在水动力荷载作用下的稳定性。

表  5  南海海况参数

Table  5.  Sea state parameters in South China Sea

类型		环境参数	一年一遇	百年一遇
波浪	有效波高/m		7.3	12.3
	有效周期/s		11.3	14.5
洋流	洋流流速/(m·s-1)		1.13	2.33

下载: 导出CSV 

| 显示表格

图 9为在一般海况条件下，管道模型20号节点在竖直方向上受力变化情况。管重为0.8 kN，当管道所受竖直荷载大于管重时，会在上升力的作用下克服自身重力上浮，管道失去稳定性。极端海况条件下的平均竖直荷载大小约为一般海况条件下的2.4倍，即在极端荷载条件下，管道更易出现失稳现象，如图 9（a），（b）所示。

图  9  管道竖直方向力-时间关系曲线

Figure  9.  Vertical force-time relationship curve of pipelines

下载: 全尺寸图片 幻灯片

图 10为在一般海况条件下，管道模型^#20节点在水平方向上受力变化情况。管道所受水平荷载为拖曳力与惯性力之和，当管道所受水平荷载大于土抗力时，管道会克服土体抗力发生侧向运动，管道失去稳定性，即在极端荷载条件下，管道更易出现失稳现象，如图 10（a），（b）所示。

图  10  管道水平方向力-时间关系

Figure  10.  Horizontal force-time relationship curve of pipelines

下载: 全尺寸图片 幻灯片

3.   结语

本文介绍了用以研究海底管道在位稳定性而搭建的水-管-土耦合分析计算平台，提供了一种分析海底管道在时域中的三维稳定性的数值计算方法。该计算平台共由管-水作用与管-土相互作用两部分程序组成，其中管-水相互作用程序采用傅里叶分析法计算作用于管道上的水动力荷载，管-土相互作用程序采用Verley模型计算管道在运动过程中受到的土抗力。随后使用ABAQUS自定义荷载子程序DLOAD和自定义单元子程序UEL将上述两部分进行了耦合，形成水-土-管分析计算平台，最后利用该计算平台，分析了在南海一般海况条件与极端海况条件下的管道在位稳定性。

同时需要说明的是，由于海床地形变化或者极端海况下海水的冲刷作用，海底长输管道的部分管段可能悬空于海底，发生悬跨现象。一方面，悬跨段的管道周围缺少对其起到约束作用的土体，影响海管稳定性；另一方面，海流流过悬跨段时，易发生由于漩涡脱落产生的涡激振动现象，引发海管共振，同样会对海底管道的稳定性产生影响^[17-18]。将在未来的研究中着重考虑该现象，使本文计算平台可以更全面地分析海底管道的在位稳定性。