DOI ：10.16030/jki.issn.1000-3665.202006027
各向异性对软土力学特性影响的离散元模拟
赵 洲1
，宋 晶
1,2,3
，刘锐鸿1 ，杨守颖1 ，李志杰
1
（1. 中山大学地球科学与工程学院，广东 广州 510275；2. 广东省地球动力作用与地质灾害重点实验
室，广东 广州 510275；3. 广东省地质过程与矿产资源探查重点实验室，广东 广州 510275）
摘要：软土预压工程中，初始和诱发各向异性对软土力学性质的影响十分显著，而现有研究缺乏对初始和诱发各向异性的统一研究方法。采用离散单元法，以颗粒长宽比作为定量评价指标，构建真实形态的颗粒模型，生成5组不同沉积角的初始各向异性试样，并进行竖直和水平两方向加载的双轴模拟实验，研 究了初始各向异性和诱发各向异性对软土力学特性影响；在细观层面，以颗粒为对象研究了颗粒接触形式和转动角度的变化规律，以接触为对象研究了配位数和接触法向各向异性的发展趋势，在此基础上探究抗剪强度指标与各向异性关系。结果表明：初始和诱发各向异性共同影响试样力学性质，当加载方向和软土沉积方向垂直时，土体有最大的峰值强度。颗粒接触形式中面面接触的比例随加载的进行逐渐增大，并影响着试样初始模量和抗剪强度，配位数和接触法向各向异性受颗粒接触形式的影响有不同的演化规律，并在加载后期趋于稳定；同时，初始各向异性试样相较各向同性试样有更大的黏聚力，诱发各向异性主要影响试样内摩擦角，进而影响试样抗剪强度。 关键词：软土；各向异性；峰值应力；抗剪强度指标；细宏观性质
中图分类号：TU411.3 文献标志码：A 文章编号：1000-3665（2021）02-0070-08
Discrete element simulation of the influence of anisotropy on the
rglt
mechanical properties of soft soil
ZHAO Zhou 1
，SONG Jing
1,2,3
，LIU Ruihong 1 ，YANG Shouying 1 ，LI Zhijie
1
（1. School of Earth Sciences and Engineering , Sun Yat-Sen University , Guangzhou , Guangdong 510275, China ；2. Guangdong Provincial Key Lab of Geodynamics and Geohazards , Guangzhou , Guangdong 510275, China ；
3. Guangdong Provincial Key Laboratory of Geological Processes and Mineral Resource Exploration ,
Guangzhou , Guangdong 510275, China ）
Abstract ：The initial and induced anisotropy has a significant effect on the mechanical properties of soft soil in preloading engineering. However, there is a lack of unified research methods for the initial and induced anisotropy. Discrete element method is adopted in this study, and the length-width ratio of particles is used as the quantitative evaluation index. Five types of initial anisotropy samples with different deposition angles are generated. The effects of initial anisotrop
y and induced anisotropy on the mechanical properties of soft soil are studied by vertical and horizontal loading. At the micro level, the contact form and rotation angle of particles are examined from the point of view of particles, and the development trend of coordination number and contact normal to anisotropy is studied from the point of view of contact. The relationship between shear strength index
收稿日期：2020-06-12；修订日期：2020-08-04
基金项目：国家自然科学基金项目（41877228；41402239；41877229）；广东省自然科学基金项目（2019A1515010554）；广州市科技计划项目
（201904010136）
第一作者：赵洲（1995-），男，硕士研究生，主要从事软土微观结构及数值模拟。E-mail ：*******************.edu 通讯作者：宋晶（1982-），女，博士，副教授，主要从事软土工程及微观结构研究。E-mail ：****************.edu
Vol. 48 No. 2水文地质工程地质
第 48 卷 第 2 期
Mar.，2021HYDROGEOLOGY & ENGINEERING GEOLOGY 2021 年 3 月
and anisotropy is explored. The results show that initial and induced anisotropy influences the mechanical properties of the samples. When the loading direction is perpendicular to the soft soil deposition direction, the soil mass has the maximum peak strength. In the particle contact form, the proportion of surface contact gradually increases with loading. The initial modulus and shear strength are influenced. The coordination number and contact normal anisotropy have different evolutions under the influence of particle contact form, and tend to be stable at the later loading stage. In terms of shear strength parameters, the initial anisotropy will increase the cohesion of the sample, and the induced anisotropy will mainly change the internal friction angle, which is consistent with the law of shear strength.
Keywords：soft soil；anisotropy；peak stress；shear strength index；mesoscopic and macroscopic characteristics
软土预压工程中，由于颗粒沉积环境与沉积条件的差异，促使颗粒形成不同方向的初始排列，从而造成结构上的各向异性，称为初始各向异性[1]。采用真空预压和堆载预压的处理方法，导致颗粒水平和竖直向的应力状态不同，进而影响土体力学行为，称为诱发各向异性[1]。因此，对软土力学性质的研究，需综合考虑初始和诱发各向异性的影响，建立细观组构与宏观力学特性之间的相互联系。
目前，研究各向异性的常用方法主要包括了土工试验和数值模拟。早期工作集中于研究各向异性土体对地基承载力的影响，Oda等[2]、Siddiquee等[3]通过室内试验发现，当荷载方向与颗粒沉积方向平行时，地基承载力会大幅提升。针对初始各向异性的研究，Zhang等[4]通过提取不同深度的天然黏土来获取各向异性试样，进行一系列固结试验，研究各向异性对压缩指数和屈服应力的影响，并从微观机理予以解释，但其理论有待完善。对于诱发各向异性的研究，Rouain等[5]在复杂应力路径的真三轴试验基础上，发现土体强度受应力路径的影响，并提出结构性土的统一本构模型。曾玲玲等[6]、刘恩龙等[7]、刘嘉英等[8]也针对土体在不同应力路径下强度、变形及破坏特性进行了探究，极大丰富了土体各向异性的基础理论。王清等[9 − 10]从微观角度研究了土体结构强度的形成机制，并指出结构体单元的动态平衡特征。对于土体颗粒内在结构信息的提取，目前常用有扫描电镜（SEM）、X射线断层成像（CT）、数字图像技术成像、热电导率和水电导率[11 − 13]等研究手段。室内试验和现场试验作为基础研究手段仍存在诸多局限，如试样的初始各向异性难以控制，制样方式和试验手段会对最终结果产生很大影响。更为关键的是，室内试验无法做到组构演化过程实时监测，从而难以通过细观机制解释宏观力学行为。
离散元法将散体介质独立为离散单元体，通过单元接触判断发生的相互作用，能较好反映散体材料细观结构方面特征。Ehsan[14]建立凸多边形的离散元颗粒，并模拟层理的不同初始角度，来探索初始各向异性对试样剪切和变形行为的影响。DAI等[15]采用三种不同的制样手段制备各向异性试样并进行三
轴剪切，分析了固有各向异性对剪切强度和体积膨胀率的影响机理，并使用组构张量对土体宏观响应进行表征。Wang等[16]通过二维和三维的DEM构建出4类沉积角试样，发现颗粒材料的胀缩性与颗粒各向异性和接触法向的演变过程有关。当前离散元研究方法多是将初始和诱发各向异性分开考虑，本文将二者结合考虑，研究初始和诱发各向异性对软土宏微观力学响应的共同影响。通过构建实际形状的软土颗粒模型，制备初始和诱发各向异性试样，进行双轴剪切试验，分析了各向异性对宏观力学性质的影响，并通过颗粒接触形式、配位数和接触法向各向异性等细观层面影响宏观力学差异的机理进行探讨，以期丰富细观土力学的理论。
1 颗粒形状统计与构建
室内试验土样取自深圳大鹏新区海域工程现场，土体呈深黑，属于高液限黏土[17]，其外形多呈片状结构，颗粒间的接触方式以边-面、面-面接触为主。在离散元程序中基本单位为标准的圆盘，无法真实反映颗粒的实际形状与接触方式。基于此，本文在考虑形状因素的基础上，建立真实颗粒以反应宏观力学性质，具体步骤如下：①采用PCAS图像处理软件对原状软黏土的SEM图像进行统计分析，如图1所示，并计算软土颗粒的形状参数；②对单个颗粒形态进行描述，选取长宽比As作为指标（定义见图2），反映颗粒轮廓的整体形态。绘制软土颗粒的长宽比累积分布
2021 年赵 洲，等：各向异性对软土力学特性影响的离散元模拟 · 71 ·
曲线，如图2所示；③基于颗粒实际形状，采用clump 模拟单颗粒模型，由多个球单元等间距排列为刚性集合体。在低围压下，不考虑颗粒内部的破碎，仅考虑clump 之间的相互作用，以及颗粒与墙体的接触，以模拟三轴试验中颗粒与橡胶膜的接触行为。
10 μm
(a) SEM 图
(b) 二值化图
图 1 软土SEM 图像Fig. 1 SEM image of soft soil
a b
A s =b /a
806040200
0.2
0.40.60.8 1.0
室内试验数值模拟
颗粒长宽比
图 2 颗粒形状和长宽比累积分布曲线
Fig. 2 Particle shape and cumulative distribution
curve of the aspect ratio
2 试样制备与加载
θθ制备试样具体方案如下：①在高100 mm 、宽50 mm 的矩形空间中，制备不同长宽比的颗粒，并按长宽比累积曲线分布于矩形空间中，初始孔隙率为0.2，与SEM 图片的面积孔隙率一致，颗粒粒径在4～8 mm 之间均匀分布，且满足边界效应。②制备初始各向异性试样，将颗粒长轴绕水平向旋转角后，限制颗粒旋转，并在一定时步内施加重力作用，此过程模拟沉积角为的土样原始沉积。③在矩形
空间高50 mm 处重新生成顶墙，形成边长为50 mm 的正方形试样，即为所需的加载试样，试样尺寸满足边界效应。卸除重力并解除颗粒旋转限制。④固结过程则向中心缓慢移动四周墙体，等向施加围压进行固结，达到所设定的初始围压并维持平衡。
固结完成后，对试样进行不同方向的加载以模拟诱发各向异性，加载方向设定为水平和竖直两个方
向，当进行水平向加载时，采用伺服机制保持上下端墙体应力恒定，左右端墙体以恒定速率相向运动，加载速率采用0.01 mm/s ，整个过程维持准静态。竖直加载方式则维持左右端墙体应力恒定，利用上下端墙体进行加载。初始围压设为100，200，400 kPa ，与室内试验一致。初始试样如图3所示。对模型局部放大，可反映出颗粒的3种接触方式：红代表面面接触，蓝代表边面接触，绿代表边边接触（图3）。
水平加载
θ
竖直加载
图 3 数值模型加载方式及颗粒接触方式
Fig. 3 Numerical model loading mode and particle contact mode
3 细观参数标定
为模拟软土颗粒之间的接触强度和黏结强度特性，选用线性接触黏结模型作为颗粒接触法则，并通过室内三轴剪切试验（固结不排水剪）进行参数标定，采用初始各向同性试样（iso ）的竖向加载进行离散元模拟，不断调整参数至模拟结果与室内试验的拟合程度达到最佳，数值模型的力学参数标定结果如表1所示。数值模拟和室内试验的应力应变曲线如图4所示。由图4可知，本文基于真实颗粒形态建立的数值试样与室内试验结果具有较好的一致性，在此基础上进行针对初始和诱发各向异性的三轴试验研究。
表 1 数值模型参数
Table 1 Parameters used in the numerical model
有效模量E /MPa
刚度比k 摩擦系数μ阻尼系数颗粒密度ρ/(kg·m −3
)
1
1
0.03
0.7
2 700
4 宏观力学结果
软土真空预压联合堆载工程中，主要存在2个方向的荷载，一类是竖向堆载，另一类为水平向吸力。数值模拟过程中，将考虑竖直方向和水平方向作为荷载方向，设置5组沉积角（0°～10°、20°～30°、40°～50°、60°～70°、80°～90°），用于模拟软土颗粒的5类沉积方向。开展三轴试验模拟，应力应变曲线如图5所示。
· 72 · 水文地质工程地质第 2 期
轴向应变 εq /%0°20°40°60°80°~90°
4020
10
20
轴向应变 εq /%604020
10
20(b) 水平加载
(a) 竖直加载
0°20°40°60°80°~90°
图 5 不同沉积角下偏应力与轴向应变关系Fig. 5 Relationship between the deviator stress and
axial strain of the different deposition angles
初始各向异性试样在两类加载方向下，偏应力随轴向应变的演化过程既具有相似规律，也存在显著差异之处。相似规律：加载初期偏应力随轴向应变的增大迅速增大，约在5%轴向应变处相对变缓；加载后期偏应力增大趋势逐渐平缓，达到峰值应力并最终维持恒定，曲线整体呈应变硬化型。显著差异：剑舞翩跹
竖直加载中（图5a ），颗粒沉积角为0°～10°时的峰值应力最大，而
水平加载中（图5b ），沉积角为80°～90 °时的峰值应力最大。即偏应力大小不仅限于受加载类型的影响，而且取决于颗粒沉积方向与加载方向的关系。当颗粒沉积方向与加载方向垂直时，试样峰值应力可达到最大值，即试样表现更大的抗剪强度。可见，软土宏观力学性质受初始各向异性和诱发各向异性的影响显著。当加载方向相同时，颗粒沉积方向与加载方向垂直的试样抗剪强度更高，而同一沉积角试样在不同方向下加载，达到的抗剪强度也不同。软土固结过程中，若全程采用真空预压，含量较高的黏土颗粒沉积角较大，土样抗剪强度不能有效增加。因此，软土工程后期多采用真空联合堆载，堆载的竖直方向与颗粒沉积角更接近90°，可以进一步提高软土抗剪强度。
图6为当轴向应变为15%时对应的试样偏应力为峰值应力统计值。当试样为横向沉积0°～10°和轴向沉积80°～90°时，加载方向对峰值应力影响最显著，而当沉积角为40°～50°时影响最小，与各向同性试样最为接近。当加载方向相同时，相对于各向同性（iso ）试样峰值应力变化原因为试样的初始各向异性，沉积角与加载方向垂直时，所产生的峰值应力会远远大于各向同性试样。而沉积角与加载方向一致时，峰
值应力则小于各向同性试样。软土固结过程前，充分搅拌或者沉积的工法更有利于形成各向异性的土体，也更利于固结后抗剪强度增加。
0~10
20~3040~5060~7080~90
iso
沉积角/(°)
差值
陈和平
竖直加载水平加载图 6 峰值应力及差值
Fig. 6 Peak stress and difference
各向异性软土的力学特性与微观组构的演化过程密切相关，朱楠等
[18]
通过室内微观试验对比不同地
区软土工程特性，表明颗粒接触方式的不同是软土工程性质差异的主要原因。在此基础上，本文针对
三轴试验中颗粒接触方式的演化过程进一步研究，发现在单一周围压力加载条件下，接触方式对力学性质的影响不仅体现在某一临界状态，而且体现在加载的全过程。长镜头理论
1201008060402000
5
10
15
轴向应变/%
试验数值
=400 kPa 3=200 kPa σ3=100 kPa
图 4 数值试验与室内试验应力应变曲线Fig. 4 Stress and strain curves of the
numerical and laboratory tests
2021 年
赵 洲，等：各向异性对软土力学特性影响的离散元模拟
· 73 ·
已从实验模拟和工程实践角度分析了沉积角和加载方向对宏观力学特征的影响，为体现剪切过程颗粒接触形式的变化特征，下面从细观角度进行分析。
5 细观机理分析
5.1 颗粒面面接触及颗粒转角
对颗粒长轴a 、短轴b 及枝向量l 进行统计，本次模拟中颗粒接触类型包括边边接触、边面接触和面面接触，其判断方式如图7所示。面面接触系数P 可定义为面面接触数占总接触数的比例，即为：
式中，C 1−
面面接触数；
C 0−
初始面面接触数；C −总接触数。
(b) 边面接触
(c) 面面接触
变形计爱在远山全集(a) 边边接触
a 1/2+a 2/2=l a 1/2+
b 2/2=l b 1/2+b 2/2=l 图 7 颗粒接触形式及判别方法
Fig. 7 Particle contact form and the discriminant method
基于离散元模拟结果，图8显示了面面接触系数P 随轴向应变的演化过程。在剪切初期，颗粒通过平移、旋转等方式调整自身位置，大量接触形式由边
边、边面接触转化为面面接触，因此颗粒接触形式中的面面接触形式迅速增加，并在偏应变约5%时达到峰值，此变化趋势与试样应力应变具有一致性，说明剪切初期面面接触形式能有效增加试样初始模量。剪切后期，颗粒接触基本稳定，面面接触系数受初始各向异性影响较大，随轴向应变呈现增长趋势，并在剪切末期逐渐稳定。
在加载过程中，不同初始各向异性对应的面面接触系数曲线形状相似但具有不同的初始斜率和峰值状态，且与应力应变曲线形态一致，说明面面接触对力学性质的影响不仅作用在初始状态[18]
，也作用于整
个加载过程。
基于离散元试验结果，对不同形状的颗粒达到峰值应力时的转角进行统计，以40°～50°沉积角试样为例，结果如图9所示。在二维平面中颗粒长宽比越小，形状更为细长，颗粒转动惯量越大，在加载过程中受阻力越大而难以发生偏转，这与Marschall 等[19]
、袁
斌等
[20]
研究结论一致。因此初始各向异性试样在二
维加载中，颗粒发生转动从而改变颗粒的接触形式，最终影响峰值应力。
5.2 配位数演化规律
离散元模拟过程中，可参考配位数反映试样在加载过程中的密实程度。配位数越高，颗粒之间的相互连接越紧密。配位数定义为
[21]
：
2
0°20°40°60°80°轴向应变 εq /%0
10
20
2
0°20°40°60°80°~
90°
轴向应变 εq /%0
vb1210
20
图 8 不同沉积角下面面接触系数与轴向应变关系Fig. 8 The relationship between surface contact coefficient and
axial strain of the different diposition angles
图 9 沉积角40°～50°试样颗粒长宽比与转角关系Fig. 9 Relationship between the aspect ratio and angle of the
sample with diposition angle 40°-50°
· 74 · 水文地质工程地质
第 2 期
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议