非饱和重塑黄土的一种简化离散元模型

郭鸿，黄文华，李军，魏鸿林

（1.陕西理工大学土木工程与建筑学院，陕西汉中723000；

2.西北农林科技大学水利水电科学研究院，陕西杨凌712100；

3.中国城市建设研究院有限公司，北京100120）

非饱和重塑黄土的一种简化离散元模型

郭鸿1，2，黄文华1，李军1，魏鸿林3

（1.陕西理工大学土木工程与建筑学院，陕西汉中723000；

2.西北农林科技大学水利水电科学研究院，陕西杨凌712100；

3.中国城市建设研究院有限公司，北京100120）

黄土是一种水敏性颗粒材料，其强度受含水率影响十分明显。基于非饱和土力学的基本理论，利用颗粒离散元方法，建立了一种非饱和黄土的简化离散元模型。在对非饱和黄土的力学特性进行离散元模拟分析后，与室内三轴试验进行了对比验证。对比模型和试验结果发现：不同含水率试样的应力应变关系受土体颗粒间摩擦系数直接影响，表现为试样含水率越高，摩擦系数越小，说明该离散元模型在一定范围内能够反映重塑黄土的土力学特性。

非饱和黄土；离散元模型；含水率；三轴试验

1 研究背景

黄土的分布十分广泛，在中国分布面积达64万km2。对黄土的力学特性研究非常多，尤其以室内试验为主。但是就微观层面的研究，尤其是基于微观结构的本构模型研究，还需要做很多工作。近年来，离散元模型越来越受到广大科研工作者的青睐，而且其在黄土力学中作用也越来越大［1-2］。

离散元法（Discrete Element Method，DEM）［3］是将每一个颗粒或块体当作一个单元，观察记录数值试验过程中单元之间相互作用与运动规律，来预测散体群行为的方法。它是一种能准确、实时获取数值试验过程中试样所有颗粒、应力、孔隙等微观信息的工具［4-6］。蒋明镜等［7-9］在离散元模拟黄土方面做了大量的研究。Lu Ning等［10-11］基于吸力应力理念对非饱和土做了全新的阐释。然而，他们的方法由于其复杂的微观参数，在工程实践中推广具有一定的局限性。

离散元所建立的土颗粒之间的胶结参数诸如土颗粒之间的黏结力，无法和实际土体相对应。并且，黄土颗粒之间由于不同颗粒之间的胶结情况随着土颗粒的形状以及胶结盐的分布等变化非常复杂。因此，笔者认为，黄土剪切破坏的根本原因并不在于颗粒之间胶结由于“拉伸”引起微观破坏，进而引起宏观的破坏；而在于胶结情况发生变化，颗粒之间的切向接触发生了变化，即不同含水率的黄土颗粒，其颗粒在相同法向应力的情况下，滑动的难易程度不同。基于此理念，本文提出了非饱和黄土的一种基于含水率为变量的切向滑动离散元模型。该模型运用最简单的线弹性模型，具有参数简单，预测效果理想等特点。

2 颗粒离散元模型

非饱和黄土由于其内部存在胶结物质，而且其具有很强的水敏性，水的作用使非饱和黄土的强度大大减小。非饱和土颗粒之间由于水表面张力的存在，颗粒间产生基质吸力，其力学示意图如图1所示。颗粒间的水存在两个曲面，一个曲面半径为r1，另一个曲面半径为r2，颗粒半径为R。

如图1（a）所示，根据受力平衡分析［11］，得到基质吸力的值为：

其中Ts为水的表面张力。假设水膜圆弧和颗粒的切点与颗粒之间接触点形成的圆弧夹角为θ（如图1（b）所示），则由几何关系可以推导出：

将（2）代入（1）中，并令

得：

容易看出，f(θ)是关于θ的减函数，这也就意味着Su也是θ的减函数。

图1 非饱和土受力分析［11］

假设沿法向接触的面积为An，沿切向接触的面积为As，则颗粒法向和切向的黏结力分别为：

考虑到以上因素，本文基于Cundall［3］提出的二维圆盘模型，用颗粒流程序PFC2D（Particle Flow Code in Two Dimensions）对非饱和黄土颗粒进行模拟，其简易离散元模型如图2所示。

图2 颗粒离散元接触模型

不考虑水的作用，法向和切向的接触力为Fn=Knδn和Fs=Ksδs，式中Kn和Ks为分别为接触颗粒的法向刚度，δn和δs分别为颗粒间法向和切向的相对位移。最大切向力表示为：

式中μ0为颗粒间的摩擦系数。

然而由于基质吸力的存在，对于颗粒切向的最大接触力则为：

令μ=μ0+Ps/δKn=μ(θ)，则式（7）变为：

由于饱和度Sr和接触角θ的关系可以表示为［1］：

东西如此，人也如此。一个人好与不好，不在模样上。俊朗的人不一定是好人，丑陋的人也不一定是坏人。这么说吧，高人不一定个高，高人是境界高，与身量无关。人附加值在哪里呢？在心。一个人，心好，就是好人，好人好气息。心坏，就是坏人，坏人坏气息。所以人也分贵贱，高人则贵，滥人则贱。君子则贵，小人则贱。

故，新的摩擦系数可以用含水率表示为：

式（8）则变为：

前人的研究也发现，在一定范围内，内摩擦角与试样含水率由存在着特定关系［12-14］。因此，笔者建立简化模型中主要考虑参数Kn、Ks和μ即可，其中Kn和Ks与土颗粒的刚度（软硬）有关，而μ则与含水率w相对应。

3 非饱和重塑黄土的三轴试验

3.1 试验材料与方案本文采用陕西泾阳地区黄土作为试验原料土，取土深度在3.0～4.0 m之间，测得天然含水率16.2%，天然干密度1.28 g/cm3，孔隙比1.17，为Q3黄土，试验基本物理性质指标见表1。

表1 黄土基本物理指标

试样采用非饱和重塑试样，为原料土风干、碾碎、过筛后加水压制而成。试样直径为39.1 mm，高为80mm，控制试样密度为天然干密度，含水率为11%、18%和26%。三轴试验采用SBL-1型应力应变控制式三轴剪切渗透试验仪，控制固结围压为100 kPa、200 kPa和300 kPa作不固结不排水（UU）三轴剪切试验，试验方案如表2所示。

表2 三轴剪切试验方案

3.2 试验结果剪切试验采用应变控制，试样装样完成后，在轴向加载系统中设定好轴向应变速率，在计算机中设置好应力、应变记录时间间隔，最后进行剪切试验。当试样轴向形变大于15%，或试样出现明显剪切带，应力值急剧下降，即可停止试验，认为试样达到破坏。

图3 三轴试验应力应变曲线

4 模拟方法及结果

本次研究基于PC环境下的颗粒流程序PFC2D，用笔者建立的离散元模型模拟二维黄土的三轴剪切试验。如图4为实际三轴试验试样的颗粒级配曲线和离散元模型所选定级配曲线。表3为模拟选定参数和试验土样实际参数。需要说明的是离散元数值试验的颗粒刚度（法向和切向）是根据实际三轴试验标定的，而墙体的刚度则取值为颗粒刚度的十分之一，模拟三轴试验的橡皮膜。

表3 离散元模型模拟参数和试验参数

图4 试样颗粒级配曲线和模型选定级配曲线

离散元模拟三轴试验时，首先根据表3的参数生成离散元数值试样，其中颗粒的接触关系满足式（11），然后用fish语言编写的围压伺服程序通过周围墙体（模拟橡皮膜）以及上下墙体（模拟上下加压板）对试样施加固结围压σ3，等循坏稳定后，通过向下移动（速度和实际三轴试验保持一致）试样上部墙体施加轴向应力σ1。最后通过程序输出σ1-σ3与ε（竖向位移与试样初始高度的比值）的关系。

如图5所示为PFC2D根据模型参数选定值模拟出来的模拟曲线与实际三轴试验的试验曲线进行的比较分析，图例中短线前面的数字代表含水率，短线后面的数字为围压，sim开头的为模拟值，其余为试验值。

从图5中可以看出，笔者建立离散元模型模拟得到的应力应变关系曲线与试验实际应力应变关系曲线能较好吻合，模拟精度在高含水率和高围压条件下尤为突出。进一步根据室内三轴试验（图3）和离散元数值模拟结果（图5）绘制莫尔圆，通过最小二乘法拟合得到内摩擦角的试验值和模拟值，如表4所示。可以看出，虽然试验得到的内摩擦角随含水率增大存在一点起伏，而模型模拟得到的内摩擦角随着含水率在11%至26%间增大而减小，但是两者总体趋势一致，且相同含水率下试验值与模拟值差值较小。

笔者建立简化颗粒离散元模型，针对不同含水率的重塑黄土，只需调整一个参数，即颗粒间的摩擦系数μ，就可以较为理想地反应土体的应力应变关系。本文的研究表明该离散元模型在一定范围内能够反映陕西泾阳地区非饱和重塑黄土的力学性质。从本文非饱和土力学的观点看，含水率较小时，θ减小，基质吸力Su增大，颗粒之间黏结的更紧密，更不容易被剪切错动。因此，含水率较低的黄土剪切强度较高，颗粒之间的摩擦系数和内摩擦角较高。

图5 不同含水率、围压下模拟曲线和试验曲线

表4 试验与模拟得到内摩擦角比较

5 结论

（1）本文简化离散元模型对非饱和重塑黄土的三轴剪切试验进行数值模拟，对比两者结果发现，模型所得应力应变关系曲线与实际试验结果吻合良好，曲线表现出的趋势一致，且在高含水率和高围压下模拟效果较为理想。

（2）不同含水率下，试验得出土样内摩擦角与模型模拟得出内摩擦角大小相近，且随着含水率增大，内摩擦角总体上呈减小趋势。

（3）笔者建立简化颗粒离散元模型，针对不同含水率的重塑黄土，只需调整一个参数，即颗粒间的摩擦系数，就可以较为理想地反应土体的应力应变关系。

［1］Jiang M J，Leroueil S，Konrad J M.Insight into shear strength functions of unsaturated granulates by DEM analy⁃sis［J］.Computers and Geotechnics，2004，31（6）：473-489.

［2］蒋明镜，胡海军.密实和松散颗粒材料等吸力三轴剪切试验离散元数值模拟［J］.中南大学学报：自然科学版，2010，41（6）：2350-2359.

［3］Cundall P A，Strack O D L.The distinct element method as a tool for research in granular media［R］.Minnesota：University of Minnesota，1978.

［4］李凡.岩土材料破损特性的颗粒流研究［J］.土木工程材料，2007，40（9）：78-81.

［5］Cheng Y P，Boltln M D，Nakata Y.Crushing and plastic deformation of soil simulated using DEM［J］.Geotech⁃nique，2004，54（2）：131-141.

［6］Bolton M D，Nakaat Y，Cheng Y P.Micro-and macro-mechanical behavior of DEM crushable materials［J］.Geo⁃technique，2008，58（6）：471-480.

［7］Jiang M J，Konrad J M，Lreoueil S.An efficient technique for generating homogeneous specimens for DEM studies［J］.Computers and Geotechnics，2003，30（7）：579-597.

［8］Jiang M J，Yu H S，Leroueil S.A simple and efficient approach to capturing bonding effect in naturally micro⁃structured sands by discrete element method［J］.International Journal for Numerical Methods in Engineering，2007，69（6）：1158-1193.

［9］蒋明镜，沈珠江，赵魁芝，等.结构性黄土湿陷性指标室内测定方法的探讨［J］.水利水运科学研究，1999，21（1）：65-71.

［10］Lu N，Lechman J，Miler K T.Experimental verification of capillary force and water retention between unevensized spheres［J］.Journal of Engineering Mechanics，2008，134（5）：385-395.

［11］Lu N，Likos W J.Unsaturated soil mechanics［M］.New York：John Wiley，2004.

［12］王浩林，白晓红，冯俊琴.压实黄土状填土抗剪强度指标的影响因素探讨［J］.岩土工程学报，2010，32（增2）：131-134.

［13］贾亮，朱彦鹏，朱鋆川.兰州马兰、离石压实黄土抗剪强度影响因素探讨［J］.岩土工程学报，2014，36（增2）：120-124.

［14］王铁行，罗少锋，刘小军.考虑含水率影响的非饱和原状黄土冻融强度试验研究［J］.岩土力学，2010，31（8）：2378-2382.

A simple discrete element model of unsaturated remodeling loess

GUO Hong1，2，HUANG Wenhua2，LI Jun1，WEI Honglin3
（1.School of Civil Engineering and Architecture，Shaanxi Sci-Tech University，Hanzhong723000，China；
2.Institute of Water Resources and Hydropower Research，Northwest A&F University，Yanglilng712100，China；
3.China Academy of Urban Construction Co.，Ltd.，Beijing100120，China）

Loess is a kind of water sensitive granular material,and its strength is obviously affected by wa⁃ter content.Based on the basic theory of unsaturated soil mechanics,a simplified discrete element model of unsaturated loess is established by using the discrete element method.After analyzing the mechanical charac⁃teristics of unsaturated loess,the triaxial test was carried out for verification.After comparison of the re⁃sults from the model and experiments，it is found that the inter-particle friction coefficient influences the strain-stress relationship with different water contents.In other words,samples with higher water content could be simulated by lower internal friction coefficient in DEM.This，to some extent，shows that the DEM model in this paper could reflect the soil mechanics properties of unsaturated remodeling loess.

unsaturated loess；DEM model；water content；triaxial shear tests

TU432

A

10.13244/j.cnki.jiwhr.2016.06.000

1672-3031（2016）06-0448-06

（责任编辑：李琳）

2015-08-26

国家自然科学基金项目（51178392）；陕西省教育厅科研计划项目（15JK1117）

郭鸿（1984-），男，汉，陕西长武人，博士，讲师，主要从事岩土离散元数值模拟研究。E-mail：aaronkwo@qq.com