真空吸蝕致巖溶塌陷的穩定性分析及其數值模擬

王文軒,夏 源

(1.桂林理工大學 環境科學與工程學院;廣西 桂林 541004;2.桂林理工大學巖溶地區水污染控制與用水安全保障協同創新中心,廣西 桂林 541004)

0 引言

地下水活動導致巖溶塌陷非常普遍,前人研究后總結出幾種致塌機制,如潛蝕致塌,真空吸蝕致塌等[1],真空吸蝕致塌是重要的致塌機理之一。隨著人類活動的加劇,生產生活中大量抽排地下水[2],導致巖溶空腔中產生氣真空值,此時若巖溶空腔上部貫通致巖基面,與覆蓋層接觸,且覆蓋層氣密性較好,導致部分處于平衡狀態的巖溶土洞發生地面塌陷[3]。這種現象被稱為“真空吸蝕”[4]。針對該原理,學者們進行的大量研究,姚春梅等建立巖溶塌陷預警系統模型對巖溶地質災害進行預警[5];熊啟華等[6]結合塌陷成因建立真空吸蝕相關力學模型,并對巖溶塌陷進行分析預測;張鑫等[7]通過理論分析和室內模型試驗,定量化的研究水位下降對巖溶塌陷的影響。鐘瑤[8]利用FLAC3D數值模擬軟件模擬水位下降時巖溶塌陷演化過程,從應力應變角度分析了巖溶塌陷的發育過程;于林弘等[9]通過改變土洞開口大小及水位降幅,對三元覆蓋型巖溶塌陷模式進行模擬;方先知等[10]利用有限元軟件對巖溶體穩定性進行分析;張曉宸等[11]運用Modflow軟件,分析了地下水位變化對巖溶塌陷的影響。雖然學者已經對真空吸蝕導致的巖溶塌陷進行了大量研究,但使用最多的是模擬連續性、均勻性和小變形的FLAC3D、GMS等軟件,基于網格方法,可以模擬出塌陷的發生,而后期的大變形破壞則不能模擬,而離散元法可以模擬非連續性、不均勻性和大變形破壞,此類研究較少,本文使用離散元軟件MatDEM對巖溶塌陷進行模擬研究。

1 真空吸蝕致巖溶塌陷分析

1.1 覆蓋層土體極限平衡公式

在實際的工程中,不同的地下水位變化導致土洞中出現不同真空度,對覆蓋層的穩定產生不同的影響。根據文獻[12],氣壓差ΔP一般情況下取值范圍在0～50 kPa之間,所以分別選取0、10、20、30、40、50 kPa進行分析,通過不同的氣壓差模擬不同水位下降深度對巖溶塌陷的影響,根據選取不同的氣壓差,本實驗將真空吸蝕分為六種工況進行模擬分析。

目前,國內對于分析覆蓋層土體在巖溶塌陷過程中的受力條件與穩定狀態進行了研究,本文通過前人[13]的方法并結合二維離散元模型進行分析,基本假定如下:(1)在真空吸蝕作用下,塌陷體剪切破裂面從拱腳處垂直發展至地表;(2)考慮到覆蓋層多為黏性土體孔隙率較小,地下水驟降條件下,外界氣體未能及時補給,洞內氣壓暫時保持不變;(3)地下水位始終位于巖基面以下,土體容重始終不變;(4)在不受外力作用下自然塌陷所形成的土洞體積較小,進行力學分析時忽略此處脫落土體。

如圖1所示,土體重度G=γSH;水位下降引起的負壓F=ΔPS;覆蓋層破壞體側面摩擦力

圖1 覆蓋層土體受力分析圖

(1)

其中K0為土體側壓力系數;γ為覆土天然容重;H為覆蓋層厚度;C為土體黏聚力;φ為內摩擦角;D為覆蓋層施加氣壓的面積;W為二維離散元模型發生側面位移的周長;為巖溶裂隙寬度,取W=3 m。

由此得到基于二維離散元模型得到的覆蓋層土體極限平衡狀態時平衡公式:

(2)

式(2)右側分子部分表示土體極限應力值,分母部分為實際受力情況。為了對巖溶塌陷進行預測評價,將塌陷點的基本數據代入數學表達式計算,當塌陷系數K>1,理論上不會產生巖溶地面塌陷;當塌陷系數K<1,理論上發生巖溶地面塌陷。具體計算結果如表1所示。

表1 六種工況下的塌陷系數K值表

由表1可知,當不考慮其他外界因素時,僅考慮土體的自重壓力及氣壓差時。氣壓差為0、10、20 kPa覆蓋層土體受到的應力達不到極限值,處于穩定狀態,巖溶地面塌陷很難發生。當氣壓差為30、40、50 kPa時覆蓋層土體受到的應力達到極限值,發生地面塌陷。

1.2 塌陷臨界氣壓分析

為得到不同覆蓋層厚度的塌陷臨界氣壓差ΔP,令式(2)中的K=1,得到覆蓋層厚度H與塌陷臨界氣壓差ΔP的關系(圖2)。如圖3所示,塌陷臨界氣壓差與覆蓋層厚度呈正相關,且隨著覆蓋層厚度的增大臨界氣壓差增大的速率增加;當覆蓋層厚度大于9.16 m時,覆蓋層土體始終保持穩定,覆蓋層厚度小于6.71 m覆蓋層土體始終處于失穩狀態。

圖2 覆蓋層厚度與塌陷臨界氣壓的關系

圖3 巖溶塌陷離散元模型

2 建立離散元模型

2.1 模型尺寸

如圖3所示,巖溶塌陷二維離散元模型長×寬為30 m×12 m,根據礦區水位地質條件,結合數值模擬經驗,建立巖層厚度為4 m,覆蓋層厚度為8 m,巖溶裂隙寬度為3 m。顆粒單元共37 470個,其中固定單元12 435個(巖土層),活動單元25 035(覆蓋層),將巖層單元顆粒設置為固定單元減小計算時間。

2.2 材料參數

設置合適的力學參數在離散元數值建模中設置非常重要,MatDEM離散元隨機堆積模型的宏觀力學參數與微觀力學參數存在解析解,即轉換公式[14]。線彈性接觸模型的五個微觀力學可以通過轉換公式,由材料的五個宏觀力學性質計算得到

一般情況下,根據轉換公式計算得到的緊密堆積模型力學性質小于理論值,為了獲取模型顆粒合適的微觀參數,通過MatDEM的MatTraining文件實現自動訓練材料,經過多次數值測試后得到模型最終的材料參數。根據文獻[15]可以得到巖溶覆蓋層宏觀力學參數值,經過材料訓練得到時顆粒單元適合顆粒單元的宏觀力學參數值如表2。

表2 覆蓋層的力學參數

2.3 數值模擬過程

在MatDEM離散元軟件的空箱子中生成一系列具有一定速度的顆粒單元,隨后進行重力沉積壓實,建立離散元模型;利用MatDEM的切割函數C-tool將模型切割為所需要的巖溶塌陷模型,將巖層設置為墻單元減少計算時間,隨后賦材料于覆蓋層平衡模型,重新計算模型顆粒的受力狀態,使模型重新平衡。

在賦予材料后進行通過過濾矩陣篩選,得到巖溶裂隙上部地表的顆粒單元,根據不同的內外壓力差賦予,篩選出來表層顆粒相應的體力;根據不同的真空值對巖溶裂隙正上方的覆蓋層施加一個短時間豎直向下的壓力來模型巖溶裂隙內外氣壓差導致的壓力,以模擬空蝕對巖溶覆蓋層的影響,隨后進行自然迭代。如圖4所示,在地表顆粒單元上施加豎向體力,模擬真空吸蝕原理對覆蓋層土體的作用。

圖4 真空吸蝕模型圖

3 塌陷數值模擬預測分析

針對六種模擬工況,利用MatDEM進行數值模擬,圖5為模擬結果,從中可以看出當氣壓差為0、10、20 kPa覆蓋層土體未發生位移,處于穩定狀態,巖溶地面塌陷未發生。當氣壓差為30、40、50 kPa時覆蓋層土體失穩,巖溶地面塌陷發生。數值模擬結果與力學解析法相符。

圖5 巖溶塌陷演化過程圖(單位:m)

3.1 位移分析

如圖5所示,在工況一、二、三情況下巖溶塌陷狀況類似,在經過自然塌陷后施加0、10、20 kPa的豎向壓力均不能造成覆蓋層發生地面塌陷,且根據單元連接圖可以看出土體顆粒未發生明顯位移,此時覆蓋層處于穩定狀況,但是施加豎向應力后,巖溶土洞開始發育;在工況四中,覆蓋層在氣壓差為30 kPa的豎向壓力作用下,巖溶裂隙處覆蓋層向上發育,部分土體脫落。但未形成較大規模,所以覆蓋層土體未發生明顯豎向位移;在工況四、五、六中均發生地面塌陷。

經過對圖5 d-f的單元檢測可以得到巖溶地面塌陷的范圍(表3)。如表3所示,當發生巖溶地面塌陷時,隨著ΔP增加,巖溶地面塌陷寬度增大。

表3 巖溶地面塌陷范圍表

3.2 應力分析

由于模擬工況較多,本文選取未發生地面沉降的工況三,和發生地面塌陷的工況五進行應力分析,圖6與圖7分別是工況三和工況五情況下巖溶覆蓋層應力狀態。

圖6 工況三覆土層應力圖(單位:Pa)

圖7 工況五覆土層應力圖(單位:Pa)

由圖6 a可知,覆蓋層土體僅發生了內部塌陷,形成巖溶土洞,土洞周圍顆粒在水位升降結束后重新排列,應力狀態也隨之重新分布,在拱效應的作用下土洞上部土體水平應力均大于相同深度土體的水平應力,土洞空腔中的塌落土體的水平應力接近于0,即該處土體未受到形成土洞所產生的水平應力;圖6 b為豎直方向應力狀態圖,整體上豎直方向上應力狀態隨著土層深度的增加而增大,土洞周圍的土體單元經過重新排列和應力重分布后,在洞趾處出現應力集中,土洞上層土體由于拱效應的影響,該處土體的豎向應力明顯低于相同深度土體的豎向應力,

由圖7 a可知,覆蓋層土體失穩,發生地面塌陷,土體水平方向應力整體上隨著土層深度的增大而增大,塌落土體經過土體單元重新排列后,水平方向應力重新分布呈拱狀,土層表面水平方向上應力較小或者為拉應力狀態;圖7 b為豎直方向應力狀態圖,整體上豎直方向上應力狀態隨著土層深度的增加而增大,塌落土體豎向應力明顯小于同一土層深度的其他土體。

4 結語

(1)基于二維離散元模型及前人的研究得到真空吸蝕致巖溶塌陷的解析公式,并利用MatDEM離散元軟件進行模擬。離散元數值模擬在模擬真空吸蝕致巖溶塌陷中可以較為準確地預測巖溶塌陷的穩定性,

(2)在覆蓋層厚度與巖溶裂隙寬度不變的情況下,發生巖溶地面塌陷的塌陷范圍與氣壓差呈正相關。

(3)巖溶裂隙寬度一定時,塌陷臨界氣壓差ΔP和覆蓋層H厚度呈正相關,且隨著覆蓋層厚度的增大塌陷臨界氣壓差增大的速率增加;當覆蓋層厚度H大于或小于一定值后覆蓋層穩定不再受氣壓差影響。