郭鴻，黃文華，李軍，魏鴻林

（1.陜西理工大學(xué)土木工程與建筑學(xué)院，陜西漢中723000；

2.西北農(nóng)林科技大學(xué)水利水電科學(xué)研究院，陜西楊凌712100；

3.中國(guó)城市建設(shè)研究院有限公司，北京100120）

非飽和重塑黃土的一種簡(jiǎn)化離散元模型

郭鴻1，2，黃文華1，李軍1，魏鴻林3

（1.陜西理工大學(xué)土木工程與建筑學(xué)院，陜西漢中723000；

2.西北農(nóng)林科技大學(xué)水利水電科學(xué)研究院，陜西楊凌712100；

3.中國(guó)城市建設(shè)研究院有限公司，北京100120）

黃土是一種水敏性顆粒材料，其強(qiáng)度受含水率影響十分明顯?；诜秋柡屯亮W(xué)的基本理論，利用顆粒離散元方法，建立了一種非飽和黃土的簡(jiǎn)化離散元模型。在對(duì)非飽和黃土的力學(xué)特性進(jìn)行離散元模擬分析后，與室內(nèi)三軸試驗(yàn)進(jìn)行了對(duì)比驗(yàn)證。對(duì)比模型和試驗(yàn)結(jié)果發(fā)現(xiàn)：不同含水率試樣的應(yīng)力應(yīng)變關(guān)系受土體顆粒間摩擦系數(shù)直接影響，表現(xiàn)為試樣含水率越高，摩擦系數(shù)越小，說(shuō)明該離散元模型在一定范圍內(nèi)能夠反映重塑黃土的土力學(xué)特性。

非飽和黃土；離散元模型；含水率；三軸試驗(yàn)

1 研究背景

黃土的分布十分廣泛，在中國(guó)分布面積達(dá)64萬(wàn)km2。對(duì)黃土的力學(xué)特性研究非常多，尤其以室內(nèi)試驗(yàn)為主。但是就微觀層面的研究，尤其是基于微觀結(jié)構(gòu)的本構(gòu)模型研究，還需要做很多工作。近年來(lái)，離散元模型越來(lái)越受到廣大科研工作者的青睞，而且其在黃土力學(xué)中作用也越來(lái)越大［1-2］。

離散元法（Discrete Element Method，DEM）［3］是將每一個(gè)顆粒或塊體當(dāng)作一個(gè)單元，觀察記錄數(shù)值試驗(yàn)過(guò)程中單元之間相互作用與運(yùn)動(dòng)規(guī)律，來(lái)預(yù)測(cè)散體群行為的方法。它是一種能準(zhǔn)確、實(shí)時(shí)獲取數(shù)值試驗(yàn)過(guò)程中試樣所有顆粒、應(yīng)力、孔隙等微觀信息的工具［4-6］。蔣明鏡等［7-9］在離散元模擬黃土方面做了大量的研究。Lu Ning等［10-11］基于吸力應(yīng)力理念對(duì)非飽和土做了全新的闡釋。然而，他們的方法由于其復(fù)雜的微觀參數(shù)，在工程實(shí)踐中推廣具有一定的局限性。

離散元所建立的土顆粒之間的膠結(jié)參數(shù)諸如土顆粒之間的黏結(jié)力，無(wú)法和實(shí)際土體相對(duì)應(yīng)。并且，黃土顆粒之間由于不同顆粒之間的膠結(jié)情況隨著土顆粒的形狀以及膠結(jié)鹽的分布等變化非常復(fù)雜。因此，筆者認(rèn)為，黃土剪切破壞的根本原因并不在于顆粒之間膠結(jié)由于“拉伸”引起微觀破壞，進(jìn)而引起宏觀的破壞；而在于膠結(jié)情況發(fā)生變化，顆粒之間的切向接觸發(fā)生了變化，即不同含水率的黃土顆粒，其顆粒在相同法向應(yīng)力的情況下，滑動(dòng)的難易程度不同?；诖死砟?，本文提出了非飽和黃土的一種基于含水率為變量的切向滑動(dòng)離散元模型。該模型運(yùn)用最簡(jiǎn)單的線彈性模型，具有參數(shù)簡(jiǎn)單，預(yù)測(cè)效果理想等特點(diǎn)。

2 顆粒離散元模型

非飽和黃土由于其內(nèi)部存在膠結(jié)物質(zhì)，而且其具有很強(qiáng)的水敏性，水的作用使非飽和黃土的強(qiáng)度大大減小。非飽和土顆粒之間由于水表面張力的存在，顆粒間產(chǎn)生基質(zhì)吸力，其力學(xué)示意圖如圖1所示。顆粒間的水存在兩個(gè)曲面，一個(gè)曲面半徑為r1，另一個(gè)曲面半徑為r2，顆粒半徑為R。

如圖1（a）所示，根據(jù)受力平衡分析［11］，得到基質(zhì)吸力的值為：

其中Ts為水的表面張力。假設(shè)水膜圓弧和顆粒的切點(diǎn)與顆粒之間接觸點(diǎn)形成的圓弧夾角為θ（如圖1（b）所示），則由幾何關(guān)系可以推導(dǎo)出：

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

得：

容易看出，f(θ)是關(guān)于θ的減函數(shù)，這也就意味著Su也是θ的減函數(shù)。

圖1 非飽和土受力分析［11］

假設(shè)沿法向接觸的面積為An，沿切向接觸的面積為As，則顆粒法向和切向的黏結(jié)力分別為：

考慮到以上因素，本文基于Cundall［3］提出的二維圓盤(pán)模型，用顆粒流程序PFC2D（Particle Flow Code in Two Dimensions）對(duì)非飽和黃土顆粒進(jìn)行模擬，其簡(jiǎn)易離散元模型如圖2所示。

圖2 顆粒離散元接觸模型

不考慮水的作用，法向和切向的接觸力為Fn=Knδn和Fs=Ksδs，式中Kn和Ks為分別為接觸顆粒的法向剛度，δn和δs分別為顆粒間法向和切向的相對(duì)位移。最大切向力表示為：

式中μ0為顆粒間的摩擦系數(shù)。

然而由于基質(zhì)吸力的存在，對(duì)于顆粒切向的最大接觸力則為：

令μ=μ0+Ps/δKn=μ(θ)，則式（7）變?yōu)椋?/p>

由于飽和度Sr和接觸角θ的關(guān)系可以表示為［1］：

東西如此，人也如此。一個(gè)人好與不好，不在模樣上?？±实娜瞬灰欢ㄊ呛萌耍舐娜艘膊灰欢ㄊ菈娜?。這么說(shuō)吧，高人不一定個(gè)高，高人是境界高，與身量無(wú)關(guān)。人附加值在哪里呢？在心。一個(gè)人，心好，就是好人，好人好氣息。心壞，就是壞人，壞人壞氣息。所以人也分貴賤，高人則貴，濫人則賤。君子則貴，小人則賤。

故，新的摩擦系數(shù)可以用含水率表示為：

式（8）則變?yōu)椋?/p>

前人的研究也發(fā)現(xiàn)，在一定范圍內(nèi)，內(nèi)摩擦角與試樣含水率由存在著特定關(guān)系［12-14］。因此，筆者建立簡(jiǎn)化模型中主要考慮參數(shù)Kn、Ks和μ即可，其中Kn和Ks與土顆粒的剛度（軟硬）有關(guān)，而μ則與含水率w相對(duì)應(yīng)。

3 非飽和重塑黃土的三軸試驗(yàn)

3.1 試驗(yàn)材料與方案本文采用陜西涇陽(yáng)地區(qū)黃土作為試驗(yàn)原料土，取土深度在3.0～4.0 m之間，測(cè)得天然含水率16.2%，天然干密度1.28 g/cm3，孔隙比1.17，為Q3黃土，試驗(yàn)基本物理性質(zhì)指標(biāo)見(jiàn)表1。

表1 黃土基本物理指標(biāo)

試樣采用非飽和重塑試樣，為原料土風(fēng)干、碾碎、過(guò)篩后加水壓制而成。試樣直徑為39.1 mm，高為80mm，控制試樣密度為天然干密度，含水率為11%、18%和26%。三軸試驗(yàn)采用SBL-1型應(yīng)力應(yīng)變控制式三軸剪切滲透試驗(yàn)儀，控制固結(jié)圍壓為100 kPa、200 kPa和300 kPa作不固結(jié)不排水（UU）三軸剪切試驗(yàn)，試驗(yàn)方案如表2所示。

表2 三軸剪切試驗(yàn)方案

3.2 試驗(yàn)結(jié)果剪切試驗(yàn)采用應(yīng)變控制，試樣裝樣完成后，在軸向加載系統(tǒng)中設(shè)定好軸向應(yīng)變速率，在計(jì)算機(jī)中設(shè)置好應(yīng)力、應(yīng)變記錄時(shí)間間隔，最后進(jìn)行剪切試驗(yàn)。當(dāng)試樣軸向形變大于15%，或試樣出現(xiàn)明顯剪切帶，應(yīng)力值急劇下降，即可停止試驗(yàn)，認(rèn)為試樣達(dá)到破壞。

圖3 三軸試驗(yàn)應(yīng)力應(yīng)變曲線

4 模擬方法及結(jié)果

本次研究基于PC環(huán)境下的顆粒流程序PFC2D，用筆者建立的離散元模型模擬二維黃土的三軸剪切試驗(yàn)。如圖4為實(shí)際三軸試驗(yàn)試樣的顆粒級(jí)配曲線和離散元模型所選定級(jí)配曲線。表3為模擬選定參數(shù)和試驗(yàn)土樣實(shí)際參數(shù)。需要說(shuō)明的是離散元數(shù)值試驗(yàn)的顆粒剛度（法向和切向）是根據(jù)實(shí)際三軸試驗(yàn)標(biāo)定的，而墻體的剛度則取值為顆粒剛度的十分之一，模擬三軸試驗(yàn)的橡皮膜。

表3 離散元模型模擬參數(shù)和試驗(yàn)參數(shù)

圖4 試樣顆粒級(jí)配曲線和模型選定級(jí)配曲線

離散元模擬三軸試驗(yàn)時(shí)，首先根據(jù)表3的參數(shù)生成離散元數(shù)值試樣，其中顆粒的接觸關(guān)系滿足式（11），然后用fish語(yǔ)言編寫(xiě)的圍壓伺服程序通過(guò)周?chē)鷫w（模擬橡皮膜）以及上下墻體（模擬上下加壓板）對(duì)試樣施加固結(jié)圍壓σ3，等循壞穩(wěn)定后，通過(guò)向下移動(dòng)（速度和實(shí)際三軸試驗(yàn)保持一致）試樣上部墻體施加軸向應(yīng)力σ1。最后通過(guò)程序輸出σ1-σ3與ε（豎向位移與試樣初始高度的比值）的關(guān)系。

如圖5所示為PFC2D根據(jù)模型參數(shù)選定值模擬出來(lái)的模擬曲線與實(shí)際三軸試驗(yàn)的試驗(yàn)曲線進(jìn)行的比較分析，圖例中短線前面的數(shù)字代表含水率，短線后面的數(shù)字為圍壓，sim開(kāi)頭的為模擬值，其余為試驗(yàn)值。

從圖5中可以看出，筆者建立離散元模型模擬得到的應(yīng)力應(yīng)變關(guān)系曲線與試驗(yàn)實(shí)際應(yīng)力應(yīng)變關(guān)系曲線能較好吻合，模擬精度在高含水率和高圍壓條件下尤為突出。進(jìn)一步根據(jù)室內(nèi)三軸試驗(yàn)（圖3）和離散元數(shù)值模擬結(jié)果（圖5）繪制莫爾圓，通過(guò)最小二乘法擬合得到內(nèi)摩擦角的試驗(yàn)值和模擬值，如表4所示?？梢钥闯?，雖然試驗(yàn)得到的內(nèi)摩擦角隨含水率增大存在一點(diǎn)起伏，而模型模擬得到的內(nèi)摩擦角隨著含水率在11%至26%間增大而減小，但是兩者總體趨勢(shì)一致，且相同含水率下試驗(yàn)值與模擬值差值較小。

筆者建立簡(jiǎn)化顆粒離散元模型，針對(duì)不同含水率的重塑黃土，只需調(diào)整一個(gè)參數(shù)，即顆粒間的摩擦系數(shù)μ，就可以較為理想地反應(yīng)土體的應(yīng)力應(yīng)變關(guān)系。本文的研究表明該離散元模型在一定范圍內(nèi)能夠反映陜西涇陽(yáng)地區(qū)非飽和重塑黃土的力學(xué)性質(zhì)。從本文非飽和土力學(xué)的觀點(diǎn)看，含水率較小時(shí)，θ減小，基質(zhì)吸力Su增大，顆粒之間黏結(jié)的更緊密，更不容易被剪切錯(cuò)動(dòng)。因此，含水率較低的黃土剪切強(qiáng)度較高，顆粒之間的摩擦系數(shù)和內(nèi)摩擦角較高。

圖5 不同含水率、圍壓下模擬曲線和試驗(yàn)曲線

表4 試驗(yàn)與模擬得到內(nèi)摩擦角比較

5 結(jié)論

（1）本文簡(jiǎn)化離散元模型對(duì)非飽和重塑黃土的三軸剪切試驗(yàn)進(jìn)行數(shù)值模擬，對(duì)比兩者結(jié)果發(fā)現(xiàn)，模型所得應(yīng)力應(yīng)變關(guān)系曲線與實(shí)際試驗(yàn)結(jié)果吻合良好，曲線表現(xiàn)出的趨勢(shì)一致，且在高含水率和高圍壓下模擬效果較為理想。

（2）不同含水率下，試驗(yàn)得出土樣內(nèi)摩擦角與模型模擬得出內(nèi)摩擦角大小相近，且隨著含水率增大，內(nèi)摩擦角總體上呈減小趨勢(shì)。

（3）筆者建立簡(jiǎn)化顆粒離散元模型，針對(duì)不同含水率的重塑黃土，只需調(diào)整一個(gè)參數(shù)，即顆粒間的摩擦系數(shù)，就可以較為理想地反應(yīng)土體的應(yīng)力應(yīng)變關(guān)系。

［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］蔣明鏡，胡海軍.密實(shí)和松散顆粒材料等吸力三軸剪切試驗(yàn)離散元數(shù)值模擬［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，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］蔣明鏡，沈珠江，趙魁芝，等.結(jié)構(gòu)性黃土濕陷性指標(biāo)室內(nèi)測(cè)定方法的探討［J］.水利水運(yùn)科學(xué)研究，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］王浩林，白曉紅，馮俊琴.壓實(shí)黃土狀填土抗剪強(qiáng)度指標(biāo)的影響因素探討［J］.巖土工程學(xué)報(bào)，2010，32（增2）：131-134.

［13］賈亮，朱彥鵬，朱鋆川.蘭州馬蘭、離石壓實(shí)黃土抗剪強(qiáng)度影響因素探討［J］.巖土工程學(xué)報(bào)，2014，36（增2）：120-124.

［14］王鐵行，羅少鋒，劉小軍.考慮含水率影響的非飽和原狀黃土凍融強(qiáng)度試驗(yàn)研究［J］.巖土力學(xué)，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

（責(zé)任編輯：李琳）

2015-08-26

國(guó)家自然科學(xué)基金項(xiàng)目（51178392）；陜西省教育廳科研計(jì)劃項(xiàng)目（15JK1117）

郭鴻（1984-），男，漢，陜西長(zhǎng)武人，博士，講師，主要從事巖土離散元數(shù)值模擬研究。E-mail：aaronkwo@qq.com