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

（1.武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072；2.水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

1 研究背景

臨界狀態(tài)是砂土等顆粒材料變形過(guò)程中達(dá)到的極限狀態(tài)。經(jīng)過(guò)幾十年的發(fā)展，Roscoe等［1］提出的臨界狀態(tài)理論［2］已成為砂土本構(gòu)框架的重要組成部分。剪脹性也是砂土、粗粒土等顆粒材料的一個(gè)重要特性。經(jīng)典的應(yīng)力剪脹理論［3］構(gòu)建了土體的變形特性與應(yīng)力狀態(tài)之間的關(guān)系，被廣泛應(yīng)用于砂土彈塑性本構(gòu)模型中，但忽略了其內(nèi)部結(jié)構(gòu)特性對(duì)剪脹的影響［4］。因此，不少學(xué)者［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)試驗(yàn)中得到了驗(yàn)證。

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

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

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

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

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

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

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

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

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

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

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

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

表1 不同圍壓試樣表

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

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

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

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

體積應(yīng)變?chǔ)舦、剪應(yīng)變?chǔ)舚用3個(gè)主應(yīng)變表示為：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 真三軸試驗(yàn)的剪脹特性

4.1 狀態(tài)相關(guān)剪脹模型 在三維應(yīng)力空間，顆粒材料的剪脹因子d定義為塑性體變?cè)隽颗c塑性剪應(yīng)變?cè)隽恐龋?/p>

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

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

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

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

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

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

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

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

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

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

式（18）所描述的剪脹模型，當(dāng)狀態(tài)參量ψ變化很小時(shí)（如加載初期），d-η曲線(xiàn)近似為直線(xiàn)。該直線(xiàn)的斜率為與應(yīng)力羅德角無(wú)關(guān)，即不同中主應(yīng)力系數(shù)對(duì)應(yīng)的d-η曲線(xiàn)初始段的斜率是相同的。而本文的離散元數(shù)值模擬結(jié)果表明，應(yīng)力-剪脹曲線(xiàn)的初始段斜率是與中主應(yīng)力系數(shù)有關(guān)的，如圖8、圖9所示。已有試驗(yàn)研究也表明應(yīng)力-剪脹曲線(xiàn)的初始直線(xiàn)段的斜率隨著應(yīng)力羅德角的變化而變化［21］，因此需要對(duì)三維剪脹模型進(jìn)一步修正。由于顆粒材料在三維應(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ù)測(cè)中主應(yīng)力系數(shù)為0.25，0.5，0.75和1.0的d-η關(guān)系曲線(xiàn)，將其與數(shù)值試驗(yàn)結(jié)果繪于圖10。由圖10（a）—（d）可知，根據(jù)式（19）預(yù)測(cè)的應(yīng)力剪脹關(guān)系與數(shù)值試驗(yàn)結(jié)果基本吻合。圖10（e）為不同b值對(duì)應(yīng)的剪脹曲線(xiàn)，可以發(fā)現(xiàn)，在初始的近似直線(xiàn)段，隨著中主應(yīng)力系數(shù)的增加，其斜率和截距均不斷增加。圖10（f）對(duì)比了在中主應(yīng)力系數(shù)為1.0時(shí)，本文模型與黃茂松等［41］采用模型的預(yù)測(cè)結(jié)果。由此可以明顯看出，本文提出的剪脹方程在較大的中主應(yīng)力系數(shù)條件下對(duì)應(yīng)力-剪脹曲線(xiàn)的拐點(diǎn)及斜率預(yù)測(cè)能力較好。為量化式（18）與式（19）的擬合結(jié)果，表3給出了二者的預(yù)測(cè)值與離散元數(shù)值模擬值的平均均方差，可以看出，隨著中主應(yīng)力系數(shù)b的增加，文獻(xiàn)［41］給出的模型的預(yù)測(cè)值與數(shù)值模擬結(jié)果的相離程度較大，而本文給出的公式計(jì)算所得的均方差則基本都保持在一個(gè)相對(duì)小的范圍內(nèi)。

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

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

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

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

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

5 剪脹模型的討論

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

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

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

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

6 結(jié)論

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

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

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

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

參 考 文 獻(xiàn)：

［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］.巖土工程學(xué)報(bào)，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］ 扈萍，黃茂松，馬少坤，等.粉細(xì)砂的真三軸試驗(yàn)與強(qiáng)度特性［J］.巖土力學(xué)，2011，32（2）：465-470.

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

［18］ 施維成，朱俊高，代國(guó)忠，等.粗粒土在π平面上的真三軸試驗(yàn)及強(qiáng)度準(zhǔn)則［J］.河海大學(xué)學(xué)報(bào)：自然科學(xué)版，2015，43（1）：11-15.

［19］ 施維成，朱俊高，代國(guó)忠，等.球應(yīng)力和偏應(yīng)力對(duì)粗粒土變形影響的真三軸試驗(yàn)研究［J］.巖土工程學(xué)報(bào)，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)下堆石體的變形和強(qiáng)度特性研究［J］.巖土力學(xué)，2012，33（10）：3006-3012.

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

［28］ 馬剛，劉嘉英，常曉林，等.堆石體在真三軸應(yīng)力狀態(tài)下的非共軸性與剪脹特性［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，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］ 鄭穎人，沈珠江，龔曉南.巖土塑性力學(xué)原理［M］.北京：中國(guó)建筑工業(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］.巖石力學(xué)與工程學(xué)報(bào)，2014，33（7）：1479-1487.