王文軒,夏 源

(1.桂林理工大學(xué) 環(huán)境科學(xué)與工程學(xué)院;廣西 桂林 541004;2.桂林理工大學(xué)巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心,廣西 桂林 541004)

0 引言

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

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

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

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

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

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

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

(1)

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

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

(2)

式(2)右側(cè)分子部分表示土體極限應(yīng)力值,分母部分為實(shí)際受力情況。為了對(duì)巖溶塌陷進(jìn)行預(yù)測評(píng)價(jià),將塌陷點(diǎn)的基本數(shù)據(jù)代入數(shù)學(xué)表達(dá)式計(jì)算,當(dāng)塌陷系數(shù)K>1,理論上不會(huì)產(chǎn)生巖溶地面塌陷;當(dāng)塌陷系數(shù)K<1,理論上發(fā)生巖溶地面塌陷。具體計(jì)算結(jié)果如表1所示。

表1 六種工況下的塌陷系數(shù)K值表

由表1可知,當(dāng)不考慮其他外界因素時(shí),僅考慮土體的自重壓力及氣壓差時(shí)。氣壓差為0、10、20 kPa覆蓋層土體受到的應(yīng)力達(dá)不到極限值,處于穩(wěn)定狀態(tài),巖溶地面塌陷很難發(fā)生。當(dāng)氣壓差為30、40、50 kPa時(shí)覆蓋層土體受到的應(yīng)力達(dá)到極限值,發(fā)生地面塌陷。

1.2 塌陷臨界氣壓分析

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

圖2 覆蓋層厚度與塌陷臨界氣壓的關(guān)系

圖3 巖溶塌陷離散元模型

2 建立離散元模型

2.1 模型尺寸

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

2.2 材料參數(shù)

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

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

表2 覆蓋層的力學(xué)參數(shù)

2.3 數(shù)值模擬過程

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

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

圖4 真空吸蝕模型圖

3 塌陷數(shù)值模擬預(yù)測分析

針對(duì)六種模擬工況,利用MatDEM進(jìn)行數(shù)值模擬,圖5為模擬結(jié)果,從中可以看出當(dāng)氣壓差為0、10、20 kPa覆蓋層土體未發(fā)生位移,處于穩(wěn)定狀態(tài),巖溶地面塌陷未發(fā)生。當(dāng)氣壓差為30、40、50 kPa時(shí)覆蓋層土體失穩(wěn),巖溶地面塌陷發(fā)生。數(shù)值模擬結(jié)果與力學(xué)解析法相符。

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

3.1 位移分析

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

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

表3 巖溶地面塌陷范圍表

3.2 應(yīng)力分析

由于模擬工況較多,本文選取未發(fā)生地面沉降的工況三,和發(fā)生地面塌陷的工況五進(jìn)行應(yīng)力分析,圖6與圖7分別是工況三和工況五情況下巖溶覆蓋層應(yīng)力狀態(tài)。

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

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

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

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

4 結(jié)語

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

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

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