考虑形状因素的岩土细观力学特性研究

Research on Mesomechanical Properties of Rock and Soil Considering Shape Factors

摘要: 本文基于PFC6.0软件的数值模拟工具，针对三轴压缩试验中的颗粒形状这一因素进行了深入研究，在引入了颗粒细长比AR、颗粒凹陷程度C这两个因素后，构建四颗粒模型来模拟真实岩块形状，并从旋转角、各向异性、滑移接触占比等细观角度进行了力学特性的研究，最后与圆形颗粒模型的数值模拟结果作对比，得到了在重点考虑变形的情况下，圆颗粒会低估位移量的结论，造成工程风险。

Abstract: Based on the numerical simulation tool of PFC6.0 software, this paper conducts an in-depth study on the factor of particle shape in triaxial compression test. After introducing the two factors of particle slenderness ratio AR and particle depression degree C, four particles are constructed. The model is used to simulate the shape of the real rock block, and the mechanical characteristics are studied from the micro-angles of rotation angle, anisotropy, and slip contact ratio. Finally, the numerical simulation results of the round particle model are compared, and the key points are obtained. Considering the deformation, round particles will underestimate the conclusion of the displacement, causing engineering risks.

1. 引言

岩堆体 [1]，是岩质山体在各种物理、化学作用下发生塌滑、剥落，形成由大小不一岩石碎块、岩屑组成的松散堆积体 [2] [3]，其性质完全不同于土或岩石，表现出强烈的非均质性和各向异性 [4] [5] [6]。

国内外的研究表明岩堆体的力学行为受诸多因素不同程度的影响，如含石量、含水率、围压、块石形状、块石的空间分布、土石强度比、试样尺寸等 [7]，李维树等 [8] 针对奉节新城的岩堆体，做了大量的原位直剪试验，建立了强度参数与含石量之间的对应关系。李晓等 [9] 在三峡库区自农庵滑坡处典型的岩堆体上，在现场进行了23个大尺度原位推剪试验和压剪试验。高春玉等 [10] 对四川盆地区红层无粘性土石混合料进行多组大型三轴试验。Xu等 [11] 基于云南丽江金沙江中游梨园电站坝址区的岩堆体边坡，研究了四组含石量(0%，30%，50%和70%)下的岩堆体强度特性。王凯 [12] 针对G314国道公格尔隧道进出口段的土石混合体地层做了现场原位试验。但就目前所用到的试验方法而言，主要采用直剪试验或水平推剪试验，难以进行更加复杂的力学试验，也难以研究其细观力学特性。

此外，随着电脑计算能力以及数值分析技术的发展，越来越多的研究者开始使用数值试验来研究岩堆体细观力学行为和变形破坏机制。徐文杰等 [13] [14] [15] 提出了一种基于数字图像处理的非均质岩土材料细观结构PFC2D数值计算模型自动生成方法，开展了一系列数值直剪试验研究工作。张振平等 [16] 生成团粒结构建立土石混合体直剪试验离散元模型，研究不同含石量下应力峰值强度的变化。三维离散元计算方面，金磊等 [17] [18] 使用不规则颗粒三维离散元精细模拟技术，考虑块石形状为球体、正方体和长方体三种情况，探究不同棱角对胶结岩堆体力学特性的影响。张强 [19] 综合运用计算机三维扫描与随机模拟技术，建立了不同块石含量和空间分布的土石混合体三维随机细观结构模型和离散元模型，采用颗粒流程序进行大型离散元三轴试验模拟，研究了块石含量和空间分布对土石混合体力学特性和变形破坏规律的影响。

纵观国内外对岩堆体的力学特性研究，数值方法已经成为了一个重要的手段，但目前的研究大多将块体简单地简化为规则几何体，如球体或方形，这一定程度上是受限于计算机的算力，但过度的简化会忽略块体表面凹凸部分能嵌固在一起形成一个整体，加大摩擦，提高整体强度的事实，忽略了颗粒形态对细观力学特性的影响。

为了研究不同形状土石混合体的细观力学特性，论文提出了一种基于真实块石形状的土石混合体二维建模方法，并展开了三轴压缩数值试验，并以此为基础开展了岩堆体的细观力学特性研究。

2. 三轴压缩数值试验

2.1. 模型建立

离散元计算方法在处理大变形、破裂过程具有天然的优势，因此本文运用PFC2D6.0模拟进行三轴压缩数值试验。PFC中的粒子为圆形刚性粒子，为了模拟复杂的颗粒形状，采用以下两个颗粒形状参数，AR (Aspect ratio)以及C (Convexity)，前者表示颗粒的细长比，后者表征颗粒的凹陷程度，如图1所示。

AR = D^Fmin/D^Fmax C = A/(A+B)

Figure 1. Schematic diagram of shape parameters

图1. 形状参数示意图

为了兼顾算力和形状参数，现采用四颗粒来表示一个clump单元，如图2所示，首先其细长比需满足真实块体的AR参数，其次根据真实块体的凹陷程度C，通过matlab建立方程可解出小球半径r。

${\theta }_1=\arccos \frac{AR-2r}{2r}$

${\theta }_{\text{2}}=\arccos \frac{\text{1}-2r}{2r}$

$\frac{AR-\left(4-\pi \right)r^2-2\left[r(AR-2r)-\pi r^2\frac{90-{\theta }_1}{180}-\frac{1}{2}{r}^2\sin 2{\theta }_1\right]-2\left[r(1-2r)-\pi r^2\frac{90-{\theta }_2}{180}-\frac{1}{2}{r}^2\sin 2{\theta }_2\right]}{AR-\left(4-\pi \right)r^2}=C$

其中

Figure 2. Particle model

图2. 颗粒模型

具体参数标定参考了Luo X D针对福建土的力学特性研究室内试验 [20]，其模型盒宽71.1 mm，高142.2 mm，颗粒粒径级配以及形状参数级配如图3所示。

最终建立整体模型如图4所示。

模型用到的具体参数如下表1所示。

2.2. 标定结果

三轴压缩数值试验和Luo X D的室内试验结果对比如下图5所示，可以看到，在偏应力–轴向应变曲线图与体积压缩应变曲线图中，标定数值试验与室内试验的结果均十分接近，可以采用作为后续细观力学性能数值试验的参数。

Figure 3. Particle shape parameters and diameter gradation diagram

图3. 颗粒形状参数及直径级配图

Figure 4. Schematic diagram of PFC2D6.0 generated model

图4. PFC2D6.0生成模型示意图

Table 1. Basic parameters of Fujian soil particles

表1. 福建土颗粒基本参数

Figure 5. Comparison of calibration results with indoor test results

图5. 标定结果与室内试验结果比较

更进一步，在改变围压至200 kPa~1500 kPa范围内时，福建土的破坏线与室内试验所得到的破坏线十分相近，如图6所示。

Figure 6. Destruction line graph of Fujian soil

图6. 福建土的破坏线图

其中横坐标 $p=\frac{{\sigma }_1+{\sigma }_3}{2}$ ，纵坐标 $q=\frac{{\sigma }_1-{\sigma }_3}{2}$

该图像的拟合直线斜率为1.239，而原试验的拟合直线为1.21，误差仅有2%，已十分接近。

3. 细观力学特性分析

在PFC三轴压缩数值试验中，每隔0.5%的轴向应变就会遍历导出所有颗粒数据及接触信息，包括颗粒的位置、速度、不平衡力、旋转角、角速度、所属分组、体积、ID，以及接触的位置、轴向应力、切向应力、轴应力的方向、连接颗粒的id，以下分析是基于这些数据所处理得到的。

3.1. 旋转角

将颗粒簇的旋转角度提取出来，每5%轴向应变之间做一个差值，可以看到每个小阶段下颗粒旋转的程度。下图7用颜色来区分顺逆方向，蓝色为顺时针，红色为逆时针，且颜色越深转动角越大。

0%~5%轴应变 5%~10%轴应变 10%~15%轴应变 15%~20%轴应变 20%~25%轴应变 25%~30%轴应变

Figure 7. Rotation angle of particle clusters at different axial strain stages

图7. 不同轴应变阶段颗粒簇旋转角度

可以看到，在前10%轴应变阶段，旋转程度较大的颗粒簇分布还较为均匀，对应偏应力应变图看的话，此阶段试样整体还未进入屈服；而在15%轴应变之后，颜色深的颗粒簇开始集中在某些条带上，象征着滑动带的产生，且这些滑动带大多与水平面呈45˚夹角。

该规律的出现与偏应力–应变图中曲线进入屈服的时刻相吻合，试样整体压缩进入屈服阶段在细观层面表现为滑动带的产生。

3.2. 各向异性

在考虑形状因素的情况下，颗粒簇是有一个朝向的，生成模型时朝向0~360˚完全随机，所以是各向同性，但在经历固结–压缩之后，由于细长比的存在，颗粒簇朝向可能不再等可能地朝向任意方向，因此有必要用一个参数来表征压缩过程中各向异性的变化。

在二维平面中，颗粒的各向异性可以根据Rothenburg [21] 提出的傅里叶级数来计算：

$f_n\simeq \frac{1}{2\pi }\left[1+a\cos 2\left(\theta -{\theta }_a\right)\right]$

其中 $a$ 表示各向异性系数， ${\theta }_a$ 表示主角度，并满足 ${\displaystyle \oint f\left(n\right)}dn=1$

而根据Oda [22] 所提出的张量计算式可以用来计算接触力的各向异性：

${\Phi }_{ij}\text{=}{\displaystyle \underset{\Theta }{\int }E\left(\Theta \right)n_i{n}_{j}d\Theta }=\frac{1}{N_c}{\displaystyle \underset{c\in N_c}{\sum }{n}_i{n}_j}$

此处n表示接触应力朝向的单位向量，N_c表示接触的数量，在二维空间下，公式可进一步写为：

$\left(\begin{array}{cc}{\Phi }_{xx}& {\Phi }_{xy}\\ {\Phi }_{yx}& {\Phi }_{yy}\end{array}\right)=\frac{1}{N_c}\left(\begin{array}{cc}{\displaystyle \underset{k=1}{\overset{N_c}{\sum }}{n}_{kx}{n}_{kx}}& {\displaystyle \underset{k=1}{\overset{N_c}{\sum }}{n}_{kx}{n}_{ky}}\\ {\displaystyle \underset{k=1}{\overset{N_c}{\sum }}{n}_{ky}{n}_{kx}}& {\displaystyle \underset{k=1}{\overset{N_c}{\sum }}{n}_{ky}{n}_{ky}}\end{array}\right)$

Oda等人指出，二维空间内各向异性参数可以由上式二阶应力张量的两个主应力差值求出，即

$a={\sigma }_1-{\sigma }_3$

最终作出在30%轴应变内各向异性系数的变化曲线图，如下图8所示。

Figure 8. Variation curve of anisotropy coefficient with axial strain

图8. 各向异性系数随轴应变变化曲线

图中轴向应变为0的状态为施加均等围压的固结后状态，可以看出固结完毕时结构已经不是各向同性的了，但系数大约只有0.11左右，大体还是保持各向同性状态的；而进入压缩阶段后，各向异性系数会突增至0.53左右，并在之后很长一个阶段内保持几乎不变的状态，说明颗粒细长比不为1的状态下，压缩后的试样内结构不是各向同性的，有一定的偏好方向。

3.3. 滑动颗粒占比

遍历所有pebble-pebble接触，统计满足 $f_n=\mu f_s$ 的接触个数，其中 $f_n$ 是法向接触力， $f_s$ 是切向接触力， $\mu$ 是摩擦系数，两颗粒若摩擦系数不同，取大者。满足上式的接触即为滑动接触，将压缩阶段不同时刻下滑动接触占所有接触的比例汇总，做出图9。

可以看到，在整个压缩阶段，颗粒与颗粒之间产生滑动的占比是几乎保持不变的，不与是否进入屈服阶段相关。

3.4. 与圆形颗粒的比较

不采用四颗粒模型，不考虑形状因素，仅用一个圆球来模拟，粒径级配以及模型基本参数均保持不变，也进行了对比数值试验，其比较如下图10所示。

Figure 9. Sliding contact ratio curve

图9. 滑动接触占比曲线图

Figure 10. Comparison of compression curves between round particles and four-particle models

图10. 圆颗粒与四颗粒模型的压缩曲线比较图

可以看出，忽略形状因素的条件下模型的偏应力应变曲线几乎没有变化，而体积压缩曲线却少了约三分之一，说明形状因素主要影响的是模型的变形能力。

通过计算得到四颗粒模型在固结完毕时的孔隙比为0.295，而圆形颗粒模型在固结完毕时的孔隙比为0.304，结合之后的体积压缩应变曲线，说明不规则的形状更易堆积得更密，且在偏应力作用下角度调整后能够堆积得愈密，相反圆形颗粒的堆积就比较松散，初始空隙较大，且更难压密。

做出圆形颗粒模型的各向异性系数在压缩过程中的变化曲线图如下图11，可以看到其各向异性也远不如不规则四颗粒模型的大，且圆颗粒在偏应力为零的固结阶段后，其各向异性系数也为零，不像四颗粒模型在固结阶段就开始出现轻微的各向异性。

Figure 11. Comparison of anisotropy between round particles and four-particle models

图11. 圆颗粒与四颗粒模型的各向异性比较图

综上比较可得，在仅考虑强度的时候，颗粒形状影响很小，可以使用圆颗粒进行简化；而当考虑变形的时候，使用圆颗粒模型计算得到的形变会变小，在用数值模型进行实际工程中的变形模拟时，很可能低估地层的收敛，造成风险。

4. 结论与展望

4.1. 结论

本文通过细长比AR与凹度C两个形状参数构建了用来模拟真实块体形状的颗粒簇，在PFC软件下进行了三轴压缩数值试验，并进行了一系列细观力学分析，得到以下结论：

1) 在整体试样压缩进入屈服阶段时，细观表现为滑移带的产生，转动较大的颗粒开始集中在某些与水平呈约45˚夹角的带上；

2) 在四周受到同样围压下，颗粒的长细比对整体各向异性影响较小；而当偏应力产生时，试样整体内部便会表现出较大的各向异性；

3) 滑动接触的数量占比与试样模型是否进入屈服阶段无关，整体保持在一个定值附近；

4) 在进行强度计算时，颗粒形状的影响很小，可采用圆颗粒进行简化计算；在进行收敛位移计算时，圆颗粒会低估形变量，造成工程风险，宜采用接近真实的颗粒形状模型进行数值模拟。

4.2. 展望

本论文模拟的岩块颗粒形状虽然考虑了长细比和凹度，但使用的形状样本仍然不足，实际工程遇到的岩块可能极薄，甚至有的凹凸会形成锁扣，且数值模拟未考虑岩块的破裂，将来可以从以上几个角度进行更深入的研究。