基于多場耦合的膨脹土邊坡非飽和降雨入滲分析

韋秉旭，陳亮勝，肖羅明，鄭 威，白玉祥

(1.長沙理工大學 交通運輸工程學院，長沙 410114； 2.湖南路橋建設集團有限責任公司，長沙 410004)

1 研究背景

降雨誘發(fā)膨脹土邊坡產(chǎn)生滑坡、崩塌等地質災害，嚴重威脅公路工程建設和運營安全，造成了大量的人員傷亡、財產(chǎn)損失和環(huán)境破壞[1]。降雨入滲引起土坡濕度場的變化是一個典型的非飽和滲流過程，隨著降雨的持續(xù)，大氣影響深度內(nèi)坡體的孔隙水壓力、含水率和重度顯著提高，導致基質吸力消散和抗剪強度降低，逐漸形成暫態(tài)飽和區(qū)[2]。此外，降雨入滲會引起膨脹土膨脹變形和強度軟化等不良性質的出現(xiàn)，加劇邊坡地質災害的形成[3]。因此，掌握降雨入滲對膨脹土邊坡非飽和滲流影響的時間和空間分布規(guī)律，能為膨脹土邊坡防治提供理論基礎與技術支持。

數(shù)值模擬因其可復現(xiàn)性、求解速度快、數(shù)據(jù)易采集、節(jié)約成本等優(yōu)點，成為廣大學者研究膨脹土邊坡非飽和滲流的重要手段。文獻[4-6]采用滲流與應力的非耦合數(shù)值方法，直接從孔隙水壓力、含水率、滲透力和暫態(tài)飽和區(qū)的時空分布規(guī)律來闡述膨脹土邊坡內(nèi)部非飽和滲流的過程，并說明了暫態(tài)飽和區(qū)與滲流發(fā)展速度、面積、時間等諸多因素有關，但目前尚未建立一個描述不同降雨歷時期間暫態(tài)飽和區(qū)變化規(guī)律的定量指標，且此方法忽略了土體變形對土骨架和孔隙水的相互作用，導致其滲透特性、水分運移通道和滲流過程與實際存在差異。有限差分軟件FLAC3D具有強大的非線性應力應變求解功能，能有效處理非飽和土滲流場與應力場相互耦合的難題。蔣中明等[7]基于內(nèi)嵌FISH語言的二次開發(fā)平臺，修正并完善了非飽和滲流的降雨入滲計算功能。汪華斌等[8]、謝強等[9]著重處理入滲邊界計算方法與非飽和區(qū)流量的兩大關鍵問題，充分驗證了有限差分軟件FLAC3D實現(xiàn)非飽和土流-固耦合數(shù)值模擬的可靠性。實際上，膨脹土的非飽和降雨入滲過程伴隨著膨脹應變和應力應變的產(chǎn)生，是典型的膨脹-滲流-應力耦合問題，如何將膨脹土的膨脹應變引入非飽和土的流-固耦合模型成為膨脹土滲流耦合分析的另一難點。

鑒于此，基于有限差分滲流理論、降雨入滲模塊的二次開發(fā)和膨脹應變的施加，采用C++和FISH語言編制程序，提出一種綜合考慮應力應變和膨脹應變的非飽和膨脹土多場耦合分析法。以某膨脹土路塹邊坡為研究對象，運用FLAC3D實現(xiàn)了降雨入滲條件下膨脹土邊坡非飽和滲流的數(shù)值模擬，并探討了不同降雨時間下膨脹應變對飽和度和暫態(tài)面積比變化規(guī)律的影響。

2 非飽和膨脹土多場耦合分析法

2.1 有限差分滲流理論

準靜態(tài)比奧理論是有限差分FLAC3D計算非飽和滲流的方法，它將多孔介質結構的飽和度、孔隙壓力和體積應變建立聯(lián)系[10]，假設非飽和滲流中不考慮氣相的影響，其流體連續(xù)性方程的表達式為

(1)

式中：s為飽和度；uw為孔隙水壓力；εv為體積應變；M為比奧模量；ζ為孔隙水運移過程中發(fā)生的流體體積改變量；n為孔隙率；t為時間；α為比奧系數(shù)。

當膨脹土邊坡經(jīng)降雨入滲處于飽和-非飽和滲流階段時，因含水率的變化脹縮現(xiàn)象明顯，總應變更新為應力應變和膨脹應變之和[11]，即

(2)

將式(2)代入連續(xù)性方程式(1)中，得到考慮膨脹應變的膨脹土連續(xù)性方程為

在降雨入滲條件下，非飽和土多孔介質水分運移的過程遵循達西定律和質量守恒定律，其達西流動方程的表達式為

qi=-Kilk(s)(uw-ρfxjgj)il。

(4)

式中：qi為單位流量向量;ρf為流體密度;Kil為滲透系數(shù)張量；Kil=ks/γw，ks為達西定律中單元的飽和滲透系數(shù)，γw為流體重度;k(s)為相對滲透系數(shù)，滿足k(s)=s2(3-2s)∈[0，1]；xj為笛卡爾坐標分量;gj為重力加速度分量。

流體質量平衡方程的表達式為

?ζ/?t=-qi,i+qv。

(5)

式中qv為流體源強度。

非飽和區(qū)的飽和度和滲透系數(shù)是隨基質吸力動態(tài)變化的參數(shù)，其值分別通過土-水特征曲線和滲透系數(shù)曲線與孔隙水壓力建立聯(lián)系，從而賦值給每一節(jié)點或單元?？赏ㄟ^Fredlund-Xing模型[12]確定基質吸力與飽和度的曲線公式，即

(6)

式中：ψ為基質吸力；ψr為殘余基質吸力；a為進氣值的倒數(shù)；m為控制殘余含水率的參數(shù)，m=3.67ln(1/si)，si為進氣值所對應飽和度；n為曲線斜率。

非飽和區(qū)單元的滲透系數(shù)k為單元飽和滲透系數(shù)與相對滲透系數(shù)的乘積，其表達式為

k=kss2(3-2s) 。

(7)

2.2 降雨入滲模塊二次開發(fā)

基于FLAC3D平臺的降雨入滲模塊二次開發(fā)主要涉及初始孔壓場的生成、非飽和區(qū)飽和度和滲透系數(shù)的設置、邊界條件的開發(fā)三大技術問題。

(1)在FLAC3D中設置流體抗拉強度，能允許負孔隙水壓力的存在，這為實現(xiàn)非飽和滲流數(shù)值模擬提供了可能。通過內(nèi)置變量zone.pp(z)可指定節(jié)點孔隙水壓力，軟件內(nèi)部運用自動插值的方法獲取土體的孔隙水壓力，即可生成與實際工程適應的初始孔壓場。

(2)通過式(6)利用內(nèi)置變量gp.sat(z)獲取非飽和區(qū)的節(jié)點飽和度，并采用距離加權的方法轉化為單元飽和度。土體的滲透系數(shù)為單元屬性參數(shù)，根據(jù)相對滲透系數(shù)公式，采用內(nèi)置變量z.fluid.prop (z)計算土體單元的相對滲透系數(shù)，并通過式(7)獲取非飽和區(qū)單元修正后的滲透系數(shù)。

(3)FLAC3D計算滲流邊界時，通常按降雨強度或飽和滲透系數(shù)的大小作為流量邊界施加于坡面，事實上，降雨入滲及出滲邊界是一個動態(tài)變化的邊界[13]。采用接觸面單元提供的intface功能獲取坡面不同位置的滲透系數(shù)，調用降雨強度與坡面滲透系數(shù)兩者之間的小值作為坡面入滲率，并實時監(jiān)控坡面節(jié)點孔隙水壓力的大小，如當節(jié)點孔隙水壓力>0時，節(jié)點附近出現(xiàn)積水，將該節(jié)點的流量邊界調整為壓力邊界。

根據(jù)上述三大技術問題的解決措施編寫FISH程序，采用FISHCALL命令調用程序，實現(xiàn)降雨入滲對土質邊坡非飽和滲流的模擬。

2.3 膨脹應變的施加

通過開發(fā)的降雨入滲模塊與修正有效應力后的應力場模塊相互耦合，可構成普通土的非飽和流-固耦合模型，但膨脹土非飽和滲流的變化與普通土的區(qū)別主要在于膨脹應變的產(chǎn)生[14]。從式(3)膨脹土的流體連續(xù)性方程可知，膨脹土是典型的膨脹變形場-非飽和滲流場-應力場多場耦合問題。目前，在FLAC3D內(nèi)嵌本構模型中，暫無本構模型適用于膨脹土脹縮性這一特殊性質。因此，需修正內(nèi)嵌本構模型來體現(xiàn)膨脹應變對膨脹土邊坡非飽和滲流的影響。

膨脹應變在各向同性彈性體中沿土體的各個方向相同，因此，膨脹應變的主應變方向相同，切應變?yōu)?，其膨脹應變分量為

(8)

根據(jù)式(2)可知：應力應變和膨脹應變遵循Hooke定律，其應變分量為

(9)

同時將式(9)變換為總應力分量形式，即

(10)

式中：σij和εij分別為總應力分量和總應變分量；σkk和εkk分別為三個方向的主應力之和和主應變之和；E為彈性模量；v為泊松比；δij為Kornecker符號，當i=j時，δij=1，否則取0。

將式(10)改寫成有限差分的總應變增量格式，即

(12)

膨脹應變的施加主要在Visual Studio 2010平臺中操作，其開發(fā)工作基于內(nèi)嵌本構模型的C++代碼，利用上述理論修正頭文件(后綴為.h)和源文件(后綴為.cpp)，并根據(jù)文獻[15]修正屈服準則和流動法則，由此生成一個動態(tài)庫文件(后綴為.dll)。其中，頭文件的修正包括私有變量和成員函數(shù)的重新定義，源文件的修正通過實際函數(shù)代替基類的成員函數(shù)來實現(xiàn)，在成員函數(shù)中，主要對Initialize()和Run()函數(shù)進行修正，其作用分別為運行初始化模型的中間變量與根據(jù)應變增量計算應力增量。

在FLAC3D計算非飽和滲流的流-固耦合分析中，調用修正后的膨脹土本構模型，即可實現(xiàn)非飽和膨脹土的膨脹-滲流-應力多場耦合分析。其原理為：基于有效應力原理建立應力場與非飽和滲流場的聯(lián)系，可以通過降雨入滲引起非飽和滲流的時空分布控制膨脹應變來影響應力應變關系的變化過程，也可以通過體現(xiàn)應力應變關系的膨脹應變來影響非飽和降雨入滲的變化過程，三者實現(xiàn)相互耦合，簡化計算流程如圖1所示。

圖1 計算流程Fig.1 Computation process

3 工程算例

根據(jù)相關工程資料，考慮南陽某膨脹土路塹邊坡為計算模型，經(jīng)長期氣候作用，距坡面深度3 m范圍內(nèi)分布著棕黃色膨脹土(上層土)，以下為褐黃色膨脹土(下層土)。圖2所示為路塹邊坡的網(wǎng)格模型，一級坡和二級坡的邊坡坡率為別為1∶1.5和 1∶1.75，中間設臺階連接，寬度為2 m，設計Ⅰ-Ⅰ、Ⅱ-Ⅱ兩個截面作為計算剖面。

圖2 網(wǎng)格模型Fig.2 Geometric model

表1 數(shù)值模擬參數(shù)Table 1 Numerical simulation parameters

計算參數(shù)：膨脹土邊坡的土-水特征曲線和飽和滲透系數(shù)分別采用壓力板和雙環(huán)滲透試驗，并分別采用式(6)和式(7)對數(shù)據(jù)進行擬合，其余計算參數(shù)根據(jù)室內(nèi)外試驗確定，參數(shù)具體取值見表1。實際工程中假設膨脹土的干密度一定，膨脹應變主要與初始含水率和含水率的變化情況有關。因此，膨脹應變隨含水率的變化關系采用無荷膨脹試驗確定[16]，其曲線關系如圖3所示，并通過MatLab確定其多元回歸的擬合公式為

式中：εw為膨脹應變；w0和w分別為初始含水率和過程含水率。

圖3 膨脹應變隨含水率的變化Fig.3 Variation of expansion strain with water content

邊界及水文條件：模型底部固定約束，頂部自由，四周水平約束。在降雨入滲數(shù)值計算過程中，路面設為不透水邊界，坡頂水平面和坡面設置為降雨入滲及出滲邊界。值得說明的是，當降雨強度小于坡面入滲能力時，形成無壓入滲階段，此階段為流量邊界；當降雨強度大于坡面入滲能力時，形成有壓入滲階段，此階段為壓力邊界。在降雨結束后，根據(jù)坡面入滲及出滲的情況重新調整，當水分由于重力作用向坡內(nèi)入滲時，坡面孔隙水壓力轉變?yōu)樨撝担藭r調整為零流量邊界；當水分向坡外溢出時，坡面孔隙水壓力轉變?yōu)檎?，此時調整為零壓力邊界。

根據(jù)當?shù)貧庀筚Y料擬定降雨強度為1.65×10-6m/s，降雨持續(xù)3 d，停雨持續(xù)5 d，其初始孔壓場如圖4所示，并采用暫態(tài)分析法每隔12 h記錄一次計算結果，實時監(jiān)測邊坡模型的變化情況。

圖4 初始孔壓場Fig.4 Initial pore water pressure

4 計算結果及分析

4.1 孔隙水壓力分布

利用本文提出的多場耦合分析法模擬非飽和膨脹土的降雨入滲過程，得到如圖5所示不同降雨歷時階段的孔隙水壓力變化。從圖5可以看出：降雨過程中，孔隙水壓力的變化與降雨時間正相關，坡面最先出現(xiàn)暫態(tài)水壓力，隨降雨時間推移，坡內(nèi)的孔隙水壓力也逐漸上升；降雨結束后，水分因重力作用會繼續(xù)向坡內(nèi)運移，濕潤鋒逐漸被均化，其上緣的孔隙水壓力減小，下緣的孔隙水壓力增加，且隨深度增加，土體固結度提高，濕潤鋒向坡內(nèi)擴展速率變慢，介于-90～0 kPa之間孔隙水壓力的區(qū)域增加。還可以看出：由于下層土的滲透能力小于上層土，濕潤鋒難以透過滲透率低的下層土，導致降雨入滲對非飽和滲流場的影響區(qū)域集中在上層土，對下層土孔隙水壓力的影響程度低，在土層交界處形成一條顯著的自上而下孔隙水壓力逐級遞減帶，具有距坡面位置越遠，受降雨入滲影響程度越小的分布特征。其數(shù)值模擬結果滿足文獻[17]和文獻[18]的非飽和膨脹土降雨入滲分析，說明本文提出的多場耦合分析法可進一步研究膨脹應變對膨脹土非飽和降雨入滲的變化規(guī)律。

圖5 孔隙水壓力分布Fig.5 Distribution of pore water pressure

4.2 膨脹應變對飽和度的影響規(guī)律

根據(jù)邊坡初始飽和度的狀況，其區(qū)域可劃分為暫態(tài)飽和區(qū)、非飽和區(qū)、初始飽和區(qū)。實際降雨入滲過程中，土體孔隙難以全部被水分填充，故將飽和度>0.9的區(qū)域劃分為飽和區(qū)[19]。

圖6為邊坡截面Ⅰ-Ⅰ和截面Ⅱ-Ⅱ處的飽和度隨邊坡深度變化的曲線關系。由圖6可知:①降雨3 d時，暫態(tài)飽和區(qū)的范圍成懸掛型向坡內(nèi)擴展，降雨入滲深度和飽和區(qū)深度分別處于2.5～3.0 m、1.8～2.3 m之間，且不考慮膨脹變形的降雨入滲深度和暫態(tài)飽和區(qū)的擴展深度要大于考慮膨脹應變的情況；②停雨5 d時，濕潤鋒上緣變?yōu)榉秋柡蛥^(qū)，降雨入滲深度介于3.5～4.0 m之間，土層交界處存在較大范圍的飽和區(qū)，暫態(tài)飽和區(qū)的面積具有顯著的時間滯后性；③截面Ⅰ-Ⅰ和Ⅱ-Ⅱ飽和度變化的差異性主要表現(xiàn)在土層交界處，且考慮膨脹應變的飽和度變化幅度小于不考慮膨脹應變的情況。

圖6 飽和度變化Fig.6 Variation of saturation

造成上述現(xiàn)象的原因是：降雨入滲改變了土體含水率的時空分布，引起含水率發(fā)生變化的中下部土體膨脹應變變化量顯著增加，產(chǎn)生濕潤鋒的影響范圍和水分運移的擴散速度受控于膨脹應變，當考慮膨脹應變時，膨脹土在降雨增濕過程中產(chǎn)生膨脹變形，因吸水膨脹土體被壓密，引起土顆粒間的孔隙尺寸減小，水分經(jīng)滲流通道的連通孔隙變窄，從而造成膨脹土的滲透性能降低，水分擴散速率變緩，抑制了降雨入滲對水分擴散的過程。在降雨結束后，坡面附近的飽和度降低，膨脹土開始出現(xiàn)失水收縮的現(xiàn)象，引起膨脹應變衰減，但其飽和度在降雨入滲過程中始終大于初始飽和度，膨脹應變繼續(xù)發(fā)揮著抑制降雨入滲對水分擴散的作用。

4.3 膨脹應變對暫態(tài)面積比的影響規(guī)律

為定量分析降雨入滲過程中膨脹土邊坡暫態(tài)飽和區(qū)面積的時程變化，定義暫態(tài)面積比η(暫態(tài)飽和區(qū)面積與初始非飽和區(qū)面積之比)來揭示膨脹土邊坡非飽和降雨入滲的變化規(guī)律，用以彌補描述暫態(tài)飽和區(qū)時程效應的不足。暫態(tài)面積比的表達式為

式中：ST(m2)和S(m2)分別為暫態(tài)飽和區(qū)面積和初始非飽和區(qū)面積；ST(pt)和S(pt)分別為暫態(tài)飽和區(qū)像素面積和初始非飽和區(qū)像素面積。運用圖像處理技術，根據(jù)兩者像素之比可直接獲取暫態(tài)面積比。

圖7所示為考慮膨脹應變與不考慮膨脹應變的暫態(tài)面積比時程變化，可以看出，隨降雨歷時的增加，暫態(tài)面積比先上升后下降，呈“凸”型分布，且暫態(tài)面積比的“凸”點發(fā)生在停雨前期，其極大值為21.5%。這是由于膨脹土邊坡在降雨期間和停雨前期，濕潤鋒呈飽和狀態(tài)向坡內(nèi)擴展，懸掛型暫態(tài)飽和區(qū)面積持續(xù)上升；隨停雨時間的持續(xù)，濕潤鋒抵達土層交界處，暫態(tài)飽和區(qū)的面積向下擴展受阻，而坡面水分不斷向坡外溢出或向坡內(nèi)慢滲，濕潤鋒被均化，暫態(tài)面積比呈不斷縮小的趨勢，且暫態(tài)面積比的降低速率明顯小于增長速率。在整個降雨歷時期間，針對考慮膨脹應變和不考慮膨脹應變的情況，暫態(tài)面積比的“凸”點形成時間分別在停雨2 d和1 d；降雨結束5 d后，暫態(tài)面積比分別變?yōu)?3.1%和10.2%，相差2.9%。其原因是膨脹應變對非飽和降雨入滲過程起抑制作用，水分擴散作用變緩，考慮膨脹應變的暫態(tài)面積比的變化速率變慢，其暫態(tài)面積比的雨后時間滯后效應長于不考慮膨脹應變時的情況。

圖7 暫態(tài)面積比的時程變化Fig.7 Change in transient area ratio over time

膨脹土邊坡非飽和降雨入滲的數(shù)值結果驗證了本文提出的多場耦合分析法的正確性和合理性。在非飽和降雨入滲過程中，一方面，因降雨增濕產(chǎn)生的膨脹變形降低了土體的滲透率，抑制了降雨入滲對水分的擴散作用；另一方面，膨脹土的膨脹應變受非飽和降雨入滲過程的控制，導致膨脹應變的發(fā)展過程不同步，造成膨脹應變變化量的大小差異幅度也不盡相同。因此，采用數(shù)值模擬的方法應綜合考慮膨脹土的脹縮性對降雨入滲過程的影響，從而進一步分析多場耦合作用下的非飽和膨脹土邊坡的破壞機制與變形機理。

5 結 論

(1)根據(jù)有限差分滲流理論和增量理論，將膨脹應變引入非飽和土的流-固耦合模型，采用FISH語言和C++語言編制程序，進行了非飽和降雨入滲模塊和膨脹土本構模型的二次開發(fā)，提出了一種適用于非飽和膨脹土的多場耦合分析法。

(2)濕潤鋒影響范圍集中在邊坡淺層，濕潤鋒隨深度增加向坡內(nèi)擴展速率變慢，雨后介于-90～0 kPa間的孔隙水壓力區(qū)域增加，且土層交界處形成一條顯著的自上而下孔隙水壓力逐級遞減帶。

(3)降雨增濕產(chǎn)生的膨脹變形造成土體滲透性能降低，抑制了降雨入滲對水分的擴散作用，導致膨脹土邊坡的降雨入滲深度、濕潤鋒擴展和均化速度以及飽和區(qū)時空分布受膨脹應變的控制。

(4)暫態(tài)面積比的時程效應呈“凸”型分布，其“凸”點發(fā)生在停雨階段，且在考慮膨脹應變后，暫態(tài)面積比的變化速率會變慢，且雨后時間滯后效應相應延長。