港里鐵礦水下開采流固耦合數(shù)值模擬與風(fēng)險分析

馬海濤 李迎佳 王云海

(1.北京科技大學(xué)土木與環(huán)境工程學(xué)院，北京100083;2.中國安全生產(chǎn)科學(xué)研究院，北京100012;3.魯中礦業(yè)有限公司，山東萊蕪，271113)

據(jù)不完全統(tǒng)計，我國目前受河流影響而滯留大量礦石儲量的冶金礦山65座，有色金屬礦25座，鐵礦山最為突出，10億t以上的礦石儲量被壓埋在河流之下未能合理開發(fā)，僅邯鄲、萊蕪地區(qū)就有5億多t鐵礦石賦存于水體之下，難以開發(fā)。因而，解決水體下安全開采問題成為冶金礦山所面臨的重大生產(chǎn)安全問題之一［1］。水下開采地質(zhì)災(zāi)害本質(zhì)上是由水、巖等固體、準(zhǔn)流體及流體構(gòu)成的復(fù)雜地質(zhì)系統(tǒng)對開挖擾動作出的響應(yīng)或反饋。響應(yīng)的方式和程度不同，災(zāi)害的類型和規(guī)模也就不同。從物理形態(tài)上硬巖、軟巖、水等是可以分開的，但是，無論是天然條件下，還是開挖環(huán)境下，它們的變形、運動往往都互相聯(lián)系、彼此影響。即應(yīng)力場和滲流場之間存在耦合效應(yīng)。因此，需要對水下礦山開采風(fēng)險進行分析，并考慮流固耦合作用的影響。

有關(guān)流固耦合的計算方法大致可以分為兩類，一類是解析解，一類是數(shù)值解。數(shù)值解法通常采用的是有限單元法或有限差分法［2-4］。有限差分法FLAC3D分析三維流－固耦合問題，應(yīng)用較為廣泛。

1 三維流固耦合計算方程

相互關(guān)聯(lián)的變形－擴散過程的增量公式是通過線性準(zhǔn)靜態(tài)耦合理論的框架結(jié)構(gòu)完成的。描述多孔介質(zhì)中的滲流需要的變量包括孔隙壓力和3個特殊的釋放矢量。變量之間通過3個方程相互聯(lián)系，3個方程分別是:①流體質(zhì)量平衡方程;②滲流力學(xué)本構(gòu)定律，體積應(yīng)變對孔隙壓力的影響是通過流體本構(gòu)定律反映出來的，孔隙壓力的變化引起了力學(xué)變形的發(fā)生;③滲流邊界和初始條件，初始條件符合一個特定的壓力場，邊界條件通常是根據(jù)孔隙壓力或者根據(jù)垂直于邊界的比流量的分量來確定。

本研究數(shù)值計算考慮了4種情況:①給定孔隙壓力;②給定垂直邊界的比流量分量;③給定滲流邊界;④給定不具有滲透性區(qū)域的邊界。

2 工程概況

港里河在魯中礦業(yè)有限公司港里1號礦體上部流經(jīng)，為大汶河支流，屬源短流急的間歇性河流，全長46 km，流域面積144 km2，歷年最大洪峰流量為180.6 m3/s。

礦區(qū)內(nèi)均有第四系地層覆蓋。礦床主要含水層為中奧陶系大理巖、灰?guī)r，含水性和透水性均較弱。礦體頂板為石炭系板巖、角巖，節(jié)理裂隙發(fā)育，性脆易碎，穩(wěn)定性差;底板為中奧陶系大理巖，裂隙、溶洞一般不發(fā)育，較為穩(wěn)定。當(dāng)其含泥質(zhì)較多時，易破碎坍塌。

根據(jù)礦山采用無底柱分段崩落法多年生產(chǎn)實踐，本區(qū)地應(yīng)力較大，接近礦體附近時巷道多變形或坍塌，不經(jīng)支護幾乎無法生產(chǎn)，回采難度較大，采準(zhǔn)巷道多采用噴錨網(wǎng)支護。因此，礦床開采時第四系地下水和河水可能通過第三系黏土層中形成的導(dǎo)水裂隙帶下滲，出現(xiàn)坍塌、突水危害。另一方面，礦體上部第三系為隔水層，如果采用充填法進行開采，限制圍巖產(chǎn)生較大的變形，則不會破壞覆巖隔水層，可以保證礦山生產(chǎn)安全。

因此，本研究在礦山地質(zhì)條件分析、上覆巖層隔水性分析基礎(chǔ)上，進行三維流固耦合的數(shù)值模擬分析，評估河下礦體開采的安全性。

3 地表隔水層安全性分析

3.1 地質(zhì)條件分析

港里礦地質(zhì)條件簡單，區(qū)內(nèi)無較大斷層、節(jié)理切割隔水層，地表水和地下水聯(lián)系較少。故頂板裂隙帶的發(fā)育高度，決定了隔水層的安全性。根據(jù)國家煤炭工業(yè)局制定頒發(fā)的《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)程》，正常地質(zhì)構(gòu)造和充分采動條件下，港里礦(中硬和軟弱覆巖)裂隙帶最大高度

式中，HL為冒落帶加裂隙帶的高度，m;M為開采礦體真厚度，m;h為開采垂高，m。

計算得港里礦開采后頂板覆巖裂隙帶高度為105 m。由于港里礦采用充填法開采，在保證充填體強度和完全接頂?shù)臈l件下，冒落帶高度可以忽略，實際發(fā)生的裂隙帶高度也遠遠小于該值。港里礦頂板上覆巖層為厚150～700 m的第三系紅板巖，具有較好隔水性。因此，從地質(zhì)條件上看，井下導(dǎo)水裂隙不會發(fā)育到地表，不會與地表水系產(chǎn)生直接聯(lián)系，滿足井下安全采礦的隔水要求。

3.2 巖石物理力學(xué)性質(zhì)試驗

通過現(xiàn)場巖石取樣室內(nèi)試驗獲得的巖石物理力學(xué)參數(shù)見表1。

表1 礦巖石力學(xué)參數(shù)Table1 Rock mechanics parameters

港里礦礦體上覆巖層變形量較大，且在采動影響下，其蠕變性明顯，即巖體變形隨著時間緩慢發(fā)展，且持續(xù)時間較長。研究結(jié)果表明，軟巖中一般存在黏塑性區(qū)，且隨時間而增大，但最終趨于穩(wěn)定。通過實驗室三軸試驗結(jié)果表明，上覆軟弱巖層——紅板巖殘余強度為2.0～3.0 MPa。當(dāng)應(yīng)力達到6 MPa時為穩(wěn)定蠕變;當(dāng)大于6 MPa時轉(zhuǎn)化為加速蠕變。

采用最小二乘法、遺傳算法和遺傳規(guī)劃法分別對蠕變本構(gòu)模型進行參數(shù)辨識和模式識別。與紅板巖蠕變試驗結(jié)果進行比對，如圖1。

圖1 蠕變模型預(yù)測結(jié)果Fig.1 Forecasting result of themodel of creep behavior

如圖1可見遺傳規(guī)劃預(yù)測曲線與實驗曲線吻合程度較高，參數(shù)取值和蠕變模型較為合理。

3.3 地面河床安全性分析

采用三維概率積分法對地表移動變形進行預(yù)測，揭示開采沉陷的相應(yīng)規(guī)律，從而判定開采對河床帶來的影響。

選取大地坐標(biāo)(55 500，17 000)為原點(0，0)建立模型，其中X為正軸指向東、Y為正軸指向北，預(yù)測分析在最大限度利用礦產(chǎn)資源的情況下地表移動變形對河流的影響，即在礦體開采結(jié)束后，開采沉陷影響范圍。預(yù)測結(jié)果如圖2。

圖2 開采后地表下沉影響范圍預(yù)測(單位:mm)Fig.2 Region of influence of ground depression after m ining(mm)

由圖2經(jīng)計算，地面河床最大下沉量在8 cm范圍內(nèi)。地表最大傾斜值小于3 mm/m，最大水平變形值小于2 mm/m。滿足《采礦設(shè)計手冊(礦床開采卷下)》和《建筑物、水體、鐵路開采規(guī)程》規(guī)定的允許變形極值，即地表傾斜3 mm/m，地表曲率0.2×10－31/ m，地表水平變形2 mm/m。

4 三維流固耦合數(shù)值模擬

4.1 計算模型及邊界條件

根據(jù)礦區(qū)地質(zhì)特征，地表第四系取8 m，第三系隔水層取186 m，鐵礦厚度取平均厚度14.35 m，傾角16°。1號礦體劃分為2個階段開采，階段高度50 m，劃分為－146 m和－196 m 2個中段開采，上向進路充填采礦法方案實際開采范圍為－90～－190 m。

建立礦區(qū)三維模型，共計79 396個三維滲流計算單元，88 501個網(wǎng)格節(jié)點，如圖3所示。

圖3 數(shù)值模擬三維模型Fig.3 3D model for numerical simulation

根據(jù)礦區(qū)巖性屬性分析，采用Mohr－Coulomb屈服準(zhǔn)則判斷材料的破壞［5］:

式中，σ1、σ3分別是最大和最小主應(yīng)力;c、分別是黏結(jié)力和內(nèi)摩擦角。當(dāng)fs＞0時，材料將發(fā)生剪切破壞。在通常應(yīng)力狀態(tài)下，巖體的抗拉強度很低，因此可根據(jù)抗拉強度準(zhǔn)則判斷巖體是否產(chǎn)生拉破壞。計算中參數(shù)選取見表1。

4.2 計算結(jié)果及分析

從整體計算過程的不平衡力變化情況來看，模型在自重應(yīng)力作用下逐漸趨于收斂，最大不平衡力變化率達到10－5以下。說明模型最終穩(wěn)定，計算數(shù)據(jù)可信［6］。

4.2.1 無充填開采模擬結(jié)果

受到礦體開采擾動的影響，地表以及頂板巖層在不同程度位移作用下，會進入塑性狀態(tài)，破壞形式主要以拉伸破壞為主。破壞程度較大。

在不采用充填采礦法的情況下，對上向進路采礦進行數(shù)值模擬，底板最大下沉17 m，地表下沉7 m。

如果不采用充填法開采，地表會發(fā)生很大的塌陷，裂隙帶直接與地表水溝通，造成嚴重的安全事故。因此，港里礦開采必須采用充填法。

4.2.2 充填開采模擬結(jié)果

(1)水平地應(yīng)力場分析。水平地應(yīng)力的分布服從由上到下逐漸增加的變化趨勢，模型底部水平地應(yīng)力達到極值。隨著開挖的進行，受到開采擾動的影響，磁鐵礦內(nèi)部應(yīng)力卸載，形成應(yīng)力降低區(qū)，圍巖有應(yīng)力集中現(xiàn)象。

(2)破壞場分析。受到礦體開采擾動的影響，地表在不同程度位移作用下，會進入塑性狀態(tài)，破壞形式主要以拉伸破壞為主。地表稍微有破壞，如圖4所示。盡管合理使用進路充填采礦法可以有效控制地應(yīng)力，在應(yīng)力集中和應(yīng)力升高區(qū)，可以通過適當(dāng)提高充填體強度或開采時間調(diào)整等來緩解地壓的影響，但仍需采取河道治理的措施。

圖4 圍巖破壞區(qū)分布Fig.4 Distribution map of wall rock failure zone

圖5 垂直位移分布Fig.5 Distribution map of vertical displacement

(3)位移場分析。圍巖的位移方向主要指向采場，由于采用上向水平充填采礦法，鐵礦開采后立刻回填，圍巖產(chǎn)生的位移量較小。礦區(qū)沉降位移不明顯，受開采擾動影響，仍會形成下沉盆地(地表移動盆地)。地表整體沉降量不大，礦區(qū)稍微有些沉降，在4～8 cm范圍。礦區(qū)水平位移最大值也主要發(fā)生在頂板巖層，最大下沉量15 cm，如圖5所示。三維數(shù)值計算結(jié)果與開采地表沉陷理論分析相符，說明預(yù)測結(jié)果數(shù)據(jù)可信。

(4)滲流場分析?？紫端畨毫χ饕性诘谒南档讓樱瑯O值約為0.1 MPa，礦體圍巖孔隙水壓力很小。滲流矢量最大值為1.962×10－2m/s。但隨著充填開采的進行，孔隙水壓力和滲流矢量在不斷增加。隨著開挖的進行，受到開采擾動的影響，第四系的水會向下流。同時，頂板的水會流向開采區(qū)域，成層狀分布。

5 結(jié)論

(1)由于受開采擾動的影響，地表仍然會發(fā)生移動變形，有一定下沉量，但礦區(qū)整體沉降量不大，在4～8 cm范圍。滿足采礦擾動引起的地表最終下沉情況和《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)程》等相應(yīng)的安全規(guī)定。

(2)礦井不受地表水害威脅，可以滿足河下安全開采的要求。圍巖產(chǎn)生的位移量較小，礦體圍巖孔隙水壓力很小，隨著充填開采的進行，孔隙水壓力和滲流矢量在不斷增加。

(3)根據(jù)三維滲流數(shù)值模擬和概率積分法計算分析，港里鐵礦1號礦體的開采，造成地表平緩、均勻的下沉和水平移動，總體下沉量不大，在港里河下開采，不會造成地表河流堤壩變形破壞，不會影響河流儲水能力。

(4)采用三維流固耦合(滲流場與力場相互作用)力學(xué)數(shù)值模擬分析研究港里鐵礦1號礦體采用上向進路充填采礦法開采過程中的滲流場變化和力學(xué)穩(wěn)定性，為礦山設(shè)計和生產(chǎn)提供了參考依據(jù)。

［1］ 何國清，楊 倫，賈鳳彩.礦山開采沉降學(xué)［M］.北京:中國礦業(yè)大學(xué)出版社.1991.

He Guoqing，Yang Lun，Jia Fengcai.Mining Subsidence Science［M］.Beijing:China University of Mining and Technology Press，1991.

［2］ 朱 軍，劉光廷，藍青松，等.飽和模型非定常滲流的數(shù)值求解方法［J］.巖土力學(xué)，2002，23(3):355-357.

Zhu Jun，Liu Guangting，Lan Qingsong，et al.Numericalmethod for saturated unsteady seepage［J］.Rock and Soil Mechanics，2002，23 (3):355-357.

［3］ 王金安，馬海濤.河下金屬礦合理開采順序的固流耦合力學(xué)分析［J］.金屬礦山，2005(S):133-134.

Wang Jinan，Ma Haitao.Numerical analysis of the scheme formining metal ore under the body of water［J］.Metal Mine，2005(S):133-134.

［4］ 馬海濤.“11.6”特別重大坍塌事故礦區(qū)采場穩(wěn)定性三維數(shù)值模擬分析［J］.中國安全生產(chǎn)科學(xué)技術(shù)，2007，3(6):68-72.

Ma Haitao.3D numerical simulation of stope stability in 11.6 accident of gypsum mine collapse［J］.Journal of Safety Science and Technology，2007，3(6):68-72.

［5］ 蔡美峰，何滿潮，劉東燕.巖石力學(xué)與工程［M］.北京:科學(xué)出版社，2002.

CaiMeifeng，He Manchao，Liu Dongyan.Rock Mechanics and Engineering［M］.Beijing:Science Press，2002.

［6］ 王金安.巖土工程數(shù)值計算方法實用教程［M］.北京:科學(xué)出版社，2010.

Wang Jinan.Practical Tutorial of Numerical Calculation Method on Geotechnical Engineering［M］.Beijing:Science Press，2010.