韩同春，邱子义，豆红强



基于颗粒离散元的抗滑桩土拱效应分析

韩同春1, 2，邱子义1, 2，豆红强3

(1. 浙江大学 滨海与城市岩土工程研究中心，浙江 杭州，310058；2. 浙江大学 软弱土与环境土工教育部重点实验室，浙江 杭州，310058；3. 福州大学 环境与资源学院，福建 福州，350116)

运用 YADE离散元软件对三维空间下的土拱效应进行研究。探讨在不同土层埋深、桩间距及土颗粒不均匀分布等影响因素下的土拱效应发展及变化。研究结果表明：浅层土中的土拱较之深层土中的土拱强度较差，破坏也有所提前；随着桩间距的增大，土拱效应的影响将会减弱，并且减弱的速度逐渐变慢；土颗粒粒径分布的不均匀性增强会使得土拱的形成提前，但是对加载过程中，桩上水平力分担比峰值的影响不大。

土拱效应；YADE离散元；三维模型；土层深度；桩间距变化；土体不均匀性

土拱效应是岩土工程中一种常见的现象，TERZAGHI[1]在1943年就通过“活动门”试验证实了土拱效应的存在。近些年来，众多学者分别从数值模拟、理论推导和室内试验等角度对土拱效应进行了深入研究[2−7]。李忠诚等[8]利用三维有限元结合ITOS 等[9]提出的土压力理论对成拱效应中桩间距的影响进行了系统的分析；杨明等[10]使用离心机对桩间的土拱效应进行了研究；向先超等[11]利用二维的离散元也进行了部分的研究。上述研究主要是建立在连续介质的基础之上，事实上土体是一种非均质各向异性的散粒体，利用有限元连续介质模型的研究必然存在一定的误差。并且，除此之外，由土拱效应的触发机理可知，土拱效应的产生过程即为桩后土体及桩间土和桩后土的水平位移发展不协调及运动土体的“楔紧”过程，土体的压密发生在三维空间中，仅仅在二维平面内考虑土拱问题与实际不符。随着离散元在岩土工程中的广泛应用[12−14]，该方法对于大变形和非连续介质模拟的适用性良好的优点越来越突出。为此，本文作者采用离散元软件对三维空间中抗滑桩土拱效应进行研究，探讨不同深度土层的成拱形式及桩间距、土颗粒粒径分布等因素对土拱效应的影响。

1 土拱效应形成机理简述

土拱的形成主要是土体具有颗粒间的黏结力和摩阻力[2]。就抗滑桩而言，土拱效应产生的机理为桩后土体在土压力或者滑坡的作用下产生位移，而抗滑桩的存在致使桩后土位移的发展受阻。与此同时，土体为典型的非均质各向异性散粒体材料，表现为桩后土体位移的发展在土压力作用下并不均匀，进而使得土体间发生“楔紧”作用。在土体自身存在抗剪强度及桩侧的摩阻力作用下，桩后土体形成以相邻两桩为拱脚的土拱，且桩后的土压力随着土拱传递到拱脚及其周围土体上(见图1，其中1为桩径，2为桩间距)。所以，土拱效应的实质是应力在土体和抗滑桩之间的转移过程，即从运动的土体转移到抗滑桩上，使桩土共同分担荷载。

图1 桩间土成拱分析示意图

2 YADE离散法基本理论简介

离散元法的主要思想是将整个研究的物体离散为一系列刚性小球，小球几何特性(直径、空间位置、排列等)由人为控制产生，并且可以赋予刚性小球一定的物理特性(重度、摩擦角等)或是化学特性。小球运动受经典力学理论控制，整个研究物体的宏观物理性质由小球的运动和相互作用描述。

就YADE离散元软件而言，它是一个内核基于C++语言编写，外部的控制代码为Python语言编写的开源代码非商业三维离散元软件，方便用户基于此软件进行二次开发和更加深入的研究，亦可以根据自己的需要调用软件中相应的模块来实现对工程问题的模拟。而且，在离散元法的计算效率上，接触面的发现算法是影响计算效率的重要因素。YADE离散元程序在接触面的发现算法中采用“边界盒法”，具有较高的效率[15]。

参考该离散元软件的教程[16],两刚性小球的接触刚度分为法向刚度N和切向刚度T。切向刚度T由法相刚度N乘上T/N得到。其中，T/N决定宏观上泊松比的取值，对于法向刚度N，用2个长度等于球体半径的弹簧链表示，结合图2得

式中：1和2分别为两接触小球的弹性模量；和在本次模拟采用的接触特性中取=21，=22。

图2 两小球法向接触刚度

对于应变的计算，YADE中将应变分为法向应变和切向应变。法向应变的计算如图3所示。其中，和分别为两小球初始球心坐标，定义：

图3 两小球接触初始状态图

法向位移N及法向应变N的计算式如下：

对于压缩大变形的情况，由式(5)和(6)可知：其应变将趋向于负无穷，但是小球间的力还保持为一个有限的值，会发生球体间彼此穿透的现象。所以，对法向应变的计算式进行如下调整：

剪切应变的计算有全量算法和增量算法2种，增量算法的剪切位移为

对于应力的计算，在已知的法向和切向位移条件下，正应力N和切应力的计算式如下：

式中：为土体的内摩擦角。

3 数值模型

参照杨明等[10−11]建立的二维离散元模型,建立三维离散元模型。根据对称性，取排桩中相邻两桩中心线土体为研究对象，所研究的土体限于桩后及两桩中间的土体。由于所研究的土体为某一深度处的土层，在上部土层的约束作用下，土体竖向的位移相对于水平向位移较小，并且，整个土层的竖向位移为零并不会导致土颗粒间的“楔紧”作用仅发生在水平方向，土拱效应的产生依然是一个三维的过程。所以，作近似处理，假定土体的位移仅仅发生在水平方向上，在矩形盒的顶部和底部施加竖直向的刚性约束。

包围土体的矩形盒表面光滑，前端开口。为保证土体因水平方向上产生的不均匀位移引发的“楔紧”作用不是人为施加，而是抗滑桩和土颗粒间及土颗粒和土颗粒间相互作用产生的结果，使用竖向加载墙对各土层施加以相同的初始水平位移。加载的速度为0.05 m/s，所施加的速度方向为沿矩形盒长度方向。矩形盒的长度为5 m，高度为3 m，宽度依据不同的桩间距调整。

抗滑桩为直径0.9 m的圆截面桩，桩土摩擦因数为0.5，忽略桩的弹性变形，视之为刚体且桩端完全固定，无水平方向的位移。同时忽略桩后土体的土压力作用以形成临空面。

采用刚性小球来模拟离散土体，若是按照实际的尺寸来模拟土颗粒，则所建立的模型中小球将多达数十万个，以目前计算机的计算速度无法满足要求。而且，本文所研究的是土拱效应在不同条件下的变化趋势，数值仅需比较相对大小，并不需与实际的工程做比较。所以，模拟中对土颗粒的尺寸进行适当的放大，土体的孔隙率为15%，小球直径为0.02 m。根据室内三轴试验和YADE离散元模拟实验之间的参数对比分析，确定土体的其他细观物理力学参数见表1。

表1 土体细观参数

依据以上边界条件及土体的各项物理指标，利用YADE开源代码离散元软件，在程序中定义符合上述物理指标的材料直接赋予所产生的刚性小球，并且，控制小球产生的数量来满足初始孔隙率的要求。建立数值模型如图4所示。

图4 离散元三维数值模型

4 土拱效应影响因素的数值分析 结果

4.1 不同深度土层的成拱形式

取抗滑桩的净间距为2.7 m，程序中虚拟时间每过1 s保存1次加载过程中各离散元小球的物理状态。保存的数据利用paraview后处理软件进行处理，即可得到离散元小球接触面上正应力力传递图。通过观察离散元小球接触面上正应力力传递图的变化，来间接得知抗滑桩间土拱的形成、发展和破坏过程。所进行的模拟步为25 000步左右，所用时间约为1.5 h。

对于整个模型所模拟的土体，在加载墙不同的推进位移情况下，如图4所建立坐标系，取不同处的正应力传递图的水平剖面，可得结果如图5所示的土颗粒小球接触面间正应力传递图(为加载墙位移)。从图5可以看出：处于不同埋置深度的土层，所成的土拱形状也不相同。随着土层深度的增加，土拱的矢高逐渐降低。对于同一处的土层水平剖面，当= 0.5 m时，随着加载墙上加载位移的发展，土拱在早期就可形成，并且在加载位移=5~15 cm时，土拱的形状变得更为清晰，在=40 cm时，两相邻抗滑桩间的应力分布图已出现不连续的现象，土拱发生破坏；当=1.5 m时，土拱在=5 cm时，土拱已然形成，在= 15 cm时，土拱形状已经很完整，在=15~40 cm的过程中，土拱已经开始产生破坏，但是破坏的程度不如=0.5 m时明显；当=3.0 m时，在整个加载过程中，土拱形状始终完整。