近年来，山区机场工程越来越广泛，这些工程通常采用以碎石料为主体的土石混合料填筑而成的高填方地基。强夯法以其经济高效、施工简单、处理效果明显等优点，常被用来处理这类高填方地基，从而提高其承载力、降低其压缩性^[1]。

强夯作为强非线性动力过程，引起的地基动力响应非常复杂，很难用解析方法进行分析，故目前国内外学者主要采用数值模拟和现场测试等方法对其进行研究^[2-4]。常用的有限元等数值方法用于研究高填方体这类典型的散体介质时具有较大的局限性，无法反映散体介质特有的细观组构变化及其与宏观力学性质的关系。而现场测试一般费用较高，能获取的可靠数据有限，尚无法基于这些数据对强夯动力响应问题进行深入分析。基于离散颗粒介质特征建立的颗粒流理论^[5]，可以有效模拟动应力作用下的颗粒流动和大变形问题，亦可以实时观察颗粒的细观组构特征及其变化，为模拟强夯过程并探索填筑体的细观行为提供了一个具有一定优势的工具。利用其进行土体动力性质研究的可行性已为众多学者所验证^[6-8]。贾敏才等^[6]利用PFC^2D基于圆盘颗粒建立了模拟砂土地基强夯加固过程的细观颗粒流模型并与室内缩尺实验结果进行了对比，验证了颗粒流方法模拟强夯过程的可行性。蒋鹏等^[7]同样利用PFC^2D基于圆盘颗粒对块石土强夯过程进行模拟。Ma等^[8]采用PFC^2D/PFC^3D对块石土的强夯过程进行了模拟分析，模型利用强夯过程中地基土孔隙率的变化作为评价强夯效果和影响深度的指标，对夯锤形状、颗粒接触模型参数、颗粒级配对强夯效果及其深度的影响进行了细致分析。目前绝大多数模拟都是基于圆盘或圆球颗粒，没有考虑非圆颗粒形状和细颗粒的影响。

高填方填料中的碎石料往往由爆破开采得到，具有形状不规则的特点。在Deluzarche等^[9]和Alaei等^[10]采用PFC建立的碎石料模型中，均考虑了长条形颗粒。Lin等^[11]专门开发了椭球形颗粒用于离散元计算。王磊等^[12]利用PFC^2D采用clump方法近似椭圆形颗粒对强夯作用下碎石料颗粒的定向排列问题进行了分析。另外，山区机场高填方地基填料通常由碎石料和土料混合组成。已有研究表明：土石混合料中土料的加入使其压缩特性与纯碎石料有较大区别。它们在强夯响应和效果方面的区别则有待进一步研究。

本文针对上述高填方填料的特点，采用PFC^3D建立了考虑长条形碎石颗粒形状的土石混合料模型进行强夯模拟，并对强夯后土体的横观各向同性性质进行宏观和细观研究。为了研究颗粒形状和细粒含量对强夯效果的影响，同时建立了椭球形颗粒模型和球形颗粒模型作为对比模型。

1 颗粒流计算模型 1.1 颗粒接触模型

采用滞回阻尼模型来模拟碎石土和砂土在强夯作用下的动力滞回响应^[6-8]。如图 1所示，滞回阻尼模型其切向接触采用一般的线性接触模型，其法向接触引入加载模量和卸载模量，用两者模量的区别来反映颗粒之间的能量耗散。定义R_h为法向接触刚度k_{n, load}和法向卸载刚度k_{n, unload}之比，0 < R_h < 1。k_{n, load}和k_{n, unload}按下式计算：

$ {k_{{\rm{n,load}}}} = \frac{{2{R_{\rm{h}}}{k_{\rm{m}}}}}{{1 + {R_{\rm{h}}}}},\;\;\;{k_{{\rm{n,unload}}}} = \frac{{2{k_{\rm{m}}}}}{{1 + {R_{\rm{h}}}}}. $

(1)

其中，k_m为k_{n, load}和k_{n, unload}的平均值，当R_h=1时此模型退化为一般的线性接触模型。对于碎石料这种无黏性土，上述滞回效应只在受压加卸载时存在。根据Ma等^[8]的建议，取R_h=0.5。颗粒单元的法向刚度k_m、切向刚度k_s取为5×10⁷ N/m，摩擦系数取为0.3。

颗粒与墙之间采用线性接触模型，法向刚度、切向刚度与颗粒之间接触刚度取相同值。

1.2 地基模型

考虑到一般工程强夯的有效加固深度及模型计算量，取8 m×8 m×8 m范围的地基模型。模型的边界采用墙单元，在底面和侧面设置5个墙单元来约束颗粒法向位移。模型生成过程如下。

1)采用膨胀法在封闭墙单元内生成相应的球形颗粒，粒径分别为0.40、0.35、0.30和0.13 m，各占总体积的50%、25%、21%、4%，颗粒密度ρ=2 600 kg/m³。

2)根据“体积等效”和“质量等效”的原则，对0.40、0.35和0.30 m球形颗粒进行替换，生成长轴方向随机分布的椭球形颗粒(clump)来模拟长条形碎石颗粒，形状和尺寸如图 2a所示。若替换前的球形颗粒半径为R、密度为ρ，则替换后小球半径r₁、大球半径r₂、和密度ρ′分别为

$ {r_1} = \sqrt[3]{{\frac{{162}}{{569}}}}R,\;\;\;{r_2} = \frac{4}{3}\sqrt[3]{{\frac{{162}}{{569}}}}R,\;\;\rho ' = \frac{{569}}{{708}}\rho . $

(2)

3)施加重力计算至平衡。图 2b为三维填筑体模型示意图。本文建立了2个对比模型，球形碎石料地基模型和纯长条形碎石料地基模型。

1.3 夯锤模型

如图 2c所示，根据实际常见的工程参数^[2]，圆柱形夯锤直径取为3 m，锤重45 t，落距取为10 m (接触地基时夯锤初速度为14.14 m/s)，夯锤在模型中心处上方落下。采用57个直径为1.0 m的球形颗粒粘结成的颗粒簇(clump)来模拟夯锤。夯锤颗粒的切向刚度、法向刚度取为颗粒刚度的10倍，即5×10⁸N/m。

1.4 强夯过程模拟

模拟强夯过程时，首先在地基上表面生成夯锤，并赋予相当于从高处自由落体运动至表面的初速度，然后进行动力计算。待夯锤速度降为0后删除夯锤，再让体系在自重作用下计算至平衡。至此，第一次夯击模拟完成。此后重复这一过程，模拟连续夯击过程。本文模型暂不考虑强夯下颗粒的破碎效应。

2 强夯数值模拟结果分析

根据强夯数值模拟结果，首先分析重锤和地基动力响应，证明模型计算的可靠性。随后分析强夯过程中夯沉量和孔隙比的变化，采用孔隙比变化量表征强夯效果。

2.1 动力响应分析

图 3给出了重锤位移随接触时间的变化曲线，可以看出重锤约在80 ms时位移不再增加，重锤位移达到1.8 m左右。

图 4为锤下3、4、5和6 m的竖向动应力时程曲线，可以看出: 1)夯锤激发的地基内部动应力有一尖峰，同时伴随一个较弱的第二峰值，第一尖峰是由于夯锤下土体在压缩波作用下的推拉运动引起的。需要指出的是，第二峰值是边界的反射作用引起的，这一反射现象在实际工程中并不存在，而是由于计算量的限制导致本文PFC模型所选取范围不够大，没有完全消除这种反射现象。后续研究可以考虑利用对称性取1/4模型进行分析，从而可以在相同的计算量下将地基模型取更大一些，并且侧壁和底板边界可考虑加入传输边界，通过这2种方法尽量消除计算中反射波的影响。2)锤下3 m处竖向动应力峰值为726 kPa，约为重锤引起的自重应力的12倍，动接触应力的增大和衰减过程大约在60~80 ms内完成，但动接触应力增大和衰减过程并不对称，其增大过程较衰减过程快，这与文[4, 6-7]的实测与计算结果基本一致。3)竖向动应力峰值沿深度方向存在明显的滞后，且峰值随深度减小，这是压缩波在地基内传播并不断产生能量耗散导致的。

图 5为夯锤下正中、左侧1.5 m和左侧3 m不同深度处峰值应力的对比，可以看出在夯锤正下方，动应力峰值明显大于同一深度其他位置处动应力峰值，且随深度增加，动应力峰值逐渐降低。这说明强夯引起的土体内的压缩波能量在振源正下方这个方向最为集中，引起的动应力峰值最大，且随着传播距离的增加能量有衰减。在夯锤左侧1.5 m处，动应力峰值也呈现出随深度增加逐渐降低的规律，这也是压缩波能量随着传播距离的增加出现衰减导致的。在夯锤左侧3 m处，动应力峰值呈现出随深度增加逐渐增加的规律，这说明由于压缩波是沿一扩散角向四周传播，该处基本未受到强夯的影响，特别是深度较小的地方。动应力峰值的规律实际上是自重应力的规律。

2.2 夯击效果分析

为研究碎石颗粒形状及细颗粒含量对强夯效果的影响，提取土石混合料模型及2个对比模型(椭球形颗粒模型与球形颗粒模型)的计算结果进行比较分析。为方便行文，将球形颗粒模型、椭球形颗粒模型与土石混合料模型分别称为模型1、模型2、模型3。

图 6a为3个模型夯沉量随夯击数的变化情况。从图中可发现，随夯击数增加，总夯沉量增加，但每击增加的夯沉量逐渐减小，这是由于每一次夯击过后土体密实硬化、土体变形模量增加所造成的。

以夯锤正下方地基深度4 m处的孔隙比变化为例分析强夯引起的土体密实度变化情况。孔隙比的计算采用的是1 m直径范围的孔隙比平均值(由PFC^3D提供的测试球技术获得)。图 6b为3个模型孔隙比随夯击数的变化情况。可以看出，随夯击数增加，总体来看孔隙比在逐次减小，说明在计算的夯击次数内，强夯密实效果显著。图中初始孔隙比表现为模型1>模型2>模型3，这是模型地基生成的原因：在模型2生成过程中，由于椭球颗粒替换圆球颗粒造成地基在自重下沉降有所增加，所以模型2初始孔隙比较模型1小，而模型3在模型2的基础上，有小粒径土颗粒的掺入，故初始孔隙比更小。

为了对比不同模型的强夯效果，参考Ma等^[8]研究建议，采用孔隙比的减少量(初始孔隙比减去当前孔隙比)来表征强夯对土体密实度的提高效果，如图 6c所示。

从图 6c中可以看出：1)采用椭球形颗粒的模型2和3比采用球形颗粒的模型1强夯效果明显更好，这说明考虑长条形颗粒，土体在相同的夯击能作用下更易密实；2)对比模型2和3，考虑土料的模型3相比不考虑土料的模型2孔隙比减少量更大，这是由于小粒径土颗粒的填充作用使得强夯时地基更容易变得密实。

由上述分析可知，长条形颗粒比球形颗粒的强夯密实效果明显，细粒含量的存在也有利于强夯的密实作用。

3 强夯后土体横观各向同性研究

目前有关天然沉积砂土的研究表明^[15]，天然沉积的砂土出现颗粒定向排列^[16]，横观各向同性的性质显著：其不同方向的强度差异明显，而变形的各向异性表现在垂直于沉积面方向的模量较大，而沉积面内的模量较小。强夯填筑形成的高填方地基是否存在类似性质有待研究。以下通过对地基土宏观取样进行不同方向的三轴剪切实验对强夯后土体的性质进行研究，并结合细观组构分析其细观原因。为研究颗粒形状与细粒的影响，对3个模型进行对比分析。

3.1 三轴剪切实验

从强夯前、后的土体中，以锤下3 m处为中心，沿水平和竖直方向分别取出直径4 m、高4 m的圆柱体，即该圆柱体轴心线方向与水平面夹角分别为0°和90°。随后将大主应力施加于轴心线方向，以模拟大主应力方向与各向同性面(水平面)夹角(大主应力方向角)为0°和90°的情况。加载过程首先对试样施加围压(模型1为300 kPa，模型2、3为100 kPa)，然后施加三轴剪切荷载。不同模型计算的主应力差与轴向应变的曲线如图 7所示。

图 7是采用模型1、2、3进行强夯前后2个不同方向试样的三轴剪切实验曲线，可见强夯前2个方向试样的实验曲线基本相同，表明强夯前地基土为各向同性材料；而强夯后2个方向的三轴剪切实验曲线由于地基土密实度增加，比强夯前在模量和强度方面都有较大提高。图 7a强夯后2个方向实验曲线基本相同，说明强夯后仍是各向同性材料。但图 7b和7c强夯后表现出横观各向同性性质。进一步观察，大主应力方向角为90°对应的曲线在模量和强度方面均高于方向角为0°所对应的曲线，这说明地基土在竖直方向(垂直于各向同性面)强度和模量高于水平方向(平行于各向同性面)。

上述结果表明，当填料颗粒各向尺寸基本相同可采用球形颗粒来模拟时，强夯前后填筑体都表现出各向同性性质。原因是强夯动力冲击下，颗粒之间在竖直方向挤密的同时，也会在水平方向互相挤密，整体来说，各个方向挤密的程度基本相同，最终仍表现为各向同性。而当填料颗粒以长条形为主可采用椭球形近似模拟时，在强夯过程中颗粒受到竖向冲击作用会定向排列，从而导致竖向挤密的程度更为明显，最终表现为横观各向同性。上述定性分析可以结合下文的细观组构分析进行验证。

3.2 细观组构分析

细观组构分析是通过土体颗粒之间细观关系来反映土体的性质^[15-18]。组构参量是对细观组构的一种度量，通常用组构参量的分布函数来反映其统计特征。以往学者多采用二维组构参量进行描述，本文采用PFC^3D进行计算，故需要采用三维组构参量进行描述。并用PFC^3D自带的FISH语言开发了以下对三维组构参数进行统计的程序模块。

配位数：指与某颗粒相接触的颗粒数目，是衡量颗粒材料密实程度的一个指标，实际采用平均配位数，通过统计一定范围内颗粒之间总的接触数然后除以范围内的颗粒总数得到。

颗粒定向：对于椭球形颗粒，指的是颗粒长轴方向的空间分布情况。通过统计一定范围内颗粒的长轴方向与竖轴的夹角γ，并按下式进行处理得到：

$ E\left( \gamma \right) = \frac{{\Delta M\left( \gamma \right)}}{{M\left[ {\sin \left( {\gamma + \Delta \gamma } \right) - \sin \left( \gamma \right)} \right]/2}}. $

(3)

其中：E(γ)为长轴方向分布函数；ΔM(γ)为长轴方向在γ至Δγ角度范围内的颗粒数目；M为总颗粒数。

粒间接触法向：指颗粒接触处颗粒外表面法线的方向。通过统计一个范围内颗粒接触法线方向与竖轴的夹角γ，并仿照式(2)处理得到。

粒间平均接触法向力：指接触法线方向相同的接触法向力的平均值。通过统计一定范围内颗粒接触法向力，然后将一个角度范围内的接触法向力取平均值得到。

3.2.1 配位数

图 8为夯锤下2 m处半径1 m的球形范围内颗粒平均配位数随夯击数的变化情况。从图中可以看出，随夯击数增加，配位数逐渐增加，但增幅逐渐减小，到后期甚至有所下降。这说明强夯冲击作用在前几次夯击中时颗粒挤密，互相之间接触增多，配位数提高。但随夯击次数增大，颗粒之间已经很密，动力冲击作用反而会使颗粒间接触有所减少，配位数降低。对比3个模型，模型1配位数最小，模型2最大，模型3介于两者之间，说明椭球形颗粒较球形颗粒接触更多。由此可见颗粒形状对平均配位数有显著影响。

3.2.2 颗粒定向

图 9a为模型2和模型3中椭球颗粒在强夯后的颗粒定向分布。由于强夯前颗粒的排列是随机分布的，故颗粒定向分布图及下文中接触法线方向分布图、平均法向接触力分布图均为圆形(均不给出)。可以看出强夯后颗粒长轴方向集中在与水平面夹角0°~30°，这说明强夯作用下，颗粒发生了较为明显的定向排列，颗粒长轴方向向水平方向倾斜。

3.2.3 接触法线方向分布

图 9b为3个模型在强夯后颗粒间接触法向分布情况。从图中可以看出，模型1在强夯后颗粒间接触法向分布图保持为圆形，而模型2、3则呈椭圆形，竖向接触多于水平向接触。

3.2.4 平均法向接触力分布

图 9c为3个模型在强夯后颗粒间平均法向接触力分布情况。从图中可以看出，模型1在强夯后颗粒间平均法向接触力分布图保持为圆形，而模型2、3则呈椭圆形，竖向接触力大于水平向接触力。

综合上述细观组构分析，由长条形颗粒组成的地基在强夯后呈现明显的横观各向同性性质，在细观角度主要是由于在强夯过程中发生了颗粒的定向排列，其颗粒长轴方向向水平方向倾斜，竖向接触多于水平向接触，平均竖向接触力大于水平向接触力。而由球形颗粒组成的地基则无此现象，因此保持各向同性性质。

4 结论

本文针对以长条形碎石为主的高填方土石混合填料，采用PFC^3D建立了考虑长条形碎石颗粒形状的土石混合料模型进行强夯模拟，并对夯后地基的横观各向同性性质进行研究，得到如下结论：

1)模型计算的动力反应与实测结果规律相同，说明本文建立的模型在模拟强夯过程方面的可靠性。采用椭球形颗粒的模型2和3比采用球形颗粒的模型1强夯效果更好，考虑土料的模型3相比不考虑土料的模型2孔隙比减少量更大。说明当高填方地基填料为土石混合料且碎石形状以长条形为主时，强夯效果相比非长条形纯碎石料和长条形纯碎石料更明显。

2)通过对强夯前后地基土取样进行不同大主应力方向的三轴剪切实验，从宏观角度证明了当颗粒呈长条形时，土体经强夯后呈现明显的横观各向同性性质，竖直方向较水平方向模量及峰值强度明显增大，而当颗粒近似圆形时，土体保持各向同性。

3)采用三维组构参量对强夯前后地基土的细观组构进行了对比分析，从细观角度说明当高填方地基以长条形颗粒为主时，在强夯过程中颗粒受到竖向冲击作用进行了充分的定向排列，从而导致竖向挤密的程度更为明显，竖向接触数量和接触力均较大，最终表现为横观各向同性。

4)实际高填方地基采用的是土石混合料，当碎石形状以长条形为主时，更接近于本文模型3，故可认为这种高填方填料在强夯后存在横观各向同性性质。