堆石料的三維應(yīng)力分?jǐn)?shù)階本構(gòu)模型

吉華，孫逸飛

(1.中國電建集團(tuán) 華東勘測設(shè)計研究院，杭州 310014；2.河海大學(xué) 土木與交通學(xué)院，南京 210098；3.德國波鴻魯爾大學(xué) 土木與環(huán)境工程學(xué)院，波鴻 44801)

堆石料常常被用于水利工程，諸如堆石壩、堤壩等的建造[1]。因此，在此類工程的設(shè)計、建造和運(yùn)營全過程中，堆石料的強(qiáng)度和變形特性顯得尤為重要。由于河流水荷載的復(fù)雜性，堆石料在其全壽命周期內(nèi)經(jīng)常會遭受不同應(yīng)力路徑的外部荷載作用。為了分析不同應(yīng)力路徑下堆石料的強(qiáng)度和變形行為，學(xué)者們[2-6]開展了一系列試驗和理論研究，并建立了相應(yīng)的本構(gòu)模型。Xiao等[5]通過對不同圍壓和初始密度條件下堆石料的強(qiáng)度和剪脹行為進(jìn)行研究，提出了狀態(tài)依賴的堆石料邊界面塑性力學(xué)本構(gòu)模型。陳生水等[7]通過研究顆粒破碎對堆石料臨界狀態(tài)特性等的影響，建立了一個可以考慮堆石料顆粒破碎特性的廣義塑性力學(xué)模型。此外，在實際工程中，堆石料等粗粒土常遭受三維荷載作用，中主應(yīng)力對其強(qiáng)度和變形行為影響較大[8]。為此，Xiao等[9]提出了基于統(tǒng)一破壞準(zhǔn)則的堆石料三維邊界面塑性本構(gòu)模型；姚仰平等[10-11]、Tian等[12]提出了基于轉(zhuǎn)換應(yīng)力法[13]的軟土和粗粒土三維統(tǒng)一硬化模型。為了考慮材料的非關(guān)聯(lián)塑性流動特性，以上這些模型在建立過程中均額外引入了不同于屈服面函數(shù)的塑性勢面函數(shù)。這些塑性勢面函數(shù)往往基于經(jīng)驗假定。那么，是否可以在不引入塑性勢的條件下建立堆石料的三維本構(gòu)模型呢？為此，Sun等[14]通過分析巖土材料的三軸試驗特性，提出了土體的應(yīng)力分?jǐn)?shù)階塑性力學(xué)方法[14-15]。Lu等[16]基于特征應(yīng)力法，提出了軟土的三維分?jǐn)?shù)階塑性本構(gòu)模型。Qu等[17]建立了適用于巖石的分?jǐn)?shù)階塑性力學(xué)模型。李海潮等[18]進(jìn)一步引入溫度效應(yīng)，拓展建立了軟巖的分?jǐn)?shù)階熱塑性模型。這些模型均處于“應(yīng)力分?jǐn)?shù)階塑性力學(xué)”框架，還有一些諸如“時間分?jǐn)?shù)階塑性力學(xué)”框架下的模型[19]，仍然采用屈服面和塑性勢面不一致的假定，這里不再贅述。不同于傳統(tǒng)塑性力學(xué)，應(yīng)力分?jǐn)?shù)階塑性力學(xué)基于非局域性的分?jǐn)?shù)階微分算子，在某一應(yīng)力點求導(dǎo)所得塑性流動方向不僅與該點的應(yīng)力狀態(tài)有關(guān)，還與到達(dá)該點的加載歷史或相對臨界狀態(tài)的距離有關(guān)。在實際建模過程中，僅需對屈服面進(jìn)行一階和分?jǐn)?shù)階求導(dǎo)，即可建立堆石料的非關(guān)聯(lián)彈塑性本構(gòu)模型。然而，目前提出的分?jǐn)?shù)階塑性力學(xué)模型[14,20]僅對粗粒土的三軸應(yīng)力-應(yīng)變特性進(jìn)行了分析，對土體真三軸應(yīng)力-應(yīng)變特性尚未進(jìn)行預(yù)測和模擬。

為解決這一問題并同時完善現(xiàn)有分?jǐn)?shù)階巖土塑性力學(xué)理論體系，有必要建立一個三維應(yīng)力分?jǐn)?shù)階塑性力學(xué)本構(gòu)模型。筆者擬通過特征應(yīng)力法[21]對現(xiàn)有分?jǐn)?shù)階塑性力學(xué)模型進(jìn)行修正，使其可以合理地預(yù)測土體的三維應(yīng)力-應(yīng)變關(guān)系。

1 基本定義

1.1 應(yīng)力-應(yīng)變的定義

(1)

(2)

除此之外，還有平均有效主應(yīng)力(p′)、偏應(yīng)力(q)、中主應(yīng)力系數(shù)(b)這幾個需要用到的定義。

p′=1/3trσ′

(3)

(4)

(5)

式中：偏應(yīng)力張量s=σ′-p′I，I為單位張量。

1.2 分?jǐn)?shù)階導(dǎo)數(shù)

Caputo的分?jǐn)?shù)階左導(dǎo)數(shù)和右導(dǎo)數(shù)分別定義為[23]

σ′>0+

0->σ′

(6b)

1.3 特征應(yīng)力法

(7)

(8)

(9)

2 三維分?jǐn)?shù)階應(yīng)力-應(yīng)變關(guān)系

增量主應(yīng)變可以分解為彈性應(yīng)變和塑性應(yīng)變兩部分，即

(10)

(11)

式中：G是剪切模量；K是體積模量，分別定義為

(12)

(13)

(14)

式中：Kp為硬化模量，其定義將在后面的本構(gòu)模型部分給出。n和m分別為加載方向和塑性流動方向，定義為[24]

(15)

(16)

式中：df和dg分別定義為

(17)

(18)

在知道屈服面的情況下，分?jǐn)?shù)階塑性力學(xué)模型無需再額外假定一個塑性勢面函數(shù)，而是通過分?jǐn)?shù)階微分直接得到一個非關(guān)聯(lián)的流動法則，從而使模型得到簡化。由文獻(xiàn)[15]可知，當(dāng)α越大時，所得到的非關(guān)聯(lián)性越明顯。具體參見圖1，當(dāng)α由1.05增加到1.3時，同一應(yīng)力水平下加載方向與流動方向的夾角變大，即非關(guān)聯(lián)程度增加。

圖1 分?jǐn)?shù)階數(shù)對非關(guān)聯(lián)性的影響Fig.1 Effect of fractional order on

3 三維分?jǐn)?shù)階塑性力學(xué)模型

3.1 屈服面

由傳統(tǒng)塑性力學(xué)可知，加載方向是屈服面的法向。選取特征應(yīng)力條件下的修正劍橋模型屈服面函數(shù)來描述土體的受荷屈服。

(19)

(20)

式中：φc是常規(guī)三軸壓縮時的臨界狀態(tài)摩擦角。如圖2所示，特征應(yīng)力法修正后的屈服面在常規(guī)主應(yīng)力空間的屈服面不再是一個圓錐，而是一個在π平面的投影，是曲邊三角形的曲面，如圖3所示。因此，可以反映土體在真三軸應(yīng)力條件下的屈服特性。

圖2 主應(yīng)力空間的三維屈服面(β=0.4)Fig.2 3D yielding surface in principal stress space

圖3 π平面上的屈服面(β=0.2)Fig.3 Yielding surface on the π plane

3.2 加載張量與流動張量

將式(19)分別代入式(17)、式(18)，進(jìn)一步求得df和dg為

(21)

(22)

將式(21)和式(22)代回至式(15)和式(16)，即可得到塑性加載張量和塑性流動張量的具體表達(dá)式

(23)

(24)

比較式(23)和式(24)可知，當(dāng)α=1時，模型退化為傳統(tǒng)關(guān)聯(lián)流動法則。由圖4可知，對屈服面進(jìn)行分?jǐn)?shù)階求導(dǎo)后，所得子午平面上塑性流動方向不再垂直于屈服面，而是與加載方向存在一個夾角。

圖4 子午平面上的屈服面Fig.4 Yielding surface on the meridian

此外，為了可以反映受荷變形的狀態(tài)依賴特性，Sun等[14-15]建議采用狀態(tài)依賴的分?jǐn)?shù)階數(shù)

α=exp(Δψ)

(25)

式中：Δ>0，是材料常數(shù)。式(25)是經(jīng)驗公式，來源于基于試驗現(xiàn)象的考慮：土體所處狀態(tài)會影響其塑性流動特性，而分?jǐn)?shù)階數(shù)決定了塑性流動方向。因此，運(yùn)用式(25)直觀地考慮了加載過程中土體所處狀態(tài)對其塑性流動方向的影響。ψ是狀態(tài)參數(shù)，可以定義為

ψ=e-ec

(26)

式中：ec是臨界狀態(tài)孔隙比，定義為[25]

(27)

式中：eΓ、λ和ξ是e-p′空間的臨界狀態(tài)參數(shù)。

3.3 硬化模量

硬化模量描述了土體的硬化特性。Sun等[14]采用Li等[25]提出的常規(guī)三軸應(yīng)力條件下的硬化模量。為了使建立的模型可以反映真三維應(yīng)力條件下的硬化特性，基于特征應(yīng)力法修正了Li等[25]的硬化模量，見式(28)。

(28)

(29)

4 模型參數(shù)識別與模型驗證

4.1 模型參數(shù)識別

模型含有11個參數(shù)，分別為：φc、λ、eΓ、ξ、Δ、β、k、h1、h2、G0以及ν。這些模型參數(shù)均可以通過三軸試驗獲得。

參數(shù)φc、λ、eΓ和ξ為臨界狀態(tài)參數(shù)。φc可以通過擬合p′-q平面臨界狀態(tài)線的斜率得到；λ、eΓ、ξ可以通過擬合e-p′空間的臨界狀態(tài)數(shù)據(jù)得到[5]。圖5所示為臨界狀態(tài)參數(shù)對模型預(yù)測結(jié)果的影響，可見，隨著φc和eΓ數(shù)值的增加，模型預(yù)測得到的同一應(yīng)變水平下的應(yīng)力比逐漸增加；隨著λ的增加，模型預(yù)測的同一應(yīng)變水平下的應(yīng)力比逐漸減小。

圖5 臨界狀態(tài)參數(shù)對模型預(yù)測的影響Fig.5 Effect of critical-state parameters on model

參數(shù)Δ決定了土體的塑性流動特性，可以由試樣狀態(tài)轉(zhuǎn)換時的剪脹關(guān)系得到，即d=0。將式(25)代入式(22)，得

(30)

因此

(31)

(32)

再將式(29)代入式(32),得

(33)

因此

(34)

參數(shù)Δ和k對模型預(yù)測的影響見圖6。分析發(fā)現(xiàn)：隨著參數(shù)Δ的增加，模型預(yù)測得到的同一應(yīng)變水平下的應(yīng)力比逐漸增加；隨著參數(shù)k的增加，模型預(yù)測得到的同一應(yīng)變水平下的應(yīng)力比逐漸降低。

圖6 Δ和k對模型預(yù)測的影響Fig.6 Effects of Δ and k on model

特征應(yīng)力參數(shù)β可由土體相關(guān)三軸壓縮(φc)和拉伸(φe)時的臨界狀態(tài)摩擦角關(guān)系式獲得[21]。β≤0.1時，模型在π平面的屈服軌跡非常接近于Matsuoka-Nakai準(zhǔn)則[26]。因此，β近似取0.1。

參數(shù)h1和h2可以通過擬合土體受荷時ε1-q關(guān)系曲線獲得。具體參見文獻(xiàn)[15,25]，這里不再贅述。圖7所示為參數(shù)h1和h2對模型預(yù)測結(jié)果的影響，分析可知：當(dāng)參數(shù)h1增加時，同一應(yīng)變水平下模型預(yù)測的應(yīng)力比增加；當(dāng)參數(shù)h2增加時，同一應(yīng)變水平下模型預(yù)測的應(yīng)力比減小。

圖7 h1和h2對模型預(yù)測的影響Fig.7 Effects of h1 and h2 on model

兩個彈性參數(shù)(G0、ν)可以通過擬合堆石料初始加載階段的應(yīng)力-應(yīng)變曲線獲得，具體參見文獻(xiàn)[27]。由于受荷后堆石料彈性變形相對塑性變形較小，因此，不再探討彈性參數(shù)對模型預(yù)測的影響。

4.2 模型驗證

選取不同初始狀態(tài)的堆石料[28]在真三軸壓縮條件下的試驗數(shù)據(jù)進(jìn)行擬合和對比，以驗證所提出三維分?jǐn)?shù)階巖土塑性力學(xué)模型的合理性。

所選堆石料為兩河口心墻土石壩筑壩料，試驗初始孔隙比為0.68。試驗時固定第三主應(yīng)力和中主應(yīng)力比不變，第二和第一主應(yīng)力按照設(shè)定中主應(yīng)力比值的大小比例增加。其他有關(guān)材料和試樣準(zhǔn)備資料見文獻(xiàn)[28]，不再贅述。根據(jù)前述方法，確定具體模型參數(shù)為：G0=90、ν=0.25、φc=46°、λ=0.11、eΓ=0.404、ξ=0.1、Δ=0.2、β=0.1、k=0.1、h1=1.2、h2=0.3。

圖8～圖10為所提出分?jǐn)?shù)階塑性力學(xué)模型對真三軸試驗數(shù)據(jù)的模擬。對比分析可以發(fā)現(xiàn)：模型可以合理地模擬堆石料在不同圍壓和中主應(yīng)力比條件下的真三維應(yīng)力-應(yīng)變特性。具體來說，隨著第一主應(yīng)變的增加，試樣的應(yīng)力比不斷增加；第二主應(yīng)變只有在中主應(yīng)力比等于0時，即常規(guī)三軸試驗條件下，數(shù)值為負(fù)，其他中主應(yīng)力比條件下均為正；第三主應(yīng)變在任何中主應(yīng)力比條件下均為負(fù)值。這些應(yīng)力-應(yīng)變特性均可被較好地模擬。

圖8 圍壓100 kPa時三維應(yīng)力-應(yīng)變關(guān)系的模型預(yù)測Fig.8 Model predictions of the 3D stress-strain behavior of rockfill at confining pressure of 100

此外，需要指出的是，中主應(yīng)力系數(shù)較小時，圖8～圖10中試樣在同一偏應(yīng)力水平所產(chǎn)生的第一主應(yīng)變大于第二、第三主應(yīng)變，但隨著中主應(yīng)力系數(shù)的增加，第一主應(yīng)變小于第二、第三主應(yīng)變；這是因為隨著中主應(yīng)力系數(shù)從0變化到1，原先的第二主應(yīng)力大小慢慢接近于第一主應(yīng)力，對側(cè)向應(yīng)變的主導(dǎo)作用逐漸變大。從而在第一主應(yīng)力方向產(chǎn)生的應(yīng)變會逐漸小于同等偏應(yīng)力水平下其他兩個主應(yīng)力方向的應(yīng)變。

圖9 圍壓150 kPa時三維應(yīng)力-應(yīng)變關(guān)系的模型預(yù)測Fig.9 Model predictions of the 3D stress-strain behavior of rockfill at confining pressure of 150

圖10 圍壓200 kPa時三維應(yīng)力-應(yīng)變關(guān)系的模型預(yù)測Fig.10 Model predictions of the 3D stress-strain behavior of rockfill at confining pressure of 200

5 結(jié)論

基于特征應(yīng)力法，修正了已提出的分?jǐn)?shù)階應(yīng)力塑性力學(xué)模型，使其可以對堆石料在一般三維應(yīng)力條件下的應(yīng)力-應(yīng)變進(jìn)行預(yù)測。并選擇堆石料的真三軸壓縮試驗數(shù)據(jù)，對模型進(jìn)了驗證。主要結(jié)論如下：

1)基于特征應(yīng)力修正的屈服面，其在主應(yīng)力空間不再是一個圓錐，而是一個可以反映土體三維屈服特性的曲面，該面在π平面的投影為曲邊三角形。

2)無需額外假定塑性勢面，只需要對已有屈服面進(jìn)行分?jǐn)?shù)階微分求導(dǎo)，便可以得到非關(guān)聯(lián)的流動法則。進(jìn)一步得出剪脹方程受分?jǐn)?shù)階求導(dǎo)階數(shù)、特征應(yīng)力參數(shù)以及中主應(yīng)力系數(shù)的影響。

3)基于特征應(yīng)力法修正了已有常規(guī)三軸應(yīng)力條件下的硬化模量，使其可以反映三維應(yīng)力條件下堆石料的硬化特性。

4)對比試驗數(shù)據(jù)，驗證了所提模型的合理性；模型可以對不同初始狀態(tài)試樣在不同中主應(yīng)力比條件下的真三軸應(yīng)力-應(yīng)變特性進(jìn)行較好的模擬。

5)建立模型僅考慮了完全飽和堆石料的變形特性，對于干燥或半干燥堆石料遇水飽和后存在的濕化或者軟化性能還有待進(jìn)一步研究。此外，模型也未考慮堆石料在填筑或者運(yùn)維過程中的顆粒破碎，對于考慮顆粒破碎影響的堆石料三維分?jǐn)?shù)階塑性模型也有待進(jìn)一步研究。