基于有限元模擬的汨羅市泥石流發生機理研究及該區域風化砂土剪切特性考察

周 飛，王 磊

(山東正元地球物理信息技術有限公司，山東 濟南 250101)

0 引言

泥石流是我國僅次于滑坡、崩塌的第三大地質災害，嚴重影響了山區居民的生命財產安全、山區資源可持續發展以及山區工程的建設和維護，每年造成直接經濟損失達數十億元[1-4]。2016年7月4日汨羅市川山坪鎮玉池山村玉明片區陳家山組后山出現了多處泥石流，導致上千群眾受災，直接受災損失高達782萬元。泥石流成災過程具有較強的偶然性和復雜性，時程變化特性顯著，野外原型觀測與模型實驗存在較大困難，因此，構建合理的數值模型來反演和再現泥石流過程的復雜物理現象逐步成為泥石流研究的有效手段?；陔x散介質的網格方法通常忽略了泥石流運動中的流體特性，僅適用于含水較少的碎屑流模擬[5-6]。基于連續介質網格法發展較為成熟，與GIS平臺以及目前通用的數字高程模型結合程度較好，便于各類泥石流本構模型與流量、侵蝕模型整合[7-8]。

經野外調查發現，在汨羅市泥石流發生的源頭區域，土體的種類構成以及其粒徑分布往往具備不均勻性，而滲透性的差異大量存在，如風化巖和堆積物，其具有高滲透性，因此與透水性能較差的巖層有一定的差異，此時不能滲入基巖的雨水在滲透和非滲透層的邊界流動，并導致流入土體和地面塌陷的現象，參與到原始泥石流沉積物中形成大規模的地質災害[9-11]。而大量的調查報告表明，管流孔道的痕跡大量出現在泥石流的源頭部，這可能是后期導致泥石流發生的重要原因。管流的特性，流動和水壓對安全系數和邊坡崩壞的影響尤為顯著，也是當前需要不斷進行深入研究的著力點。當前大量的研究基于有限元法的理論基礎，通過滲流分析對存在管流的坡面的地下水流向進行了討論[12-14]。研究結果表明：管道中粗粒部分和孔道的堵塞引起的局部水壓變化可能導致坡面塌陷，但孔道本身的結構和管道中泥沙的運動過程的結論尚未統一。

本研究中，暫時不對管流現象的產生與否進行討論，利用有限元數值模擬，從壓力水頭，總水頭差，流速等角度對邊坡破壞機理進行詳細分析探討，對汨羅市玉池山村泥石流發生進行了有效再現和破壞機理分析。同時通過室內剪切試驗對上部軟弱風化土層的物理力學性質進行了分析，該分析對泥石流源頭滑坡安全率理論計算具有重要作用，同時對該區域地質災害預警和防治具有重要意義。

1 有限元滲透模擬

1.1 概述

本節通過有限元滲流模擬，以把握滲流速度和地下水水頭，分析崩塌原因。在數值模擬中，把風化土層厚度、邊坡坡度、表層條件作為變量。

1.2 泥石流特征整理

以汨羅市玉池山村泥石流的斜坡為模型，調查發現該斜坡具有以下特點：

地質條件：發生泥石流的地區以中生代白堊紀形成的花崗巖為主?；◢弾r硬巖暴露在附近的山區，砂土化的風化堆積物稀薄地分布在山谷中。

地形：中游的河道為約5°的緩坡，裸露出花崗巖基巖(圖1a)。泥石流末端區部分堆積物主要由砂土組成，即使在2°～3°以下的緩坡，也發生了250～400 m的流動距離(圖1b)。

源頭區域狀況：泥石流源頭部崩塌平均深度為1.5 m或更小，坡度多為25°～40°，尤其集中在30°～35°范圍內。在風化層與基巖的交界處，有粗粒沉積物分布，說明滲透性較強(圖1c)。

圖1 區域現場拍攝圖像Fig.1 Field images in the area

1.3 數值方法、模型和案例

根據1.2泥石流特征，建立邊坡模型，采用有限元法進行二維滲流數值模擬。

1)土層構成和初期條件設定。通過對泥石流發生邊坡的野外勘察，發現表層土沒有或者很薄。其次，由于發現風化的砂土層與花崗巖基巖的交界處夾雜一層具有高滲透性的粗粒層，推測為以往泥石流的堆積物。數值模擬將分別在風化砂土層上有/無表土層的情況、風化砂土層部分分布表土層的情況、風化層與基巖之間有/無高滲透層的情況結合進行分析。關于非飽和滲透性，圖2中表示了體積含水率θ與負壓水頭ψ的關系、體積含水率θ與比滲透系數Kr的關系。各土層的滲透特性見表1。采用線彈性模型進行數值分析，風化砂土的彈性模量設為10 MPa，泊松比為0.3，條件為排水。由圖2和表1可知，在深度方向的吸力分布設置為0，并且未設置初始地下水位。

圖2 非飽和滲透特性Fig.2 Unsaturated permeability characteristics

表1 各土層特性Table 1 Characteristics of each soil layer

2)坡度。參考發生泥石流的狀況，將斜面坡度設置為25°，中游坡度設置為5°，不易發生泥石流的條件。

3)降雨量與邊界條件。將表2中從2016年7月4日湖南省汨羅市觀測到的降雨量作為降雨條件。在這種降雨條件下，泥石流產生于上午7：00至8：00之間。由圖3可見，給模型的左右兩邊設置一定水頭(右邊20.0 m，左邊至表土下)，模型底部為排水面。在本數值分析中，采用恒定流量滲透邊界。為了便于與真實降雨變化情況進行對比分析，默認模擬開始時間與真實降雨開始時間一致，即為凌晨1：00，之后進行時間累計分析。

表2 2016年7月4日汨羅市山區降雨Table 2 Rainfall in the mountainous area of Miluo City on July 4，2016

根據現場勘察情況，對表3和圖3中所示的5種情況進行了數值分析。

圖3 邊坡模型Fig.3 Slope model

表3 有限元數值分析案例Table 3 Finite element numerical analysis case

2.4 數值模擬結果

分別對表3中五類工況進行分析，邊坡模型為長55 m、地下水位線右側20 m處。圖4為各次降雨開始后5.5 h的上午6：30邊坡內水頭分布情況。由圖4a可見，如果有表層土壤但沒有高滲透層，即使給以大量降雨，壓力水頭也不會隨著時間的推移而發生劇烈變化。

另一方面，表土完全不存在的工況b中，圖4b中，壓力水頭為0，即地下水位不規則地出現在邊坡中，這種情況已經無法用達西法則來解釋說明。這種狀況是從凌晨2：00的降雨引起的，無法解釋7：00到8：00泥石流的發生。

圖4 有限元模擬結果Fig.4 Finite element simulation result

在缺少部分表土的工況c和工況e情況下，殘留表土在這兩種情況下都會抑制滲透，斜坡內的滲流相對穩定。在工況c和工況d中，無論有無表土，在高滲透層中都會發生比其他情況更快的滲流，壓力水頭隨也會時間增加而增加。特別是圖5中工況d時，4：30以后，高滲透性區域的壓力水頭逐漸增大。在分析初期，2：00時水頭分布的不自然被看作是過渡現象，之后在3：00時開始變得穩定。

圖5 工況d的經時變化Fig.5 Time variation in case d

表4為風化砂土層內部的流速、水頭、安全率的計算結果。實際流速為觀測1、2點的平均流速除以有效孔隙率的值。工況c和工況d的實際流速達到0.06 cm/s。本文采用經典的Taylor(1948)[15]極限流速ν計算公式(1)，得到ν=1.78 cm/s。

(1)

其中：g為重力加速度；γs為土的浮重度；A為受到滲流作用土的橫截面積；γw為水重度；Gs為土粒相對密度，取值2.62；d為粒徑，D10= 0.003 cm。

可以推斷出風化砂土層中的流速足夠緩慢和穩定，不會因流速而導致滲透破壞。圖4中兩點(1點和2點)流土現象安全系數F用公式(2)計算[16]。壓力水頭最大的點為第2點，距第2點水平距離為5～10 m處上游點為第1點，計算結果見表4。

(2)

其中：e為孔隙比，0.8；ha為壓力水頭差；D為位置水頭差。

由表4可知，只有在工況d中，沒有表土，砂土層厚1.0 m，交界處有高滲透層，F才略低于1.00。第一次強降雨后的4：30分，F≥1.00，但第二次強降雨后的6：30，F<1.00，與泥石流發生時間一致。雖然沒有因流速而發生滲流破壞，但認為由于壓力水頭的增加，邊坡穩定性有所下降。

表4 觀測點1和2的流速、水頭、安全率Table 4 Flow velocity， water head and safety rate at observation points 1 and 2

3 風化花崗巖剪切特性考察

對于泥石流源頭多處的邊坡崩塌，本文數值分析從水頭變化等角度對邊坡穩定性進行了討論，但仍需要從力學性質角度對邊坡構成土層進行力學分析以及進下一步的邊坡安全率計算等。本節通過直剪試驗對汨羅市泥石流發生區域風化花崗巖滑動過程中的剪切特性進行了初步判定。

3.1 野外調查和土樣采集

土樣采集點位于汨羅市川山坪鎮玉池山村玉明片區陳家山組后山邊坡，坡度約30°。表層花崗巖風化強烈，地層組成以砂土為主。土樣采集圖片見圖6，圖7為原狀土采集土筒的尺寸。

圖7 原狀土采集土筒Fig.7 Undisturbed soil collection tube

3.2 剪切儀概況

剪切儀見圖8，土樣盒直徑6 cm，高2 cm，內部光滑，分為上下兩層。剪切前進行超過20 min的壓密，試驗開始后維持同一垂直應力σ不變，保持0.2 mm/min的剪切速度，每隔1 min記錄剪切應力τ，垂直位移ΔH和剪切位移δ。當剪切位移達到7 mm時終止試驗(剪切儀剪切位移上限)。

圖8 室內剪切儀Fig.8 Laboratory shear equipment

3.3 風化層物理性質

風化層物理性質詳見表5，風化土層原狀土土樣抽取3個進行剪切試驗。

表5 風化層物理性質Table 1 Physical property of weathered layer

3.4 剪切試驗結果

圖9中展現了風化土在自然含水比(Sr=37%)下，剪切應力τ、垂直位移ΔH和剪切位移δ關系。剪切過程中剪切應力τ沒有出現明顯峰值后下降的現象，到達一定值后趨于穩定。

圖9 剪切應力-剪切位移關系(a)，垂直位移(b)-剪切位移關系Fig.7 Relationship between shear stress and displacement (a)， vertical displacement and shear displacement (b)

從體積變化角度來看，最主要展現了剪切過程中土樣體積收縮的特性。剪切應力τ和垂直應力σ的關系見圖10。

圖10 剪切應力和垂直應力關系Fig.10 Relationship between shear stress and vertical stress

4 有限元模擬和剪切試驗結果考察

在有限元分析結果中，每種工況下都得到了隨著降雨的持續而出現的壓頭升高情況，但在工況d中，坡面無表土，且透水性有差異，即基底與砂層之間夾雜高滲透層，隨著降雨量和時間的增加，觀測點2的壓力水頭逐漸增大。本案例中，在計算安全系數時，在上午6：30，即降雨后5.5 h，得到F<1.00，推測有流土現象發生，可能由于以下原因造成：

首先，如果基底的巖層非透水，而表層土又不能抑制滲透，那么強降雨時，約1.0 m厚的薄砂土層會很快飽和。此外，由于表土沒有覆蓋整個坡面，且與基底的邊界附近分布著粒徑相對較大的高滲透層，隨著降雨的不斷進行，導致壓力水頭升高。到了一定程度，壓力水頭達到所謂的流土狀態，土顆粒失去凝聚力，從而進一步導致崩壞發生。

在后續風化花崗巖土樣剪切試驗考察部分，初步得到了風化土層物理參數及其剪切特性，為今后邊坡穩定性分析和理論研究提供重要參數支持。

5 結論

本研究通過數值分析在不同土層分布條件下，對汨羅市玉池山村泥石流發生進行了有效再現和破壞機理分析，有限元數值模擬和剪切試驗得出的結論如下：

1) 對2016年7月4日汨羅市玉池山村泥石流災害發生源頭邊坡進行數值分析，孔隙水壓隨土層構成不同而變化較大，當砂土層和不透水層之間存在粗粒高滲透層時，觀察到孔隙水壓急劇增加，造成坡穩定性下降。

2) 在無表層土和存在粗粒高滲透層的情況下，在第二次強降雨后，孔隙水壓上升導致安全率小于1，邊坡失穩。

3) 作為數值模擬水壓變化的補充，進行了部分土的力學性質實驗。在剪切試驗中，風化花崗巖展現出剪切收縮特性，并初步得到了該風化土層物理力學性質，為下一步邊坡穩定性分析奠定力學基礎。

本研究對該區域地質災害預警和防治具有重要意義，對于夾雜高滲透層導致管流和邊坡流土發生的情況，今后研究中仍需邊坡模型試驗對數值模擬進一步對比分析。