土-石混合体系指块石和土体互相混合的岩土介质，其结构特征主要包括粒径、粒度分布、含石量、定向性及形状特征等^[1]。天然的土-石混合体赋存于残坡积物、崩坡积物和冲洪积物这样的天然堆积地层中，是完整新鲜岩体经过长期内、外动力地质作用的产物；人工的土-石混合体来自于工程建设产生的杂填土，块石形状更加复杂，分布更加不均匀。本次研究的土-石混合体，以块石形状规则且分布相对均匀、充填土体较密实为研究对象，主要突出块石和土体弹性模量的高低差异和强度差异。

目前已有学者对粒度分布和含石量进行过大量研究^[2-3]，而对定向性的研究^[4]却甚少。Wen等^[5]研究发现堆积层滑带土中粗颗粒定向性对强度的影响不可忽略。在深部滑移的堆积层滑坡中，土-石混合体在局部发生强烈剪切变形，从远离剪切带到靠近剪切带，块石定向性特征逐渐凸显，长轴方向趋向剪切变形方向。在隧道等水平线型工程中，纵向跨度范围内，土-石混合体中块石排列方向将随着冲积方向和沉积作用程度的不同而发生变化。在竖桩竖井等垂直线型工程中，含石量也随着沉积规律而有所变化。因此，考虑块石的含量及定向性对土-石混合体材料宏观力学性能产生的影响，在工程实践中十分必要^{[2, 6-8]}。

现有研究土-石混合体的方法主要有室内、原位试验法，以及基于有限元、有限差分或离散元的数值模拟方法。一些学者通过开展室内大型试验如直剪^[9]、三轴^[10]和水平推剪^[11]等试验，研究分析土-石混合体的变形和强度性质，基于此建立弹塑性本构模型^[12]。但是，大型室内试验具有经济和时间成本投入大、试样不可重复、试验过程不便于控制、试样内部细观特征难以观察等诸多不足，因而其发展受到制约。为了弥补室内试验的不足，徐文杰等^[13-16]、周伟等^[17]、Barbero等^[18]利用宏、细观力学特征探究了土-石混合体及堆石体的变形破坏过程。土-石混合体数值模拟试验中试样细观模型的构建是数值模拟试验至关重要的环节。目前构建方法主要有两种:一种是基于图像处理技术的真实模型生成方法^[16]，即选择某一特定揭露的土-石混合体区域，采用拍照或CT扫描方法并结合图像处理技术获得其细观结构，进而建立数值试验模型；另一种是基于统计分析原理的随机模型生成方法^[17]，即在对土-石混合体颗粒粒径进行统计分析的基础上，采用随机性颗粒生成方法建立试样的宏、细观模型。基于图像处理技术的建模方法只适用于局部特定区域，而基于随机模型生成的建模方法更加符合土-石混合体的实际空间特征和分布特点，能够适应多维、多尺度建模的需求。

笔者将建立二维椭圆块石集合体的细观随机模型，并进行数值模拟双轴压缩试验。该方法具有可重复性高、研究对象突出可控、细观变化明显的特点，可以研究土-石颗粒定向性特征对不同含石量的土-石混合体强度的影响，揭示土-石混合体局部与整体的作用机制。

1.   数值模拟试验

1.1   几何模型生成方法

笔者提出的二维椭圆块石集合体的随机模型(2D Elliptic Block Collection Stochastic Model, 简称EBCS^2D)，用于生成以椭圆颗粒为块石，其余范围为均质土的土-石混合体二维几何模型。

其中，等效椭圆模型一直是各研究^{[2, 7, 19-20]}中用于表征块石几何形态的重要模型。它将不规则的多边形块石简化为椭圆，其长轴的长度表示块石的最大可视粒径，长轴的方向表示块石的定向方向，长短轴之比表示圆扁形状特征。本次研究采用同样的处理方式来简化建模过程，排除不规则多边形块石对计算结果带来的离散性，突出对定向性特征的研究。

EBCS^2D模型通过调整椭圆形颗粒的级配、长轴长度、长轴方向、长短轴比来实现对土-石混合体的结构特征(如含石量、最大可视粒径、优势倾角等)的控制。实现步骤如下:

(1) 圆形颗粒填充  在PFC^2D程序中采用普通的接触模型，在与试样同大小的矩形区域，随机生成一定级配互不重叠的一系列的圆形颗粒，如图 1-a所示，然后导出圆形颗粒的中心坐标及半径。

图  1  EBCS^2D方法模型建立过程示意图

a.圆形颗粒填充；b.圆形颗粒筛选；c.椭圆形颗粒生成；d.混合模型生成

Figure  1.  Diagram of modeling procedure of EBCS^2D method

下载: 全尺寸图片 幻灯片

(2) 圆形颗粒筛选  通过Matlab编程，读取步骤(1)导出的数据。按照某一土-石阈值，筛选出将用作生成块石颗粒的圆形颗粒，其直径为块石颗粒的最大可视直径。图 1-b中蓝色圆形颗粒即为所选中用于生成块石的基体。

(3) 椭圆形颗粒生成  通过Matlab编程，将步骤(2)所筛选的圆作为椭圆形块石颗粒的外边界，以其圆心为椭圆中心，直径为长轴，按照一定的长短轴比、长轴倾角，绘制该圆所在范围内的椭圆。所生成的椭圆形颗粒即为土-石混合体中的块石颗粒，即图 1-c中黑色线条勾勒的椭圆。

(4) 混合模型生成  椭圆形颗粒集合与试样矩形区域的差集，即为土基质、土-石混合体二维几何模型的建成。经过网格划分得到图 1-d所示的计算模型。

1.1.1   含石量与长短轴之比的控制方法

土-石混合体几何结构模型的含石量，可通过椭圆形颗粒的数目、形状和大小来控制。在步骤(2)中，土-石阈值为d_min，最大粒径为d_max，则作为生成块石的圆形基体直径在[d_min, d_max]之间，即最大可视粒径的范围。由于步骤(1)中的PFC^2D颗粒模型中存在空隙，则用作块石生成域的圆形基体总面积S_ball为：

式中：μ为PFC^2D模型的孔隙率；W为试样的宽度；H为试样的高度；D为[d_min, d_max]粒径区间内圆形基体的总面积占试样总面积之比。

步骤(3)中，椭圆形颗粒长短轴之比P服从正态分布N(P₀, σ²)，其中σ²为方差，ρ₀为椭圆形颗粒长短轴之比的均值。在预算含石量时，为简化计算将所有椭圆形颗粒的长短轴之比都取均值P₀，这对所生成模型的含石量影响不大。根据椭圆与其外接圆的关系，有∑S_ellipse=P₀∑S_ball，可得(体积)含石量ABP(areal block proportion)：

μ、D在步骤(2)中通过PFC^2D程序设置，P₀在步骤(3)中通过Matlab代码设置，则土-石混合体模型的含石量及其块石的圆扁形状特性都得以控制生成。

1.1.2   定向性特征的控制方法

该方法以椭圆形颗粒长轴方向为块石颗粒的定向方向，定义竖直向上为标准0°方向，顺时针为正，椭圆形颗粒倾角θ在(-90°, 90°]范围内。在无优势倾角的模型中，倾角θ服从(-90°, 90°]范围的均匀分布；而在有优势倾角的模型中，倾角θ服从正态分布N(θ₀, α²)，定义均值θ₀为优势倾角, α²为方差。

椭圆形颗粒在互不相交的圆形基体内生成，椭圆形块石以圆形基体的中心为中心任意改变自身倾角，这样实现了含石量与定向性特征的独立控制，保持了块石集合体相对位置不变。

但是由于块石是在预留的圆形基体内生成的椭圆，块石之间的组合接触关系受到了限制。该方法仅能生成图 2-a所示的关系，而对图 2-b，c所示的关系则无法生成，这将在以后的研究中做进一步改进提高。

图  2  块石组合接触关系

Figure  2.  Contact relationship of block combination

下载: 全尺寸图片 幻灯片

1.2   数值模拟试验过程与参数

依照上述二维椭圆形块石集合体的随机模型生成方法，保持土-石混合体的部分结构参数(优势倾角分布方差α²、长短轴之比分布均值P₀、长短轴之比分布方差σ²)固定不变，仅对优势倾角均值θ₀和含石量ABP，取间隔变化的量(表 1)进行土-石混合体结构几何模型随机生成，并采用有限差分法进行数值模拟试验研究。

表  1  土-石混合体结构参数组合

Table  1.  Parameters of soil-rock mixture structure

块石定向性特征			块石形状特征		含石量ABP/%
θ0/(°)	α2		P0	σ2
无定向性			1.5	0.01	40.1829.0623.17
—	—
有定向性
0, 15, 30, 45, 60, 75, 90	0

下载: 导出CSV 

| 显示表格

所生成的土-石混合体二维几何模型，经过投影拉伸在FLAC^3D中建立相应的三维数值模拟试样，如图 3所示。模型尺寸为200 mm×100 mm×1 mm，经过均匀网格剖分，六边形网格边长均为1 mm。无定向性特征以及有定向性特征优势倾角45°的两种典型模型网格剖分如图 4所示。双轴压缩数值模拟试验的边界条件为:y轴方向位移固定为零；x方向底面的z方向位移固定为零；侧向荷载恒定200 kPa，垂直于x方向两侧表面向内施加压力；轴向荷载通过恒定速度边界条件施加，z轴方向速度为-0.5×10^-8 m/s。土体和块石的屈服面模型均采用Mohr-Coulomb模型；应力-应变关系均为线弹性。土体和块石的力学参数在各个不同的结构模型中保持一致，详见表 2。

图  3  几何模型

Figure  3.  Geometric model

下载: 全尺寸图片 幻灯片

图  4  典型模型网格剖分

含石量40.18%: a.无定向性；d.优势倾角45°
含石量29.06%: b.无定向性；e.优势倾角45°
含石量23.17%: c.无定向性；f.优势倾角45°

Figure  4.  Mesh in typical models

下载: 全尺寸图片 幻灯片

表  2  土-石混合体中块石和土体力学参数

Table  2.  Parameters of soil-rock mixture block and soil

名称	密度/(kg·m-3)	弹性模量E/MPa	泊松比μ	黏聚力c/kPa	内摩擦角φ/(°)
块石	2 580	1 040	0.2	500	38
土体	1 820	10	0.3	43	24

下载: 导出CSV 

| 显示表格

定向性特征(无，优势倾角0°~90°)与含石量(φ_B分别为40.18%、29.06%、23.17%)交叉组合试验24组，详见图 3。据文献[2]，选择一种典型的土-石阈值，取阈值粒径d_min=5 mm，最大粒径d_max=20 mm。

1.3   数值模拟试验结果

在数值模拟试验加载过程中，其应力-应变关系如图 5所示。各定向性特征土-石混合体的应力-应变曲线都在加载初期呈现阶梯状近线性上升。这一阶梯状曲线可能是速度边界条件下裂缝分布不均与扩展造成的。当加载进一步持续，土-石混合体的轴向应力不再增加时即表示试样达到了峰值抗压强度。

图  5  含石量40.18%试样的应力-应变图

Ⅰ、Ⅱ、Ⅲ为轴向应变的3个阶段

Figure  5.  Stress-strain curve of 40.18% block proportion

下载: 全尺寸图片 幻灯片

图 6和图 7分别表示了Ⅰ、Ⅱ、Ⅲ 3个应变阶段有定向性特征和无定向性特征的土-石混合体的最大剪应变率分布云图，表示的是剪裂缝的扩展滑移速率，图中网格的变形放大系数为4.0。

图  6  含石量40.18%无定向性特征试样加载过程最大剪应变率云图

Figure  6.  Nephogram of maximum shear strain rate in the soil-rock mixture of 40.18% block proportion with no preferred orientation

下载: 全尺寸图片 幻灯片

图  7  含石量40.18%有定向性特征θ₀=75°试样加载过程最大剪应变率云图

Figure  7.  Nephogram of maximum shear strain rate in the soil-rock mixture of 40.18% block proportion with preferred orientation θ₀=75°

下载: 全尺寸图片 幻灯片

在数值模拟试验中，土-石混合体试样均以剪切破坏形式达到极限强度，主剪切破坏面明显。从加载过程来看，整个试样范围内，局部剪切变形在初期零星分布。随着加载的持续，主剪切破坏面逐渐显现，其周围的剪切变形加剧，显著高于其他局部剪切变形。剪裂缝从局部扩展逐渐贯通，直至试样破坏。破坏试样的剪裂缝总体与轴向成30°~45°夹角，且仅发生于块石之间的土体中，局部形态受块石影响而发生弯曲。

无定向性特征的土-石混合体的抗压强度随着含石量升高而升高，这一结论与其他研究者^{[2, 15]}的结论吻合。以侧向压力200 kPa下不同定向性特征的土-石混合体为例，其三级含石量40.18%、29.06%、23.17%的抗压强度如图 8所示，分别代表较低、中等、较高3种含石量的土-石混合体。

图  8  不同含石量不同定向性特征下的抗压强度

Figure  8.  Compressive strength of the soil-rock mixtures with different block proportion and orientation

下载: 全尺寸图片 幻灯片

基于摩尔-库伦理论，侧向压力200 kPa下，本次研究采用的岩土强度参数如下：块石的抗压强度为2 891.06 kPa，土体的抗压强度为606.67 kPa(图 8中红线所示)。可见本次研究所选含石量的土-石混合体强度远小于块石强度，而接近土体强度。较低含石量的土-石混合体，在任何定向性特征下，其强度都略小于土体；中等含石量的土-石混合体，在优势倾角为15°、30°、45°的定向性特征下，其强度略小于土体强度，而在其他条件下，其强度均大于土体强度；较高含石量的土-石混合体，在任何定向性特征下，其强度都大于土体强度。

具有相同含石量的土-石混合体，其定向性特征对试样抗压强度的影响有显著不同。图 9为具有定向性特征试样的抗压强度相对于无定向性特征试样的变化率。含石量23.17%试样的抗压强度变化率在0.6%以内，含石量29.06%和40.18%试样的抗压强度变化率的大小最大分别达到2.7%和3.6%。从整体变化幅度来看，含石量越大，其抗压强度受到块石定向排列的影响越显著。较低含石量的土-石混合体，其抗压强度随着优势倾角增大而缓慢增大；而中等含石量和较高含石量的土-石混合体，随着优势倾角均值的增大，其抗压强度均呈现先减小后增大的规律，其中，最小抗压强度都在优势倾角30°~45°的结构特征下。

图  9  不同含石量定向性特征对试样抗压强度的影响

Figure  9.  Impact on the compressive strength of the soil-rock mixtures with different block proportion and orientation

下载: 全尺寸图片 幻灯片

2.   土-石相互作用机制

为探究土-石混合体的整体强度受块石影响的机制，笔者提出了土-石相互作用单元(土-石单元)的概念，它包括块石及其周围包裹的土体，代表受块石影响的局部变形破坏单元。笔者认为土-石混合体的整体强度取决于土-石相互作用单元局部作用的效果和这些局部作用效果的叠加。

通过分析模拟试验结果和应变云图发现，土-石单元在压缩和剪切两种局部变形形式下，呈现两种不同的作用效果，分别能够降低和提高局部土体强度。

观察发现，除了主剪切破坏面外，土-石界面是局部变形最集中的地方，多见张拉剥离和剪切滑移两种形式。这主要是土体和石块弹性模量的显著差异，二者容易发生不耦合变形造成的。这使土-石界面成为判断评价土-石单元局部作用效果的关键。

对于以局部压缩变形为主的土-石单元，与最小主应力方向垂直的土-石界面发生张拉剥离。该作用形成的局部张拉损伤裂隙随块石方向变化而发生开度和长度的变化，不同程度地降低局部土体强度。如图 10所示，局部压缩变形强度的变化顺序为：a < b < c。

图  10  土-石相互作用单元(土-石单元)

a，b，c.压缩变形；d，e，f.剪切变形

Figure  10.  Soil-rock interaction unit

下载: 全尺寸图片 幻灯片

对于以局部剪切变形为主的土-石单元，土-石界面上拉压应力交替分布，块石起到阻碍剪切变形的作用，与剪切变形同向的土-石界面发生剪切滑移。土-石单元受块石方向的影响，块石方向与主剪切面夹角越大，对局部土体强度的提高越显著。如图 10所示, 局部剪切变形强度的变化顺序为：d>e>f。

以上两种局部变形形式下的土-石单元，对局部土体强度呈现相反的作用效果，且以不同方式受到块石方向的影响。在本次模拟试验的范围内，这两种土-石单元的局部作用在整体中所占的份额与块石的疏密程度相关。

当块石之间间距较大、块石之间的影响不明显时，土-石界面的变形破坏可以无阻碍地发生，为压缩变形下的作用效果提供了有力条件。同时，稀疏的块石导致主剪切破坏面仅经过数量很少的土-石单元，限制了剪切变形的土-石单元所发挥的作用(图 11)。与上述块石孤立作用相对的是块石团簇作用，即当块石之间间距较小、块石之间的影响不可忽略时，块石的土-石界面不能自由地变形破坏，而是受到相邻块石的约束，从而削弱了压缩变形的土-石单元的作用效果。另外，密集的块石团簇增加了主剪切破坏面穿过土-石单元的机会，使得更多的土-石单元在剪切变形下发挥提高局部土体强度的作用(图 12)。

图  11  块石孤立作用

Figure  11.  Isolated rock effect

下载: 全尺寸图片 幻灯片

图  12  块石团簇作用

Figure  12.  Clustered rock effect

下载: 全尺寸图片 幻灯片

以上作用机制可以解释为本次模拟试验中，不同含石量下，土-石混合体强度受定向性影响的内在原因。当含石量较低时，块石孤立作用主导，土-石混合体的强度受到压缩变形的土-石单元控制，任何定向性特征下都低于土体，且有定向性特征的土-石混合体强度随着块石优势倾角的增大而增大；当含石量较高时，块石团簇作用主导，土-石混合体的强度受到剪切变形的土-石单元控制，任何定向性特征下都高于土体，且有定向性特征的土-石混合体强度随着块石优势倾角与主剪切面夹角的减小而减小。当含石量中等时，块石孤立作用和团簇作用难分高下，土-石混合体的强度，在主剪切破坏面受到剪切变形的土-石单元影响，其他则受到压缩变形的土-石单元影响，不同定向性特征下，可能高于土体也可能低于土体。

3.   各向异性屈服面的适用性

Oda等^[21]通过对砂岩中颗粒定向性特征的研究，提出了扩展的Drucker-Prager屈服面方程，基于其系列研究提出用各向异性组构张量F_ij来描述其中颗粒的定向性特征，其主方向表达式如下:

其中

式中：M是颗粒个数；θ(k)是第k个颗粒的倾角；0≤Δ≤1为表征倾角密度的参量^[22]，当Δ=1时表示颗粒倾角θ全部一致。

考虑颗粒定向性特征的扩展的Drucker-Prager屈服面方程^[21]如下：

式中：α, c为试验所得的参数，c一般取零,

式中：δ_ij, δ_ik, δ_lj分别为不同方向的克罗内克函数；σ_ij为球应力张量；S_ij, S_kl分别为不同方向的偏斜应力张量；F_ij, F_ik分别为不同方向的各向异性组构张量。

各参数含义如下:

式中：λ，μ，Δ为试验相关参数。

将第一节数值模拟试验中块石的定向性特征参数、侧向压力和峰值强度，代入上述宏观各向异性屈服面方程，分别拟合得到三种含石量下的参数μ、λ、α(表 3)。拟合结果如图 13所示，由表 3和图 13可以看出，随着优势倾角的变化，抗压强度变化规律基本吻合，决定系数均接近1，拟合效果较好，说明这一各向异性屈服面方程在一定程度上适用于土-石混合体。

图  13  数值模拟结果与拟合曲线

Figure  13.  Numerical simulated results and the fitting curve

下载: 全尺寸图片 幻灯片

表  3  拟合参数结果

Table  3.  The results of fitting parameters

含石量φB/%	μ	λ	α	决定系数
40.18	0.786 62	0.997 74	0.245 24	0.999 88
标准差	0.042 21	0.066 38	0.000 89
29.06	0.641 80	0.769 78	0.236 39	0.999 97
标准差	0.014 01	0.013 23	0.000 69
23.17	0.065 14	0.081 56	0.232 17	1.000 00
标准差	0.138 54	0.146 25	0.000 53

下载: 导出CSV 

| 显示表格

由表 3可见，随着含石量的降低，参数μ、λ、α不同程度地减小，其中μ和λ减小显著。张佩等^[4]用抛物线拟合得到的参数A和参数B也随着含石量有单调变化趋势。该规律可以在一定程度上反映土-石混合体强度对其定向性特征的敏感性，以及含石量对该敏感性的影响，而这些参数与土-石混合体含石量相关联的物理意义还有待探明。

4.   堆积层边坡模型

因为堆积层边坡失稳十分普遍，且岩土体性质在不同条件下变化复杂，堆积层边坡一直备受关注^[23-25]。本次研究以理想堆积层边坡为例，以坡体中块石的定向性特征为研究重点，建立理想的土-石混合体的堆积层边坡模型。模型的坡高30 m，坡面倾角30°。堆积层土-石混合体的含石量30%，岩土材料参数与前文一致，块石粒径范围50~150 mm。分别建立块石无定向分布(图 14-a)和块石定向分布(图 14-b)2种特征的堆积层边坡模型，后者定向分布的块石长轴方向与坡面方向一致。由强度折减法求解得出后者的稳定性系数低于前者，说明在块石顺坡向分布的堆积层边坡中，滑动面更容易在土体中延伸，土-石混合体的块石定向性特征在一定程度上影响边坡稳定性。更多的关于块石定向性、含石量、粒径、边坡倾角等对堆积层边坡稳定性的影响还有待进一步研究。

图  14  堆积层边坡模型

Figure  14.  Debris slope model

下载: 全尺寸图片 幻灯片

4.   结论

(1) 本次研究建立的二维椭圆块石集合体的随机模型(EBCS^2D)能在合理简化模型的基础上，控制含石量、定向性、长短轴之比这些结构特征变量，生成研究所需的几何模型。它在研究土-石混合体结构特征对其强度影响的建模过程中，具有较好的适用性。

(2) 数值模拟结果显示，含石量不同，块石定向性对土-石混合体总体性能的影响程度也不同。土-石混合体的含石量越大，自身强度越大，受块石定向性特征的影响越显著。

(3) 土-石混合体的整体强度由土-石单元的局部作用和不同分布的局部作用效果叠加决定。土-石相互作用单元，在压缩和剪切局部变形条件下，对其周边土体强度作用的效果相反。整体上，块石的疏密程度直接影响压缩和剪切两种土-石单元局部作用效果的数目和强弱。块石孤立和团簇两种作用的叠加解释了不同含石量的土-石混合体强度受块石影响不同的原因。

(4) 各向异性屈服面方程的计算结果与数值模拟结果具有较高相关性，说明该屈服面模型的参数μ和λ能用于表征土-石混合体中块石对其强度的影响。