張廣超，陶廣哲，曹志國，趙勇強(qiáng)，徐祝賀，李 友，左 昊

（1.山東科技大學(xué) 能源與礦業(yè)工程學(xué)院，山東 青島 266590；2.煤炭開采水資源保護(hù)與利用國家重點實驗室，北京 102209；3.北京低碳清潔能源研究院，北京 102211）

0 引 言

近年來，隨著我國煤炭開發(fā)戰(zhàn)略西移，煤炭大規(guī)模開采與水資源保護(hù)間的矛盾極為突出。 為此，國家能源投資集團(tuán)有限責(zé)任公司提出煤礦地下水庫技術(shù)體系［1－2］，利用煤炭開采形成的采空區(qū)巖體空隙儲水，為我國西部煤炭開采與水資源保護(hù)利用協(xié)調(diào)提供重要技術(shù)途徑［3－4］。 在地下水庫建設(shè)過程中，煤炭開采覆巖破壞形態(tài)及發(fā)育高度直接決定著地下水庫的儲水庫容能力［5－7］。 因此，研究并揭示煤炭開采后覆巖破壞形態(tài)對于地下水庫高效建設(shè)具有重要意義。

國內(nèi)外科技工作者對地下水庫建設(shè)中采場覆巖破壞規(guī)律進(jìn)行了大量研究。 如曹志國等［8］提出了采動覆巖導(dǎo)水裂隙的主通道分布模型，并闡明了水體流動力學(xué)特性；李全生等［9］分別基于單一煤層開采與煤層群開采條件下的導(dǎo)水?dāng)嗔褞Ц叨阮A(yù)計方法，對地下水庫保水技術(shù)的適應(yīng)性進(jìn)行了評價；鞠金峰等［10］通過對導(dǎo)水裂隙帶內(nèi)空隙的定量計算，確定了地下水庫極限庫容的計算方法；上述研究多采用理論分析、現(xiàn)場原位觀測和室內(nèi)實驗的方法對導(dǎo)水?dāng)嗔褞Оl(fā)育高度或形態(tài)進(jìn)行分析預(yù)測及基于上述探測結(jié)果的地下水庫庫容計算。 事實上，采場覆巖破壞規(guī)律受諸多因素影響，理論計算方法多基于大量假設(shè)，而現(xiàn)場實測方法費(fèi)時費(fèi)力且實測結(jié)果受客觀因素影響，因此，數(shù)值模擬成為分析采場覆巖破壞特征的簡單、有效方法。 但在以往的數(shù)值模擬研究中，忽略了實際工程中冒落矸石漸進(jìn)壓實過程而誘發(fā)的應(yīng)變硬化行為［11－15］，因而無法科學(xué)地指導(dǎo)地下水庫庫容計算。

基于李家壕煤礦31108 工作面地質(zhì)生產(chǎn)條件，探討采空區(qū)垮落巖體壓實過程中的應(yīng)變硬化力學(xué)特性，并基于FLAC3D內(nèi)置雙屈服本構(gòu)模型尋找其在數(shù)值模擬中的仿真實現(xiàn)方法，模擬分析采場覆巖破壞形態(tài)及應(yīng)力位移場特征，為地下水庫設(shè)計提供理論支撐。

1 工程地質(zhì)概況

李家壕煤礦位于內(nèi)蒙古鄂爾多斯市東勝煤田中南部，設(shè)計生產(chǎn)能力600 萬t／a。 當(dāng)前該礦主采2－2煤和3－1 煤。 2－2 煤厚度1.4～5.3 m，平均厚度2.03 m，大部分可采，煤層結(jié)構(gòu)簡單，頂板以粉砂巖與細(xì)粒砂巖為主，底板以砂質(zhì)泥巖為主；3－1 煤厚度2.5 ～6.3 m，平均4.0 m，煤層結(jié)構(gòu)簡單，頂板以砂質(zhì)泥巖和粉砂巖為主，底板主要是砂質(zhì)泥巖。 3－1 煤埋深204～254 m，地層中主要含水層為基巖孔隙裂隙含水層以及第四系松散含水層。 試驗工作面為李家壕煤礦31108 工作面，主采3－1 煤層，工作面內(nèi)TL32鉆孔揭示煤層柱狀如圖1 所示。

圖1 31108 工作面鉆孔柱狀Fig.1 Borehole column of panel No.31108

31108 工作面頂板砂質(zhì)泥巖的孔隙率為4.54%～30.00%，含水率為0.13%～5.42%，吸水率為1.08% ～12.28%，抗壓強(qiáng)度為2.4 ～80 MPa，平均41 MPa，抗拉強(qiáng)度為0.41 ～3.24 MPa，黏聚力為1.2 MPa，內(nèi)摩擦角為20°，泊松比為0.3，軟化系數(shù)為0.12～0.73，巖體碎脹系數(shù)取1.3。

2 垮落巖體應(yīng)變硬化力學(xué)行為仿真

2.1 垮落巖體應(yīng)變硬化力學(xué)特性

在煤炭開采過程中，巖層自下而上運(yùn)移破壞，形成垮落帶、斷裂帶和彎曲下沉帶。 破碎巖體由松散體向承載體轉(zhuǎn)化的力學(xué)過程，即應(yīng)變硬化行為，將對覆巖變形應(yīng)力分布發(fā)生顯著影響［16－18］。 Salamon 將冒落塊體作為顆粒物質(zhì)，提出的垮落巖體應(yīng)力－應(yīng)變關(guān)系被廣泛采用，其表達(dá)式如下：

式中：σcap為垮落巖體受到的垂直載荷；ε為在應(yīng)力σcap作用下垮落巖體的體積應(yīng)變；εmax為產(chǎn)生的最大體積應(yīng)變；E0為垮落巖體的初始彈性模量。εmax和E0的取值取決于垮落巖體的碎脹系數(shù)與強(qiáng)度：

式中：Kp為垮落巖體的碎脹系數(shù)；σc為巖體的抗壓強(qiáng)度。 將式（2）、式（3）代入（1）可得，垮落巖體應(yīng)變硬化力學(xué)特性如下：

根據(jù)式（4）可以推導(dǎo)得出采空區(qū)垮落巖體應(yīng)力狀態(tài)與碎脹系數(shù)、巖石強(qiáng)度的關(guān)系，如圖2 所示。

由圖2a 可知，隨著應(yīng)變ε 的增加，垮落巖體的應(yīng)力先緩慢增大，當(dāng)達(dá)到某一臨界應(yīng)變后，應(yīng)力呈近似指數(shù)型快速增長；且隨著碎脹系數(shù)的增加，應(yīng)力呈現(xiàn)指數(shù)型增長所對應(yīng)的臨界應(yīng)變亦隨之增大。由圖2b 可知，隨著應(yīng)變增加，應(yīng)力先緩慢增加后呈快速增長趨勢，且?guī)r塊強(qiáng)度越高，最終應(yīng)力恢復(fù)值越高。

圖2 碎脹系數(shù)、巖體強(qiáng)度對垮落巖體力學(xué)特性的影響Fig.2 Effect of bulking coefficient and rock mass strength on mechanical property of caved rock masses

2.2 碎脹巖體力學(xué)特性在數(shù)值模擬中的實現(xiàn)

FLAC3D內(nèi)置雙屈服模型能較準(zhǔn)確地描述垮落巖體壓實過程中的應(yīng)力恢復(fù)行為，并被國內(nèi)外專家學(xué)者使用［19－20］。 在雙屈服模型使用過程中，需要蓋帽壓力和材料特性兩大類參數(shù)，其中，材料參數(shù)包括密度、體積模量、剪切模量、內(nèi)摩擦角與剪脹角。 蓋帽壓力可通過Salamon 經(jīng)驗公式確定，材料參數(shù)可采用反演－試錯方法確定。

3 導(dǎo)水?dāng)嗔褞Оl(fā)育規(guī)律模擬分析

為驗證上述仿真計算方法的可行性，以李家壕煤礦31108 工作面地質(zhì)條件為例，對比分析采場覆巖位移場和破壞場，確定覆巖破壞狀態(tài)。

3.1 數(shù)值模型建立

根據(jù)31108 工作面地質(zhì)生產(chǎn)條件，建立數(shù)值模型如圖3 所示。

圖3 三維計算模型Fig.3 Three dimension simulation model

模型尺寸為500 m×250 m×260 m。 模型X方向和Y方向邊界施加水平位移約束，模型底部邊界速度限定為0。 沿X、Y方向施加水平應(yīng)力，側(cè)壓系數(shù)設(shè)定為1.2。 Mohr－coulomb 模型用于模擬頂?shù)装鍘r層，其參數(shù)在室內(nèi)實驗基礎(chǔ)上通過Hoek－Brown強(qiáng)度準(zhǔn)則得出；雙屈服模型用于模擬冒落矸石。

數(shù)值模擬過程為：初始應(yīng)力計算→31108 工作面分步開挖（50、100、150、200 m）→采空區(qū)雙屈服模型定義（方案一）、未充填采空區(qū)（方案二）→運(yùn)算平衡。

3.2 垮落巖體力學(xué)參數(shù)確定

根據(jù)31108 工作面地質(zhì)條件，垮落矸石的碎脹系數(shù)取1.3，單軸抗壓強(qiáng)度取41 MPa，代入式（2）—式（3），得到最大應(yīng)變εmax為0.23，彈性模量E0為66；代入式（4），得到垮落巖體應(yīng)力分布如圖2 所示。 通過建立1 m×1 m×1 m 單元體，模型上部施加1×10－5m／s 速度，模型四周及底部位移限制為0，獲取模型平衡運(yùn)算過程中的應(yīng)力－應(yīng)變曲線。 通過調(diào)試材料參數(shù)并擬合，最終獲得的材料參數(shù)見表1。

圖4 冒落巖體應(yīng)力－應(yīng)變關(guān)系Fig.4 Stress－strain relationships of caved rock masses

表1 雙屈服模型中的材料特性Table 1 Materials’ parameters for double－yield model

3.3 模擬結(jié)果分析

3.3.1 位移分布特征

31108 工作面推進(jìn)至不同位置時位移變化曲線如圖5 所示。

圖5a 為工作面正常開挖時頂板位移變化曲線。由圖可知：①頂板下沉位移曲線呈近似倒梯形分布形態(tài)，工作面開挖區(qū)域頂板下沉量最大，未開區(qū)域頂板下沉量逐漸降低并趨于0。 ②隨著工作面開挖尺寸增大，頂板下沉量亦逐漸增大，工作面推進(jìn)50、100、150、200 m 時，頂板最大下沉量依次為134、212、270、320 mm，增加了1.38 倍。

圖5b 為定義雙屈服模型后頂板位移變化曲線。 由于采空區(qū)應(yīng)變硬化模型作用，采空區(qū)頂?shù)装暹吔缬蔁o約束轉(zhuǎn)化為應(yīng)力約束，從而極大限制了頂板下沉，見表2。 工作面推進(jìn)至不同位置時，頂板垂直位移分別為124、195、237、256 mm，相比未充填采空區(qū)時的頂板最大位移依次降低了10、17、33、64 mm。

圖5 采空區(qū)頂板位移變化曲線Fig.5 Displacement curve of roof on goaf

表2 不同推進(jìn)距離時采空區(qū)頂板位移對比Table 2 Comparison of roof displacement with various retreating distance

3.3.2 塑性破壞特征分析

圖6 為工作面正常開挖時頂板塑性破壞特征。 當(dāng)工作面推進(jìn)50 m 時，頂板局部出現(xiàn)垮落，根據(jù)拉剪破壞單元分布，可以確定垮落帶高度約為8.5 m，導(dǎo)水裂隙帶高度約為34 m；當(dāng)工作面推進(jìn)到100 m 時，垮落帶和斷裂帶高度均顯著增大，“兩帶”高度分別為15 m 和50 m；隨著工作面繼續(xù)推進(jìn)至150 m 和200 m 時，垮落帶橫向范圍增大但整體高度基本穩(wěn)定，最終垮落帶高度為18 m，斷裂帶繼續(xù)向上發(fā)育，高度依次增長為66 m 和75 m。

圖6 不同推進(jìn)距離時塑性區(qū)分布特征Fig.6 Distribution characteristics of the plastic zone with various retreating distance

圖7 為定義采空區(qū)雙屈服模型時頂板塑性破壞特征。 當(dāng)考慮落矸石碎脹特性時，塑性區(qū)破壞形態(tài)變化較小，見表3。 當(dāng)考慮冒落巖體碎脹特性時，垮落帶高度最終確定為15 m，減少約16%；斷裂帶高度確定為65 m，減少約13.3%。 在31108 工作面現(xiàn)場探測實踐中，觀測結(jié)果表明“兩帶”發(fā)育高度為57.7～65.9 m，而當(dāng)考慮垮落巖體碎脹特性時，測得的導(dǎo)水裂隙帶高度約為65 m，與現(xiàn)場實測結(jié)果較好地吻合。

圖7 考慮垮落巖體力學(xué)特性時塑性區(qū)分布特征Fig.7 Distribution characteristics of plastic zone in consideration of mechanical properties of caved rock masses

表3 不同推進(jìn)距離時“兩帶”發(fā)育高度對比Table 3 Comparison of development height of two zones with various retreating distance 單位：m

綜上分析，在模擬分析過程中，考慮冒落巖體碎脹力學(xué)特性時，由于冒落矸石對采空區(qū)的完全充實、壓實特性，頂板巖層的運(yùn)動范圍和程度明顯減低，得到結(jié)果更加接近現(xiàn)場實際。 因此，采用垮落帶力學(xué)特性及雙屈服模型，實現(xiàn)了對垮落巖體碎脹特性的仿真模擬。

4 結(jié) 論

1）隨著應(yīng)變的增加，垮落巖體的應(yīng)力先緩慢增大，當(dāng)達(dá)到某一臨界應(yīng)變后，應(yīng)力呈近似指數(shù)型快速增長；且隨著碎脹系數(shù)的增加，應(yīng)力呈現(xiàn)指數(shù)型增長所對應(yīng)的臨界應(yīng)變亦隨之增大。

2）基于FLAC3D內(nèi)置的雙屈服模型，提出了采空區(qū)垮落巖體應(yīng)變硬化特性的仿真計算方法，并提出了蓋帽壓力與材料特性的反演計算方法。

3）通過數(shù)值模擬對比分析可知，當(dāng)考慮冒落巖體碎脹力學(xué)特性時，由于冒落矸石對采空區(qū)的完全充實、壓實特性，頂板巖層的運(yùn)動范圍和程度明顯降低，得到結(jié)果更加接近現(xiàn)場實際。