基于数值模拟方法对船山村滑坡流固耦合计算分析

Fluid Structure Coupling Calculation and Analysis of Chuanshan Village Landslide Based on Numerical Simulation Method

摘要: 本文基于实际监测资料及勘查资料，获得拜殿乡船山村滑坡的降雨量、孔隙水压力及土壤含水率的参数，滑坡的变形发展受地下水位波动变化明显，地下水位上升通常是滑坡发生和发展的最重要影响因素，降雨对滑坡的影响主要通过地下水及其渗流场来发生作用。本文采用Geo-studio数值模拟软件，对于降雨条件下该滑坡土壤入渗规律进行模拟，分析滑坡在不同渗流场条件下的稳定性，根据数值模拟结果来判断滑坡发生预警级别。

Abstract: This paper is based on the actual monitoring data and exploration data. The parameters of rainfall, pore water pressure and soil moisture content of Chuanshan village landslide in Baidian township are obtained. The deformation development of landslide is obviously affected by the fluctuation of groundwater level. The rise of groundwater level is usually the most important factor affecting the occurrence and development of landslide. The influence of rainfall on landslide mainly occurs through groundwater and its seepage field. In this paper, geo studio numerical simulation software is used. The soil infiltration law of the landslide under rainfall conditions is simulated. The stability of landslide under different seepage field conditions is analyzed. According to the numerical simulation results, the early warning level of landslide is judged.

1. 滑坡地质环境条件

1.1. 地形地貌

该滑坡点位于衡阳市南岳区寿岳乡镇船山村枫树湾组，该滑坡所在区域属于低山地貌。该区域滑坡群断续延长1500 m，由多个滑坡体组成，一般长度50~100 m。

1.2. 地层岩性

根据地表调查和钻探揭露情况，滑坡区域地层结构简单，其上覆土层主要为第四系全新统残坡积(Q₄^del)碎石土，下伏基岩为燕山期(γ₅²)全风化~中风化花岗岩。其中：碎石土呈黄褐色，硬塑，土质不均匀，主要由粉质粘土及大量中砂、强风化花岗岩碎块组成，局部夹有块石，层厚7.60~12.70 m，分布广泛。

1.3. 水文地质条件

滑坡区地下水类型主要为第四系松散岩类孔隙水及花岗岩风化裂隙水两大类，富水量贫乏。地下水补给来源为大气降水与地表水，雨季时地表水补给地下水，旱季时地下水补给地表水。地下水流向受地形控制，一般从山脊、山坡流向溪沟，基本与地表径流方向一致。地下水的排泄方式有两种：一是沿基岩裂隙面与坡洪积层接触带呈线状、带状分散排泄，雨季或大雨过后，沿山坡、冲沟处是排泄场所；二是在冲沟下部地形低洼地段或坡洪积层陡坎处，顺基岩裂隙面以下降泉的方式进行集中排泄。

2. 滑坡体基本特征

2.1. 滑坡形态特征

该滑坡纵长200 m，横宽约250 m，厚约10~12 m，面积约5 × 10⁴ m²，体积约6 × 10⁵ m³，滑坡属于中层中型牵引式土质滑坡 [1]。

2.2. 滑体

根据钻探揭露，该滑坡体主要物质由碎石土及全风化花岗岩组成，平均厚度约11.7 m，主要由粉质粘土及大量中砂、花岗岩碎块组成，局部夹有块石。

2.3. 滑带(面)

岩性主要为全风化花岗岩，灰褐色，原岩结构及构造已基本破坏，岩芯已风化呈砂土状，土体呈软塑状，饱和，该层广泛分布，层厚较薄，揭露厚度为0.30~0.50 m。

2.4. 滑床

滑床岩层主要为强风化花岗岩层，基岩风化裂隙较发育，岩芯破碎，透水性一般，降水入渗后造成滑床部位水量相对富集，湿润滑带土，降低滑带土的工程力学性质。

3. 滑坡监测数据分析

3.1. 监测平台站点

该滑坡共布设监测点位12个，布设18件套设备。在ZK16附近安装一体化雨量自动监测站1套，在ZK15和ZK16中间位置安装一体化土壤含水率监测仪1套，在ZK20号钻孔位置安装一体化孔隙水压力自动监测站1套 [2]。

3.2. 降雨数据

在2017年5月1日至2017年7月31日92天间，监测平台无数据天数为13天，累计生产数据1508条。根据图1所示，降雨数据的产生数量基本由降雨决定，降雨量大的时候数据生产数量也相应多，由此可将降雨分为5月1日~24日、6月1日~7日、6月13日~7月2日、7月22~31日4个周期，并将其作为后期分析的基础。

Figure 1. Distribution map of rainfall data

图1. 降雨数据量分布图

对比湖南省气象台数据(表1)，可见两者变化趋势非常一致，但对比统计累计降雨量、该点数据值存在少18%的差值。

Table 1. Statistical cumulative rainfall in June 2017

表1. 2017年6月统计累计降雨量

3.3. 土壤含水率数据

自2017年5月1日至2017年7月31日，该点仅5月1~4日无数据，平均每日6条数据，共计547条、1641个数据。根据图2所示，1 m处：5月22日~6月12日、6月27日~7月31日数据下降到0.1附近波动；1.5 m处在6月25日~7月22日较大波动；0.5 m处在31%附近波动，该点数据需结合设备进一步探讨。

Figure 2. Soil moisture distribution

图2. 土壤含水率分布

3.4. 平均水渗压数据

自2017年5月1日至2017年7月31日，该点缺5月1~4日数据，共计产生537条数据，平均每日6条。从图3中可见，该点水渗压峰值与降水峰值响应情况良好。

3.5. 数据质量分析

1) 监测降雨数据质量最优。

2) 位移类数据变化较小，GNSS数据不完整；地表位裂缝移、深部测斜数据变化较小。

3) 土壤含水率、孔隙水渗压数据与降雨呈现一定相关性，土壤含水率数据需核实 [3]。

Figure 3. Pore water seepage pressure distribution

图3. 孔隙水渗压分布

4. 滑坡流固耦合计算分析

本文基于实际监测资料及相关勘查报告，获得拜殿乡船山村滑坡的降雨量、孔隙水压力、土壤含水率的参数，采用Geo-studio数值模拟软件，对于降雨条件下拜殿乡H2滑坡土壤入渗规律进行模拟，分析滑坡在不同渗流场条件下的稳定性。

4.1. 计算方法

地下水在滑坡内运动，由于滑坡形成过程复杂，导致孔隙的大小、形状都很复杂，水流在滑坡内运动相对紊乱，有些地方杂乱无章甚至不连续。为简化研究过程，依照已有土体内渗流分析方法，使用与真实水流属于同一流体的、并且充满整个含水层的假想流体来代替仅在孔隙中运动的真实水流，通过对这一假想流体的研究来达到了解真实水流渗透规律的目的。假想流体设定性质如下：

1) 它通过任意断面的流量与真实水流通过同一断面的流量相等；

2) 它在某断面上的压力或水头应等于真实水流的压力或水头；

3) 它在任意土体或岩体体积内所受到的阻力应等于真实水流所受到的阻力。

满足上述条件的假想水流称为渗透水流，一般简称渗流。它的基本表征量有两个，即流速及水头，都是空间坐标(x, y, z)和时间t的函数。

目前应用数值模拟进行滑坡地下水渗流分析的方法最为普遍，而最广泛的渗流场的数值模拟是有限元法的应用。有限元法的计算主要是解决渗流场中的水头函数，确定渗流场的自由面和渗流参数。有限元法即把渗流微分方程和边界条件按变分原理转变成一个泛函极值问题，首先把连续体或研究领域离散成有限元，并形成一组代数方程组，在计算机上解决。

有限元分析是用较简单的问题代替复杂的问题然后再求解，非常适合渗流这种复杂的问题。本滑坡的渗流计算基于饱和-非饱和理论方法，采用Geo-studio系列软件中的SEEP/W模块进行模拟。

降雨入渗是水在土体中流动的过程，其流动特征分为水在饱和土体中的流动和在非饱和土体中的流动。非饱和区土壤水流动和饱和区的地下水流动是互相联系的，将两者统一起来研究比较方便，即所谓饱和非饱和流动问题。

饱和土中水的流动通常用Darcy定律来表示，Darcy假设土体中水的流速与其水力梯度成正比：

$v_w=-k\frac{\partial h_w}{\partial y}$ (4-1)

式中：v_w为水的流速，k为土体的渗透系数， $\frac{\partial h_w}{\partial y}$ 为y方向的水力梯度。

降雨入渗时部分岩土体由非饱和状态转变为饱和，坡体内孔隙水压力逐渐增加，非饱和区孔隙水的运动和饱和区孔隙水的运动是相互关联的，将两者统一起来就是所谓的饱和–非饱和渗流问题。

采用水头h作为控制方程的因变量，对于各向异性的二维饱和与非饱和渗流控制方程为：

$\frac{\partial }{\partial x}\left(k_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_y\frac{\partial h}{\partial y}\right)=m_w{\rho }_{w}g\frac{\partial h}{\partial t}$ (4-2)

式中： $k_x$ 、 $k_y$ 分别为水平和垂直方向的饱和渗透系数； ${\rho }_w$ 为水的密度； $g$ 为重力加速度； $m_w$ 为比水容重，定义为体积含水量 ${\theta }_w$ 对基质吸力 $\left(u_a-u_w\right)$ 偏导数的负值：

$m_w=-\frac{\partial {\theta }_w}{\partial \left(u_a-u_w\right)}$ (4-3)

降雨引起的滑坡体内暂态渗流场浸润线的变化，其边界条件主要包含定水头边界及定流量边界两类：

水头边界：

${k\frac{\partial h}{\partial n}|}_{{\Gamma }_1}=h\left(x,y,t\right)$ (4-4)

流量边界：

${k\frac{\partial h}{\partial n}|}_{{\Gamma }_2}=q\left(x,y,t\right)$ (4-5)

k为渗透系数张量，n为边界面单位法向量。

当水流过土体中的时候，一部分水要驻留在土体结构中，驻留的水量是孔隙水压力和土体结构特征的函数。对渗流分析来说，定义水量驻留部分的体积和总体积的比值为体积水含量，用公式表示为：

$\theta =\frac{V_W}{V}$ (4-6)

式中： $V_W$ ―土体单元中水的体积；

V―土体的总体积。

当饱和度是100%的时候，体积水含量等于土壤的孔隙率。体积水含量的改变依赖于应力状态的改变和土体的性质。在渗流计算中，假定总应力是不变的，且孔隙气压力也保持不变，这意味着体积水含量的改变仅依赖于孔隙水压力的改变，用公式可表示为：

$\partial \theta =m_w\partial u_w$ (4-7)

式中： $u_w$ —为孔隙水压力。

由于渗透系数(水力传导率)是表示土体导水能力的一个参数，水力传导系数依赖于体积水含量(含水率)，而含水率又是孔隙水压力的函数，渗透系数是含水率的函数，因此渗透系数是孔隙水压力的间接函数。

本文采用Geo-Studio系列软件中的SEEP/W模块对该滑坡进行渗流分析。SEEP/W模块具有饱和−非饱和渗流稳定渗流计算、饱和−非饱和非稳定渗流计算等功能。不同降雨强度下的滑坡稳定性定量计算基于Geo-studio有限元分析软件进行地下水位数值模拟并采用极限平衡法计算滑坡稳定性。首先采用SEEP/W模块计算出不同工况下的浸润线；然后，将渗流计算获得的地下水位变动曲线与传递系数法相耦合，计算滑坡的稳定性系数。

4.2. 数值建模及参数选取

4.2.1. 模型的建立

分析滑坡的物质组成、地层结构、变形机制和主要影响因素，把握其变形破坏的基本规律和主控因素的基础上，建立拜殿乡船山村滑坡合理的地质模型，选择滑坡的主滑剖面作为该滑坡的计算模型。

根据船山村滑坡的地质模型选取主剖面2-2’剖面进行地质建模。根据上述勘查资料将滑坡从上到下分为碎石土(滑体)、全风化花岗岩(滑体)、全风化花岗岩(滑带)、强–中风化花岗岩(基岩) 4个部分，计算模型见图4所示。

Figure 4. Calculation model diagram

图4. 计算模型图

计算剖面采用四边形或三角形对全局进行网格划分，全局单元的大概尺寸为1 m，计算剖面共划分2688个节点，2499个单元。

4.2.2. 计算工况

基于船山村滑坡自2017年12月28日至2018年12月28日一个水文年的降雨监测数据，根据国家气象局颁布的降雨强度划分标准，在一个水文年内该区域不同日降雨强度I的发生概率如表2所示。

Table 2. Rainfall intensity and occurrence probability

表2. 降雨强度与发生概率

根据表2所示，船山村滑坡在未来很有可能在24 h内遭受暴雨、大雨、中雨以及小雨的作用。根据一个水文年内的降雨监测数据，船山村滑坡连续降雨天数分别为2、3、4、5、6、8、12天，连续降雨最长为12天，累计降雨量达195 mm。此次计算工况结合降雨监测数据及考虑降雨强度和降雨时长综合确定，最终选择不同的日降雨强度为5 mm/d、10 mm/d、20 mm/d、30 mm/d、40 mm/d、50 mm/d及60 mm/d，历时分别为：3天、5天、7天、10天以及12天，累计降雨量从15 mm到480 mm作为计算工况。

4.3. 船山村滑坡渗流场变化特征

4.3.1. 不同降雨过程中滑坡渗流场的动态响应特征分析

本文采用Geo-studio中SEEP/W模块进行地下水渗流研究，分别对降雨3 d、5 d、7 d、10 d、12 d典型工况进行地下水浸润线模拟，结果如图5~9所示。

Figure 5. Variation of groundwater phreatic line under 10 mm rainfall intensity for 3 consecutive days

图5. 连续3天10 mm降雨强度工况地下水浸润线变化

Figure 6. Variation of groundwater phreatic line under 40 mm rainfall intensity for 5 consecutive days

图6. 连续5天40 mm降雨强度工况地下水浸润线变化

Figure 7. Variation of groundwater phreatic line under 30 mm rainfall intensity for 7 consecutive days

图7. 连续7天30 mm降雨强度工况地下水浸润线变化

Figure 8. Variation of groundwater phreatic line under 30 mm rainfall intensity for 10 consecutive days

图8. 连续10天30 mm降雨强度工况地下水浸润线变化

Figure 9. Variation of groundwater phreatic line under 30 mm rainfall intensity for 12 consecutive days

图9. 连续12天30 mm降雨强度工况地下水浸润线变化

在Geostudio中SEEP/W模块进行地下水渗流模拟计算，数值模拟结果表明：开始降雨后，坡体内地下水浸润线出现明显抬升，抬升幅度较大，降雨结束后随着地下水排出坡体，浸润线开始缓慢下降。计算结果见表3所示，分析表明，地下水浸润线的抬升，同时受到降雨时长和日降雨量的影响，降雨时长的增加或者累计降雨量的增加，地下水位都会出现明显的抬升。

Table 3. Water level rise under calculation conditions

表3. 计算工况水位线上升值

4.3.2. 不同渗流场特征下的滑坡稳定性分析

将各工况渗流场的模拟结果导入到Geostudio中Slope/W模块进行相应工况下的滑坡稳定性计算。稳定性计算剖面继承用于进行地下水渗流模拟计算所建立的模型，采用Morgenstern-Price法，在稳定性分析中，既能考虑滑坡体应力变化对稳定性的影响，又能真实的反映危险滑动面的分布情况，还能以安全系数的形式来表征坡体稳定性的情况。

各典型工况稳定性系数随时间变化曲线如图10~14，由稳定性系数随时间变化情况可得：在未降雨条件下，随着地下水排出坡体，滑坡稳定性系数缓慢上升，开始降雨之后，地下水位开始抬升，稳定性系数随着地下水位的抬升开始下降，降雨结束之后，坡体内的地下水开始排出，稳定性系数随着地下水位的下降呈现出上升趋势。根据上述分析说明拜殿乡船山村滑坡稳定性对降雨响应明显，降雨对拜殿乡船山村滑坡的稳定性影响较大。

Figure 10. Response law of 3 D rainfall landslide stability to seepage field

图10. 3 d降雨滑坡稳定性对渗流场的响应规律

Figure 11. Response law of 5 D rainfall landslide stability to seepage field

图11. 5 d降雨滑坡稳定性对渗流场的响应规律

Figure 12. Response law of 7 D rainfall landslide stability to seepage field

图12. 7 d降雨滑坡稳定性对渗流场的响应规律

Figure 13. Response law of 10 D rainfall landslide stability to seepage field

图13. 10 d降雨滑坡稳定性对渗流场的响应规律

Figure 14. Response law of 12 D rainfall landslide stability to seepage field

图14. 12 d降雨滑坡稳定性对渗流场的响应规律

5. 滑坡预警模型建立

本项目选取1年的监测周期，分多次降雨过程，结合数理统计和数值模拟计算的手段，从降雨量和稳定性系数两方面构建单指标的预警判据。结合地质灾害四级预警机制，建立滑坡危险性预警判据，确定不同工况下的滑坡危险性预警等级。

5.1. 降雨条件下滑坡稳定性响应规律

对不同降雨条件下船山村滑坡地下水渗流模拟以及稳定性数值计算的结果进行分析(详见图15)，发现随着累计降雨量的增加，滑坡的稳定性系数不断减小。同时也出现相同累计降雨量条件下，降雨时长不同滑坡的稳定性不同的现象，如累计降雨量为150 mm时，降雨历时3天，日降雨量为50 mm/d滑坡最小稳定性系数为1.137；历时5天，日降雨量为30 mm/d滑坡最小稳定性系数为1.125。两工况间最小稳定系数相差0.012，在日降雨量小而降雨历时长的工况条件下，滑坡最小稳定性系数小。

不同降雨时长和累计降雨量会引起滑坡稳定性的变化不同。当连续降雨3天时，累计降雨量从15变化到180 mm，滑坡稳定性系数的变化幅度较小为0.026；连续降雨5天时，累计降雨量从25变化到300 mm，稳定性系数降幅为0.058；连续降雨7天时，累计降雨量从35变化到420 mm，稳定性系数降幅为0.112；连续降雨10天时，累计降雨量从50变化到400 mm稳定性系数降幅为0.153；连续降雨12天时，累计降雨量从60变化到480 mm，稳定性系数变化幅度较大为0.186。

上述分析说明拜殿乡船山村滑坡的稳定性不仅受累计降雨量的影响，滑坡稳定性对降雨时长的响应也非常明显。因此，该滑坡的预警模型的建立应综合考虑累计降雨量和降雨时长的影响。

Figure 15. Variation characteristics of landslide stability coefficient under different rainfall conditions

图15. 不同降雨条件下滑坡稳定性系数变化特征

5.2. 滑坡预警模型

本文通过分析滑坡专业监测数据和上述计算的结果，最终确定船山村滑坡主要从降雨量累计降雨量以及降雨时长两方面建立滑坡预警的指标体系。分别对相同降雨历时不同日降雨量工况下滑坡稳定性计算的结果进行非线性曲线拟合，根据拟合结果与稳定性系数和累计降雨量的对应关系，分别建立不同降雨历时的滑坡预警模型(图16~19)。

比如降雨历时3 d时，当累计降雨量在0~50 mm之间，滑坡处于基本稳定状态；在50 mm~180 mm之间，滑坡最小稳定性系数小于1.15，滑坡预警级别进入注意级。降雨历时12 d时，当累计降雨量小于200 mm时，滑坡预警级别为注意级；累计降雨量在200~270 mm之间时，滑坡最小稳定性系数在1.05~1.10之间，为警示级；累计降雨量在270~365 mm之间时，最小稳定性系数在1.05~1.10之间，为警戒级；累计降雨量在365~500 mm之间时，最小稳定性系数在1.05~1.10之间，为警报级。详见表4。

Table 4. Landslide early warning level and index system under different rainfall conditions

表4. 不同降雨工况下滑坡预警等级及指标体系

Figure 16. 5 D rainfall landslide early warning model

图16. 5 D降雨滑坡预警模型

Figure 17. 7 D rainfall landslide early warning model

图17. 7 D降雨滑坡预警模型

Figure 18. 10 D rainfall landslide early warning model

图18. 10 D降雨滑坡预警模型

Figure 19. 12 D rainfall landslide early warning model

图19. 12 D降雨滑坡预警模型

5.3. 船山滑坡监测平台预警综合模型工作流程

根据上文的Geo-studio数值模拟结果，判断滑坡发生预警等级，标准见表5，由此作为该监测点的预警级别。

Table 5. Early warning level of landslide

表5. 滑坡发生预警等级

6. 对今后工作的建议

在未来模型的使用过程中应注重参数及统计方法的修正，根据监测点在湖南省不同滑坡地带分区、不同地形条件、不同地质环境进一步完善降雨诱发机制，注重滑坡中期及活动期监测，以提高模型预测的科学程度和精确度。

在系统平台的应用中，一方面需要改善平台数据的实时更新情况，提高基础数据质量；另一方面应注意软件的与系统平台的对接，不断完善算法，实现软件的稳定、准确运算。