0.   引言

海底管道是海洋油气能源运输的生命线,其安全运营是保障海洋油气资源开发的关键之一。深海管道在安装过程中通常直接铺设在海床表面,在管道自重、波浪、洋流等荷载和施工船舶扰动等动态铺设效应的共同作用下贯入海床至一定深度。管道动态运动一方面引起土体不排水抗剪强度降低,造成土体重塑,另一方面使得海床表面形成沟壕,两者共同作用造成管道贯入深度大于准静态计算值。管道埋置深度的准确评估对深海管道的设计至关重要,特别是在管道的在位稳定性、疲劳寿命^[1]、热胀屈曲和轴向移动^[2-3]等几个方面。

对于管道在不排水黏土中的静态竖向贯入过程的模型试验^[4]已有大量相关研究。受土体参数限制,对于更为一般的问题,往往采用经典塑性理论^[5-7]或小应变有限元分析理论^[8-9]进行分析。近年来,随着数值模拟技术的高速发展,大变形分析开始得到广泛运用。当前,最常见的大变形模拟方法包括基于有限元分析的任意欧拉-拉格朗日法（ALE法）^[10]、耦合欧拉-拉格朗日法（CEL法）^[11-12]和小应变下网格重划分与差值方法（RITSS法）^[2],但采用上述方法的研究对管道贯入深度多局限于一倍管径内。Hawlader等^[13]建立了计算流体动力学模型和CEL模型对单调贯入至5倍直径的管道运动过程进行模拟,真实反映了由管道贯入引起的极端土体变形,但由于网格重划分涉及到复杂的应力映射,该模型未模拟土体的坍塌回流问题。在管道贯入海床的过程中,管周土体将发生严重重塑并坍塌至管道上部。这一过程中土体自由界面表面演变复杂,对数值模拟带来了极大的挑战。因此,开发一种稳定高效的数值模拟方法是解决此类问题的关键。

本文基于连续极限分析法（SLA）对管道在竖向循环运动过程中的土体大变形问题进行了研究。首先分析了单调和循环位移下管土相互作用特性,并与已有的模型试验结果进行了对比验证。之后模拟了不同复杂加载路径下管道与土体的循环竖向作用,探索了土体抗力和破坏机构演变规律。该工作以期为复杂荷载条件下管-土相互作用简化模型的建立提供探索性依据,同时推广SLA方法在极端大变形岩土工程问题中的应用。

1.   数值模型

1.1   SLA法介绍

该方法在早期通常用于解决结构工程领域的一些简单问题,例如框架和平板等,此后被Kong等^[14-15]进一步拓展使其适用于土结相互作用的分析。SLA方法以牛津大学开发的有限元极限分析软件OxLim为核心,使用极限上限分析求解塑性变形区域的速度场,并将其在单个时间步内进行积分从而更新模型的瞬时位移场。通过重复这一小应变极限分析步骤,最终获得大变形问题的求解。

在SLA数值分析中,土体本构模型通常采用改进后的Tresca模型。Einav等^[16]在大量研究基础上提出通过如下经验公式考虑应变软化和应变率效应对土体不排水抗剪强度的影响：

式中,s_u0为土体的初始不排水抗剪强度,一般情况下由s_{u0 =} s_um + kz确定,s_um为泥面强度,k为强度沿深度变化梯度,z为泥面以下深度；第二项应变软化参数用于考虑由塑性应变$\xi$累积造成的土体强度衰减,第三项率效应参数则用于说明应变速率对土体强度的影响。残余强度系数${\delta }_{\text{rem}}$对应为灵敏度$S_{\text{t}}$的倒数,表示土体完全重塑后强度与初始强度的比值,延展参数${\xi }_{\text{95}}$则与土样强度降低到其峰值强度5%时的累积剪应变相对应。其中,海洋软土的${\delta }_{\text{rem}}$变化范围一般为0.2～0.5,${\xi }_{\text{95}}$为10～50。黏度系数$\mu$用于分析土体强度随应变率增大而增强的关系,${\dot{\gamma }}_{\text{max}}$为土体最大剪应变速率,${\dot{\gamma }}_{\text{ref}}$为土体参考剪应变速率,对于海洋软土而言,μ取值范围通常为0.05～0.2,${\dot{\gamma }}_{\text{ref}}$一般取3×10⁻⁶/s^[16]。

1.2   管土模型及参数选取

本文基于平面应变假定进行分析,图1为管道-土体相互作用模型及初始网格划分,管道直径为D,其属性设置为平面应变刚体,整个圆周划由60条直线段连接而成；管道在竖直方向上可自由移动,仅对其旋转自由度进行约束。土体的半宽、半深分别用B和L表示,底边和右侧边界设置为固定约束,表面设置为自由边界。由于该计算模型具有对称性,本文对其右半部分进行考虑,并对对称轴上的节点的横向移动进行约束。模型的最小单元尺寸为0.0005D²,每次分析的管道增量位移δd设置为0.005D。

图  1  计算模型

Figure  1.  SLA model

下载: 全尺寸图片 幻灯片

土体浮重度为${\gamma }^{\prime }$,管土界面切向最大容许剪应力为${\tau }_{\mathrm{int}}$=$\alpha s_{\text{u}}$,其中$\alpha$表示界面黏聚力系数（$\alpha$=0为完全光滑,$\alpha$=1为完全粗糙）,s_u为土体不排水抗剪强度。管土法向最大拉应力为T,在贯入过程中假定管土界面不能承受拉应力（T=0）,而上拔过程中管道底侧承受无穷拉应力（T=∞）。

2.   单调加载

本节采用SLA方法对较深埋深管道的竖向贯入过程进行数值模拟,Lee^[17]在常重力下针对高岭土地基进行了此类问题的小比尺物理模型试验,从中观察到土体由整体破坏到局部破坏的破坏机构演变。本节将模型管道视作刚体,由海床表面贯入至4D处,假定管土界面完全光滑,故本算例中$\alpha$=0,具体土体参数见表1所示。

表  1  竖向贯入模型试验参数

Table  1.  Parameters of model tests on vertical loading

参数	管道直径/m	土体强度$s_{\text{u0}}$/kPa	有效重度${\gamma }^{\prime }$/(kN·m-3)	埋深w/D	加载速率vp/(m·s-1)	灵敏度St	延展参数${\xi }_{\text{95}}$
Lee[17]	0.03	1.1+4.5z	6	7	0.003	1.7	10
Hu等[18]	1.0	5.19z	—	3	0.001	—	—

下载: 导出CSV 

| 显示表格

图2为管道从海床表面贯入至4D过程的荷载-位移关系曲线。结果表明：不考虑应变软化和率效应的理想土（S_t = 1,$\mu$= 0）在整个加载过程中与模型试验结果存在明显差异,考虑土体的软化而不考虑率效应的算例（S_t = 1.7,$\mu$= 0）可以更合理地预测了深埋处的阻力,但在深度小于3D时仍存在计算结果偏小的情况；当考虑率效应（$\mu$= 0）时土体抗力计算结果出现明显的增大,在埋深小于1D时与模型试验结果高度吻合。

图  2  深埋管道受荷响应

Figure  2.  Loading response of pipe during deep penetration

下载: 全尺寸图片 幻灯片

图3为管道贯入至一定深度后土体速度场和无量纲化剪切应变率分布情况。随着管道贯入深度逐渐增大,最终剪切破坏带集中分布在管道周围,土体破坏机构转变为局部破坏模式。由于该模型试验是在常重力下以较小比例进行,土体重量的影响远小于现场,因此在模型试验和数值模拟中均未发现管道上方土体的回填现象。

图  3  土体破坏机构（w/D = 4）

Figure  3.  Failure mechanism of soils

下载: 全尺寸图片 幻灯片

3.   循环加载

本节对复杂加载路径下管道与土体的循环竖向作用进行分析。首先对循环位移下管土相互作用特性进行模拟,并将模拟结果与Hu等^[18]提供的离心试验结果进行了比较,验证了SLA模型的有效性；之后针对几个关键参数开展参数分析,对复杂加载路径下土体抗力和破坏机构演变规律进行探究。

3.1   验证

采用SLA方法对Hu等^[18]离心机试验的最初几个加载循环进行模拟。该试验中,管道从泥面处贯入至初始埋深3D,之后进行一次振幅为3D的上拔/再贯入循环。管-土界面张拉力在上拔过程中设为无穷大,贯入过程中为零,其他具体土体参数见表1所示。注意,由于Hu等^[18]给出的初始土体强度在泥面处大小为0,为保证模拟过程数值稳定性,将土体初始强度设置为s_u0=1+4.8z,这一处理使SLA模型中3D深度处的初始土体强度与离心试验保持一致。离心机试验中未提供的参数S_t和${\xi }_{\text{95}}$分别取3.3和10,土体有效重度${\gamma }^{\prime }$取6 kN/m³,与实际深海海床土体重度一致。

图4（a）将数值模拟结果（$\mu$=0.05）与离心机试验得到的土体抗力进行比较。在初始贯入阶段两种方法的计算结果吻合程度较好,在下行到达3D深度、阻力到达最大值前贯入阻力与埋深基本成线性增加关系；稍有不同的是,在开始贯入前离心试验中的管道所受的土体抗力为非零值,该抗力是否源于浮力尚不明确。在上拔阶段,由于假设土体为刚塑性材料且管-土界面张拉力在上拔过程中为无穷大,管道在向上运动初始阶段土体抗力会立刻达到峰值,而在离心试验结果中则需要大约0.3D的距离以达到峰值阻力,但在管道上行至2.4D深度后重新逼近离心试验结果。在再贯入阶段,数值模拟结果与试验曲线吻合较好,土体抗力较初次贯入时有明显减小,表明土体内部发生了严重的重塑现象。图4（b）显示了率效应对土体抗力的影响,此处仅给出初始贯入过程和上拔阶段的结果。在初始贯入阶段,$\mu$=0.05和0.1的数值模拟结果均与离心机试验结果高度吻合；而在上拔过程中土体抗力随$\mu$的增大而显著增加,特别是在埋置深度大于1D时,与离心机试验结果存在较大差距。

图  4  离心试验与SLA模拟竖向循环运动结果比较

Figure  4.  Comparison between centrifuge tests and SLA simulations

下载: 全尺寸图片 幻灯片

图5展示了竖向循环加载试验中不同阶段的土体变形以及应变软化系数和率效应系数分布情况。此处,率效应系数与最大剪切应变率${\dot{\gamma }}_{\text{max}}$相关,也间接反映土体剪切变形程度。管道上拔后管周剪切带附近土体强度接近于完全重塑土强度,远离剪切带上方的土体则经历不同程度应变软化,因此,在再贯入过程中的竖向土体抗力远小于初始贯入过程中的竖向土体抗力（图4（a））。

图  5  竖向循环运动下土体变形

Figure  5.  Deformations of soils during cyclic vertical loading

下载: 全尺寸图片 幻灯片

3.2   参数研究

本节对管道与土体的循环竖向作用中循环振幅w_ult和加载路径的影响进行探究。土体初始强度设置为s_u0=1+2z,土体重度为3 kN/m³。考虑土体的软化而不考虑率效应,残余强度系数${\delta }_{\text{rem}}$和延性参数ξ₉₅分别取0.3和20,其他参数与3.1节中一致。

图6为管道竖向贯入至一定深度（w_ult）后上拔至海床表面过程中土体抗力的变化规律。对4组不同深度进行分析：w_ult/D=0.5,1.0,1.5,2.0。在加载过程中4种工况下土体抗力变化具有相似的规律。值得注意的是,w_ult/D=0.5和1两种工况下管道上拔至泥面处时土体抗力接近于0,这与海洋立管的典型管土相互作用模型吻合,而其它两种工况则表现出大于0的土体抗力。

图  6  不同振幅下竖向循环受荷响应

Figure  6.  Penetration and extraction of a pipe with different amplitudes

下载: 全尺寸图片 幻灯片

图7为非横幅循环加载路径下的土体抗力-埋深关系曲线,上述w_ult/D=1的算例在图中作为参照,将该曲线定义为骨干曲线。在管道每一次再贯入过程的初始阶段,受先期管道运动的影响,土体抗力最初与骨干曲线值相比较小,当管道继续向下贯入,超过历史最大埋深后,土体抗力与骨干曲线逐渐逼近,这种现象与Randolph等^[19]的研究结论一致。此外,图7中前2次循环的管道上拔承载力与骨干曲线值相比偏大,直到第3次循环（0.2≤w/D≤1.0）其值低于骨干曲线,表明管道移动距离越大,土体应变软化越严重,对管道的上拔承载力的削弱作用越明显。值得注意的是,在第三次上拔过程的最后阶段（w/D≤0.2）,循环贯入工况下管道承受的土体抗力明显大于直接贯入的情况,且当管道上升至泥面处时土体抗力为非零值,该现象与图8中两种工况下软化系数相吻合。观察图8可以发现,在循环贯入—上拔试验中管道底部存在土体黏结的情况,而在单次贯入—上拔的工况下管道几乎完全脱离土体。

图  7  非恒幅竖向循环受荷响应

Figure  7.  Retreating/encroaching behaviours during cyclic vertical loading

下载: 全尺寸图片 幻灯片

图  8  土体软化系数

Figure  8.  Softening factor

下载: 全尺寸图片 幻灯片

图9（a）为管道小振幅（0.1D）竖向循环加载路径下的土体抗力-埋深关系曲线。管道在初始状态贯入至0.5D后在0.4D和0.5D之间进行十次垂直振动循环,然后再贯入至1D深度,最后上拔至海床表面。此外,图中同样给出管道直接贯入到1D然后上拔至泥面的骨干曲线以作对比。图9（a）表明,在循环运动过程中土体竖向抗力随循环次数增加而逐渐下降,尽管这一弱化速率随振次逐渐降低,当管道结束循环再次向下贯入时,土体抗力在初始阶段远小于骨干曲线值,随着深度增加逐渐逼近骨干曲线。图9（b）为管道在较大振幅（0.5D）竖向循环加载路径下的土体抗力-埋深关系曲线。在此工况下管道首先被直接贯入至2D,处于完全掩埋状态,然后在1.5D～2D进行5次垂直振动循环,管道在结束循环后最终贯入至3D处,其现象与图9（a）相似。图9（a）和9（b）涵盖了全埋管道和部分埋入管道的贯入/上拔荷载特性,展示了SLA在复杂加载路径下模拟循环管土相互作用的能力。

图  9  非恒幅竖向循环受荷响应

Figure  9.  Loading responses of a pipe subjectal to vertical oscillations

下载: 全尺寸图片 幻灯片

4.   结论

本文采用连续极限分析（SLA）方法对深海管道的竖向管土大变形相互作用进行了研究,针对不同复杂位移路径和土体变形的工况进行了详细分析,得到以下3点结论。

（1）SLA法通过求解一系列连续小变形问题实现对极端大变形问题的模拟,该过程中通过网格重划分与场变量映射避免了网格失真问题。本文中SLA法对管道位移过程中涉及的破坏机构演变、土体沟槽、土坡坍塌回流和上拔破土行为均能进行有效模拟,且所得管土抗力与模型试验结果高度吻合。

（2）SLA因无需构建应力应变关系从而可采用较长的时间步长分析大变形问题,因此具有较高的数值稳定性和计算效率。例如,本文中竖向循环贯入涉及9倍直径（图4（a））的管道位移和土体极端大变形,利用SLA法在常用计算机上可在2 h内完成模拟。

（3）管道在竖向循环运动过程中土体发生应变软化使得管道竖向抗力峰值随循环次数增加而逐渐减小,但埋深突破历史最大埋深时逐渐逼近骨干曲线,表现出典型的迟滞现象,与以往模型试验结果高度吻合。采用SLA法可标定Randolph等^[19]所提出的管-土相互作用非线性弹簧模型,以进行简化三维模拟分析管道疲劳问题。

本文首次针对管道自泥面贯入后的竖向循环管土相互作用展开了研究,介绍了连续极限分析SLA法并得到了模型试验的验证,同时开展典型复杂工况下的模拟展示了SLA在模拟极端土体大变形问题方面的特点和优势。由于文章篇幅限制,未进一步开展更细致的参数分析以建立简化的管-土竖向作用弹簧模型。此外,本文因采用总应力法,未考虑管道贯入过程中所产生的超静孔压,也是本文的局限性之一。