土水特征曲線的一個滯回模型

陸業(yè)奇，劉斯宏，傅中志

(河海大學水工結構研究所，江蘇南京 210098)

非飽和土是多相混合體，通常認為是由土骨架、孔隙水和孔隙氣三相組成.Fredlund等[1]認為:應將收縮膜(液-氣分界面)作為獨立的第四相，其狀態(tài)對非飽和土的強度、變形等特性均具有顯著的影響;建立非飽和土平衡條件時，需要同時考慮土體的整體平衡和收縮膜的局部平衡，故建議采用凈應力和基質吸力2個應力變量來描述非飽和土的應力狀態(tài).Fredlund的雙應力變量物理意義明確，且在試驗中可以獨立控制，在20世紀90年代非飽和土本構模型研究中得到了廣泛的運用[2-3].但基于凈應力和基質吸力的本構模型在數(shù)值計算中有諸多不便:非飽和土數(shù)值計算中通常假定孔隙氣壓力為常量，故凈應力與基質吸力分別對應總應力與負孔壓，無法實現(xiàn)非飽和狀態(tài)與飽和狀態(tài)之間的連續(xù)過渡[4].正因為如此，許多學者近年來重新審視了Bishop的廣義有效應力原理[5-6]，并提出了一系列基于廣義有效應力與基質吸力的本構模型[7-8].

運用Bishop建議的有效應力必然涉及有效應力組合系數(shù) χ的選擇.原則上 χ可以作為飽和度Sw的函數(shù)，亦可取作基質吸力s的函數(shù)[9].Houlsby[10]通過能量分析指出:若將Bishop有效應力公式(1)中的 χ取為Sw，則有效應力與應變共軛，s與體積含水率nSw共軛，這為采用廣義有效應力原理提供了理論基礎.

采用Sw作為 χ還可以直接反應出土體含水量對材料強度、變形等力學特性的影響[9].此外，對土體的總應力平衡方程和孔隙水的連續(xù)性方程進行有限元離散后可以發(fā)現(xiàn)，當 χ=Sw時，單元的勁度矩陣是對稱的，這對于提高非飽和土問題的求解效率極具意義.鑒于上述原因，非飽和土計算中廣泛采用飽和度作為有效應力組合系數(shù)[11].

從數(shù)值計算的角度來說，采用飽和度作為有效應力組合系數(shù)必須建立飽和度與基質吸力之間的關系，同時進行力學和水力學分析(通常是耦合的)，從而逐步完成初邊值問題的求解.可見，除材料的力學本構模型外，還必須建立描述飽和度與基質吸力關系的水力學本構模型.土力學研究歷史中，對力學本構模型進行了大量的研究，建立了彈塑性、亞塑性、邊界面塑性等豐富的建模體系，并針對土體的不同變形特性提出了大量的本構模型.相比而言，水力學本構模型得到的重視程度遠低于力學本構模型，在處理具體問題時通常僅采用脫濕曲線或者吸濕曲線，而忽略滯回特性，在涉及干濕循環(huán)過程時，得到的解答一般不能反應土體的實際含水量變化.

近年來，部分學者吸收了力學本構模型建模過程中的基本概念和方法，建立了全面考慮吸濕曲線、脫濕曲線和掃描曲線的數(shù)學模型.Wei等[12]基于熱力學2個基本定律，認為毛細松弛和毛細滯回都是一種耗散過程，可以借助內變量來表征這一過程，并提出一個用內變量來考慮毛細滯回的模型.Li[13]基于邊界面塑性理論，定義了給定飽和度時掃描曲線與邊界曲線(主吸濕曲線和主脫濕曲線)斜率之間的關系，并通過積分建立了掃描曲線的表達式.劉艷等[14]對Wei和Li等的工作進行了分析，將不可恢復的含水率作為內變量，提出了一個修正的簡化模型，并運用試驗對該模型進行了檢驗.

本文以Li模型為例，研究了邊界面水力學本構模型的特點，并指出Li模型預測的掃描曲線易于出界的不足，提出了克服該不足的方法，并建立了一個新的水力學模型.由于以飽和度作為狀態(tài)變量的數(shù)值計算中需根據(jù)基質吸力確定飽和度，故本文模型將體積含水量作為因變量，以基質吸力作為自變量.

1 Li滯回模型簡介

Li將土體的主吸濕曲線(fd(，Sr)=0)和主脫濕曲線(fw(，Sr)=0)作為2條邊界曲線，限定所有的掃描曲線(f(s，Sr)=0)必須位于兩者之間，且當掃描曲線接近邊界曲線時，兩者斜率逐漸趨同，以保證掃描曲線與邊界曲線融合.為建立掃描曲線斜率與邊界曲線斜率之間的關系，Li將干濕轉換點 α作為記憶量，并在趨近的邊界曲線上定義 α的映射點，如圖1所示.為數(shù)學運算方便，Li對基質吸力s取對數(shù)(s*=ln s)，并在半對數(shù)圖s*～Sr中分析掃描曲線與邊界曲線斜率之間的關系，如圖2所示.

圖1 水力狀態(tài)變量定義Fig.1 Definitions of hydraulic state variables

圖2 吸濕/脫濕路徑和水力狀態(tài)變量的變化Fig.2 Update of wetting/drying paths and hydraulic state variables

現(xiàn)記掃描曲線上的飽和度增量為d Sr，相應的基質吸力增量為d s*，并記增量d s*在趨近的邊界曲線上的映射增量為d*.Li基于對試驗結果的觀察，認為始終大于等于1，且該比值取決于 s*同 α*和*之間的距離.當 s*趨近于 α*時可近似為無窮大;當s*趨近于*時，該比值趨近于1.由此Li給出了式(2)，以表征 d s*，d*的關系:

式中β是該模型唯一參數(shù)，可以通過試驗數(shù)據(jù)進行標定.

現(xiàn)將式(2)從干濕轉化點α*到當前吸力狀態(tài)s*進行如下積分:

可以得到封閉形式的掃描曲線表達式:

式中“+”用于脫濕過程，“-”用于吸濕過程.對于每一個給定的飽和度Sr，式(4)右端均為已知，由此即可計算相應的對數(shù)吸力值s*，并給出任意干濕過程中的掃描曲線.

筆者對一些試驗進行模擬時發(fā)現(xiàn)，Li模型能夠較好地重現(xiàn)掃描曲線，但亦存在不合理之處，即預測的掃描曲線可穿越邊界，如圖3所示.脫濕掃描曲線在A點之后與主脫濕曲線(MDC)的斜率近似相同，但掃描曲線穿越了主吸濕曲線(MWC)，這與引入邊界面的初衷相悖.實際上，式(2)僅保證了脫濕過程中掃描曲線不穿越MDC，僅通過斜率條件無法完全控制掃描曲線的形態(tài).在實際運用中，一旦預測的掃描曲線出界，其后續(xù)的干濕循環(huán)過程不可能再得到合理的模擬.Li還指出，該模型只能直接得到給定飽和度增量情況下吸力的變化;在非飽和土有限元計算中，通常需根據(jù)吸力變化求解飽和度增量，這可以通過迭代計算[13]確定，但這對于提高數(shù)值計算的效率是不利的.

圖3 預測掃描線出界Fig.3 Predicted scanning curve out of boundary

2 一個新的滯回模型

2.1 基本模型的建立

筆者參照Li的建模思路，通過掃描曲線體積含水率增量與邊界曲線體積含水率增量之間的關系，來建立掃描曲線表達式.假定土體在(θ0，s)點發(fā)生干濕轉換，該記憶量在趨近邊界曲線上的投影記為(，s)，如圖4所示.現(xiàn)設土體狀態(tài)由(θ，s)變化為(θ+dθ，s+d s)，相應的投影點由(，s)變化為(+d，s+d s)，則給定吸力增量d s時，必然有dθ≤d，為滿足該條件，可設

圖4 模型狀態(tài)變量定義(吸力以水頭表示)Fig.4 Definitions of state variables(suction marked with water head)

對式(5)兩端積分可得

式中θ0為干濕轉換記憶點的體積含水率.式(6)可以顯式地表達如下:

式中“-”用于脫濕過程，“+”用于吸濕過程.

這樣就建立了吸力變化驅動體積含水率變化的滯回模型，該模型與Li模型相似，僅含有一個參數(shù)β.值得指出的是式(7)建立的模型并未解決掃描曲線出界的問題，下面對Li的基本模型進行修正.

2.2 出界問題的處理

如前所述，僅通過掃描曲線與邊界曲線的斜率關系不足以對掃描曲線的范圍進行全面的控制，建模過程中宜同時使用2條邊界曲線共同確定掃描曲線.現(xiàn)將掃描曲線的體積含水率θ分成θ1和θ2:其中θ1為脫濕和吸濕過程中相應的邊界曲線對掃描曲線的貢獻，可由(7)式得到;θ2為另一條邊界曲線對掃描曲線的貢獻.在脫濕過程中可令

類似地，在吸濕過程中令

式中:θWT(s)，θDR(s)——給定吸力下，主吸濕曲線和主脫濕曲線上的體積含水率.本文采用Feng等[15]建議的函數(shù)，即

式中:b，d——材料參數(shù)，對于干燥曲線和浸潤曲線分別取不同的值;θsat——飽和體積含水率;θirr——殘余體積含水率.

為將掃描曲線嚴格地控制在2條邊界曲線之間，將掃描曲線的含水率θ表示成θ1和 θ2的線性組合，即

基本模型的數(shù)值試驗表明，出界問題主要取決于干濕轉換時的體積含水率 θ0，為此取下述形式的組合系數(shù):

式中 μ為控制掃描曲線的另一個參數(shù).容易看出，通過這樣的線性組合，掃描曲線必然會落在2條邊界曲線之間，從而避免了掃描曲線出界的現(xiàn)象.

3 模型驗證

a.用本文模型預測caribou粉土(土樣的材料參數(shù)見表1).在脫濕掃描曲線中 β=4，μ=3;吸濕掃描曲線中β=3，μ=5.從圖5可以看到，修正后的模型預測的掃描曲線不再出現(xiàn)初始階段的直線段，也沒有出界的現(xiàn)象，并且能較好地擬合試驗點.

圖5 模型預測caribou粉土與試驗點比較Fig.5 Comparison of model predictions and experimental results for caribou silt loam

b.利用本文模型預測砂土(土樣的材料參數(shù)參見表1).對于脫濕掃描曲線 β=4，μ=1.5;吸濕掃描曲線 β=4，μ=3.從圖6可以看到，模型在預測砂土時，也可以克服上述直線段和出界2個問題.除了在吸濕過程起始點飽和度較高時，模型預測曲線與試驗點存在著一些偏差，其余的試驗點都能與模型較好地吻合.

圖6 模型預測砂土與試驗點比較Fig.6 Comparison of model predictions and experimental results for sand

表1 土樣材料參數(shù)Table 1 Material parameters of soil samples

4 結 語

本文根據(jù)Li的邊界面建模思路，以基質吸力為自變量，以體積含水率為因變量，建立了一個可以考慮復雜干濕循環(huán)過程的土水特征曲線模型.本文模型中引入2個參數(shù)控制掃描曲線的形態(tài)，解決了Li模型可能出界的問題，并且能夠較好地模擬不同土類的試驗結果.

本文模型可以用于非飽和土力學、水力學耦合本構模型的建立;也可以用于非飽和土降雨入滲、心墻壩庫水位下降等實際工程問題的求解.

[1]FREDLUNDDG ，RAHARDJOH.Soil mechanics for unsaturated soils[M].New York:John Wiley&Sons，1993:66-90.

[2]ALONSO E E，GENSA，JOSA A.A constitutive model for partially saturated soils[J].Géotechnique，1990，40(3):405-430.

[3]WHEELERS J，SIVAKUMARV.An elastic-plastic critical state framework for unsaturated soil[J].Géotechnique，1995，45(1):35-53.

[4]SHENG Dai-chao，GENSA，F(xiàn)REDLUNDDG ，et al.Unsaturated soils:from constitutivemodelling to numerical algorithms[J].Computers and Geotechnics，2008 ，35:810-814.

[5]NUTHM，LALOUI L.Effective stress concept in unsaturated soils:clarification and validation of a unified framework[J].International Journal for Numerical and Analytical Methods in Geomechanics，2008，32:771-801.

[6]LALOUI L，NUTH M.On the use of the generalized effective stress in the constitutive modeling of unsaturated soils[J].Computers and Geotechnics，2009 ，36:20-23.

[7]SUN De-an，MATSUOKA H，CUI Hong-bin，et al.Three-dimensional elastic-plastic model for unsaturated compacted soils with different initial densities[J].International Journal for Numerical and Analytical Methods in Geomechanics，2003，27:1079-1098.

[8]RUSSELL R，KHALILIN.A unified bounding surface plasticity model for unsaturated soils[J].International Journal for Numerical and Analytical Methods in Geomechanics，2006，30:181-212.

[9]GENSA，SANCHEZM ，SHENG Dai-chao.On constitutivemodeling of unsaturated soils[J].Acta Geotechnical，2006，1:137-147.

[10]HOULSBY G T.Thework input to an unsaturated granular material[J].Géotechnique，1997，47(1):193-196.

[11]ABAQUSTheory Manual Version 6.8[R].West Lafayette:Simulia Corp，2008:629-701.

[12]WEI Chang-fu，DEWOOLKARM M.Formulation of capillary with internal state variables[J].Water Resour Res，2006 ，42:W07405.

[13]LI Xiang-song.Modelling of hysteresis response for arbitrary wetting/drying paths[J].Computers and Geotechnics，2005，32:133-137.

[14]劉艷，趙成剛.土水特征曲線滯后模型的研究[J].巖土工程學報，2008，30(3):399-405.(LIU Yan，ZHAO Cheng-gang.Hysteresis model for soil-water characteristic curves[J].Chinese Journal of Geotechnical Engineering，2008，30(3):399-405.(in Chinese))

[15]FENGM，F(xiàn)REDLUNDDG.Hysteretic in fluenceassociated with thermal conductivity sensor measurements[C]//Proceedings from Theory to the Practice of Unsaturated Soil Mechanics:The 52nd Canadian Geotechnical Conference and Unsaturated Soil Group.Canada，Regina:RDGINA SASD，1999，14:214-220.

[16]TOPPG C.Soil water hysteresis in silt loam and clay loam soils[J].Water Resour Res，1971 ，7(4):914-920.

[17]POULOVASSILISA.A hysteresis of porewater in granular porous bodies[J].Soil Sci，1970，109(1):5-12.