蔣明鏡 ，盧國文 ，李 濤

(1.天津大學建筑工程學院，天津 300072；2.水利工程仿真與安全國家重點實驗室(天津大學)，天津300072；3.同濟大學土木工程防災國家重點實驗室，上海 200092；4.同濟大學土木工程學院，上海 200092)

在經(jīng)典土力學中，土體本構模型能很好地反映土體的力學性質，并能將土體的強度問題和變形問題融合，以劍橋模型為代表的彈塑性模型的發(fā)展標志著人們認識土體特性的一次飛躍.但這一時期多數(shù)本構模型都是針對飽和重塑土發(fā)展起來的，在實際工程中遇到的原狀土大都具有結構性[1].其不僅與受力歷史和土體密度有關，還取決于土粒間膠結物的力學特性.

在過去幾十年中，研究人員利用多種方法對結構性黃土進行了數(shù)學建模，早期發(fā)展的是微觀機制模型，其中應用較多的一種為微結構模型，其實質是對土體顆粒及孔隙的排列、形狀、接觸關系的類型劃分和數(shù)學建模.Rowe[2]認為宏觀塑性應變是在滑移面上所產生的滑移的累加，并在此基礎上建立了微觀力學模型.沈珠江等[3]率先把損傷理論引入到土體的本構模型研究中，認為天然土的結構破損是從原狀土逐漸向擾動土(損傷土)變化的過程.隨后建立了彈塑性損傷模型[4]、非線性損傷模型[5]、堆砌體模型[6]，并在此建模思路上開始向巖土破損力學[7]過渡，建立并完善了二元介質類模型.在結構性黃土研究方面，沈珠江[8]較早提出了一個可以適用于黃土的損傷力學模型.Desai[9]于1974年提出擾動狀態(tài)概念(disturbed state concept)，此后 Desai 等[10]、吳剛[11]在擾動狀態(tài)概念的基礎上建立一系列描述巖土材料力學特性的本構模型，統(tǒng)稱為擾動狀態(tài)模型.與上述建模思路不同，許多研究人員提出表達土體結構性強弱的宏觀定量化參數(shù)來建立能夠描述土體變形和強度規(guī)律的數(shù)學模型.近年來許多研究人員在臨界狀態(tài)土力學基礎上建立適用于結構性土的本構模型[12-13]，稱之為臨界狀態(tài)模型.在黃土本構模型研究方面，陳正漢等[14]將該關系式引入到彈塑性模型中，給出了濕陷變形的計算方法，得到了濕陷初始面，但模型參數(shù)確定需要進行較多的試驗，且只能計算增濕到飽和狀態(tài)的變形.Jiang等[15]基于沈珠江提出的廣義吸力的概念提出了可描述結構性黏土的逐漸破損過程的結構性黏土的結構吸力模型.謝定義等[16]提出了結構性定量化指標——綜合結構勢.

鑒于傳統(tǒng)方法建立的本構模型中的一些關鍵變量缺乏微細觀力學機制的支撐，基于宏微觀土力學的研究思路[17]，蔣明鏡等[18]通過二維離散元法驗證了膠結破損規(guī)律，并建立了天然結構性土的本構模型.張伏光等[19]通過三維離散元法建立了膠結砂土的本構模型.筆者認為，結構性黃土的研究需要同時考慮黃土的非飽和性.然而非飽和土本構模型非常復雜，基質吸力對土體抗剪強度的影響呈非線性，另外，基質吸力對非飽和土壓縮性的影響以及土水特性曲線的引入都加大了建模難度.為了簡化建模過程，李廣信等[20]曾對非飽和土的實用化模型展開過探索和研究.

為此，本文以巖土破損力學及臨界狀態(tài)土力學為框架，基于實用化宏微觀土力學研究思路，通過建立非飽和吸應力與飽和度的關系式，將臨界狀態(tài)線、擴展橢圓屈服面和臨界 Lade-Duncan強度準則推廣到非飽和土體，采用三維離散元模擬分析表征結構性損傷的膠結破損規(guī)律及其參數(shù)建議公式，然后將其引入到硬化規(guī)律中，建立了基于經(jīng)典彈塑性模型的非飽和黃土三維本構模型，并通過黃土室內試驗進行了驗證.

1 本構模型

1.1 有效應力原理

對于飽和土，太沙基提出了有效應力的概念并在此基礎上建立了古典土力學.對于非飽和土，有效應力(沈珠江[21]稱為廣義有效應力)為

為了模型的實用性，本文僅引入飽和度，不引入吸力和土水特性曲線，用飽和度計算吸應力ps.本文對不同含水量結構性黃土離散元試樣進行了離散元模擬分析，采用冪函數(shù)描述吸應力和初始飽和度的關系以反映初始飽和度對吸應力的影響，如式(2)所示.

式中：ξ表示土體非飽和程度；cw和nw為擬合參數(shù)；Sr為土體的飽和度.

如果考慮某些試驗中飽和重塑土也存在黏聚力(如式(3)所示)，可以將吸應力計算式延伸為

式中ps0為飽和土的吸應力，根據(jù)飽和重塑土或結構性土臨界狀態(tài)強度線求取.

1.2 臨界狀態(tài)強度面

圖 1給出了結構性黃土離散元試樣真三軸試驗π平面強度(數(shù)據(jù)點)與 Lade-Duncan準則對比[22](曲線)，該模擬在軟件 PFC中實現(xiàn)，主要步驟為離散元試樣制備、預壓、真三軸試驗模擬.在使用 Lade-Duncan準則時，臨界狀態(tài)應力比Mcr取相應飽和度常規(guī)三軸試驗結果，峰值強度為Mf=qf/p，qf為相應飽和度常規(guī)三軸試驗峰值強度包面上與p′相對應的縱坐標(偏應力)值.從圖 1中可以看出，在使用有效應力的前提下，Lade-Duncan準則能夠較好反映非飽和結構性黃土離散元試樣的π平面強度.因此，采用有效應力概念后，Lade-Duncan準則可以非常方便地將應用范圍從無黏性土擴展到有黏聚力的非飽和土體.

鑒于此，本文在建立非飽和結構性黃土本構模型時采用Lade-Duncan準則作為臨界狀態(tài)強度.其強度包面可表示為

式中：θσ為應力羅德角；M(θσ)為隨θσ變化的土體強度線斜率.從圖中也可以看出，Lade-Duncan考慮了中主應力與羅德角的影響，不同主應力系數(shù)b或羅德角θσ條件下的M(θσ)可由中主應力系數(shù)為0(三軸壓縮條件下θσ＝-π/6)時的強度線斜率Mcr求得，具體計算方法見文獻[23].

1.3 屈服函數(shù)

橢圓屈服面在應用于結構性土體時，屈服面“干側”不能很好地描述室內試驗情況[24].依沈珠江[25]關于破壞準則和屈服函數(shù)的總結，采用封閉型函數(shù)拓展橢圓屈服面，即

式中：ns為形狀參數(shù)；η′為有效應力比；為結構性土屈服面在p軸上的截距.

將式(6)變形可得式(7)，即

圖 2給出了不同ns的擴展橢圓屈服面形狀，從圖中可以看出，當ns＝2時，式(7)退化為修正劍橋模型的橢圓屈服面；當ns＜2時，屈服面峰值左移；當ns＞2時，屈服面峰值右移.一般說來，對于重塑和結構性土體，應取ns＞2.

圖2 擴展橢圓屈服面形狀Fig.2 Shape of the extended elliptical yield surface

1.4 破損參數(shù)

鑒于室內試驗無法直接獲取土體破損演化規(guī)律，一般先假設后通過試驗的方法間接確定，缺乏明確的物理意義.本文在巖土破損力學的框架內，結構性土顆粒可看作是由無膠結顆粒集合體和膠結顆粒集合體組成，根據(jù)均質化理論，引入體積破損率λ概念，有如下關系[26].

結構性土的屈服面與對應的重塑土的屈服面大小不同，但形狀相似[18].如圖 3所示，圖中pc為重塑土屈服面在p軸上的截距.

從原點引出任意一條直線分別交重塑土和結構性土屈服面于A、B點，兩點的應力狀態(tài)分別可以表示為.宏觀上結構性土強度組成和微觀上無膠結和膠結部分應力分擔之間存在以下關系[25].

在應變分擔二元介質模型中有

綜合式(10)、(11)可得

式中：λs為基于修正彈塑性方法定義的結構破損參數(shù)，其物理意義為代表性單元內無膠結顆粒的局部應力與代表性單元的平均應力的比值，λs會有一個初值(λs初值不能為零)，表示結構性的強弱，隨著結構性土的受荷或增濕過程，會不斷增加，其極限值為1，此時結構性土完全演化為重塑土.

基于宏微觀土力學的研究思路，建立非飽和結構性黃土本構模型的關鍵為確定合理的破損參數(shù)演變規(guī)律，其需要從微觀尺度對土體結構性破損演變規(guī)律進行定量分析，從而建立其與宏觀力學參量間的關系.

基于二維離散元模擬結果，孫渝剛[26]選擇建立大主應變與破損參數(shù)的關系，劉靜德[27]選擇建立等效塑性應變(包含剪切應變和體應變的復合作用)與破損參數(shù)的關系.由于結構破損參數(shù)與塑性應變的關系式中需要同時考慮塑性體應變和塑性偏應變的影響，用某種形式的等效塑性應變來描述結構破損參數(shù)演化更具有普適性.因此，本文建議將結構破損參數(shù)的演化建立在等效塑性應變之上，以上為一個關鍵點.對于非飽和結構性黃土而言，怎么考慮含水量變化對結構破損參數(shù)演化規(guī)律的影響是結構性黃本構模型構建的另一個關鍵點.

通過室內試驗猜想假設[25]和非飽和結構性黃土三維離散元模擬分析可知，該離散元試樣的結構破損演化速率跟含水量相關，含水量越低，結構破損速率越慢.因此，本文定義等效塑性應變與結構性土的結構屈服應力(等向壓縮試驗測得)的比值為等效塑性應變系數(shù)，即有為等效塑性應變.為了拓展模型的適用性，增加參數(shù)k將等效應塑形應變定義為當k＝1時，退化為.根據(jù)室內研究結構屈服應力σy′可以采用式(13)求解.

式中cy1與cy2為模型參數(shù).

經(jīng)過加荷增濕等試樣離散元模擬，并分析試驗過程中的土體膠結破壞過程，本文建議采用以下函數(shù)表征土體結構膠結破損規(guī)律.

令Ep＝0時的λs為其初始值λs0，則可得到cb表達式，將其代入式(14)得

式中：νs＝ca，表示結構破損演化的速率；λs0為結構破損參數(shù)初始值.

1.5 流動法則

基于修正劍橋模型剪脹率-應力關系，劉靜德[27]參考 Li等[28]研究思路，建立了考慮土體狀態(tài)參數(shù)ψ的剪脹率-應力關系，即

式中：Md為反映土體由剪縮到剪脹的特征狀態(tài)應力比；ψ＝e-ecr為狀態(tài)參數(shù)，ecr為臨界狀態(tài)孔隙比；n為模型參數(shù).特別注意，雖然本文采用的剪脹率-應力關系表達式與前人相同，但是具體實施時需要采用非飽和土的有效應力比和臨界狀態(tài)孔隙比，比如有效應力比η′＝q/p′(式中，)考慮了非飽和土的吸力造成的黏聚力對塑性流動的影響，非飽和土的臨界狀態(tài)孔隙比也需要考慮吸力對臨界狀態(tài)壓縮曲線的影響.

對剪脹率-應力關系式(16)進行積分即可得到塑性勢函數(shù).

1.6 硬化規(guī)律

在修正劍橋模型中，飽和土的正常固結線(NCL)假設為直線.

式中：N為參考應力(1 kPa)對應的孔隙比；λ為土體等向壓縮試驗所得的壓縮指數(shù).

根據(jù) Hu等[29]研究表明，式(21)能很好地描述非飽和土的正常固結線.

式中ah和bh為模型參數(shù).

根據(jù)式(21)，可求出非飽和土的cp′(ξ)[30]，即

為了考慮土體狀態(tài)參數(shù)ψ對硬化規(guī)律的影響，劉靜德[27]在式(22)中引入了一個乘數(shù)，cp′(0)表達式改進為

式中：Mb＝Mexp(-mψ)為土的峰值應力比[28]，m為模型參數(shù)；cp＝(λ-k)/(1+N).

式(12)、(15)、(21)～(23)聯(lián)立，可得非飽和黃土的硬化規(guī)律表達式為

式中cp′即為cp′(ξ)，表達式經(jīng)聯(lián)立即可求得，由于較為繁瑣，此處不再給出.

1.7 彈塑性矩陣

根據(jù)Hook定律，應力-應變增量關系表示為

根據(jù)屈服函數(shù)f＝0，可得一致性條件為

將式(26)代入式(25)，可得塑性因子，即

其中

式中A為塑性硬化模量.將塑性因子帶入式(25)可得

式中

2 模型參數(shù)確定

下面將模型參數(shù)分為以下幾類，給出了所有模型參數(shù)的物理意義，并介紹了采用室內試驗確定模型參數(shù)的方法.

(1)有效應力參數(shù)：ps0、cw、nw.

有效應力參數(shù)與土體非飽和性相關，用來計算非飽和土的吸應力，從而計算土體的有效應力，ps0一般應取為0，即飽和土的吸應力為零(黏聚力為零).

通過不同初始飽和度的重塑土等含水量三軸試驗確定臨界狀態(tài)強度包線(p-q平面)，各狀態(tài)強度包線與橫軸的截距即為不同初始飽和度的吸應力，然后通過擬合獲得吸應力與初始飽和度的關系獲得有效應力參數(shù)ps0、cw、nw.

(2)彈性參數(shù)：κ、ν.

κ由飽和重塑土等向壓縮回彈試驗結果在半對數(shù)e-lnp′平面內的回彈曲線求得；泊松比ν可按經(jīng)驗選取.

(3)屈服面形狀參數(shù)：ns.

ns可通過擬合不同應力路徑下的屈服點，根據(jù)屈服點的位置確定屈服面而獲得，一般建議取ns≥2，不方便確定時可以直接選定ns＝2～3.

(4)臨界狀態(tài)參數(shù)：Mcr、λ、N、ah、bh、m、n.

Mcr、λ、N為修正劍橋模型固有參數(shù)，ah、bh是為反映非飽和土正常固結線上移而引入的參數(shù).

λ可由飽和重塑土等向壓縮試驗結果在e-lnp′平面內的壓縮曲線求得；N為飽和重塑土等向壓縮試驗的正常固結曲線上參考應力(1kPa)所對應的孔隙比，由常規(guī)三軸試驗的臨界狀態(tài)應力比求得，可對不同含水量試驗結果取平均值.Mcr可依據(jù)不同初始飽和度重塑土等向壓縮試驗，ah、bh通過擬合不同飽和度(ξ＝1-Sr)在不同平均有效應力下的孔隙比比值eN/esN的關系獲得.

引入?yún)?shù)m、n反映土體狀態(tài)參數(shù)ψ對硬化規(guī)律的影響，求解ψ需知道臨界狀態(tài)孔隙比ecr.m、n、ecr確定方法可參見文獻[27]，特別注意其主要應用于砂土，對于粉土和黏土(包括黃土)，一般取m＝n＝0.

(5)結構性參數(shù)：cy1、cy2、νs、k.

通過不同初始飽和度的結構性土等向壓縮曲線可以獲得結構屈服應力σy′與初始飽和度的關系，從而通過擬合獲得cy1、cy2.破損參數(shù)νs可以通過試算結構性土的等向壓縮線確定，參數(shù)k可通過不同應力比的等應力比壓縮試驗確定，也可以假設為 1，建議0≤k≤1.

3 模型驗證

為驗證模型的正確性及有效性，選取天然非飽和結構性黃土及人工制備結構性黃土進行模擬，并與室內試驗進行對比分析.

圖4給出了Q3天然非飽和結構性黃土(數(shù)據(jù)來源：文獻[31]，取土地點：西安理工大學曲江校區(qū)附近)三軸試驗偏應力-軸向應變關系和體應變-軸向應變關系實測結果與本文離散元模擬結果(不同初始飽和度下的等含水量(排氣不排水)及相應離散元試樣).采用前文提到的參數(shù)確定方法，模擬使用的參數(shù)見表1.初始飽和度Sr＝0.362編號為#1，Sr＝0.459編號為#2，三軸試驗圍壓為 100kPa、200kPa、300kPa和 400kPa.其中，cw和nw為根據(jù)不同飽和度三軸試驗臨界狀態(tài)強度線的橫軸截距(吸應力)計算；ah和bh需要根據(jù)不同飽和度重塑土等向壓縮曲線計算，本文根據(jù)不同飽和度結構性黃土壓縮曲線估算.

從圖 4中可以看出，由于引入了參數(shù)反映飽和度對吸應力ps和結構屈服應力σy′的影響，模型能夠基本反映不同飽和度下結構性黃土在三軸試驗條件下的力學特征，隨著飽和度的增加，相同圍壓下的偏應力下降.由于結構性的影響，試樣在低圍壓下可以表現(xiàn)出較高的剪切模量，試樣傾向于弱硬化甚至軟化，隨著圍壓的增加試樣在剪切開始就可能出現(xiàn)初始屈服，試樣呈現(xiàn)硬化.

圖 5給出了人工制備結構性黃土(數(shù)據(jù)來源：非飽和人工制備結構性黃土[32]的側限壓縮試驗和17.7%含水量黃土等含水量三軸試驗的試驗結果和模擬結果.采用前文提到的參數(shù)確定方法，模擬使用的參數(shù)見表 1，試驗編號為#3，三軸試驗圍壓 5 kPa、25kPa、200kPa 和 300kPa.其中，cw、nw、ah、bh選取方法與前述一致，但因為僅有一個飽和度試驗，選定nw＝1.

從圖 5中可以看出，模型能夠基本反映不同飽和度人工制備結構性黃土的壓縮特征，對于飽和人工結構性黃土，土體在較低壓力下屈服，對于 17.7%含水量結構性黃土，試樣在較高壓力下屈服，屈服以后試樣的壓縮線斜率變陡.17.7%含水量的平均有效應力較大，因此彈性模量比飽和土模擬結果大，這與試驗結果相符.對于三軸試驗，與天然非飽和結構性黃土試驗結論相同，試樣在低圍壓下呈現(xiàn)出快速到達峰值后逐步軟化，在高圍壓下試樣較早屈服后逐漸硬化，這些主要特點都能夠為模型所反映.

圖4 非飽和結構性黃土等含水量三軸試驗與模擬結果Fig.4 Predicted and experimental results in triaxial compression test on natural loess with constant water content

表1 模擬使用參數(shù)Tab.1 Model parameters used in verification tests

圖5 非飽和人工制備結構性黃土側限及17.7%含水量三軸試驗與模擬結果Fig.5 Predicted and experimental results in oedometer and triaxial compression tests on artificial loess with 17.7% water content

4 結 論

本文在巖土破損力學及臨界狀態(tài)土力學框架內，通過膠結材料微觀力學理論及三維離散元模擬分析，給出了表征結構性損傷的膠結破損規(guī)律，將其引入到硬化規(guī)律中，建立了基于經(jīng)典彈塑性模型的非飽和黃土本構模型.主要得到以下結論.

(1)通過建立非飽和土吸應力與飽和度的關系得到了非飽和土有效應力表達式，可以方便地將臨界狀態(tài)強度線、擴展橢圓屈服面和臨界狀態(tài) Lade-Duncan準則推廣到非飽和土體.

(2)定義等效塑性應變系數(shù)以考慮含水量對結構破損演化參數(shù)規(guī)律的影響，膠結破損參數(shù)表達式物理意義明確，通過離散元驗證分析，能夠很好地反映土體結構性損傷微觀機理.

(3)本文模型參數(shù)均可通過室內試驗確定，Q3天然非飽和結構性黃土三軸試驗及人工制備結構性黃土側限壓縮試驗和三軸試驗結果與本文模擬結果較為吻合，表明該模型是合理可行的.