碎石桩单桩受荷模型试验的离散单元法数值模拟
材料，采用二维离散单元与有限差分耦合数值方法建立了软土地基中碎石桩单桩竖向受荷模型，对碎石桩单桩受荷变形破坏及桩土相互作用全过程进行了模拟，通过荷载-沉降曲线、桩体及土体变形场应力场讨论了碎石桩单桩承载破坏机制.数值模拟结果与室内模型试验实测的荷载-沉降曲线、桩体鼓胀变形吻合良好，验证了本文数值模型的合理性.利用离散单元法建立的碎石桩模型无需复杂的本构模型假设便能较好地反映其鼓胀变形和失稳特性，同时讨论了模型的不足与有待进一步研究的問题.
关键词：复合地基;模型试验;离散单元法;数值模拟;鼓胀
中图分类号：TU44 文献标志码：A
Abstract：A 2D DEM-FDM coupled numerical modelling scheme was proposed to simulate a single stone column in soft clay. The discrete granular materials which firm the stone column can be better simulated by the DEM model， which is beyond the capability of continuum approaches. In the proposed numerical model， the stone column was represented by discrete rigid blocks based on the Voronoi tessellation， the surrounding clay was represented by the continuum Mohr-coulomb material. The settlement，bugling deformation，failure process and the involved complicated interactions between the stone column and surrounding soil were well captured by the proposed numerical model. The accuracy of the proposed model was verified using experimental model test results. The consistency between the numerical simulation and model test results shows a good potential for the proposed modelling scheme to further study the load-deformation and failure behavior of soft soils reinforced by stone columns. The limitations of the proposed model were also discussed in this study.
Key words：composite foundation;model test;discrete element method;numerical simulation;bulging
采用碎石桩对软弱土地基进行加固形成复合地基，能够有效提高地基承载力、减小沉降变形，并且具有施工快速简便、造价经济等优点.大量工程实践应用证明碎石桩复合地基是一种非常适合我国国情的地基处理方法.碎石桩复合地基的承载力取决于桩体置换率和布置，以及碎石桩桩体自身的承载力.碎石桩与其桩周土体形成复合地基，协调变形、共同承载，桩土的复杂相互作用决定了碎石桩承载破坏的机制，国内外许多研究者都进行过较为深入的研究[1-5].由于碎石桩体内部填料之间并无粘结强度，离散碎石的承载力来源于碎石颗粒之间的摩擦效应，所以桩体承载力受到桩周土的侧向约束力的影响.大量现场及室内试验证明当碎石桩埋置深度超过一定范围时（通常认为是4倍桩径）时，绝大多数碎石桩破坏是由桩体径向鼓胀变形引起的[6].故现在散体材料桩单桩极限承载力的理论及计算方法，均建立在碎石桩桩体上部发生鼓胀从而破坏失效的机制上[6-9].
室内及现场试验往往难以实时监测到碎石桩体与桩周土体共同受力变形的关系，很多学者因此采用数值模型对碎石桩复合地基变形及承载特性进行了研究[10-12].大部分数值模型基于连续介质方法，难以准确描述碎石的细观接触关系以及变形不连续性导致的桩土复杂相互作用.而离散元法则以颗粒或块体为单元，单元之间没有变形协调的约束，能很好地模拟碎石及土体大变形和位移非连续的过程[13].本文将碎石桩视为离散块体集合，桩周软土视为理想弹塑性材料，采用离散元数值模型对碎石桩单桩受荷室内模型试验的变形及破坏全过程进行模拟，研究碎石桩单桩受荷作用下承载及失效机制.
1 室内模型试验
如图1所示，在湖南大学地基试验室进行了多组碎石桩单桩受荷的模型试验[14].试验所填筑软土采用湘江枯水期的淤泥质软土.其原状土的不排水抗剪强度cu ≤ 10 kPa，天然含水率在30 % ～ 60 %之间，局部可达70 %以上.
模型试验地槽软土地基填筑深度为1 000 mm，填筑软土时分层填筑，每层200 mm高软土，共5次填筑完毕.填筑过程中，每层填筑完毕后采用5 kg的砝码从50 cm高度自由落下对软土层进行压实挤密，然后采用微型触探仪检验软基表面各点处的密实度和均匀性.模型地槽软土全部填筑完毕后静置较长时间，采用塑料薄膜覆盖包裹软土表面，防止水分蒸发，尽量满足正常固结条件.
2 碎石桩的离散单元数值模型
2.1 碎石离散单元数值模型
在离散单元法数值模型中，碎石等离散颗粒集合体可以假设为一系列刚性或弹性块体球体的集合.集合体的受力、运动及变形等一系列静力或动力过程则通过块体或球体之间的力学接触模型来模拟.二维离散块体力学问题的离散单元法程序UDEC可以用来模拟多边形块体集合体的力学行为.UDEC程序的基本组成对象是块体和接触，分别用刚体或者弹塑性力学关系和接触关系加以描述.每个块体的边界均设置成为接触，在力学求解过程中块体视为独立的运动对象，运动过程中通过接触与其他块体发生相互作用.单个块体又可进一步划分成有限差分网格，并赋予某种特定的连续介质本构模型，或者考虑为不可变形的刚性块体（图 2）.在UDEC中，碎石桩桩体可以通过大量离散凸多边形块体集合进行模拟，块体之间接触力及位移通过对相邻块体运动方程求解获得，而凸多边形块体在其接触部位能够相互滑动或脱离，从而模拟碎石桩鼓胀变形特性[15-16].块体与变形网格耦合计算的方法则是将块体与有限差分网格交界面上的力或应力与位移进行交互迭代，并在同一个时步内将迭代结果返回各自模型根据各自算法进行下一时步受力与变形计算.有限差分网格的节点力作用在刚性块体边界上，在一个时步内块体通过牛顿第二定律计算位移.块体角点位移则作用在有限差分网格节点上，作为节点位移通过本构关系计算有限差分网格内应力应变，下一个时步如此循环.
2.2 碎石桩单桩荷载数值模型
碎石桩桩体变形特性是通过块体接触的接触刚度（法向刚度kn及切向刚度ks）进行设置.桩体强度特性则通过接触的强度参数（摩擦角？准c、黏聚力cc及抗拉强度σtc）设置.其中碎石的摩擦强度一部分来源于块体摩擦角，一部分则来源于多边形块体之间的咬合摩擦.当块体接触内力未超过接触强度时，模型通过块体弹性变形模量和接触刚度产生弹性变形;而当接触内力超过其强度时，接触便发生破坏丧失其抗拉强度及黏聚力，此时破坏接触相邻的块体则发生滑移或者脱离，块体接触摩擦和块体咬合摩擦仍能提供一定抗剪强度[17].但不同于连续介质方法，离散单元法采用的是块体及接触的微观力学参数.室内试验通常无法直接提供这些微观参数，所以必须经过多次标定校核后才能获得相对应的微观参数.通过离散单元法模拟碎石桩体时，可首先根据碎石桩体的宏观力学特性比如压缩模量、摩擦角等预先选择微观参数，再进行多组数值试验（比如单轴、三轴、侧限压缩试验等）将模拟结果与试验获得宏观参数进行比对.根据对比修改微观参数，使数值模拟结果与实验室结果趋于一致.
本文碎石桩桩体模型基于UDEC程序建立，选取了与模型试验相近的Voronoi块体尺寸，接触的抗拉强度及黏聚力均默认为零.接觸摩擦角的取值可根据试验获得的碎石宏观摩擦角减去由于数值块体镶嵌咬合引起的摩擦角.块体咬合摩擦强度则可通过设置零接触摩擦强度的碎石块体数值试件进行三轴试验获得.由于多边形块体咬合摩擦角较大，所以本文模型接触摩擦角取了较小值.通过标定后的微观接触参数及对应的宏观力学特性见表 1.
模型中软土则采用摩尔库伦模型进行模拟，力学参数均从室内试验获得，列于表 2.
数值模型如图3所示，碎石桩首先通过Voronoi块体随机离散方法利用随机分布凸多面体块体集合进行模拟.模型简化为轴对称平面问题，模型边界均为简支边界.数值模型加载前在重力作用下达到其初始平衡状态.可以看到由于桩体土体刚度密度差异，碎石桩的设置给初始应力场造成了轻微扰动.在达到初始平衡后，桩顶荷载的施加与室内模型试验一致，通过桩顶刚性加载板施加持续位移产生.
3 碎石桩单桩荷载下变形失稳全过程数值模拟
图4为单桩荷载试验荷载-沉降曲线的数值模拟结果与模型试验结果对比，两者在沉降达到40 mm之前比较吻合.在沉降达到15 mm之前曲线初始段基本呈线性;随着荷载增大，沉降增长率逐渐增大，当位移达到20 mm时，非线性沉降明显发展;当沉降达到40 mm时，数值模拟的曲线出现了应力跌落，预示着数值碎石桩模型发生了垮塌失稳.此时数值模型中碎石桩内个别块体由于产生了较大接触变形，并在接近破坏状态时由于周围接触约束力在有限时步内无法抵消所累计的动能，所以发生个别关键块体离开碎石桩体，嵌入土体当中，从而导致数值桩体的垮塌，桩顶应力跌落.数值荷载-沉降曲线峰值因此可对应数值模型碎石桩的单桩极限承载力（89 kPa），与模型试验施加的最大荷载（91 kPa）非常接近.
图5为数值模型桩身侧限鼓胀变形与模型试验实测值的对比.桩体的径向应变可以表示其鼓胀变形.随着沉降的增加，碎石桩的鼓胀变形量逐渐增加，桩体不同深度的鼓胀量均有增加，沉降20 mm之前，模型试验实测鼓胀量在桩体上部较数值模拟结果偏大.
由于数值碎石桩采用紧密镶嵌的多边形块体，在荷载初始阶段产生了过高的咬合效应.当沉降达到20 mm以后，数值桩体的最大鼓胀变形与实测结果逐渐接近.沉降达到40 mm，数值桩体径向应变为7 %，模型试验实测值为10 %;沉降超过40 mm数值桩体破坏后径向应变为16 %，模型试验最后一级荷载阶段实测值为14 %.桩体的下部非鼓胀段径向变形数值与实测结果非常接近.与连续介质数值模型或者理论方法结果不同，离散单元法数值桩体径向应变沿桩长分布曲线很不光滑，这是由于每个凸多边形块体的运动并不连续，块体角点与土体网格接触处产生较大应力，而块体边与土体网格接触处应力较小，从而引起桩侧土压力分布具有一定离散性（图6）.
桩体上部鼓胀较大的部位对应着较大的侧向水平土压力.沉降达到40 mm时，鼓胀段承受最大侧向水平土压力为35 kPa.当桩体破坏后，由于块体嵌入桩周土体程度不一致，导致侧限水平土压力分布更不均匀.模型试验没有测得桩周水平土压力，但考虑到碎石颗粒的离散性，真实侧向水平土压力分布很有可能如数值结果一样不均匀分布，与通常假设为线性分布的被动土压力有较大差别.
图7、图8分别为数值模型不同沉降时对应的土体水平及竖直应力场分布.随着桩顶荷载的增加及桩身的变形，土体中的竖直应力随着桩土界面的剪应力传递也开始增加.即便桩土接触只设置了较小的剪切强度，但由于局部块体与土体的嵌入咬合使得桩体竖向荷载仍然能够通过剪应力传递到桩周土体当中.相对初始竖直应力场，沉降达到20 mm时模型底部竖直应力提高了一倍.在桩体宏观破坏之前，土体竖向应力是持续增加的.当沉降达到40 mm时，离散块体构成的桩体已经产生宏观破坏，导致出现了桩顶竖向应力跌落，此时下部桩体由于上部桩体竖向应力跌落实际处于卸载状态，从而导致底部竖向应力反而减小.在桩体宏观破坏之前，土体竖向应力是持续增加的.从土体水平应力场分布来看，更明显观察到最大土体水平应力发生在桩体碎石鼓胀最强烈的地方.水平应力因为直接受到碎石桩鼓胀变形影响，扰动较竖直应力场更剧烈，与初始应力场相比，鼓胀最剧烈处水平应力提高了将近7倍.模型所设置的软土刚度特别低，实际案例中水平应力提高幅度会更大.
图9显示了数值模型在不同沉降下的节点位移矢量分布，可以看到桩体位移随着深度逐步减小.在荷载初期，土体靠近桩体附近由于桩土剪切也会产生一定沉降;桩体破坏时，桩体鼓胀段周围土体产生较大水平位移，并可明显看到部分土体由于水平推移而产生的竖直上抬趋势.图10显示的是数值模型在不同沉降下桩周土体水平位移场分布，与鼓胀变形分布比较一致.在桩体破坏之前，两侧土体最大水平位移均大致发生在桩顶以下1.5倍桩径处;桩体破坏后，由于桩体内部出现了宏观剪切带，两侧土体最大位移位置不再一致.图 11为不同沉降下桩周土体剪切应变场分布，因为土体剪切应变大小直接揭示了桩土共同作用的强弱，从图中可以很明显观察到桩土作用强烈的区域集中在桩顶下2.5倍桩径以内的区域.
结合桩周土体在不同桩顶沉降下的塑性区发展（图 12）与图11所示土体剪切应变的发展，揭示了桩周土体的破坏过程.在荷载初期，由于桩身整体位移的差异，桩体上部发生了较大的压缩，从而导致桩体上部周围的土体发生较大剪切变形而产生塑性屈服;随着荷载增加，桩身整体都产生压缩，土体塑性区域往深部发展，并且由于桩身鼓胀变形的影响，部分桩侧土体塑性区域（图12中沉降20 mm时）往两侧开展;直到土体塑性区域扩展到土体表面或者土体边界，形成贯通的塑性破坏区域（被动土压力状态），桩体发生破坏丧失其承载力.最终破坏的模型见图12所示，此时桩体上部碎石块体间发生了较大的位移并导致部分碎石重新排列.由于侧向鼓胀变形，块体之间发生分离，产生了更多的空隙，从而减少了块体间的相互接触和咬合，与荷载初期相比桩体变成了一个更不稳定的离散系统.如果要继续维持这个离散系统的稳定则需要更大的桩周约束力，但此时土体已形成了整体破坏无法提供额外所需的约束力，桩体因此破壞并导致了图4所示的应力跌落.
4 结 语
本文将碎石桩视为离散凸多边形块体的集合，桩周软土视为理想弹塑性材料，采用二维离散单元与有限差分耦合数值方法建立了软土地基中碎石桩单桩竖向荷载模型，对软土地基中碎石桩单桩受荷变形破坏及桩土相互作用的全过程进行了数值模拟，通过荷载-沉降曲线、桩体及土体变形场应力场讨论了碎石桩单桩承载破坏机制.
碎石桩体单桩受荷主要破坏是由于鼓胀变形段碎石块体间发生了较大的相对位移.碎石块体在移动及重新排列过程中，发生镶嵌及分离，减少了块体间的相互接触和咬合，碎石桩体的不稳定性急剧提高.若桩周土体无法提供所需的额外约束力，桩体则将破坏并丧失其承载能力.
数值模拟结果和室内模型试验实测的荷载-沉降曲线、桩体鼓胀变形吻合良好，说明了本文数值模型的合理性.离散单元法建立的碎石桩模型无需复杂的本构模型假设便能较好地反应其鼓胀变形和失稳特性.与连续介质方法相比，本文建立的数值模型能较好地模拟碎石及土体大变形和位移非连续的过程，可进一步研究更多形式碎石桩承载特性，比如垫层效应、群桩效应及加筋碎石桩等.离散单元法细观接触参数的标定和选取方法则仍需进一步研究.
参考文献
[1] 陈建峰，韩杰. 夯扩碎石桩群桩承载性状研究[J]. 中国公路学报，2010，23（1）：26—31.
CHEN J F，HAN J. Research on bearing behavior of rammed aggregate pier group[J]. China Journal of Highway and Transport，2010，23（1）：26—31.（In Chinese）
[2] 赵明华，张玲，邹新军，等. 土工格室-碎石桩双向增强复合地基研究进展[J]. 中国公路学报，2009，22（1）：1—10.
ZHAO M H，ZHANG L，ZOU X J，et al. Research progress in two-direction reinforced composite foundation formed by geocell reinforced mattress and gravel piles [J]. China Journal of Highway and Transport，2009，22（1）：1—10. （In Chinese）
[3] 刘杰，张可能. 散体材料桩复合地基极限承载力计算[J]. 岩土力学，2002，23（2）：204—207.
LIU J，ZHANG K N. Calculation of ultimate bearing capacity of composite foundation on discrete material pile[J]. Rock and Soil Mechanics，2002，23（2）：204—207. （In Chinese）
[4] 龚晓南. 广义复合地基理论及工程应用[J]. 岩土工程学报，2007，29（1）：1—13.
GONG X N. Generalized composite foundation theory and engineering application[J]. Chinese Journal of Geotechnical Engineering，2007，29（1）：1—13. （In Chinese）
[5] 刘杰，赵明华，何杰. 碎石桩复合地基承载及变形性状研究[J].湖南大学学报（自然科学版），2007，34（5）：15—19.
LIU J，ZHAO M H，HE J. Research on bearing and deformation characters of composite foundation with granular columns [J]. Journal of Hunan University （Natural Sciences），2007，34（5）：15—19. （In Chinese）
[6] BRAUNS J. Die anfangstraglast von schottersaulen im bindigen untergrund[J]. Die Bautechnik，1978，55（8）：263—271.
BRAUNS J. The initial load of graval pile in the clay foundation[J]. The Construction Technology， 1978， 55（8）： 263—271. （In Germany）
[7] HUGHER J，WITHERS N J. Reinforcing of soft cohesive soils with stone columns[J]. Ground Engineering，1974，7（3）：42—49.
[8] 陈昌富，李其泽.基于圆孔扩张理论顶部加箍碎石桩承载力计算[J].湖南大学学报（自然科学版），2011，38（10）：7—12.
CHEN C F，LI Q Z. Calculation of bearing capacity of geosynthetic-encased gravel pile based on cavity expansion theory [J]. Journal of Hunan University （Natural Sciences），2011，38（10）：7—12. （In Chinese）
[9] WONG H Y. Field instrumentation of vibroflotation foundation[J]. Field Instrumentation in Geotechnical Engineering，1974，23（4）：475—487.
[10] 蒋敏敏，肖昭然，蔡正银. 高速公路碎石桩复合地基加固数值模拟[J]. 公路，2012，31（1）：8—12.
JIANG M M，XIAO Z R，CAI Z Y. Numerical simulation of stone columns reinforced foundation of expressway[J]. Highway，2012，31（1）：8—12. （In Chinese）
[11] 陈建峰，韩杰. 夯扩碎石桩单桩载荷试验数值模拟[J]. 岩土工程学报，2009，31（9）：1366—1370.
CHEN J F，HAN J. Numerical modeling of loading tests on a rammed aggregate pier[J]. Chinese Journal of Geotechnical Engineering，2009，31（9）：1366—1370. （In Chinese）
[12] 刘洋，闫鸿翔，汪成林. 复合振冲碎石桩加固机理及施工过程数值模拟[J]. 工程科学学报，2015，37（3）：390—397.
LIU Y，YAN H X，WANG C L. Reinforcement mechanism of vibro-stone columns with additional wicks and numerical analysis[J]. Chinese Journal of Engineering，2015，37（3）：390—397. （In Chinese）
[13] 蒋明镜，肖俞，陈双林，等. 砂土中单桩竖向抗压承载机制的离散元分析[J]. 岩土力学，2010，31（S2）：366—372.
JIANG M J，XIAO Y，CHEN S L，et al. Discrete element analysis of bearing mechanism of single pile in sand under vertical load[J]. Rock and Soil Mechanics，2010，31（S2）：366—372. （In Chinese）
[14] 赵明华，顾美湘，张玲，等. 竖向土工加筋体对碎石桩承载变形影响的模型试验研究[J]. 岩土工程学报，2014，36（9）：1587—1593.
ZHAO M H，GU M X，ZHANG L，et al. Model tests on influence of vertical geosynthetic-encasement on performance of stone columns[J]. Chinese Journal of Geotechnical Engineering，2014，36（9）：1587—1593. （In Chinese）
[15] HANY，CUNDALL P. Verification of two-dimensional LBM-DEM coupling approach and its application in modelng episodicsand production in borehole[J]. Petroleum，2016，3（2）：179—189.
[16] CUNDALL P A，HART R D. Numerical modeling of discontinue[C]//Comprehensive Rock Engineering. Oxford：Pergamon Press，1993：193—213.
[17] 譚鑫，HEINZ K. 含层理构造的非均质片麻岩巴西劈裂试验及离散单元法数值模拟研究[J]. 岩石力学与工程学报，2014，33（5）：938—946.
TAN X，HEINZ K. Brazilian split tests and numerical simulation by discrete element method for heterogeneous gneiss with bedding structure[J]. Chinese Journal of Rock Mechanics and Engineering，2014，33（5）：938—946. （In Chinese）
上一篇：论固流耦合渗透固结数学模型的建立