基于离散元与有限元耦合的滑坡作用下输气管道动力响应分析

徐玲，彭星煜，鲁广航

(西南石油大学石油与天然气工程学院，成都 610500)

油气长输管道具有输送距离长、沿途地形地貌多样、外部载荷复杂等特点[1-2]。在沿途地质灾害破坏形式中，滑坡是主要形式也是破坏性极大的一种[3-4]。滑坡对油气管道的安全运行以及诱发管道失效造成的次生危害对自然环境和人类生命财产具有潜在威胁[5-6]，造成的经济损失不可控估量。因此能准确分析滑坡作用下油气管道的动力响应对保障管道安全运行以及维修维护具有重要意义。

国内外学者针对滑坡灾害对油气管道的影响进行了一系列研究，目前研究手段主要包括解析解法、数值模拟以及物理模型三种。Philips等[7]、O’rourke等[8]、Chan[9]充分考虑了管道与滑坡土体间的相对滑动，运用土弹簧模型来表示管道与滑坡体的接触作用，得到纵向、横向两种类型滑坡作用下管道受力分布情况。陈利琼等[10]、王学龙[11]基于数值模拟分别采用CAESAR Ⅱ和ANSYS对埋地输气管道在纵穿和横穿滑坡段下进行了应力分析与对比，并分析不同参数对管道应力应变的影响；朱勇等[12]、Daiyan等[13]通过物理实验建立横向滑坡模型，并结合弹性地基梁构建有限元模型，对比验证了有限元模型的有效性，得出滑坡段管道应力分布情况。王荣有等[14]、王金安等[15]采用理论推导与数值模拟方法结合研究滑坡体的演化规律。当前大多数学者普遍采用数值模拟方法，管道与滑坡土体相互作用模型的选择主要采用非线性接触模型，但该模型未能充分考虑滑坡土体颗粒的离散性，不能准确反映出管道与滑坡体在接触过程中的颗粒与管道、颗粒与颗粒之间的相互作用关系。模拟方法上采用单一的离散元或有限元方法，不能准确地模拟出颗粒与几何结构间的接触作用。

基于此，现提出采用离散元与有限元(discrete element method-finite element method，DEM-FEM)耦合分析法建立管道滑坡模型，从小尺度分析滑坡土体颗粒与管道的接触作用，得到管道在滑坡作用下的受力情况，使模拟过程尽可能地符合工程实际，并通过具体工程案例进行研究分析。

1 离散元-有限元耦合

对于离散元与有限元的耦合方式，根据已有研究，主要有单向耦合和双向耦合两种。通过对比分析，选择单项耦合方式，其示意图如图1所示。首先，利用离散元软件模拟计算出颗粒与几何结构之间产生的作用力；然后，将计算结果作为一种荷载条件导入到有限元软件中；最后，结合其他约束条件对结构进行力学分析。这种方式忽略了几何结构变形对颗粒的位置以及其与几何体的接触产生的影响，因此计算效率方面大大提高，联合仿真也比较容易实现。而双向耦合由于研究案例较少，在技术运用层面上还不够成熟且操作过程较为复杂和困难，对计算机的性能要求极高，因此仅考虑单向耦合。

图1 单向耦合示意图

2 离散元基本方程

2.1 本构方程

离散元本构方程表征作用在两实体上的接触力与相对位移之间的关系。总接触力由法向分量和切向分量组成。力与位移关系表达式为

(1)

2.2 运动方程

运动方程由平移运动方程和旋转运动方程组成，表达式为

(2)

3 工程案例

联络线输气管道起于普光首站，经平昌输气站，终于元坝站，管线全长202 km，管材为X70直缝埋弧焊钢管，规格为Ф1 016 mm×17.5 mm，埋深为2 m，现实际运行压力为7.5 MPa。由于省道409公路升级改造，在斜坡前缘进行开挖，加之施工期间当地降雨较多且较为集中，坡体产生变形滑移，造成管道滑坡事故，现场滑坡情况如图2所示。

图2 滑坡现场航拍图

3.1 基本假设

在对实际工况进行模拟时，结合现场数据采集情况，需对模型做出合理的简化与假设。

(1)在建立管道滑坡模型时，一般认为非滑坡区的宽度应大于滑坡区总宽度的1/4。

(2)假定滑移面为一平面，滑坡坡度为一定值，通过设置挡墙模拟边坡，添加挡墙移动速度来模拟边坡失稳产生的滑坡效果。

(3)假定滑坡体内只含有土体颗粒而不含其他石块等成分，滑坡过程仅考虑土颗粒对管道产生冲击作用，颗粒形状采用不规则的圆球组合模型。

(4)由于管道采用焊接方式连接，为降低计算难度，将焊缝的力学性能视为与母材一致。

3.2 滑坡计算模型建立

3.2.1 管道滑坡模型

根据滑坡输气管道现场情况，以及表1、表2基本模型参数，通过Solidwork建立滑坡输气管道三维模型，其俯视图与侧视图如图3所示。为保证管道在非滑坡区域内与基岩接触充分，实现良好的嵌固，并综合考虑仿真计算效率，建立管道模型总长为32 m，滑坡体宽度为10 m，基岩两侧的非滑坡区域分别为12 m。在滑坡体前人为设置挡墙模拟边坡防护工程，通过设定挡墙移动速度来模拟因边坡失稳后颗粒因失去支撑而产生滑坡的工况。

表1 管道特性参数

表2 滑坡基岩与挡墙几何参数

图3 管道滑坡模型图

3.2.2 土体颗粒模型

实际情况下滑坡土体颗粒形状不规则，颗粒直径大小也不尽相同，综合考虑模型大小、计算效率等因素，采用简单颗粒与复杂颗粒组合模型来模拟滑坡土体颗粒，颗粒模型如图4所示，以便更加符合真实的滑坡体土体性质。颗粒参数参照文献[16]的微观参数标定实验结果数据，如表3所示。采用系统中“Hertz-Mindlin with JKR”模型[17-18]来表征颗粒与几何结构之间以及颗粒间的接触特性。

图4 颗粒模型图

表3 颗粒特性参数

3.3 离散元模拟分析

3.3.1 离散元参数设置

由滑坡体几何模型估算填充满滑坡体所需颗粒，在EDEM软件[19-20]中设置颗粒生成速率为20 000颗/s，生成的颗粒沿重力方向下落对滑坡体进行填充，颗粒填充过程总时长设为51 s。通过试算确定时间步长，在颗粒填充阶段操作时间步设为Rayleigh时间步[21]的30%，在滑坡模拟阶段操作时间步设为Rayleigh时间步的35%。数据存储间隔在颗粒填充阶段设置为1～2 s，在滑坡阶段设置为0.5～1 s。通过虚拟颗粒工厂生成约102万颗直径为20 mm、位置均随机分布的颗粒填充滑坡体。在此基础上，设置挡墙移动速度为沿Z轴正方向3 m/s，滑坡过程模拟15 s，其他设置均保持不变，再选择8线程提交运算。

3.3.2 离散元模拟结果

由离散元模拟得到滑坡过程中土体颗粒与管道不同时刻的接触状态如图5所示，可以看出在滑坡发生3 s后管道前方土体颗粒与管道失去接触，出现一定的临空，此时管道后方土体无法及时移动至管道前方的空隙处，导致管道后部土体向下滑移并挤压管道。随着滑坡的继续，管道前方的临空面逐渐增大，并且逐渐向重力方向偏转，而管道后部土体颗粒因管道外表面的阻挡无法继续滑移，持续推挤管道，直到滑坡结束，管道前方外表面已经完全临空。

图5 滑坡过程中不同时刻颗粒-管道三维接触状态

通过EDEM后处理得到管道在滑坡过程中受到的合力以及各方向受力时程变化曲线如图6所示，可以看出在滑坡发生后管道所受合力出现先增大后减少的趋势，管道在轴向(X方向)上受力较小，可忽略其影响，在重力方向(-Y方向)上由于管道上方的颗粒的减少，管道受到该方向上的作用力随时间逐渐减小，直至管道部分暴露出埋土外，管道受力保持稳定。此外，管道在迎坡方向(Z方向)上受力变化趋势与管道所受合力变化趋势一致。由图可以看出，管道所受最大合力与迎坡方向受到的最大作用力均发生在第3秒，最大合力为 53 686.7 N，重力方向最大受力为53 686.7 N，迎坡方向最大受力为33 822.7 N，管道在滑坡作用下所受最大推力方向与重力方向夹角为36°。

图6 滑坡过程中管道受力时程曲线

3.4 有限元模拟分析

3.4.1 本构模型与参数设置

基岩本构模型选择目前普遍使用的Mohr-Coulomb模型，管材本构模型采用Ramberg-Osgood模型。设置基岩与管道间为摩擦接触，根据研究摩擦因数设为0.3，将管道在第3秒时刻各单元受到的最大滑坡推力结果文件通过插件导入ANSYS，软件自动识别出滑坡推力作用在管道上的空间分布，管道内压设置为7.5 MPa，并且考虑重力场的作用，对基岩底部添加固定约束，其作用面选择滑坡区域管道外表面，完成求解设置后提交运算。

3.4.2 有限元模拟结果

在滑坡作用下，由于滑坡土体颗粒对管道的冲击，由图7(a)可以看出，在滑坡体主滑线中部受到最大等效应力为237.53 MPa；在滑坡过程中，滑坡体与基岩在轴向会形成相对运动，在滑坡边界基岩会对管道产生较大的剪切作用，因此在该处管道受到剪切应力最大，由图7(b)可得最大剪切应力为60.555 MPa；同时由图7(c)可以看出，管道在滑坡体主滑线附近的变形位移最大为12.495 mm，处于非滑坡区域内管道变形位移为零，由此表明非滑坡区管道在基岩中处于嵌固状态，因此可以说明两侧非滑坡区管道长度设为11 m长能够有效实现嵌固作用。由图7(d)滑坡段管道的等效应力分布可知，在滑坡管段中部和滑坡周界处的应力最大。

图7 管道应力分析

3.5 模型验证

通过离散元与有限元的耦合得到滑坡过程中管道受到的最大滑坡推力与应力分布情况，结合管道应力监测平台，平台设有3个监测点，分别位于滑坡两边界处(监测点1、3)和滑坡管道中部(监测点2)。对比3处管道2020年4月4—9日的应力监测数据，如图8所示，在同一滑坡推力下，滑坡管道中部受到的应力最大监测值为224.52 MPa，滑坡边界处的应力最大监测值分别为277.91 MPa和286.84 MPa；有限元模拟值分别为237.53 MPa和299.79 MPa(左右边界对称)。由此可见，现场监测数据和有限元模拟结果得到的应力大小趋势基本一致，且3个监测点的相对误差在7.30%以内，在可接受的范围中，因此验证该模型具有一定的可靠性和参考意义。

图8 监测平台实时监测数据

4 结论

综上分析，可得出以下结论。

(1)通过离散元和有限元单向耦合方式，基于EDEM模拟土体颗粒滑坡过程，得到管道在各方向的受力时程变化曲线，将第3秒时刻管道所受滑坡最大推力作为载荷导入ANSYS中，模拟得到管道等效应力、剪切应力以及位移的分布云图。应力在滑坡管段中部和滑坡周界处的最大，变形位移在滑坡体主滑线附近最大。

(2)将数值模拟得到的滑坡管道中部和滑坡周界处的最大等效应力值与对应位置的监测点数据进行对比，相对误差在7.30%以内，因此表明模型具有一定的可靠性。

(3)根据滑坡段管道等效应力分布，输气管道若不可绕避滑坡体，则建议加强对滑坡管道中部和滑坡周界处，以及环焊缝连接处的实时监测，避免应力集中造成管道加速破坏。