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

周 飛，王 磊

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

0 引言

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

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

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

1 有限元滲透模擬

1.1 概述

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

1.2 泥石流特征整理

以汨羅市玉池山村泥石流的斜坡為模型，調(diào)查發(fā)現(xiàn)該斜坡具有以下特點：

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

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

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

圖1 區(qū)域現(xiàn)場拍攝圖像Fig.1 Field images in the area

1.3 數(shù)值方法、模型和案例

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

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

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

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

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

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

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

根據(jù)現(xiàn)場勘察情況，對表3和圖3中所示的5種情況進行了數(shù)值分析。

圖3 邊坡模型Fig.3 Slope model

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

2.4 數(shù)值模擬結(jié)果

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

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

圖4 有限元模擬結(jié)果Fig.4 Finite element simulation result

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

圖5 工況d的經(jīng)時變化Fig.5 Time variation in case d

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

(1)

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

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

(2)

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

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

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

3 風(fēng)化花崗巖剪切特性考察

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

3.1 野外調(diào)查和土樣采集

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

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

3.2 剪切儀概況

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

圖8 室內(nèi)剪切儀Fig.8 Laboratory shear equipment

3.3 風(fēng)化層物理性質(zhì)

風(fēng)化層物理性質(zhì)詳見表5，風(fēng)化土層原狀土土樣抽取3個進行剪切試驗。

表5 風(fēng)化層物理性質(zhì)Table 1 Physical property of weathered layer

3.4 剪切試驗結(jié)果

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

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

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

圖10 剪切應(yīng)力和垂直應(yīng)力關(guān)系Fig.10 Relationship between shear stress and vertical stress

4 有限元模擬和剪切試驗結(jié)果考察

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

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

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

5 結(jié)論

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

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

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

3) 作為數(shù)值模擬水壓變化的補充，進行了部分土的力學(xué)性質(zhì)實驗。在剪切試驗中，風(fēng)化花崗巖展現(xiàn)出剪切收縮特性，并初步得到了該風(fēng)化土層物理力學(xué)性質(zhì)，為下一步邊坡穩(wěn)定性分析奠定力學(xué)基礎(chǔ)。

本研究對該區(qū)域地質(zhì)災(zāi)害預(yù)警和防治具有重要意義，對于夾雜高滲透層導(dǎo)致管流和邊坡流土發(fā)生的情況，今后研究中仍需邊坡模型試驗對數(shù)值模擬進一步對比分析。