基于非圆颗粒Trapdoor试验模拟土拱效应

赵志国，王炳龙，韩高孝，唐 斌，梁 杰

（1.同济大学道路与交通工程教育部重点实验室，上海201804；2.中水电成都建设投资有限公司，成都610212；3.中国水利水电第七工程局有限公司，成都610081）

土拱现象在岩土工程实践中较为常见［1-2］，例如桩基工程中针对土拱效应的研究，对解决桩基合理间距［3-6］、桩基承载力以及应力传递特性［7-9］等问题具有重要意义。另外，如何合理确定砂土地层中隧道土拱的范围及土拱效应的稳定性是决定隧道土压力荷载是否合理的关键因素［10］。实践表明在岩土工程中充分利用拱效应，不仅有助于预防工程事故的发生，而且可以得到良好的工程效益。太沙基［11］利用砂土的Trap-door 试验证实了土拱效应的存在，在试验的基础上推导了松动土压力的解析表达式；武军［10］基于椭球理论、陈若曦［12］考虑土拱作用下主应力偏转、李文博［13］基于侧压力系数分别对太沙基松动土压力公式进行了修正；Adachi［14］等利用三维砂土Trap-door试验模拟隧道施工过程支护压力的变化；Fin［15］运用连续介质理论给出了Trap-door 平移或旋转时竖向应力的平面应变弹性解；Evans C H［16］给出了无黏性、均匀、各向同性土体中Trap-door 向上和向下运动时的平面应变塑性解；Koutsabelouis［17］利用有限元分析不同高宽比条件下Trap-door下降过程中Trap-door上竖向应力的变化。

虽然诸多学者针对土拱效应问题提出了不同的计算模型，但实际上土颗粒并不是连续的，应用有限元、有限差分等方法不能有效地模拟土颗粒的大变形、流动等诸多问题。离散元方法克服了传统连续介质力学模型的宏观连续性假设，可从细观层面上对土的工程特性进行数值模拟，并通过细观参数的研究来分析宏观力学行为。

韩高孝［18］利用PFC2D中二维圆盘单元模拟砂土中的土拱效应，而文献［19］指出黏性土颗粒可采用一系列的二维圆盘单元来模拟，直接采用二维圆盘单元来模拟砂土颗粒是不适当的。基于以上原因，本文从散粒体和细观角度，应用非圆颗粒对砂土Trap-door试验过程中的土拱效应进行分析，以期更真实地反应砂土中土拱效应的形成机理。

2 非圆颗粒Trapdoor试验DEM模拟方法

2.1 非圆颗粒模拟

PFC2D中提供了模拟非圆颗粒的CLUMP方法用于非圆颗粒的二次开发，所谓CLUMP就是通过纯圆颗粒的粘结组成任意形状颗粒，颗粒之间可以重叠也可以不重叠，CLUMP单元在受力过程中不会被破坏，从而提高了非圆颗粒的计算效率。

本文即利用“团颗粒”方法，开发形状近似椭圆的“椭圆团”颗粒模拟砂土。模拟椭圆颗粒时，首先生成纯圆颗粒，然后，按照质量等效和面积等效原则进行转化［20］，具体思路如图1所示。

图1 初始纯圆颗粒与“椭圆团”颗粒转化示意图Fig.1 The conversion chart of the initial round particles to the elliptical particle

2.1.1 面积等效原则

假设纯圆颗粒半径为R，“椭圆团”颗粒由主圆颗粒和两个副圆颗粒组成，主圆颗粒半径为rb，副圆颗粒半径为ra，“椭圆团”颗粒长轴长为A，短轴长为B，颗粒长短轴比为Se，则满足

初始纯圆颗粒与开发的“椭圆团”颗粒面积相等，则满足

式中：S纯圆为初始纯圆颗粒面积；S主圆为“椭圆团”颗粒主圆颗粒面积；S副圆为“椭圆团”颗粒副圆颗粒面积；S叠合为“椭圆团”颗粒主圆颗粒与副圆颗粒叠合面积。

2.1.2 质量等效原则

“椭圆团”颗粒的质量是组成团颗粒的各个圆盘颗粒的质量之和，则满足

式中：ρ纯圆为初始纯圆颗粒密度；ρ′为“椭圆团”颗粒密度。

2.2 Trapdoor试验DEM模拟

2.2.1 数值试样制备

整个数值模型中选用三面墙体模拟模型箱，模型箱大小为0.4 m×0.4 m，其中底板中心处可以移动的墙体用以模拟Trap-door，墙体刚度的取值十分关键，过大边界效应明显，过小则造成颗粒逃逸，本文取颗粒刚度的10倍，颗粒细观参数取值参照文献［18］的结果。

颗粒利用落雨法制备，即：①在0.4 m×0.4 m 的模型箱中生成纯圆颗粒，初始孔隙比取0.45；②编写FISH程序找到每个纯圆颗粒的中心位置及其颗粒尺寸，删除纯圆颗粒，利用FISH二次开发程序按照质量等效原则和体积等效原则生成“椭圆团”颗粒；③每隔一定运算步消除颗粒中的位移场，在重力作用下颗粒下沉，删除指定高度以上的颗粒；④反复循环，直至颗粒之间平均不平衡力与平均接触力之比小于1%即认为颗粒制备完成，颗粒达到平衡状态；等待试验，最终生成Trap-door模型如图2和图3所示。

2.2.2 Trapdoor试验DEM模拟方法

图2 Trap-door DEM模型Fig.2 The model of the Trap-door DEM

图3 Trap-door DEM局部放大图Fig.3 An enlarged view of the Trap-door DEM

为了模拟非圆颗粒土拱效应，按照控制单一变量的原则，本文选择了以下11种方案进行模拟，其中方案1、2、3和4保证颗粒形状、颗粒摩擦系数和Trap-door宽度不变的前提下，控制Trap-door位移，用于研究Trap-door位移对土拱效应的影响；方案5和6保证颗粒摩擦系数、Trap-door宽度和Trap-door位移不变的前提下，控制颗粒形状变化，研究颗粒长短轴比即颗粒形状对土拱效应的影响；方案7、8和9 保证颗粒形状、Trap-door宽度和Trap-door位移不变的前提下，改变颗粒摩擦系数，用于研究颗粒摩擦特性对土拱效应的影响；方案10和11保证颗粒形状、颗粒摩擦系数和Trap-door位移不变的前提下，改变Trap-door宽度用于模拟Trap-door大小对土拱效应的影响，主要计算参数和方案如表1所示。

表1 计算参数及其方案表Tab.1 Calculation parameters and program

整个试验过程中通过移动模型箱底板上的6 号墙模拟Trap-door，墙体移动速度为0.005 m·s-1，试验过程可以认为是准静态的，即土体内部颗粒有足够的时间进行应力和位移调整。

3 结果分析

为了从宏观上描述Trap-door试验过程中的土拱效应，通过编写相应的FISH 程序记录模型中横向和纵向每隔4 cm共48个监测点不同时刻的水平位移和竖向位移，同时在模型箱底板上布置直径为8 cm的5个测量圆监测试验过程中底板应力变化，监测点布置情况如图4所示。

图4 Trap-door试验监测点布置示意图（单位：cm）Fig.4 Arrangement of monitoring points

3.1 Trap-door位移对土拱效应的影响分析

以颗粒长短轴比为1.4、颗粒摩擦系数为1和Trap-door宽度为50 mm为例分析在不同Trap-door位移条件下“椭圆团”颗粒土拱效应。

从图5可知，由于Trap-door 的下移，造成应力释放，模型中不同埋深处的颗粒相应地产生沉降；随着Trap-door的下移，距离Trap-door越远的截面沉降越小，沉降曲线呈两头小、中间大的特点；在模型箱底面距离模型箱中心8 cm外沉降很小，主要影响区域宽度为16 cm，这一宽度约为3.2倍的Trap-door宽度。由图6可知，当Trap-door位移为8.15 mm时F截面处的位移几乎相等，即出现了所谓的等沉现象，这也是土拱效应的宏观表现。

图5 不同Trap-door位移条件下A截面位移曲线图Fig.5 Displacement curve of A cross-section in different displacement of Trap-door

图6 不同Trap-door位移条件下F截面位移曲线图Fig.6 Displacement curve of F cross-section in different displacement of Trap-door

由图7可知，当Trap-door 位移达17 mm 时，模型中颗粒最大位移为8.472 mm，小于Trap-door 位移；更为有趣的是，截面D以上的所有颗粒沉降几乎相等，形成所谓的等沉面，这便为“椭圆团”颗粒中土拱效应的宏观表现，而这一截面高度即为土拱效应影响的上限，上限高度约为16 cm，这一高度约为3.2倍的Trapdoor宽度。

通过布设在模型箱底板的测量圆获得Trap-door 在移动过程中底板应力变化值，由图8可知，随着Trap-door 的下移，Trap-door 位移从0.849 mm 增大到17 mm 的过程中Trap-door 上的应力由202.4 N 变为952.07 N；并且，随着Trap-door的下移，距离模型中央约5 cm范围内底板应力大幅度减小，距离模型中央两侧各4～15 cm，模型箱底板应力有所增加，这一宽度约为2.2倍Trap-door 宽度；由此说明模型箱底板应力在不断调整，这个过程正是土拱效应形成的过程。这一现象与文献［21］运用砂土通过室内试验获得规律是一致的，这也从侧面验证了本文数值试验模型的正确性。

图7 Trap-door位移为17 mm条件下各截面位移变化曲线图Fig.7 Displacement curve of different cross-section under the condition of Trap-door displacement equal to 17 mm

图8 不同Trap-door位移条件下底板应力变化曲线图Fig.8 Stress curve of floor in different displacement of Trap-door

3.2 颗粒长短轴比对土拱效应的影响分析

取Trap-door位移为17 mm时A截面颗粒竖向位移变化分析，由图9可知，随着颗粒长短轴比的增大，A截面颗粒竖向位移由10.66 mm减小到8.472 mm，由图10可知，底板中心应力变化值随着颗粒长短轴比的增大而增大，由718.8 N增大到1 008.4 N，土拱效应随着颗粒长短轴比的增大而更加明显。

图9 不同长短轴比A截面颗粒位移变化曲线图Fig.9 Displacement curve of A cross-section in different ratio of axis of particle

图10 不同长短轴比底板应力变化曲线图Fig.10 Stress curve of floor in different ratio of axis of particle

3.3 颗粒摩擦系数对土拱效应的影响分析

取Trap-door位移为17 mm时，A截面颗粒竖向位移变化分析，由图11可知，随着颗粒摩擦系数的增大，A截面颗粒竖向位移由12.36 mm减小到3.03 mm，由图12可知，底板中心应力变化值随着颗粒摩擦系数的增大而增大，由1 008.41增大至1 258.95 N，综合可知，土拱效应随着“椭圆团”颗粒间摩擦系数的增大而更加明显。

图11 不同摩擦系数A截面颗粒位移变化曲线图Fig.11 Displacement curve of A cross-section in different frictional coefficient

图12 不同摩擦系数底板应力变化曲线图Fig.12 Stress curve of floor in different frictional coefficient

3.4 Trap-door宽度对土拱效应的影响

取Trap-door 位移为17 mm 时，A截面颗粒竖向位移变化分析，由图13 可知，随着Trap-door 宽度的增大，A截面颗粒竖向位移由5.07 mm增大到8.47 mm；由图14可知，底板中心应力变化值随着Trap-door宽度的增大而减小，由1 163 N减小至1 008.4 N，综合可知，土拱效应随着Trap-door宽度的增大而减弱，这主要是由于应力释放的增大所造成的。

图13 不同Trap-door宽度A截面颗粒位移变化曲线图Fig.13 Displacement curve of A cross-section in different width of Trapdoor

图14 不同Trap-door宽度底板应力变化曲线图Fig.14 Stress curve of floor in different width of Trapdoor

4 结论

本文通过应用PFC2D中CLUMP方法二次开发的“椭圆团”颗粒模拟砂土在Trapdoor试验过程中的土拱效应，并分析了Trapdoor位移、颗粒形状、颗粒摩擦系数和Trapdoor宽度对土拱效应的影响，最后得到如下结论：

1）颗粒长短轴比为1.4、颗粒摩擦系数为1、Trap-door宽度为50 mm时，“椭圆团”颗粒中土拱影响宽度约为3.2倍的Trap-door宽度，高度约为3.2倍的Trap-door宽度；

2）由于土拱效应，Trap-door位移从0.849 mm增大到17 mm的过程中，模型箱底板应力在不断调整；

3）增大颗粒长短轴比、颗粒摩擦系数可以增强土拱效应；

4）减小Trap-door宽度即减小应力释放可以增强土拱效应。

[1] 贾海莉,王成华,李江洪.关于土拱效应的几个问题[J].西南交通大学学报,2003,38(4):398-402.

[2] 杨雪强,何世秀,庄心善.土木工程中的土拱效应[J].湖北工学院学报,1994,9(1):1-7.

[3] 杨明,姚令侃,王广军.桩间土拱效应离心模型试验及数值模拟研究[J].岩土力学,2008,29(3):817-822.

[4] 蒋良潍,黄润秋,蒋忠信.黏性土桩间土拱效应计算与桩间距分析[J].岩土力学,2006,27(3):446-450

[5] 赵明华,廖彬彬,刘思思.基于拱效应的边坡抗滑桩桩间距计算[J].岩土力学,2010,31(4):1211-1216.

[6] 冯君,吕和林,王成华.普氏理论在确定抗滑桩间距中的应用[J].中国铁道科学,2003,24(6):79-81.

[7] BRITTON E J,NAUGHTON P J.The arching phenomena observed in experimental trap door model tests[J].Geo-Frontiers, 2011:788-797

[8] 王漾,周萌,宫全美,等.高速铁路CFG桩网复合结构设计参数分析[J].华东交通大学学报,2010,26(6):17-22.

[9] 元翔,宫全美,王炳龙.高速铁路CFG桩复合地基处理方案数值模拟[J].华东交通大学学报,2008,25(1):32-36.

[10] 武军,廖少明.基于颗粒流椭球体理论的隧道极限松动区与松动土压力[J].岩土工程学报,2013,35(4):714-721.

[11] TERZAGHI K.Theoretical soil mechanics[M].New York:Wiley,1943.

[12] 陈若曦,朱斌,陈云敏,等.基于主应力轴旋转理论的修正Terzaghi松动土压力[J].岩土力学,2010,31(5):1402-1406.

[13] 李文博,陶连金,蔡东明,等.地铁隧道竖向土压力计算公式探讨与改进[J].铁道建筑,2013(3):78-81.

[14] ADACHI T, KIMURA M, KISHIDA K.Experiment study on the distribution of earth pressure and surface settlement through three-dimensional Trapdoor tests[J].Tunnling and underground space technology,2003,18:171-183.

[15] FINN.Boundary value problems of mechanics[J].Journal of the Soil Mechanics and Foundation Division,ASCE, 1963,89(SM5):39-72.

[16] EVANS C H.An examination of arching in granular soils[D].Cambridge:Massachusetts Institute of Technology, 1983.

[17] Koutsabelouis N C,Griffiths D V.Numerical modelling of Trap door problem[J].Geotechnique,1989,39(1):77-89.

[18] 韩高孝,宫全美,周顺华.摩擦型岩土材料土拱效应微观机制颗粒流模拟分析[J].岩土力学,2013,34(6):1791-1798.

[19] 胡欣雨,张子新.不同地层条件泥水盾构开挖面失稳状态颗粒流模拟方法研究[J].岩土力学,2013,32(11):2258-2267.

[20] 史旦达.单调与循环加荷条件下砂土力学性质细观模拟[D].上海:同济大学,2007.

[21] 加瑞.盾构隧道垂直土压力松动效应研究[D].南京:河海大学,2007.