基于SPH方法的埋地PE燃气管道挖掘破坏数值模拟研究

唐宇峰,史君林，李 涛

(1.四川轻化工大学 机械工程学院,四川 宜宾 644000；2.过程装备与控制工程四川省高校重点实验室，四川 宜宾 644000)

0 引言

随城镇化快速发展，诸如道路施工、地基开挖、基础夯土等第3方作用导致埋地管道损伤或破坏事故屡见不鲜，造成严重经济损失。其中，挖掘破坏导致管道损伤与破坏事故最为严重，且比例最大。目前，针对第3方挖掘造成管道破坏研究方法主要为有限元法：蓝朝逊等[1]将管道破坏分为接触、屈服、挖裂和挖穿4个阶段，并采用有限元法分析管道弯头在不同阶段力学响应；李军等[2]采用有限元分析基于应力与应变2种失效准则在管道失效过程中适用性；周立国等[3]采用ABAQUS分析挖掘载荷作用下燃气管道力学特征及动力响应；Neacsa[4]模拟分析不同挖掘及影响因素对管道损失影响。在挖掘作用下，管道受土体、斗齿等多方面综合作用力，出现裂纹扩展或破坏大变形，而基于网格的有限元法研究大变形问题时会产生网格畸变，即使采用网格重构技术，也存在计算费用昂贵，精度低甚至计算终止等问题。

光滑粒子流体动力学(Smoothed Particle Hydrodynamics,SPH)是1种新型纯拉格朗日无网格计算方法，在诸如流体、岩土、爆炸与冲击等计算力学领域得到广泛应用[5-11]，粒子间无网格连接，避免因大变形导致的网格畸变，以表征大变形过程实时运动状态。但在许多领域SPH仍处于起步阶段，尤其在挖齿-土体-管道等多方接触或冲击作用等领域。因此，本文基于SPH理论，建立挖齿-土体-管道间相互接触和冲击力学模型，分析土体和管道在挖齿作用下动态响应。研究结果可为挖齿-土体-管道等多方接触或冲击领域研究提供新的计算方法。

1 SPH理论基础及冲击方程离散

1.1 SPH理论基础

SPH法关键步骤为对任意宏观变量函数f(x)的积分近似[12]，以及对连续积分表达式的粒子近似。任意宏观变量函数f(x)的积分近似如式(1)所示：

(1)

(2)

连续积分表达式的粒子近似是将相关连续积分表达式转化为支持域内所有与粒子叠加求和的离散化方程，如式(3)所示：

(3)

(4)

在式(1)～(4)中，光滑函数Wij的选择对离散近似起重要作用。本文选取3次光滑核函数，如式(5)所示：

(5)

式中：R为点x和x′处2粒子间相对距离，R=r/h=|x-x′|/h，h为光滑长度，m；αd在1，2，3维空间下分别等于1/h，15/7πh2，3/2πh3。

1.2 基于弹性及弹塑性动力学理论的冲击方程SPH离散

基于SPH基本理论，对质量守恒方程、动量守恒方程离散得到SPH形式下连续性方程及动量方程，如式(6)～(7)所示：

(6)

(7)

式中：vi与vj分别表示粒子i与j速度，m/s；α，β代表坐标方向；σαβ为总应力张量，N；Fα表示外力引起加速度，如重力、齿轮齿条相互作用力及固壁边界作用力等;∏ij为人工黏度[13]，防止粒子接近发生非物理穿透，如式(8)所示：

(8)

(9)

(10)

(11)

(12)

vij=vi-vj

(13)

xij=xi-xj

(14)

式中：vij代表粒子i与粒子j速度差，m/s；xij代表粒子i与粒子j距离,m。

同时，引入XSPH速度平均方法，如式(15)所示：

(15)

式中：ε为修正因子，取0.3。

PE管道为具有黏弹性的聚合物，其弹性模量随时间增长而减少。由于挖齿与管道冲击速度极快，可忽略时间影响，将模型简化为分段弹性模型。对于土体模型，采用基于D-P屈服准则的弹塑性本构模型，如式(16)～(17)所示：

(16)

(17)

(18)

(19)

2 挖齿-土体-管道接触模型

定义挖齿、土体、管道为不同部件。其中，挖齿作为刚体不考虑其变形，土体假设为连续均匀弹塑性介质，管体因受力时间短而假设为分段弹性介质。在SPH方法中，粒子一般通过核估计计算相互作用，如果不额外定义特殊接触特性，当各部件相互分离时，接触面处粒子将产生与实际情况不符的剪切及拉伸力。因此，本文作以下3个假设：1)将土体、管体内部视为连续体，挖齿-土体-管道接触面处视为离散体。2)当接触面处不同部件粒子间距为各部件粒子平均间距1/2时，部件间粒子视为相互接触。3)当部件间粒子接触时，给部件间粒子施加1个具有弹簧特性的相互作用力，该力与各部件变形及弹性模量有关，绝对值随接触粒子间距离减小而增大，且施加给各部件作用力大小相等，方向相反。当部件间粒子间距大于部件粒子平均间距1/2时，表示部件未接触，即作用力为0；当部件接触时，该力在2维条件下如式(20)所示：

(20)

基于SPH采用FORTRAN编制埋地燃气管道挖掘破坏数值模拟程序。其中，密度采用连续密度法求解，粒子信息更新采用跳蛙法，计算流程如图1所示。

图1 基于SPH的埋地燃气管道挖掘破坏分析计算流程

3 数值算例

本文以挖齿挖掘薄层土覆盖管道为研究对象，管道外径0.104 m，壁厚0.008 m，土体长1 m，高0.24 m，挖齿截面为扇形，圆心角35°。挖齿视为刚体，并以1 m/s速度垂直向下运动；土体弹性模量20 MPa，泊松比0.25，黏聚力12.38 KPa，摩擦角20°；管体弹性模量1.115 GPa，泊松比0.3。模型共离散为1 986个粒子，其中土粒子1 262个，初始间距0.017 m；管粒子480个，初始间距0.005 m；挖齿粒子110个，初始间距0.01 m；边界粒子132个，初始间距0.012 5 m。时间步长4×10-7s，共计算8×105个时间步，即0.32 s。PE管道塑性阶段应力应变见表1。挖齿-土体-管道模型初始状态如图2所示。模型分别在0.208，0.288，0.32 s发生接触、挖裂和挖穿，模型发生接触、挖裂、挖穿时粒子分布如图3所示。

图2 挖齿-土体-管道模型粒子初始状态

表1 PE管塑性阶段应力应变

由图3可知，SPH法可呈现挖掘过程挖齿-土体-管道实时动态及破坏状态。在挖齿直接接触管道前，管道变形较小；当挖齿直接作用管道后，管道以接触点为中心产生较大变形直至断裂。

图3 模型发生接触、挖裂、挖穿时粒子分布

进一步采用椭圆度对管道变形定量评价，如式(21)所示：

(21)

式中：Dmax、Dmin分别为管道最大外径与最小外径，m。挖齿挖掘土体及管道全过程，管道椭圆度随挖掘时间变化关系如图4所示。

图4 管道椭圆度随挖掘时间变化关系

由图4可知，尽管在0.208 s前挖齿尚未直接接触管道，但部分挖齿作用力由土体传递给管道，仍会使管道产生一定椭圆度变化，但此时椭圆度不足0.1，不会对管道产生破坏；当挖齿与管道直接接触后，管道椭圆度迅速增大并且逐渐挖裂，挖裂后管道抵抗变形能力大幅降低，直至挖穿时椭圆度约为0.45。

各挖掘阶段管道环向应力及环向应变分布如图5～6所示。由图5～6可知，管道环向应力随挖齿向下作用而逐渐增大，管道由接触至挖穿过程中，最大应力应变发生在挖齿与管道接触部分，即管道正上方；当管道发生挖裂时，最大应力约35 MPa，等效应变约0.045；当管道挖穿时，环向应力约65 MPa，等效应变约0.08。此外，当管道贯通后，管道上方约45°位置会出现较大应力应变，易发生进一步破坏。

图5 各挖掘阶段管道环向应力分布

图6 各挖掘阶段管道环向应变分布

采用有限元法对相同模型模拟，以验证SPH法准确性。但有限元法挖裂后网格趋于畸变，因此，以挖裂初环向应力为例进行对比，如图7所示。有限元方法在挖裂阶段得到环向最大应力约35 MPa，与SPH法结果一致。

图7 挖裂阶段管道环向应力分布(有限元)

4 结论

1)当挖齿作用于土体而尚未作用于管道时，管道椭圆度受一定影响，但不会产生断裂或破坏；当挖齿直接作用于管道且继续向管道方向挖掘时，管道会发生局部断裂直至挖穿。

2)挖齿与管道直接接触点为管道发生破坏或断裂第1危险点；当管道继续大变形后，其左右两侧受拉力作用可能产生断裂破坏，为断裂第2危险点。

3)模型离散粒子数量及计算域等相关因素对计算精度产生一定影响。