黄土单向固结试验微观非连续变形分析

郭龙骁，张常亮，杨德广，李同录

(长安大学 地质工程与测绘学院， 西安 710054)

黄土单向固结试验微观非连续变形分析

郭龙骁，张常亮，杨德广，李同录

(长安大学 地质工程与测绘学院， 西安 710054)

黄土骨架颗粒的架空结构体系控制着其主要力学行为，应用非连续变形分析方法可从微观揭示其力学行为的本质。采取陕西省泾阳南塬的L5(Q2)黄土原状样，进行各种土力学试验，测定其物理力学性质参数；采用光学显微镜和图像处理软件得到土样的微观结构图像，描绘颗粒形态，建立该土样的非连续变形分析微观结构模型，并模拟了土样在1.0，12.5，25.0，50.0，100.0，200.0，400.0 kPa荷载作用下的压缩变形。结果得到：土样在压缩过程中，微观结构变化图像及压缩曲线和室内固结试验结果较为接近，说明采用非连续变形分析方法模拟黄土的微观变形特性是可行的。研究结果为分析黄土的微观物理力学性质提供了一种思路。

黄土；骨架颗粒的架空结构；固结试验；微观结构；非连续变形分析

1 研究背景

黄土是一种典型的结构性土，由骨架颗粒的点接触和面胶结构成的架空结构体系控制着其主要力学行为。黄土宏观上的变形和破坏表现为微观上的局部破坏和连通性破坏。从微观结构入手可揭示其宏观变形和破坏等力学行为的本质。离散单元法为土体微观研究提供了一种新思路，目前广泛采用其模拟砂土和黏性土的微观结构，其中颗粒流[1-3](Particle Flow Code)和非连续变形分析方法是较常用的离散单元法。

近年来，在利用颗粒流模拟黄土变形特性方面有一些开拓性工作，蒋明镜等[4-5]针对黄土大孔隙和胶结特性，应用离散元方法生成不同含水率的黄土试样，研究试样的一维湿陷特性，其结果为认识黄土复杂力学特性和建立其本构理论提供了基础。但是该方法有其不足之处，颗粒流假定单元颗粒为圆形且为刚体，颗粒形态单一，与黄土颗粒的实际形态相差大，只有点接触，在初始模型中颗粒容易滑落而形成密实的结构，难以形成黄土特有的架空结构。

非连续变形分析方法(Discontinuous Deformation Analysis，简称DDA)由石根华博士等[6]于20世纪80年代提出，可模拟散体系统，允许单元自身有弹性形变，单元形状可以是任意凸形或凹形。该方法主要用于岩石块体分析，在边坡及地下工程中得到发展和应用[7-9]，一些学者试图将其扩展应用到土力学领域，如张国新等[10]用DDA分析了砂土在平面应变条件下变形特性的微观机理。其数值试验将土的颗粒简化成由横断面为11条边构成的等边棱柱体，与颗粒流的模型类似，没有考虑颗粒的实际形态。郭培玺等[11]应用DDA模拟重力作用下粗黏料在承模筒内的下落过程并以此建立粗黏料数值模拟试验模型，说明DDA数值模拟应用于粗粒料力学特性的研究是可行的。其采用满足级配要求的任意多边形颗粒作为试验模型，不足之处在于只考虑了颗粒级配和试验结果一致，没有考虑颗粒的实际形态。

本文采取陕西省泾阳泾河南岸L5(Q2)黄土原状样，进行基本物理力学性质测试、固结和三轴固结不排水试验，测定了其基本物理力学性质指标，将其用于模型参数的标定和数值模拟结果的对比。用光学显微镜和图像处理软件获取其微观结构图像，描绘其颗粒形态，建立该土样的DDA微观结构模型，模拟土样在各级荷载作用下的变形，得到土体微观结构动态变化过程和压缩曲线，并与单向固结试验结果进行对比分析。

2 L5黄土试样的物理力学特性

图1 试样的粒径分布曲线Fig.1 Curves of particle size distribution of specimen

图2 试样的e-p曲线Fig.2 Curve of e-p of the specimen

对土样进行室内试验，测得其基本物理力学参数如下：相对密度为2.66；重度为17.7 kN/m3；天然含水率11.3%；饱和度46.5%；孔隙比0.687；塑限19.3%(搓条法)；液限28.8%(碟式液限仪法)；塑性指数9.5。其粒度曲线如图1，可见该土样粒度成分以粉粒为主，其次为黏粒。土较致密，不具有湿陷性。 将原状土样制成环刀样，做饱和土样的单向固结试验，竖向荷载按1.0，12.5，25.0，50.0，100.0，200.0，300.0，400.0，500.0，600.0，700.0，800.0 kPa逐级增加，测定土样在各级竖向荷载作用下产生的变形，计算各级荷载下相应的孔隙比，得到图2所示的e-p曲线。

同时对该土样还进行了饱和状态的三轴固结排水试验。图3(a)为固结围压分别在100～500 kPa时的应力应变关系曲线，可知土样在>300 kPa时呈应变硬化特点，在<300 kPa时呈理想弹塑性特点；图3(b)为试验所得总应力路径和有效应力路径曲线，2种曲线接近，孔压较小。取轴向应变达到15%时的剪应力为破坏强度，在图3(b)中可得到土体的破坏主应力迹线Kf′, 由其斜率b和倾角α，通过式(1)和式(2)，可计算得到土体的黏聚力为4.0 kPa，内摩擦角12.8°。

(1)

(2)

3 微结构图像提取及其DDA模型

利用自主改装的微观结构成像系统采集黄土微观结构的信息，如图4所示，该系统由光学显微镜、环形无影灯、数码相机、数据线、图像处理软件组成。利用该系统可将置于光学显微镜下的黄土试样的微观结构图通过数码相机和数据线直接导入图像处理系统中并进行存储。

图4 黄土微观结构成像系统Fig.4 Imaging system of loess microstructure

图像获取前，需先将黄土试样切割成1.5 cm×1.5 cm×4 cm的柱状体，并在中间部位刻一道1 mm深的浅槽，将试样从中间掰开，尽量使断面平整，不扰动其掰开后的原始状态(图5)，将断面置于显微镜下观察，显微镜放大倍率最大可达500倍，能够清晰、直观地观察研究黄土颗粒组成、孔隙特征。图6是通过黄土微观结构成像系统拍摄的全景深图，基于光学显微镜所获取的黄土微观图片建立DDA分析模型。

图5 黄土试样制备及微结构观察断面Fig.5 Preparation of loess specimen and a section for microstructure observation

图6 黄土微观结构图像Fig.6 Microstructure image of loess specimen

取图6白色框内土体作为分析对象，应用图像处理软件CAD提取土颗粒几何数据，作为DDA分析的单元，可生成图7所示的颗粒几何模型图。

图7 DDA前处理子程序建立的颗粒模型Fig.7 DDA particle model built with pretreatment subroutine

DDA以单元位移和应变作为基本未知量，利用势能最小原理将单元刚度、质量、荷载、惯性力等子矩阵，以及单元间接触和约束子矩阵统一成一个整体系统矩阵，建立系统的平衡方程组，引入边界条件求解。该方法中，每个块体是一个弹性单元，各块体单元可以发生弹性变形，可以具有不同的弹性模量和泊松比，但每个单元是常应变和常应力的连续介质，这一点类似于有限元；而块体之间为弹性接触，可以彼此脱离或自由运动，这一点类似于离散元。DDA通过求解隐式的方程组获得各单元的位移和应变，进一步求解应力，求解方法类似有限元。这些特点都使得将该方法可用到具有任意形状颗粒的黄土微观结构分析中。

参考固结试验设备，模型左、下、右3个面施加固定约束板作为土样盒，顶部设置可移动加载板，在加载板上施加20个均匀点荷载模拟固结试验的压力，同时设置2个观测点用于竖向变形测量，在土样盒的四周设置8个固定点限制其线位移和角位移。试验用的土样盒及压力盖都是DDA的块体单元，土样盒和加载板内部的弹性模量为土颗粒的20倍，近于刚体，结合室内试验在环刀内部涂凡士林的方法，对土样盒内壁做光滑处理，因此将其内壁的摩擦角设为0°，黏聚力设为0 kPa。由于该数值模拟采用准静态模型，边界条件对试验结果影响不大。

模型主要几何参数与物理参数取值为：模型规范为0.60 mm×0.40 mm，其中土样为0.28 mm×0.60 mm；土颗粒单元总数为604；土颗粒的总面积为9.98×10-2mm2；孔隙面积为6.82×10-2mm2；模型初始孔隙比为0.683，与实测的土样孔隙比0.686较为接近。

在微观模拟中，土颗粒的密度取值不包括孔隙部分，可取其相对密度值。根据对该土样的粒度分析结果(图1)，其粉粒及粉砂粒质量百分比占72%，黏粒占28%。黄土中粉粒及砂粒是其骨架，黏粒存在于孔隙中或空隙间。黄土的骨架颗粒主要由长石和石英构成，因此土粒的弹性模量和泊松比参考完整岩石(如花岗岩)的取值。根据三轴试验测得该土样的有效内摩擦角为12.8°，黏聚力为4.0 kPa。由于摩擦系数是无量纲的数，摩擦力与法向接触力的作用点面积无关，因此颗粒间的摩擦角直接取实测的土的摩擦角。黏聚力与法向接触力无关，但是黏聚力是一个与接触面积成正比的物理量，可将三轴试验得到的黏聚力换算成模型的接触黏聚力，根据单位面积接触点的数量换算出粒间接触黏聚力为0.8 Pa。

DDA物理模型需提供颗粒物理力学参数如表1所示。

表1 模型物理力学参数Table 1 Physical and mechanical parameters of the model

DDA的块体内部采用线弹性本构关系，不考虑颗粒的破碎，块体间的接触面上采用摩尔-库伦准则，接触有3种状态：张开、滑动和锁定。当块体相互接触后，根据其受力的大小和方向判断施加或去掉切向弹簧或法向弹簧。DDA可以通过测量点的数据来判断模型整体是否稳定，在本数值模拟中如果2个测量点的位移在一定时间内变化不大，即可认为土样在该级荷载下稳定；然后施加下一级荷载，荷载分为12.5，25.0，50.0，100.0，200.0，400.0 kPa共6级。为了保证土体颗粒与压力盖接触良好，在施加第1级荷载之前先施加了1.0 kPa的预应力，以施加预应力稳定后的土体作为模型初始状态，然后施加其它6级荷载稳定。

4 模拟结果

由于荷载增量小时，变形不容易观察到，因此提取其中变形明显的4幅图像，如图8所示。

图8 各级荷载稳定后的模型Fig.8 Images of the model under different levels of loadings at stable state

可以看出，在低应力状态下(25 kPa)由骨架颗粒构成的架空结构很少被破坏，架空结构存在且稳定，能够承受一定的压力。而随着压力升高(100～400 kPa)，架空结构逐渐被破坏，粒间孔隙逐渐被压密，土体的孔隙比减小且以粒间孔隙为主。从微观上反映出黄土固结过程中的结构变化。

模拟微结构微缩模型的单向固结过程，相应变形微小，荷载作用下的总变形量只有0.015 mm左右，占总高度的5%(图8)，仅凭肉眼观察，各模型图之间的变化似乎不是很明显，但是对比首尾模型图可以很明显地看出孔隙被压密。提取2个测量点的位移，用其平均值作为压缩量，得出图9所示的压缩量与迭代次数关系曲线。由图9可见，在各级应力下，固结变形随着荷载的增加而增加；在某一级荷载上，随着迭代次数增长，变形量变化不明显，说明其收敛很快，当收敛于某一值后，压缩位移基本保持不变。利用各级荷载作用下产生的变形计算对应的孔隙比，可得到图10所示的e-p曲线。

图9 迭代次数-压缩量曲线Fig.9 Curve of compression deformation vs. iteration number

图10 DDA模拟与室内试验对比e-p曲线Fig.10 Comparison of e-p curves between DDA simulation and laboratory test

与单向固结试验结果相比可见，2条曲线较为相似，孔隙比相差≤0.01，但是随着压力增加，DDA模拟得到的孔隙比与实际试验结果逐渐产生了偏差，在400.0 kPa时DDA模拟得到的孔隙比是0.600，而试验情况是0.617，在400.0 kPa以内两者的误差<3%。

5 结论与讨论

本文采用黄土微观结构成像系统提取L5黄土原状样的微观图像，描绘黄土的颗粒和孔隙，真实反映黄土颗粒组成及各类孔隙的特点，依此建立的DDA黄土颗粒模型相较于理想模型更能体现黄土的结构特性。通过原状黄土单向固结模拟试验前后的对比，可以看出压缩过程中土颗粒发生运动，黄土架空结构被破坏，粒间孔隙被压密。模拟得到的e-p曲线与室内试验基本吻合，试验结果从微观上说明土体的压缩是由于受压后土体架空结构被破坏，部分土颗粒进入大孔隙，部分粒间孔隙闭合，总体积减小所致。

模拟试验的结果与固结试验的结果存在一定的偏差，主要原因：①受其计算能力的限制，无法做到与室内试验一样的比例模型，存在尺寸效应；②参数标定有可能存在偏差，宏观变形和强度参数用于微观情况需要进一步研究；③由于仅考虑了土、气两相介质，土体中的胶结、接触-胶结等连接方式不易体现，凝块等颗粒形态也不易表现。将毛细水对土颗粒的作用引入DDA程序中，从三相角度考虑黄土微观水土作用有待于进一步研究。

总之，本文只是提供了利用DDA研究黄土力学行为的一个简单思路，大量具体问题还需要深入研究。

[1] 周 健，池毓蔚，池 永，等.砂土双轴试验的颗粒流模拟[J].岩土工程学报，2000，22(6)：701-704.

[2] 周 健, 池 永. 颗粒流方法及 PFC2D 程序[J]. 岩土力学, 2000, 21(3): 271-274.

[3] 廖雄华, 周 健, 徐建平, 等. 粘性土室内平面应变试验的颗粒流模拟[J]. 水利学报, 2002, (12): 11-17.

[4] 蒋明镜，胡海军，彭建兵. 结构性黄土一维湿陷特性的离散元数值模拟[J].岩土力学,2013,34(4):1121-1130.

[5] 蒋明镜，李 涛，胡海军. 结构性黄土双轴压缩试验的离散元数值仿真分析[J].岩土工程学报,2013,35(增2):242-246.

[6] SHI G H, GOODMAN R E. Two Dimensional Discontinuous Deformation analysis[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1985, 9(6): 541-556.

[7] 邬爱清, 丁秀丽, 卢 波, 等. DDA 方法块体稳定性验证及其在岩质边坡稳定性分析中的应用[J]. 岩石力学与工程学报, 2008, 27(4): 664-672.

[8] 邬爱清，丁秀丽，李会中，等. 非连续变形分析方法模拟千将坪滑坡启动与滑坡全过程[J]. 岩石力学与工程学报,2006,25(7):1297-1303.

[9] 王书法，李树忱，李术才，等. 节理岩质边坡变形的DDA模拟[J]. 岩土力学,2002,23(3):352-354.

[10]张国新，李广信，郭瑞平. 不连续变形分析与土的应力应变关系[J]. 清华大学学报(自然科学版),2000,40(8):102-105.

[11]郭培玺，林绍忠. 粗黏料力学特性的DDA数值模拟[J]. 长江科学院院报，2008，25(1):58-60，69.

(编辑：姜小兰)

Discontinuous Deformation Analysis with Microstructural Modelfor Axial Odometer Test of Loess

GUO Long-xiao,ZHANG Chang-liang,YANG De-guang,LI Tong-lu

(School of Geology Engineering and Geomatics, Chang’an University, Xi’an 710054， China)

Mechanical behavior of loess is controlled by overhead structure of skeleton particle, and microstructure of loess could reveal its mechanical behavior by Discontinuous Deformation Analysis(DDA). In order to simulate the characteristics of consolidation of microscopic loess, we carried out test on undisturbed L5(Q2) loess specimens from southern Jingyang plateau of Shaanxi province, and measured the physical and mechanical parameters of loess. Then we obtained microstructure image and particle morphology of the specimens by using optical microscope and image processing software. On this basis, we established soil microstructure model and simulated the compression deformation under load of 1.0, 12.5, 25.0, 50.0, 100.0, 200.0, 100.0 kPa, respectively. Finally, we obtained microstructure image and compression curve of the model in the compression process, which are close to the results of indoor test. Results show that DDA is feasible for simulating the characteristics of the microscopic deformation of loess and this paper provides an idea for analyzing the microscopic physical and mechanical properties of loess.

loess; overhead structure of skeleton particle; consolidation tests; microscopic structure; discontinuous deformation analysis

2015-12-23；

2016-05-09

陕西省自然科学基础研究计划项目(2016JM4018)；中央高校基本科研业务专项基金(310826161022)

郭龙骁(1990-)，男，陕西西安人，硕士研究生，主要从事黄土微观模拟方面的研究，(电话)13323138969(电子信箱)fudenun@163.com。

10.11988/ckyyb.20151103

2017,34(3):80-84

TU444

A

1001-5485(2017)03-0080-05