路基邊坡雨水浸透變形模擬研究

馮振華，王利超，韓玉龍

摘要：為了減少由于雨水滲透引起的路基邊坡局部變形和失穩現象，在概述模型實驗現象的基礎上，采用飽和及非飽和滲流分析方法以及彈塑性有限元分析方法，模擬路基邊坡模型中水位的變化和滑動破壞發生前的初始變形，建立漸進式邊坡失穩的預測方法，通過實驗進行再現分析，以明確邊坡初始變形機理。研究結果表明：1）在飽和及非飽和滲流模擬中，當滲流系數設置值比滲透實驗結果大時，模型實驗的再現性良好;2）在彈塑性有限元分析中，通過邊坡水位升高導致變形系數（剛度）減小，進而能夠分析出模型實驗中邊坡的初始變形。研究過程為分析降雨期間路基邊坡穩定性提供了新方法，其結果可為工程施工過程中避免發生大規模崩塌現象采取適當措施提供依據，為施工標準的設計和制定提供重要參考。

關鍵詞：巖土工程;降雨模型實驗;初始變形;飽和及非飽和滲流分析;彈塑性有限元分析

中圖分類號：TU45文獻標識碼：ADOI： 10.7535/hbgykj.2021yx05004

Simulation research on rainwater infiltration deformation of

roadbed slope

FENG Zhenhua，WANG Lichao，HAN Yulong

（Henan Transportation Planning and Design Institute Company Limited，Zhengzhou，Henan 450000，China）

Abstract：In order to reduce the local deformation and instability of subgrade slopes caused by rainwater infiltration，on the basis of an overview of the model experiment phenomena，the saturated and unsaturated seepage analysis and the elastoplastic finite element analysis were used to simulate the change of water level in the roadbed slope model and the initial deformation before sliding failure.Then the progressive slope instability prediction method was established.Finally，through the reconstruction analysis of the rules with experiments，the initial deformation mechanism of the slope was clarified.The results show that：1） in saturated and unsaturated seepage simulations，when the seepage coefficient is set larger than the result of the seepage test，the reproducibility of the model experiment is good;2） in the elastoplastic finite element analysis，the deformation coefficient （stiffness） can be reduced by increasing the water level of the slope，thereby the initial deformation of the slope in the model experiment can be analyzed.The study provides a new way for analyzing the stability of roadbed slopes during rainfall，and the results can provide a reference for taking appropriate construction measures to prevent the occurrence of large-scale collapse，which is also beneficial for formulating future construction standards.

Keywords：geotechnical engineering;rainfall model experiment;initial deformation;aturated and unsaturated seepage analysis;lastoplastic finite element analysis

河南地區每年7—9月為雨季，近年來大部分地區中雨、大雨以及暴雨日數呈增加趨勢，自然災害頻繁。極端降水會引起路基土體結構變形、坍塌，對公路路基影響巨大，嚴重影響交通運輸安全。

豫東地區土質多為粉質砂土，降雨過程中邊坡土體結構破壞的狀況主要通過模型實驗進行分析，確定邊坡的穩定性[1-3]。對降雨過程中邊坡水位變化和崩塌機理的研究發現，粉質砂土邊坡的穩定性分析大多基于完全剛塑性體的極限平衡假設。但在全尺寸模型[4]和離心機模型[5]實驗中，使用由粉質砂土模型邊坡進行降雨滑坡實驗，隨著降雨造成的邊坡中的水位和孔隙水壓力的增加，邊坡會在某時間點發生局部變形，之后破壞逐漸向邊坡上部延伸，最終導致土體大規模滑動，即漸進式邊坡失穩。此類研究中的邊坡破壞形態并非極限平衡法中所假設的瞬間滑動。在漸進式破壞發生前，坡面還出現了裂紋和局部變形。考慮到實際邊坡狀況，需要采取適當的施工方案以便在初期變形前防止土體大規模滑動崩塌發生。因此，闡明初始變形的機理，建立邊坡漸進式破壞的預測方法顯得尤為重要。

近年來SPH法等粒子分析方法在對邊坡漸進式失穩的分析中被證明是有效的，并得到了坍塌滑動過程中的影響因素分布的結論[6-8]。但在粒子分析方法中，需要的參數更多并且難以從現場實驗勘測和室內土力學實驗中獲得。因此，通過有限元法簡單地再現邊坡漸進式失穩具有更大的實際應用價值。

本研究對模型邊坡實施了人工降雨實驗，記錄了降雨過程中坡體的飽和度和孔隙壓力變化，并進一步通過飽和及非飽和滲流模擬再現了模型邊坡中的水位變化，并嘗試通過彈塑性有限元法再現粉質砂土邊坡模型中水位上升引起的邊坡初始變形來研究其整體變形機理。在彈塑性有限元分析中，主要針對巖土材料的變形系數展開研究，并假設邊坡地下水位升高時，變形系數隨有效應力的減小和應變的增加而產生折減。

1降雨邊坡模型實驗

1.1模型概要

邊坡模型的概略圖如圖1所示。

水平支撐地層厚為70 mm，在其上制作了坡度比為1∶1.5，高600 mm、寬600 mm的邊坡模型，并在土箱壁上涂抹硅脂，從而減小摩擦。模型中使用的土體材料取自河南省道323高家溝大橋附近的粉質砂土（土粒相對密度s=2.683 g/cm3，平均粒徑D50=0.165 mm，不均勻系數Uc=2.59，曲率系數U'c=1.32，細粒含量（質量分數）Fc=10.8%，最大干密度dmax=1.533 g/cm3，最優含水率w = 19.5%），制作方法為普羅克特壓實法，干密度d /最大干密度dmax×100%=95%，w=19.5%。以邊坡體壓實度Dc=75%，w=13.0%為目標，每次鋪設50 mm厚度土層，然后用振動器壓實后再進行下一層土層的鋪設。研究中，在坡箱底部放置一層平均直徑約10 mm的礫石材料，以防止由于水位快速上升發生滲透崩塌。

在邊坡模型的坡頂安裝位移計，在坡體中安裝孔隙水壓計和土壤水分儀用于測量土體含水量變化，儀器安裝布置如圖1所示。降雨量設置為40 mm/h，降雨持續450 min，其后采用灑水降雨。在支撐土層的3個位置安裝壓力計，以測量邊坡水位。

1.2數據分析

土體壓力計讀數隨降雨時間的變化如圖2所示，壓力計值表示了邊坡中地下水位高度。降雨持續450 min后直至模型邊坡崩塌過程中，隨著降雨不斷進行，壓力計值不斷增大，說明邊坡模型內土體水位持續上升;降雨551 min后，坡面土體開始發生變形。實驗過程發現，壓力計值和孔隙水壓力值增加不同步（如圖3所示），這是因為壓力計值是在測點周圍的土體高度飽和后其讀數才開始上升，而孔隙水壓力在降雨過程中土體連續滲透作用下就會產生，并不是在土體高度飽和狀態下才會發生。由于壓力計讀數與孔隙壓力上升過程的趨勢基本一致，因此可判斷壓力計讀數是合理的。模型邊坡從輕微變形到坍塌的過程如圖4所示。圖中①處，在深度約200 mm的地方坡體發生了初始裂縫和輕微變形;②處，地下水開始從模型坡體產生裂縫的地方滲出;③處，一定時間后，滲水點上方的斜坡表面開始出現裂縫;④處，以坡體頂部裂縫為邊界的滑塌持續發生。另外從圖4中可以觀察到坡面下端的塌陷土變得泥濘。

自坡面有微小初期變形開始，邊坡截面形狀隨時間的變化如圖5所示。在初期變形發生10 min后邊坡下端發生了較大的變形，并且發生了小規模的崩塌。隨著崩壞規模逐漸變大，在坡面發生變形60 min后崩塌范圍延伸到了邊坡的頂部，最終形成大規模圓弧狀滑動。

圖3中亦顯示了飽和度Sr和孔隙壓力增量Δuw隨著時間的變化，飽和度是根據土壤水分儀獲得的體積含水量以及與制作模型時的相同ρd（ρd=1.15 g/cm3）計算出。Δuw為澆灌過程中孔隙水壓的增量，其初始值為澆灌開始時的孔隙水壓值。

飽和度在開始灑水降雨的初期階段，邊坡底端（M14和M15）土體的Sr值有所上升，之后此處Sr值約為0.8，高于其他區域。這是因為浸潤面在斜坡底部較快到達防滲層（在本研究中防滲層為支撐上部邊坡的堅實土層），隨后形成水位使飽和度Sr值增加。從Δuw的變化可以看出，在降雨75 min左右，斜坡末端（P07）的Δuw上升，隨著降雨的滲透，Δuw趨向于從邊坡底端向模型中心（從P07到P01）方向上升。在模型底部飽和度Sr=0.8時，邊坡中間部分（M12）的Sr值增大，當中部的Sr值和邊坡底端的Sr值大致相同時，發生初始變形。由此可知，造成邊坡初期變形的原因之一是由于邊坡底端和中部水位升高，而導致邊坡抗剪切強度下降[7]。

2灑水降雨時水位上升的模擬

2.1數值模型和條件

數值模型截面如圖6所示，數值模型邊界條件如圖7所示。

數值模型的尺寸與模型實驗尺寸相同。滲流數值分析邊界條件：模型上部作為降雨入滲邊界，支撐層的左端作為排水邊界，模型的最底部和右緣均假定為不透水邊界。降雨條件設置與實驗相同。數值分析的初始條件是當安裝在邊坡和支撐層邊界之間的壓力值發生變化（降雨250 min），即支撐層的水位達到70 mm時，相當于給支撐層底面施加70 mm的壓力水頭。與模型實驗土體相同壓實度（Dc）試樣的浸透排水實驗結果如圖8所示。

滲流分析中使用了水分特征曲線Van Genuchten模型[9]以及比滲透系數模型Mualen[10]。

θ=1+（α×ψ）n-m×（θs-θr）+θr，（1）

Kr=Se1/21-（1-Se1/m）m2，（2）

m=1-1/n，（3）

式中：θ為體積含水率;θr為殘留體積含水率;θs為飽和時體積含水率;ψ為基質吸力，kPa;Se為有效飽和度;Kr為比滲透系數。

2.2數值分析

灑水降雨期間壓力計高度hm隨時間的變化如圖9所示。這里表示了使用與模型實驗干密度相同的試樣進行定水頭滲透實驗得到的滲透系數k（k=klab）和使用與實驗結果整合后逆向計算得到的滲透系數（k=kmod）的對比結果。

采用定水頭滲透實驗結果的分析中，Ma01的數值結果始終大于實驗結果，Ma02在t=400 min左右的再現性較高，之后數值結果遠大于實驗結果。在t=400 min之前，Ma03的數值結果略低于實驗結果，此后數值解比Ma02中的實驗結果大很多。是因為所提供的雨水沒有及時從模型邊坡中排出，從而導致模型中的水位被過高計算。通過使用表1所示kmod數值分析，可以很好再現實驗結果（圖9中k=kmod）。

滲透坍塌部分的土體會從坡面流出。現推測是由于土體壓實不均勻在坡體中產生了細小過水通道進而導致局部排水，這些過水通道的出現，導致整個模型實驗的滲透系數大于滲透實驗結果。另外，考慮到底部礫石材料的影響，使邊坡底部滲流特性發生變化，數值分析結果與實驗值相比大幅減小，難以取得整合性，從而對模型邊坡整體的滲透系數產生較大影響。圖10為坡面發生變形時，通過滲透系數值分析得到的流速矢量。

模型降雨過程中，在邊坡發生變形區域附近，流速矢量在邊坡表面變形位置的下方開始變大，且流速矢量沿水平向外變化，說明邊坡底部發生初期變形可能是受滲透壓作用的影響。

3彈塑性有限元對初期變形的再現模擬

3.1數值模型和條件

圖11為數值模型截面的網格視圖以及邊界條件。

由圖可知，邊界條件是左右兩端的垂直方向為自由端，水平方向為固定端。此外，分析條件將邊坡水位設定為初期變形時從壓力計上得到的數值。在此水位狀態下，先進行自重應力分析，然后根據強度折減法確定安全系數和變形行為。巖土材料的本構方程遵從Mohr-Coulomb準則。非飽和土的有效應力根據Bishop[11]提出：

σ′=（σ-ua）+xg（ua-uw）=

σnet+xgs，（4）

式中：σ′為有效應力;σ為全應力;ua為孔隙氣壓;uw為孔隙水壓;s為吸力;x為Bishop有效應力系數，0

數值模擬中使用的物理參數如表2所示。

邊坡發生初期變形前后的狀態如圖12所示。

可以觀察到邊坡底部設置的礫石向末端方向移動。此時坡面表層的位移量非常小，水平方向的位移量為0.5～1.5 mm。在本研究中，通過考慮微應變到小應變區域的變形系數來再現圖12所示A點的位移實驗結果。

巖土材料的變形系數隨著有效應力的減小和應變的增大而減小并表現出強烈的非線性[12-13]。變形系數的大小對應力和應變分析很重要，因為會影響初期模擬計算。本研究中數值模型的微應變至小應變的變形系數可通過常規三軸壓縮實驗獲得。本研究采用經驗公式（6）和公式（7）[14]計算了邊坡模型和支撐地層彈性區域的變形系數。

Gmax=1 400×f（e）×σ′m0.4，（6）

f（e）=（2.17-e）2/（1+e），（7）

式中：Gmax為微應變區的剪切模量;f（e）為孔隙比函數;σ′m為平均有效應力;e為孔隙比，由邊坡和下部支撐土層的ρd和ρs計算得出。

3.2數值模擬結果分析

剪切模量變化差值ΔGmax與坡面變形量的關系如圖13所示。

圖10和圖12中坡面變形，A點產生的水平位移為ux。在0.7Gmax左右時實驗與數值模擬結果基本一致，表示可以通過降低無水位狀態時的Gmax值來預測邊坡初期變形狀態。非飽和狀態下，Gmax的增減除了受滲透壓的影響外，也需考慮由于飽和度增加吸力降低導致有效應力減少的可能性[15]。

采用0.7Gmax得到的邊坡底部變形的數值分析結果如圖14所示。

數值模擬的邊坡底部發生變形，表明數值分析結果與實驗結果一致。在上述模擬中，假設變形系數減小的區域是整個模型邊坡，為了確定模型的平均變形系數，式（6）中的垂直有效應力取值為模型中心所受垂直應力（H=300 mm）。

水平位移分布的模擬結果如圖15所示。發生初期變形后持續降雨，650 min后發生大面積崩塌。發生大面積崩塌時模型內水位相比圖15的水位仍處于上升狀態。從圖15可以看出，模擬結果發現發生大變形的區域與圖16所示實驗中的最終破壞區域基本一致。因此，在初期變形發生時，邊坡變形系數（剛度）下降引起的微小變形極有可能形成潛在的滑移面，從而導致大面積崩塌。邊坡出現初始變形后進一步降雨，地下水位上升引起的滲透壓作用使得有效應力減小，從而導致邊坡逐漸失穩。為了防止土體結構逐漸失穩而導致大面積崩塌，建議采取坡體排水等措施來防止邊坡發生初期變形。由于邊坡的初始微小變形很可能與潛在滑動面有關，因此，如果檢測到這種微小變形，就必須考慮邊坡發生崩塌或滑坡的可能性。

4結論

基于粉質砂土模型邊坡降雨實驗，通過土體飽和度和孔隙水壓力的變化分析了邊坡失穩的過程，通過飽和及非飽和滲流模擬再現邊坡水位變化過程，采用彈塑性有限元法再現邊坡初始變形狀態，具體結論如下。

1）通過對滲流過程中地下水位上升進行分析，并使用彈塑性有限元法對模型實驗隨水位上升的初期變形進行數值模擬，發現：必須使用比滲透實驗得到的滲透系數更大的數值進行計算，才能在數值分析中再現降雨過程中的邊坡水位變化情況。

2）在模擬邊坡變形系數降低的彈塑性有限元分析中，通過將彈性區域的變形系數降低30%得到與實驗值相一致的數值解，由此可以看出，可以采用降低變形系數的數值模擬方法來預測邊坡初期變形狀態。

3）當前研究結果僅限于模型實驗的再現且數值模擬基于邊坡整體變形系數降低，而實際邊坡破壞時，在滲透壓作用下只有邊坡底部局部變形系數的降低。

在今后研究中，需要對降雨實驗和實際因降雨引起滑坡的案例進行更詳細的數值分析，在此基礎上建立更準確的邊坡崩塌預測方法，提出更有效的滑坡防治措施。

參考文獻/References：

[1]ORENSE R P.Slope failures triggered by heavy rainfall[J].Philippine Engineering Journal，2004，25（2）：73-90.

[2]TOHARI A，NISHIGAKI M，KOMATSU M.Laboratory rainfall-induced slope failure with moisture content measurement[J].Journal of Geotechnical and Geoenvironmental Engineering，2007，133（5）：575-587.

[3]張明禮，溫智，董建華，等.考慮降雨作用的多年凍土區不同地表土質活動層水熱過程差異分析[J].巖土力學，2020，41（5）：1549-1559.

ZHANG Mingli，WEN Zhi，DONG Jianhua，et al.Response of hydrothermal activity in different types of soil at ground surface to rainfall in permafrost region[J].Rock and Soil Mechanics，2020，41（5）：1549-1559.

[4]羅中寶.液壓離心地震模擬振動臺模型的建立及其控制策略研究[D].哈爾濱：哈爾濱工業大學，2015.

LUO Zhongbao.Modelling and Control Strategr Study of Hydraulic Centrifuge Shaker Used in Earthquake Simulation[D]Harbin：Harbin Institute of Technology，2015.

[5]LING H，LING H I.Centrifuge model simulations of rainfall-induced slope instability[J].Journal of Geotechnical and Geoenvironmental Engineering，2012，138（9）：1151-1157.

[6]胡嫚，吳飛，汪時機，等.基于光滑粒子法邊坡失穩影響范圍的模擬[J].重慶大學學報，2019，42（5）：56-65.

HU Man，WU Fei，WANG Shiji，et al.Modeling of influenced area after slope failure based on smoothed particle hydrodynamics （SPH）[J].Journal of Chongqing University，2019，42（5）：56-65.

[7]包一丁.排土場邊坡穩定性評價及潛在滑坡災害范圍預測——以朱家包包排土場為例[D].長春：吉林大學，2020.

BAO Yiding.Stability Analysis and Landslide Failure Potential Forecasting of the Waste Dump：Zhujiabaobao Waste Dump Case Study[D]Changchun：Jilin University，2020.

[8]BUI H H，FUKAGAWA R，SAKO K，et al.Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics（SPH）[J].Géotechnique，2011，61（7）：565-574.

[9]VAN GENUCHTEN M T.A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal，1980，44（5）：892-898.

[10]MUALEM Y.A new model for predicting the hydraulic conductivity of unsaturated porous media[J].Water Resources Research，1976，12（3）：513-522.

[11]BISHOP A W，BLIGHT G E.Some aspects of effective stress in saturated and partly saturated soils[J].Geotechnique，1963，13（3）：177-197.

[12]TATSUOKA F.Laboratory stress-strain tests for developments in geotechnical engineering research and practice[C]//Proceedings of the Fifth International Symposium on Deformation Characteristics of Geomaterials.Amsterdam：IOS Press，2011：3-50.

[13]TATSUOKA F，IWASAKI T，FUKUSHIMA S，et al.Stress conditions and stress histories affecting shear modulus and damping of sand under cyclic loading[J].Soils and Foundations，1979，19（2）：29-43.

[14]IWASAKI T，TATSUOKA F.Effects of grain size and grading on dynamic shear moduli of sands[J].Soils and Foundations，1977，17（3）：19-35.

[15]KAWAJIRI S，SHIBUYA S，KATO S，et al.Effects of matric suction on elastic shear modulus for three unsaturated soils[C]//Proceedings of the 5th International Conference on Unsaturated Soils.[S.l.]：[S.n.]，2010：271-275.