考慮水力耦合作用的降雨誘發土質滑坡過程數值模擬研究

曹 磊, 吳 清, 趙 廣

(1.成都市城市安全與應急管理研究院,四川 成都 610023;2.四川省公路規劃勘察設計研究院有限公司,四川 成都 610041;3.中國科學院地理科學與資源研究所,北京 100101)

0 引言

降雨誘發土質滑坡是我國西南山區最為普遍的災害之一,也是災害防治的重點[1-2]。降雨誘發土質滑坡的過程可劃分為邊坡失穩階段和滑坡運動階段,這2個階段之間存在緊密關聯。針對上述2個階段,有研究利用室內實驗[3-4]、物理模型[5-6]、數值模擬[7-8]等方法揭示其形成機理。對于土質滑坡而言,降雨入滲是坡體從非飽和狀態到飽和狀態的過程,該過程中坡體內部物理力學性質發生顯著改變,從而導致坡體失穩破壞。以抗剪強度為例[2,5],降雨入滲導致坡體內部含水率增加與基質吸力下降,使得土體抗剪強度下降,進而造成邊坡失穩。

為定量描述降雨入滲導致土體抗剪強度的改變,Lu等[9]將廣義有效應力的概念引入淺層邊坡穩定性計算中,并廣泛應用于局部飽和條件下的邊坡穩定性分析。國內外將非飽和土入滲理論引入淺層邊坡穩定性研究并取得長足進展[10-11]。當邊坡失穩后,如何定量描述滑坡運動過程是評價降雨誘發滑坡危害的重要前提。一方面,通過物理實驗可測量直接反映滑坡運動特征的滑坡厚度和滑坡速度等[12-13]參數,隨著測量技術的發展,坡體內部孔隙水壓力、有效壓應力以及地震動信號等參數也間接被用于分析滑坡運動特征[14-15];另一方面,數值計算方法成為滑坡運動過程定量評價常用方法之一[16-17],如Savage等[18]首次提出基于深度平均理論的滑坡物理模型,并成功應用于顆粒流、滑坡、泥石流等物質的運動模擬中[17-18];此外,侵蝕效應[19]、剪脹效應[20]等顯著影響滑坡運動特征的現象也被考慮。盡管針對降雨誘發土質滑坡過程單一階段的研究已較為深入,但有關降雨誘發土質滑坡過程多階段的耦合研究較為缺乏,可同時模擬邊坡失穩和滑坡運動過程的數值算法還有待建立。

因此,本文以降雨誘發土質滑坡過程為研究對象,建立可反映邊坡失穩與滑坡運動2個階段演化機理的物理模型,明確2個階段轉化銜接的關鍵因子,提出適用于2個階段物理模型同步計算的數值方法。搭建考慮水力耦合作用下從邊坡失穩到滑坡運動全過程的數值模擬流程,實現降雨誘發土質滑坡過程的定量描述,在此基礎上開展數學模型的參數敏感性分析,旨在為降雨誘發土質滑坡提供一種考慮水力耦合作用的可定量化評價方法,為此類滑坡防災減災提供理論與技術支撐。

1 考慮水力耦合作用的降雨誘發滑坡耦合數學模型

1.1 邊坡失穩數學模型

水力耦合作用對土質邊坡失穩階段的影響主要體現在降雨入滲下邊坡內部土體含水率與吸應力的改變。為描述這一過程,本文采用基質勢形式的Richards模型[8],模型中涉及的土—水特征曲線和非飽和滲透系數等采用Van-Genuchten公式[21]擬合。模型方程表達如下:

(1)

(2)

(3)

式中:

t—土壤吸應力變化的時間,s;

z—土壤垂向厚度,m;

ψ—土體基質吸力,Pa;

θ—土體含水率;

θr—土體殘余含水率;

θs—土體飽和含水率;

C—比水容重,C=dθ/dh,其中,h表示邊坡厚度,m;

Se—土體飽和度,Se=(θ-θr)/(θs-θr);

α,n,m—與邊坡物理特征相關的擬合參數;

K—非飽和滲透系數;

Ks—飽和滲透系數。

考慮非飽和條件下的降雨誘發土質邊坡穩定安全系數通常使用基于無限邊坡假設的一維極限平衡模型進行評價[9],方程表達如下:

(4)

式中:

FoS—邊坡穩定安全系數;

φ—邊坡內摩擦角,°;

β—邊坡斜面傾角,°;

cp—邊坡內部有效黏聚力,Pa;

ρs—邊坡固相密度,kg/m3;

ρf—邊坡孔隙水密度,kg/m3。

1.2 滑坡運動數學模型

降雨誘發土質邊坡失穩后通常處于松散及臨近飽和狀態,變形破壞時由于剪縮作用會產生超孔隙水壓力,可降低滑坡底部有效摩擦阻力,提升滑坡流動性;相反,對于密實飽和土體,變形破壞時由于剪脹作用導致孔隙水壓力消散,可增加滑坡底部有效摩擦阻力,降低滑坡流動性。因此,為考慮水力耦合作用對滑坡運動階段的影響,本文采用可考慮剪脹效應的Iverson-George模型[20],該模型由質量守恒方程、動量守恒方程、體積含水率方程和孔隙水壓力方程組成,方程表達如下:

(5)

(6)

(7)

(8)

(9)

式中:

u,v—滑坡沿x和y方向的運動速度,m/s;

gx,gy,gz—重力加速度沿x,y和z方向的分量,m/s2;

pb—滑坡內部孔隙水壓力,Pa;

kap—滑坡內部側向土壓力系數,通?？稍O為1[17];

D—滑坡剪脹率,D=-2K(pb-ρfgzh)/μh;

μ—滑坡內部孔隙水黏度,Pa·s;

c—滑坡固相體積分數,c=1-θ;

ρ—滑坡密度,kg/m3,ρ=ρfθ+ρsc;

ω—滑坡固相壓縮率;

κ—滑坡剪脹角,°,κ=atan(c-ceq),其中,ceq表示滑坡固相臨界體積分數;

χ,ξ—無量綱系數,χ=(ρs+3ρf)/4ρ,ξ=3/(2ωh)+gzρf(ρ-ρf)/4ρ。

上述模型將剪脹理論耦合至滑坡運動模型,基于滑坡所處狀態(固相體積分數大小)來判斷其內部發生剪脹或剪縮效應,當滑坡固相體積分數大于臨界值時表明滑坡處于剪脹狀態,反之則處于剪縮狀態,從而判定孔隙水壓力消散或增加。

1.3 數學模型耦合

通過上述數學模型分析可見,含水率是同時影響邊坡失穩和滑坡運動2個階段的關鍵參數,可以作為實現階段耦合的銜接因子,如圖1。通過邊坡失穩模型可獲得邊坡內部含水率、基質吸力、穩定安全系數、失穩體積等變量狀態。當滑坡開始運動后,將邊坡失穩體積和含水率作為初始條件代入滑坡運動模型,可獲得滑坡運動過程的厚度、速度等變量狀態。值得注意的是,由于Iverson-George模型是由Navier-Stokes方程基于深度平均理論簡化而來,無法考慮滑坡垂直方向上的孔壓分布,故滑坡物質含水率可取垂直方向上的平均值代入。

圖1 考慮水力耦合作用的降雨誘發土質滑坡全過程物理模型

2 數學模型求解算法

由于上述數學模型方程的結構特征存在差異,且包含變量的演化方向也有所不同,故分別采取相應的算法進行數值求解。

針對Richards模型方程,為保持計算穩定性,采用隱式格式離散,具體離散形式如下:

(10)

式中:

k—時間迭代次數;

i—垂向計算網格位置;

ψi—非飽和區基質勢在垂向i網格處取值,Pa;

Ci—比水容重在垂向i網格處取值;

Δt—單次迭代時間步長;

Δz—垂直方向網格長度。

可見,方程(10)屬于典型的三對角矩陣形式,利用Tridiagonal Matrix Algorithm(TDMA)方法[22]求解。

針對Iverson-George模型方程,同樣采用隱式格式離散,具體離散格式如下:

(11)

式中:

j—平面計算網格位置;

Δx—水平x方向網格長度;

Δy—水平y方向網格長度;

U—變量矩陣;

F,G—通量矩陣;

S—源項矩陣。

可見,方程(11)屬于非線性雙曲格式,利用Finite Volume Method(FVM)方法求解[17]。值得注意的是,為保證2個階段模型同步計算時間的一致性,時間步長取決于2個階段中偏小的一方。綜上所述,搭建考慮水力耦合作用的降雨誘發土質滑坡過程數值流程,如圖2。

圖2 考慮水力耦合作用的降雨誘發土質滑坡過程數值計算流程

3 數學模型及數值算法驗證

為驗證所用各階段物理模型及數值算法的適用性,首先,對楊詩秀等[23]開展的室內一維降雨入滲非飽和土壤試驗進行模擬,該試驗選擇60cm高的砂壤土柱,分2種工況開展,一種將土柱水平放置,一端表面以飽和含水率作為進水條件(第一類邊界條件),經150min后測定土柱內部含水率分布;一種將土柱垂直放置,頂端以0.02427cm/min的供水強度噴灑在土柱表面(第二類邊界條件),經180 min后測定土柱內部含水率分布。土柱初始含水率θ為0.025,飽和含水率θs為0.4,殘余含水率θr為0.01,滲透系數K=1.42θ10.24為土壤飽和度函數,擬合參數取值α=0.025,n=1.391。網格大小為1cm條件下的數值結果如圖3中黑色虛線,表明模型數值數據與試驗實測數據是一致的。

圖3 2種條件下入滲監測數據與計算數據對比

其次,對2015年深圳光明新區紅坳村土質滑坡運動過程開展模擬?，F場調查顯示,該滑坡發生時處于飽和松散狀態,受剪脹效應影響具有滑動迅速、滑距長等特點[24],計算模擬參數通過相關文獻資料[24-25]獲取,見表1。網格大小為5m條件下的模擬結果,如圖4。圖4給出了滑坡運動4個時刻的厚度分布,時間為0.041667h時刻線條代表實際滑坡波及范圍,滑坡堆積邊緣處存在幾棟被波及的建筑。結果表明計算滑坡堆積范圍與實際滑坡堆積范圍吻合度較好。

表1 模型計算所需參數取值

圖4 深圳滑坡運動過程不同時刻滑坡形態分布

4 數學模型參數敏感性分析

為深入研究水力耦合作用在降雨誘發土質滑坡失穩及運動過程中的影響,利用上述耦合物理模型開展全過程模擬及關鍵參數敏感性分析。由于缺乏降雨誘發土質滑坡從邊坡失穩到滑坡運動全過程的實測資料,本文采取一維數值模擬方案進行研究。具體方案設置,如圖5。初始邊坡位于坡度39°的斜坡上,邊坡水平長5m,斜坡水平長10m,邊坡初始含水率0.1,飽和含水率0.4,殘余含水率0.05,飽和滲透系數7×10-9m/s,雨強50mm/h,滑坡固相臨界體積分數0.64,滑面摩擦系數0.65,其余模型參數與上述模擬案例一致。

圖5 降雨誘發土質滑坡初始條件示意圖

邊坡失穩過程計算結果,如圖6。隨著降雨入滲,坡體內部含水率由上到下逐漸提升,相應的坡體穩定性系數隨之下降,當降雨至5h,坡面向下1m處發生失穩破壞,坡體失穩部分趨于飽和,平均含水率為0.38。將上述計算結果代入滑坡運動模型模擬滑坡運動過程(如圖7),由于滑坡初始固相體積分數小于臨界固相體積分數,滑坡初始滑動階段處于剪縮狀態,孔隙水壓力增大,從而減少滑坡基底摩擦阻力,增加滑坡流動性,滑坡最大速度可達10m/s。當滑坡進入平面后,摩擦阻力增大,導致滑坡減速堆積,最大水平滑距為32m。

圖6 降雨入滲下邊坡含水率及穩定安全系數變化

圖7 滑坡運動過程中不同時刻滑坡形態

在上述案例的基礎上,選擇滲透系數及土體含水率作為研究水力耦合作用的關鍵因子開展敏感性分析。該參數對降雨入滲過程的影響分析已有較多研究成果[5-6,26],因此本文僅針對滑坡運動過程進行研究。在其余模型參數相同的條件下選取K=7×10-8,7×10-9,7×10-10進行模擬,在相同時間內的滑坡形態,如圖8(a)。從圖8(a)可知,滲透系數較小,水力耦合作用對滑坡運動能力影響明顯,這是由于當滲透系數較大時滑坡內部孔隙水壓力耗散迅速,導致滑坡底部摩擦阻力減小幅度變小,從而減少了滑坡運動距離。同樣,在其余模型參數相同的條件下選取不同滑坡固相體積分數c=0.61,0.62,0.64進行模擬,在相同時間內的滑坡形態,如圖8(b)。從圖8(b)可知,滑坡固相體積分數小于臨界體積分數時,滑坡運動能力越強,滑坡固相體積分數的細微差異可對滑坡運動能力帶來顯著變化。綜上可知,滑坡土體物質特征(如滲透系數、含水率等)均對滑坡運動能力存在著關鍵影響作用。

圖8 不同條件下的滑坡運動形態

5 結論

針對降雨誘發土質滑坡形成過程,將其劃分為2個階段:邊坡失穩階段和滑坡運動階段。通過分析水力耦合作用對降雨誘發土質滑坡的影響機理,采用Ricahrds模型與Iverson-George模型相耦合,以及針對方程結構特征采用不同數值方法求解,搭建降雨誘發土質滑坡過程數值模擬流程,實現該過程的定量計算及參數敏感性分析,得到以下結論。

(1)水力耦合作用是影響降雨誘發土質滑坡過程的重要因素,且對不同階段影響的方式存在差異。邊坡失穩階段,主要通過入滲影響邊坡含水率變化;滑坡運動階段,主要通過剪脹影響內部孔隙水壓力進而改變滑坡基底摩擦阻力。

(2)通過模擬物理實驗及實際案例,驗證了構建的物理模型及數值算法的適用性。數值模擬分析表明,土體滲透系數及含水率均對滑坡運動能力存在顯著影響。同一初始條件下,土體滲透系數越小,滑坡運動能力越強,土體趨于飽和時含水率越大,滑坡運動能力越強,相比于滲透系數而言土體含水率對滑坡運動影響程度較大。