0.   引言

采用碎石桩对软土地层进行加固形成碎石桩复合地基,能够有效提高其地基承载力,减小沉降变形,并且具有施工快速简便、造价经济等优点。大量工程实践应用证明了碎石桩复合地基是一种非常适合中国国情的地基处理方法^[1]。国内外许多研究者采用理论和试验等方法对碎石桩复合地基的变形和承载原理进行过较为深入的研究^[2-6]。

由于涉及到离散碎石和软土两种特性迥异的材料之间的相互作用,理论解析方法难以完全揭示碎石桩的承载机制,因此很多学者也采用数值模型对碎石桩复合地基中的桩土相互作用展开模拟研究。陈建峰等^[7]建立数值模型精确模拟了碎石桩单桩夯扩和载荷试验加载过程,其模拟和实测的夯扩碎石桩桩顶、桩端荷载-沉降曲线均吻合良好。蒋敏敏等^[8]通过数值模拟分析了高速公路碎石桩复合地基在桩体施工、路堤填筑、运行期全过程和地震动荷载作用等问题。王复明等^[9]通过有限元模型研究了高聚物碎石桩承载特性。但目前大部分碎石桩数值模型均基于连续介质方法,难以准确描述碎石颗粒细观接触关系以及变形不连续性导致的桩-土复杂相互作用^[10]。离散单元法则以颗粒或块体为基本单元,单元之间没有变形协调的约束,能很好地模拟碎石及土体大变形和非连续位移的过程^[11]。Indraratna等^[12]提出了一种离散-连续介质耦合数值计算格式,利用二维PFC圆形颗粒模拟碎石桩,用FLAC有限差分网格模拟桩周土体,建立耦合数值模型成功地模拟了碎石桩单桩荷载试验。基于类似建模思路,谭鑫等采用基于二维Voronoi块体^[13]离散单元数值模型的方法,对散体材料单桩^{[11, 14]}以及复合单桩^[15]展开了竖向荷载下变形承载全过程的数值模拟研究,并据此提出了碎石桩单桩承载的典型变形及失稳模式。Gu等^[16-17]则基于PFC^3D采用球形颗粒模拟了土工格栅包裹碎石桩的无侧限压缩试验,其模拟结果与模型试验结果获得了比较好的吻合度。

采用离散单元方法建立的碎石桩数值模型无需复杂的本构模型及假设便能较好地反映其鼓胀变形和失稳特性,为碎石桩复合地基数值模型提供了一条有效途径。但是上述二维模型或者球形颗粒假设忽略了碎石颗粒本身三维形态对碎石桩力学特性的影响。较多研究^[18-20]指出颗粒材料的力学特性受到其离散颗粒形状、大小以及粒径级配的影响。因此本文在考虑碎石颗粒三维形态的基础上,将碎石桩视为刚性凸多面体集合,桩周软土视为理想弹塑性材料,采用连续-非连续介质耦合数值模型对碎石桩单桩受荷室内模型试验的变形及破坏全过程进行模拟。通过数值模拟结果有效地揭示了碎石桩宏观变形失效与其细观颗粒受力运动之间的联系。

1.   软土中碎石桩三维数值模型建立

1.1   三维碎石颗粒模型及接触

碎石颗粒的真实形态和级配直接影响到碎石桩体的力学特性,但数值模型不可能也不需要真实模拟每一颗不同形态的碎石。本文采取的方法是预先建立包含有限个典型碎石颗粒三维形态的模板库,在生成碎石材料时从库中随机选取三维颗粒形态,按照级配曲线指定颗粒大小混合后生成碎石堆积体。首先选取一定数量的典型真实碎石颗粒进行三维形态的简化和提取,在CAD软件中将三维碎石颗粒转换为凸多面体,然后导入PFC^3D程序建立碎石形态模板库。考虑到碎石桩体最可能受到的应力状态区间,可以认为碎石颗粒本身不发生破坏,因此直接采取刚性多面体模拟单个碎石颗粒,碎石桩体的变形由颗粒接触和颗粒移动产生。

颗粒间的接触模型如图1所示,在结合颗粒形态和粒径级配的基础上采用简单的线性接触模型即可较好地模拟碎石材料的力学特性。碎石颗粒间的变形和力的传递通过刚性块体接触的接触刚度（法向刚度k_n及切向刚度k_s）进行计算。接触的强度则通过强度参数（摩擦系数μ、黏结及抗拉强度）来控制。其中碎石摩擦强度一部分取决于颗粒接触之间的摩擦系数,一部分则通过不规则三维颗粒之间的咬合力提供。由于已经考虑了三维颗粒的几何形态,因此接触摩擦系数可以采用碎石材料的真实摩擦系数。散体碎石之间没有黏结材料,因此黏结和抗拉强度均可以设置为零。

图  1  三维碎石颗粒及其接触模型

Figure  1.  3D gravel particles and contact model

下载: 全尺寸图片 幻灯片

当块体接触内力未超过其接触强度时,模型通过颗粒间接触刚度产生弹性变形；而当接触内力超过其强度时,接触便发生破坏,颗粒将产生滑移、转动或者脱离。但不同于连续介质方法,离散单元法采用的是接触的微观力学参数。室内试验通常无法直接提供这些微观参数,所以仍然需要通过数次标定校核后才能获得相对应的微观参数。通过离散单元法模拟碎石桩体时,可首先根据碎石桩体的宏观力学特性比如压缩模量、摩擦角等预先选择微观参数,再进行多组数值试验（比如单轴、三轴、侧限压缩试验等）将模拟结果与试验获得的宏观力学特性进行比对。根据对比调整微观参数,使数值模拟结果与实验室结果趋于一致。有部分研究^{[13, 18, 20]}提供了一些半解析半经验的宏细观参数标定方法,可以用于初选参数,本文不再赘述。

1.2   模型试验与数值模型

本文选取了文献[3]记录的模型试验作为本文提出数值方法的模拟案例,来验证模型及方法的有效性。在图2所示的直径为500 mm,高度为500 mm的圆桶型模型箱中采用粒径在2～10 mm之间的细小碎石建立了单根模型碎石桩。其中2～5 mm粒径碎石含量占总质量37.5%,5～10 mm粒径碎石含量为59.5%,仅3%含量的粒径小于2 mm。碎石桩桩径100 mm,高度500 mm,利用预埋的薄壁钢管分批灌入并分层击实碎石料直至埋设高度。桩周填土是液限为33.8%,塑限为19.6%,含水率为22.8%的中塑黏土,其十字板剪切强度值为48 kPa。在模型碎石桩顶部安装了与桩径相同的刚性加载板,通过固定于加载架上的千斤顶进行竖向加载。采用多级加载方式加载（每级荷载0.7 kN）,直至模型碎石桩破坏。加载过程中记录了每个加载稳定阶段后的沉降,并在图2所示压力计位置进行了土压力测量。试验结束后,将熔化的石蜡注入破坏的模型桩体中,待石蜡凝固后小心去除桩周土体后,可直接观察到变形的桩体轮廓。

图  2  模型试验设置

Figure  2.  Setup of model tests

下载: 全尺寸图片 幻灯片

采用1.1节所述的三维碎石颗粒模型建立了与模型试验相同尺寸的数值模型如图3所示。其中PFC^3D碎石桩体采用了大概十余种典型的颗粒形态按照模型试验采用的颗粒尺寸与粒径分布随机混和而成。桩周土体则采用服从莫尔库仑破坏准则的FLAC^3D有限差分网格进行模拟。在PFC^3D桩体和FLAC^3D土体之间设置了Wall单元,分别作为耦合边界用来传递离散单元颗粒与有限差分网格之间的力及位移。土体边界除了自由顶面外均施加法向约束,桩体上下两端通过与桩体相同直径的Wall来施加底部约束以及顶部荷载。数值模型在重力的作用下达到初始平衡后才进行单桩加载。颗粒宏细观参数通过一系列三轴压缩试验模拟标定后列于表1,其中细观参数为模型输入参数,宏观参数为三轴压缩试验模拟获得的材料特性。土体参数列于表2。

图  3  数值模型设置

Figure  3.  Numerical model

下载: 全尺寸图片 幻灯片

表  1  碎石桩体的材料参数

Table  1.  Material parameters of stone column

块体及微观接触参数（输入）							宏观表象参数（计算）
接触杨氏模量Ec/Pa	接触刚度比ks/ks	接触摩擦系数μ	碎石粒径d/mm	相对密实度	碎石数量	块体重度γ/(kN·m-3)	压缩模量Eo/MPa	弹性模量Ec/MPa	泊松比$\nu$	表观摩擦角φ/(°)
2.0×109	2.0	0.55	3~12	0.8	8130	25	75	45	0.35	42

下载: 导出CSV 

| 显示表格

表  2  软土的材料参数

Table  2.  Material parameters of soils

cu/kPa	φ/(°)	E/MPa	$\nu$	$\rho$/(g·cm-3)
48	0	6.0	0.4	1700

下载: 导出CSV 

| 显示表格

2.   数值模拟结果与验证

图4为数值模拟的单桩荷载-沉降曲线（p-s曲线）与模型试验结果的对比,两者吻合较好。在单桩加载初期沉降达到2.5 mm之前p-s曲线基本呈线性增加,桩顶应力在较小的沉降量下即迅速增加到约500 kPa；之后随着荷载增大,桩体出现软化,沉降增长率逐渐增大,当沉降超过5 mm时,非线性变形明显发展；当桩顶应力超过600 kPa后沉降迅速发展,沉降超过25 mm后p-s曲线几乎竖直,表明模型碎石桩丧失了承载能力。数值模拟的p-s曲线显示的单桩承载力与模型试验施加的最大荷载（1073 kPa）非常接近,两者均略大于采用Branus公式计算的单桩极限承载力。

图  4  单桩载荷试验的p-s曲线

Figure  4.  p-s curves of loading tests on a single stone column

下载: 全尺寸图片 幻灯片

图5为数值模型在不同沉降阶段的桩体鼓胀变形（水平位移）及桩周土体应力沿桩身分布曲线图。其中,桩侧土体水平应力分布曲线与模型试验实测值对比见图5（b）,两者在最大鼓胀量和鼓胀长度上非常接近。图5（a）显示随着沉降的增加碎石桩全桩鼓胀变形量均逐渐增加,但鼓胀最剧烈部分限制在桩体上部两倍桩径的埋深范围内。桩体破坏后最大鼓胀变形量为5 mm左右,对应径向应变约为10%。由于上部加载板和下部桩体对鼓胀变形段内碎石颗粒的摩擦约束,鼓胀变形呈现出上下部小中间凸出的形态。模型桩侧土体的水平应力和竖直应力分布与桩体鼓胀变形分布较为一致。其中桩侧应力在鼓胀剧烈区域相应增加明显,桩体破坏阶段桩侧土体水平应力最大值接近200 kPa。由于鼓胀变形剧烈范围内的桩体压缩量较大,桩土相对位移也最为明显,因此部分竖直荷载通过桩土之间的剪切位移传递到桩周土体中。可以看到无论是在反映整体碎石桩整体受力变形的p-s曲线还是反映局部桩土相互作用的水平应力分布曲线上,数值模拟结果与模型试验监测结果均吻合较好,从而验证了本文所提出数值模型的合理性。

图  5  不同沉降下桩体鼓胀变形及桩周土体应力分布曲线

Figure  5.  Distribution of column bulging deformation and soil stresses under varying settlements

下载: 全尺寸图片 幻灯片

不同于连续介质数值模型计算的应力应变曲线,离散单元法建立的碎石桩体由于考虑了碎石颗粒的形状和离散特性,能够反映离散碎石颗粒的非连续运动特征,块体凸起点与土体网格接触处产生局部应力集中,导致p-s曲线和土体变形应力分布曲线并不光滑。

3.   碎石桩单桩承载过程中桩周土体受力变形特征

图6为数值模型在不同桩顶沉降时对应的桩周土体水平位移及应力场分布。其中碎石颗粒颜色分布代表着相应颗粒沉降量大小,可以看到桩体压缩主要产生在桩体上部大约2.5倍桩径以内的范围。该段区域内碎石产生剧烈鼓胀变形,从而引起桩周土体产生水平位移。土体水平应力在桩体鼓胀导致的水平挤压下增长剧烈,与初始应力场相比,鼓胀最剧烈处水平应力提高了将近10倍。

图  6  不同沉降下土体水平位移及应力分布

Figure  6.  Horizontal displacement and stress contours of surrounding soils under varying settlements

下载: 全尺寸图片 幻灯片

图7为数值模型在不同桩顶沉降时对应的桩周土体竖直位移及应力场分布。随着桩顶荷载的增加及桩体的压缩变形,桩周土体在靠近桩体压缩段也出现了较大向下的竖直位移,但位移量明显小于碎石颗粒的沉降量。因此可以判断在桩体压缩强烈的2.5倍桩径范围内,桩体和桩周土体产生了较大相对位移。图7（b）显示在该区域内桩周土体产生了较大竖直应力增量,增加的竖直应力主要由桩侧摩阻力及变形后的桩体与土体挤压产生,并通过土体剪应力以一定扩散角度传递到远离桩体的土体区域。在加载最后阶段,桩体的鼓胀变形也引起了靠近地表的部分桩周土体出现隆起和抬升。

图  7  不同沉降下土体竖直位移及应力分布

Figure  7.  Vertical displacement and stress contours of surrounding soils under varying settlements

下载: 全尺寸图片 幻灯片

图8为不同沉降下桩周土体剪切应变分布,土体剪切应变大小直接反映桩体变形对土体的扰动强弱。最先受到桩体扰动的为紧邻桩体压缩剧烈区域附近的桩周土体,扰动原因为桩土界面的剪切变形以及碎石颗粒刺入土体产生的应力集中。随着桩体沉降量的增加,土体剪切扰动区域在桩体水平鼓胀推动下向四周扩展。由于模型箱边界尺寸的约束,最终没有观察到贯通地表的联系剪切滑动面,但根据剪切应力发展的趋势在图中标出了潜在的剪切滑移面。

图  8  不同沉降下桩周土体剪应变分布

Figure  8.  Shear strain contours of surrounding soil under varying settlements

下载: 全尺寸图片 幻灯片

4.   碎石桩单桩承载过程中桩体受力变形特征

通过在DEM桩体内设置应力监测球（图9（a））,将球体范围内碎石接触力转换成平均应力分布曲线示于图9,可大致反映桩体宏观应力随着沉降的发展过程。图9（b）所示桩体内水平应力表现出随深度衰减的趋势,这与桩周土体水平应力的分布模式比较一致,颗粒接触力与单元应力取得了较好的一致性,说明离散单元颗粒与有限差分法网格耦合效果较好。但是桩顶处桩体水平应力与相应位置土体水平应力差距较大（图5（b））,桩体内应力显著高于桩侧土体应力。其原因为桩顶碎石颗粒受到加载板的强烈摩擦约束,因此尽管靠近地表处土体侧向约束力不强,但也不会发生鼓胀破坏,鼓胀破坏区域通常位于一倍桩径以下。

图  9  不同沉降下桩土平均应力分布

Figure  9.  Average stresses of stone column under varying settlements

下载: 全尺寸图片 幻灯片

图9（c）所示的桩体竖直应力也随埋深发生一定衰减,其桩端阻力显著低于桩顶应力。说明了作用在桩顶的一部分外荷载在桩体上部变形剧烈的区域通过水平力和桩侧摩阻力的形式传递到桩侧土体中。碎石桩由于鼓胀区域变形剧烈,即便在单桩荷载条件下也将与桩周土体发生较强相互作用,这一点有利于充分调用土体自身承载力,符合复合地基的设计原理。但应该注意,图9中桩体应力由一定数量颗粒接触力换算而成,其代表颗粒集合体的平均受力状态,与常规意义上土体中一点的应力状态概念需要有所区别。

通过离散单元法构建的碎石桩体模型,可以从细观尺度直接观察碎石颗粒的受力及位移。图10显示了桩体碎石颗粒在不同沉降下接触力的分布,接触力在图中显示为圆柱体,圆柱体的粗细与接触力大小相联系。由于桩体内包含了近万个碎石颗粒,颗粒接触总数更是数万个的量级,因此在图10（a）中仅显示50 N以上的较强接触力便于观察。这些接触力彼此联结并交叉分散形成了桩体内部独特的力链网络结构。可以看到在桩体顶部由于直接承受竖向荷载产生了数量最少但是强度最高的接触力链,这些强接触力链以一种树根状的形式向桩体下部生长并传递荷载。在外荷载向下的传递过程中,这些强力链逐渐转化成数量更多但是强度更低的弱力链。随着桩体沉降的增加,由强弱力链相互联结形成的树根状力链网络也逐渐向桩体底部发展。这与通过监测球提取的桩体内上大下小的桩体平均应力分布模式一致（图9）。

图  10  桩体内部力链网络

Figure  10.  Force chain networks inside stone column

下载: 全尺寸图片 幻灯片

图10（a）中显示的力链网络在桩体底部消失,实际上是因为底部力链网络更多由弱接触力组成因而没有显示。将桩体同一部位放大后按不同接触力区间分别显示在图10（b）中,可以观察到力链网络结构具有明显分级现象。其中大于200 N以上的接触力数量很少（图中区域115个）,但形成了力链网络中的主干,主干强力链方向以竖直为主。接触力大小位于100～200 N的中等接触力数量显著增加（达到365个）,形成了力链网络中的次主干,次主干力链的方向出现了一定偏转但大致还是沿竖直方向。小于100 N的弱接触力数量庞大（5049个）,与强接触力的个数出现了数量级的差异,形成了力链网络中的“毛细血管”,并且这些弱力链大部分与竖直方向呈较大夹角,一部分甚至出现近似水平方向。

力链网络结构的分布与形态揭示了碎石桩体承载和传力的机制：桩顶竖直荷载通过较强的力链主干、次主干传递到桩体内部；而数量众多的弱力链则将水平接触力比较均匀地传递给桩周土体,从而支撑其强力链主干。当桩周土体无法提供足够水平约束力时,鼓胀区域碎石颗粒的水平位移将无法限制,碎石颗粒弱力链无法支撑持续增大的强力链主干,从而导致主干力链网络的坍塌和桩体的宏观破坏。

与强力链分布区域一致,图11桩体内碎石颗粒的位移矢量在整个加载过程中均集中在桩体上部两倍桩径范围内。在加载过程后段可以观察到颗粒位移矢量从竖直方向往水平方向出现一定偏转,表示桩体上部出现了较大鼓胀（侧向）变形。破坏后的桩体轮廓切面与模型试验观察到的桩体变形形态也呈现较好的一致性（图12）。

图  11  不同沉降下碎石位移矢量

Figure  11.  Displacement vectors of gravels under varying settlements

下载: 全尺寸图片 幻灯片

图  12  最终破坏桩体轮廓

Figure  12.  Failure stone column at final loading stage

下载: 全尺寸图片 幻灯片

5.   结语

本文在考虑碎石颗粒三维形态的基础上,将碎石桩视为刚性凸多面体集合,桩周软土视为理想弹塑性材料,采用连续-非连续介质耦合数值模型对碎石桩单桩荷载模型试验的变形及破坏全过程进行了模拟。数值模拟获得的单桩荷载-沉降曲线以及鼓胀变形曲线与模型试验结果显示出较好的吻合程度,验证了本文数值模型的合理性。

通过荷载-沉降曲线、土体应力应变场以及桩体碎石接触力链网络讨论了碎石桩单桩承载破坏过程,并有效揭示了碎石桩宏观变形失效与其细观颗粒受力运动的联系。碎石颗粒相互接触形成的力链网络结构分布与形态能解释碎石桩体承载和传力的机制：荷载通过碎石桩内的强力链主干、次主干传递到桩体内部,数量众多的弱力链则将水平接触力均匀地传递给桩周土体并支撑其强力链主干。当桩周土体无法提供足够水平约束力时,鼓胀区域碎石颗粒的水平位移将无法限制,颗粒系统弱力链无法支撑持续增大的强力链,从而导致主干力链的坍塌和桩体的宏观破坏。

本文建立的数值模型能较好地模拟碎石的非连续变形过程以及由此引起的桩土相互作用,可进一步研究更多形式碎石桩承载特性,比如垫层效应、群桩效应及加筋碎石桩等。