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

趙思同， 高國耀， 周鳳璽

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

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

由于土體中的應(yīng)力場(chǎng)和孔隙流體壓力分布影響著土骨架的變形，與此同時(shí)，土體孔隙特征的改變將會(huì)引起流體滲透相關(guān)參數(shù)的變化，進(jìn)而影響土體中滲流場(chǎng)的分布?？梢姡秋柡屯恋膽?yīng)力場(chǎng)、滲流場(chǎng)和位移場(chǎng)是相互影響，相互作用的。因此，很有必要采用滲流-變形(流-固)耦合的方法對(duì)非飽和土的力學(xué)行為進(jìn)行研究。

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

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

1 基本方程

分別以s，l，g表示非飽和土中的固相骨架、孔隙液相和孔隙氣相?？紤]土體的孔隙率為n，孔隙被液相和氣相共同填充，體積分?jǐn)?shù)分別為φl和φg，則有n=φl+φg?？紫兑合嘤梢簯B(tài)水和溶解的氣體兩種組分組成，氣相由水蒸氣和干燥氣體組成。非飽和土各相以及各組分構(gòu)成如圖1和表1所示。

圖1 非飽和多孔介質(zhì)的示意圖

表1 非飽和多孔介質(zhì)中的各相組成

1.1 物理方程

非飽和土有效應(yīng)力公式為：

(1)

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

(2)

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

(3)

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

1.2 質(zhì)量守恒方程

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

(4)

1.2.1 水分質(zhì)量守恒

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

(5)

(6)

(7)

(8)

基質(zhì)吸力Pm與飽和度的關(guān)系可表示為:

Pm=Pg-Pl=PdSe-1/?

(9)

考慮土體骨架位移速率的散度與體積應(yīng)變率的關(guān)系[15]:

(10)

由固相質(zhì)量守恒方程可得:

(11)

引入平均相對(duì)速度[16]:

(12)

(13)

根據(jù)廣義達(dá)西定律，液相和氣相的對(duì)流通量分別表示為[15]:

(14)

qg=-Kg(Pg+ρgg)

(15)

式中的各參數(shù)表達(dá)式如下[15]：

1.2.2 干燥氣體質(zhì)量守恒

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

(16)

(17)

1.3 動(dòng)量守恒方程

多孔介質(zhì)的動(dòng)量守恒方程為：

·σ+f=0

(18)

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

根據(jù)有效應(yīng)力原理式(1)和物理方程式(3)，結(jié)合方程式(18)可得到三相耦合的動(dòng)量守恒方程：

(19)

2 一維非穩(wěn)態(tài)問題的控制方程

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

(20)

ρa(bǔ)Sg?tu,i-ρa(bǔ)Kg(Pg,ii+ρgg)-

(21)

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

(22)

考慮一維問題，并應(yīng)用Laplace積分變換即可得到頻域內(nèi)問題的控制方程為：

(23)

(24)

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

(25)

(26)

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

3 數(shù)值算例

表2 非飽和土的邊界條件

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

3.1 Liakopoulos試驗(yàn)的數(shù)值模擬

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

圖2 水壓力試驗(yàn)結(jié)果與模擬結(jié)果對(duì)比

表3 Del Monte砂土的材料物理力學(xué)參數(shù)

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

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

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

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

對(duì)比Liakopoulos的試驗(yàn)數(shù)據(jù)結(jié)果可以得知，本文提出的滲流與變形耦合分析的數(shù)值模型在理論數(shù)值分析上基本反映了非飽和土體內(nèi)液相與氣相的存在狀態(tài)，同時(shí)也較為準(zhǔn)確地反映了氣-水兩相流動(dòng)變化與土骨架位移變形之間的相互影響關(guān)系。

3.2 參數(shù)分析

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

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

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

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

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

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

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

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

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

圖11 不同滲透系數(shù)下飽和度的變化

4 結(jié) 論

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

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

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