氣液兩相流原理在邊坡穩(wěn)定性分析中的有效性研究

董順德,孫芝玲,何玉鵬,馬文禮

(1.青海省環(huán)境地質(zhì)勘查局，青海西寧 810001；2. 青海省環(huán)境地質(zhì)重點實驗室，青海西寧 810001；3.青海省水利水電科學研究院有限公司，青海西寧 810001)

近年來，極限平衡等分析方法已被廣泛應用于雨水滲透導致邊坡崩塌的機理研究中[1-4]。降雨使得土壤飽和度增加，致使坡體質(zhì)量增加，從而導致邊坡滑動力增大，出現(xiàn)邊坡失穩(wěn)現(xiàn)象[5]。另一方面，由于坡體中水壓上升和基質(zhì)吸力減小，導致邊坡抗剪強度下降，促使滑坡發(fā)生。因此對雨水滲流行為和孔隙水壓變化的研究是評估邊坡穩(wěn)定性的必要條件。雨水一般通過地表的非飽和區(qū)進入邊坡中，而非飽和土的滲透系數(shù)與其體積含水率或基質(zhì)吸力之間的關系是非線性的，應用單純的數(shù)值解析法對其研究較為困難。有學者利用有限差分法(FLAC3D)將飽和土流-固耦合計算原理擴展到非飽和土中[6]。還有學者采用有限元法進行了飽和-非飽和滲流分析[7-8]。有限元法主要是對滲透到坡體的雨水流量進行估算，是計算降雨期間孔隙水壓力隨時間變化的較為有效的方法[9]。

在對降雨邊坡穩(wěn)定性研究中，缺乏對孔隙氣壓造成影響的分析。針對降雨對非飽和土影響的問題，有學者研究孔隙氣壓對邊坡造成破壞的影響?；A實驗證實，在土中發(fā)生空氣堵塞時，可通過氣壓的增加導致地表出現(xiàn)局部破壞的現(xiàn)象觀察氣壓變動與地表破壞的關系[10-11]；CHO[12]研究表明，雨水的滲透和水位的上升會引起氣壓上升和壓裂現(xiàn)象，堤防內(nèi)部的壓縮空氣會對堤防整體造成一定程度的損傷；SAITO等[13]針對孔隙氣壓對雨水滲透的影響進行了室內(nèi)試驗研究，分別對兩相流和單相流分析法進行了研究，分析此類方法對雨水滲透過程的適用性。室內(nèi)試驗的結(jié)果證實了通過氣液兩相流分析可以重現(xiàn)孔隙氣壓的上升和排水量增加的現(xiàn)象，但其室內(nèi)模型試驗忽略了對飽和度、基質(zhì)吸力和滲透系數(shù)之間關系的研究。

筆者對氣液兩相流理論在降雨邊坡穩(wěn)定性分析中的應用進行了研究。首先，利用圓柱形容器進行一維滲流模型實驗，采用不考慮孔隙氣壓影響的單相流分析和考慮孔隙氣壓影響的氣液兩相流分析，分別對一維滲流模型實驗進行了分析。對比結(jié)果證實了兩相流分析法的適用性與合理性。針對均勻坡度的風化巖粉質(zhì)黏土邊坡進行了人工降雨實驗，通過氣液兩相流理論對邊坡降雨滲透行為進行計算分析。采用氣液兩相流法進行應力-滲透耦合計算，模擬2012年青海省黃南州同仁縣隆務鎮(zhèn)西山發(fā)生的滑坡，分析降雨誘發(fā)邊坡失穩(wěn)的原因。

1 應力-滲透流耦合分析概述

1.1 有限差分法的特征

有限差分法和有限元法是典型的將目標視為連續(xù)體的數(shù)值分析方法。由于有限差分法能利用結(jié)構(gòu)網(wǎng)格的拓撲優(yōu)勢輕松擴大模板，構(gòu)造出高精度格式，在基于顯式的應力-變形的時間增量計算中，它比有限元法更具優(yōu)勢。本研究使用了有限差分法解析軟件FLAC，此軟件具有分析多孔介質(zhì)滲透的功能，可以與應力-變形分開，從而進行獨立的滲流分析，也可以進行應力-滲流耦合分析。

1.2 單相流分析的基本方程式

進行單相流分析時，不考慮孔隙氣壓的影響，僅考慮水的滲透條件，符合達西定律：

(1)

其中：qi為速度矢量；kij為飽和滲透系數(shù)；P為水壓；gk為重力加速度；ρf為水的密度；κ為比滲透系數(shù)(飽和滲透系數(shù)與飽和滲透系數(shù)的比),用飽和度s表示：

κ=s2(3-2s)。

(2)

流體的平衡方程如式(3)所示：

(3)

其中:qv為水的流入和流出量大小，而分析中，水被視為不可壓縮狀態(tài)物質(zhì)。

1.3 兩相流分析的基本方程

兩相流是指氣相和液相流體的流動，2種流體流動都遵循達西定律，氣相和液相相關的達西定律分別如式(4)和式(5)所示：

(4)

(5)

其中：μ為黏性系數(shù)；P為孔隙壓力；w表示液體；a表示氣體。式(6)和式(7)中的比滲透系數(shù)κr可以用Van Genuchten模型來表示[14]。

(6)

(7)

其中：a，b，c均為常數(shù)；Se為有效飽和度。在本研究中，主要采用基于Van Genuchten模型的水特性曲線：

(8)

流體的平衡方程可表示為

(9)

(10)

式中：Sw為飽和度；ζ為單位體積內(nèi)流體的體積變化；qv為流入和流出量。流體的本構(gòu)方程如下：

(11)

(12)

其中:Kw和Kg為體積彈性系數(shù)。此外在本研究中，不考慮空氣在水中的溶解作用。

2 一維滲流模型實驗的再現(xiàn)分析

2.1 實驗方法

一維滲流模型的實驗裝置如圖1所示。

圖1 一維滲流模型實驗裝置Fig.1 Experiment device of 1D seepage model

試樣是直徑為280 mm，高為500 mm的圓柱體，孔隙壓力表和張力計分別測量其中心和底部的正、負孔隙壓力。采用含水比為10%左右的粉質(zhì)黏土，土體顆粒密度為2.74 g/cm3，孔隙率為0.4。通過壓實達到預定的濕密度制備模型試樣。為了得到規(guī)定的濕密度，對每層分30次壓實，分成8層構(gòu)建模型。模型完成后，整個模型的初始飽和度將達到預計的60%左右，之后進行變水頭實驗，得出的滲透系數(shù)k=1.8×10-5m/s。

調(diào)整降雨強度，并從試樣頂部供給。為了便于水從頂部均勻地供給，在試樣頂部設置了6 cm的沙層。在砂層與試樣的邊界處設置了4個直徑為2 mm的排水孔，使無法滲入粉質(zhì)黏土的水從砂層中流出，防止積水對試樣造成過大壓力。

2.2 數(shù)值解析

一維滲流解析數(shù)值模型如圖2所示。

圖2 一維滲流解析數(shù)值模型Fig.2 Analytical numerical model of 1D seepage

為了再現(xiàn)一維滲流模型實驗，將模型尺寸設定為寬280 mm，高500 mm，滲透系數(shù)設定為k=1.8×10-5m/s。中部(No.1)和下部(No.2)的設計與模型實驗中的孔隙水壓測量點相同，并計算中部和下部的孔隙水壓。設定降雨強度為50 mm/h，在底面為不排水條件下進行了一維滲流模型實驗和數(shù)值解析。在氣液兩相流分析中，需要確定土壤水分參數(shù)，即須確定式(6)-式(8)中的a，b，c和P0值。設定b=1/2，c=1/3[15]，對于a和P0的取值，分別用張力計測量在恒定干密度壓實條件下不同飽和度試樣的基質(zhì)吸力，確定飽和度與基質(zhì)吸力的關系。隨后，通過改變式(10)中a和P0的數(shù)值，與實驗結(jié)果進行擬合，擬合結(jié)果顯示a=0.47，P0=8.5 kPa。

2.3 模型實驗和數(shù)值模擬結(jié)果

圖3分別為模型實驗、單相流分析和兩相流分析中孔隙壓力的變化。從實驗結(jié)果看，模型中心No.1測點的孔隙壓力首先開始增大，其次底部測量點No.2測點的孔隙壓力增大。這表明，滲透是從頂部向底部進行的。張力計測得的孔隙水負壓會很快歸零，說明在非飽和狀態(tài)下，孔隙壓力表測量得到的是正孔隙壓力。導致這個結(jié)果的原因是處于相同深度的孔隙壓力表和張力計的受壓部分沒有設置在相同位置。開始通水后約250 min，No.2測點的孔隙水壓力升高，No.1測點的孔隙壓力在稍作延遲后也有所增加。這是因為前面的滲透面到達模型底部后，下部飽和，水位從模型底部上升。之后2個測量點在整個模型的飽和狀態(tài)下，靜水壓力趨于一致。

圖3 模型實驗和數(shù)值模擬結(jié)果Fig.3 Model experiment and numerical simulation results

從單相流分析結(jié)果可以看出，孔隙水壓先從No.2測點開始增大，出現(xiàn)了孔隙水壓力即說明土體飽和。因此，認為在單相流分析中不能考慮孔隙空氣的運動，水以一定的速度平穩(wěn)滲透，從解析結(jié)果可以看出試樣下部比試樣中部更早地達到了飽和狀態(tài)。在兩相流分析中，No.1測點和No.2測點的孔隙壓力在大致相同時間內(nèi)由負變正，但由于有微小的誤差，No.2測點測得的結(jié)果較快。因此，認為土體飽和是從試樣的下部向上部進行的。與單相流分析相比，No.1測點和No.2測點的孔隙壓力變?yōu)檎龎核钑r間相同，其原因是在兩相流分析中，試樣在飽和層從底部上升前，中間土層的飽和度較高，使其飽和所需的水量變少。

在模型實驗中，孔隙水壓呈階梯式上升，并趨于靜水壓力值。在這2種情況下，孔隙水壓繼續(xù)緩慢增加，直至靜水壓力值。盡管實驗中使用的試樣是由不均勻性的土壤顆粒組成的，而從數(shù)值模擬角度認為模型是一個均質(zhì)連續(xù)體，但是，在兩相流分析中孔隙水壓的上升趨勢和實驗結(jié)果相同，只是No.1測點的孔隙水壓比No.2測點先開始上升。對于單相流分析的結(jié)果，試樣下部先達到飽和的時間要早于實驗結(jié)果和兩相流分析結(jié)果，無法再現(xiàn)水的滲透被孔隙空氣阻擋的情況，并且也很難捕捉非飽和狀態(tài)下的孔隙壓力變化情況。

由以上分析可知，在50 mm/h降雨強度下，降雨強度相對于土壤滲透性較大時，兩相流分析可以考慮孔隙空氣對水滲透的影響。反之，當降雨強度相對于土壤滲透率較小時，認為孔隙空氣的影響較小，可以通過單相流分析來評價滲流行為。由于單相流分析無法捕捉到非飽和狀態(tài)下的孔隙壓力變化情況，而通過兩相流分析可以得到非飽和狀態(tài)下孔隙壓力和土中空氣的影響，所以兩相流更適合評估降雨滲透導致的邊坡穩(wěn)定性。

3 邊坡滲流行為的室內(nèi)斜面模型

3.1 斜面模型實驗

斜面模型實驗如圖4所示，斜面模型的形狀和儀器的位置如圖5所示。實驗材料使用的是從同仁縣隆務鎮(zhèn)西山Ⅲ號滑坡現(xiàn)場采取的表層粉質(zhì)黏土。邊坡模型制作時試樣的初始含水量為10%，干密度為1.68 g/cm3，黏聚力為12.7 kPa，內(nèi)摩擦角為36.9°。坡面模型分8塊壓實，為了達到預定的濕密度，每段都分為3層并用木槌進行壓實。為了確認坡面和坡體重的滲透行為，在No.1—No.3測點設置了張力計，No.2和No.3測點分別安裝了孔隙壓力計和土壤水分儀，并從模型頂部調(diào)整噴水量模擬降雨。為驗證不同降雨強度對雨水的滲透行為的影響，進行了3種不同模擬降雨強度條件下的實驗，即在45°傾角下降雨強度分別為20，50和100 mm/h。

圖4 斜面模型實驗Fig.4 Slope model experiment

3.2 邊坡模型實驗再現(xiàn)數(shù)值模型

筆者根據(jù)圖5所示的模型進行數(shù)值模擬。在相同實驗條件下，在同一位置測量孔隙水壓。在數(shù)值模型中，通過對比斜面模型試驗中測得的孔隙水壓和吸力值與數(shù)值模擬得到的孔隙壓力和吸力值，對模型進行了驗證。應用兩相流分析需要的土壤水分參數(shù)b=1/2，c=1/3。在模擬中，應用了表示非飽和時的水分特征曲線的式(8)。P0是土壤的特征毛細管壓力，它隨土壤類型的不同而變化，a和P0的取值是通過下述實驗確定的。式(13)是表現(xiàn)Van Genuchten水分特征曲線的模型[14]。式中的ψ表示基質(zhì)吸力，Se表示有效飽和度，

(13)

為了確定式(9)中的a和P0的取值，在與3.1固結(jié)條件相同情形下，將風化泥巖壓實到體積為1 000 cm3的模具中，安裝張力計和土壤濕度計，從模具上方供水，測量土壤含水量和基質(zhì)吸力，直至內(nèi)部飽和。結(jié)果如圖6所示，4組的結(jié)果幾乎相同，以圖6中Ec5_1&吸力1的值為代表，與SWRC Fit中的土壤水分特征曲線的非線性回歸程序擬合，得出如圖7所示的表現(xiàn)水分特征曲線模型。從圖7得出式(13)中的參數(shù)α=0.125，n=1.49。反轉(zhuǎn)橫軸和縱軸，將Van Genuchten模型式(8)轉(zhuǎn)化為如圖8所示的曲線，得出a=0.33，P0=8.0 kPa。

圖6 體積含水率和基質(zhì)吸力關系Fig.6 Relationship between volumetric water content and matrix suction

圖7 式(13)水分特征曲線Fig.7 Water characteristic curve of Eq.(13)

圖8 飽和度和基質(zhì)吸力關系Fig.8 Relationship between saturation and matrix suction

3.3 實驗結(jié)果和數(shù)值模擬結(jié)果

圖9表示了實驗和數(shù)值模擬過程中孔隙水壓的變化情況。實驗開始時孔隙水壓大小不同的原因是由于模型制作時壓實過程產(chǎn)生的振動，導致模型中的孔隙水向下移動。通過調(diào)整模型中的初始飽和度，使分析實驗在降雨開始后立即獲得與實驗中相同的孔隙水壓值。在降雨開始后，因為No.2測點比No.1測點更靠近坡面，No.2測點的孔隙水壓立即開始增加，隨后No.1測點的孔隙水壓才有上升趨勢。當浸潤面到達計測點時，孔隙水壓迅速增加，在飽和度維持一段時間后，模型下方上升的水面到另一測點，孔隙水壓都是先增加一次，然后在短時間內(nèi)又開始增加，呈現(xiàn)階梯變化。數(shù)值分析表明，孔隙水壓并不像實驗中那樣呈階梯式增加，而是連續(xù)增加。這可能是因為和第2節(jié)的一次元滲透流模型實驗相同，無法充分再現(xiàn)現(xiàn)階段土壤孔隙中的水滲透行為。但當浸潤面已經(jīng)達到計測點時，孔隙水壓突然上升的行為與實驗中的是一樣的。其次，孔隙水壓開始增加后，即浸潤面到達計測點的順序和時間與實驗結(jié)果差別也不大。此外，降雨強度為100 mm/h的孔隙水壓結(jié)果如圖9 c)所示，與其他2種情況不同，計測點No.2測點的孔隙水壓也同樣比No.3測點的孔隙水壓增加得快。從這些結(jié)果看，雖然無法準確再現(xiàn)實際邊坡滲透的復雜行為，但總體上可以再現(xiàn)滲流行為。

圖9 孔隙水壓變化Fig.9 Changes of pore pressures

4 基于暴雨引起的滑坡數(shù)值模擬

4.1 邊坡崩塌現(xiàn)場的描述

筆者通過氣液兩相流分析對2012年青海省黃南州同仁縣隆務鎮(zhèn)西山發(fā)生的滑坡進行了模擬。滑坡全貌如圖10所示。

圖10 黃南州同仁縣隆務鎮(zhèn)西山滑坡全貌Fig.10 General view of Xishan Landslide in Longwu Town, Tongren County, Huangnan prefecture

以西山Ⅲ號滑坡為研究對象，滑坡巖性以強風化泥巖堆積的粉質(zhì)黏土為主，土質(zhì)不均勻，泥巖風化較嚴重，局部含少量礫石。強風化泥巖節(jié)理裂隙發(fā)育，風化程度較高，且形成的粉質(zhì)黏土不但自身強度低，而且對含水量變化極為敏感，遇水后易軟化，抗剪強度大大降低，易形成軟弱滑動面(帶)，為易滑地層?；麦w平均坡度約35°，由淺層滑坡形成，滑體軸向長為280～310 m，橫向?qū)挒?0～180 m，平均厚度為21.5 m。建議采用抗滑樁和樁間擋土墻的組合形式進行邊坡防治。

4.2 基于兩相流分析的滲透行為評價

圖11表示滑坡的數(shù)值模型。由于滑坡是強風化泥巖土層的表層塌陷，因此設定作為基層土的弱風化巖是彈性體，而表層的強風化泥巖土是按照Mohr-Coulomb斷裂標準的彈塑性體。在邊界條件方面，模型的右側(cè)和底部設計為不排水，而左側(cè)和坡面設計為允許排水。底面為垂直和水平方向約束，兩側(cè)為水平方向約束。弱風化巖土層和強風化巖土層的初始飽和度分別設定為1.0和0.3，降雨強度設定為100 mm/h。分析開始時，對模型進行自重模擬，將計算后的穩(wěn)定狀態(tài)定義為初始應力狀態(tài)。

圖11 滑坡再現(xiàn)數(shù)值模型Fig.11 Numerical model of landslide reproduction

模擬結(jié)果顯示：降雨開始2 h后，坡面附近飽和度增加，可以確認到坡腳孔隙水壓開始轉(zhuǎn)為正值。降雨開始3 h后，邊坡上部孔隙水壓轉(zhuǎn)為正值，強風化泥巖土層孔隙氣壓增大。降雨開始4 h后，上部的孔隙壓力繼續(xù)增大，孔隙水壓增加的位置孔隙氣壓也在增加。降雨開始6 h后，孔隙水壓力的增加領域的面積逐漸擴大，在剪切應變增大的地方，逐漸失穩(wěn)，上部坡面淺部最大剪應變增大，可以確定為表層滑坡。結(jié)果顯示，塌方發(fā)生在坡面，該處塌方最大深度約3 m，現(xiàn)場滑坡發(fā)生的位置與數(shù)值模擬確定的劃坡位置高度一致。

數(shù)值模擬結(jié)果表明，滑坡形態(tài)和滑坡位置與實際滑坡高度相似。邊坡滑動的機理是由于雨水的滲透，孔隙水壓力和孔隙壓力增大，而有效應力減小。所以表層滑坡不一定是由地下水位上升導致滑坡面有效應力下降引起的，在非飽和條件下也可能發(fā)生剪切應變，即在降雨量達到足以抬高地下水之前就有發(fā)生滑坡的可能性。

5 結(jié) 論

為驗證在孔隙氣壓影響下坡面兩相流模擬的有效性，筆者開展了一維滲流實驗和斜面模型滲流實驗，并與數(shù)值解析結(jié)果進行了對比。結(jié)論如下：

1)一維滲流模型實驗結(jié)果表明，單相流未考慮孔隙氣壓的影響，可能導致水的滲透速度比實驗結(jié)果快。在氣液兩相流分析中，實驗中各測點的孔隙水壓上升趨勢相似；

2)模型實驗結(jié)果顯示，一旦滲透面到達計測點，孔隙水壓就會突然上升。此外，數(shù)值模擬中滲透面到達計測點的時間也與實驗無明顯差異；

3)采用兩相流應力-滲流耦合分析邊坡失穩(wěn)，結(jié)果表明，滑坡發(fā)生位置和形態(tài)與實際滑坡高度相似，滑動面附近的孔隙水壓和孔隙氣壓增加是造成滑坡的重要原因之一。

本研究對氣液兩相流模擬方法在邊坡穩(wěn)定性分析中的應用進行了基礎研究?，F(xiàn)階段還沒有充分再現(xiàn)土壤孔隙中的雨水滲流行為。在數(shù)值分析和模擬實驗中均使用了均質(zhì)材料，但由于實際的自然邊坡多為非均質(zhì)體，因此，邊坡土層的非均質(zhì)性研究是今后的重要課題之一。