丁靖洋

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

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

丁靖洋

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

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

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

0 前言

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

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

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

1 工程背景

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

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

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

2 三維建模工作

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

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

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

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

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

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

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

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

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

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

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

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

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

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

孔隙介質(zhì)中的應(yīng)力應(yīng)變理論基于Hook定律。在OpenGeoSys中，采用返回映射算法計(jì)算彈塑性應(yīng)力應(yīng)變?cè)隽俊?/p>

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

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

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

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

3.3 非飽和滲流理論

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

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

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

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

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

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

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

4.1 力學(xué)參數(shù)

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

4.2 滲流參數(shù)

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

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

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

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

圖2 數(shù)值模擬中（a）土水特征曲線(xiàn)；（b）相對(duì)滲透率曲線(xiàn)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ù)值模擬中，對(duì)三層土分別賦予初值，第一層初始孔隙吸力為48000Pa, 表明第一層土的初始有效飽和度約為0.09；第二層初始孔隙吸力為9000Pa，表明第二層土的初始有效飽和度為0.23；第三層初始孔隙吸力為110000Pa，代表基巖層初始有效飽和度為0.01。

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

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

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

本研究工作基于單元體積積分計(jì)算模型中的透水量，通過(guò)定義計(jì)算時(shí)間步長(zhǎng)，得到模型表層的入滲流速。通過(guò)GMESH軟件得到模型表層面積為3722.32m2，進(jìn)一步計(jì)算得到降雨強(qiáng)度（mm）。

5.2.1 降雨過(guò)程中的流固耦合

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

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

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

圖3 降雨入滲過(guò)程中的滲流云圖(a)降雨100mm時(shí)有效飽和度云圖 (b)降雨277mm時(shí)流速云圖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 降雨過(guò)程中的應(yīng)變

降雨入滲過(guò)程中堆渣體的應(yīng)變?cè)茍D如4所示。圖4(a)為降雨100mm時(shí)堆渣體應(yīng)變?cè)茍D，第一層土的左右和上部坡腳處體應(yīng)變最大。當(dāng)降雨達(dá)到277 mm時(shí)，如圖4(b)所示，山體左右和上部坡面的應(yīng)變?cè)龃蟆?/p>

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

5.2.3 降雨過(guò)程中的危險(xiǎn)區(qū)域判別

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

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

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

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

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

需要說(shuō)明的是，本文數(shù)值模擬中的參數(shù)多借助于工程經(jīng)驗(yàn)、國(guó)內(nèi)外學(xué)者的研究成果等，下一步的工作重點(diǎn)是，在現(xiàn)場(chǎng)實(shí)際監(jiān)測(cè)數(shù)據(jù)的基礎(chǔ)上，合理反演各種力學(xué)、滲流參數(shù)，以獲得符合北京山區(qū)地質(zhì)特征的流固耦合參數(shù)，并進(jìn)一步預(yù)測(cè)潛在的危險(xiǎn)區(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]. 巖土力學(xué)，22(4)：447-450.

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

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

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

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

徐晗，朱以文，蔡元奇，等，2005. 降雨入滲條件下非飽和土邊坡穩(wěn)定分析[J]. 巖土力學(xué)，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)秀人才資助項(xiàng)目（2016400617931G171）

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