基于數值模擬的松散層下采煤水砂突涌預測研究

汪 鑫，王偉中，饒 傳，翟 封

（1.中煤科工集團重慶研究院有限公司，重慶 400039；2.新疆地質災害防治重點實驗室(新疆工程學院)，烏魯木齊 830000）

華東、華北許多礦區淺部煤層上覆巖層厚度在0～60 m，其覆巖具有厚度薄、強度低、風化程度高、裂隙發育等特征，導致淺部煤層開采面臨突水潰砂的威脅[1]。目前由于缺乏對相關地質條件下水砂突涌機制的研究，傳統理論方法已經很難適應對近松散層煤層開采水砂突涌發生的評價，本文采用數值模擬建立地質模型對水砂突涌進行預測[2-4]。

1 研究區概況

袁二煤礦72310 工作面72 煤層厚1.8～4.6 m，煤層傾角9°～24°，淺部煤層露頭為隱伏露頭，上覆蓋層為新生界松散層，厚度為202.5～309.4 m，平均262 m。新生界松散層分為3 個隔水層和4 個含水層，其中第三隔水層（以下簡稱“三隔”）厚度43.1～67.75 m，平均厚度51.72 m，分布穩定，隔絕了第四含水層（以下簡稱“四含”）與上覆第三含水層及以上含水層的水力聯系。

“四含”厚度8.4～9.9 m，平均9.18 m，巖土體性質為砂礫、粘土夾礫石、粗砂及粘土質砂等，其間夾有0～4 層薄層狀粘土。單位涌水量為0.002 12～0.015 14 L/(s.m)，K=0.000 7～0.11 m/d，為弱富水性含水層。

2 “四含”土體特征與滲透破壞類型

2.1 “四含”土體特征

2.1.1 “四含”巖性分析

“四含”巖性為砂礫、細砂、粘土夾礫石、粗砂及粘土質砂等組成，其間夾有0～4 層薄層狀粘土、砂質粘土、鈣質粘土。

圖1 含粘土細砂，“四含”

1）“四含”粘土為褐黃色，局部含少許粉砂且分布不均，含少量鐵錳質結核，具粘性和可塑性，半固結。

2）“四含”細砂為棕褐- 棕黃色，細粒，分選性中等，半固結，致密。

3）“四含”礫石深灰色，局部夾褐黃色，成分以灰巖為主，次為長石石英砂巖，礫徑3～8 cm，次棱角狀、次圓狀，堅硬，松散。

2.1.2 “四含”粗細顆粒含量分析及顆粒級配

“四含”細砂中細粒含量約占32%，砂粒含量約占68%；“四含”礫石中細粒含量約占17%，砂粒含量約占49%，礫粒含量約占34%。工作面“四含”細砂Cu=12.8，Cc=1.02；礫石Cu=120，Cc=2.89，故“四含”砂層屬于級配良好的土，滲透性差。

表1 “四含”土體顆粒含量與級配

2.2 “四含”土體破壞類型

根據《煤礦防治水細則解讀》：“厚度大于0.25 m的粉砂細層，砂土的孔隙率大于43%，含水率高，土的組成中粉砂含量大于75%”[7]和“四含”土體顆粒含量與級配，“四含”土的組成中粉砂含量小于75%，不滿足《煤礦防治水細則解讀》中涌水潰砂條件。

根據“四含”土體顆粒含量與級配細砂中細粒含量32%，“四含”細砂屬于不均勻不連續的過渡型土；礫石中細粒含量17%，屬于不均勻不連續管涌土。根據“四含”巖土體性質，不易發生潰砂。

3 采動裂隙平面分布及臨界水力坡度確定

3.1 采動裂隙平面分布

3.1.1 導水裂隙帶高度計算

本次理論計算導水裂隙帶公式依據為2017 年國家煤礦安全監察局《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規程》[5]和《煤礦防治水細則解讀》[7]。

導水裂隙帶計算公式：

式中：H為導水裂隙帶高度，m；M為工作面煤層采厚，m。

計算結果如下：

表2 理論計算導水裂隙帶高度

根據2 個公式預測導水裂隙帶高度為36.22～63.09 m。

3.1.2 導水裂隙帶高度數值模擬

FLAC3D數值模擬判斷導水裂隙帶高度的判據有塑性區和應力2 種。應力判據時根據該點的應力狀態來判斷其是否破壞；塑性區判據是塑性區發育高度可以反映導水裂隙帶發育高度。

1）參數選取。依據巖石力學參數數據，并根據相鄰工作面導水裂隙帶實測值進行參數校正，綜合確定工作面巖體力學參數。

表3 巖土力學參數

2）建立模型。根據工作面煤層底板等高線圖及鉆孔、勘探線，建立了三維地質模型。X、Y、Z 方向為走向、傾向、鉛直，分別為850 m×450 m×400 m，采空部分使用零單元模擬。

圖2 地質模型

3）模型運行。結合模擬過程中應力、位移等特征分析，當工作面回采130、220 m 時，導水裂隙帶波及基巖風化帶和“四含”；工作面全部回采完塑性區出現頂板破壞產生的剪應力，剪切破碎帶鉛直范圍和預測裂隙帶高度值相近，約51 m，裂采比11.09。

圖3 全回采塑性區圖

圖4 全回采應力云圖

3.1.3 采動裂隙帶分布

地表下沉盆地的形狀和位置取決于采空區的形狀和煤層傾角。在地表下沉盆地內，不同位置的移動和變形性質且大小各異，超充分采動時形成的地表下沉盆地可劃分為3 個區域：中性區域、壓縮區域、拉伸區域[8]。

針對工作面實際情況及“兩帶”發育情況繪制出工作面“兩帶”平剖面分布圖?！皟蓭А庇绕淅靺^域波及基巖風化帶與“四含”為礦井涌水與潰砂的提供通道。

圖5 工作面采動裂隙平面分布圖

3.2 臨界水力坡度確定

考慮到砂層具有一定的粘結程度，采用扎馬林公式[6]：

式中：γs為土粒容重，N/m3；Γ為水的容重，N/m3；N為孔隙率，n=e/(1+e)；E 為孔隙比，取0.6。

根據公式計算，確定“四含”潰砂臨界水力坡度Jcr=1.28。

4 數值模擬預測實際水力坡度

4.1 建立模型

4.1.1 概念模型

由于工作面的邊界是人為劃定的邊界，為非自然邊界。根據鉆探、物探、水文地質資料，斷層未對松散層“四含”造成影響，即采區邊界設置為流量邊界。模型頂部為砂質粘土類隔水層（三隔），底部邊界為煤層底板以下50 m 泥巖隔水層，在模型中均概化為隔水邊界。

根據模擬區內施工鉆孔資料（勘探線），去除或修正異常點相關數據，得到區內各模型層的底板標高離散點數據，在此基礎上采用反距離加權插值方法生成地質模型。

圖6 采區地層模型

4.1.2 數學模型

研究區域內地下水為非穩定流運動，采用三維滲流數學模型進行有限差分法求解[9]。

式中：H為水頭高度，H；K為滲透系數，m/d；S為彈性釋水率；W為補給量；t為時間變量。

4.1.3 幾何模型

根據Visual MODFLOW 的要求，采用等間距有限差分的離散方法對含水介質進行剖分，根據本區滲透介質空間分布特點，為了盡可能真實地反映巖層中地下水的滲流狀況，在滿足模擬精度的前提下，建立網格單元為50×50×4，模型單元格總數10 000。

4.2 模型校驗

模型校驗是檢驗所建計算機數值模型在模擬期受各種因素影響下地下水動力場時空分布與實測值時空分布的一致性，從而確定水文地質參數。

參數反演采用已回采的相鄰工作面放水試驗進行模型識別與調試。以本次放水試驗數據為基礎，經過反復調整與試算，使計算的水位值與實測的水位值之差最小，從而達到數值仿真的目的。模型調試后，觀1 孔地下水位與時間擬合曲線如下。

圖7 曲線反映了工作面疏放水時造成觀1 孔水位隨時間持續下降約1 m。基于上述數據導入excel軟件，生成觀1 孔水位歷史擬合曲線圖（圖8），對觀1孔在相鄰工作面回采實際水位與時間進行線性回歸，回歸方程為y = -0.0393x+2.9679，擬合系數R2=0.967 1；再對軟件模擬水位與時間進行線性回歸，回歸方程為y=-0.034 3x+2.955，擬合系數R2=0.886 8。2 條擬合曲線趨勢一致，故數值模擬值與實測值趨勢基本一致，該模型可用于模擬“四含”流場變化規律。

圖7 觀1 孔模擬與實測水位歷史曲線

圖8 觀1 水位歷史擬合曲線圖

4.3 模型運行

結合Visual MODFLOW 模擬結果和采動裂隙平面分布成果，確定導水裂隙帶影響范圍內的水力坡度。

根據圖9、圖10，通過對導水裂隙帶范圍內“四含”水位變化進行涌水潰砂判別：

表4 “四含”涌水潰砂判別依據

圖9 回采130、220 m 時導水裂隙帶范圍內水位變化

圖10 全回采時導水裂隙帶范圍內水位變化

當回采130、220 m 時，預測水力坡度為0.4、0.64，均小于臨界水力坡度；工作面482 m 全部回采時，預測水力坡度0.8，仍然小于臨界水力坡度。說明工作面回采過程中水力坡度不斷增大，直至全部回采結束達到最大值0.8，仍小于臨界水力坡度1.28，工作面回采不易潰砂。

5 結 論

1）根據扎馬林公式計算工作面臨界水力坡度為1.28；

2）根據Moflow 模擬結果，工作面回采期間實際水力坡度呈現增大的趨勢；

3）預測工作面全部回采結束后，導水裂隙帶影響范圍內實際水力坡度達到最大值為0.8 m，仍小于臨界水力坡度1.28，說明工作面回采不易潰砂。