粗粒土大型静止侧压力系数测定试验的颗粒流模拟

蒋明杰，栗书亚，吉恩跃，张小勇，朱俊高

(1.广西大学工程防灾与结构安全教育部重点实验室，南宁 530004；2.广西大学广西防灾减灾与工程安全重点实验室，南宁 530004；3.南京水利科学研究院，南京 210024；4.河海大学岩土力学与堤坝工程教育部重点实验室，南京 210098)

粗粒土因具有抗剪强度高、压实性好的特点，常作为土石坝和道路路基的重要填筑材料。静止侧压力系数K0作为评价粗粒土基本力学性质的重要参数，研究K0对岩土工程有重要的理论意义。由于粗粒土的最大粒径较大，通过大型静止侧压力系数试验仪准确地测定其K0值耗时长、费用昂贵。因此，通过建立相关数值仿真模型，对颗粒细观参数的调整来探究K0演化规律，可规避试验中设备以及人为误差。

离散元软件PFC3D由于能从细观角度模拟土体力学特性，在土力学试验的模拟中已有一定的应用。例如，邢炜杰等[1]分析了细观参数的变化对黏性土的应力-应变关系曲线的影响，得出边界条件、孔隙率、黏结强度、摩擦因数等细观参数对黏土宏观力学性质的影响规律；陈亚东等[2]提出了砂土的孔隙率、内摩擦角及压缩模量等宏观力学参数相匹配的细观参数的确定方法，并给出了细观结构参数与宏观力学参数的函数关系；田世雄等[3]对不同粗粒含量砂卵石土宏观及细观力学特性开展了研究，提出了不同粗粒含量砂卵石土的接触模量、摩擦因数等颗粒离散元细观参数；杨舒涵等[4]采用离散单元法进行了颗粒材料的真三轴加载路径数值试验来研究粒间摩擦对颗粒材料宏观力学特性的影响；周玉县等[5]基于平行黏结模型，建立了岩土颗粒材料的摩擦因数、颗粒黏结强度等细观参数与宏观强度参数的标定关系。

虽然离散单元法已被广泛应用于土体力学特性的研究[6]，但极少从细观角度对粗粒土静止侧压力系数进行模拟分析。因此，现利用三维颗粒流分析程序(PFC3D)基于大型静止侧压力系数试验仪建立仿真模型，探讨孔隙率、颗粒摩擦因数、切向刚度、法向刚度等细观参数对粗粒土K0的影响规律，建立粗粒土K0与孔隙率和摩擦因数等参数的函数关系。

1 离散元模型建立及验证

1.1 离散元模型的建立

建立的K0试验数值模型的边界由一个刚性圆桶和两个作为加压板的刚性方形板组成，模型尺寸与大型K0试验仪的压力室[7]尺寸一致，直径400 mm，高300 mm，如图1所示。其中刚性圆筒与下部加压板在试验过程中保持固定，通过对上部加压板设定法向移动，对散粒体进行轴向加载，从而使散粒体处于单向压缩状态。

在K0试验模拟过程中，为了上部加压板产生稳定的压力μn，给模型边界(刚性墙)施加一定的速度来产生恒定的轴压，数学表达式为

μn=G(σm-σr)

(1)

(2)

式中：G为伺服调整系数；σm和σr分别为墙上的应力和目标应力；fwx和fwy分别为墙与颗粒在x、y方向的接触力；A为边界墙面积。

通过伺服机制，不断调整上压板的运动速度使墙体平均应力达到指定值并保持恒定，保证模型在加载过程中的稳定性。

通过FISH函数计算获得加载过程中试样的应力状态，表达式为

(3)

K0=σh/σv

(4)

式中：σv、σh分别为试样竖向应力以及水平应力；F1、F2、F3分别为刚性墙1、2及3所受力。由式(4)可得到实时土体静止侧压力系数。

1.2 细观参数标定

选取文献[8-9]试样D10～D60的K0试验数据作为依据，标定堆料K0数值模拟细观参数，级配参数如表1所示。模型所受最大竖向应力为4 MPa，颗粒密度取2 650 kg/m3，为准确模拟K0试验，结合相关参考文献[10-11]，颗粒数目超过30 000个。在细观参数标定过程中，通过反复调整细观参数(300～400组)，最终标定出堆石料的ks、kn、τc、σc及μ共5个细观参数，如表2所示。根据表2细观参数模拟出来的K0-σ′v关系曲线及试验所得K0-σ′v关系曲线一同绘于图2。

r为刚性墙1水平截面半径；h为试样实时高度图1 离散元模型示意图Fig.1 Schematic diagram of discrete element model

表1 试样级配参数Table 1 Grading parameters of samples

表2 数值模拟材料细观参数Table 2 Numerical simulation of material parameters

由图2可知，D10～D60数值试验K0-σ′v关系曲线趋势与室内试验结果较为接近，相同竖向应力下数值试验所得K0与室内试验结果相差不到14%，因此可以认为表2细观参数与实际情况接近。

2 离散元模拟结果分析

2.1 位移云图变化

模型加载过程中颗粒位移变化情况如图3所示，颗粒颜色由青到红表示位移逐渐增大。由图3可知，随着竖向应力的增加，颗粒受到上部压力逐步向下部移动，且随着压力的增加位移逐渐增大。该现象从细观角度解释了粗粒土体变形机理：在不考虑颗粒破碎的条件下，土体变形量主要是由楔入和咬合颗粒易位，颗粒骨架结构的改变组成。

2.2 颗粒间接触力分布

不同竖向应力对应的模型应力链分布，如图4所示。力链的粗细代表颗粒间接触力的大小，力链的个数代表了接触力的数目。随着竖向应力的增大，应力链逐渐扩展至密实，并逐步向试样下部发展，反映到宏观力学上是，粗粒土随着竖向应力的增加，由疏松的状态逐步向较为密实的状态转变。

2.3 孔隙率对粗粒土K0的影响

为分析孔隙率对粗粒土K0影响，设计初始孔隙率n为0.3、0.36、0.4和0.5共4种情况下的K0试验，采用D60试样级配以及表2中细观参数，并将各孔隙率试样的K0-σ′v曲线绘于图5。由图5可知，在相同竖向应力下，试样孔隙率越大，土体密实程度越高，其静止侧压力系数越小。另外，根据各试样数值试验结果整理得σ′v分别在500、1 000、1 500、2 000及3 000 kPa时，对应不同初始孔隙率试样的K0值及各σ′v相对应的孔隙率n，如图6所示。由图6可知，当试样所受竖向应力一样时，试样K0随着孔隙率n的增大逐渐增大，但增大趋势趋于平缓。

图2 数值试验以及室内实验K0-σ′v关系曲线对比Fig.2 Comparison of relationship K0-σ′v curves between numerical experiments and laboratory experiments

根据实验结果，发现K0与n的关系可表示为

(5)

式(5)中：A1和A2为材料参数。

将加载时不同竖向应力下拟合所得堆石料材料参数A1和A2以及决定系数R2分别列于表3，所得拟合曲线绘于图6。由图6可知，拟合曲线与对应试验数据吻合度较高，最大误差仅4.03%，决定系数R2最小值也达到了0.97，显然，式(5)拟合效果较好。因此，在相同σ′v下堆石料静止侧压力系数K0与孔隙率n关系可用式(5)表示。

图3 颗粒位移云图变化Fig.3 Changes in particle displacement nephogram

图4 颗粒间接触力随竖向力小应力变化Fig.4 Variation of contact force between particles with vertical stress

图5 不同孔隙率下各试样的K0-σ′v曲线Fig.5 K0-σ′v curves of specimens with different porosity

表3 各试样不同竖向应力下的A1和A2Table 3 A1 and A1 of each sample under different vertical pressure

图6 不同竖向应力下试样的K0-n曲线Fig.6 K0-n curves of each specimen under different vertical stress

2.4 摩擦因数对粗粒土K0的影响

为研究粗粒土表面粗糙程度对粗粒土静止侧压力因数的影响，本节设计了摩擦因数μ分别为0.4、0.5、0.6及0.74共4种情况下的K0试验，采用D60试样级配以及表2中细观参数，并将不同摩擦因数试样的K0-σ′v曲线绘于图7。由图7可知，在相同竖向应力下，试样摩擦因数μ越大，其静止侧压力系数越小，说明粗粒土表面粗糙程度对其静止侧压力系数有较大影响，土体表面越粗糙，静止侧压力系数越低。

根据各试样数值试验结果整理得σ′v分别在500、1 000、1 500、2 000及3 000 kPa时，对应不同μ试样的K0值，如图8所示。由图8可知，当试样所受竖向应力一样时，试样K0与其摩擦因数近似呈线性负相关关系，即

K0=B1-B2μ

(6)

式(6)中：B1和B2为材料参数。将加卸载时不同竖向应力下拟合所得堆石料材料参数B1和B2以及决定系数R2分别列于表4，所得拟合曲线绘于图8。

图7 不同摩擦因数试样的K0-σ′v曲线Fig.7 K0-σ′v curves of specimens with different friction coefficients

由图8可知，与对应试验值相比，式(6)预测值误差较小，最大误差只有2.25%，决定系数R2最小值也达到了0.952，显然，式(6)能较好地描述K0与摩擦因数μ的关系。

图8 不同竖向应力下各试样K0-σ′v关系曲线Fig.8 K0-σ′v curves of specimens with different vertical stress

表4 各试样不同竖向压力下的B1和B2Table 4 B1和B2 of each specimens under different vertical pressure

2.5 接触刚度对粗粒土K0影响

为分析接触刚度对粗粒土K0影响，采用D60试样级配，设计了切向接触刚度ks为1×107、2×107、4×107、6×107以及1×108N/m共5种情况下的K0试验，其他细观参数采用表2中所列参数。不同切向接触刚度试样的K0-σ′v曲线如图9所示。由图9可知，切向接触刚度对粗粒土静止侧压力系数有较大影响，在相同竖向应力下，随着切向接触刚度的增大，粗粒土静止侧压力系数逐渐增大，且随着切向接触刚度的增大，粗粒土静止侧压力系数增大趋势趋于平缓。显然，母岩强度对粗粒土静止侧压力系数有较大影响。

为研究法向接触刚度对粗粒土静止侧压力系数的影响，采用D60试样级配，设计了法向接触刚度kn为1×107、1×108、5×108及1×109N/m共4种情况下的K0试验，其他细观参数采用表2中所列参数。不同法向接触刚度试样的K0-σ′v曲线如图10所示。

由图10可知，在相同竖向应力下，试样法向接触刚度越大，其静止侧压力系数越小，且法向接触刚度小于1×108N/m时，静止侧压力系数随法向接触刚度的变化较大，而法向接触刚度大于1×107N/m时，法向接触刚度对静止侧压力系数影响较小。

图9 不同切向接触刚度试样的K0-σ′v曲线Fig.9 K0-σ′v curves of specimens with different tangential contact stiffness

图10 不同法向接触刚度试样的K0-σ′v曲线Fig.10 K0-σ′v curves of specimens with different normal contact stiffness

3 结 论

通过介绍粗粒土静止侧压力系数K0试验数值模型的建模过程，并根据数值试验结果研究了孔隙率、颗粒摩擦因数、切向刚度、法向刚度等宏细观参数对粗粒土静止侧压力系数的影响，得出以下结论。

(1)随着竖向应力增大，粗粒土的颗粒位移和力链逐渐增大。同一竖向应力下试样K0随着孔隙率的增大逐渐增大，但增大趋势趋于平缓，证明了粗粒土土体密实程度对静止侧压力系数有较大的影响。

(2)颗粒摩擦因数对粗粒土静止侧压力系数有较大影响，两者近似成线性负相关关系，粗糙程度越高，静止侧压力系数值越低。表明粗糙度对粗粒土静止土压力系数有较大影响。

(3)在相同竖向应力下，试样法向接触刚度越大，或切向接触刚度越大，其静止侧压力系数K0越小。

(4)结合数值试验结果，分别建立了在加载过程中粗粒土孔隙率和摩擦系数与静止侧压力系数K0的关系式。