基于離散元的顆粒材料三維臨界狀態(tài)與剪脹特性研究

劉嘉英 ，馬 剛 ，周 偉，常曉林

（1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072；2.水工巖石力學教育部重點實驗室，湖北 武漢 430072）

1 研究背景

臨界狀態(tài)是砂土等顆粒材料變形過程中達到的極限狀態(tài)。經(jīng)過幾十年的發(fā)展，Roscoe等［1］提出的臨界狀態(tài)理論［2］已成為砂土本構(gòu)框架的重要組成部分。剪脹性也是砂土、粗粒土等顆粒材料的一個重要特性。經(jīng)典的應(yīng)力剪脹理論［3］構(gòu)建了土體的變形特性與應(yīng)力狀態(tài)之間的關(guān)系，被廣泛應(yīng)用于砂土彈塑性本構(gòu)模型中，但忽略了其內(nèi)部結(jié)構(gòu)特性對剪脹的影響［4］。因此，不少學者［5-9］引入不同形式的狀態(tài)參量，提出了狀態(tài)相關(guān)的剪脹模型，并建立了相關(guān)的本構(gòu)關(guān)系。如Li等［7-9］引入了與臨界狀態(tài)孔隙比相關(guān)的狀態(tài)參量，所提出的剪脹模型在Toyoura砂的室內(nèi)試驗中得到了驗證。

目前常采用常規(guī)三軸試驗研究砂土等顆粒材料的力學特性，大多數(shù)的本構(gòu)模型、經(jīng)典臨界狀態(tài)理論等也建于此基礎(chǔ)上。而在實際工程中，由于各種因素的影響，砂土等材料常處于三向不等的應(yīng)力狀態(tài)，屬于典型的三維問題。真三軸試驗是研究土體三向不等應(yīng)力狀態(tài)的重要手段，Lade等［10-15］分別對Monterey砂、Santa Monica Beach砂、Nevada砂進行了真三軸試驗，研究了強度、變形、剪切帶以及各向異性等問題。扈萍等［16］、張敏等［17］、施維成等［18-19］研究了中主應(yīng)力系數(shù)對粉細砂、福建標準砂及粗粒土應(yīng)力-應(yīng)變關(guān)系及強度指標的影響。然而受試驗條件的限制，真三軸物理試驗?zāi)軌蚣虞d到臨界狀態(tài)且不產(chǎn)生局部化變形的很少［20］。此外，離散單元法（DEM）［21-26］、連續(xù)-離散耦合方法（FDEM）［27-29］等非連續(xù)介質(zhì)力學方法，使得從宏細觀兩個層面研究了砂土、粗粒土等顆粒材料在三維應(yīng)力條件下的變形和強度特性成為可能。

在三維應(yīng)力條件下，顆粒材料的臨界狀態(tài)和剪脹特性有別于軸對稱的應(yīng)力條件，然而受限于試驗條件等因素，所涉及的研究不多且存在一些不同看法。關(guān)于顆粒材料在三維應(yīng)力條件下臨界狀態(tài)線的唯一性，學術(shù)界沒有一致的結(jié)論，如Zhao等［30］根據(jù)離散元數(shù)值試驗結(jié)果，認為顆粒材料在三維應(yīng)力條件下臨界狀態(tài)線有唯一性，但Huang等［31］同樣采用離散元方法，卻得出不同中主應(yīng)力系數(shù)對顆粒材料臨界狀態(tài)線有所影響的結(jié)論。而在理論上［32］、室內(nèi)試驗［18］中，關(guān)于臨界狀態(tài)線的唯一性亦有不同的論述。此外，砂土等顆粒材料在三維應(yīng)力條件下的剪脹特性亦較為復雜，不同的學者采用不同的方法，對顆粒材料的剪脹特性進行了多種數(shù)學描述。施維成等［18-19］、Xiao等［33-34］根據(jù)粗粒土真三軸試驗結(jié)果，研究了三維空間的應(yīng)力剪脹關(guān)系，并提出了多個剪脹模型。馬剛等［28-29］采用連續(xù)-離散耦合的數(shù)值分析方法，研究了堆石體在真三軸應(yīng)力路徑下的非共軸性和剪脹特性。總體而言，盡管在顆粒材料三維應(yīng)力條件下的臨界狀態(tài)和剪脹特性方面有了一些研究成果，但由于缺乏足夠的試驗數(shù)據(jù)，仍需要作進一步的研究和探討。

砂土等顆粒材料是由大量離散的固體顆粒相互作用組成的復雜體系，采用細觀數(shù)值模擬方法（如離散單元法DEM）可以再現(xiàn)顆粒材料的宏觀力學特性，并揭示其復雜力學特性背后的細觀機理。在細觀數(shù)值模擬中，在對顆粒的形狀、級配和接觸特性等做了一些簡化，但其能準確控制試樣的級配、排水等條件，精確實現(xiàn)預定的加載條件和加載路徑，避免試樣在加載過程中受外界因素的影響。采用周期性邊界控制試樣的加載，能夠最大程度地避免剪切帶的產(chǎn)生，使試樣整體能夠在剪切過程中達到臨界狀態(tài)。

本文將基于顆粒離散元方法，采用開源軟件LIGGGHTS（LAMMPS Improved For General Granular and Granular Heat Transfer Simulations）［35］進行多組顆粒材料的真三軸數(shù)值試驗，研究了中主應(yīng)力對顆粒材料剪脹特性和臨界狀態(tài)的影響。根據(jù)離散元數(shù)值試驗數(shù)據(jù)，分析三維應(yīng)力條件下顆粒材料臨界狀態(tài)的唯一性，對Li等［8］提出的狀態(tài)相關(guān)剪脹方程進行驗證，并通過角隅函數(shù)修正三維情況下的模型參數(shù)，為建立合理的顆粒材料本構(gòu)模型提供一定的參考。

2 真三軸數(shù)值試驗

2.1 周期性邊界應(yīng)力控制 進行DEM模擬時，不同的邊界條件會對顆粒材料的力學響應(yīng)有一定的影響。剛性邊界、柔性邊界以及周期性邊界為較為常用的邊界條件。本文研究需要將試樣加載到臨界狀態(tài)，并盡量避免應(yīng)變局部化現(xiàn)象的發(fā)生。周期性邊界條件下，邊界溢出顆粒將循環(huán)出現(xiàn)在對稱的邊界上。采用周期性邊界進行數(shù)值模擬，能夠反映局部推廣到整體的性質(zhì)，減小試樣在加載過程中受邊界的影響，并且使試樣在空間中保持較為均勻的應(yīng)變場［21，31］。因此本文采用周期性邊界進行真三軸數(shù)值試驗。

采用周期性邊界進行應(yīng)力控制時，顆粒集合體的應(yīng)力張量可以表示為：式中：V為體積單元的體積；Np為體積單元內(nèi)顆粒的數(shù)量；為顆粒p的應(yīng)力張量；Vp為顆粒p的體積。

在對周期性邊界進行應(yīng)力控制時，Thornton等［36］、Huang等［31］采用控制加載過程中應(yīng)變率的方法，使顆粒集合體的應(yīng)力狀態(tài)（由式（1）計算所得）達到預定的值以滿足真三軸應(yīng)力路徑。本文應(yīng)用Huang等［31］的算法進行周期性邊界應(yīng)力控制，加載過程中各個方向的應(yīng)變率采用如下方法控制：

式中：ε˙為試樣加載的應(yīng)變率；σ*為試樣在某方向需要達到的應(yīng)力值；σ為當前的應(yīng)力值；g為增益參數(shù)；ε˙max為加載過程中允許的最大應(yīng)變率。通過試算得到合適的ε˙max，使應(yīng)變率ε˙在加載過程中合理變化，從而使顆粒集合體的應(yīng)力狀態(tài)符合預定的應(yīng)力路徑。

2.2 數(shù)值試樣與參數(shù) 制樣時，首先在立方體空間隨機生成19 436個無接觸的圓球顆粒，粒徑分布如圖1所示，平均粒徑d50為6.25 mm。為了避免由制樣產(chǎn)生的初始各向異性，首先在試樣的各個方向采用位移控制等速地壓縮試樣直至目標大小，然后給試樣施加三向等壓應(yīng)力直至達到預定的圍壓值。本文數(shù)值試驗考慮了6個級別的圍壓（0.5、1.0、2.0、4.0、8.0和16.0 MPa），其對應(yīng)的初始孔隙比、平均配位數(shù)如表1所示。在數(shù)值試驗中不考慮重力，采用位移等速壓縮進行制樣，與室內(nèi)試驗有所區(qū)別，但這種制樣方法能減小制樣過程產(chǎn)生的各向異性，生成初始各向同性數(shù)值試樣，從而避免加載過程中試樣的原生各向異性的影響。

顆粒間接觸模型采用Hertz-Mindlin模型，并引入了轉(zhuǎn)動阻抗［37］，即在顆粒動力平衡系統(tǒng)中添加一個轉(zhuǎn)動阻力矩，采用轉(zhuǎn)動摩擦系數(shù)控制材料的抗轉(zhuǎn)動特性。本文數(shù)值試驗采用的細觀參數(shù)如表2所示。

圖1 數(shù)值試樣與粒徑分布

表1 不同圍壓試樣表

表2 數(shù)值模擬細觀參數(shù)

2.3 真三軸加載路徑 采用中主應(yīng)力系數(shù)b或應(yīng)力羅德角θσ反映3個主應(yīng)力之間的關(guān)系：

式中：中主應(yīng)力系數(shù)b的取值范圍為0≤b≤1。當b=0時，σ2=σ3，應(yīng)力羅德角θσ=0，試樣處于三軸壓縮狀態(tài)；當b=1時，σ2=σ1，應(yīng)力羅德角試樣處于三軸拉伸狀態(tài)。本文進行的真三軸數(shù)值試驗的b值分別取為0.0，0.25，0.5，0.75和1.0。

采用應(yīng)力不變量廣義剪應(yīng)力q、平均靜水壓力p描述顆粒材料的應(yīng)力狀態(tài)：

體積應(yīng)變εv、剪應(yīng)變εq用3個主應(yīng)變表示為：

本文通過控制周期性邊界的移動速率來控制試驗各個方向應(yīng)力的大小。進行等p等b數(shù)值試驗時，施加在試樣上的中主應(yīng)力σ2和小主應(yīng)力σ3分別為：

進行等σ3等b數(shù)值試驗時，小主應(yīng)力σ3不變，施加在試樣上的中主應(yīng)力σ2為：

2.4 宏觀力學響應(yīng) 初始圍壓為0.5 MPa的等p等b數(shù)值試驗結(jié)果如圖2、圖3所示，各試樣在加載過程中保持平均應(yīng)力p=0.5 MPa。

圖2 不同加載路徑的廣義剪應(yīng)力演化曲線

圖3 不同加載路徑的體積應(yīng)變演化曲線

圖2為廣義剪應(yīng)力q隨大主應(yīng)變ε1的演化曲線。中主應(yīng)力系數(shù)的影響體現(xiàn)在，隨著b的增加，峰值廣義剪應(yīng)力減小，其對應(yīng)的大主應(yīng)變和臨界值也逐漸減小。圖3為體積應(yīng)變εv隨大主應(yīng)變ε1的演化曲線，體積應(yīng)變以剪縮為正，剪脹為負。由于試樣比較密實，且不可破碎，在加載初期未出現(xiàn)明顯的剪縮。加載過程中不同中主應(yīng)力系數(shù)對應(yīng)的試樣均產(chǎn)生了明顯的剪脹，且b值越大，剪脹越明顯。綜合圖2和圖3，當加載到軸向應(yīng)變?yōu)?0%～40%時，試樣的應(yīng)力狀態(tài)、體積應(yīng)變均不再變化，可以認為試樣均達到了臨界狀態(tài)。

3 復雜應(yīng)力路徑下的臨界狀態(tài)

3.1 臨界狀態(tài)線 Li等［38］對有效靜水壓力進行歸一化，提出了砂土等無黏性顆粒材料的臨界狀態(tài)線表達式：

式中：ec、p′c分別為臨界狀態(tài)下的孔隙比與有效靜水壓力；pa為標準大氣壓；eΓ、λ、ξ為模型參數(shù)。根據(jù)式（12）擬合不同初始圍壓下等p等b、等σ3等b應(yīng)力路徑試驗得到的臨界狀態(tài)點，得到的臨界狀態(tài)線如圖4所示。擬合參數(shù)值eΓ=0.699，λ=0.0002，ξ=0.971。由圖4可見，不同中主應(yīng)力系數(shù)對應(yīng)的臨界點有所離散，但離散程度較低，基本落在擬合臨界線附近。而在對式（12）進行數(shù)據(jù)擬合時，其相關(guān)系數(shù)達到0.96，因此，本文認為，在不同的中主應(yīng)力系數(shù)的真三軸試驗，其臨界狀態(tài)線CSL是唯一的。

對于真三軸加載路徑下臨界狀態(tài)線的唯一性問題，目前的研究存在一定分歧。Zhao等［30］認為臨界狀態(tài)線是唯一的，與中主應(yīng)力系數(shù)無關(guān)，而Huang等［31］則認為臨界狀態(tài)線與中主應(yīng)力系數(shù)相關(guān)。將Zhao等［30］、Huang等［31］的臨界狀態(tài)數(shù)據(jù)點繪于圖5中。二者數(shù)值模擬所得的臨界狀態(tài)數(shù)據(jù)點在e-lgp平面內(nèi)的臨界狀態(tài)線附近也僅有較小的離散，二者結(jié)論的差異可能是因為其臨界狀態(tài)線的坐標取值范圍不同而造成的。Huang等［31］對于不同中主應(yīng)力系數(shù)的臨界狀態(tài)數(shù)據(jù)較大的離散性，可能是因為其計算所得的臨界狀態(tài)孔隙比變化范圍較小，從而使同一圍壓下對應(yīng)的孔隙比的隨機性差異放大所致。而在圖5中，在孔隙比變化較大的坐標范圍下，Huang等［31］的不同的臨界狀態(tài)線也基本重合。Li等［32］也通過理論證明了臨界狀態(tài)線（CSL）的唯一性，因此可認為本文對臨界狀態(tài)線的唯一性的認識是合理的。本文的后續(xù)工作也是基于臨界狀態(tài)線唯一性的前提。

圖4 臨界狀態(tài)線

圖5 其他文獻中的臨界狀態(tài)數(shù)據(jù)

3.2 臨界應(yīng)力比 為了反映中主應(yīng)力對應(yīng)力比的的影響，一般將臨界應(yīng)力比M表示為應(yīng)力羅德角θσ的函數(shù)

式中Mc為三軸壓縮條件(b=0)對應(yīng)的臨界狀態(tài)應(yīng)力比。

角隅函數(shù)g(θσ)有多種形式，本文數(shù)值試驗中試樣的摩擦角較?。é眨?2°），因此，采用形式較為簡單的角隅函數(shù)即可滿足不同應(yīng)力羅德角對應(yīng)的應(yīng)力比關(guān)系［39］

式中K為三軸拉伸與壓縮條件下的臨界狀態(tài)應(yīng)力比的比值。

將數(shù)值試驗得到的臨界應(yīng)力比與式（13）預測的臨界應(yīng)力比繪于π平面內(nèi)，如圖6所示?？梢钥闯觯鶕?jù)公式計算所得的臨界應(yīng)力比與數(shù)值試驗結(jié)果基本一致。關(guān)于顆粒材料在三維條件下臨界狀態(tài)線以及臨界狀態(tài)應(yīng)力比的規(guī)律，與以往的試驗研究結(jié)果相似，也進一步驗證了本文離散元數(shù)值模擬結(jié)果的合理性。

4 真三軸試驗的剪脹特性

4.1 狀態(tài)相關(guān)剪脹模型 在三維應(yīng)力空間，顆粒材料的剪脹因子d定義為塑性體變增量與塑性剪應(yīng)變增量之比：

砂土等無黏性顆粒材料的剪脹特性不僅與應(yīng)力比有關(guān)，還與材料的內(nèi)部結(jié)構(gòu)有關(guān)。Been等［40］將砂土的臨界狀態(tài)作為參考狀態(tài)，定義了狀態(tài)參量ψ，其值為當前有效靜水壓力下的孔隙比e與臨界孔隙比ec之差：

圖6 π平面內(nèi)臨界狀態(tài)應(yīng)力比

圖7 狀態(tài)參數(shù)量示意

圖7為狀態(tài)參量在e-p平面內(nèi)的示意，A點在臨界狀態(tài)線下方，ψ＜0，表示砂土的狀態(tài)較為緊密，剪切時會發(fā)生剪脹；B點在臨界狀態(tài)線上方，ψ＞0，表示砂土的狀態(tài)較松，剪切時會發(fā)生剪縮。

在經(jīng)典劍橋模型中，應(yīng)力剪脹關(guān)系表示為d=M-η，其中M為臨界狀態(tài)應(yīng)力比，η為當前應(yīng)力比。此模型未考慮砂土的內(nèi)部狀態(tài)。Li等［8］提出了基于狀態(tài)相關(guān)的砂土剪脹模型，將剪脹因子d與應(yīng)力比η的關(guān)系表示為：

式中d0和m為模型參數(shù)，可從常規(guī)三軸試驗結(jié)果得到。當M=0d0=M時式（17）還原為劍橋模型的剪脹方程。

采用圍壓0.5 MPa，b=0條件下的等p等b數(shù)值試驗數(shù)據(jù)，根據(jù)式（17）進行擬合，得到m=4.570，d0=0.619。將式（17）預測的d-η關(guān)系曲線繪于圖8。由于本文數(shù)值試驗采用的試樣較為密實，因此d-η曲線存在一個明顯的轉(zhuǎn)折點，即剪脹因子在峰值應(yīng)力比時達到最小值，此時試樣的剪脹速率最快。在該點之前，d-η曲線大致呈直線型式；達到該點以后，應(yīng)力比逐漸減小，剪脹因子逐漸增加直至為0，此時應(yīng)力比和剪脹因子均不再變化，達到臨界狀態(tài)。除了少數(shù)數(shù)據(jù)點偏離預測值外，式（17）較為準確地預測了數(shù)值試驗中剪脹因子d與應(yīng)力比η的關(guān)系曲線，因此可應(yīng)用狀態(tài)相關(guān)剪脹理論描述顆粒材料的剪脹特性。

4.2 剪脹方程的三維化 砂土的實際應(yīng)力條件并非軸對稱，是典型的三維問題。Li［9］在原始模型的基礎(chǔ)上，將剪脹模型三維化，即將式（18）中的臨界應(yīng)力比表示為Lode角的函數(shù)。黃茂松等［41］在構(gòu)建松砂靜態(tài)液化本構(gòu)模型時也采用了類似的型式：

圖8 三軸壓縮條件下的剪脹-應(yīng)力比關(guān)系

圖9 b=1.0條件下式（18）預測的應(yīng)力-剪脹曲線與數(shù)值實驗結(jié)果對比

式（18）所描述的剪脹模型，當狀態(tài)參量ψ變化很小時（如加載初期），d-η曲線近似為直線。該直線的斜率為與應(yīng)力羅德角無關(guān)，即不同中主應(yīng)力系數(shù)對應(yīng)的d-η曲線初始段的斜率是相同的。而本文的離散元數(shù)值模擬結(jié)果表明，應(yīng)力-剪脹曲線的初始段斜率是與中主應(yīng)力系數(shù)有關(guān)的，如圖8、圖9所示。已有試驗研究也表明應(yīng)力-剪脹曲線的初始直線段的斜率隨著應(yīng)力羅德角的變化而變化［21］，因此需要對三維剪脹模型進一步修正。由于顆粒材料在三維應(yīng)力條件下的應(yīng)力-剪脹關(guān)系跟應(yīng)力羅德角有關(guān)，因此本文擬根據(jù)式（17），結(jié)合數(shù)據(jù)擬合的結(jié)果，將其中的模型參數(shù)d0和m表示為與角隅函數(shù)g(θσ)有關(guān)的形式，即：

式中d0與m為b=0.0即三軸壓縮條件下的模型參數(shù)值。

采用式（19）預測中主應(yīng)力系數(shù)為0.25，0.5，0.75和1.0的d-η關(guān)系曲線，將其與數(shù)值試驗結(jié)果繪于圖10。由圖10（a）—（d）可知，根據(jù)式（19）預測的應(yīng)力剪脹關(guān)系與數(shù)值試驗結(jié)果基本吻合。圖10（e）為不同b值對應(yīng)的剪脹曲線，可以發(fā)現(xiàn)，在初始的近似直線段，隨著中主應(yīng)力系數(shù)的增加，其斜率和截距均不斷增加。圖10（f）對比了在中主應(yīng)力系數(shù)為1.0時，本文模型與黃茂松等［41］采用模型的預測結(jié)果。由此可以明顯看出，本文提出的剪脹方程在較大的中主應(yīng)力系數(shù)條件下對應(yīng)力-剪脹曲線的拐點及斜率預測能力較好。為量化式（18）與式（19）的擬合結(jié)果，表3給出了二者的預測值與離散元數(shù)值模擬值的平均均方差，可以看出，隨著中主應(yīng)力系數(shù)b的增加，文獻［41］給出的模型的預測值與數(shù)值模擬結(jié)果的相離程度較大，而本文給出的公式計算所得的均方差則基本都保持在一個相對小的范圍內(nèi)。

表3 公式預測與離散元數(shù)值模擬的平均均方差

圖10 本文模型預測的三維剪脹-應(yīng)力比關(guān)系

4.3 三維剪脹模型的驗證 采用不同應(yīng)力加載路徑的數(shù)值試驗數(shù)據(jù)驗證本文提出的三維剪脹模型。圖11為初始圍壓為2 MPa的等p等b應(yīng)力路徑對應(yīng)的剪脹曲線，圖12為圍壓為0.5 MPa的等σ3等b應(yīng)力路徑對應(yīng)的應(yīng)力剪脹曲線。由圖11、圖12可知，式（19）預測的d-η關(guān)系曲線與離散元數(shù)值試驗結(jié)果非常接近，且不同中主應(yīng)力系數(shù)對應(yīng)的應(yīng)力-剪脹曲線規(guī)律與圖10類似。因此，本文提出的剪脹方程在一定程度上可以描述顆粒材料在三維應(yīng)力條件下的剪脹特性。

圖12 等σ3等b實驗對應(yīng)的剪脹-應(yīng)力比關(guān)系

圖13 Xiao等［21］的試驗數(shù)據(jù)

5 剪脹模型的討論

本文提出的基于臨界狀態(tài)相關(guān)的三維剪脹模型，是基于離散元數(shù)值試驗建立的。該剪脹模型在初始加載直線段的規(guī)律較好。Xiao等［21］進行了堆石體真三軸室內(nèi)試驗，并擬合了剪脹與應(yīng)力比之間的直線關(guān)系。將Xiao等［21］的實驗數(shù)據(jù)繪于與本文一致的坐標中，如圖13所示。可以發(fā)現(xiàn)直線的斜率以及截距隨中主應(yīng)力系數(shù)的變化趨勢與本文提出的剪脹模型基本一致。由于Xiao等［21］的室內(nèi)試驗沒有出現(xiàn)應(yīng)變軟化段，故其可以采用直線較好地描述堆石體的剪脹性。

本文提出的剪脹方程，是在Li等［8］提出的狀態(tài)相關(guān)剪脹模型的基礎(chǔ)上，將模型參數(shù)表示為角隅函數(shù)的型式。該公式可以較好地描述離散元模擬條件下的應(yīng)力-剪脹關(guān)系，可以反映中主應(yīng)力系數(shù)對顆粒材料剪脹特性影響的一般規(guī)律，但是其適用性以及適用范圍等仍然需要更多的理論分析以及物理試驗驗證。

現(xiàn)如今關(guān)于砂土、粗粒土等材料的真三軸試驗很多，然而能真正達到臨界狀態(tài)的真三軸試驗數(shù)據(jù)并不多，因此本文沒有采用更多的試驗數(shù)據(jù)進行驗證。筆者擬開展真三軸物理試驗，采用松砂、密砂兩種特征的土對本文提出的公式進行驗證。

采用離散元細觀數(shù)值模擬探討顆粒材料的宏觀臨界狀態(tài)與剪脹特性，是本文的一個嘗試，希望能夠從數(shù)值模擬中得到一定規(guī)律以對顆粒材料的性質(zhì)有更深入的認知。顆粒材料的剪脹性從本質(zhì)上是顆粒間、顆粒與孔隙等結(jié)構(gòu)之間的相互影響和共同作用的結(jié)果。本文進行的數(shù)值試驗及提出的剪脹方程，亦可從微觀的角度進行探討和解釋，如組構(gòu)量、各向異性等，進一步探討剪脹的產(chǎn)生機理，從而加深對土體變形本質(zhì)的認識。而在離散元數(shù)值模擬中，對于顆粒形狀、顆粒尖角等因素的考量也是必要的，本文采用的轉(zhuǎn)動阻抗還不足以模擬這種效應(yīng)，因此對復雜顆粒形狀的離散元模擬以及其對應(yīng)力-剪脹關(guān)系的影響也將是筆者下一步要展開的工作。

6 結(jié)論

本文采用顆粒離散元方法，進行了顆粒材料的真三軸數(shù)值試驗，分析了顆粒材料在真三軸應(yīng)力狀態(tài)下的臨界狀態(tài)和剪脹特性，主要結(jié)論如下：

（1）在等p等b真三軸應(yīng)力路徑加載試驗中，應(yīng)力應(yīng)變規(guī)律符合已有的物理試驗和數(shù)值模擬結(jié)果。隨著中主應(yīng)力系數(shù)的增大，廣義剪應(yīng)力的峰值和臨界值均減小，體積應(yīng)變在初期表現(xiàn)出明顯的剪脹，到加載結(jié)束時基本保持不變。

（2）試樣加載結(jié)束時各試樣均達到了臨界狀態(tài)。擬合e-lgp平面臨界狀態(tài)線，不同中主應(yīng)力系數(shù)對應(yīng)的e-lgp數(shù)據(jù)點落在唯一的臨界狀態(tài)線附近；q-p平面的臨界狀態(tài)應(yīng)力比可以采用角隅函數(shù)近似描述。

（3）狀態(tài)相關(guān)剪脹理論可以較好地預測密實顆粒集合體的剪脹特性。通過引入角隅函數(shù)修正模型參數(shù)，提出了一個三維條件下狀態(tài)相關(guān)的剪脹模型，可以近似反映真三軸應(yīng)力狀態(tài)下的應(yīng)力剪脹關(guān)系。離散元數(shù)值結(jié)果可以從一定程度上反映顆粒材料的應(yīng)力-剪脹規(guī)律，但該公式仍需更多的理論研究和物理試驗驗證。

參 考 文 獻：

［1］ ROSCOE K H，SCHOFIELD A N，WROTH C P.On the yielding of soils［J］.Geotechnique，1958，8（1）：22-53.

［2］ SCHOFIELD A，WROTH P.Critical State Soil Mechanics［M］.US：Mc Graw Hill.1968.

［3］ ROWE P W.The stress-dilatancy relations for static equilibrium of an assembly of particles in contact［J］.Proc R Soc London Ser A，1962，269（1339）：500-527.

［4］ 蔡正銀，李相菘.砂土的剪脹理論及其本構(gòu)模型的發(fā)展［J］.巖土工程學報，2007，29（8）：1122-1128.

［5］ MANZARI M T，DAFALIAS Y F.A critical state two-surface plasticity model for sands［J］.Geotechnique，1997，47（2）：255-272.

［6］ WAN R G，GUO P J.A simple constitutive model for granular soils：modified stress-dilatancy approach［J］.Computers and Geotechnics，1998，22（2）：109-133.

［7］ LI X S，DAFALIAS Y F，WANG Z L.State-dependent dilatancy in critical-state constitutive modelling of sand［J］.Canadian Geotechnical Journal，1999，36（4）：599-611.

［8］ LI X S，DAFALIAS Y F.Dilatancy for cohesionless soils［J］.Geotechnique，2000，50（4）：449-460.

［9］ LI X S.A sand model with state-dependent dilatancy［J］.Geotechnique，2002，52（3）：173-186.

［10］ LADE P V，DUNCAN J M.Cubical triaxial tests on cohesionless soil［J］.Journal of Geotechnical and Geoenvi?ronmental Engineering，1973，101（10）：793-812.

［11］ WANG Q，LADE P V.Shear banding in true triaxial tests and its effect on failure in sand［J］.Journal of Engi?neering Mechanics，2001，127（8）：754-761.

［12］ LADE P V，WANG Q.Analysis of shear banding in true triaxial tests on sand［J］.Journal of Engineering Mechan?ics，2001，127（8）：762-768.

［13］ ABELEV A V，LADE P V.Effects of cross anisotropy on three-dimensional behavior of sand.I：Stress-strain be?havior and shear banding［J］.Journal of Engineering Mechanics，2003，129（2）：160-166.

［14］ LADE P V，ABELEV A V.Effects of cross anisotropy on three-dimensional behavior of sand.II：Volume change behavior and failure［J］.Journal of Engineering Mechanics，2003，129（2）：167-174.

［15］ RODRIGUEZ N M，LADE P V.True triaxial tests on cross-anisotropic deposits of fine Nevada sand［J］.Interna?tional Journal of Geomechanics，2013，13（6）：779-793.

［16］ 扈萍，黃茂松，馬少坤，等.粉細砂的真三軸試驗與強度特性［J］.巖土力學，2011，32（2）：465-470.

［17］ 張敏，許成順，杜修力，等.中主應(yīng)力系數(shù)及應(yīng)力路徑對砂土剪切特性影響的真三軸試驗研究［J］.水利學報，2015，46（9）：1072-1079.

［18］ 施維成，朱俊高，代國忠，等.粗粒土在π平面上的真三軸試驗及強度準則［J］.河海大學學報：自然科學版，2015，43（1）：11-15.

［19］ 施維成，朱俊高，代國忠，等.球應(yīng)力和偏應(yīng)力對粗粒土變形影響的真三軸試驗研究［J］.巖土工程學報，2015，37（5）：776-783.

［20］ XIAO Y，SUN Y，LIU H，et al.Critical state behaviors of a coarse granular soil under generalized stress condi?tions［J］.Granular Matter，2016，18（2）：1-13.

［21］ THORNTON C.Numerical simulations of deviatoric shear deformation of granular media［J］.Géotechnique，2000，50（1）：43-53.

［22］ NG T T.Macro-and micro-behaviors of granular materials under different sample preparation methods and stress paths［J］.International Journal of Solids and Structures，2004，41（21）：5871-5884.

［23］ MAHMUD SAZZAD M，SUZUKI K，MODARESSI-FARAHMAND-RAZAVI A.Macro-micro responses of gran?ular materials under different b values using DEM［J］.International Journal of Geomechanics，2012，12（3）：220-228.

［24］ ZHOU W，LIU J，MA G，et al.Three-dimensional DEM investigation of critical state and dilatancy behaviors of granular materials［J］.Acta Geotechnica，2017，12（3）：527-540.

［25］ ZHOU W，YANG L，MA G，et al.Macro-micro responses of crushable granular materials in simulated true triaxi?al tests［J］.Granular Matter，2015，17（4）：497-509.

［26］ 周偉，謝婷蜓，馬剛，等.基于顆粒流程序的真三軸應(yīng)力狀態(tài)下堆石體的變形和強度特性研究［J］.巖土力學，2012，33（10）：3006-3012.

［27］ 周偉，劉東，馬剛，等.基于隨機散粒體模型的堆石體真三軸數(shù)值試驗研究［J］.巖土工程學報，2012，34（4）：748-755.

［28］ 馬剛，劉嘉英，常曉林，等.堆石體在真三軸應(yīng)力狀態(tài)下的非共軸性與剪脹特性［J］.中南大學學報：自然科學版，2016，47（5）：1697-1707.

［29］ MA G，CHANG X L，ZHOU W，et al.Mechanical response of rockfills in a simulated true triaxial test：A com?bined FDEM study［J］.Geomechanics and Engineering，2014，7（3）：317-333.

［30］ ZHAO J，GUO N.Unique critical state characteristics in granular media considering fabric anisotropy［J］.Géotechnique，2013，63（8）：695-704.

［31］ HUANG X，HANLEY K J，O’SULLIVAN C，et al.DEM analysis of the influence of the intermediate stress ratio on the critical-state behaviour of granular materials［J］.Granular Matter，2014，16（5）：641-655.

［32］ LI X S，DAFALIAS Y F.Anisotropic critical state theory：role of fabric［J］.Journal of Engineering Mechanics，2011，138（3）：263-275.

［33］ XIAO Y，LIU H L，ZHU J G，et al.Dilatancy equation of rockfill material under the true triaxial stress condition［J］.Science China Technological Sciences，2011，54（1）：175-184.

［34］ XIAO Y，LIU H，SUN Y，et al.Stress-dilatancy behaviors of coarse granular soils in three-dimensional stress space［J］.Engineering Geology，2015，195：104-110.

［35］ KLOSS C，GONIVA C.Liggghts-a new open source discrete element simulation software［C］//Proc.of The 5th In?ternational Conference on Discrete Element Methods.2010：25-26.

［36］ THORNTON C，ZHANG L.On the evolution of stress and microstructure during general 3D deviatoric straining of granular media［J］.Geotechnique，2010，60（5）：333-341.

［37］ AI J，CHEN J F，ROTTER J M，et al.Assessment of rolling resistance models in discrete element simulations［J］.Powder Technology，2011，206（3）：269-282.

［38］ LI X S，WANG Y.Linear representation of steady-state line for sand［J］.Journal of Geotechnical and Geoenvi?ronmental Engineering，1998，124（12）：1215-1217.

［39］ 鄭穎人，沈珠江，龔曉南.巖土塑性力學原理［M］.北京：中國建筑工業(yè)出版社，2002.

［40］ BEEN K，JEFFERIES M G.A state parameter for sands［J］.Géotechnique，1985，35（2）：99-112.

［41］ 黃茂松，曲勰，呂璽琳.基于狀態(tài)相關(guān)本構(gòu)模型的松砂靜態(tài)液化失穩(wěn)數(shù)值分析［J］.巖石力學與工程學報，2014，33（7）：1479-1487.