OGS软件在地质灾害监测预警系统中的应用

丁靖洋

（北京市地质工程勘察院，北京 100048）

OGS软件在地质灾害监测预警系统中的应用

丁靖洋

（北京市地质工程勘察院，北京 100048）

地质灾害对交通网络线性工程的安全性构成重大威胁。本文依托于“京张铁路（居庸关隧道进口段）地质安全监测示范工程”项目，以一处堆渣体为研究对象，基于有效应力原理、Mohr-Coulomb屈服准则、Richard-flow非饱和渗流方程，采用OGS软件对降雨入渗条件下的流固耦合过程进行数值模拟。研究结果表明，随着降雨强度的不断增加，堆渣体上部土体的有效饱和度逐渐增大，地灾发生区域开始出现并相互贯通，这将为地质安全预警提供依据。

地质灾害预警；降雨入渗；流固耦合；OGS数值模拟

0 前言

北京地区以中低山地貌为主，其中山区面积约占全市面积的62%，且沟谷分布广泛，地形相对高差较大，岩体较为破碎，地质构造复杂，新构造活动强烈；北京降水多集中在7—8月，且时常有暴雨、大暴雨发生；加之随着经济社会的发展，人类活动对地质体环境的影响日趋明显。因此北京山区内的地质灾害，尤其泥石流、滑坡、崩塌、采空塌陷等较为严重。据统计，北京地区包含泥石流591条、滑坡15个、崩滑塌34550个、采空塌陷162个（董桂芝等，1997）。

作为“京津冀一体化”交通网络建设的重要组成部分，京张城际铁路的地质安全性尤其引人注意。由于交通网络线性工程的特点，地质灾害的发生将对其造成巨大损失。其中京张城际铁路（北京段）穿越军都山，沿途存在多种不良地质灾害，为了确保京张城际铁路在建设、运行中的安全，对沿线周边不良地质灾害体进行实时监测并及时预警具有重要意义。

京张铁路（居庸关隧道进口段）地质安全监测示范工程，是京津冀协同发展交通网络地质安全监测预警系统建设中的试验研究项目。本文在该示范工程试验区选址的基础上，以一处堆渣体为分析对象，通过对三维地质数据的处理建立了有限元计算模型，基于有效应力原理、非饱和渗透原理和Mohr-Coulomb塑性理论，对堆渣体位移与降雨量关系、渗流-应力耦合过程和潜在危险区域进行了研究。

1 工程背景

试验区位于北京市昌平区姚店一处采石场，面积约0.5km2，地面高程约285m，山体自然坡度一般为25°～59°，区内整体地形呈两侧高，中间低，地形起伏，横向“V”型冲沟发育，均汇入中部低洼处的纵向河谷中，局部成陡崖岩墙状。

通过野外调查发现，试验区内发育有一条泥石流沟，沟谷长度较长，规模较大，为V型沟。沟内主要物源为采石场开采过程中产生的废石、废渣和大量松散坡积物，称之为“堆渣体”。

为了便于地质灾害的现场监测预警，试验区内布设有GNSS、静力水准、SAA、测斜仪、土压力计、地下水动态自动监测仪等设备，24小时自动记录数据，监测内容包括降雨量、水平/垂直位移、土压力、土壤含水率和地下水渗透率等。

2 三维建模工作

2.1 三维地质体展示模型

借助于Creater建模软件平台，在对堆渣体钻孔数据进行筛选及处理的基础上，开展三维地质体展示模型的构建工作。堆渣体建模分为3个地层、10个剖面，其中从上到下地层分别为坡积物层、坡积碎石及碎石层和基岩层。

2.2 三维地质体模型网格剖分

建模分析工作以OpenGeoSys (OGS) 数值计算软件为平台：OpenGeoSys是针对孔隙介质中热-水-力-化学（THMC）多场耦合过程的国际著名开放性数值模拟平台，其提供地质和水力方面的多场耦合有限元计算框架平台，软件程序允许C++语言进行嵌入编程。

基于三维展示模型获得的堆渣体高程数据，导出了数据散点形式的不规则三角形文件(TINs)格式。通过自行编写的C++语言程序，读取TINs格式文件数据，从而建立了GMSH网格文件（.Geo）。为了减少数值计算工作量，模型减少了底层基岩的厚度。

通过GMSH软件自动划分了四面体网格，基于Frontal算法，有限单元网格划分如图1所示，数值计算单元共24498个。

图1 数值模拟有限单元网格模型Fig.1 Finite Element Mesh of the model

3 非饱和渗流-应力耦合理论

3.1 有效应力原理

降雨诱发地质灾害的过程本质上为非饱和孔隙介质入渗的过程，即渗流—应力耦合过程。在描述非饱和孔隙介质入渗过程的应力分量理论中，最为广泛应用的理论为有效应力原理（Terzaghi，1923）。

在降雨诱发地质灾害的过程中，堆渣体中的孔隙压力随降雨入渗而不断变化。Bolzon 等（1996）通过热力学理论对非饱和土中的能量方程进行了验证，从而得到了非饱和有效应力公式为：

式中：S为孔隙吸力，uw为孔隙水压力，uα为孔隙气压力，Se为有效饱和度，δij为克罗内克（Kroenker delta）符号。

本次研究工作中，堆渣体坡积物以碎石土为主，孔隙率较大，故忽略孔隙气压力uα，只考虑孔隙水压力uw对渗流-应力耦合过程的影响。因此，有效应力公式为

需要说明的是，Bolzon推导的有效应力公式引入了有效饱和度Se这一参数，表征了降雨入渗对非饱和孔隙介质的影响。当饱和度为1时，Bolzon有效应力公式可退化为太沙基公式。

孔隙介质中的应力应变理论基于Hook定律。在OpenGeoSys中，采用返回映射算法计算弹塑性应力应变增量。

3.2 考虑入渗影响的Mohr-Coulomb理论

Mohr-coulomb屈服准则广泛应用于描述岩土体的塑性力学行为，其表达式为

式中：τ为剪切应力，σb为主应力，φ和c分别为内摩擦角和粘聚力。

若考虑降雨入渗对岩土体屈服准则的影响，则引入有效应力作为Mohr-Coulomb准则的主应力，将式（2）代入式（3），可得

3.3 非饱和渗流理论

若忽略物质溶解，则孔隙介质中的流体流动应遵循质量守恒定律。以达西尺度的代表性体积单元（RVE）为例，RVE内流体质量的变化量与通过RVE边界的流体质量应相等。对不可变形孔隙介质，在无源、无汇的情况下，流体质量守恒关系表达式为

式中：下标w代表孔隙介质中的某流体相，ρw为孔隙介质中流体的密度，Se表示孔隙介质中流体相的有效饱和度，n为孔隙介质孔隙度，vw为流体流过RVE的流速，nu为在RVE边界上沿流动方向的单位向量，∂U为RVE的边界。

Richard方程（Richards，1931）是描述非饱和孔隙介质流动的经典方程。对式（5）进行积分，由链式微分法则，引入达西方程，考虑孔隙介质形变对流体质量平衡的影响，则Richard渗流本构方程为：

式中：孔隙吸力 用毛细管压力等效，等于负孔隙水压力，即uw=-S，k为孔隙介质固有渗透率，krel为非饱和到饱和过程中的相对渗透率，μw为流体粘性系数。毛细管压力与有效饱和度的关系（通常称为土水特征曲线，Water Retention Curve, WRC）将通过Van Genuchten模型确定（Genuchten，1980），即

式中：a和m为相关力学参数。

相对渗透率公式采用幂函数形式如：

4 数值模拟参数

4.1 力学参数

堆渣体自上而下分别为坡积物层、坡积碎石及碎石层和基岩层，根据现场勘察资料、室内土工试验及参考《工程地质手册》，并借助于国内外学者的相关研究成果（王新刚等，2013；胡凯衡等，2014；黄龙波等，2016；陈善雄等，2001；徐晗等，2005），确定数值模拟所需的力学参数（表1）。

4.2 渗流参数

对于渗流参数，堆渣体非饱和到饱和过程中的相对渗透率采用公式（8）计算。数值模拟中的渗流参数见表2。

本次研究工作中的土水特征曲线（WRC）引用了文献（Cho，2016）中的数据，通过提取数据，得到了土水特征曲线如图2(a)，相对渗透率计算结果如图2(b)所示。

表1 数值模拟力学参数Tab.1 Mechanics parameters for the numerical modeling

表2 数值模拟渗流参数Tab.2 Seepage parameters for the numerical modeling

图2 数值模拟中（a）土水特征曲线；（b）相对渗透率曲线Fig.2 (a)water retention curve; (b)relative permeability curve for the numerical modeling

5 降雨诱发地质灾害的数值模拟研究

5.1 初值及边界条件

参考相关研究成果（王新刚等，2013；胡凯衡等，2014；黄龙波等，2016；陈善雄等，2001；徐晗等，2005；Cho，2016），本次数值模拟中，对三层土分别赋予初值，第一层初始孔隙吸力为48000Pa, 表明第一层土的初始有效饱和度约为0.09；第二层初始孔隙吸力为9000Pa，表明第二层土的初始有效饱和度为0.23；第三层初始孔隙吸力为110000Pa，代表基岩层初始有效饱和度为0.01。

设置环绕模型的曲面在X和Y方向不可移动，即固定曲面在X和Y方向的位移。

为研究降雨强度与堆渣体位移的关系，将降雨强度转化为模型中的透水量，并假设第一层土的上表面持续降雨入渗，以此为关系变量进行数值模拟计算。

5.2 数值计算结果分析

本研究工作基于单元体积积分计算模型中的透水量，通过定义计算时间步长，得到模型表层的入渗流速。通过GMESH软件得到模型表层面积为3722.32m2，进一步计算得到降雨强度（mm）。

5.2.1 降雨过程中的流固耦合

降雨入渗过程中的渗流云图如3所示，图3(a)为降雨强度达到100mm时的有效饱和度分布云图。在降雨入渗过程中，堆渣体左右和上方边界处的有效饱和度逐渐增大，山体中部的有效饱和度增量较小。图3(b)显示了第二层土的孔隙水压力分布云图和第一层土中的流速分布云图。由图3(b)可知，在降雨277mm时，第一层土中山体左右和上方边界处的流速较大，且由四周汇向中部。

由图3中可知，由于堆渣体左右和上部区域土层较薄，中部区域土层较厚，由渗透压力梯度可知，土层厚度越薄，渗透压力梯度越大，由图3(b)中可知，降雨277 mm 时周边区域的渗透驱动力为孔隙水重力。在较大的渗透压力作用下，孔隙水逐渐流向中部。

综上可知，在孔隙水压力的作用下，堆渣体左右和上部边界区域将发生较大位移。

图3 降雨入渗过程中的渗流云图(a)降雨100mm时有效饱和度云图 (b)降雨277mm时流速云图Fig.3 Seepage diagrams during rainfall inf i ltration(a) Diagram of the effective saturation when rainfall is 100mm(b) Diagram of the fl ow rate when rainfall is 277mm

5.2.2 降雨过程中的应变

降雨入渗过程中堆渣体的应变云图如4所示。图4(a)为降雨100mm时堆渣体应变云图，第一层土的左右和上部坡脚处体应变最大。当降雨达到277 mm时，如图4(b)所示，山体左右和上部坡面的应变增大。

综上所述，随着降雨强度的增大，堆渣体左右坡脚和上部坡面处的体应变不断增加，并逐渐相互贯通汇成一个较大的滑动带，且图3、图4相互印证，表明降雨的持续入渗是导致堆渣体产生滑移的最重要原因，这有助于为地质灾害预警提供依据。

5.2.3 降雨过程中的危险区域判别

根据文献（黄龙波等，2016；Cho，2016），基于最大剪应力准则，计算堆渣体的剪切强度：

式中：ccritical、φcritical和fcritical分别为与极限剪应力强度相关的粘聚力、内摩擦角和安全系数，τcritical和σs为剪切强度和压缩强度。取安全系数为1.3，根据表2计算得到第一层土的剪切强度为8.3kPa，第二层土的剪切强度为30.4kPa。

图4 降雨入渗过程中的应变云图(a)降雨100mm时的应变云图 (b)降雨277mm时的应变云图Fig.4 Strain diagrams during rainfall inf i ltration(a) Strain diagram when rainfall is 100mm (b) Strain diagram when rainfall is 277mm

图5(a)-(c)显示了降雨强度在137mm、189mm和277 mm 时的地质灾害区域分布图。图5(a)中显示，在降雨累积137mm时堆渣体出现了初始灾害，随着累积降雨量的增加，灾害区域开始扩展并逐渐汇集成片，见图5(b)和(c)。地质灾害区域主要发生于陡坡区域，从坡脚向坡面延伸。

本文的数值模拟研究中，不仅考虑了土体孔隙压力变化导致的有效应力变化，还考虑了堆渣体自重影响下的有效应力变化，图5的模拟结果表明，基于Mohr-Coulomb屈服准则来判定地质灾害危险区域具有可行性和合理性。

图5 降雨入渗过程中的灾害分布区域(a)降雨137mm时的灾害区域 (b)降雨189mm时的灾害区域(c)降雨277mm时的灾害区域Fig.5 Geology disaster area during rainfall inf i ltration(a) Geology disaster area when rainfall is 137mm (b) Geology disaster area when rainfall is 189mm(c) Geology disaster area when rainfall is 277mm

6 结论

本文在三维地质体展示模型的基础上进行了有限元网格剖分，并借助于OGS多物理场数值计算平台，基于有效应力原理、Mohr-Coulomb屈服准则、Richard-flow非饱和渗流方程，对堆渣体在降雨入渗过程中的渗流-应力耦合过程进行了数值模拟研究。研究结果表明：堆渣体左右及上部土层较薄的区域在较大的渗透压力梯度作用下易产生较大变形，且该处坡面较陡，雨水入渗过程中最易诱发地质灾害，其中的主要驱动力有渗透压力、孔隙水自重和堆渣体土层自重。随着降雨量增大，有效饱和度增加，受孔隙压力变化的影响，潜在危险区域逐渐贯通并将发生较大规模的地质灾害。

需要说明的是，本文数值模拟中的参数多借助于工程经验、国内外学者的研究成果等，下一步的工作重点是，在现场实际监测数据的基础上，合理反演各种力学、渗流参数，以获得符合北京山区地质特征的流固耦合参数，并进一步预测潜在的危险区域。

Bolzon G, Schrefler B, Zienkiewicz O, 1996. Elastoplastic soil constitutive laws generalized to partially saturated states [J].Géotechnique, 46(2): 279-289.

Cho S E, 2016. Stability analysis of unsaturated soil slopes considering water-air flow caused by rainfall infiltration[J].Engineering Geology, 211: 184-197.

Genuchten M T V, 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of Americal Journal, 44(44):892-898.

Richards L A, 1931. Capillary conduction of liquids through porous mediums [J]. Journal of Applied Physics, 1(5): 318-333.

Terzaghi K V, 1923. Die berechnung der durchlassigkeitsziffer des tones aus dem verlauf der hydrodynamischen spannungserscheinungen [J]. Sitzungsberichte der Akademie der Wissenschaften in Wien, Mathematisch-Naturwissenschaftliche Klasse, Abteilung IIa, 132: 125-138.

陈善雄，陈守义，2001. 考虑降雨的非饱和土边坡稳定性分析方法[J]. 岩土力学，22(4)：447-450.

董桂芝，纪玉杰，单青生，1997. 北京市地质灾害现状分析[J].北京地质，(1)：30-33.

胡凯衡，马超，2014. 泥石流启动临界土体含水量及其预警应用[J]. 地球科学与环境学报，(2)：73-80.

黄龙波，汪洪星，谈云志，等，2016. 降雨对岩土边坡稳定影响的实例分析[J]. 三峡大学学报（自然科学版），38(4)：40-45.

王新刚，胡斌，刘强，等，2013. 碎石土大型直剪研究与边坡稳定性分析[J]. 长江科学院院报，30(6)：63-67.

徐晗，朱以文，蔡元奇，等，2005. 降雨入渗条件下非饱和土边坡稳定分析[J]. 岩土力学，26(12)：1957-1962.

The Application of OpenGeoSys Software in the Geological Disaster Monitoring and Forecast System

DING Jingyang

(Beijing Institute of Geological and Prospecting Engineering, Beijing 100048, China)

Geological disaster is of great treat for the linear project of transportation network. Based on the demonstration project of “Geological Safety Monitoring of Beijing-Zhangjiakou Railway (entrance of Juyong Pass tunnel )”, the fl uid-solid coupling process of a slag heap under the condition of rainfall inf i ltration is researched with the help of OpenGeoSys software, also the theories of effective stress, yield criterion of Mohr-Coulomb and Richard-f l ow equations. The results show that, for the slag heap, the effective saturation of the upper soil increases with the increase of rainfall intensity, also the area of geological disaster is appeared and joined. It provides a basis for the geological safety forecast.

Geological disaster forecast; Rainfall inf i ltration; Fluid-solid coupling; OGS numerical simulation

A

1007-1903（2017）03-0081-06

10.3969/j.issn.1007-1903.2017.03.0016

北京市优秀人才资助项目（2016400617931G171）

丁靖洋（1985- ），男，博士，工程师，主要从事岩土工程研究。E-mail：ak47djy@sina.com