桩间三维土拱效应颗粒流数值模拟及其演化规律

王桂林，冉龙宝，张亮

(重庆大学 a.土木工程学院；b.库区环境地质灾害防治国家地方联合工程研究中心，重庆 400045)

土拱效应的研究与发展已经有100多年的历史,Roberts在1884年发现和提出“粮仓效应”，可看做拱效应的开端。土拱效应最初由Terzaghi[1]提出，并通过活动门实验证明了土拱效应在岩土领域的存在。近年来，很多学者采用现场模拟试验[2]、合理拱形理论分析[3-4]和数值分析[5-7]等方法将土拱效应理论应用到岩土工程领域并取得了很好的效果。其中，抗滑桩间土拱效应尤为明显。

随着计算机技术的普及与发展，运用数值模拟进行抗滑桩间土拱效应作用机理研究取得了长足发展。目前的研究多建立在已有合理拱轴线理论[8]上，分析各因素对桩间土拱力学效应的影响。不少学者利用有限差分软件和离散元软件从二维角度分析了桩间距、土体内摩擦角、粘聚力、孔隙比等因素对抗滑桩间土拱效应的影响。Chen等[9]学者使用Flac有限差分软件对群桩在侧向土体运动作用下的成拱机理进行了探讨，解释了应力从土传递到桩的过程。张建华等[10]建立起抗滑桩结构的有限差分模型分析了桩间距对桩间土拱效应的的影响，研究结果对抗滑桩间距的设计有一定的参考价值。向先超等[11]利用颗粒流方法，研究了抗滑桩截面大小、间距、桩土相对变形速度和土体颗粒粒度组成对土拱效应的影响，并揭示了这些因素对土拱的极限承载能力、残余承载能力和桩土荷载分担比的影响规律。Sun等[12]采用二维颗粒流数值模型，研究了相邻悬臂式抗滑桩间土拱的形成和破坏过程。以上研究均将桩间土拱效应简化为二维平面问题，但实际上桩间土拱既有水平土拱的作用，又有竖向摩擦拱的存在，三维土拱效应明显。近年来，对桩间土拱效应的三维特征虽有研究[13-17]，但对三维土拱效应的发展及土拱厚度演化规律的研究不够深入。利用颗粒流数值模拟软件PFC3D，分析抗滑桩间土的三维土拱效应的形成发展及破坏过程，并对土拱厚度的演化规律进行深入研究。

1 基本模型及参数标定

1.1 抗滑桩三维模型

模拟抗滑桩桩间土拱的三维效应，将桩视为刚性墙体，桩后黏性填土视为离散圆球颗粒并将颗粒间链接设为接触黏结(contact bond)，利用墙体来限制颗粒运动作为边界约束条件(图1)，利用加载墙的运动来模拟土拱的形成演化过程。

模型参照董捷等[15]、王振强[18]等的桩间三维土拱效应物理模型试验，长、宽、高取为3.2 m×2.0 m×2.0 m，其中，墙后填土区尺寸为2 m×2 m×1.8 m，抗滑桩截面尺寸为0.3 m×0.3 m，桩间距取1 m，模型如图1所示。模拟时，按照指定孔隙度生成填土范围内的圆球颗粒并在自重作用下压实稳定。文中模型宏观参数取值详见表1。

图1 抗滑桩三维模型简图Fig.1 Three-dimensional model sketch of anti-slide

重度/(kN·m-3)弹性模量/kPa内摩擦角/(°)粘聚力/kPa202×1032010

1.2 参数标定

模拟常规三轴试验试件尺寸为高10 cm、半径2.5 cm的圆柱形，模型上下设置边界墙(wall)单元模拟加压过程，侧壁用圆柱形墙体单元(cylinder wall)模拟橡胶套筒并通过伺服保持恒定围压，如图2所示。颗粒间接触为设为接触黏结(contact bond)，通过分别施加50、75、100 kPa的围压得到材料的应力-应变关系曲线，如图3所示，对应的偏应力峰值强度分别为84、118、146 kPa。材料在3种围压下的摩尔强度包络线[19-20]如图4所示。

由摩尔应力强度包络曲线图可得到数值模型的抗剪强度参数粘聚力和内摩擦角分别为10 kPa和22°。最终确定模型的细观参数如表2所示。

图2 三轴试验试件模型Fig.2 Triaxial test specimen

图3 偏应力-应变曲线Fig.3 Deviation stress-strain

图4 摩尔强度包络线Fig.4 The envelope line of Mohr

弹性模量Ec/Pa刚度比(kn/ks)摩擦系数最大最小半径比2×1062.00.201.66法向链接强度/Pa切向链接强度/Pa孔隙度/n密度/(g·m-3)8×1038×1030.292 000

2 三维土拱的形成演化规律

采用三轴剪切模拟实验得到的细观参数，建立抗滑桩间土拱效应三维数值模型，取颗粒最小半径为0.014 m，共使用约17×104个球形颗粒模拟黏性土颗粒(图1)。模拟过程中，通过对加载墙施加大小为 0.01 m/s方向水平(y负方向)的速度来模拟下滑力对桩后土体的作用，撤去桩间挡土板计算至土体破坏。

2.1 土拱的形成与发展

土拱效应可解释为桩土产生相对位移引起的桩后土体应力重分布现象，故可用土体颗粒位移量的差异来表示土拱的演化过程。付海平等[21]在桩承式路堤土拱效应的研究中验证了桩土相对位移方法验证土拱效应演化的可行性。本研究利用PFC内置fish语言，获得桩后土体产生的最大位移d，并将其分成15等份，绘制桩后土颗粒相对位移三维等值图。

撤去桩间挡土板并给加载墙施加水平方向的速度。由于加载墙速度恒定，因此，加载墙位移与计算时步成正比关系。计算100时步(局部受压阶段)时，由图5(a)、图6(a)可以看到，临空面附近土体开始产生偏向桩外侧的位移且桩底土体位移要多于上部，这是由于抗滑桩底部应力较大，撤去挡土板后临空面位置出现较大的应力集中，颗粒之间相互作用力较大，桩间临空面底部颗粒首先产生相对较大的变形位移，同时，桩背侧土体由于桩的约束，变形几乎为零，桩间临空面与桩后土体变形存在差异。可以发现，此时桩间并没有形成弧链桩的土拱形状，表明土颗粒仍处于挤压致密的过程；继续运算至500时步时，根据图5(b)、图6(b)可以看到桩间土颗粒与桩背土位移之间出现了较大的分异。具体表现为自抗滑桩内侧边缘起形成较为直观的弧链状位移分层。同时，通过图6(b)可看到相对位移的差别沿高度也有很好的辨识度，认为此时桩间土拱已经开始形成；计算至2 000时步时，通过图5(c)、图6(c)可观察到相对位移为(2/15～3/15)d的土颗粒分层明显向土体内部延伸，计算至4 000时步时，土颗粒相对位移分组与2 000时步时相比没有太大变化，表明此段时间内土拱处于稳定发展阶段。

图6 抗滑桩桩间对称面(x=1 m)颗粒相对位移Fig.6 Relative displacement of particles on symmetrical plane (x=1 m) between anti-slide

2.2 土拱的破坏

继续计算至14 000时步时，由图7(a)、(b)可以看到，桩间弧链状土拱层不再完整，从桩间对称面上可以看到临空面附近出现较大的断裂，表明土体已超过极限平衡状态，原来形成的稳定土拱遭到破坏。

图7 计算14 000时步时的颗粒位移Fig.7 Displacement of particles ( timesteps=14 000 )

继续运算，直至桩后土体出现垮塌，如图8所示，可以看到数值模拟得到的桩后土体破坏情况与物理模型实验所得结果相同，均表现为中下部土体率先开裂并垮塌形成空腔，顶部出现悬链状裂缝并伴有土体下陷垮塌的迹象。产生这一现象的原因是，在下滑推力不断增加的过程中，中下部土体产生的应力相对较大，而不同位置处相同土体的极限抗剪强度是相同的，因而中下部土体将率先出现开裂。临空面附近土体由于缺少约束，在自重作用下将不断地垮塌，最终形成贯通的滑动面。

图8 物理模型实验与数值模拟实验破坏结果对比Fig.8 Comparison of damage results between physical model experiment and numerical simulation

为了更直观地分析三维土拱的演化过程，在桩后高度z为0.5、1.0、1.5 m处分别设置应力监测球，监测桩后y方向的应力随加载墙位移的变化情况，与此同时，用抗滑桩承受的土压力与加载墙承受的压力之比表示荷载传递效率，荷载传递效率越高，则桩间三维土拱效应越显著。监测结果见图9。

图9 桩后y方向应力及荷载传递效率Fig.9 y-direction stress and load transfer efficiency behind

图9中桩后不同深度处y方向的应力监测显示，填土底部(z=0.5 m)处的应力监测在加载墙位移达到12 mm后开始急剧减少，该位置形成的土拱此刻已达极限状态；抗滑桩中部(z=1.0 m)处的应力在加载墙位移达18 mm左右时开始下降，抗滑桩上部(z=1.5 m)在加载墙位移达25 mm时才开始下降，表明土拱的破坏是从抗滑桩底部开始，逐渐向上发展的。

2.3 土拱效应阶段划分

根据图9所示的桩后y方向应力及荷载传递效率监测结果，结合土拱的形成演化过程分析，将土拱的形成发展分为3阶段：土拱初步形成阶段、土拱稳定发展阶段和土拱破坏阶段。

土拱初步形成阶段：桩土相对位移较小，桩后土体局部受压，随着桩土相对位移不断扩大，桩背侧土体受到阻拦，且范围不断扩大，临近土体相互锲紧，并向后侧土体发展，形成具有较高承载力的且沿桩身均有分布的拱形结构，将此阶段定义为土拱初步形成阶段。此阶段的特点是桩后土体应力增加较快，荷载传递效率呈增加的趋势。此阶段出现在加载初期，加载墙位移不足2 mm时。

土拱稳定发展阶段：土拱雏形形成后，随着加载墙位移的继续增加，拱后土体密实度增加，土拱作用进一步增强。此阶段的特点是桩后应力增长变缓，荷载传递效率基本保持不变。该阶段持续时间较长，出现在加载墙位移达约2～15 mm期间。

土拱破坏阶段：土拱稳定阶段发展到后期，随着加载墙位移不断增大，土拱进一步挤密，达到土体抗剪强度，土拱发生破坏。此阶段的特点是桩后应力突然减小，荷载传递效率呈现降低的趋势。该阶段出现在加载墙位移达到约15 mm后。

3 土拱厚度的演化规律

3.1 土拱厚度确定

土拱效应的产生与桩土间的相对位移密不可分，拱身部分土体在锲紧作用下密实度较高，而临空面附近由于没有约束，土颗粒位移较大，拱身和拱前土体可能会产生较大的相对位移，发展到后期甚至出现裂纹或形成贯通的滑动面，故可用土体颗粒相对位移来确定土拱前缘。

取土拱稳定发展阶段(s=5 mm)的土拱效应进行分析。现选取距桩底部z=0.1 m处的水平截面，用相对位移确定土拱范围的前缘。土拱前缘开始于抗滑桩内侧边缘，向内部土体延伸，如图10所示。土拱前缘线大致呈悬链状，左右对称。

图10 由相对位移确定土拱前缘Fig.10 Determination of leading edge of soil arch

悬臂式抗滑桩土拱效应的空间特征十分显著，按照土拱的作用面将其分成水平拱、竖向拱和临空面拱。其中，水平面拱为大主应力偏转形成的大主应力拱并认为拱轴线上主应力处处相等，即拱轴线也是大主应力等值线。基于此假定，可以利用大主应力等值线确定土拱后缘。通过导出颗粒流数值模型中的应力张量，利用MATLAB求解最大主应力并生成二维等值线图，如图11所示。可以发现最大主应力等值线在临空面附近发生了偏转，等值线呈现出一系列的弧形状排列，且靠近临空面时等值线较密集，远离临空面达到一定距离后，等值线变得稀疏直至不具有明显特征。以桩后应力集中位置为起点，参考等值线的稀疏分布及规律，向两边延伸形成的光滑弧形拱作为土拱的后缘。

图11 最大主应力等值线图确定土拱后缘Fig.11 Determination of trailing edge of soil arch by contour map of maximum principal

将由上述方法确定土拱的前缘和后缘综合绘制在图12上，可以看到土拱前缘后缘形状基本相似，形成与桩宽大致相等(0.3 m)的等厚土拱。

图12 土拱厚度

在桩间对称面上z=0.1 m位置处设置一排监测点监测各点的x方向和y方向应力情况。当土拱效应发挥时，应力沿拱轴线方向传递到桩后，在桩间对称面上表现为x方向应力突变，在土体内部，由于加载墙的移动，y方向应力基本保持不变，靠近土拱效应发挥区域时，应力发生偏转，y方向应力将突然变小，利用这一特性可以检验前述方法的合理性。

由图13可以发现，x方向应力变化趋势减缓发生在距离临空面约0.3 m处，y方向应力变化趋势减缓发生在距离临空面约0.8 m处，与图12确定的土拱范围基本一致，说明采用前述方法确定土拱效应的拱厚度是可行的。

图13 z=0.1桩间对称面上应力监测Fig.13 Stress monitoring on symmetrical plane between

3.2 土拱厚度随深度的变化规律

在土拱稳定发展阶段，计算时步为4 000时步时，(加载墙位移约5 mm)分别取距桩底部0.1、0.5、1.0、1.5 m处的水平截面按照3.1节的方法确定土拱厚度，如图14所示。

图14 加载墙位移5 mm时不同深度土拱厚度Fig.14 Thickness of soil arch at different depths when loading wall displacement 5

根据图14可知，土拱厚度是随深度变化的。为了更明显地显示土拱厚度随深度的变化情况，通过分析多个不同深度水平剖面的土拱厚度分布，作出桩间对称面土拱厚度随深度变化图，如图15所示。在桩底位置时，土拱厚度相对较薄，这是因为桩底位置除了抗滑桩对土体有阻挡约束作用外，部分应力通过下层土体传递到桩前土(模型中由底部固定边界摩擦承担)，因此，土拱效应相对较弱；土体更容易产生较大位移，土拱效应更加明显。桩顶位置由于土体应力较小，颗粒产生的位移量更小，土拱厚度相对较薄。由于采取的模型中设置了加载墙以便模拟土拱的破坏，使得桩土相对位移量比自然状态下的边坡桩土相对位移更大，土拱的厚度与实际情况相比也要偏大。

图15 桩间对称面上土拱沿深度变化图Fig.15 The change of soil arch along the depth on the symmetrical surface between

3.3 土拱厚度随时间的变化规律

由于数值模型设置的加载墙一直在运动中，表现为桩后下滑力一直处在变化中，故土拱厚度随时间的变化规律更多的表现为土拱厚度在演化过程中的变化规律。

取计算时步为100、500、2 000、4 000、14 000分析桩后土拱效应随时间的演化。以抗滑桩中部(z=1.0 m)水平截面为例进行分析，结果如图16所示。

图16 z=1.0 m不同时刻土拱厚度Fig.16 Soil arch thickness at different time(z=1.0 m)

图16的结果表明，土拱的厚度在土拱的形成演化过程中也是在变化的。在最开始计算的100时步时，桩后土体处于局部受压状态，桩间土体位移相对较小，土拱效应较弱，形成的土拱仅在临空面附近0.35 m范围内，且厚度较薄；计算至500时步时，桩间土体处于土拱初步形成阶段，桩间土体相对位移增大，土拱的前后缘均向土体内部延伸，土拱的影响范围扩展至临空面附近0.8 m范围内，土拱厚度增加；计算至2 000时步时，桩间土体处于土拱稳定发展阶段，土拱前缘向土体内部产生较小的延伸，土拱厚度变薄，继续计算至4 000时步时，土拱厚度进一步变薄，土拱后缘朝临空面附近靠近，产生这一变化的原因是，土拱稳定发展过程中，土拱形状基本确定，前缘经锲紧作用压密，而后侧土体还有压缩的空间，随着加载墙位移的增加，土拱影响范围内的土体不断压密导致最大主应力等值线向临空面附近靠近，土拱厚度变薄；继续计算至14 000时步时，土拱处于破坏状态，由图16(b)可看到，此时的土拱范围与土拱稳定发展阶段的大不相同，前缘土体已经破坏，土拱前缘已不存在。

综上所述，在土拱的形成演化过程中，土拱的厚度具有明显的空间与时间变化规律。目前，对于土拱厚度的认知公认的是假定土拱厚度等于桩宽或者为桩宽的一半。这种做法显然具有一定局限性，首先，桩后不同深度处或者不同阶段时土体受力状态不同，所形成的土拱厚度也不同。当土拱处于稳定阶段时，土拱的厚度与桩宽相差不大，表明在精度要求不高的情况下，将桩宽作为土拱厚度也是可行的。

4 结论

依据已有的桩间三维土拱效应物理模型试验，选取合理的模型尺寸与参数进行三维颗粒流模拟研究，得出以下结论：

1)把抗滑桩当作刚性桩时，撤去抗滑桩间挡土板后桩底的土体首先产生位移，之后不断向土体内部和桩顶发展；土拱的破坏首先出现在桩底位置，并逐渐向桩顶发展。

2)土拱厚度随深度是变化的，具体表现为桩底位置土拱厚度较小，中部土拱厚度最大，桩顶位置较薄。

3)土拱厚度在滑体推力施加过程中是变化的。随着加载墙位移的增加，土拱厚度经历了由薄变厚再变薄的变化过程。