复合型滑坡固液耦合过程数值模拟分析——以无山坪滑坡为例

张 晗高 杨李 滨李 军吴伟乐

1.长安大学地质工程与测绘学院,陕西 西安 710054;

2.中国地质科学院地质力学研究所,北京 100081;

3.自然资源部活动构造与地质安全重点实验室,北京 100081

全球气候变化造成了极端天气事件的增多,进而导致特大型地质灾害频发。强降雨条件已成为诱发滑坡及复合型地质灾害的最主要原因之一(孙广忠,1988;殷跃平,1998;李滨等,2016;高杨等,2017;闫金凯等,2020)。复合型地质灾害是指包含两种或两种以上基本地质灾害类型,并可引发堰塞坝和涌浪等次生效应的链式地质灾害(Varnes,1978;Yin et al.,2009;殷跃平等,2010,2021;Hungr et al.,2014)。在中国西南山区,高位崩滑体失稳启动后,经过势动能转化不断碰撞解体形成碎屑流,往往在水动力作用下极易向泥石流或水石流转化,形成复合型地质灾害,对人类的生命和财产安全造成巨大的威胁。针对水动力作用造成滑坡远程运动理论,国内外学者开展了大量研究。滑坡碎屑流在水动力条件下向泥石流的转换现象已经被大量地质工作者通过野外调查和试验证明(Plafker and Ericksen,1978;Voight et al.,1983;崔鹏,1991;崔鹏和关君蔚,1993;Iverson,1997;Takarada et al.,1999;余国安,2022)。固体和流体不同的相互作用不仅可以划分滑坡或泥石流的类型,还使这类灾害具有独特的破坏性。碎屑流-泥石流复合型滑坡的流变试验和数值研究表明不同液体体积分数、液体性质、固液速度差和固体颗粒尺寸是影响固液相互作用,从而影响滑坡运动距离的关键因素(Gao et al.,2022a)。流体存在会增加孔隙水压力,有效降低颗粒间摩擦并减少能量损耗,提高了滑坡的运动特性。对比流体拖曳力,具有黏性和不可压缩性的液体存在使得运动滑体的边界层保有较高的孔隙水压力,在高速剪切、振动作用下形成超孔隙水压力,导致滑体摩擦阻力下降,增加滑坡运动距离(Hutchinson and Bhandari,1971;Evans et al.,2001;Sassa et al.,2004;殷跃平等,2010;Iverson and George,2016)。Strom and Abdrakhmatov(2018)在进行大量滑坡调查后发现,滑坡高速运动造成的饱水下垫层液化是其侧向扩离运动及覆盖范围增大的重要原因。碎屑流向泥石流转化形成的复合型滑坡在世界范围内多次发生,给人类生命财产安全带来巨大威胁,例如:1970年 秘 鲁Yugay滑 坡(Plafker and Ericksen,1978),2002年 俄 罗 斯Kolka滑 坡(Evans et al.,2009),2006年菲律宾Leyte滑坡(Evans et al.,2007),2008年四川文家沟滑坡(余斌等,2010),2014年重庆咸池水库滑坡(卫童瑶,2021;Gao et al.,2022a;Li et al.,2022),2019年贵州鸡场滑坡(Gao et al.,2020;李壮等,2020;高浩源,2021)。

为开展高位滑坡后破坏成灾范围和成灾风险的定量评估,进行科学避险,数值模拟技术成为反演分析滑坡运动过程和动力学机制的最为高效的手段之一。目前研究两相复合型滑坡的数值模拟方法包括欧拉-欧拉模型(Reynolds,1976;Iverson and Denlinger,2001;Pitman and Le,2005;Pailha and Pouliquen,2009;George and Iverson,2011;Bouchut et al.,2015,2016)、拉格朗日-拉格朗日模型(DEM-SPH等;Tan and Chen,2017;徐文 杰,2020)、欧 拉-拉 格 朗 日 模 型(SPHDCDEM、DEM-CFD等;Shan and Zhao,2014;Jing et al.,2019)。欧拉-拉格朗日模型较其他模型能够更好地还原固液两相流的运动过程,但由于现有数值模型没有充分考虑拖曳效应及技术限制,还未实现液体和固体颗粒的共同运动。

2014年9月2日晚,重庆市奉节县竹园镇无山村所在山体发生高位滑坡(以下简称无山坪滑坡),滑坡体内房屋全部倒塌,共1504间房屋被毁坏。由于及时的降雨警报,无山村所有居民在山体滑坡发生前20分钟被疏散。滑坡沿岩层层面产生整体位移,前缘滑体大面积滑塌,滑塌的岩土体混和暴雨形成泥石流,沿沟道冲入并堵塞岔河,属典型的滑坡转化为泥石流的复合型滑坡。文中采用基于两相耦合模型和SPH算法技术结合的滑坡后破坏数值模拟平台(LPF3D),针对无山坪滑坡,采用现场调查、遥感影像对比和数值模拟等方法对该滑坡的水力学特征进行了初步探讨,分析了滑坡基本情况和成灾机理,模拟并重现了滑坡运动堆积过程,为类似固液耦合复合型滑坡的预测和反演提供帮助。

研究区位于重庆市东北部奉节县竹园镇,距离奉节主县城约80 km,中心地理坐标为109°01′17″E,30°29′19″N。研究区滑坡前后缘公路交通条件较好,与乡镇、奉节县城相连。滑坡发生后公路损坏严重。地形地貌属构造剥蚀斜坡、中低山地貌,地形总体呈南高北低。滑坡后缘最高点海拔为1330 m,前缘最低点位于岔河河谷,海拔约800 m,相对高差550 m;地形坡度一般为15°～35°,多陡坎。山体表面发育多条冲沟,滑体失稳下滑后,沿运动路径中的冲沟滑动,并在前缘沿东西向岔河内堆积。

研究区位于新华夏系第三沉降带四川盆地东端,地处川东坳褶带、大巴山南缘弧形褶皱带及鄂湘黔隆褶带接合部。地质构造以褶皱为主,断裂欠发育。滑坡处于红岩向斜轴部近西端(图1)。由于早期的构造作用和风化作用,部分岩体破碎,有利于雨水的渗透和地下水的储存。

图1 研究区地质构造图Fig.1 Geological structural map of the study area

滑坡区域覆盖第四系松散堆积物,主要由黏性土、砂泥岩碎块石和少量生活及建筑垃圾组成。基岩为侏罗系中统新田沟组(J2x)的砂岩和泥岩,上部以砂岩为主,中部砂、泥岩互层或混杂沉积成岩,下部砂岩较泥岩发育。从无山坪滑坡的后缘到前缘,地层的倾角逐渐减小到零;在滑坡的前端甚至出现倒置的地层。该区砂岩强度较高,容易沿着岩层的倾角发生缓慢变形,从而在斜坡表面逐渐形成许多卸荷裂缝,这些裂缝加速了地表水和雨水的渗入(殷跃平和胡瑞林,2004;刘新荣等,2008;王志俭等,2008;刘小红等,2015)。砂岩和泥岩之间的薄弱界面和广泛拉张裂缝的发展有利于雨水的渗透,促进了滑坡形成(Li et al.,2022)。

研究区多年年平均降雨量由川鄂交界处的2000～2100 mm及北部边缘的1500～1600 mm向长江河谷降低为1100 mm左右,其中3—8月降雨量占全年的68%以上。滑坡发生一周内,重庆市奉节县竹园镇遭到百年一遇最大集中降雨,奉节县气象局数据显示奉节县竹园镇累计降雨量可达266 mm。短时强降雨是包括无山坪滑坡在内的许多地质灾害的主要诱发因素。

2014年8月31至9月1日,重庆市奉节县竹园镇累计降雨量达到266 mm,数小时内诱使无山坪滑坡发生整体位移,并出现大面积滑塌,滑塌的岩土体同地表径流混和形成泥石流,沿沟槽冲入岔河,堵塞岔河形成两处堰塞湖。无山坪滑坡目前仍残余大量潜在不稳定体,平台前部与右侧共发育有五条沟道,坡体内发育有大量拉张裂隙,若遇强降雨,滑坡极有可能发生再次下滑,威胁滑坡下方岔河两岸居民及基础设施安全。

通过现场调查、卫星和无人机影像分析(图2),滑坡失稳后发生解体,沿N20°W和N8°E两个方向高位剪出,约4×106m3滑体残留在滑源区,2.83×106m3滑体向下游流态化运动堆积。强降雨在研究区内造成了丰富的地表径流,滑坡失稳启动后,滑体同丰富的地表径流混合形成泥石流复合灾害,泥石流沿着无山坪平台前部的五条沟道向下游流动,汇入坡脚处东西向的主河道。根据无山坪滑坡的堆积特征,从总体上将该滑坡-泥石流分为滑源区、流通区和堆积区(图3)。

图2 无山坪滑坡滑前滑后遥感对比图Fig.2 Image comparison before and after the Wushanping landslide

(1)滑源区

滑源区的平面形态呈“舌型”,顺坡向滑坡体总体长约500 m,前源宽约820 m,后缘宽约360 m,平均厚度约19 m。平面面积约39.2×104m2,体积约6.83×106m3。滑坡运动过程中,滑体沿与岩层之间的滑面旋转和滑动,导致后缘垂直错位超过30 m,滑后可见砂、泥岩互层的基岩出露,与水平夹角约30°(图4a)。滑源区中部地形较高,在滑动过程中形成了许多长10～20 m的裂缝,伴随着不均匀的下沉。随处可见无山村被毁道路与房屋(图4b)。

(2)流通区

流通区主要位于高程1000 m以下的沟道处,斜坡基岩主要以第四系堆积物和砂泥岩为主。滑源区岩体失稳后,在巨大的势能和雨水拖曳的作用下,快速顺沟道向下运动(沟道位置见图3a)。流通区顺坡向长约400 m,宽约700 m,面积约28×104m2。滑坡剪出口位于一个陡峭斜坡的自由面,坡度为40°～50°,距离坡脚的高度为150～200 m(图4c、4d)。

图3 无山坪滑坡平剖面图Fig.3 Profile and plan of the Wushanping landslide

(3)堆积区

滑前在堆积区有一条梅溪河的支流岔河流过山谷,方向为160°。由于暴雨的缘故,滑体在陡坡上运动解体转化为泥石流,大量流固混合物沿沟道冲向河流。最终形成顺滑向670 m、顺河长250 m、平均厚度14 m、体积约为2.3×106m3的堰塞体,将岔河堵塞。此外,大量滑坡碎屑继续向下流动约1000 m。泥石流和堰塞湖淹没了原河谷两边的一些房屋(图4e、4f)。

图4 滑源区、流通区及堆积区现场调查照片Fig.4 Site photos of the slide source area, propagation area and accumulation area

针对此类两相流复合型滑坡,文中基于SPH方法和两相耦合动力学理论研发的LPF3D(Landslides post failure 3D)数值计算方法,能够处理具有自由表面、变形边界、运动交界以及极大变形的问题,具有计算速度快和计算精度高的适用性。流体采用基于连续介质的SPH方法计算,SPH粒子承载着质量、密度、速度、加速度等物理量,粒子间通过核函数相互作用,适用于大变形情况。固体颗粒体基于颗粒动力学模型将传统SPH方法改造成适用于离散颗粒的SPH方法,采用理想弹塑性本构进行求解,使得SPH粒子不仅承载质量、密度、速度、加速度等各种物理量,同时承载着颗粒的粒径、体积分数等颗粒属性(Cui et al.,2021;Chen and Yan,2021)。

3.1 控制方程

不同状态下滑体材料的控制方程由质量和动量守恒原理描述,见公式(1)和(2)。

连续方程:

运动方程:

公式中:α和β分别是x、y两个方向上的分量;ρ为材料密度;vα为速度在x方向上的分量;xα和xβ分别为x在x和y两个方向上分量;t为时间;f代表其他外力(如重力、相间作用力);d/dt为方程全导数。总张量σαβ通常分为各向同性压力P和剪应力τ两部分。表示为:

其中δαβ是克罗内克函数,当α=β时,δαβ=1;当α≠β时,δαβ=0。

3.2 计算本构

3.2.1 流体本构

该本构模型以计算纯水流体和浓缩流体为主,流体密度ρ与压强P关系为:

其中P0=ρ0c20/γ,为初始压强;ρ0为流体初始密度;γ与流体的可压缩性相关,γ=7;c0为初始声速,为保证流体可压缩性,一般取c0=(10～40)vmax,vmax为流体最大速度。

对于牛顿流体,剪应力ταβ(公式5)与剪应变率εαβ(公式6)成比例,比例系数为黏性系数μ;对于非牛顿流体,黏性系数μ为剪应变率函数。

其中·v为速度的散度,其他参数同公式(1)与公式(3)。

3.2.2 固体本构

静压力Ps是使用平均应力的标准定义从颗粒构成方程中直接计算的,表示为:

其中σxx、σyy和σzz是应力张量σγγ在x、y和z方向上的分量。当颗粒处于准静态状态时,颗粒整体几乎不产生大变形,主要以弹性变形为主,在此阶段,根据线弹性模型(胡克定律)计算颗粒的应力-应变关系。

其中σ·αβ为增量 形 式 的 应 力 分 量,G为剪 切 模量,K为弹性体积模量,E为弹性模量(杨氏模量),υ为泊松比,e·αβ为偏剪切应变率张量,ε·αβ为应变率张量。

ταβ为偏差应力分量;为第二不变量;为应变率张量,其定义与摩擦系数μ(I)相对应,这也意味着体积分数对惯性常数的单调性相关。

此处, 实验和数值模拟表明在函数μ(I)中,在非常低的惯性常数I(准静态)下,最小值μp逐渐增加到I增加时的有限值μ2。I0是一个常数,公式(4)中的参数取决于材料特性。例如,典型值时

μp=tan(21°),μ2=tan(33°),I0=0.28(MiDi, 2004)。

惯性常数表示惯性时间尺度(d2ρp/Ps)0.5与宏观变形时间尺度之间的比率;d表示颗粒粒径;ρp表示颗粒密度。

3.3 相间作用

大量研究表明,两相流运动中颗粒体和流体的相间作用对滑体运动范围至关重要(Davies,1990;Iverson et al.,2010;Tayyebi et al.,2022;Gao et al.,2022b)。计算固液相间拖曳力时,采用Gidaspow Gidaspow,提出的公式(1994),即对于密相的计算采用Ergun方程以及对于稀相的计算采用Wen-Yu方程(Ergun,1952;Wen and Yu,1966):

其中β为流体与固体间的动量传递系数;CD为曳力系数;φP为固体内摩擦角;μf为流体黏度;dP为固体颗粒粒径;ρf为流体密度;vf、vP分别代表流体和固体的速度;αf为流体体积分数。

曳力系数CD为:

相对雷诺数ReP定义为:

为消除两个方程间的不连续性,引入松弛因子φfP对过渡区域中的动量交换系数进行光滑,

因此,动量交换系数β可以表示为:

可得作用于单位质量颗粒上的曳力R′fP为:

其中参数含义同公式(16)。

3.4 边界条件

3.4.1 法向边界力

SPH作为无网格方法,不能像网格方法那样直接将界面力施加在边界点上。边界力求解采用强洪夫等(2011)提出的罚函数方法求解颗粒与边界之间的法向边界力。

3.4.2 切向边界力

滑体与基底的切向力根据滑体材料性质分别选择,流体和颗粒体分别采用库伦摩擦模型和层流黏滞模型(Hungr, 1995)。

固体-库仑摩擦模型:

公式中:T为基底剪切阻力;γ为重度;Hi为流体厚度;θ为运动路径坡角;ac=v2/R为离心加速度,大小取决于运动路径的曲率半径R,v为运动速度;g是重力加速度,取9.8 m·s-2;ru为孔隙水压力系数;φ为摩擦角。

流体-层流模型:

公式中:T为基底剪切阻力,A为面积,v为流体速度,μ为动态黏滞系数,H为流体深度。

4.1 模型建立及参数选取

为定量描述研究区地形地貌特征,采用1∶2000的三维等高线数据建立研究区地质模型,在滑源区建立滑体与流体模型。此次模拟不考虑前期滑体失稳过程,假设滑坡由降雨持续影响并突然爆发,目的是为了模拟流体与岩土体颗粒的运动堆积过程。

滑体材料物理力学参数采用实际参数。滑坡岩性主要为侏罗系中风化砂岩,岩石的物理力学参数为ρ=2240 kg/m3,粒径为0.1 m。流体采用泥浆参数,ρ=1200 kg/m3,黏滞系数为0.2 Pa·s。为探讨水动力对固液两相流滑坡后破坏过程的影响,此次模拟共设计四种不同工况(表1)。①工况Ⅰ:纯固体颗粒下滑,不考虑孔隙水压力和相间作用力。②工况Ⅱ:考虑孔隙水压力对下滑运动过程中固体颗粒的影响。③工况Ⅲ:考虑运动过程中流体对固体颗粒持续拖曳的影响。④工况Ⅳ:同时考虑流体孔隙水压力与拖曳效应的影响。

表1 无山坪滑坡LPF模拟参数Table 1 LPF simulation parameters of the wushanping landslide

4.2 结果分析

分别取颗粒运动停止时刻为计算终止时间。不同工况下无山坪滑坡模拟结果如图5所示。工况Ⅰ:纯固体颗粒下滑,大部分运动颗粒停留在沟道内或沟道出口处,均未抵达下方河道,最大运动速度为23 m/s,最远运动距离为900 m。工况Ⅱ:滑体底部孔隙水压力存在降低了摩擦,使固体颗粒运动距离较工况Ⅰ更远,大部分颗粒停留在沟槽出口位置,最大运动速度为26 m/s,最远运动距离为1150 m。工况Ⅲ:仅考虑流体拖曳作用时,固体颗粒与流体耦合向下流动,在各沟道出口处均形成了不同程度的冲积扇,具有明显的流态化特征,最大运动速度为31.5 m/s,最远运动距离为1250 m。工况Ⅳ:孔隙水压力减阻和流体拖曳力增程共同作用时,滑坡的运动距离更远,致灾范围更广,与无山坪滑坡灾后实际堆积距离、堆积范围最为接近,最大运动速度为34 m/s,最远运动距离为1300 m。

图5 四种工况下的堆积结果图Fig.5 Accumulations of the Wupingshan landslide under four working conditions

通过对固液耦合作用的对比,结合堆积模拟结果,可以发现工况Ⅳ得到的计算结果与实际最为接近。因此选用工况Ⅳ代表无山坪滑坡泥石流后破坏过程,并对其运动堆积特性进行论述。

文中选取0～40 s和180 s六个关键节点来分析无山坪滑坡破坏后的运动过程。工况Ⅳ下无山坪滑坡泥石流随时间变化的运动过程如图6所示,棕色颗粒表示固体颗粒,蓝色颗粒表示流体颗粒。整个运动过程中固液整体的运动速度如图7所示,在180 s时,颗粒的平均速度已无法影响滑坡的堆积状态,故180 s时运动已停止。从图中可以看出,在0～10 s之间滑坡高位启动,固体颗粒与流体混合顺多处沟道快速向下流动。10～20 s,滑坡前缘抵达下方处,大约在16 s左右,速度达到最高值34 m/s。20～40 s,由于泥浆的拖曳作用,后部固体颗粒沿着沟道继续运动不断汇入堆积区,堆积厚度不断增大。在40 s左右,大部分颗粒运动到河道位置,并在此处淤积,平均速度明显下降。在180 s左右,整体颗粒接近于静止运动,计算停止。

图6 工况Ⅳ下运动过程图Fig.6 Diagrams showing the fluid-solid coupled movement of the Wushanping landslide under working condition Ⅳ

图7 工况Ⅳ下运动速度图Fig.7 Velocity diagrams of the Wupingshan landslide under working condition Ⅳ

滑坡运动过程中固体堆积厚度变化等值线图如图8所示,图中红色线代表真实滑坡,红黄蓝三种色系体现固体颗粒不同堆积厚度变化情况。滑坡运动全过程,岩土体材料混合流体顺多条沟道向下流动,逐渐在各沟道处堆积。颗粒停止运动时,滑坡前缘在河道处呈流态化堆积,最大运动距离达1300 m,最大堆积厚度约21.5 m,堆积形态与现场调查结果较为一致。

图8 工况Ⅳ下滑坡堆积厚度图Fig.8 Diagrams showing the deposition thickness with time of the Wupingshan landslide under working condition Ⅳ

流体与固体颗粒间的速度差是决定拖曳力大小的关键因素,为进一步了解无山坪滑坡运动过程中拖曳力及孔隙水对固体颗粒运动的影响,选取②号沟道处两个速度对比曲线图作为对比。四种工况下滑坡前缘固体颗粒的速度对比曲线如图9a所示,水和相间作用力的影响下,固体颗粒的运动速度有明显的增大,且孔隙水作用大于相间作用。同时有流体参与运动时,滑坡运动时间更长,运动距离更远。工况Ⅲ、工况Ⅳ下滑坡前缘固体颗粒及其周围流体的速度对比曲线如图9b所示,滑动和加速阶段,工况Ⅲ和工况Ⅳ的流体速度相同,但固体颗粒速度有明显不同,故滑坡启动和加速主要受孔隙水压力的影响。滑坡减速阶段,两种工况下固体颗粒速度随流体速度减小,且变化趋势基本一致,故相间作用力在滑坡减速阶段起主要作用。

图9 无山坪滑坡速度曲线图Fig.9 Velocity change of the front-edge granules under different working conditions

文中用自主研发的LPF3D软件,以重庆奉节无山坪滑坡为例,模拟再现了无山坪滑坡运动堆积过程,探讨分析了远程滑坡后破坏运动过程中的水动力作用,认为流体会在孔隙水压力减阻和拖曳力增程两方面影响高位滑坡的远程运动,取得了以下认识。

(1)通过四种工况对无山坪滑坡泥石流的模拟结果分析,孔隙水压力和拖曳力的共同作用会使滑坡的运动速度、运动距离增大。表明水动力在滑坡运动过程中的动力学效应主要分为液化效应和拖曳效应两种,两种效应的共同作用会使碎屑流转化为泥石流,导致远程成灾。

(2)基于SPH的多相耦合计算模型,将流体状态方程、弹塑性本构方程和相间作用力进行搭建结合,模拟结果较好地再现了无山坪滑坡的运动堆积。表明基于多相耦合理论的LPF3D计算方法能够更真实地还原滑坡后破坏过程,对复合型滑坡的模拟具有较好适用性。该方法对于计算自由表面流、大变形和损伤破坏问题具有较大优势,避免了网格法在计算时遇到的界面追踪困难、网格扭曲和缠绕、计算量大等问题。

(3)模拟结果显示:无山坪滑坡全程最大运动速度为34 m/s,最大堆积厚度为21.5 m,最大堆积面积为0.12 km2,最远运动距离为1300 m。在前期强降雨影响下,无山坪滑坡滑体材料充分饱水,强度降低。同时降雨入渗补给使滑面处孔隙水无法即时排出,导致孔隙水压力增加,摩擦系数降低,滑坡高速启动,并在运动过程中与流体混合形成泥石流。流固两相间拖曳力的存在进一步促进滑坡颗粒流化运动,提高了其运动速度,使得滑坡远程成灾。

猜你喜欢 滑体泥石流滑坡 滑坡推力隐式解与显式解对比分析——以河北某膨胀土滑坡为例河北地质(2021年1期)2021-07-21基于遥感数据的灾后滑坡信息快速提取方法中国地质灾害与防治学报(2020年2期)2020-05-09泥石流杂文月刊(2018年21期)2019-01-05滑坡稳定性分析及处治方案中国公路(2017年18期)2018-01-23西南某水电站滑坡稳定性及其失稳影响研究资源环境与工程(2017年4期)2017-09-03露天矿区滑坡压煤后减少二次剥离的回采方法中国地质灾害与防治学报(2017年2期)2017-07-18“民谣泥石流”花粥:唱出自己海峡姐妹(2017年6期)2017-06-24泥石流环球时报(2017-06-14)2017-06-14露天矿反铲挖掘机处理滑体的方式中国科技纵横(2017年5期)2017-05-12浅谈公路滑坡治理北方交通(2016年12期)2017-01-15

推荐访问:滑坡 耦合 为例