馬田田，韋昌富，陳 盼，魏厚振，伊盼盼

（1.中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071；2. 武昌理工學院，武漢 430223）

1 引 言

在分析非飽和土問題時，通常將水力學和力學性質分開考慮。但試驗現(xiàn)象表明，非飽和土中滲流與變形存在復雜的耦合作用[1－2]。滲流過程影響土體的強度與變形，例如反復的干、濕循環(huán)會使非飽和土的強度降低[3]。同時，土體的變形又會改變非飽和土的土-水特性，從而影響孔隙水的滲流過程[4]。因此，有效地模擬任意含水率變化路徑下非飽和土的滲流、變形及強度變化之間的耦合作用是亟需解決的關鍵問題。要解決該關鍵問題首先需要建立兩個本構方程，即應力-應變關系和土-水特征關系及其之間的相互影響。

Alonso等[5]在修正劍橋模型的基礎上，提出BBM（Barcelona Basic Model）[6]模型。該模型采用雙應力變量，引入 LC（Loading-Collapse）屈服面概念來考慮吸力對屈服應力的硬化作用，能夠模擬非飽和土吸水濕陷現(xiàn)象。采用 SI、SD屈服面來描述飽和度的塑性屈服。Wheeler等[7]用吸力控制的非飽和土等向壓縮和三軸試驗數據驗證了BBM模型的合理性。Wheeler[8]、Nuth[9]、孫德安[10]、繆林昌[11]和盛岱超等[12]在BBM模型基礎上，采用有效應力代替凈應力作為應力狀態(tài)變量，也用LC屈服面描述非飽和土在不同應力下的吸水濕陷現(xiàn)象。

Wheeler等[7]指出，BBM等非耦合本構模型不能預測在常吸力下飽和度隨著剪切而發(fā)生變化的特性，也不能反映飽和度隨變形的發(fā)展增大的這種耦合特性。因此，不能準確地預測不排水條件下非飽和土的水力和力學性狀。BBM模型及其以此為基礎發(fā)展的模型，以塑性體變?yōu)橛不瘏?，因此，不能反映當塑性體變不變時，吸濕即基質吸力減小使屈服面縮小的現(xiàn)象[8]。

在描述非飽和土的本構行為時，目前最流行的做法是基于“LC曲線”提出的模型，用來反映屈服應力隨吸力的硬化作用。為了反映土體塑性變形中的非飽和效應，LC曲線（加載-崩塌）已成為許多傳統(tǒng)非飽和土本構模型的核心組成部分[2,5－7,13]。

但研究發(fā)現(xiàn)，基于LC屈服面的模型有以下局限性[14]：①該類模型在含水率循環(huán)變化時不適用，如在飽和與非飽和狀態(tài)的轉換區(qū)域附近，如果土體處于脫濕過程且基質吸力小于進氣值，可視為飽和土；如果土體處于吸濕過程，則表現(xiàn)出明顯的非飽和特性；②LC曲線只是描述屈服強度隨吸力增加而增大的性質，而沒有引入飽和度作為基本變量，因此，不能考慮前期含水率及其變化歷史的影響[8]。由于存在毛細滯回效應，即使在吸力相同的情況下，飽和度并不相同，土體表現(xiàn)出完全不同的力學特性；③模型中沒有考慮毛細循環(huán)滯回特性；④沒有考慮飽和度及其循環(huán)變化和土體變形及強度變化之間的耦合效應[8]。

目前主要有兩種方式考慮吸力對前期固結壓力的硬化作用，韋昌富[14]、Khalili[15]和 Tamagnini[16]等在原本的硬化函數基礎上乘以一個受吸力影響的系數，此種方法可以描述吸力的影響與體積硬化作用兩者的耦合效應，其中韋昌富[14]是將吸力改為飽和度，未考慮吸力的影響；Jommi[17]、Li[18]則認為，吸力對土體應變硬化有一個附加作用，將吸力的影響與體變硬化的作用分別考慮，總的硬化函數等于兩者之和。

Loret等[19]分析試驗數據提出以式（1）表示吸力對前期固結壓力的硬化作用，

式中：pc為非飽和土的屈服應力；pc0為飽和土的初始屈服應力；X為硬化參數；ξ為系數；pcm(Sc)、pca(Sc)為吸力對前期屈服應力的貢獻，Sc為基質吸力； pcm(Sc)=1為體應變硬化與吸力影響之間是解耦的， pca(Sc)=1為體應變硬化與吸力影響是耦合的。但上述模型都沒有考慮飽和度和基質吸力同時對初始屈服應力的貢獻。

Wheeler[20]、Mualem[21]通過試驗發(fā)現(xiàn)，飽和度對非飽和土的應力-應變關系有很大的影響，建議在非飽和土彈塑性本構模型中引入循環(huán)土-水特征曲線，來描述含水率循環(huán)變化對非飽和土變形的影響。

在描述土-水特征曲線時，所見模型主要分為4種：經驗模型、區(qū)域模型[21]、理性外推模型[22]和邊界面模型[23]。盡管這些模型不同程度地考慮了毛細滯回現(xiàn)象，但很少反映土體變形對持水特征的影響。繆林昌[24]、方祥衛(wèi)[25]、龔壁衛(wèi)等[26]通過試驗研究了土體的密實狀態(tài)對土-水特征關系的影響，Gallipoli[27]通過試驗提出了土體體變與 van Genuchten模型[28]參數的關系，因此，可以考慮體變對土-水特征曲線的影響?；谝陨洗嬖诘膯栴}，本文將吸力和飽和度同時作為硬化參數，在同一個框架中建立了塑性變形和毛細滯回耦合的本構模型。

2 非飽和土毛細滯回與骨架變形耦合本構模型

2.1 應力狀態(tài)變量

應力狀態(tài)變量采用Bishop[29]有效應力σi′j和基質吸力Sc：

在三軸應力狀態(tài)下：

式中：σij為總應力張量；σi′j為有效應力張量；Sr為飽和度；ua、uw分別為孔隙氣壓力和孔隙水壓力；δij為單位張量。式（1）是將Bishop有效應力[29]中的參數χ取為飽和度Sr[15]。Houlsby[30]提出與此應力狀態(tài)變量（σi′j，Sc）功共軛的應變?yōu)椋é舏j，Sr），其中εij為應變張量。p′為平均有效應力；q為廣義剪應力；p為平均總應力。

2.2 基質吸力和飽和度對屈服應力的作用

在達到相同吸力的條件下，分別通過吸濕和脫濕之后的土樣含水率的大小、分布可能并不相同，如圖1所示。從非飽和土的CT掃描圖[31－32]可以看出，相同含水率條件下，孔隙水的分布狀態(tài)是不同的，經歷吸濕過程的土樣比經歷脫濕過程的含水率分布更加均勻，因此，非飽和土的力學性質也不同。通過測量兩種狀態(tài)試樣的壓縮波傳播時發(fā)現(xiàn)，前者的壓縮波速明顯比后者小[31,33－34]。

圖1 非飽和土的CT掃描圖（飽和度均為92%時分別經過吸濕和脫濕而成的土樣）Fig.1 CT images of unsaturated soil(degree of saturation is 92%, after drying and wetting respectively)

從上面的分析可知，由于毛細滯回的影響，非飽和土體的土-水狀態(tài)與土體所經歷的水力狀態(tài)密切相關，而土-水狀態(tài)的變化必將影響土體的強度與變形。因此，必須要同時考慮基質吸力和飽和度對前期固結壓力的影響。通過試驗[35]發(fā)現(xiàn)，屈服應力隨著吸力的增加、飽和度的減小而增大，因此，假設用式（6）來考慮吸力和飽和度對屈服應力的硬化作用，

采用類似Alonso[5]提出的壓縮曲線的斜率與吸力的關系，定義h函數為

式中：r為非飽和土處于殘余飽和度時的強度與飽和土屈服強度的比值；m為屈服應力隨飽和度和吸力增長的速度；為最大的塑性體積應變；為殘余飽和度；Sb為非飽和土的進氣值； Sc/Sb為吸力比。〈 〉為Macauley大括號，定義為 x = xH(x)，其中H(x)為Heaviside公式。

式（7）滿足以下條件：①當 Sr= 100%和Sc=0時，h=1.0；②隨著Sr的減小，Sc的增大而增大，并趨于一個穩(wěn)定值；③隨著的增大，飽和度和基質吸力對強度的影響逐漸減弱，當趨于最大值即時，飽和度和基質吸力不再對強度有影響。

本文采用修正劍橋模型[36]，如圖2所示，屈服函數為

式中：f為屈服函數；M為臨界狀態(tài)線的斜率。

圖2 修正劍橋模型的屈服面Fig.2 Yield surfaces of modified Cam-clay model

圖中飽和土的屈服應力pc0是塑性體變的函數。非飽和土的屈服應力pc是塑性體變、飽和度、基質吸力的函數。

彈性變形和塑性變形分別為

式中：λ、κ為初始壓縮曲線的斜率和回彈斜率；υ為土體比容積。

采用相關流動法則

式中：G為土體的剪切模量。硬化準則為

從前面的推導得出，吸力和飽和度對屈服應力的硬化作用如圖3所示。圖中實線為脫濕時基質吸力對前期屈服應力的影響，虛線為吸濕時基質吸力對前期屈服應力的影響，圖上的數據點為Wheeler[7]所做的試驗。從圖上可以看出，當基質吸力相同時，由于存在滯回效應，脫濕對應的飽和度比吸濕時高，所以相應的屈服應力較吸濕時小，而且含水率分布較不均勻，這點與前面提出的飽和度的作用相對應。

圖3 模型得出的LC曲線與試驗數據的對比Fig.3 Comparison between simulated LC curves and tested data

從上面的分析可以看出，盡管沒有將基質吸力當作獨立變量來考慮，但推導得出的LC曲線與試驗數據擬合結果良好。

2.3 變形對土-水特征關系的影響

Wei等[37]嚴格以內變量熱動力學和多孔介質理論為基礎，提出了一個簡單實用的毛細滯回內變量模型。該模型基于多相、多孔介質理論，將毛細循環(huán)現(xiàn)象視為一個非飽和土的能量耗散過程。在循環(huán)邊界曲線給定的情況下，該模型只需一個參數就能確定任意含水率變化路徑下非飽和土的土-水特征關系。該模型改進后可以考慮體變對土水特征關系的影響。

該模型以Feng等[38]提出的土-水特征曲線為邊界線

式中：b、d、b0、d0、α為材料參數，與土體的土-水狀態(tài)相關，在脫濕和吸濕時取不同值。

當不考慮變形的影響時，基質吸力與飽和度變化關系為

通過試驗發(fā)現(xiàn)，土體的孔隙比對土-水特征曲線有重要影響[39]，孔隙比減小時，進氣值增大，土-水特征曲線向右移動。

考慮變形對土-水特征曲線的影響，有

式中：Vw、Vv分別為水和孔隙的體積；等號右邊第1項表示基質吸力對飽和度的影響，第2項描述變形對飽和度變化的影響。忽略彈性變形的影響，式（17）可以變?yōu)?/p>

式中：n為孔隙率；第2項表示毛細滯回與變形的耦合效應。

3 模型驗證

3.1 各向同性壓縮和吸濕試驗

采用孫德安[40]的各向同性加壓濕陷試驗來驗證所提出的模型。試樣采用屬于粉質黏土的珍珠黏土（pearl clay），非飽和擊實樣。此黏土的液限為49%，塑限為22%，土粒相對密度為2.71，初始含水率為26%。圖4為模型模擬的結果與試驗數據對比圖。應力路徑為從 A點（ pnet=20 kPa，Sc=147 kPa）保持吸力為常數加載到 B點（pnet=98 kPa， Sc= 147 kPa），然后保持凈應力不變，進行吸濕試驗，到達 C點（ pnet=98 kPa,Sc=0 kPa），最后土樣在 0吸力下各向同性壓縮至196 kPa。其中pnet為平均凈應力， pnet=p-ua。圖上從A點到B點和C點到D點為各向同性壓縮，B點到C點為吸濕路徑，當吸力減小時，屈服應力降低。類似地，另外一條路徑從 pnet= 196 kPa 開始吸濕。從圖上可以看出，盡管吸濕時有效應力降低，但屈服應力也在減小，土體的塑性變形大于回彈體變，因此，發(fā)生了濕陷現(xiàn)象。

圖4 凈應力與比容積的關系Fig.4 Relation between net stress and specific volume

圖5 基質吸力與濕陷變形的關系和SWCC曲線Fig.5 Relation between matric suction and deformation and SWCC

圖6為當 Sc= 147 kPa時各向同性壓縮之后，在時進行吸濕。圖中是飽和度隨平均凈應力的變化關系，因為有塑性應變發(fā)生，飽和度隨著凈應力增大而略微增大。這是由于塑性變形使孔隙比減小造成的。

圖6 飽和度隨凈應力的變化關系Fig.6 Variation of saturation with net stress

圖7為土樣在給定的平均凈應力下經過吸濕和脫濕之后的比容積變化，從圖可以看出，土體在吸濕階段會產生彈性膨脹，而在脫濕階段則表現(xiàn)出明顯塑性收縮現(xiàn)象，這與Sharma[3]得出的結論相一致。

圖7 模型對吸濕、脫濕循環(huán)中比容積變化的模擬Fig.7 Specific volume change obtained by simulation under of wetting and drying cycles

3.2 常凈應力和常吸力下的剪切試驗

采用孫德安[1]所做的非飽和土三軸試驗結果進行模擬。根據試驗結果，得出pearl clay的力學參數為 λ=0.12，κ=0.03， pc0= 70 kPa，r=4.5，m= 3.65，=0.15，G= 8 MPa。

土 樣 在 pnet= 196 kPa 、 Sc= 147 kPa時 ，q= 0 → 296 kPa加載至破壞。圖8～11為試驗結果和模型模擬的對比。圖8為體變εv與主應變ε1之間的關系。圖9、11為主應力比 σ1/σ3與主應變ε1、ε3之間的關系。圖 10為 σ1/σ3與飽和度Sr之間的關系，圖上初始的直線段表示土樣為彈性狀態(tài)，與前面忽略彈性體變對飽和度的影響假設相一致。從圖9、11中可以看出，模型得出第一主應變較小，第三主應變較大，與試驗數據有偏差，但總體上的趨勢是相一致的。圖 10中的直線段為初始的彈性變形，對飽和度沒有影響。

圖8 第一主應變與體變的試驗結果與模型預測Fig.8 Test and prediction results of maximum principal strain and volume strain

圖9 第一主應變與主應力比的試驗結果與模型預測Fig.9 Test and prediction results of maximum principal strain and ratio of principal stresses

圖10 飽和度與主應力比的試驗結果與模型預測Fig.10 Test and prediction results of saturation and ratio of principal stresses

圖11 第三主應變與主應力比的試驗結果與模型預測Fig.11 Test and prediction results of minor principal strain and ratio of principal stresses

3.3 體變對土-水特征曲線的影響

土體發(fā)生體變時，孔隙比降低，進氣值增大，因此，使得土-水特征曲線向外移動。模型所得結果如圖12所示，從圖可以看出，體變對吸濕的影響較脫濕大。這是因為體積壓縮導致孔隙變小，大孔隙比例減小，使得大孔隙對小孔隙的阻隔作用減弱，形成的滯回圈變小。

圖12 體變對土-水特征曲線的影響Fig.12 Effect of volume strain on soil-water characteristic curves

4 結 論

本文在修正劍橋模型的基礎上，考慮了非飽和土的硬化作用建立了毛細滯回與變形耦合的本構模型。該模型摒棄了傳統(tǒng)模型的LC曲線，在修正劍橋模型和毛細滯回內變量模型的基礎上增加了兩個參數考慮基質吸力和飽和度對屈服應力的硬化作用。兩個參數可通過試驗得出的LC曲線得到。其他參數可通過常規(guī)的力學和水力學試驗得到。通過與試驗結果對比，證明此模型可以模擬基本的非飽和土力學與水力學行為，包括濕陷、強度隨基質吸力的增加而增大等。而且同時考慮了體變對土-水特征曲線的影響，毛細滯回及其含水率變化歷史對變形的影響。

[1]SUN D A, SHENG D, SLOAN S W. Elastoplastic modelling of hydraulic and stress-strain behaviour of unsaturated soils[J]. Mechanics of Materials, 2007,39(3): 212－221.

[2]WHEELER S J. Inclusion of specific water volume within an elastoplastic model for unsaturated soils[J].Canadian Geotechnical Journal, 1996, 33(1): 42－57.

[3]SHARMA R S. Mechanical behavior of an unsaturated highly expansive clays[D]. Oxford, UK: University of Oxford. 1998.

[4]CHARLES W W N, PANG Y W. Influence of stress state on soil-water characteristics and slope stability[J].Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(2): 157－166.

[5]ALONSO E E, GENS A, JOSA A. A constitutive model for partially saturated soils[J]. Geotechnique, 1990, 40(3):405－430.

[6]BOLZON G, SCHREFLER B A, ZIENKIEWICZ O C.Elastoplastic soil constitutive laws generalized to partially saturated states[J]. Géotechnique, 1996, 46(2): 279－289.

[7]WHEELER S J, SIVAKUMAR V. An elasto-plastic critical state framework for unsaturated soil[J].Geotechnique, 1995, 45(1): 35－53.

[8]WHEELER S J, SHARMA R S, BUISSON M S R.Coupling of hydraulic hysteresis and stress-strain behaviour in unsaturated soils[J]. Geotechnique, 2003,53(1): 41－54.

[9]NUTH M, LALOUI L. Advances in modelling hysteretic water retention curve in deformable soils[J]. Computers and Geotechnics, 2008, 35(6): 835－844.

[10]SUN D A, LI X. An elastoplastic hydro-mechanical model for unsaturated soils[C]//Proceedings of the 7th International Conference on Computational Science, Part III: ICCS. Beijing: Springer-Verlag, 2007: 1138－1145.

[11]繆林昌. 非飽和土的本構模型研究[J]. 巖土力學, 2007,28(5): 855－860.MIAO Lin-chang. Research of constitutive model of unsaturated soils[J]. Rock and Soil Mechanics, 2007,28(5): 855－860.

[12]SHENG D, GENS A, FREDLUND D G, et al.Unsaturated soils: From constitutive modelling to numerical algorithms[J]. Computers and Geotechnics,2008, 35(6): 810－824.

[13]CUI Y J, DELAGE P. Yielding and plastic behaviour of an unsaturated compacted silt[J]. Geotechnique, 1996,46(2): 291－311.

[14]韋昌富, 李幻, 王吉利. 考慮彈塑性變形和毛細循環(huán)滯回的非飽和土本構模型[C]//第一屆全國巖土本構理論研討會論文集. 北京: 北京航空航天大學, 2008: 259－266.

[15]KHALILI N, GEISER F, BLIGHT G E. Effective stress in unsaturated soils: Review with new evidence[J].International Journal of Geomechanics, 2004, 4(2):115－126.

[16]TAMAGNINI R. An extended Cam-clay model for unsaturated soils with hydraulic hysteresis[J].Geotechnique, 2004, 54(3): 223－228.

[17]JOMMI C. Remarks on the constitutive modelling of unsaturated soils[C]//Experimental Evidence and Theoretical Approaches in Unsaturated Soils. Rotterdam:A. A. Balkma, 2000: 139－153.

[18]LI X S. Thermodynamics-based constitutive framework for unsaturated soils 2: A basic triaxial model[J].Geotechnique, 2007, 57(5): 423－435.

[19]LORET B, KHALILI N. A three-phase model for unsaturated soils[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2000, 24(11): 893－927.

[20]WHEELER S J. An alternative framework for unsaturated soil behavior[J]. Geotechnique, 1995, 41(2): 257－260.

[21]MUALEM Y. A conceptual model of hysteresis[J].Water Resources Research, 1984, 10(3): 514－520.

[22]PARLANGE J Y. Capillary hysteresis and the relationship between drying and wetting curves[J]. Water Resources Research, 1976, 12(2): 224－228.

[23]LI X S. Modelling of hysteresis response for arbitrary wetting/drying paths[J]. Computers and Geotechnics,2005, 32(2): 133－137.

[24]MIAO L C, JING F, HOUSTON S L. Soil-water characteristic curve of remolded expansive soils[J].Geotechnical Special Publication, 2006, 147(1): 997－1004.

[25]方祥位, 陳正漢, 申春妮, 等. 剪切對非飽和土土-水特征曲線影響的探討[J]. 巖土力學, 2004, 25(9): 1451－1454.FANG Xiang-wei, CHEN Zheng-han, SHEN Chun-ni.et al. A study of effect of shear on soil-water characteristic curve of an unsaturated soil[J]. Rock and Soil Mechanics, 2004, 25(9): 1451－1454.

[26]龔壁衛(wèi), 吳宏偉, 王斌. 應力狀態(tài)對膨脹土SWCC的影響研究[J]. 巖土力學, 2004, 25(12): 1915－1918.GONG Bi-wei, WU Hong-wei, WANG Bin. Influence of stress states on soil-water characteristics of expansive soils[J]. Rock and Soil Mechanics, 2004, 25(12): 1915－1918.

[27]GALLIPOLI D, WHEELER S J, KARSTUNEN M.Modelling the variation of degree of saturation in a deformable unsaturated soil[J]. Géotechnique, 2003,53(1): 105－112.

[28]GENUCHTEN M T V. A closed-form equation for predicting the hydraulic conductivity of unsatutated soils[J]. Soil Science Society of America Journal, 1980,44(5): 892－898.

[29]BISHOP A W. The principle of effective stress[J].Teknisk Ukeblad, 1959, 106(39): 113－143.

[30]HOULSBY G T. The work input to an unsaturated granular material[J]. Geotechnique, 1997, 47(1): 193－196.

[31]CADORET T, MARION D, ZINSZNER B. Influence of frequency and fluid distribution on elastic wave velocities in partially saturated limestones[J]. Journal of Geophysical Research, 1995, 100(B6): 9789－9803.

[32]CADORET T, MAVKO G, ZINSZNER B. Fluid distribution effect on sonic attenuation in partially saturated limestones[J]. Geophysics, 1998, 63(1): 154－160.

[33]KNIGHT R, NOLEN-HOEKSEMA R. A laboratory study of the dependence of elastic wave velocities on pore scale fluid distribution[J]. Geophysical Research Letters,1990, 17(10): 1529－1532.

[34]RAVALEC L M, GUEGUEN Y, CHELIDZE T. Elastic wave velocities in partially saturated rocks: Saturation hysteresis[J]. Journal of Geophysical Research, 1996,101(B1): 837－844.

[35]SUN D A, MATSUOKA H, CUI H B, et al.Three-dimensional elasto-plastic model for unsaturated compacted soils with different initial densities[J].International Journal for Analytical and Numerical Methods in Geomechanics, 2003, 27(12): 2079－2098.

[36]ROSCOE K H, BURLAND J B. On the generalized stress-srain behavior of wet clay[C]//Engineering Plasticity. Cambridge: Cambridge University Press, 1968:535－609.

[37]WEI C F, DEWOOLKAR M M. Formulation of capillary hysteresis with internal state variables[J]. Water Resources Research, 2006, 42(7): W07405. 1－16.

[38]FENG M, FRELUND D G. Hysteresis influence associated with thermal conductivity sensor measurements[C]//The 52nd Canadian Geotechnical Conference and Unsaturated Soil Group, Proceedings from Theory to the Practice of Unsaturated Soil Mechanics, Regina, Saskatchewan, Canada: [s. n.], 1999:651－657.

[39]VANAPALLI S K, FREDLUND D G, PUFAHL D E.The influence of soil structure and stress history on the soil-water characteristics of a compacted till[J].Geotechnique, 1999, 49(2): 143－158.

[40]SUN D A, SHENG D, XU Y. Collapse behavior of unsaturated compacted soil with different initial densities[J]. Canadian Geotechnical Journal, 2007,44(6): 673－686.