OGS軟件在地質(zhì)災(zāi)害監(jiān)測預(yù)警系統(tǒng)中的應(yīng)用

丁靖洋

（北京市地質(zhì)工程勘察院，北京 100048）

OGS軟件在地質(zhì)災(zāi)害監(jiān)測預(yù)警系統(tǒng)中的應(yīng)用

丁靖洋

（北京市地質(zhì)工程勘察院，北京 100048）

地質(zhì)災(zāi)害對交通網(wǎng)絡(luò)線性工程的安全性構(gòu)成重大威脅。本文依托于“京張鐵路（居庸關(guān)隧道進口段）地質(zhì)安全監(jiān)測示范工程”項目，以一處堆渣體為研究對象，基于有效應(yīng)力原理、Mohr-Coulomb屈服準則、Richard-flow非飽和滲流方程，采用OGS軟件對降雨入滲條件下的流固耦合過程進行數(shù)值模擬。研究結(jié)果表明，隨著降雨強度的不斷增加，堆渣體上部土體的有效飽和度逐漸增大，地災(zāi)發(fā)生區(qū)域開始出現(xiàn)并相互貫通，這將為地質(zhì)安全預(yù)警提供依據(jù)。

地質(zhì)災(zāi)害預(yù)警；降雨入滲；流固耦合；OGS數(shù)值模擬

0 前言

北京地區(qū)以中低山地貌為主，其中山區(qū)面積約占全市面積的62%，且溝谷分布廣泛，地形相對高差較大，巖體較為破碎，地質(zhì)構(gòu)造復(fù)雜，新構(gòu)造活動強烈；北京降水多集中在7—8月，且時常有暴雨、大暴雨發(fā)生；加之隨著經(jīng)濟社會的發(fā)展，人類活動對地質(zhì)體環(huán)境的影響日趨明顯。因此北京山區(qū)內(nèi)的地質(zhì)災(zāi)害，尤其泥石流、滑坡、崩塌、采空塌陷等較為嚴重。據(jù)統(tǒng)計，北京地區(qū)包含泥石流591條、滑坡15個、崩滑塌34550個、采空塌陷162個（董桂芝等，1997）。

作為“京津冀一體化”交通網(wǎng)絡(luò)建設(shè)的重要組成部分，京張城際鐵路的地質(zhì)安全性尤其引人注意。由于交通網(wǎng)絡(luò)線性工程的特點，地質(zhì)災(zāi)害的發(fā)生將對其造成巨大損失。其中京張城際鐵路（北京段）穿越軍都山，沿途存在多種不良地質(zhì)災(zāi)害，為了確保京張城際鐵路在建設(shè)、運行中的安全，對沿線周邊不良地質(zhì)災(zāi)害體進行實時監(jiān)測并及時預(yù)警具有重要意義。

京張鐵路（居庸關(guān)隧道進口段）地質(zhì)安全監(jiān)測示范工程，是京津冀協(xié)同發(fā)展交通網(wǎng)絡(luò)地質(zhì)安全監(jiān)測預(yù)警系統(tǒng)建設(shè)中的試驗研究項目。本文在該示范工程試驗區(qū)選址的基礎(chǔ)上，以一處堆渣體為分析對象，通過對三維地質(zhì)數(shù)據(jù)的處理建立了有限元計算模型，基于有效應(yīng)力原理、非飽和滲透原理和Mohr-Coulomb塑性理論，對堆渣體位移與降雨量關(guān)系、滲流-應(yīng)力耦合過程和潛在危險區(qū)域進行了研究。

1 工程背景

試驗區(qū)位于北京市昌平區(qū)姚店一處采石場，面積約0.5km2，地面高程約285m，山體自然坡度一般為25°～59°，區(qū)內(nèi)整體地形呈兩側(cè)高，中間低，地形起伏，橫向“V”型沖溝發(fā)育，均匯入中部低洼處的縱向河谷中，局部成陡崖巖墻狀。

通過野外調(diào)查發(fā)現(xiàn)，試驗區(qū)內(nèi)發(fā)育有一條泥石流溝，溝谷長度較長，規(guī)模較大，為V型溝。溝內(nèi)主要物源為采石場開采過程中產(chǎn)生的廢石、廢渣和大量松散坡積物，稱之為“堆渣體”。

為了便于地質(zhì)災(zāi)害的現(xiàn)場監(jiān)測預(yù)警，試驗區(qū)內(nèi)布設(shè)有GNSS、靜力水準、SAA、測斜儀、土壓力計、地下水動態(tài)自動監(jiān)測儀等設(shè)備，24小時自動記錄數(shù)據(jù)，監(jiān)測內(nèi)容包括降雨量、水平/垂直位移、土壓力、土壤含水率和地下水滲透率等。

2 三維建模工作

2.1 三維地質(zhì)體展示模型

借助于Creater建模軟件平臺，在對堆渣體鉆孔數(shù)據(jù)進行篩選及處理的基礎(chǔ)上，開展三維地質(zhì)體展示模型的構(gòu)建工作。堆渣體建模分為3個地層、10個剖面，其中從上到下地層分別為坡積物層、坡積碎石及碎石層和基巖層。

2.2 三維地質(zhì)體模型網(wǎng)格剖分

建模分析工作以O(shè)penGeoSys (OGS) 數(shù)值計算軟件為平臺：OpenGeoSys是針對孔隙介質(zhì)中熱-水-力-化學（THMC）多場耦合過程的國際著名開放性數(shù)值模擬平臺，其提供地質(zhì)和水力方面的多場耦合有限元計算框架平臺，軟件程序允許C++語言進行嵌入編程。

基于三維展示模型獲得的堆渣體高程數(shù)據(jù)，導(dǎo)出了數(shù)據(jù)散點形式的不規(guī)則三角形文件(TINs)格式。通過自行編寫的C++語言程序，讀取TINs格式文件數(shù)據(jù)，從而建立了GMSH網(wǎng)格文件（.Geo）。為了減少數(shù)值計算工作量，模型減少了底層基巖的厚度。

通過GMSH軟件自動劃分了四面體網(wǎng)格，基于Frontal算法，有限單元網(wǎng)格劃分如圖1所示，數(shù)值計算單元共24498個。

圖1 數(shù)值模擬有限單元網(wǎng)格模型Fig.1 Finite Element Mesh of the model

3 非飽和滲流-應(yīng)力耦合理論

3.1 有效應(yīng)力原理

降雨誘發(fā)地質(zhì)災(zāi)害的過程本質(zhì)上為非飽和孔隙介質(zhì)入滲的過程，即滲流—應(yīng)力耦合過程。在描述非飽和孔隙介質(zhì)入滲過程的應(yīng)力分量理論中，最為廣泛應(yīng)用的理論為有效應(yīng)力原理（Terzaghi，1923）。

在降雨誘發(fā)地質(zhì)災(zāi)害的過程中，堆渣體中的孔隙壓力隨降雨入滲而不斷變化。Bolzon 等（1996）通過熱力學理論對非飽和土中的能量方程進行了驗證，從而得到了非飽和有效應(yīng)力公式為：

式中：S為孔隙吸力，uw為孔隙水壓力，uα為孔隙氣壓力，Se為有效飽和度，δij為克羅內(nèi)克（Kroenker delta）符號。

本次研究工作中，堆渣體坡積物以碎石土為主，孔隙率較大，故忽略孔隙氣壓力uα，只考慮孔隙水壓力uw對滲流-應(yīng)力耦合過程的影響。因此，有效應(yīng)力公式為

需要說明的是，Bolzon推導(dǎo)的有效應(yīng)力公式引入了有效飽和度Se這一參數(shù)，表征了降雨入滲對非飽和孔隙介質(zhì)的影響。當飽和度為1時，Bolzon有效應(yīng)力公式可退化為太沙基公式。

孔隙介質(zhì)中的應(yīng)力應(yīng)變理論基于Hook定律。在OpenGeoSys中，采用返回映射算法計算彈塑性應(yīng)力應(yīng)變增量。

3.2 考慮入滲影響的Mohr-Coulomb理論

Mohr-coulomb屈服準則廣泛應(yīng)用于描述巖土體的塑性力學行為，其表達式為

式中：τ為剪切應(yīng)力，σb為主應(yīng)力，φ和c分別為內(nèi)摩擦角和粘聚力。

若考慮降雨入滲對巖土體屈服準則的影響，則引入有效應(yīng)力作為Mohr-Coulomb準則的主應(yīng)力，將式（2）代入式（3），可得

3.3 非飽和滲流理論

若忽略物質(zhì)溶解，則孔隙介質(zhì)中的流體流動應(yīng)遵循質(zhì)量守恒定律。以達西尺度的代表性體積單元（RVE）為例，RVE內(nèi)流體質(zhì)量的變化量與通過RVE邊界的流體質(zhì)量應(yīng)相等。對不可變形孔隙介質(zhì)，在無源、無匯的情況下，流體質(zhì)量守恒關(guān)系表達式為

式中：下標w代表孔隙介質(zhì)中的某流體相，ρw為孔隙介質(zhì)中流體的密度，Se表示孔隙介質(zhì)中流體相的有效飽和度，n為孔隙介質(zhì)孔隙度，vw為流體流過RVE的流速，nu為在RVE邊界上沿流動方向的單位向量，?U為RVE的邊界。

Richard方程（Richards，1931）是描述非飽和孔隙介質(zhì)流動的經(jīng)典方程。對式（5）進行積分，由鏈式微分法則，引入達西方程，考慮孔隙介質(zhì)形變對流體質(zhì)量平衡的影響，則Richard滲流本構(gòu)方程為：

式中：孔隙吸力 用毛細管壓力等效，等于負孔隙水壓力，即uw=-S，k為孔隙介質(zhì)固有滲透率，krel為非飽和到飽和過程中的相對滲透率，μw為流體粘性系數(shù)。毛細管壓力與有效飽和度的關(guān)系（通常稱為土水特征曲線，Water Retention Curve, WRC）將通過Van Genuchten模型確定（Genuchten，1980），即

式中：a和m為相關(guān)力學參數(shù)。

相對滲透率公式采用冪函數(shù)形式如：

4 數(shù)值模擬參數(shù)

4.1 力學參數(shù)

堆渣體自上而下分別為坡積物層、坡積碎石及碎石層和基巖層，根據(jù)現(xiàn)場勘察資料、室內(nèi)土工試驗及參考《工程地質(zhì)手冊》，并借助于國內(nèi)外學者的相關(guān)研究成果（王新剛等，2013；胡凱衡等，2014；黃龍波等，2016；陳善雄等，2001；徐晗等，2005），確定數(shù)值模擬所需的力學參數(shù)（表1）。

4.2 滲流參數(shù)

對于滲流參數(shù)，堆渣體非飽和到飽和過程中的相對滲透率采用公式（8）計算。數(shù)值模擬中的滲流參數(shù)見表2。

本次研究工作中的土水特征曲線（WRC）引用了文獻（Cho，2016）中的數(shù)據(jù)，通過提取數(shù)據(jù)，得到了土水特征曲線如圖2(a)，相對滲透率計算結(jié)果如圖2(b)所示。

表1 數(shù)值模擬力學參數(shù)Tab.1 Mechanics parameters for the numerical modeling

表2 數(shù)值模擬滲流參數(shù)Tab.2 Seepage parameters for the numerical modeling

圖2 數(shù)值模擬中（a）土水特征曲線；（b）相對滲透率曲線Fig.2 (a)water retention curve; (b)relative permeability curve for the numerical modeling

5 降雨誘發(fā)地質(zhì)災(zāi)害的數(shù)值模擬研究

5.1 初值及邊界條件

參考相關(guān)研究成果（王新剛等，2013；胡凱衡等，2014；黃龍波等，2016；陳善雄等，2001；徐晗等，2005；Cho，2016），本次數(shù)值模擬中，對三層土分別賦予初值，第一層初始孔隙吸力為48000Pa, 表明第一層土的初始有效飽和度約為0.09；第二層初始孔隙吸力為9000Pa，表明第二層土的初始有效飽和度為0.23；第三層初始孔隙吸力為110000Pa，代表基巖層初始有效飽和度為0.01。

設(shè)置環(huán)繞模型的曲面在X和Y方向不可移動，即固定曲面在X和Y方向的位移。

為研究降雨強度與堆渣體位移的關(guān)系，將降雨強度轉(zhuǎn)化為模型中的透水量，并假設(shè)第一層土的上表面持續(xù)降雨入滲，以此為關(guān)系變量進行數(shù)值模擬計算。

5.2 數(shù)值計算結(jié)果分析

本研究工作基于單元體積積分計算模型中的透水量，通過定義計算時間步長，得到模型表層的入滲流速。通過GMESH軟件得到模型表層面積為3722.32m2，進一步計算得到降雨強度（mm）。

5.2.1 降雨過程中的流固耦合

降雨入滲過程中的滲流云圖如3所示，圖3(a)為降雨強度達到100mm時的有效飽和度分布云圖。在降雨入滲過程中，堆渣體左右和上方邊界處的有效飽和度逐漸增大，山體中部的有效飽和度增量較小。圖3(b)顯示了第二層土的孔隙水壓力分布云圖和第一層土中的流速分布云圖。由圖3(b)可知，在降雨277mm時，第一層土中山體左右和上方邊界處的流速較大，且由四周匯向中部。

由圖3中可知，由于堆渣體左右和上部區(qū)域土層較薄，中部區(qū)域土層較厚，由滲透壓力梯度可知，土層厚度越薄，滲透壓力梯度越大，由圖3(b)中可知，降雨277 mm 時周邊區(qū)域的滲透驅(qū)動力為孔隙水重力。在較大的滲透壓力作用下，孔隙水逐漸流向中部。

綜上可知，在孔隙水壓力的作用下，堆渣體左右和上部邊界區(qū)域?qū)l(fā)生較大位移。

圖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 降雨過程中的應(yīng)變

降雨入滲過程中堆渣體的應(yīng)變云圖如4所示。圖4(a)為降雨100mm時堆渣體應(yīng)變云圖，第一層土的左右和上部坡腳處體應(yīng)變最大。當降雨達到277 mm時，如圖4(b)所示，山體左右和上部坡面的應(yīng)變增大。

綜上所述，隨著降雨強度的增大，堆渣體左右坡腳和上部坡面處的體應(yīng)變不斷增加，并逐漸相互貫通匯成一個較大的滑動帶，且圖3、圖4相互印證，表明降雨的持續(xù)入滲是導(dǎo)致堆渣體產(chǎn)生滑移的最重要原因，這有助于為地質(zhì)災(zāi)害預(yù)警提供依據(jù)。

5.2.3 降雨過程中的危險區(qū)域判別

根據(jù)文獻（黃龍波等，2016；Cho，2016），基于最大剪應(yīng)力準則，計算堆渣體的剪切強度：

式中：ccritical、φcritical和fcritical分別為與極限剪應(yīng)力強度相關(guān)的粘聚力、內(nèi)摩擦角和安全系數(shù)，τcritical和σs為剪切強度和壓縮強度。取安全系數(shù)為1.3，根據(jù)表2計算得到第一層土的剪切強度為8.3kPa，第二層土的剪切強度為30.4kPa。

圖4 降雨入滲過程中的應(yīng)變云圖(a)降雨100mm時的應(yīng)變云圖 (b)降雨277mm時的應(yīng)變云圖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 時的地質(zhì)災(zāi)害區(qū)域分布圖。圖5(a)中顯示，在降雨累積137mm時堆渣體出現(xiàn)了初始災(zāi)害，隨著累積降雨量的增加，災(zāi)害區(qū)域開始擴展并逐漸匯集成片，見圖5(b)和(c)。地質(zhì)災(zāi)害區(qū)域主要發(fā)生于陡坡區(qū)域，從坡腳向坡面延伸。

本文的數(shù)值模擬研究中，不僅考慮了土體孔隙壓力變化導(dǎo)致的有效應(yīng)力變化，還考慮了堆渣體自重影響下的有效應(yīng)力變化，圖5的模擬結(jié)果表明，基于Mohr-Coulomb屈服準則來判定地質(zhì)災(zāi)害危險區(qū)域具有可行性和合理性。

圖5 降雨入滲過程中的災(zāi)害分布區(qū)域(a)降雨137mm時的災(zāi)害區(qū)域 (b)降雨189mm時的災(zāi)害區(qū)域(c)降雨277mm時的災(zāi)害區(qū)域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 結(jié)論

本文在三維地質(zhì)體展示模型的基礎(chǔ)上進行了有限元網(wǎng)格剖分，并借助于OGS多物理場數(shù)值計算平臺，基于有效應(yīng)力原理、Mohr-Coulomb屈服準則、Richard-flow非飽和滲流方程，對堆渣體在降雨入滲過程中的滲流-應(yīng)力耦合過程進行了數(shù)值模擬研究。研究結(jié)果表明：堆渣體左右及上部土層較薄的區(qū)域在較大的滲透壓力梯度作用下易產(chǎn)生較大變形，且該處坡面較陡，雨水入滲過程中最易誘發(fā)地質(zhì)災(zāi)害，其中的主要驅(qū)動力有滲透壓力、孔隙水自重和堆渣體土層自重。隨著降雨量增大，有效飽和度增加，受孔隙壓力變化的影響，潛在危險區(qū)域逐漸貫通并將發(fā)生較大規(guī)模的地質(zhì)災(zāi)害。

需要說明的是，本文數(shù)值模擬中的參數(shù)多借助于工程經(jīng)驗、國內(nèi)外學者的研究成果等，下一步的工作重點是，在現(xiàn)場實際監(jiān)測數(shù)據(jù)的基礎(chǔ)上，合理反演各種力學、滲流參數(shù)，以獲得符合北京山區(qū)地質(zhì)特征的流固耦合參數(shù)，并進一步預(yù)測潛在的危險區(qū)域。

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. 考慮降雨的非飽和土邊坡穩(wěn)定性分析方法[J]. 巖土力學，22(4)：447-450.

董桂芝，紀玉杰，單青生，1997. 北京市地質(zhì)災(zāi)害現(xiàn)狀分析[J].北京地質(zhì)，(1)：30-33.

胡凱衡，馬超，2014. 泥石流啟動臨界土體含水量及其預(yù)警應(yīng)用[J]. 地球科學與環(huán)境學報，(2)：73-80.

黃龍波，汪洪星，談云志，等，2016. 降雨對巖土邊坡穩(wěn)定影響的實例分析[J]. 三峽大學學報（自然科學版），38(4)：40-45.

王新剛，胡斌，劉強，等，2013. 碎石土大型直剪研究與邊坡穩(wěn)定性分析[J]. 長江科學院院報，30(6)：63-67.

徐晗，朱以文，蔡元奇，等，2005. 降雨入滲條件下非飽和土邊坡穩(wěn)定分析[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

北京市優(yōu)秀人才資助項目（2016400617931G171）

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