黄土颗粒结构特征及其对剪切行为的影响*

苑伟娜 范 文 邓龙胜 李维婉 周彦妍

(长安大学， 地质工程与测绘学院， 西部矿产资源与地质工程教育部重点实验室， 西安 710054， 中国)

0 引 言

土的结构是土颗粒或集合体在空间上的排列以及它们之间的联结，分为颗粒骨架、孔隙和联结3个方面的特征。黄土是典型的结构性土，微结构控制着黄土的物理力学性质，进而影响着工程土体的稳定性(彭建兵等， 2014； 邵生俊等， 2014； Ng et al.，2017； Luo et al.，2018)。研究黄土微结构及其对力学行为的影响是黄土力学的重要内容(谢定义， 2001)。

黄土微结构一直是国内外学者研究的热点问题之一，随着科学观测技术的不断提高，黄土微结构的研究成果在逐步深化。许多学者先后借助激光粒度仪、压汞仪、显微镜，扫描电镜等手段，对黄土的微结构进行了大量的研究工作，揭示了黄土的颗粒、孔隙和胶结特征(张宗祜， 1964； 高国瑞， 1980； Pye， 1995)，实现了黄土的二维微结构定量化研究，也提出了“结构强度”、“综合结构势”等结构性表征参数(谢定义等， 1999)。近年来，一些学者开始采用计算机断层扫描(CT)技术开展黄土三维微结构相关的研究工作(蒲毅彬等， 2000； 朱元青等， 2009； Li Y R et al.，2018； Wei et al.，2019; Li P et al.，2020； Yu et al.，2020； Wei et al.，2020)，取得了一系列成果，但黄土三维微结构量化及建模方面的研究工作仍在探索阶段。

国内外学者也关注到黄土微结构对其力学性质的影响并展开了相关研究。Matalucci et al. (1970)最早通过室内直剪和三轴试验研究了黄土颗粒排列对剪切强度的影响。一些学者通过对比分析原状黄土试样与重塑黄土试样剪切试验的结果来间接研究结构对黄土力学性质的影响及相应的机理(Wen et al.，2014； Xu et al.，2016； Ng et al.，2017)。许多学者也对黄土试样在载荷作用下的结构演化展开了研究(Hu et al.，2001； 雷胜友等， 2004； Yuan et al.，2019)。这些研究工作主要基于试验方法从不同角度探讨了黄土微结构对力学行为的控制作用。但是黄土的微结构特征是颗粒结构、孔隙结构和胶结三方面综合的结果，各方面都包含多个结构要素，试验的方法难以控制单个结构要素对于力学性质的影响，而深入分析各个结构要素的影响是揭示黄土宏观力学行为微观机制的基础。数值模拟方法具有这方面的优势，其中：颗粒流离散元数值模拟方法在岩石、砂和土石混合体等地质材料的结构及力学性质方面的研究中得到了广泛的应用。李识博等也采用该方法研究黄土的力学性质(李识博等， 2013； Jiang et al.，2017； 同霄等， 2019)，蔡玮彬等(2014)借助该方法分析了黄土的胶结特征，这些研究也验证了颗粒流离散元数值方法对于黄土微结构和力学性质研究的适用性。因此，采用颗粒流离散元数值方法研究黄土骨架颗粒的各个结构要素对力学性质的影响是必要和可行的。

本文拟通过X射线CT技术观测黄土微结构，借助Avizo软件定量分析黄土的三维颗粒结构特征，结合黄土的颗粒结构特征及室内三轴试验建立黄土的颗粒流离散元数值计算模型，研究黄土颗粒结构的各个要素对其剪切特性的影响，研究结果有助于揭示黄土宏观力学行为的微观机制，丰富黄土微结构及黄土力学的研究成果。

1 试样与试验

陕西泾阳南塬塬边黄土滑坡灾害频发，塬边斜坡体主要由中更新世(Q2)黄土构成(雷祥义， 1995)。本文试验所用黄土试样为陕西泾阳南塬Q2黄土，取样深度约20m，现场获取30cm×30cm×30cm的块状原状黄土试样。依据《土工试验方法标准》(GB/T50123-2019)测得试样的基本物理性质指标，试样的天然含水量为15.2%，天然密度为 1.79 g·cm-3，比重瓶法测得土粒密度为 2.69 g·cm-3，计算得到试样的孔隙度为42%。采用液塑限联合测定仪测得试样的液限为31.8%，塑限为17.0%。

采用X射线CT设备分析黄土试样的三维颗粒结构特征。从块状原状黄土试样中，加工出约直径2mm，高2mm的圆柱试样，保持试样横截面平行于黄土地层。采用中国科学院地质与地球物理研究所的ZEISS X radia 520 Versa微米CT设备对黄土圆柱试样进行扫描，该仪器的空间分辨率为700 nm，本文扫描试验使用的分辨率为1μm，试验时穿透电压为50kV，电流为80 μA。

为了验证黄土试样的颗粒流数值计算模型，对原状黄土试样进行了三轴剪切试验。三轴试验所用原状试样的尺寸为直径39.1mm，高80mm，试验采用的仪器为GDS非饱和土三轴剪切仪，剪切方法为固结排水剪，剪切速率为0.02%/min，围压为100 kPa，剪切应变达15%时停止测试。

2 黄土颗粒结构特征

采用Avizo软件对黄土试样的CT扫描图像进行处理和定量分析(图1)。首先通过阈值分割的方法提取黄土试样的颗粒骨架，依据试样的孔隙度确定阈值，阈值分割后得到试样的二值化图像(图1b)。然后采用Avizo软件中基于分水岭的算法实现颗粒骨架中颗粒之间的分割(图1c)。最后对分割后的单个颗粒进行统计和量化分析。颗粒尺寸、形态和排列是颗粒结构特征的3个基本方面(Assallay， 1998)，本文从这3个方面来分析黄土试样的颗粒结构特征。

图1 CT图像及处理结果

颗粒尺寸采用颗粒的等效直径(EqD)来表征(式(1))，考虑到扫描试验的分辨率，文中主要关注等效直径大于1μm的颗粒。试样的统计结果见图2a，结果表明，试样颗粒的等效直径小于70μm，平均粒径为27.6μm，粒径分布满足Weibull概率分布，概率分布的比例参数λ为30.37μm，形状参数k为2.33。

EqD=(6V/π)(1/3)

(1)

式中：V为颗粒的体积。

描述颗粒形态的参数较多(Blott et al.，2008)，本文采用式(2)和式(3)表示的球度来表征颗粒的形态(Fonseca， 2011)。颗粒的球度值介于0～1之间，取1时颗粒为球形，数值越大，颗粒的形状越接近于球形或者各向等轴。试样的三维颗粒球度统计结果见图2b，球度值分布符合Gamma分布。

S3D=(36πV2/A3)(1/3)

(2)

S2D=(4πA2D/P2)(1/2)

(3)

其中：S3D为三维球度；V为颗粒的体积；A为颗粒的表面积；S2D为二维球度；A2D为颗粒的面积；P为颗粒的周长。

三维颗粒排列通过球坐标系下颗粒最大Feret直径的方向角(φ，θ)分布来表征。如图2c所示，坐标系的Z轴垂直于地层，颗粒方向角在坐标系上半空间分成31个区间，每个区间用颗粒最大Feret直径与单位球面的交点表示，则图中的31个区间表示成31个点。图2d为图2c沿Z轴负方向的俯视图，图中给出了黄土试样颗粒方向角在31个区间的统计结果，各个区间的统计频率分别用相应点的直径大小表示。图2d中的统计结果表明，沿坐标系的方位角θ方向，颗粒方向角的分布相对均匀； 随着坐标系的极角φ增大，颗粒方向角的分布频率有增大的趋势，在φ接近90°时(颗粒方向平行于地层)，颗粒方向角的分布频率最大(图中的点直径最大)，φ介于60°～90°之间的颗粒占60%。这表明颗粒的排列表现出横观各向同性的特征，并且大部分颗粒最大Feret直径方向与地层近于平行。

图2 颗粒结构特征统计结果

3 建立数值计算模型

考虑到黄土微结构的横观各向同性排列特征，以及数值模拟方法的计算效率，文中借助PFC2D软件来分析颗粒结构对黄土剪切行为的影响，采用椭圆形的clump单元来表征颗粒，颗粒之间的接触模型为线性接触黏结模型。如图3a中的放大区域所示，数值试样中颗粒的尺寸分布按照图2a中试样的统计结果生成，颗粒的长轴方向角随机均匀分布，颗粒的形状取相同的球度值。已有的研究表明三维颗粒形态参数与二维参数结果有一定的差别(Fonseca， 2011； Alshibli et al.，2015)，因此，球度值采用二维球度统计结果的峰值。

PFC软件中试样的加载可以通过墙体或边界上的颗粒单元来实现。由于室内三轴剪切试验的围压通过水压来施加，是一种柔性加载方式。为了客观地模拟试验的加载条件，文中通过边界上的球颗粒对试样进行加载。按上述的结构特征生成黄土数值试样后，在试样的两侧边分别生成边界球(图3a中黄色边界球)，边界球之间的接触设置较大的强度参数，以保证边界的完整性和加载的连续性，试验时通过对两侧边边界球施加力来伺服控制围压，实现柔性加载。在试样的上下边也生成边界球(图3a中绿色边界球)，通过控制上下边的边界球速度来完成轴向加载。

图3 剪切试验前后试样

为了建立能够表征黄土试样剪切性质的数值模型，完成数值试样建立和加载设置后，参考黄土颗粒流离散元数值计算相关的文献(李识博等， 2013； 同霄等， 2019)，选取PFC数值计算模型微观力学参数的初值，经过调整和试算，并与黄土三轴剪切试验的结果进行对比验证，最后确定计算模型的微观力学参数，见表1，数值试验与三轴试验的结果对比见图3和图4。结果表明，数值模拟计算所得的偏应力-轴向应变关系曲线形态与三轴试验结果基本一致，表现出轻微应变硬化现象(图4)，在100kPa围压下数值计算和试验得到的抗剪强度分别为352.7kPa和333.6 kPa，两种方法得到的试样形态也基本一致，为塑性剪切破坏，侧向有轻微鼓胀现象(图3)。

表1 数值模型微观力学参数

图4 偏应力-轴向应变关系曲线对比

4 黄土颗粒结构对剪切行为的影响

4.1 颗粒尺寸的影响

黄土试样的颗粒等效直径分布满足Weibull概率分布(图2a)，如图5a所示，改变形状参数k的值，概率分布曲线呈现出不同的形状，k值越大，粒径分布相对越均匀。为了分析颗粒尺寸分布对黄土剪切特性的影响，保持图3a中试样的颗粒形状和排列方式，分别按照k=2.3、k=6.4、k=9和均匀分布时4种颗粒尺寸分布生成试样，采用表1中的力学参数，进行数值剪切试验，试验结果见图5b和图5c。结果表明，试样颗粒的粒径分布越均匀，应变硬化现象越显著，试样的抗剪强度越高(图5b)，但各个试样的变形破坏形态没有明显的差异(图5c)。这主要是因为试样的抗剪强度由黏结强度和摩擦强度组成，计算时所有试样取相同的颗粒接触模型参数，黏结强度差异小，但试样的粒径分布越均匀，颗粒之间的咬合摩擦越大，具有的抗剪强度越高。并且在微观上，各个试样变形破坏的本质都是颗粒不断调整最终形成塑性剪切带的过程，因此不同粒径分布的试样表现出的宏观变形破坏形态没有明显的差异。

图5 不同颗粒尺寸分布条件下试样的数值剪切试验结果

4.2 颗粒形态的影响

为了分析颗粒形态对剪切行为的影响，保持图3a中试样的颗粒尺寸分布和排列方式，分别取颗粒球度值1、0.94、0.87和0.8生成4个试样，采用表1中的力学参数，进行数值剪切试验，计算结果见图6。结果表明，试样的变形破坏方式相同，颗粒球度值越小，试样的变形模量越小，应变硬化现象越显著。球度值越大，抗剪强度也相对越大，但圆形颗粒试样的抗剪强度反而较低。文中建立的模型颗粒表面粗糙度基本相同，在此情况下，如图6a所示，颗粒的球度值越小，颗粒越细长(长短轴比越大)，在受力剪切过程中，颗粒旋转或滚动形成定向排列剪切带的难度增加，因而应变硬化现象越显著。与椭圆形颗粒相比，圆形颗粒之间的咬合力相对较小，因此圆形颗粒试样具有较低的剪切强度。

图6 不同颗粒球度值试样的数值剪切试验结果

4.3 颗粒排列的影响

数值模型中，clump单元长轴方向即为颗粒最大Feret直径方向，颗粒排列通过颗粒方向角(α)的分布来表征。考虑颗粒排列对试样剪切性质的影响时，试样颗粒形状与图3a中的试样一致，颗粒尺寸均匀分布，分别按照图7a所示的3种排列方式生成试样，进行数值剪切试验。其中排列方式1的试样颗粒方向集中分布在与轴向加载方向近平行的方向，排列方式2的试样颗粒方向随机均匀分布，排列方式3的试样颗粒方向集中分布在与轴向加载方向近垂直的方向。

试样的数值剪切计算结果见图7b和图7c，结果表明，排列方式对于试样的剪切特性影响较显著，排列方式3的试样有较大的峰值强度和应力降，有贯通的剪切破坏面。排列方式1和2的试样应力-应变关系曲线表现为应变硬化型，变形破坏形态为塑性剪切，其中排列方式2的试样具有相对大的抗剪强度。排列方式3是相对稳定的排列方式，随轴向载荷增加，试样呈现出较大的抗剪强度，但由于颗粒排列相对整齐一致，试样一旦发生破坏易形成贯通的破坏面，导致偏应力-轴向应变关系曲线出现应力降。排列方式1为相对不稳定的排列，载荷作用下结构容易失稳，试样具有较小的抗剪强度。

图7 不同颗粒排列方式下试样的数值剪切试验结果

4.4 结果讨论

岩土体具有空间变异性，室内试验采用试样的体积远小于岩土体地层的体积，相对均匀的同一地层，不同位置处试样的组成、结构及物理力学性质也可能存在差异(倪万魁等， 2001)，为尽量减小变异性，文中常规物理指标测试、CT扫描和三轴剪切试验选用取自同一位置的30cm×30cm×30cm的块状原状黄土试样。CT扫描试验的黄土试样尺寸为毫米级，试样颗粒的平均粒径约27.6μm，且文中借助Avizo软件分析了30多万个颗粒的结构特征，统计结果对于表征黄土试样颗粒的三维结构特征具有一定的代表性，而黄土结构空间变异性有待进一步研究探讨。陕西泾阳南塬地区黄土物理力学性质的相关研究成果很多，因此文中未对黄土试样的剪切性质展开详细论述，主要基于100kPa围压下试样的三轴剪切试验结果对建立的数值计算模型进行验证，并且试样的应力-应变关系曲线及抗剪强度等特征与该地区同时代同含水量原状黄土试样的三轴剪切试验结果基本一致(段钊等， 2016； 王烁， 2018)。

文中基于黄土颗粒结构特征和剪切试验结果建立的数值计算模型，探讨了颗粒尺寸分布、形状和排列对于剪切行为的影响，并对相应的微观机理进行了分析。研究结果也揭示了颗粒3种结构要素影响程度的差异，颗粒尺寸分布的影响主要表现在抗剪强度上，颗粒形状能够影响到试样的应力-应变关系曲线形态，颗粒排列的不同能够引起试样破坏机制的改变。但文中主要针对颗粒结构的各个要素对剪切行为的影响进行了初步的定性研究，量化各个结构要素的影响以及对宏观剪切行为的相对贡献是进一步研究的方向。此外，文中CT扫描采用的分辨率为1μm，所以主要关注了试样中粗颗粒(大于1μm)的结构特征及对剪切行为的影响，细颗粒(小于1μm)的结构需在纳米级分辨率下展开更小尺度的研究。黄土中的细颗粒主要为黏土矿物和盐分等，结构复杂且在黄土试样中主要以胶结物的形式存在，因此可将细颗粒的作用划分到颗粒接触特征中来表征，并深入探讨颗粒接触特征对黄土剪切行为的影响。

5 结 论

本文通过试验和数值模拟相结合的方法，研究了黄土的三维颗粒结构特征，分析了颗粒尺寸、颗粒形态和颗粒排列对黄土剪切行为的影响，得到如下结论：

(1)陕西泾阳Q2黄土试样的颗粒等效直径小于70μm，平均值为27.6μm，满足Weibull概率分布，颗粒球度分布符合Gamma分布。在球坐标系下，颗粒最大Feret直径方向角沿方位角方向的频率分布相对均匀，随极角的增大频率呈现增大的趋势， 60%颗粒的方向角集中分布在极角为60°～90°的方向，颗粒的排列表现出横观各向同性的特征。

(2)黄土颗粒的等效直径分布越均匀，试样的应变硬化现象越显著，抗剪强度越高，但试样变形破坏形态没有明显差异。颗粒球度值越小，试样的应变硬化现象越显著，但抗剪强度相对越低，并且球形颗粒的试样具有较低的抗剪强度值，试样变形破坏形态也没有明显差异。颗粒排列方式的影响表现为：颗粒方向集中分布在垂直于轴向加载方向的试样有较大的峰值强度和应力降，有贯通的剪切破坏面，颗粒方向集中分布在平行于轴向加载方向的试样具有较低的强度，应力-应变关系曲线表现为应变硬化型，变形破坏形态为塑性剪切。