地質災害流固耦合模型研究

丁靖洋

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

0 前言

我國是地質災害較為頻發的國家之一，改革開放以來，隨著大規模基礎設施的建設，地質災害也進入頻發期，且治理費用高昂，動輒幾千萬至上億元（王貢先，2004；王貢先，2005；黃潤秋，2007），而其中降雨是地質災害發生的最重要誘因之一（劉艷輝等，2011）。鐵路交通因其具有線性工程的特點，沿途一旦發生地質災害，將會對其運營造成嚴重威脅，因此針對線性交通工程進行地質安全監測、預警具有十分重要的意義。

“京張鐵路（居庸關隧道進口段）地質安全監測示范工程”是“京津冀協同發展交通網絡地質安全監測預警系統”建設的重要組成部分，旨在對“京張客專”沿途的地質環境進行監測預警，并進一步為構建完整的京津冀交通網絡地質安全監測預警系統提供經驗和技術支持。

“京張鐵路（居庸關隧道進口段）地質安全監測示范工程”的研究思路是：（1）在“京張客專”沿線選擇一處突發性地質災害較發育的試驗區域；（2）在試驗區域內布設一系列監測設備，為地質安全風險評估提供依據；（3）通過勘探手段獲得試驗區域內地質信息，建立三維地質體展示模型；（4）開展室內土力學試驗，獲取力學、滲流等一系列參數；（5）構建流固耦合本構方程，對試驗區域進行數值模擬研究，并與現場監測數據進行對比分析，驗證方程及模擬結果的可靠性、合理性。

1 監測設備埋設

試驗區位于“京張客專”居庸關隧道進口段，區內發育有一條“V”型沖溝，溝內含有大量物源體，在暴雨工況下極易發生泥石流、滑坡等地質災害。

根據勘察資料，該試驗區30m深度范圍內的地層分為新近沉積層、一般第四紀沉積層、長城系高于莊組基巖。

現場埋設的監測設備有：全球導航衛星系統（Global Navigation Satellite System，GNSS）、陣列式位移計（Shape Accel Array，SAA）、精力水準、土壓力計、孔隙水壓力計、土壤水分傳感器等。由于設備調試、數據接收等原因，本文數值模擬選用的現場監測數據為土壤含水率、地下水平位移。

（1）土壤水分傳感器：現場布設2個土壤含水率測點，設備型號為CSF11型土壤水分傳感器，如圖1(a)所示。

（2）地下水平位移監測（SAA）：在試驗區內布設1個地下水平位移監測點，SAA埋設總長為8m，如圖1(b)所示。

2 流固耦合本構模型構建

基于有效應力原理（Terzaghi，1923）、非飽和滲流理論（Richards，1931）及Mohr-coulomb屈服準則，本文作者已推導獲得了流固耦合本構方程（丁靖洋，2017），即：

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

進一步，為了考慮土體變形對孔隙率的影響，由孔隙比與體積單元變形關系可知：

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

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

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

則孔隙率與應變的關系為：

式（8）中：n0為孔隙介質的初始孔隙率，為第一應變量。

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

3 數值計算模型

3.1 三維離散模型

在現場勘察數據的基礎上，借助于OpenGeoSys（OGS）軟件及GMSH前處理軟件，建立了試驗區域的三維地質體模型（丁靖洋，2017），其中第一、二層土體的四面體單元尺寸為3m3，底層基巖單元尺寸為12m3，共計24498個（圖2）。

圖2 單元網格模型Fig.2 Element of the model

3.2 數值計算模型參數

本次研究過程中，主要針對第一層和第二層土體的降雨入滲變形過程進行模擬，第三層地層為基巖，力學強度較高，滲透率較小，因此，數值計算中的關鍵地層為第一和第二層，第三層假設為線彈性屬性。

數值模擬，流體設定為水，密度為1000kg/m3，動力粘性系數為0.001Pa·s。對試驗區現場獲得的土樣進行室內土力學試驗，確定模型中各地層的參數見表1。

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

基于現場土樣，通過試驗測定了第一層和第二層土體的土水特征曲線（戚國慶等，2004）如圖3所示：

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

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

3.3 初始條件及邊界條件

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

圖4所示為現場監測設備記錄的2017年6月21日—6月26日，共計約120小時內的土壤含水率數據，經過換算，得到有效飽和度隨時間的變化曲線。

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

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

4 模擬結果

4.1 位移監測點模擬

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

圖5 現場監測點位置及x、y位移數據Fig.5 The location of SAA equipment and x/y-displacement data

數值模擬計算獲得的該測點x、y向位移如圖6所示。對比現場監測數據與數值模擬結果可知，該本構模型能夠在一定程度上很好地反映實測結果。在SAA埋設處，地表位移較大，隨著深度的增加，位移呈減小趨勢；而數值計算結果能夠很好反映位移隨深度的變化規律，且對于潛層松散土體，即在不考慮構造因素、易形成地質災害物源體的區域，基于非飽和滲流理論的模型計算結果能夠較好地匹配降雨情況下的位移數據。

同時應該看到，由于OGS軟件的運算基礎是連續介質有限元理論，其對于介質的均勻性假設以及網格尺度造成的計算誤差，造成了模擬結果與現場實測數據存在一定的差別。但整體上，本次數值計算能夠有效反映降雨入滲過程引起的地層變形情況。

4.2 有效飽和度與位移模擬

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

圖6 現場監測與數值模擬位移數據Fig.6 Displacement data of monitoring and numerical calculation

上述現象表明，有效飽和度增大、基質吸力減小，是導致土體產生較大位移的重要原因：在降雨條件下，地層將被軟化，尤其是在地層分界處，由于力學性質發生變化，極易產生軟弱結構面；在滲流場-應力場的耦合作用下，孔隙水壓力的變化導致了有效應力的改變，并最終導致滑坡、泥石流等地質災害的發生。

圖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 結論

本文依托于“京張鐵路（居庸關隧道進口段）地質安全監測示范工程”，借助于OpenGeoSys（OGS）軟件對試驗區的降雨入滲過程（流固耦合）進行了數值模擬研究，并對現場監測設備采集到的數據進行了擬合分析，主要成果及結論為：

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

（2）開展室內土力學試驗，獲得耦合模型所需力學、滲流參數；

（3）在現場采集到的土壤含水率數據的基礎上，開展數值模擬研究，并將位移模擬結果與現場實測數據進行對比、擬合，計算結果表明該研究方法的可行性、合理性；

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