0.   引言

近年来，随着中国城市建设和相关基础设施的快速发展，山区公路或铁路修建及运行过程中滑坡问题日益突出。抗滑桩以其抗滑能力强和便于施工等特点，广泛应用于滑坡治理或边坡防护工程中。桩间水平土拱是抗滑桩加固边坡工程中一种常见现象，表现为土颗粒间产生互相“楔紧”作用，实质上是土体内部的一种应力转移现象^[1]。土拱效应在抗滑桩发挥其支护功能中起到重要作用，查明土拱作用的规律对边坡防护的抗滑桩工程具有重要指导意义。

土拱形成过程比较复杂，有关土拱效应作用机理及演变规律的研究方兴未艾。Low等^[2]通过模型试验和理论分析，研究了桩基软土地基中路堤土拱效应。芮瑞等^[3]探讨了砂填料桩承式路堤拱效应传力机制，揭示了填料内部存在的初始三角形松动滑移面。钟卫等^[4]通过引入能够表征土体相对变形和不均匀变形程度的无量纲参数，对土拱动态演化过程进行分析，将土拱发育分为弹性形成、塑性发展和最终破坏3个阶段。蒋明镜等^[5]利用离散元模拟活动门试验对净砂与胶结砂土中土拱进行研究，结果表明两者破坏形态相似。尽管土拱效应的相关研究日益深入，但目前研究主要侧重于自重作用下竖向土拱效应，对于侧向推力作用下水平土拱效应研究仍较缺乏，亟待改进。

为了深入探究水平土拱成拱机理，吕庆等^[6]采用平面应变有限元分析模型，通过数值计算结果中最大主压应力迹线，拟合了土拱轴线的描述方程，并表明桩间距对土拱形状的影响最为显著。Bosscher等^[7]对砂土边坡土拱效应进行了室内试验，结果表明，随着桩间距增大，荷载分配给桩的比例逐渐减小，即土拱效应荷载传递效率逐渐减小。进一步地，向先超等^[8]基于二维离散元方法研究了部分宏-微观参数（抗滑桩桩间距、截面尺寸、颗粒粒度组成）对土拱效应形成过程的影响，并揭示了这些因素对土拱效应性能影响的变化规律。

在土拱效应理论研究方面，尽管有研究探讨了黏性土和无黏性土在桩间土拱效应方面的影响^[9]，但现有研究多数尚未清晰区分密砂和松砂条件下桩间土拱效应。实际上，这两者的土拱形成机理存在差异。与此同时，目前大部分关土拱效应的结论皆基于联合拱模型得出，暂未见专门针对单独拱模型的研究^[10]。

鉴于此，本文采用离散元方法，从密砂与松砂的角度出发，采用颗粒流软件建立了同时包含摩擦拱和端承拱的计算模型，深入探究砂土质边坡抗滑桩加固时桩间水平土拱机理及其演变规律，并讨论了宏-微观参数对土拱形成过程与土拱效应性能的影响。

1.   离散元模型

近年来，离散元方法（DEM）已在各类岩土工程领域的数值分析中得到了广泛尝试和应用，尤其侧重于模拟颗粒材料的力学行为。鉴于本文研究对象为净砂，不考虑粒间胶结作用，应选取无胶结接触模型，可采用简单接触模型和完整接触模型^[11]。参考文献[12]的研究方法，本文采用二维颗粒流软件PFC^2D的线性接触模型（linear model）模拟密砂与松砂中的土拱现象。

1.1   成样方法：分层压缩法

在DEM模型中，尤其在包含大量颗粒的情况下，生成的试样须保持一定的均匀性，即相邻两个颗粒几乎不重叠且没有较大空隙分隔。为了避免不均匀性和各向异性，本文采用分层压缩法（multi-layer method）^[13]生成试样。分层压缩法的主要实现过程如下：①扩大竖直方向颗粒生成区域并生成颗粒；②向下移动顶部墙体单元压缩颗粒并使颗粒移动到指定的第一层区域；③执行足够的循环时步以实现在指定区域内颗粒平衡并生成第一层均匀试样；④将压缩过程划分为指定的层数，重复上述步骤，以达到指定的孔隙率和试样尺寸。如图 1所示，简要说明了利用分层压缩法生成一个4层初始试样的流程图。若想获取更为均匀的试样，可采用考虑欠压准则的分层欠压法（multi-layer with undercompaction method）^[13]。

图  1  分层压缩法成样步骤

Figure  1.  Procedure of multi-layer method for generating specimen

下载: 全尺寸图片 幻灯片

最后，为了使模拟过程接近真实的物理状态，采用伺服机制进一步模拟地下一定深度的水平自重应力。所谓伺服，就是通过模型边界的调整，使得颗粒体系间的接触达到所需要的状态，然后再在其基础上开展加载分析。在利用分层压缩法生成初始试样后，通过循环计算命令使得颗粒平衡（图 1所示的步骤⑥），然后在初始试样运算至平衡状态的基础上开展伺服并在完成后进行位移和速度清零，最后获得所需要的试样。

1.2   微观与宏观参数的标定

在DEM中，特定试样的微观参数标定是通过数值双轴试验和Mohr-Coulomb破坏准则来确定的。简而言之，Mohr-Coulomb破坏准则中的宏观抗剪强度参数是通过测量数值双轴试验在不同恒定围压情况下的抗压强度间接获得的。数值双轴试验的主要实现过程如下：①通过分层压缩法生成初始试样；②利用伺服机制以保证试样围压为一恒定值；③通过顶部与底部墙体单元压缩试样施加轴向应力直至达到目标轴向应变（本文采用15%作为最大轴向应变）；④记录相应的围压和轴向应力，然后改变围压，并重复上述步骤。

对于DEM模型，选取的微观参数是基于Li等^[12]的研究，具体参数如表 1所示。考虑到滚动阻力和实际砂粒呈棱角状，试样与拱脚的摩擦系数均取较大值以模拟实际情况，而左右边界墙的摩擦系数设置为0是为了消除边界效应。鉴于刚度比对试样强度特性的影响较为有限，刚度比取值1.0^[14]。通过在颗粒生成阶段分别赋予摩擦系数0和实际摩擦系数来改变生成试样的孔隙度，进而分别生成密砂和松砂试样。对于密砂和松砂，均在数值双轴试验中施加了3个围压（50，75，100 kPa）来标定参数。图 2展示了在50 kPa围压条件下，密砂和松砂受剪时的应力-体变-轴向应变全过程。从图 2中可以看出，随着轴向应变的增加，松砂的强度逐渐增大，而密砂的强度达到一定值后，随着轴向应变的继续增加强度反而减小，曲线呈现一个明显的峰值。密砂和松砂的轴向应力最终都趋于某一相同的恒定值，这一恒定强度通常称为残余强度。同时，松砂受剪时体积缩小，而密砂一开始表现为缩小，后面持续膨胀。砂土剪切性状表明，松砂表现为典型的应变硬化和剪缩特性，而密砂表现为典型的应变软化和剪胀特性。75和100 kPa围压的情况与50 kPa相似，故不再赘述。如图 2所示，通过Mohr-Coulomb破坏准则分别得到了密砂和松砂的内摩擦角是36°和22°。由于线性接触模型不含有抗转动能力，故得到的摩擦角相比实际砂土较低。

表  1  DEM模型的微观参数

Table  1.  Mesoscopic parameters of DEM model

参数	法向刚度/ (N·m-1)	切向刚度/ (N·m-1)	摩擦系数	颗粒密度/ (kg·m-3)	最大粒径/ mm	最小粒径/ mm	平均粒径/ mm	孔隙度
砂土(ball)	1.0×107	1.0×107	1.0	2500	1.2	0.5	0.78	0.25
边界(wall)	1.0×1012	1.0×1012	0	—	—	—	—	—
拱脚(wall)	1.0×109	1.0×109	5.0	—	—	—	—	—

下载: 导出CSV 

| 显示表格

图  2  50 kPa围压下的砂土剪切性状及参数标定

Figure  2.  Shear behaviours of sands under confining stress of 50 kPa and parameter calibration

下载: 全尺寸图片 幻灯片

1.3   计算模型

基于文献[15]的研究方法，将土拱问题简化为二维平面应变问题。鉴于PFC^2D软件中圆盘颗粒为单位厚度，水平土拱模型选取地表下一定深度的单位厚度土层作为研究对象，并假定：①该单位厚度土层的位移限定在滑动方向上；②拱脚水平位移为零，即忽略拱脚的弹性侧向变形，将拱脚视为完全刚性。

如图 3所示，建立抗滑桩加固边坡中的水平土拱简化二维土拱模型示意图。考虑模型的对称性，两个相邻的方桩（拱脚）和中心范围内的土体作为数值模型的分析对象。基于Lai等^[16]的研究以及考虑后续产生相对位移的影响，模型的高度设置为0.5(s-a)+ 1.7a+3a以保证形成完整的土拱结构。基于Tang等^[17]的室内模型试验，本文选取抗滑桩宽度a为5 cm。同时，依据尹小涛等^[18]的研究，本文选取的最小颗粒粒径与试样模型长度比小于0.01，能够保证DEM数值试验具有较好的稳定性。计算模型中不考虑桩前土体抗力作用，去除了模型下部边界的位移约束以形成自由面。为了保证形成的土拱结构具有稳定性，最大加载位移Δ设置为2a^[16]。同时，为了研究土拱形成过程中的颗粒分布，在利用分层压缩法生成试样的同时对每一层颗粒进行分组染色，效果如图 4所示。

图  3  三维土拱效应问题的平面应变模型简化示意图

Figure  3.  Simplified schematic diagram of soil arching from three-dimensional to plane strain model

下载: 全尺寸图片 幻灯片

图  4  简化二维土拱模型

Figure  4.  Simplified two-dimensional soil arching model

下载: 全尺寸图片 幻灯片

由于联合拱的分布区域、极限承载力和塑性区分布均可近似为两单独拱情况的叠加^[10]，将端承拱和摩擦拱分别单独考虑后叠加。如图 4所示，土拱模型能够同时考虑地下一定深度处受围压作用的端承拱和摩擦拱。本文采用拱跨与拱脚宽度的比值（拱跨比）进行后续参数研究，其定义为

式中：a为桩宽；s为桩间距。常用于评价土拱效应性能的指标是桩土应力比和荷载传递效率。鉴于本文中抗滑桩尺寸是固定的，消除了桩截面尺寸的影响，桩土应力比与荷载传递效率在表示土拱效应作用发挥程度方面是一致的，故只分析荷载传递效率${E_{\text{A}}}$：

式中：${F_{\text{p}}}$为由于土拱效应传递到拱脚的力；W为总推力。根据图 4所示的双拱模型，土拱结构可分为端承拱和摩擦拱。因此，式（2）中的${F_{\text{p}}}$可以表示为

式中：${N_{\text{l}}}$，${N_{\text{r}}}$和${f_{\text{l}}}$，${f_{\text{r}}}$分别为端承拱的左右拱脚处法向阻力和摩擦拱的左右拱脚处切向阻力。为了模拟土体试样受到的水平位移，采用接近于实际横向土体扩散过程的位移控制式加载方式来模拟土拱形成演化过程。对顶部光滑墙体单元施加缓慢且均匀的速度以产生土体相对位移，而顶部墙体单元施加于土体试样的力就代表总推力。因此，式（2）可进一步改写为

式中：${W_{{\text{top}}}}$为顶部墙体单元施加于土体试样的总推力。

2.   水平土拱作用机理与发育规律

为了全面了解砂土水平土拱效应，采用位移控制式加载方法，并通过DEM模拟来观察土拱的形成过程，揭示其作用机理。鉴于篇幅所限，主要探讨拱跨比$\alpha = 3$的情况，而关于拱跨比的影响将在第3节进行详细分析。

2.1   应力重分布

水平土拱的宏观行为由其细观结构和力学特性决定，土体中的应力分布规律反映了土拱效应的作用机理。通过在土体试样区域内设置测量圆，统计记录其应力张量，可以得到主应力的方向和大小分布。大主应力旋转角度记为θ：

式中：${\boldsymbol{\sigma} _{xx}}$，${\boldsymbol{\sigma} _{yy}}$和${\boldsymbol{\sigma} _{xy}}$分别为水平、竖向和切向的平均应力张量。如图 5所示，通过提取的主应力大小和旋转角度θ来绘制应力十字分布图，其中应力十字的颜色表示第二应力不变量的相对大小。在拱脚上侧一定区域内，大主应力明显向跨中偏转，即大主应力方向由竖直方向旋转到水平方向，表明“应力拱”的形成。

图  5  应力十字分布

Figure  5.  Distribution of stress cross

下载: 全尺寸图片 幻灯片

以未加载Δ=0和加载位移Δ=2a为例，分析加载前后的应力分布变化，结果如图 5所示。在未加载Δ=0状态下，密砂的应力分布较为均匀，而松砂的应力分布略有不规则，这主要是由于松砂的颗粒分布不均匀导致的局部应力集中。随着加载位移增加，在加载位移Δ=2a状态下，大主应力方向逐渐形成“应力拱”（红色虚线），拱顶处主应力方向偏转90°，大主应力变为水平方向，小主应力为竖直方向。两侧拱脚处，大主应力近似竖直，小主应力近似水平。同时，由于去除了模型下部边界的位移约束而形成自由面，桩前颗粒位移分布呈现为拱形，即“位移拱”（绿色虚线），且“位移拱”方向与“应力拱”方向相反。

此外，为了量化拱跨中线上的应力分布，在拱脚中心连线以上20 cm范围内均匀布置了20个测量圆，以获得土拱结构上部、内部与下部的应力分布情况。通过编写Fish函数提取每个测量圆的水平应力${\sigma _x}$与竖向应力${\sigma _y}$并导出数据，密砂和松砂在未加载Δ=0和加载位移Δ=2a状态下的应力分布如图 6所示。在Δ=0的时候，密砂的${\sigma _x}$和${\sigma _y}$分布较为均匀而松砂却杂乱无章，这与前述应力十字分布图的分析一致。当Δ=2a时，无论${\sigma _x}$还是${\sigma _y}$分布，在密砂和松砂中都表现出相似的规律性：在远离拱脚的区域，${\sigma _y}$大于${\sigma _x}$，随着与拱脚的距离减小，${\sigma _x}$先增大后减小，而${\sigma _y}$则是基本呈现减小的趋势，并在靠近拱脚的一定范围内，${\sigma _x}$开始大于${\sigma _y}$。由于受到拱脚的阻挡，在一定距离处${\sigma _y}$和${\sigma _x}$发生偏转，土拱效应发挥作用，土体水平位移受到限制，荷载通过土拱向两侧拱脚传递。具体的${\sigma _y}$与${\sigma _x}$随水平位移的演化过程将在下文具体阐述。

图  6  拱跨中线上σ_x与σ_y分布

Figure  6.  Distribution of σ_x and σ_y along midline of arching span

下载: 全尺寸图片 幻灯片

2.2   力链及高应力颗粒分布

在砂土类颗粒物质中，砂土颗粒间相互接触会形成传递荷载的通道，即接触力链。通过分析接触力链的变化情况，可以从细观层面研究加载位移增加过程中砂土骨架的受力情况。接触力链可以表征接触的大小和方向，力链的粗细表示相应接触力的相对大小，而力链的密集程度反映了接触的数量。如图 7所示，给出了Δ=0和Δ=2a状态下的密砂和松砂的力链分布情况。当Δ=0时，密砂的接触分布比较均匀而松砂表现为局部接触较大，并且在密砂中力链更粗更密集。当Δ=2a时，它们的表现比较相似，力链网络逐渐发生重构，强力链逐渐沿着大主应力方向汇聚，此时密砂和松砂的最大接触力分别为2.59，2.32 kN。

图  7  力链及高应力颗粒分布特征

Figure  7.  Distribution features of contact force chains and high-stress particles

下载: 全尺寸图片 幻灯片

尽管接触力链可以在一定程度上反映土拱效应，但从系统分析颗粒的应力分布特征出发，可以更深入地揭示土拱效应的作用机理。在PFC^2D中通过编写FISH函数筛选出高应力颗粒，其基本的实现步骤可以总结为：①使用ball.stress命令获取所有颗粒的应力状态；②使用tensor.prin命令获取所有颗粒的主应力；③使用comp.x命令获取所有颗粒的大主应力；④根据所有颗粒的大主应力和颗粒总数计算全局平均主应力；⑤根据单个颗粒的大主应力和全局平均主应力的大小来筛选高应力颗粒。基于以上步骤筛选出的高应力颗粒如图 7所示，其分布特征与强力链的分布基本相同。此外，高应力颗粒的下边界可以直观描述土拱拱腹的形态，并可用于量化土拱结构的净矢高，即拱腹与拱脚间高差。图 7中红色虚线表明，松砂中的土拱净矢高较密砂更高。

2.3   细观形成过程

为研究土体水平位移与土拱形成过程的关系，选取加载位移Δ分别为0.2a，0.6a，a，2a四种情况进行分析。采用分层压缩法的同时将试样分为20层进行分组染色，以直观反映加载前后颗粒位移的变化。此外，在拱脚上侧区域布置测量圆监测应力张量的变化，获取土拱形成过程中的应力分布规律。

DEM模拟结果如图 8所示。图 8中左侧为表征应力重分布的应力十字和筛选出的高应力颗粒，右侧显示了颗粒位移分布。分析图中应力分布和颗粒位移结果发现，与松砂相比，土拱在密砂中的形成更快，且形成的土拱更加稳定。在Δ=0.2a时，应力在拱脚上方范围内立即出现明显偏转，竖向应力从拱跨中线向两侧逐渐增大，在两侧拱脚处应力张量达到最大；而横向应力从两侧向拱跨中线逐渐增大，在中线处大主应力方向已近乎水平。此时密砂中的高应力颗粒已在拱脚上方形成拱形连接，但松砂中的高应力颗粒分布相对较松散。砂土颗粒间抗剪强度逐渐发挥，初始土拱形成。在Δ=0.6a时，随着加载位移的增加，主应力旋转变得更加均匀，且颗粒位移较小，说明土拱形态正处于稳定发展阶段。此时，密砂中的高应力颗粒已经紧密连接成拱形，表明密砂中已形成了完全稳定的土拱，土拱效应作用充分发挥。在Δ=a时，土拱的挡土作用开始减弱，可以明显看出拱脚下方颗粒从拱脚间挤出，此时试样开始发生塑性变形，颗粒位移呈现中间大、两侧小的分布特征。结合应力偏转的区域分析可知，较大的位移差异主要位于“应力拱”下方，而“应力拱”上方近乎均匀移动。因此，通过颗粒位移的分布可以确定土拱效应的影响范围。与密砂相比，松砂中的土拱在Δ=0.6a~2a期间稳定发展并持续发挥作用，其形成过程更为缓慢。

图  8  土拱形成过程中的应力重分布和颗粒位移

Figure  8.  Stress redistribution and particle displacement during formation of soil arching

下载: 全尺寸图片 幻灯片

鉴于水平应力${\sigma _x}$相较于竖向应力${\sigma _y}$更能反映土拱的形成过程，结合拱跨中线上${\sigma _x}$分布进一步说明土拱效应与加载位移的关系，具体如图 9所示。在Δ=0.2a时，密砂的${\sigma _x}$值比较小，峰值不明显，而松砂中几乎未出现应力偏转。当Δ=0.6a和Δ=a时，密砂表现出显著的峰值现象并保持稳定，说明密砂中已经形成稳定的土拱结构，而松砂则开始逐步呈现弱峰值。随着加载位移持续增加，试样逐渐产生塑性变形，作用在试样上的推力持续增大。在Δ=2a时，虽然两者的${\sigma _x}$均有所增大，但应力分布特征变化较小，这表明荷载传递效率基本保持稳定，已经形成稳定的荷载传递结构。

图  9  不同加载位移时拱跨中线上σ_x分布

Figure  9.  Distribution of σ_x along midline of arching span at different loading displacements

下载: 全尺寸图片 幻灯片

尽管密砂和松砂的${\sigma _x}$大体分布趋势相似，但两者峰值出现的位置不同。从图 9中可以看出，密砂的峰值一般出现在距离拱脚中心连线以上8 cm处，而松砂则在11 cm处。基于${\sigma _x}$的峰值位置，可以量化土拱结构的矢高，即拱顶与拱脚间高差。结合上文通过高应力颗粒描述的净矢高，可以进一步得出砂土水平土拱的土拱厚度的演化规律。然而，土拱厚度的分析并不是本文的研究重点，故不进行深入讨论。

2.4   组构分析与配位数

组构分析是散体力学中的重要分析手段，主要是利用颗粒间接触进行统计分析力的方向。由于砂土类散体材料在模型试验中内部细观组构演变规律与其宏观强度及变形特性密切相关，可以通过接触组构的变化来表征模型试验中的土拱效应。统计收集图 4中布置的测量圆接触信息，以角度5°为间隔进行组构分析。从图 10可以看出，密砂和松砂在Δ=0时的分布趋势存在较大差异，且密砂分布较松砂更加均匀。当Δ=2a时，通过接触组构的变化可从细观角度揭示，土拱效应导致两种砂土体内部接触力方向出现了类似的偏转规律，但偏转大小不完全一致。组构分析能够反映接触方向和大小的规律性变化，间接揭示了水平土拱的作用机理。

图  10  密砂与松砂的组构各向异性

Figure  10.  Fabric anisotropy of dense and loose sands

下载: 全尺寸图片 幻灯片

进一步地，考虑到拱脚处的复杂应力状态，对拱脚处的配位数进行了具体分析。配位数C_n是研究颗粒体系接触特征的重要指标，与结构系统的稳定性有关，其定义为每个颗粒的平均活动接触数，计算公式如下:

式中：N_b为对测量区域中有质心的物体计数；$n_{\text{c}}^{(b)}$为颗粒b的活动接触数。在拱脚两侧各设置半径为0.5a的测量圆，并取两个测量圆的平均值作为拱脚处的配位数。图 11显示了拱跨比$\alpha = 3$，5情况下密砂和松砂的配位数分布。对密砂而言，整个加载过程配位数比较稳定，而在松砂中，配位数在前三步增长速度较快，之后趋于稳定。结合上文对土拱细观形成过程的分析可知，密砂中土拱形成较快，所需相对位移较小，而松砂中土拱的形成则较为缓慢。

图  11  拱脚处的配位数

Figure  11.  Coordination number at arching foot

下载: 全尺寸图片 幻灯片

值得强调的是，若考虑抗转动系数，试样的力学配位数会随着抗转动系数的增大而减小，而颗粒材料的组构各向异性程度则会随之增大^[19]。鉴于两种砂的土拱形成过程存在较大差异，因此密砂和松砂的土拱效应性能也将呈现较大区别。第3节将通过参数研究深入分析土拱效应的演变规律。

3.   土拱效应演变规律分析与比较

尽管许多学者将观察到的土体的应力偏转及荷载传递归因于土拱效应，但以往的研究对于砂土水平土拱效应的动态演化关注还相对较少。在上文分析土拱形成全过程的基础上，进一步揭示宏-微观参数对土拱形成过程的影响，并分析土拱效应的演变规律。

3.1   相对位移速度的影响

将加载速度分为10组，其中1^-1对应加载速度2a·s^-1，10^-1对应2a/10·s^-1，以此类推。为消除其他影响因素，除加载速度外的其他模拟参数保持不变。通过比较密砂和松砂在10种不同加载速度下的推力动态发展特征，发现当加载速度低于2a/3·s^-1时，曲线整体表现出相同的趋势，产生的推力水平相仿，此时相对位移速度的影响变得不那么显著。需要说明的是，由于未考虑抗转动能力，调整摩擦系数只能提高试样的峰值强度，而几乎不影响其残余强度^[20]，相应的发展规律也会在土拱推力分布趋势中得到反映。

由于篇幅限制，选取3组加载速度2a/5·s^-1，a·s^-1，2a·s^-1作进一步分析。图 12描述了3组加载速度下产生的推力W_top以及因土拱效应传递到拱脚的荷载F_p随加载位移的动态演变规律，这与文献[21]的结论一致。可以发现，无论加载速度大小，W_top的变化均与F_p基本一致。考虑到加载速度低于2a/3·s^-1时产生的力较为稳定，本文着重分析加载速度为2a/5·s^-1的情况。如图 12所示，发生桩土相对位移时土拱效应随即产生，根据力-位移曲线特征，可将密砂中土拱发育分为3个演化阶段：①S1近似线性增长并迅速强化达到峰值，说明密砂中土拱在较小的相对位移条件下形成，土拱荷载传递效率快速提高；②S2为峰值过后的衰减阶段，由于密砂的应变软化特性，土体强度持续减小，最后趋于砂土的残余强度，进而形成稳定的荷载传递结构；③S3曲线较平缓，随着加载位移的持续增加，土体不断被压实，土拱结构持续调整，土拱效应作用充分发挥。上述土体水平成拱过程的动态演化清楚地表明，相对位移对土拱特性产生了根本性的变化。不难发现，密砂的3个演化阶段与其剪切性状，即应变软化现象相似。同样，松砂的变化曲线也可以观察到与应变硬化相似的演化阶段，其中转折阶段S2对应变为转折点P2。

图  12  3组加载速度下的W_top和F_p演变

Figure  12.  Evolution of W_top and F_p under three loading velocities

下载: 全尺寸图片 幻灯片

对比密砂和松砂的演变规律可得，由于密砂S1阶段需要的相对位移更小，其成拱速度快于松砂。进一步可以发现，除了密砂中的S2阶段，两种砂在各阶段的力均随加载位移增大而逐步增强，这主要因为加载位移的增加会影响试样内部孔隙率的变化，从而提高试样的剪切强度，说明在加载过程中土拱的强度并非恒定的，而是处于一个动态的被强化过程。

为了进一步阐释加载速度对土拱的影响，图 13描述了加载速度为2a/5·s^-1时的端承拱与摩擦拱的关系。可以看出，与端承拱的法向阻力[N_r，N_l]相比，摩擦拱的切向阻力[f_r，f_l]几乎可以忽略不计。这表明[N_r，N_l]的影响对土拱的稳定性具有决定性作用，即端承拱对土拱效应起主导作用。将图 13，12中的W_top和F_p曲线对比可知，[W_top，F_p]和[N_r，N_l]几乎呈现相同的发展趋势。同时，从曲线的重合可以看出，[N_r，N_l]和[f_r，f_l]服从对称分布，最后形成的土拱是一个对称拱结构。通过对3种典型加载速度的分析得出，无论加载速度大小，[W_top，F_p，N_r，N_l]的发展趋势一致，且[f_r，f_l]相对[N_r，N_l]总是可以忽略不计。

图  13  加载速度为2a/5·s^-1时的端承拱和摩擦拱受力

Figure  13.  Forces of end-bearing arching and friction arching at loading velocity of 2a/5·s^-1

下载: 全尺寸图片 幻灯片

此外，作为土拱效应性能的常用评价准则，进一步研究了荷载传递效率E_A随加载位移的演变规律，结果如图 14所示。总体上，E_A与[W_top，F_p，N_r，N_l]呈现相同的发展趋势，且S1、S2(P2)、S3也与E_A的演变过程基本一致。土拱在S1阶段达到极限承载能力，而后维持一个稳定的残余承载能力，这与砂土的剪切性状也保持一致。

图  14  加载速度为2a/5·s^-1时荷载传递效率E_A演变

Figure  14.  Evolution of efficacy E_A at loading velocity of 2a/5·s^-1

下载: 全尺寸图片 幻灯片

为定量分析不同状态下的土拱效应荷载传递效率E_A，本文采用小提琴图的形式对Δ=0~2a全过程的E_A进行了分析。小提琴图可以直观反映每一个位移量对应的E_A基本分布位置和形状，方框内红点代表平均值，方框底部和顶部为第一和第三四分位数（Q1和Q3）。如图 15所示，以拱跨比$\alpha = 3$情况为例，统计得到的全过程E_A平均值为0.7397，并且E_A的偏离程度较弱，主要分布在Q1和Q3的中段，集中在平均值附近。小提琴图分析不仅可以方便地获取均值等关键信息，还可以显示序列中的异常值，有利于数据清洗。考虑到所有数据波动范围和平均值大多在±5%误差区间内，因此可以合理利用小提琴图评估土拱的荷载传递效率E_A。

图  15  荷载传递效率E_A小提琴图及其在不同加载速度下E_A的变化曲线

Figure  15.  Violin plot of arching efficiency E_A and its variation curve under different loading velocities

下载: 全尺寸图片 幻灯片

基于上文分析，采用小提琴图统计获得的荷载传递效率E_A平均值来表征土拱形成发展全过程土拱效应性能，并据此分析了不同加载速度下的土拱荷载传递效率，结果如图 15所示。结果表明，两种砂土的土拱荷载传递效率均随加载速度的降低而略有提高，并且在加载速度低于2a/5·s^-1时基本不变。当加载速度从2a·s^-1改变到2a/10·s^-1时，土拱荷载传递效率在松砂中仅提高了约2.5%，而在密砂中仅提高约2%，可以发现加载速度对土拱荷载传递效率E_A的影响较小。

鉴于加载速度为2a/5·s^-1时计算成本较低且能保证稳定可靠的结果，选用此加载速度进行后续的参数分析。

3.2   拱跨比的影响

合理拱跨的确定是抗滑桩工程设计的重要内容之一，拱跨是决定土拱效应的关键指标，拱跨过大可能导致土拱效应失效，过小又易造成施工成本增加。鉴于本文拱脚尺寸不变，采用拱跨比来研究拱跨大小的影响，选取3~7的拱跨比进行了分析。图 16的结果表明，拱跨比对土拱荷载传递效率E_A影响显著，E_A随着拱跨比的增大近似线性减小，且相较于密砂，松砂的下降更为明显。同时，分析了加载速度为2a·s^-1时的情况进行对比。在加载速度为2a/5·s^-1和2a·s^-1时，两种砂的曲线变化近似平行，加载速度越快，土拱荷载传递效率E_A越低。但是，加载速度对E_A的影响很小，这与图 15的结论相同。

图  16  拱跨比对荷载传递效率E_A和端承拱承担比(N/F_p)的影响

Figure  16.  Influences of arching span ratio on arching efficacy E_A and efficacy of end-bearing arching (N/F_p)

下载: 全尺寸图片 幻灯片

此外，图 16还展示了拱跨比对端承拱承担比（N/F_p）的影响。可以看出，与E_A不同，无论加载速度及砂土的种类如何，端承拱承担比随拱跨比增加仅略有下降，且始终维持在90%以上，这也说明了端承拱在土拱效应中起主导作用。同时，密砂中端承拱承担比大于松砂，这从侧面说明密砂中土拱效应性能优于松砂。

3.3   摩擦系数的影响

鉴于摩擦系数对土体宏观强度有较大影响，本节主要分析摩擦系数的影响。如图 17所示，在50 kPa围压条件下，摩擦系数分别为0.5和1.0的密砂和松砂应力-应变曲线存在明显差异，无论摩擦系数取0.5或1.0，其对砂土的剪切性状特性都有显著影响，尤其是密砂峰值强度相差近50 kPa。然而，摩擦系数对砂土的残余强度影响较小。这表明，当摩擦系数取0.5或1.0时，砂土残余强度相对不敏感于摩擦系数变化，但是其峰值强度会受到较大影响。若考虑抗转动能力，则对峰值强度和残余强度影响均较大。

图  17  50 kPa围压下的不同摩擦系数的砂土剪切性状

Figure  17.  Shear behaviours of sands with different friction coefficients under confining stress of 50 kPa

下载: 全尺寸图片 幻灯片

鉴于[W_top，F_p，N_r，N_l]分布趋势一致，此处重点分析了推力W_top和荷载传递效率E_A。如图 18所示，摩擦系数对推力W_top演化发展趋势的影响较小，但由于摩擦系数改变了砂土的剪切特性，进而影响了推力W_top的大小。随着摩擦系数减小至原来的一半，推力W_top显著降低。与图 14的结果对比可得，虽然W_top的大小变化较大，但摩擦系数对E_A分布趋势的改变量较小。

图  18  不同摩擦系数的推力W_top和土拱荷载传递效率E_A演变

Figure  18.  Evolution of thrust W_top and arching efficacy E_A for different friction coefficients

下载: 全尺寸图片 幻灯片

4.   结论

本文进行了一系列砂土边坡桩间水平土拱效应DEM模拟，从端承拱和摩擦拱角度着重研究了密砂和松砂的土拱机理与演变规律，得到以下4点结论。

（1）水平土拱的宏观行为由其力学特性和微观结构决定，土体中的应力偏转、高应力颗粒分布规律及组构分析都能够反映土拱成拱过程和作用机理。

（2）土拱的形成和演变与相对位移及其速度的大小密不可分，根据相对位移的大小，两种砂土的[W_top，F_p，N_r，N_l，E_A]发展均经历了3个演化阶段，分别对应于密砂应变软化和松砂应变硬化的剪切性状。

（3）相比[f_r，f_l]，[N_r，N_l]的作用更大，且端承拱承担比（N/F_p）都在90%~95%，端承拱在土拱效应中占主导地位。

（4）对荷载传递效率E_A而言，拱跨比的影响最大，加载速度和摩擦系数的影响较小。