软土地基连续强夯置换碎石墩的数值分析①

谢新宇，徐玉胜，吴 健，3，江辉煌，吴勇华

（1.浙江大学 教育部软弱土与环境土工重点实验室，浙江 杭州 310027；2.铁道部科学研究院深圳研究设计院，广东 深圳 518034；3.中国汉嘉设计集团，浙江 杭州 310005）

0 引言

随着经济发展和城市化进程的深入，沿海发达地区城市用地紧张逐渐加剧，填海造陆工程不断出现。对于大面积填海工程，由于新近冲填土在变形、稳定、强度方面达不到设计要求，需对这类地基进行处理。强夯法由于效果显著、施工便捷、周期短等优点而被广泛使用。

强夯是土体短时间受到巨大冲击载荷作用的一种复杂非线性动态响应过程，具有非常明显的动力特性，而且碰撞区土体都迅速超越弹性阶段而进入塑性流动阶段，并有出现撕裂等各种形式的破坏。Chow［1］等分析了强夯地基土的应力和位移，但是由于未考虑侧向应力和应变而与实际情况存在较大出入；孔令伟等［2－3］在考虑夯锤自重的基础上得到了强夯边界接触应力与沉降在变换域中的解析式并确定了接触应力时间；张冬丽等［4］将有限元方法用于黄土地基震陷预测。任锡泰［5］对强夯的地震效应进行测试与分析；田水等［6］考虑瞬态大变形的非线性动态显式有限元，分析了强夯对地基土的冲击碰撞过程；蔡袁强［7］分析了夯击过程中土体的塑性应变，预测加固区空间分布；蒋鹏等［8］采用大变形几何非线性有限元法模拟了强夯对地基土的冲击碰撞过程，但在有限元求解中对非线性控制方程进行了线性拟合。在试验研究方面，王铁宏［9］首次采用16000 kN·m能级对某港口罐区的碎石回填场地进行了强夯试验研究，分析了锤重、落距等组合因素对有效加固深度等的影响。

碎石受到重锤冲击荷载进入土体，在碎石的挤压以及冲击波的作用下，土体短时间内吸收冲击能产生塑性应变，土体硬度也随之迅速减小，土体产生塑性破坏，失去“承担”荷载的能力；而其它较远的土体颗粒则受到冲击荷载产生的应力波作用，颗粒之间可能重新排列，土体被挤密压实。这是一个非常复杂的过程，特别是在连续荷载作用下不断受到冲击荷载作用的土体性状变化更加复杂。本文借助于有限元软件，采用土体损伤失效模型分析在连续荷载作用下每一击作用下土体失效范围及碎石墩的形成过程，同时还分析不同夯击能、夯锤直径和垫层厚度对强夯置换效果的影响。以期对以后的强夯工程在碎石工程量等方面可以提供指导。

1 有限元模型的建立

1.1 基本原理

土体在冲击荷载作用下部分区域吸收能量产生较大的塑性应变，甚至完全破坏，此时通常采用土体损伤模型，即认为在冲击波以及碎石直接作用下土体单元硬度（强度）逐渐衰减，失去硬度（强度）的土体单元因不能承担荷载而将从原来位置移除，其基本原理如图1所示。

图1 土体损伤失效模型Fig.1 Failure model of soil mass.

图1中，AB段为线弹性阶段；BC段为弹塑性硬化阶段；从C点开始，虚线延续土体硬化模型，而实线CF则表示土体在冲击荷载作用下的失效过程。D值表征土体硬度（强度）失效率，在C点土体能承担最大荷载σy0，D＝0。F点土体完全失效，D＝1。可由公式（1）表示：

土体损伤参考Hillerborg［10］提出的单位破坏能Gf，如公式（2）所示。

D为表征土体硬度失效的变量，在土体损伤产生破坏过程中其值不断发生变化。其计算可以由式（3）表示。表征土体损伤过程有很多判断法则，各个法则都有自己的判断值。

式中：dj为土体损伤模型中各自判断指标；dmult为各个指标乘积对损伤判断的表征值：

ABAQUS软件提供了弹性、弹塑性、超弹性以及黏弹塑性等模型，本文土体与碎石层采用Drucker－Prager弹塑性模型；重锤采用弹性模型。本文采用总应力法进行研究，未考虑孔隙水压力的作用。

1.2 计算模型和单元类型

模型采用轴对称形式。由于拟静力分析法无法反映夯击过程中波的特性，本文采用动力显式分析方法对夯击过程进行模拟。单元类型采用3节点线性轴对称单元（CAX3）和4节点线性轴对称单元（CAX4）。在土体塑性较大产生破坏区域网格加密，而较远部分适当增大单元尺寸。重锤与地面接触瞬时的的速度可由下式计算得到：

重锤与碎石之间接触采用摩擦接触形式，其函数为罚函数。模型尺寸以及地质参数选取见表1。

表1 模型计算参数

1.3 模型网格与边界条件

强夯动力显式算法的精度对网格依赖性较大。如式（2）所示，土体耗散能、失效判别都与网格大小直接有关。Kuhlemeyer等［11］认为，采用有限差分法时单元节点间的距离Δx应小于波长λmax的1／8；Kausel［12］认为网格大小要满足Δx≤λmin／6；而Costantino等［13］把最小网格定义为Δx≤λmin／4.4。根据不同计算精度的要求，一般可选取

式中：λT为与所需考虑的截止频率所对应的波长；Δx为单元的长边尺寸。本算例中土体产生破坏以及附近区域网格最大尺寸取为0.005m，较远区域则以等比例方式逐渐增大。

强夯过程中能量以应力波的形式在土体中传播。有限元模型受尺寸限制，在模型边界上应力波将产生反射和折射等现象，反射应力波又会二次对破坏区域及附近土体产生作用。而实际地基是半无限空间，为了减小这些反射应力波的应力的影响，文中采用有限元与无限元相结合的方式来设置模型边界条件。

2 室内模型实验与有限元结果分析

室内模型试验模拟了连续强夯25次土体失效、碎石进入土体过程。每次夯击以后夯锤下部部分碎石进入土层。模型实验测试了每次夯击后所形成碎石墩的高度（图4）。

在数值模拟时，每一次夯击前网格都需要重新建立，前一次强夯形成的碎石墩将成为后一击的初始状态。表征土体失效的指标有SDEG参数、网格失效图，表征土体变形有竖向位移图。SDEG表征土体硬度（强度）衰减程度，本计算中假定土体硬度（强度）衰减到达最大值99%，即D≥0.99时土体完全失效，对应的网格将从原位移除，相当于实际中的土体被碎石侵入，原来土体不再发挥任何作用。为了节省篇幅，文中仅选取了前面两次夯击以后土体的SDEG云图以及竖向位移图。其它每次夯击过程基本一致。

图2 第一次夯击后SDEG云图、失效网格图和竖直方向位移云图Fig.2 The SDEG nephogram （a），failure mash（b）and displacement nephogram in vertical direction（c）of the soil samples after first dynamic compaction.

图2（a）、（b）分别为第一次夯击以后土体SDEG云图以及失效网格图。从图中可以看出，在冲击荷载作用下，由于碎石硬度较土体大，冲击能主要通过碎石传递给土体，碎石直接侵入土体，而土体则产生损伤破坏，部分吸收的能量大于塑性应变能的土体此时已经不能承担荷载。图2（b）中网格空白部分代表失效网格部分，在冲击荷载作用下，这一部分土体失效，网格单元从模型中移除。失效区域成锅底型，深度与上部夯击能、土体失效准则等有关，文中第一击后碎石侵入深度大概在0.8m左右。图2（c）为在竖直方向位移云图，重锤在重力作用下将碎石挤入土体，形成凹坑，在凹坑附近碎石层受到侧向荷载向上发生位移。模拟结果与实际情况也相符合。

图3 第二次夯击后的SDEG云图，失效网格图和垂直方向位移云图Fig.3 The SDEG nephogram （a），failure mash（b）and displacement nephogram in vertical direction（c）of the soil samples after second dynamic compaction.

图3（a）、（b）、（c）为第二次夯击后土体SDEG云图、失效网格图及竖直方向位移云图。第一次夯击后碎石墩形状（图2（b）），在第二次建模时适当做了简化，呈倒梯形。第二次夯击以后碎石墩附近土体吸收能量产生破坏，如图3（b），失效土体主要在碎石墩附近区域。

最后两击土体失效网格仅在碎石墩与碎石层接触处有很小一部分，可以判断在此夯击能作用下碎石墩基本成型，即为收锤标准。

模型箱夯击次数与每次夯击后碎石墩高度关系曲线如图4所示。

图4 夯击次数与碎石墩深度Fig.4 Compaction numbers with stone column depths.

3 主要影响因素分析

由上节的室内模型实验结果与数值模拟的对比，我们得到了数值模拟时主要参数的选择方法。为分析实际工程中不同工况下夯击次数与碎石墩深度之间的关系，本文分析了不同夯击能、垫层厚度、夯锤直径时的碎石墩深度随夯击次数的变化规律。

3.1 不同夯锤直径

不同夯锤直径下夯击次数与最大碎石墩高度关系曲线如图5所示。每次夯击能相同，均为8000 kN·m，夯锤重45t，其它土层参数如表1所示。本算例中土层一、土层二和碎石层厚度分别为10、10和5m。下面算例中土层一、二厚度均相同。

图5 夯锤直径与最大夯击次数关系曲线Fig.5 Relationship curve of hammer diameters with maximum hammer numbers.

通过计算发现，夯锤直径与高度比值较小时，如夯锤直径在1m时，第一次夯击后土体侧向损伤失效的比较多，而碎石墩向下发展的较少，现场实际数据也是相同现象；当夯锤直径为4m时，此时由于直径较大，夯锤与土体接触面积较大，强夯置换效果不明显；在直径为1.5m时可以夯击的次数最大，夯击深度也最大。

图6表示直径为1.5～4m时夯锤直径与碎石墩高度关系曲线。在直径为1m时，效果很不明显；直径为4m时碎石墩只有0.3m左右；在直径为1.5m时碎石墩高度最大，说明在此夯击能下夯锤最优直径为1.5m。

图6 夯锤直径与碎石墩高度关系曲线Fig.6 Curve of relationship between hammer diameters and stone column depths.

3.2 不同垫层厚度

在实际工程中垫层厚度对强夯置换效果影响较大，垫层过厚时夯锤夯击能基本被上层碎石层吸收，对下层需要处理的地基影响较小；垫层过薄时夯锤击穿垫层，也起不到强夯置换效果。因此一定夯击能下需要合适的垫层厚度。本算例主要分析垫层厚度对置换效果的影响。垫层厚度变化范围为3～11 m。垫层厚度与最大夯击次数关系如图7所示。

图7 垫层厚度与碎石墩高度关系曲线Fig.7 Curve of relationship between cushion thicknesses and stone column depths.

由图8可以看出垫层厚度增加碎石墩高度迅速减小。当然若增加夯击能，则效果也会有所改善。算例中垫层厚度在5m时效果较好，在3m时夯锤击穿垫层，这在实际过程中可能造成丢锤。

图8 夯击能量与碎石墩高度关系曲线Fig.8 Curve of relationship between compaction energy and stone column depths.

3.3 不同夯击能

实际工程中垫层较薄则不能用较大夯击能，需要先小能量夯击，待碎石墩形成一定高度再加大夯击能量，垫层较厚时则需要大能量夯击，不同夯击能下碎石墩发展过程如图8所示。

在垫层为5m情况下，夯击能增加则碎石墩高度也增加。在10000kN·m时最大；同时通过计算也发现，在10000kN·m时第一击夯坑深度较大，在实际中有可能导致夯锤丢失。故在数值分析中先用8000kN·m夯击，再增加到10000kN·m进行夯击。

3.4 不同位置土体应力变化

为分析比较土体内部不同位置在受到冲击荷载作用过程中应力变化分布规律，本文选取了土体内部几个代表性的点。图9为各点分布情况：A0～A9为距离轴线3m处同一水平位置但不同深度的点，A0为碎石层顶面点，A3、A5、A7、A9分别距离顶面3m、5m、7m和9m处的点；H0～H4分别为距离碎石层顶面8m处同一水平线上不同点，H0位于中轴线上，H1、H2、H4分别距离中轴线1m、2m、4 m处的点。

图9 各计算点分布图Fig.9 Distribution of calculation points.

图10为各点应力随时间发展曲线。可以看出：（1）每次夯击过程中，应力S22峰值到达的时间随深度增加而推后。如图10（a），总的夯击时间为0.2 s，但3m、5m、7m、9m 处峰值到达时间分别为0.02s、0.04s、0.07s和0.11s左右；（2）本算例中，应力S22在3m处峰值最大，即夯击能量在传递过程中因为土颗粒阻力而逐渐耗散。夯击加固土体有一定区域；（3）在同一个水平面上夯锤底部一定区域内，如图10（b）所示，不同点峰值发生的时间基本一致，但是峰值大小靠近中轴线处最大，距离中轴线越远值越小。

图10 各计算点S22应力与时间关系曲线Fig.10 Relation curves for stress S22and time at different calculation points.

4 结论

（1）土体单元采用损伤力学以及土体失效准则，可以模拟土体失效区域、碎石墩形成过程，对强夯碎石墩工程量预测具有重要的参考作用。

（2）一定夯击能下，碎石墩高度不会无限增大。即一定落距和锤重条件下，仅增加夯击次数并不经济。

（3）随着深度增加，应力峰值到达时间推后，应力最大值在一定条件下只能达到某个深度。即夯击加固存在一个有效深度。

（4）夯击能一定，存在一个最优夯锤直径，在此直径下碎石墩高度达到最大；夯击能一定，垫层厚度过大或过小置换效果都不好，需要确定一个最合理的垫层厚度；同样，夯击能增大碎石墩高度也增大，但是在垫层厚度一定时，夯击能过大可能导致丢锤事故发生，需要先小能量夯击几次再逐渐增大夯击能。

［1］Chow Y K，Yong D M.Dynamic compaction analysis［J］.Journal of Geotechnical Engineering，ASCE，1992，118（8）：1141－1157.

［2］孔令伟，袁建新.强夯的边界接触应力与沉降特性研究［J］.岩土工程学报，1998，20（2）：86－92.

［3］孔令伟，袁建新.强夯作用下成层地基的表面接触应力分布特征［J］.力学学报，1999，31（2）：250－256.

［4］张冬丽，王兰民，王玉华.有限元分析方法在黄土地基震陷预测中的应用［J］.西北地震学报，2002，24（4）：208－216.

［5］任锡泰.强夯施工地震效应测试与分析［J］.西北地震学报，2002，24（4）：352－355.

［6］田水，王钊.强夯加固机理的非线性动态冲击显式有限元数值分析［J］.工程地质学报，2006，14（05）：694－698.

［7］蔡袁强，陈超，徐长节.强夯加固回填土地基的三维数值模拟［J］.岩土力学，2007，28（6）：1108－1112.

［8］蒋鹏，李荣强，孔德坊.强夯振动影响的数值分析［J］.地下空间，2001，21（5）：544－548.

［9］王铁宏，水伟厚，王亚凌.高填方碎石土地基16000kN·m能级强夯试验研究［A］∥第十届土力学及岩土工程学术会议论文集（下）［C］.重庆：重庆大学出版社，2007，43－48.

［10］Hillerborg A.Results of three comparative test series for determining the fracture energy Gf of concrete［J］.Material and Structures，1985，18（107）：407－413.

［11］Kuhlemeyer R L，Lysmer J.Finite element method accuracy for wave propagation problems［J］.Journal of Soil Mechanics and Foundation Engineering，ASCE，1973，99（SM5）：421－427.

［12］Kausel E.Wave propagation in anisotropic layered media［J］.International Journal for Numerical Methods in Engineering，1986，23（8）：1567－1578.

［13］Costantino C J，Miller C A，et al.CARES：Computer Analysis for Rapid evaluation of structures（Version 1.2）［M］.Washington DC.US Nuclear Regulatory Commission，1995.