覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型研究

熊啟華，高 旭，徐 慶，王芮瓊，李祖春，陶 良

(1.湖北省地質(zhì)環(huán)境總站，湖北 武漢 430039;2.中國(guó)地質(zhì)大學(xué)(武漢)工程學(xué)院，湖北 武漢 430074;3.武漢市測(cè)繪研究院，湖北 武漢 430022)

覆蓋型巖溶塌陷具有表觀突發(fā)性和內(nèi)部漸進(jìn)性的特點(diǎn)，塌陷的發(fā)生是在覆蓋層內(nèi)土洞逐漸演化擴(kuò)展至極限尺寸后，頂板強(qiáng)度不足以支撐其自重而突然失穩(wěn)垮塌。因此，預(yù)測(cè)不同水動(dòng)力環(huán)境下巖溶覆蓋層土洞的發(fā)育規(guī)模和極限規(guī)模尺寸對(duì)于評(píng)價(jià)覆蓋層穩(wěn)定性意義重大。然而受制于土洞發(fā)育的隱蔽性和物探解譯的多解性，直接通過(guò)技術(shù)手段查明土洞發(fā)育尺寸的難度較大或精度不高。鑒于此，從土洞發(fā)育物理力學(xué)機(jī)理出發(fā)來(lái)建立一種定量預(yù)測(cè)土洞發(fā)育規(guī)模和穩(wěn)定性評(píng)價(jià)的數(shù)值模型(即覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型)并加以驗(yàn)證是本文的研究目標(biāo)。

目前，以地下水變動(dòng)為主要誘發(fā)因素的巖溶覆蓋層中土洞發(fā)育的機(jī)理主要有潛蝕論和真空吸蝕論兩種[1]。潛蝕論認(rèn)為在巖溶溶隙開(kāi)口附近的水力梯度超過(guò)覆蓋層土體臨界水力梯度而發(fā)生潛蝕或剝蝕，形成土洞并逐漸擴(kuò)展；真空吸蝕論認(rèn)為地下水水位在覆蓋層與溶腔開(kāi)口交界處上下波動(dòng)或快速抽吸巖溶水形成負(fù)壓而使土層剝蝕擴(kuò)洞。

而目前針對(duì)覆蓋型巖溶塌陷機(jī)理的研究方法主要有解析法、物理模型試驗(yàn)法和數(shù)值模擬法。其中，解析法多是基于普氏平衡拱理論探索臨界土洞的規(guī)模大小[2]或采用極限平衡理論分析給定幾何尺寸下土洞頂板的穩(wěn)定性狀態(tài)[3]。雖然利用物理模型試驗(yàn)方法來(lái)研究覆蓋型巖溶塌陷機(jī)理的學(xué)者較多[4-5]，但定量研究土洞動(dòng)態(tài)發(fā)育特征的物理模型試驗(yàn)較少，僅見(jiàn)有：吳慶華等[6]采用恒壓取樣技術(shù)獲取巖溶塌陷物理模型試驗(yàn)過(guò)程中砂土漏失質(zhì)量，用來(lái)定量評(píng)估土洞的發(fā)育過(guò)程；張少波等[7]通過(guò)在物理模型土層中間隔一定距離布設(shè)孔隙水壓力自動(dòng)監(jiān)測(cè)點(diǎn)，通過(guò)監(jiān)測(cè)土洞發(fā)育過(guò)程中造成的不同位置處孔隙水壓力變化的差異來(lái)判定土洞的發(fā)育階段。然而，目前物理模型試驗(yàn)很難判定不同發(fā)育階段土洞所對(duì)應(yīng)的穩(wěn)定性狀態(tài)。鑒于此，部分學(xué)者利用數(shù)值模擬手段來(lái)評(píng)估土洞的穩(wěn)定性，如：孫金輝[8]結(jié)合物理模型試驗(yàn)并基于FLAC3D建立了不同規(guī)模尺寸土洞開(kāi)挖數(shù)值模型，探索覆蓋層地表沉降和內(nèi)部應(yīng)力-應(yīng)變過(guò)程；陳冬琴等[9]建立了巖溶塌陷地下水-力學(xué)耦合數(shù)值模型，并基于此對(duì)土-巖交界處溶洞塌陷進(jìn)行了地下水水位升降及降雨入滲作用對(duì)巖溶塌陷的影響分析。

上述研究基本都是人為設(shè)定土洞尺寸再進(jìn)行數(shù)值模擬，并不是自適應(yīng)的土洞擴(kuò)展演化過(guò)程。鑒于此，本文以潛蝕論為基本出發(fā)點(diǎn)，建立了覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型，并用于土洞發(fā)育規(guī)模預(yù)測(cè)。首先，以武漢市烽火村巖溶塌陷實(shí)例為分析對(duì)象，概化其巖溶塌陷工程地質(zhì)模型，并建立滲流-力學(xué)單向耦合控制方程，同時(shí)結(jié)合土洞剝蝕判據(jù)編制覆蓋型巖溶塌陷動(dòng)態(tài)演化計(jì)算程序；然后，基于此程序預(yù)測(cè)在不同水動(dòng)力環(huán)境條件下覆蓋層中潛在土洞的發(fā)育規(guī)模，并提出臨界致塌水位差和極限土洞尺寸；最后，基于普氏平衡拱理論對(duì)預(yù)測(cè)結(jié)果進(jìn)行驗(yàn)證分析。

1 覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型

1.1 覆蓋型巖溶塌陷工程地質(zhì)模型概化

圖1 烽火村覆蓋型巖溶塌陷地質(zhì)剖面圖

據(jù)此，將烽火村覆蓋型巖溶塌陷概化為如圖2所示的工程地質(zhì)模型，其覆蓋層總厚度為30 m，黏砂層厚度比為1∶5，模型尺寸為100 m×35 m，其中基巖厚度為5 m，網(wǎng)格剖分為1 m×1 m單元。

圖2 烽火村覆蓋型巖溶塌陷工程地質(zhì)概化模型

該巖溶塌陷是長(zhǎng)期人工抽取地下水誘發(fā)的，定性分析認(rèn)為覆蓋層中孔隙水位常常高于巖溶水位，導(dǎo)致存在孔隙水向巖溶水補(bǔ)給滲流過(guò)程中當(dāng)水力梯度超過(guò)某一臨界值時(shí)出現(xiàn)砂土潛蝕或剝蝕現(xiàn)象。土洞剝蝕類似地下工程開(kāi)挖卸荷，必然伴隨著洞周應(yīng)力重分布、地表沉降、塑性區(qū)發(fā)展等一系列力學(xué)響應(yīng)；加之，地下水滲流場(chǎng)的存在不僅是土洞剝蝕的直接動(dòng)力源，同時(shí)也造成對(duì)土洞的滲流荷載。因此，在滲流荷載與卸荷荷載聯(lián)合作用下，當(dāng)土洞的應(yīng)力、應(yīng)變超過(guò)其極限狀態(tài)時(shí)則產(chǎn)生巖溶塌陷。

鑒于此，覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型就是建立起能夠描述地下水流場(chǎng)中土洞擴(kuò)展的滲流-力學(xué)單向耦合控制方程(此處單向耦合是指巖土體的應(yīng)力-應(yīng)變不改變其滲透性質(zhì)，而只有滲透力的單方面力學(xué)作用)，且通過(guò)土洞剝蝕判據(jù)和土洞應(yīng)力-應(yīng)變狀態(tài)來(lái)控制土洞擴(kuò)展演化進(jìn)程。下面闡述滲流-力學(xué)單向耦合控制方程式及其所用參數(shù)和邊界條件，以及土洞動(dòng)態(tài)演化的計(jì)算流程。

1.2 滲流-力學(xué)單向耦合控制方程式及其邊界條件

數(shù)值模型不考慮剝離后土體的運(yùn)動(dòng)狀態(tài)，且假定剝離前土洞巖土體介質(zhì)處于小變形狀態(tài)，故采用二維等效連續(xù)介質(zhì)有限元數(shù)值模型模擬基巖和覆蓋層地下水流場(chǎng)，以及在地下水滲流作用下土洞剝蝕過(guò)程中變形和塑性區(qū)的發(fā)展過(guò)程，故模型的力平衡方程如下：

(1)

(2)

式中：γg為巖土層重度(kN/m3)；VA為被剝蝕單元體積(m3)；[B]為關(guān)于應(yīng)變與位移之間關(guān)系的幾何矩陣；[N]為典型四節(jié)點(diǎn)等參單元形函數(shù)。

另一方面，作用于非剝蝕單元節(jié)點(diǎn)上的滲透力則是與地下水滲流梯度有關(guān)[10]，其滲流荷載的數(shù)學(xué)表達(dá)式為

(3)

式中：γw為地下水的單位重度(kN/m3)，取10 kN/m3；Ω為土洞外滲流區(qū)域；H表示總水頭(m)，其與壓力水頭p(m)的關(guān)系式為H=Z+p，其中Z表示節(jié)點(diǎn)豎向高度(m)。

一方面，考慮到研究區(qū)覆蓋層內(nèi)土洞埋深一般在30 m以淺，所處圍壓環(huán)境較小，巖土體材料一般不會(huì)出現(xiàn)應(yīng)變硬化特性，但由于采用巖土體應(yīng)變軟化彈塑性本構(gòu)模型來(lái)描述其達(dá)到峰值強(qiáng)度后應(yīng)力-應(yīng)變關(guān)系的參數(shù)難以測(cè)定，故綜合考慮采用巖土體屈服極限不隨應(yīng)力狀態(tài)改變而改變的理想彈塑性模型來(lái)描述土洞的非線性應(yīng)力-應(yīng)變過(guò)程；另一方面，考慮到巖土體材料屬于摩擦型材料，故采用最常用的Mohr-Coulomb屈服準(zhǔn)則。本模型不考慮巖土體塑性變形過(guò)程中的剪脹效應(yīng)，因此采用非關(guān)聯(lián)流動(dòng)法則來(lái)描述單元進(jìn)入塑性狀態(tài)后的應(yīng)變規(guī)律，其塑性勢(shì)函數(shù)中的剪脹角參數(shù)設(shè)為0。

為了展示土洞演化中存在的塑性區(qū)分布情況，定義塑性應(yīng)變不變量為[11]

(4)

考慮到孔隙水位以上土層處于非飽和狀態(tài)，采用變飽和地下水穩(wěn)定流控制方程來(lái)統(tǒng)一描述整個(gè)模型中的滲流場(chǎng)：

(5)

式中:p為壓力水頭(m)；K(p)為巖土層滲透系數(shù)(m/d)，對(duì)處于飽和狀態(tài)(p≥0)的巖土層滲透系數(shù)取飽和滲透系數(shù)，而對(duì)處于孔隙水位以上的非飽和巖土層滲透系數(shù)則是壓力水頭的函數(shù)，采用最簡(jiǎn)單的指數(shù)型模型表征為

(6)

式中：Ks為巖土層飽和滲透系數(shù)(m/d)；αm為孔隙尺寸分布參數(shù)(m-1)。

如圖2所示，模型力學(xué)邊界條件設(shè)置為左右邊界約束水平位移；而底邊界同時(shí)約束水平和垂直位移；上邊界為自由邊界。對(duì)于水力邊界條件，BC段和DE段設(shè)置為總水頭H1的定水頭邊界，代表弱承壓的孔隙含水層水位；而AB段、AF段、EF段則設(shè)置為代表巖溶水位H2的可變定水頭邊界，即通過(guò)同時(shí)調(diào)整此三段水頭來(lái)模擬不同的巖溶水動(dòng)力環(huán)境條件；上邊界CD段為零流量邊界。

據(jù)文獻(xiàn)[14]可知，能使類似地層土洞出現(xiàn)剝蝕的臨界水力梯度(Ibs)取0.5。通過(guò)歸納總結(jié)相關(guān)文獻(xiàn)[15]和《武漢市巖溶塌陷調(diào)查子項(xiàng)目成果報(bào)告》(編號(hào)為1212011220189)確定數(shù)值模型所用的滲流和力學(xué)計(jì)算參數(shù)，見(jiàn)表1。

表1 數(shù)值模型所用的計(jì)算參數(shù)

上述公式(5)和(6)所代表的滲流場(chǎng)計(jì)算則采用變飽和滲流計(jì)算程序vsaft2[16]，得出的水頭場(chǎng)根據(jù)公式(3)計(jì)算出作用在節(jié)點(diǎn)上的滲透荷載;然后修改文獻(xiàn)[17]中第6.4.5節(jié)關(guān)于組裝節(jié)點(diǎn)荷載的子程序LOADPS，將滲透荷載和根據(jù)公式(2)計(jì)算的土洞剝蝕后卸荷荷載組裝成總體荷載矢量；最后由文獻(xiàn)[17]提供的非線性彈塑性力學(xué)計(jì)算程序完成計(jì)算。

1.3 土洞演化的計(jì)算流程

結(jié)合巖溶塌陷工程地質(zhì)概化模型和滲流-應(yīng)力單向耦合控制方程，可得到土洞動(dòng)態(tài)演化計(jì)算流程(見(jiàn)圖3)如下：

圖3 土洞動(dòng)態(tài)演化的計(jì)算流程

(1) 首先，在溶隙開(kāi)口部位設(shè)置一個(gè)幾何較小的初始土洞，土洞和溶隙內(nèi)都填充擾動(dòng)土，而其他各巖土層都按原生地層幾何參數(shù)賦值。

(2) 在給定的水力邊界條件下按照公式(5)求解滲流場(chǎng)，進(jìn)而得到水力梯度空間分布。一方面，利用該水力梯度空間分布情況來(lái)搜索出與土洞緊鄰的土層單元中水力梯度I大于或等于臨界水力梯度Ibs(即I≥Ibs)的單元，并納入到下一次滲流場(chǎng)計(jì)算中的擾動(dòng)土單元，同時(shí)作為下一次應(yīng)力場(chǎng)計(jì)算的剝蝕單元；另一方面，基于水力梯度空間分布，即可根據(jù)公式(3)計(jì)算出作用在模型中每個(gè)單元的滲透力，并結(jié)合土洞單元?jiǎng)兾g卸荷作用，計(jì)算模型的力學(xué)響應(yīng)(包括變形和塑性區(qū)等)。

(3) 根據(jù)新形成的土洞再次進(jìn)行滲流場(chǎng)和應(yīng)力場(chǎng)的計(jì)算，并重復(fù)上述步驟，直到滿足如下兩個(gè)終止條件之一：①緊鄰?fù)炼吹耐翆訂卧μ荻菼全都小于臨界水力梯度Ibs(即I

2 不同水動(dòng)力環(huán)境下土洞發(fā)育規(guī)模預(yù)測(cè)

本文設(shè)置模型孔隙水位H1=32 m用于模擬覆蓋層的承壓水環(huán)境且保持不變，而通過(guò)降低巖溶水位H2用于模擬孔隙水位與巖溶水位之差ΔH(見(jiàn)圖2)穩(wěn)定在4 m、6 m、8 m、10 m、12 m、14 m、16 m這七種水動(dòng)力環(huán)境，借以探索不同水位差環(huán)境持續(xù)作用下覆蓋層中潛在的穩(wěn)定土洞發(fā)育規(guī)模，以及臨界致塌水位差和極限土洞規(guī)模。

先分析在ΔH=16 m水動(dòng)力環(huán)境持續(xù)作用下土洞逐漸擴(kuò)展過(guò)程中模型滲流場(chǎng)、水力梯度、塑性區(qū)以及位移場(chǎng)情況。設(shè)置地表沉降監(jiān)測(cè)點(diǎn)位于如圖2所示藍(lán)色圓點(diǎn)處，并記錄在ΔH=16 m水動(dòng)力環(huán)境持續(xù)作用下土洞演化過(guò)程中地表監(jiān)測(cè)沉降量，見(jiàn)圖4。

圖4 ΔH=16 m水動(dòng)力環(huán)境持續(xù)作用下土洞演化過(guò)程中地表監(jiān)測(cè)沉降量(uz)變化曲線

由圖4可見(jiàn)，土洞在剝蝕次數(shù)(Nslo)從第4次到第5次時(shí)，地表監(jiān)測(cè)沉降量出現(xiàn)劇烈突變，說(shuō)明在第5次土洞剝蝕后覆蓋層出現(xiàn)地面塌陷，而第4次土洞剝蝕后形成的土洞接近極限土洞規(guī)模大小。鑒于此，本次只展示前4次土洞演化過(guò)程中滲流場(chǎng)和力學(xué)響應(yīng)的變化規(guī)律,見(jiàn)圖5和圖6。

由圖5可見(jiàn)：土洞演化過(guò)程中滲流場(chǎng)的流線主要向土洞中心匯聚，孔隙水向巖溶水補(bǔ)給，隨著土洞逐漸擴(kuò)展演化，越靠近土洞總水頭越低[見(jiàn)圖5(a)]；土洞演化過(guò)程中水力梯度小于臨界水力梯度的范圍基本不變(即覆蓋層中半圓形紅色圈定范圍)，說(shuō)明孔隙水位與巖溶水位之差ΔH一旦固定時(shí)，土洞擴(kuò)展范圍也基本確定[見(jiàn)圖5(b)]。

圖5 ΔH=16 m水動(dòng)力環(huán)境持續(xù)作用下土洞演化過(guò)程中滲流場(chǎng)和水力梯度的變化規(guī)律

由圖6可以進(jìn)一步地直觀看出：土洞剝蝕后規(guī)模大小由最初的洞高1 m和洞寬3 m演變?yōu)槎锤? m和洞寬9 m的極限土洞規(guī)模，且隨著土洞的擴(kuò)展演化，土洞周圍砂土層中塑性區(qū)范圍和塑性應(yīng)變量逐漸增大，直到第4次剝蝕后塑性區(qū)擴(kuò)展到砂土層頂部，即代表此時(shí)覆蓋層出現(xiàn)大面積塑性破壞[見(jiàn)圖6(a)]；初始土洞形成后無(wú)論地表還是土洞周圍位移都較小，而當(dāng)土洞逐漸擴(kuò)展，地表和土洞周邊位移都顯著增大[見(jiàn)圖6(b)]，但從圖6(b)中代表位移矢量方向的紅色箭頭線可知，地表主要以水平位移為主，而土洞中心線附近則以垂直位移為主,這能夠較好地解釋實(shí)際塌陷坑周邊出現(xiàn)拉裂縫且塌陷坑中心豎直沉降量最大這一現(xiàn)象。

圖6 ΔH=16 m水動(dòng)力環(huán)境持續(xù)作用下土洞演化過(guò)程中塑性區(qū)和位移場(chǎng)的變化規(guī)律

以上僅從ΔH=16 m水動(dòng)力環(huán)境條件下土洞演化過(guò)程及塌陷機(jī)理進(jìn)行闡述，而在不同水位差水動(dòng)力環(huán)境下土洞演化過(guò)程如圖7所示，圖中僅展示了模型x=40 m到x=60 m、z=3 m到z=13 m的局部塑性區(qū)發(fā)展情況[見(jiàn)圖7(a)至圖7(l)]及其對(duì)應(yīng)穩(wěn)定土洞規(guī)模和地表監(jiān)測(cè)沉降曲線[見(jiàn)圖7(m)]。其中,w代表土洞洞寬(m);h代表土洞洞高(m);uz代表地表中心監(jiān)測(cè)沉降量(m)。

圖7 不同水位差水動(dòng)力環(huán)境下土洞演化過(guò)程中塑性區(qū)及其對(duì)應(yīng)穩(wěn)定土洞規(guī)模和地表監(jiān)測(cè)沉降曲線

由圖7(a)至圖7(l)可見(jiàn):當(dāng)水位差ΔH在4～8 m時(shí)，土洞保持初始土洞大小并未擴(kuò)展，說(shuō)明在此水動(dòng)力環(huán)境下砂性土層水力梯度I超過(guò)臨界水力梯度(Ibs=0.5)的范圍較小，沒(méi)有土層單元被剝蝕；而當(dāng)水位差ΔH穩(wěn)定在10 m時(shí)，土洞發(fā)生剝蝕而擴(kuò)展至洞高3 m和洞寬5 m，塑性變形也有所增加;同樣地，當(dāng)水位差ΔH穩(wěn)定在12 m和14 m時(shí)，土洞剝蝕后形成的穩(wěn)定土洞規(guī)模相應(yīng)變大；當(dāng)水位差ΔH為14 m水動(dòng)力環(huán)境下土洞最終并未出現(xiàn)塌陷現(xiàn)象，即沒(méi)有出現(xiàn)類似如圖4所示的地表監(jiān)測(cè)沉降量突變或計(jì)算不收斂的情況，這說(shuō)明水位差ΔH需要達(dá)到16 m左右(即塌陷臨界水位差為16 m),才可能出現(xiàn)最終塌陷現(xiàn)象。

由圖7(m)可以看出：土洞洞寬曲線斜率總體大于洞高曲線斜率，說(shuō)明土洞演化過(guò)程中洞寬擴(kuò)展速率大于洞高擴(kuò)展速率；地表監(jiān)測(cè)沉降量隨著水位差ΔH增加而增加，地表沉降速率大致也可分為三個(gè)階段，即①當(dāng)ΔH=4～8 m時(shí)地表沉降速率最小，②當(dāng)ΔH=8～12 m時(shí)地表沉降速率居中，③當(dāng)ΔH=12～16 m時(shí)地表沉降速率最大，說(shuō)明地表沉降速率會(huì)隨水位差所處范圍不同而不同，水位差越大，地表沉降速率越大。

3 土洞預(yù)測(cè)結(jié)果驗(yàn)證

(7)

本文基于覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型預(yù)測(cè)的土洞極限洞高相對(duì)普氏平衡拱理論求解的土洞極限洞高偏小，其原因在于普氏平衡拱理論的致塌力只有自重作用，而本文數(shù)值模型中考慮了滲透力和自重的聯(lián)合作用，使得在相對(duì)較小的土洞規(guī)模下仍會(huì)產(chǎn)生失穩(wěn)破壞。

4 結(jié) 論

本文結(jié)合彈塑性理論和變飽和滲流理論，以臨界水力梯度為土洞剝蝕判據(jù)，提出了覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型，探索了在不同孔隙水位與巖溶水位之差水動(dòng)力環(huán)境下覆蓋層中潛在土洞的發(fā)育規(guī)模，并加以驗(yàn)證，得出如下結(jié)論：

(1) 在孔隙水位與巖溶水位之差ΔH小于8 m的水動(dòng)力環(huán)境持續(xù)作用下，土洞將保持較小規(guī)模而不會(huì)進(jìn)一步擴(kuò)展;而當(dāng)水位差ΔH超過(guò)8 m后將存在土洞剝蝕現(xiàn)象，且水位差越大，最終形成的土洞規(guī)模越大。

(2) 當(dāng)水位差ΔH達(dá)到16 m的臨界水位差時(shí)，土洞將逐步剝蝕至產(chǎn)生塌陷，其極限土洞規(guī)模為洞高5 m和洞寬9 m，比基于普氏平衡拱理論預(yù)測(cè)的土洞規(guī)模偏小，其誤差為15%。

(3) 覆蓋型巖溶塌陷動(dòng)態(tài)演化數(shù)值模型能夠較好地模擬土洞擴(kuò)展演化及其失穩(wěn)塌陷過(guò)程，可定量解釋覆蓋型巖溶塌陷的潛蝕—滲壓—自重致塌模式。