基于离散-连续耦合法的埋地管道受荷过程模拟分析

侯超群,程一朋,李永鑫,孙志彬,蒋逍文

(合肥工业大学 汽车与交通工程学院,安徽 合肥 230009)

随着我国城市化不断推进,在城市道路下方布设管道以提升城市空间利用率的方式日趋多见[1]。由于管道多敷设于地下,属隐蔽工程,因此其变形、损坏不易发现,并具有危害性大等特点。近年来,埋地管道在荷载作用下的损坏现象频发,轻则影响管道的正常使用,重则导致严重的公共安全问题。埋地管道所承受的荷载作用主要通过管周土体进行传递,因此探讨埋地管道在荷载作用下管周土体承载和传力的机制,对保证地下管网的安全运营具有重要意义。

国内外诸多学者采用数值模拟方法对荷载作用下埋地管道的受力情况及管土相互作用开展了较为深入的研究[2-5]。数值模拟可采用的方法包括有限元法、有限差分法和离散元法等。有限元法和有限差分法基于连续介质力学，难以准确详细地描述细观力学行为；而离散元法虽然能够从细观尺度研究构筑物的力学行为，但当模型颗粒数目较多时，计算耗时较长，对计算机的算力也有较高的要求。连续-离散耦合的方法能够将离散元法和有限差分法的优势结合,发挥出良好的应用价值。文献[6-8]采用连续-离散耦合方法对静载下岩土结构的变形、颗粒位移场、土压力分布进行细观分析,并验证了耦合模拟的有效性;文献[9]采用三维离散-连续耦合法对预制桩沉桩过程进行模拟分析,研究了沉桩过程中桩周土颗粒的运动趋势和振动频率对助沉效果的影响;文献[10]采用连续-非连续介质耦合数值模型从碎石桩失效与细观颗粒受力运动的联系角度对碎石桩受荷的变形及破坏全过程进行了模拟分析;文献[11-12]基于土体连续-离散耦合算法,证明了土体连续-离散耦合方法可以从细观尺度描述重点研究区域的土体和结构的破坏过程以及盾构时隧道掌子面的细观颗粒特性。以上文献研究表明,采用连续-离散耦合算法可以从细观尺度研究重点区域或结构的力、位移的变化规律以及破坏机制过程,并能达到良好的效果。但现有文献中对埋地管道采用离散-连续界面耦合法的研究却鲜有报道。

以往对于荷载作用下埋地管道的研究中,多数学者采用连续介质数值模拟等方法来研究埋地管道的受力情况,难以从细观角度全面具体地揭示荷载作用下埋地管道的管周土体的动态作用过程。本文基于离散-连续界面耦合法,建立三维埋地管道数值模型,从宏、细观尺度结合对埋地管道及管周土体受荷过程中的受力情况和荷载传递机理进行分析研究。

1 埋地管道三维耦合数值模型建立

1.1 耦合原理

离散-连续界面耦合法,即对模型内部关键细观研究区域采用离散单元法以便于细观研究,在模型非关键研究区域采用连续介质单元法以提高计算效率。针对本文研究对象,对管周土体等核心细观研究区域采用离散单元,其余部分(原状土基及管道)采用有限差分法模拟,两者之间通过设置耦合界面元素实现速度与力的相互交换传输[13-14],具体如图1所示。

图1 耦合计算原理

由FLAC模型求解得到的速度通过耦合界面传递至PFC模型,由此PFC模型得到响应更新,其对耦合界面部位产生的作用力返回至FLAC模型,作为FLAC模型力学响应更新的条件,以此循环运算。本文中的离散-连续界面耦合法通过在耦合界面处设置墙单元Wall来实现PFC与FLAC之间的数据交换传输运算。

附着于FLAC模型元素表面(单元Zone表面或结构单元面)上的PFC模型组件Wall是由三角形拼接而成,如图2所示。C为某一球颗粒单元与三角形墙面(Wall单元)的接触点,PC为墙面上距C点最近的点。从PC到三角形的顶点采用重心插值法外推。

图2 耦合中力的传递

图2中：xi表示Wall角点及其附着部位FLAC3D的Zone节点的坐标;Ai定义为3个表征三角形的面积(在耦合变量解译过程中,1个Wall单元进一步分解为3个三角形)。定义ri为PC与xi的距离。施加在每一个网格节点的力为Fi(i=1,2,3),施加在接触点上的总接触力为F,接触点处因黏结产生的弯矩为Mb。由于接触部位的变形同时存在拉伸、剪切及其扭转作用,PC与C两点在空间的位置可能存在不共点的情形,因此作用在接触平面的总力矩M为:

M=Mb+(rC-rPC)×F

(1)

由此,FLAC3D模型元素在Wall附着部位即角点xi处,承受着来自于接触耦合作用PC处的集中力Fi和弯矩Mi,其值计算可采用加权法,计算公式为:

(2)

Mi=ri×Fi

(3)

当点PC外推触发三角形区域的节点和结构单元更新,使其网格点和节点能够直接添加刚度,即分别对FLAC3D和PFC3D设置单元相关参数,通过三角形区域介质实现力与速度信息的相互传递计算。

1.2 模型建立

本文基于PFC3D6.0程序软件平台,通过加载FLAC3D有限差分计算模块来实现埋地管道受荷模拟的构建。在建模过程中,考虑了管道埋置过程。其模拟过程如下:首先建立完整的连续网格模型;然后对模型进行开挖,并放置管道;通过在管道周围填充颗粒,完成管道埋置;最后,在地表施加上部荷载,进行加载模拟。

管道回填土采用PFC3D的ball颗粒模拟,管道及原状土基采用FLAC3D的实体单元模拟。在填充颗粒之前,在管道外侧及原状土地基表面通过命令流“Wall-Zone”生成耦合墙Wall单元,作为耦合计算数据交换的媒介。填土的区域为梯形,该区域需采用geometry命令进行定义,再采用ball distribute range geometry-space命令填充颗粒,并删除管道区域内的颗粒。颗粒生成完成后,颗粒与单元网格之间自动通过耦合墙Wall单元进行计算。

模型的顶面自由,不设约束条件,其余边界均施加法向约束条件。模型在重力作用下达到初始平衡后对位移清零,然后再进行管道埋置与荷载施加。荷载施加宽度与管道宽度相同,通过设置Wall单元以便对颗粒施加上部荷载。埋地管道管径选取D=1.5 m,选取空管工况,即管内无压力。

所建立的耦合模型如图3所示(单位为m)。

图3 埋地管道模型

1.3 模型宏、细观参数选取

颗粒细观参数见表1所列，模型宏观参数见表2所列。

表1 颗粒细观参数

表2 模型宏观参数

为保证耦合模型在宏观、细观上力学特性的一致性，颗粒细观参数是与三轴数据标定后获得的。

分析埋地管道力学性能时,考虑管周回填土与原状土基性质差异,对管周回填土与原状土基选取不同的变形模量和密实度。由于管道埋深和外部荷载因素(荷载施加于管道的正上方,作用宽度与管道宽度相同)对埋地管道力学性能和受荷作用过程的影响,各种数值模拟工况方案见表3所列。

表3 数值模拟计算工况方案

2 埋地管道受荷过程受力变形特征

通过建立离散-连续耦合模型,考虑不同工况下的数值模拟来分析管道及管周土体位移、管道径向变形及管周土压力的情况。

2.1 竖向位移分析

数值模型在不同荷载大小作用下(埋深H=3D)对应的管周土体及管道竖向位移分布如图4所示(单位为m)。

图4 不同荷载作用下对应的竖向位移分布图

从图4可以看出，管周土体竖向位移呈现“U”形分层分布,即中间位移大、两侧位移小。这是由于在竖向受压条件下管道收缩,管道顶部土体产生大幅位移量,同时管顶土体在向下方移动时产生的拖曳力对周边土体起牵拉作用,使管侧土体产生小幅位移量。随着施加荷载的增大,在地表产生大变形的范围也越大,同时大变形在地下的影响深度也越深,管顶产生的变形也愈大,对管道的影响程度也愈大。

2.2 管周土压力分析

外部荷载对埋地管道的影响经由管周土体传递,管周土体既作为传递荷载的媒介,自身重力又作为荷载以压力的形式直接作用在管道上,管周土压力的分布情况将对管道的损坏有直接影响,因此有必要研究在外荷载作用下管周土压力的分布规律。本文根据模拟结果提取得到管周应力,经换算和分解得到作用在管道上的竖直和水平方向的土压力,如图5所示。

图5 不同工况下管周土压力分布

从图5a可以看出,P1工况下的管顶与管底处管周土压力分布规律与P2或P3工况有差别,说明外部荷载的大小对土压力的规律有影响,造成管周土压力重分布。随着外部荷载的增大,管顶、管侧和管底处土压力峰值均增大。在P1工况下,管顶土压力变化趋势为先增大后减小,极值出现在管顶;但在P2和P3工况下,管顶土压力从管顶向管侧先缓慢增加再急剧减小,波谷在管顶处,波峰在管肩处附近。在P1、P2和P3工况下,管底土压力从管底部向管腋部呈现轻微波动,然后再急剧减少,在管底处受到地基反力作用和沟槽底部土体受力向内侧挤压使土压力产生波动。在P1、P2和P3工况下,水平土压力均呈现先增大后减小,在管侧处达到极大值,在顶部和底部达到最小。

由图5b可知,H1、H2和H3工况下的管周土压力分布规律与在P1、P2和P3工况下相似,不再一一叙述。可以看出管周土压力随着埋深的不断增加而增大,并且增大的幅度也越来越大;且随着埋深的增大,管顶土压力也会产生重分布。

综上可知,在管道上方总荷载(外部荷载与土体自重荷载之和)达到一定数值后,使管顶受压产生一定的压缩变形,管道上方土体产生“压力拱效应”,使土体荷载传递方式发生变化,从而使管周土压力重分布。

2.3 管道径向位移分析

埋地管道在荷载作用下发生变形,这种变形不仅会影响管周土压力的分布,还会造成管道的损毁,影响管道的正常使用。不同工况下管道径向位移变化如图6所示。

由图6可知,埋地管道在外部荷载作用下竖向发生挤压,水平向发生侧向鼓胀,管体截面在压力作用下由原先的圆形变化为椭圆状,其变形程度可以由椭圆度的计算公式来衡量。椭圆度的计算公式[15]为:

(4)

其中:Dmax为最大管径长;Dmin为最小管径长。

经由(4)式计算可得,在P1、P2、P3和H1、H2、H3工况下,管道变形的椭圆度分别为2.05%、3.18%、5.26%和1.42%、3.18%、6.15%。可知随着荷载的增大,管道椭圆变形的程度越大,并且管顶与管底竖向的相对位移差也增大,易造成管道压缩损坏。在实际工程中,应禁止车辆超载等外荷载过大而造成埋地管道损坏。随着埋深的增大,由于管道上部土体重力的增大,管道椭圆变形的程度也愈大。

3 埋地管道受荷动态过程分析

3.1 水平位移分析

在固定荷载作用下(埋深H=3D，荷载P2=0.3 MPa)，耦合模型不同时步的水平方向位移分布如图7所示。

从图7可以看出，管周土体位移和原状土基土的位移在耦合界面处取得了较好的连续性与一致性，表明离散元颗粒与有限差分法网格产生良好的耦合效果。

在荷载作用下，由于管道的强度高于土体的强度，使管道上方的土体沿管道中心对称面向两侧斜下方呈“八”字形放射状位移(见图7a)。随着计算时步的增加，沟槽内土体与原状土基产生连续的位移，并该“八”字形位移范围不断下移与扩大，逐步与管道相接触产生影响；在荷载的作用下，管道竖直向受压缩，产生水平向向两侧鼓胀变形，挤压致使管侧土体向斜下方位移运动，并与外部荷载产生的向斜下方的位移作用相结合，产生形似“蝴蝶”状的水平位移范围，位移范围逐渐向原状土基两侧下部扩散；最后水平位移大部分传至原状土基。

从图7中土体水平位移趋势可以看出，埋地管道在受荷载作用过程中，管道两侧回填土及原状土基能够承担一部分外部荷载。

图7 固定荷载作用下不同计算时步的模型水平位移图

3.2 细观接触力链分析

在颗粒自身重力或外部荷载的作用下,颗粒物质体系中的相邻颗粒间引发接触,接触产生的网络结构成为颗粒间传递作用力的途径,该传力路径称为接触力链[16]。

为了研究荷载作用过程中管周土颗粒的受力情况,选择不同时步的土颗粒间的接触力链进行分析。

通过耦合法构建的埋地管道模型,可以从细观尺度直接观察埋地管周土颗粒的受力情况。在固定荷载作用下(H=3D,P2=0.3 MPa),不同计算时步的管周土颗粒间接触力链的分布如图8所示(单位为N)。

接触力在图8中以柱体显示,圆柱体的颜色与接触力大小相匹配。由于土体细化成颗粒的数量众多，并且接触力梯度的间隔较大,为了便于观察分析，在图8中显示以500 N为尺度间隔的接触力。

图8 固定荷载作用下不同计算时步的管周土体接触力链图

从图8可以看出：地表土体由于直接承受外部荷载产生了少数量的强接触力链，且这些接触力链以一种树根状的形式向下方生长扩散并传递荷载；在外荷载向下的传递过程中，这些强力链逐渐向下延伸直至管顶，并向四周扩散形成数量更多强度稍低的弱力链，以使荷载向两侧土体分散；由强、弱力链相互联结形成的力链网络逐步向管顶、管侧及管底发展延伸，使管道产生变形。随着管顶受力作用产生竖直压缩变形，使其管顶小范围内土体松动(图8d中红虚线范围内，接触力链显示为弱力链)，从而管道两肩侧土体受力产生强力链，并形成"压力拱效应"，使管顶受力有所减弱，管肩处受力增大，这与竖向最大位移出现在管顶、竖向最大压力出现在管肩处的结论相符。

3.3 位移矢量分析

为研究荷载作用过程中管周土颗粒位移变化趋势,选择不同时步的土颗粒位移矢量场进行分析。管周土颗粒在荷载作用下位移矢量分布如图9所示,其中箭头为土颗粒的位移方向。

从图9可以看出,地表土体由于直接接触竖向荷载,在地表荷载作用区域内土体瞬时产生了向下的位移。由于原状土基的模量与密实度比管周填土的大,原状土基附近填土的位移发生隆起。随着计算时步的进行,荷载作用影响的范围越来越大,沟槽内管道两侧的土体位移方向由与竖直方向有一定的夹角慢慢趋于竖直方向,位移逐渐向原状土基扩散。当在荷载作用下,回填土的密实度达到一定值后,其沟槽内的土体整体均竖直向下移动。由于沟槽壁对回填土的下沉产生向上的摩阻力,使得接触沟槽壁的土体位移量会稍减小。最后,管道上方的土体绕过管道向斜下方运动,且具有明显的方向趋势。因为管道有一定的变形吸收作用,所以管底的土体位移量小于管顶的位移量,产生管顶与管底的相对位移差。

图9 固定荷载作用下不同计算时步的管周土体位移矢量图

4 结 论

本文基于离散-连续界面耦合法建立了三维PFC-FLAC耦合埋地管道数值模型,对埋地管道受荷载作用下管周土体及管道的受力变形全过程进行了模拟,得到主要结论如下:

(1) 在不同荷载作用下,体位“U”形分布,即中间大、两侧小。外部荷载的大小和管道埋深对土压力的分布规律有影响,产生“压力拱效应”,引起管周土压力重分布,使竖直土压力最大值出现在管肩处。

(2) 荷载作用过程中，管周土体的接触力链向两侧斜下方延展，这表明管道两侧回填土及原状土基在受荷过程中能够承担一部分土体自重和外部荷载。管周土体的接触力链分布形态与管顶土压力的分布规律相互印证。

(3) 通过模型水平向位移图、管周土体的接触力链网络和位移矢量图讨论了埋地管道的受荷作用动态全过程，并有效揭示了管道变形与细观土颗粒受力位移运动之间的联系以及管周土体的承载传力机制。