地質(zhì)災(zāi)害流固耦合模型研究

丁靖洋

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

0 前言

我國(guó)是地質(zhì)災(zāi)害較為頻發(fā)的國(guó)家之一，改革開放以來，隨著大規(guī)?；A(chǔ)設(shè)施的建設(shè)，地質(zhì)災(zāi)害也進(jìn)入頻發(fā)期，且治理費(fèi)用高昂，動(dòng)輒幾千萬至上億元（王貢先，2004；王貢先，2005；黃潤(rùn)秋，2007），而其中降雨是地質(zhì)災(zāi)害發(fā)生的最重要誘因之一（劉艷輝等，2011）。鐵路交通因其具有線性工程的特點(diǎn)，沿途一旦發(fā)生地質(zhì)災(zāi)害，將會(huì)對(duì)其運(yùn)營(yíng)造成嚴(yán)重威脅，因此針對(duì)線性交通工程進(jìn)行地質(zhì)安全監(jiān)測(cè)、預(yù)警具有十分重要的意義。

“京張鐵路（居庸關(guān)隧道進(jìn)口段）地質(zhì)安全監(jiān)測(cè)示范工程”是“京津冀協(xié)同發(fā)展交通網(wǎng)絡(luò)地質(zhì)安全監(jiān)測(cè)預(yù)警系統(tǒng)”建設(shè)的重要組成部分，旨在對(duì)“京張客?！毖赝镜牡刭|(zhì)環(huán)境進(jìn)行監(jiān)測(cè)預(yù)警，并進(jìn)一步為構(gòu)建完整的京津冀交通網(wǎng)絡(luò)地質(zhì)安全監(jiān)測(cè)預(yù)警系統(tǒng)提供經(jīng)驗(yàn)和技術(shù)支持。

“京張鐵路（居庸關(guān)隧道進(jìn)口段）地質(zhì)安全監(jiān)測(cè)示范工程”的研究思路是：（1）在“京張客?！毖鼐€選擇一處突發(fā)性地質(zhì)災(zāi)害較發(fā)育的試驗(yàn)區(qū)域；（2）在試驗(yàn)區(qū)域內(nèi)布設(shè)一系列監(jiān)測(cè)設(shè)備，為地質(zhì)安全風(fēng)險(xiǎn)評(píng)估提供依據(jù)；（3）通過勘探手段獲得試驗(yàn)區(qū)域內(nèi)地質(zhì)信息，建立三維地質(zhì)體展示模型；（4）開展室內(nèi)土力學(xué)試驗(yàn)，獲取力學(xué)、滲流等一系列參數(shù)；（5）構(gòu)建流固耦合本構(gòu)方程，對(duì)試驗(yàn)區(qū)域進(jìn)行數(shù)值模擬研究，并與現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)進(jìn)行對(duì)比分析，驗(yàn)證方程及模擬結(jié)果的可靠性、合理性。

1 監(jiān)測(cè)設(shè)備埋設(shè)

試驗(yàn)區(qū)位于“京張客?！本佑龟P(guān)隧道進(jìn)口段，區(qū)內(nèi)發(fā)育有一條“V”型沖溝，溝內(nèi)含有大量物源體，在暴雨工況下極易發(fā)生泥石流、滑坡等地質(zhì)災(zāi)害。

根據(jù)勘察資料，該試驗(yàn)區(qū)30m深度范圍內(nèi)的地層分為新近沉積層、一般第四紀(jì)沉積層、長(zhǎng)城系高于莊組基巖。

現(xiàn)場(chǎng)埋設(shè)的監(jiān)測(cè)設(shè)備有：全球?qū)Ш叫l(wèi)星系統(tǒng)（Global Navigation Satellite System，GNSS）、陣列式位移計(jì)（Shape Accel Array，SAA）、精力水準(zhǔn)、土壓力計(jì)、孔隙水壓力計(jì)、土壤水分傳感器等。由于設(shè)備調(diào)試、數(shù)據(jù)接收等原因，本文數(shù)值模擬選用的現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)為土壤含水率、地下水平位移。

（1）土壤水分傳感器：現(xiàn)場(chǎng)布設(shè)2個(gè)土壤含水率測(cè)點(diǎn)，設(shè)備型號(hào)為CSF11型土壤水分傳感器，如圖1(a)所示。

（2）地下水平位移監(jiān)測(cè)（SAA）：在試驗(yàn)區(qū)內(nèi)布設(shè)1個(gè)地下水平位移監(jiān)測(cè)點(diǎn)，SAA埋設(shè)總長(zhǎng)為8m，如圖1(b)所示。

2 流固耦合本構(gòu)模型構(gòu)建

基于有效應(yīng)力原理（Terzaghi，1923）、非飽和滲流理論（Richards，1931）及Mohr-coulomb屈服準(zhǔn)則，本文作者已推導(dǎo)獲得了流固耦合本構(gòu)方程（丁靖洋，2017），即：

式（1）—式（3）中：n為孔隙介質(zhì)孔隙率，ρw為流體密度，Se為流體相的有效飽和度，uw為孔隙水壓力，k為固有滲透率，krel為非飽和到飽和過程中的相對(duì)滲透率，μw為流體粘性系數(shù)，u，為流體速度；為有效應(yīng)力，σij為主應(yīng)力，δij為克羅內(nèi)克符號(hào)；τ為剪應(yīng)力，?為內(nèi)摩擦角，c為粘聚力。

進(jìn)一步，為了考慮土體變形對(duì)孔隙率的影響，由孔隙比與體積單元變形關(guān)系可知：

式中：e，v，vs，v0，vs0和e0分別為變形過程中的孔隙比、單元體積、單元固體體積、初始單元體積、初始單元固體體積和初始孔隙比。式（4）中，表示單元體積變化量，可用第一應(yīng)變量表示為：

忽略單元變形過程中的固體介質(zhì)變形，則有：

將式（5）、（6）帶入式（4），得到：

則孔隙率與應(yīng)變的關(guān)系為：

式（8）中：n0為孔隙介質(zhì)的初始孔隙率，為第一應(yīng)變量。

綜合式（1）、（2）、（3）、（8），獲得一種新的考慮孔隙率演化的非飽和滲流-應(yīng)力耦合模型。

3 數(shù)值計(jì)算模型

3.1 三維離散模型

在現(xiàn)場(chǎng)勘察數(shù)據(jù)的基礎(chǔ)上，借助于OpenGeoSys（OGS）軟件及GMSH前處理軟件，建立了試驗(yàn)區(qū)域的三維地質(zhì)體模型（丁靖洋，2017），其中第一、二層土體的四面體單元尺寸為3m3，底層基巖單元尺寸為12m3，共計(jì)24498個(gè)（圖2）。

圖2 單元網(wǎng)格模型Fig.2 Element of the model

3.2 數(shù)值計(jì)算模型參數(shù)

本次研究過程中，主要針對(duì)第一層和第二層土體的降雨入滲變形過程進(jìn)行模擬，第三層地層為基巖，力學(xué)強(qiáng)度較高，滲透率較小，因此，數(shù)值計(jì)算中的關(guān)鍵地層為第一和第二層，第三層假設(shè)為線彈性屬性。

數(shù)值模擬，流體設(shè)定為水，密度為1000kg/m3，動(dòng)力粘性系數(shù)為0.001Pa·s。對(duì)試驗(yàn)區(qū)現(xiàn)場(chǎng)獲得的土樣進(jìn)行室內(nèi)土力學(xué)試驗(yàn)，確定模型中各地層的參數(shù)見表1。

表1 各地層模擬參數(shù)Tab.1 Parameters for the numerical modeling

基于現(xiàn)場(chǎng)土樣，通過試驗(yàn)測(cè)定了第一層和第二層土體的土水特征曲線（戚國(guó)慶等，2004）如圖3所示：

圖3 土水特征曲線及擬合參數(shù)Fig.3 Water retention curve and fitting parameters

圖3 中采用的擬合方程為Van Genuchten模型（Genuchten，1980；李少龍等，2006），即式（9）。計(jì)算獲得的第一層土體的a=0.098，m=0.385，第二層土體的a=0.212，m=0.378。

3.3 初始條件及邊界條件

根據(jù)地質(zhì)勘查結(jié)果及圖3中的參數(shù)，假定第三層基巖的含水率接近飽和，則設(shè)定第一層、第二層、第三層土體的初始飽和度分別為0.28、0.5、0.93。

圖4所示為現(xiàn)場(chǎng)監(jiān)測(cè)設(shè)備記錄的2017年6月21日—6月26日，共計(jì)約120小時(shí)內(nèi)的土壤含水率數(shù)據(jù)，經(jīng)過換算，得到有效飽和度隨時(shí)間的變化曲線。

數(shù)值計(jì)算中：設(shè)定模型四周為固定位移邊界，并固定模型底部垂直方向的位移；基于圖3擬合結(jié)果，將有效飽和度轉(zhuǎn)化為等效基質(zhì)吸力，作為模型上表面的入滲邊界條件，以此模擬降雨入滲過程引起的地表含水率和吸力變化；模型其他表面設(shè)定為不入滲邊界，在OGS軟件中表示無流體流入及流出。

圖4 地表有效飽和度隨時(shí)間變化曲線Fig.4 Effective saturation curve with time

4 模擬結(jié)果

4.1 位移監(jiān)測(cè)點(diǎn)模擬

現(xiàn)場(chǎng)布設(shè)有一個(gè)地下水平位移監(jiān)測(cè)點(diǎn)（SAA），實(shí)時(shí)記錄降雨過程中的測(cè)點(diǎn)位移變化情況；通過對(duì)數(shù)據(jù)進(jìn)行篩選、分析，選取6月23日、24日兩天持續(xù)降雨的監(jiān)測(cè)數(shù)據(jù)。由于測(cè)點(diǎn)測(cè)定的位移主要發(fā)生在x和y兩個(gè)方向，而z向位移基本可忽略不計(jì)，因此獲得位移隨SAA埋深的變化曲線如圖5所示。

圖5 現(xiàn)場(chǎng)監(jiān)測(cè)點(diǎn)位置及x、y位移數(shù)據(jù)Fig.5 The location of SAA equipment and x/y-displacement data

數(shù)值模擬計(jì)算獲得的該測(cè)點(diǎn)x、y向位移如圖6所示。對(duì)比現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)與數(shù)值模擬結(jié)果可知，該本構(gòu)模型能夠在一定程度上很好地反映實(shí)測(cè)結(jié)果。在SAA埋設(shè)處，地表位移較大，隨著深度的增加，位移呈減小趨勢(shì)；而數(shù)值計(jì)算結(jié)果能夠很好反映位移隨深度的變化規(guī)律，且對(duì)于潛層松散土體，即在不考慮構(gòu)造因素、易形成地質(zhì)災(zāi)害物源體的區(qū)域，基于非飽和滲流理論的模型計(jì)算結(jié)果能夠較好地匹配降雨情況下的位移數(shù)據(jù)。

同時(shí)應(yīng)該看到，由于OGS軟件的運(yùn)算基礎(chǔ)是連續(xù)介質(zhì)有限元理論，其對(duì)于介質(zhì)的均勻性假設(shè)以及網(wǎng)格尺度造成的計(jì)算誤差，造成了模擬結(jié)果與現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)存在一定的差別。但整體上，本次數(shù)值計(jì)算能夠有效反映降雨入滲過程引起的地層變形情況。

4.2 有效飽和度與位移模擬

本文分別模擬計(jì)算了20h、40h的位移及有效飽和度分布云圖，如圖7所示。由圖4可知，20h以內(nèi)，土壤含水率監(jiān)測(cè)位置處未達(dá)到飽和狀態(tài)；在20～40h，測(cè)點(diǎn)處持續(xù)呈飽和狀態(tài)。由圖7(b)可知，20h時(shí)模型左上部坡面最先達(dá)到飽和，且該處產(chǎn)生較大位移（圖7(a)），這是由于該處第一、第二層土體較薄，滲透壓力梯度較大，在滲透力驅(qū)動(dòng)下孔隙水將逐漸向模型中部匯集；至40h，隨著持續(xù)的降雨入滲，模型左上部出現(xiàn)大范圍飽和區(qū)域，相對(duì)應(yīng)的位移也在逐漸增大（圖7(c)），發(fā)生地質(zhì)災(zāi)害的風(fēng)險(xiǎn)增加。

圖6 現(xiàn)場(chǎng)監(jiān)測(cè)與數(shù)值模擬位移數(shù)據(jù)Fig.6 Displacement data of monitoring and numerical calculation

上述現(xiàn)象表明，有效飽和度增大、基質(zhì)吸力減小，是導(dǎo)致土體產(chǎn)生較大位移的重要原因：在降雨條件下，地層將被軟化，尤其是在地層分界處，由于力學(xué)性質(zhì)發(fā)生變化，極易產(chǎn)生軟弱結(jié)構(gòu)面；在滲流場(chǎng)-應(yīng)力場(chǎng)的耦合作用下，孔隙水壓力的變化導(dǎo)致了有效應(yīng)力的改變，并最終導(dǎo)致滑坡、泥石流等地質(zhì)災(zāi)害的發(fā)生。

圖7 位移及有效飽和度分布云圖Fig.7 Contour plots of displacement and effective saturation(a) Displacement in 20h (b) Effective saturation in 20h (c)Displacement in 40h (d) Effective saturation in 40h

5 結(jié)論

本文依托于“京張鐵路（居庸關(guān)隧道進(jìn)口段）地質(zhì)安全監(jiān)測(cè)示范工程”，借助于OpenGeoSys（OGS）軟件對(duì)試驗(yàn)區(qū)的降雨入滲過程（流固耦合）進(jìn)行了數(shù)值模擬研究，并對(duì)現(xiàn)場(chǎng)監(jiān)測(cè)設(shè)備采集到的數(shù)據(jù)進(jìn)行了擬合分析，主要成果及結(jié)論為：

（1）推導(dǎo)獲得了孔隙率隨應(yīng)變的演化方程，在有效應(yīng)力原理、非飽和滲流理論、Mohr-coulomb屈服準(zhǔn)則的基礎(chǔ)上，建立了一種新的非飽和滲流-應(yīng)力耦合模型；

（2）開展室內(nèi)土力學(xué)試驗(yàn)，獲得耦合模型所需力學(xué)、滲流參數(shù)；

（3）在現(xiàn)場(chǎng)采集到的土壤含水率數(shù)據(jù)的基礎(chǔ)上，開展數(shù)值模擬研究，并將位移模擬結(jié)果與現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比、擬合，計(jì)算結(jié)果表明該研究方法的可行性、合理性；

（4）降雨入滲過程，即有效飽和度逐漸增大、位移逐漸增加的過程，是地質(zhì)災(zāi)害發(fā)生的重要誘因。