非飽和土一維滲流與變形耦合問題的數值分析

趙思同， 高國耀， 周鳳璽

(1.甘肅省交通規劃勘察設計院股份有限公司， 甘肅 蘭州 730050；2.蘭州理工大學 土木工程學院， 甘肅 蘭州 730050)

與飽和土不同，非飽和土中不僅有固相和液相，而且還有一定量的氣相存在，它廣泛存在于干旱、半干旱地區和地下水位以上的土體中。對非飽和土滲流和變形問題的研究在解決工程實踐中經常遇到的地基沉降變形、滲透破壞等有著十分重要的意義，如地基滲流與變形，降雨引起的邊坡失穩、垃圾填埋的處理等[1-2]。

由于土體中的應力場和孔隙流體壓力分布影響著土骨架的變形，與此同時，土體孔隙特征的改變將會引起流體滲透相關參數的變化，進而影響土體中滲流場的分布。可見，非飽和土的應力場、滲流場和位移場是相互影響，相互作用的。因此，很有必要采用滲流-變形(流-固)耦合的方法對非飽和土的力學行為進行研究。

國內外學者關于非飽和土的滲流與變形問題研究已有許多的理論成果。徐炎兵等[3]基于多孔介質力學原理，建立能模擬非飽和土兩相流動與變形耦合問題的理論模型，利用Galerkin法對控制方程進行離散，得到控制方程的有限元計算格式，并驗證了所提出的分析方法在模擬非飽和土滲流以及變形問題時的有效性，從而為定量研究飽和-非飽和滲流以及變形問題提供了一條有效途徑；吳禮舟等[4]通過Fourier積分變化，考慮了任意初始的孔隙水壓力分布和流量邊界條件，研究表明變形與滲流耦合對非飽和土中水的入滲產生重要的影響，進一步分析了耦合解的各個參數的影響；唐凱等[5]運用非飽和土滲流理論對邊坡進行滲流場和應力場耦合分析；Olivella等[6]、Thomas等[7]分別從水-氣-熱三場耦合的角度分析了多孔介質孔隙中各個組分的滲流與變形問題；Schrefler[8]在Zienkiewicz等[9]的理論模型基礎上，分析了非飽和土孔隙介質中濕相與非濕相兩相滲流與變形問題；秦冰等[10]基于混合物理論，推導了在氣壓力變化條件下水與水蒸氣相變平衡所滿足的限制方程；胡冉[11]開展了非飽和土水力全耦合本構模型及數值模擬方法等內容的研究。然而由于非飽和土中流-固耦合問題的復雜性，到目前為止，無論是數學模型的描述還是數值求解方法的有效性方面尚不完善。

基于多孔介質理論，考慮固-液-氣三相耦合的情況，在耦合過程中考慮了孔隙中溶液、溶解在水中的氣相、氣相中的水分、干燥氣體等組分，建立了非飽和土多孔介質中固相、液相、氣相三相耦合的數學模型。考慮一維非穩態問題，利用Laplace積分變換，選取孔隙氣壓力、孔隙水壓力和土顆粒位移以及它們的一階導數作為狀態變量，得到了問題的狀態方程組。在給定的邊界條件下，結合打靶法和數值Laplace逆變換求解了該耦合的非線性變系數微分方程組，并通過與已有的試驗結果相比較，驗證了模型的有效性。該研究結果將對非飽和土滲流與變形耦合問題理論與數值分析研究提供一定的參考價值。

1 基本方程

分別以s，l，g表示非飽和土中的固相骨架、孔隙液相和孔隙氣相。考慮土體的孔隙率為n，孔隙被液相和氣相共同填充，體積分數分別為φl和φg，則有n=φl+φg。孔隙液相由液態水和溶解的氣體兩種組分組成，氣相由水蒸氣和干燥氣體組成。非飽和土各相以及各組分構成如圖1和表1所示。

圖1 非飽和多孔介質的示意圖

表1 非飽和多孔介質中的各相組成

1.1 物理方程

非飽和土有效應力公式為：

(1)

φl=nSl，φg=nSg，Sg+Sl=1

(2)

基于彈性小變形假設，利用幾何方程可得到土體的物理模型為：

(3)

式中：G和K分別為土體的剪切模量和體積模量；ui為土骨架位移分量。

1.2 質量守恒方程

組分α(α=w,a分別代表水分、干燥氣體)在π相(π=l,g分別表示液相和氣相)中的質量守恒方程可表示為[12-14]：

(4)

1.2.1 水分質量守恒

由于水分主要存在于液相和氣相中，忽略固體顆粒中的水分含量，結合式(1)和式(4)，可得:

(5)

(6)

(7)

(8)

基質吸力Pm與飽和度的關系可表示為:

Pm=Pg-Pl=PdSe-1/?

(9)

考慮土體骨架位移速率的散度與體積應變率的關系[15]:

(10)

由固相質量守恒方程可得:

(11)

引入平均相對速度[16]:

(12)

(13)

根據廣義達西定律，液相和氣相的對流通量分別表示為[15]:

(14)

qg=-Kg(Pg+ρgg)

(15)

式中的各參數表達式如下[15]：

1.2.2 干燥氣體質量守恒

假設氣體僅存在于氣相中，不考慮固相中的干燥氣體。結合式(5)和式(4)可得:

(16)

(17)

1.3 動量守恒方程

多孔介質的動量守恒方程為：

·σ+f=0

(18)

式中:σ為總應力；f=g[(1-n)ρs+nSlρl+nSgρg]為體積力。

根據有效應力原理式(1)和物理方程式(3)，結合方程式(18)可得到三相耦合的動量守恒方程：

(19)

2 一維非穩態問題的控制方程

通過對基本方程式(8)—式(19)進行理論推導，最終得到的以孔隙水壓力Pl、孔隙氣壓力Pg和骨架位移u為基本未知數的非飽和滲流-變形耦合問題的控制方程：

(20)

ρaSg?tu,i-ρaKg(Pg,ii+ρgg)-

(21)

g[(1-n)ρs+nSlρl+nSgρg]=0

(22)

考慮一維問題，并應用Laplace積分變換即可得到頻域內問題的控制方程為：

(23)

(24)

g[(1-n)ρs+nSlρl+nSgρg]=0

(25)

(26)

方程式(25)是一個多物理場耦合的變系數偏微分方程組，一般很難獲得其解析解答，而打靶法是一種求解強非線性微分方程的有效方法[17-18]。求解時使用Runge-Kutta和Newton-Raphson方法不斷調整初值參數，使得初值參數的解適用于終點邊界條件，從而得到原邊值問題的頻域上的精確解。為了獲得時間域上的結果，基于矩方法的穩定化算法，通過Hausdoff矩方法將其進行數值逆Laplace變換[19]。

3 數值算例

表2 非飽和土的邊界條件

其中，孔隙水壓力和孔隙氣壓力的值為相對大氣壓力。

3.1 Liakopoulos試驗的數值模擬

根據Liakopoulos[20]等提供的土樣參數，選取如表3所示的非飽和土物理力學參數進行了數值模擬。圖2給出了在重力作用下土體中水壓力模擬值與Liakopoulos的試驗值的對比結果，可以看出兩者變化基本一致。

圖2 水壓力試驗結果與模擬結果對比

表3 Del Monte砂土的材料物理力學參數

在與Liakopoulos試驗相同的邊界條件下，圖3—圖5分別給出了土體不同位置處孔隙水壓力，孔隙氣壓力和位移隨時間的變化情況。由圖3可知土體底部的孔隙水壓力隨著時間的增大而逐漸減小，大概在120 min后，孔隙水壓力的變化逐漸趨于緩和，這是因為隨著時間的增長，最后土柱在底部排水作用下孔隙水壓力逐漸趨于穩定。由于氣體可以自由進入土柱頂部邊界，使得土柱試樣下部的氣壓變化率明顯慢于試樣上部的砂樣。從圖4可以看出，隨著時間的不斷增大，孔隙氣壓力的變化率不斷減小，最后孔隙氣壓力逐漸趨于緩和，這是因為隨著時間的增大，土柱試樣頂部的氣體不斷進入土柱內，使得氣壓與大氣壓的值逐漸接近。由圖5可以看出土柱試樣在不同高度處的豎向位移隨著時間的變化逐漸增大，但是大概在120 min后，豎向變形逐漸趨于緩和，這是因為隨著時間的變化，土柱試樣在重力作用下豎向固結逐漸趨于穩定。

圖3 不同位置處孔隙水壓力的變化

圖4 不同位置處孔隙氣壓力的變化

圖5 不同位置處豎向位移的變化

對比Liakopoulos的試驗數據結果可以得知，本文提出的滲流與變形耦合分析的數值模型在理論數值分析上基本反映了非飽和土體內液相與氣相的存在狀態，同時也較為準確地反映了氣-水兩相流動變化與土骨架位移變形之間的相互影響關系。

3.2 參數分析

為了進一步研究非飽和土中水-氣-土骨架顆粒三相耦合的過程，對不同孔隙率、孔隙水壓力、滲透系數等條件下的孔隙氣壓力的變化規律進行了參數分析。

為了分析土體孔隙率對非飽和土體中流-固耦合的影響，圖6—圖8分別給出了土體上下邊界的孔隙氣壓力均為0 kPa，土體上邊界孔隙水壓力為0 kPa，下邊界孔隙水壓力分別為-20 kPa，土體孔隙率分別為0.3、0.4、0.5條件下，土體不同位置處孔隙氣壓力、孔隙水壓力以及飽和度隨時間的變化情況。從圖中可以看出土體中的孔隙氣壓力的分布趨勢基本保持不變，但是隨著土體中孔隙率的增加，土體中的孔隙氣壓力逐漸降低。隨著土體孔隙率的增大，土體中孔隙水壓力隨時間而逐漸增大，當土體的孔隙率較大時土體孔隙中的水分遷移越快，即土體底部在相同時間下土體水分飽和度越大。

圖6 不同孔隙率下孔隙氣壓力的變化

圖7 不同孔隙率下孔隙水壓力的變化

圖8 不同孔隙率下飽和度的變化

當土體孔隙率為n=0.3，上下邊界的孔隙氣壓力均為0 kPa，土體上邊界孔隙水壓力為0 kPa時，圖9給出了下邊界的孔隙水壓力分別為-8 kPa、-9 kPa、-10 kPa的條件下，土體中孔隙氣壓力的分布情況。從圖9中可以明顯地看出，隨著下邊界孔隙水壓力的減小即土體飽和度的減小，土體中孔隙氣壓力的分布趨勢基本保持不變，即孔隙氣壓力在土體內部會增大使其大于土體邊界處的孔隙氣壓力。隨著土體孔隙中水分飽和度的減小，相同時間條件下土體孔隙中的孔隙氣壓力隨之減小。

圖9 不同孔隙水壓力下孔隙氣壓力的變化

當土體孔隙率為n=0.3，上下邊界的孔隙氣壓力均為0 kPa，土體上邊界孔隙水壓力為0 kPa，土體下邊界的孔隙水壓力為25 kPa時，圖10、圖11分別給出了土體的滲透系數分別為8.6×10-6m/s、5.6×10-6m/s、2.6×10-6m/s的條件下，土體中孔隙水壓力和飽和度的變化情況。由圖10、圖11可以明顯地看出，隨著土體滲透系數的減小，在相同時間條件下土體中的孔隙水壓力而逐漸減小，即隨著土體滲透系數的減小土體中水分的遷移能力而逐漸減弱，即在相同時間條件下土體底部水分的飽和度隨著滲透系數的增大而隨之逐漸增大。

圖10 不同滲透系數下孔隙水壓力的變化

圖11 不同滲透系數下飽和度的變化

4 結 論

基于多孔介質理論，建立了等溫條件下非飽和土多孔介質中水分、氣體、骨架位移等多場耦合的數學模型。在求解非線性耦合方程時選取孔隙水壓力、孔隙氣壓力、土骨架位移、時間以及它們的一階導數變化作為狀態變量，采用打靶法進行數值分析。通過與已有的試驗結果相比較，驗證了模型的有效性。結果表明：

(1) 非飽和土體中水、氣兩相的滲流變化是相互影響的，并且與土體的變形相互作用。

(2) 土體孔隙率的變化對相關物理量的變化有著顯著影響。隨著孔隙率的增大使得土體氣相和液相滲透系數的增加，從而使土體中孔隙氣壓力降低。