向量式有限元法在管土相互作用中的应用

王 飞，李效民，马芳俊，郭海燕

（1.山东科技大学 土木工程与建筑学院，山东 青岛266590；2.山东科技大学 山东省土木工程防灾减灾重点实验室，山东 青岛266590；3.中国海洋大学 工程学院，山东 青岛266100）

0 引 言

随着深海油气资源的开发，钢悬链线立管（SCR）已经成为开发深海油气资源的主要立管形式。钢悬链线立管触地点承受较大的弯曲应力，与海床接触区域易发生疲劳破坏。因此对于SCR 触地段与海床间的相互作用分析越来越受到人们的重视。

目前，国内外专家对立管与海床间的相互作用进行了一些研究，主要有有限元法、有限差分法或通过实验研究的方法。Pesce 等[1]将海床等效为线性弹簧模型，分析了触地点处海床硬度对立管运动的影响。2H 海洋工程公司发起了工业联合开发计划（STRIDE JIP）[2]，对管土相互作用进行了全尺寸模型试验。You[3]通过有限差分法对管土相互作用进行了分析。Aubeny 等[4]根据小尺度模型实验提出了基于P-y 曲线的管土作用模型以及各个阶段的经验公式。Hu 等[5]通过试验表明立管反复作用在海床上会导致土体刚度衰减。郭海燕等[6]应用ANSYS 中的接触单元建立模型对立管-海床的作用进行了分析。王坤鹏等[7]用ABAQUS 创建触地单元分析管土相互作用。Wang 等[8]通过室内大比尺模型试验研究触地区的埋置深度对管内应力的影响。白兴兰等[9]基于改进的程序CABLE3D RSI 对非线性管土作用下钢悬链线立管触地区疲劳进行了分析。李敢等[10]用ANSYS 对管土相互作用引起的疲劳进行了分析。

尽管人们对于管土相互作用已经做了大量工作，但为保证立管在位运行的安全性，在设计分析过程中运用多种数值计算方法和模型进行验证是极其必要的。由丁承先教授[11]提出的向量式有限元法，又称为有限质点法，以向量力学作为分析的理论基础，并与数值计算相结合，是求解结构的大变形、大变位、弹塑性、碰撞、倒塌等非线性或不连续性力学行为的新方法。本文将其应用于钢悬链线立管触地段与海床土体相互作用的分析研究中。相较于传统分析方法，向量式有限元是对立管的力学行为建立计算理论，计算中无需集成结构的刚度矩阵，可方便更改单元和边界条件，基于中央差分的显式积分求解方案无需迭代求解非线性方程组，从而能够降低立管分析的难度，同时分析初始时刻便可以将荷载直接加载到立管上，即使是对静力问题的分析，亦是通过逐步施加外部载荷并不断改变立管的位置和形态直至达到最终结果，因此可以更好地模拟真实海洋环境下立管的运动变形过程，使分析相对准确。

本文基于向量式有限元法采用弹簧模型模拟海床，将悬链线立管触地段离散成一系列质点的集合，质点间用只承受内力而无质量的平面弯曲杆件元连接，采用一系列途径单元描述立管的受力变形过程，通过虚拟的逆向运动计算立管的纯变形和内力。采用牛顿第二定律来表达质点的运动行为，并用中央差分的显式积分法求解质点各时刻的控制方程，编制相应Matlab 求解程序，验证该方法的可行性并分析管土相互作用。

1 钢悬链线立管触地段管土相互作用分析

1.1 弹簧模型

钢悬链线立管触地段与海床相互作用可用一系列弹簧模拟海床，从而可以将P-y 曲线[4]简化为三种弹簧模型，如图1 所示分别为线性弹簧模型、非线性弹簧模型和线性折断弹簧模型。线性弹簧模型假设海床土体是线性的，立管触地段所受土体抗力与土体刚度和入土深度成正比；非线性弹簧模型考虑了土体的最大土体抗力与土体最大吸引力；线性截断弹簧模型考虑到管土分离后，不存在土体吸力。

图1 弹簧模型Fig.1 Spring model

海床初始嵌入骨架曲线为：

其中：Su0为海床表面抗剪强度，Sg为抗剪强度变化梯度，a 和b 随着沟槽形状和管道的粗糙度变化，y为沟槽深度，D 为管道直径。

定义最大土体吸力与极限承载力的比值：

1.2 向量式有限元分析模型

如图2 所示，钢悬链立管触地段管土模型可以简化为立管放置在一系列弹簧上，采用向量式有限元法进行分析。立管左端自由，右端铰接，立管触地段整体坐标系定义如下：取自由端为坐标原点，x 轴沿立管触地段轴线，y 轴为水深方向。采用N+1（从左向右编号依次为1，2，…，N，N+1）个等间距的质点模拟立管的直线状态，质点之间利用平面弯曲杆单元相连接，单元的质量平均分配到两端的质点上，立管的初始张力通过施加单元的初始轴力来实现。对触地段与海床土体的静力相互作用展开研究，通过对触地段左端逐步施加竖向位移代表钢悬链立管悬垂段的变形运动过程，可以使1 号节点在R 时间内运动到指定点，并在R 时间内将重力、浮力及土支撑力以斜坡函数加载的方式作用在各质点上，总共作用T 时间，使管道达到平衡状态。

图2 管土相互作用模型及坐标示意图Fig.2 Diagram of the riser deformation and coordinate

1.3 质点控制方程

对于平面问题，质点有三个位移分量，亦即两个平移量和一个转动量，考虑阻尼的影响，质点控制方程式为：

其中：Mi为i 质点的质量，为i 质点加速度，分别为i 质点所受的外力、内力和阻尼力。其展开为：

其中：mi是质点i 的质量；Ii是其对z 轴的质量惯性矩；βi是转角；ζi为立管的阻尼系数；分别为质点所受到的沿x、y 轴的外力及绕z 轴的外力偶矩，分别为质点所受到的沿x、y 轴的内力及绕z 轴的内力偶矩。

1.4 控制方程的求解

采用中央差分的显式时间积分法，质点的速度和加速度可表示为：

代入（4）式得：

1.5 单元内力计算

取一个标准途径单元tn-1≤t≤tn，杆件αβ 在tn时刻的两端点位置向量是，方向向量为：

以tn-1时刻单元内力和变形作为基础架构，通过虚拟的逆向运动计算立管下一时刻内力和变形纯增量，如图3（a）所示，然后根据材料力学的理论求得杆件元αβ 在tn时刻单元内力为：

最后通过虚拟的正向运动如图3（b）回到原先tn时刻位置。两个节点内力分量分别转换成域坐标分量和，弯矩不需转换，将各单元节点力作用在各质点，再对质点内力进行集成，得到：

当质点是αβ的起始点时τ=1，否则τ=2。

1.6 质点所受外力

立管质点所受外力包括重力、浮力以及海床土体支撑力。表达式可以写为：

其中：ω 为立管湿重，k 为土体刚度。

2 数值求解

基于上文所述理论及计算公式，采用Matlab编写求解程序， 分析管土相互作用， 并与用ABAQUS[12]软件分析的结果和You[3]基于有限差分法编制的程序进行对比分析，计算流程如图4所示。本文模型参数采用STRIDE JIP 实验[2]中的模型数据，如表1 所示。

图3 单元虚拟运动示意图Fig.3 Fictitious motion of the element

图4 向量式有限元法计算流程Fig.4 Flow chart of pipe/seafloor interaction analysis using the VFIFE

表1 模型参数Tab.1 Model parameters

根据文献[5]取a=6.73，b=0.15，沟槽深度y=D（D 为立管外径），则土体极限承载力为Pmax=3.96 kN/m，取m=0.5，则最大吸力为：Pt=0.5×Pmax=1.98 kN/m。在钢悬链立管触地段左端施加竖向位移D，来代替悬垂段运动对立管触地段的影响。因为主要变形发生在靠近自由端的部分，故图像只取靠近自由端的30 m。

图5 是分别采用线性弹簧模型、非线性弹簧模型以及线性截断弹簧模型通过向量式有限元法(VIFIFE)、有限差分法(FDM)和ABAQUS 得到的位形变化对比图。从图中可以看出三种方法得到的三种弹簧模型下的位形结果高度吻合，管道在靠近自由端的位置处向上位移，没有渗入海床，管土发生分离，随着远离自由端管道入土深度逐渐加大，直到达到最深处，然后逐渐变小趋于稳定。图6 是分别在三种模型下用三种方法模拟得到的管道弯曲应力分布图，弯曲应力在管道触地点急剧增大到最大值，随后迅速变小，在远离自由端后趋于平稳，数值很小，几乎可以忽略。从图5-6 两对比图可以看出，三种方法吻合良好，证明向量式有限元法在管土相互作用模拟分析中是可行的。

图5 三种弹簧模型位形对比图Fig.5 Comparison of configurations for three spring model

图6 三种弹簧模型弯曲应力对比图Fig.6 Comparison of bending stresses for three spring model

图7 线性截断弹簧模型下立管变形过程Fig.7 Riser deformation based on tension cut-off spring model

图8 线性截断弹簧模型下立管弯曲应力变化过程Fig.8 The variation of bending stresses based on tension cut-off spring model

图7 是采用线性截断弹簧模型模拟海床时立管触地段位形变化过程图。触地段初始状态为静态直线状态，然后在触地段左端逐步施加竖向位移使管道从水平位置移动到D 位置，最终管道达到静止状态。因此向量式有限元法分析过程是一个动态的实现过程，是对立管力学行为过程的模拟，更加符合立管的真实变化过程；图8 是与之对应的弯曲应力变化过程。相比于基于数学模型的有限差分法及其他有限元软件，向量式有限元法基于物理模型，计算中无需集成结构的刚度矩阵，也无需迭代求解，可方便更改单元和边界条件，从而简化了立管分析的难度。即使是对于静力位形及内力的计算亦是通过逐步动态实现的过程，变形及受力符合立管实际受力情况，所得结果相对精确。

3 参数分析

3.1 弹簧模型的影响

图9 是采用向量式有限元法模拟三种不同弹簧模型得到的位形图。从图中可以看出线性弹簧模型沟槽最深处最接近自由端，因为线性弹簧模型是线性的，没有考虑管土分离，所以立管一直承受海床土体吸力，沟槽最深处距离自由端最近，所取得的最大弯曲应力值也是三种模型中最大的，如图10所示。相比之下，非线性弹簧模型虽然也未考虑管土分离，但假定了最大土体吸力以及最大抗力，所以沟槽最深处离自由端远些；线性截断弹簧模型考虑了土体最大抗力以及管土分离，当土体吸力达到最大值时，立管不再受吸力，所以沟槽最深处距自由端最远，其最大弯曲应力也是三者中最小的。可以看出线性截断弹簧模型更符合海床土体实际情况。

图9 不同模型的立管位形Fig.9 The configuration of the riser with different spring model

图10 不同模型的立管弯曲应力Fig.10 The bending stress of the riser with different spring model

图11 不同吸力系数的立管位形Fig.11 The configuration of the riser with different soil suction force coefficient

图12 不同吸力系数的弯曲应力Fig.12 The bending stress of the riser with soil suction force coefficient

3.2 海床土体吸力系数的影响

采用更符合实际情况的线性截断弹簧模型模拟海床，在立管触地段自由端施加向上位移D，改变吸力系数m，分别取0、0.25、0.5、0.75 和1，得到不同立管位形和应力，如图11-12 所示。从图中可以看出随着吸力系数m 的增加，触地段所形成的沟槽最深处越靠近自由端，同时弯曲应力也相应变大。说明海床土体吸力增大了立管触地段的弯曲应力。

3.3 海床土体刚度的影响

在立管触地段自由端施加向上位移D，改变海床土体刚度，分别取100、200、300、400 和500 kN/m/m，得到不同刚度下立管位形和应力，如图13-14 所示。随着海床土体刚度的增大，沟槽最深处距离自由端越远，相应的弯曲应力也逐渐变小。当土体刚度小时，其对应力的影响很明显，但当刚度变大时，对最大应力影响变小。说明土体刚度变大会使沟槽最深处远离自由端，应力变小，土体刚度变大达到一定值后，应力值应该趋于稳定，接近刚性海床。

图13 不同海床土体刚度的立管位形Fig.13 The configuration of the riser with different soils stiffness

图14 不同海床土体刚度的弯曲应力Fig.14 The bending stress of the riser with different soils stiffness

4 结 论

本文采用弹簧模型模拟海床，将向量式有限元法用于钢悬链线立管触地段与海床土体相互作用的模拟中，并与ABAQUS 及有限元差分法模拟的结果作对比，并分析了不同参数对立管触地段与海床相互作用的静力影响。通过本文研究可以得出以下结论：

（1）三种方法模拟的结果吻合，表明向量式有限元法在管土相互作用的模拟中是可行的；

（2）即使是对于管土相互作用的静力问题分析，向量式有限元亦是动态实现的过程，即由静止开始到运动再到静止结束的过程。分析过程中可以任意改变荷载，可以更好地模拟真实海洋环境下立管的运动行为，采用中央差分的显式时间积分求解，不必形成刚度矩阵，避免了迭代求解和收敛问题，计算中容易改变单元与边界条件，计算程序中只存在质点位置计算和单元内力计算两个循环，整体步骤呈系统化；能够直接求解立管的变形和内力，避免了由变形计算内力的复杂过程和可能的误差；

（3）向量式有限元法可以有效地分析不同吸力系数和海床刚度时的管土相互作用，能够很好地契合三种模型，但还需进一步验证该方法在循环作用下土体刚度的衰减P-y 曲线模型中的应用。