考慮圍巖應變軟化及采空區接觸的深部煤層底板破壞分析

尹尚先，孟浩鵬，錢雙彬

(1.華北科技學院 安全工程學院，北京 101601；2.華北科技學院 理學院，北京 101601)

隨著時代發展，數值方法在巖土工程領域的應用逐漸普及，其中FLAC3D有限差分數值方法優點眾多，在采礦領域得到廣泛認可與應用，眾多學者借助FLAC3D圍繞煤炭安全開采進行了大量研究。董書寧等[1]在改造奧灰頂部巖層段的判別準則研究中，利用數值計算分析了采深、采高和采寬等尺寸效應對底板破壞帶的影響；劉偉韜等[2]通過正交試驗設計，進行了7 個主控因素對底板破壞深度的影響研究，并對主控因素進行排序；劉新民[3]采用現場和模擬結合的方法進行沿空留巷對底板破壞深度的研究，認為無煤柱式的沿空留巷開采技術不會對底板破壞深度造成較大影響；朱斯陶[4]、朱廣安[5]、田雨桐[6]等的研究表明，數值模擬是研究采動影響下斷層活化規律的重要手段之一，可實現采動影響下對斷層多方位的定量分析；朱慶偉[7]、甘智慧[8]等運用數值方法對采動影響下覆巖結構演化和地面沉降進行了研究。上述研究成果在數值計算中一般采用摩爾-庫倫本構模型，該本構模型為理想彈塑性模型，對圍巖塑性屈服后的狀態無法準確描述，并且對拉格朗日法遵循連續介質假設而導致網格發生大變形但節點不接觸的固有缺陷[9]未進行深入研究，這2 個問題可能拉大數值模擬同實際情況之間的差距。大變形條件下不接觸的固有缺陷，導致無法模擬采空區頂板垮落后頂底板接觸的應力傳遞現象，即“采空區不接觸”現象。觀察模擬工作面回采結果發現：頂底板接觸后并不進行接觸計算；頂底板應力基本處于泄壓狀態。真實情況是采空區在周期來壓過程中，基本頂發生周期性垮落與底板接觸，采空區的應力恢復隨時間變化呈指數函數關系[10]。因此，利用FLAC3D研究煤層回采后底板破壞時，應對圍巖應變軟化和采空區接觸進行考慮。

筆者將針對FLAC3D模擬工作面回采中本構模型的選擇和采空區不接觸的固有缺陷進行研究，以河北開平煤田林西礦2023 工作面底板導水裂隙帶實測為工程背景，建立考慮應變軟化和采空區接觸的工作面回采模型，并結合力學分析對模擬結果進行解釋，達到提高工作面回采數值計算準確性的目的。

1 工程背景及底板破壞深度實測

1.1 礦區概況

林西礦位于河北開平煤田東南翼，地質構造以褶皺為主，開平主向斜穿過井田東側深部，西部有杜軍莊背斜和黑鴨子向斜；井田內大型斷層較少，小型斷層較發育；井田內揭露地層由老到新為：奧陶系、石炭系、二疊系和第四系。石炭系、二疊系為含煤地層，煤系基底為奧陶系灰巖(簡稱奧灰)。林西礦目前主要生產水平為11 水平(?850 m)和12 水平(?1 000 m)。

12 煤位于二疊系下統趙各莊組，煤層底板至奧陶系自上而下為：石炭系上統開平組，上以K6 灰巖頂界面與二疊系下統趙各莊組分界，下以唐山組K3 灰巖頂界面與開平組分界；石炭系上統唐山組，上以K3 灰巖頂界面與開平組分界，下以G 層鋁鐵質泥巖的底界面與奧陶系石灰巖呈平行不整合接觸。趙各莊組為主要含煤組，厚度33.55～61.20 m，平均48.35 m，含11、12 兩層可采煤層；本組巖性頂部主要為黑色泥巖和灰色砂巖及褐灰色砂巖，其中砂巖向深部逐漸發展為黑色泥巖；中部及底部主要為灰色細?粗砂巖，淺部間有3～4 層砂礫巖，向深部礫巖直徑逐漸變小，砂巖粒度也逐漸變細且大部分砂巖為泥巖所替代。開平組層厚55.09～91.04 m，平均76.81 m，主要巖性為黑色粉砂質泥巖，砂巖次之，其中砂巖比例由淺向深逐漸減少。唐山組層厚64.44～80.55 m，平均70.18 m；地層巖性除K1、K2、K3 灰巖和G 層鋁鐵質泥巖外，主要為黑色-深灰色泥巖和灰色砂巖。

根據林西礦深部ZK7 奧灰水位觀測孔鉆孔水位及12 煤層底板至奧灰頂界面間距計算得到：當12 煤層底板高程小于?936.3 m 時，工作面回采期間突水系數超過0.06 MPa/m。按照《煤礦防治水細則》，如果突水系數超限問題不解決，深部區域將無法進行安全帶壓回采，嚴重影響礦井采掘接替，威脅礦井生存。因此，有必要對12 煤層底板采動導水破壞帶深度進行實測及底板破壞分析等工作。

1.2 工作面概況

林西礦深部2023 工作面開采12 煤層，位于林西井田杜軍莊背斜構造塊內，煤層走向變化較大(N11°EN36°E)，煤層厚度0.8～2.7 m，平均2.0 m；煤層傾角17°～21°，平均20°；走向長約661 m，傾斜長約93 m；地面高程30 m，開采高程?842.8～?884.0 m。直接頂為炭質泥巖，厚約3.13 m；基本頂為泥巖，厚約2.08 m；直接底為泥巖，厚約0.7 m，老底為粉砂巖，厚0.5～4.6 m。

1.3 底板破壞深度實測

為觀測林西礦2023 工作面底板采動導水破壞帶發育情況，在2023 工作面西側2023-2 工作面回風巷設計D01-1、D01-2 鉆孔，鉆孔鉆至2023 工作面法向向下30 m 左右位置，如圖1 所示，隨后進行分段壓水試驗，記錄不同位置流量穩定后的漏失量，繪制鉆孔的漏失量曲線(圖2)。觀察鉆孔的漏失量曲線發現，在垂距24 m 以后，漏失量下降明顯，注水壓力在漏失量穩定后回彈明顯且接近初始注水壓力，說明該段裂隙不發育且貫通性差。在壓水試驗基礎上，為直觀了解底板破壞情況以及對壓水試驗結果進行驗證，對D01-2鉆孔進行鉆孔成像，獲取成像數據如圖3 所示，圖中距離均已換算為距煤層底板垂深。

圖1 2023 工作面底板破壞深度實測鉆孔布置Fig.1 Layout of boreholes for the failure depth in the floor of working face 2023

圖2 2 個鉆孔漏失量曲線Fig.2 Two boreholes leakage curves

圖3 D01-2 鉆孔不同深度成像Fig.3 D01-2 borehole imaging at different depths

根據D01-2 鉆孔成像結果顯示，圖3a 中14.2 m處巖性為灰色細砂巖，處于底板破壞帶邊緣，該處巖層完整性較好，但存細微層狀裂隙；圖3b、圖3c 中15.9～19.3 m 處，巖性為細砂巖，這2 處巖層及其之間的巖石破碎嚴重，縱橫裂隙發育明顯；圖3d 中26.6 m處，巖性為細砂巖，該處巖層及以下巖層完整，未見明顯裂隙。認為D01-2 鉆孔漏失量曲線在30 m 處的回彈是由于壓水試驗中的封堵裝置壓力過大致使原生裂隙張開導致。

綜合D01-1、D01-2 孔壓水試驗實測結果和鉆孔成像數據，并參考鄰近趙各莊礦近似埋深的1237、2137 工作面正常底板實測采動導水破壞帶深度23、25 m[11]，最終確定林西礦2023 工作面正常底板采動導水破壞帶深度為24 m。

2 底板破壞力學分析

運用朗肯土壓力理論對煤層底板破壞進行定性分析[12-13]，圖4 表示具有半無限平面的煤層開采走向剖面圖，將底板塑性破壞劃分為Ⅰ區(主動區)、Ⅱ區(過渡區)、Ⅲ區(被動區)3 個區域。圖4 表示回采過程中煤層底板受力狀態的應力圓與底板強度包線之間的關系。

圖4 煤層底板破壞分區Fig.4 Coal seam floor failure zones

在原巖應力狀態下，煤層下深度為z處單元體的應力為豎向應力σv(σz)、最大水平應力σH、最小水平應力 σh，取σH=σ1(最大主應力)、σv=σ3(最小主應力)的情況，用圖5 中的應力圓①表示，此狀態下應力圓距強度包絡線較遠，底板處于彈性平衡。隨工作面的推進，底板地應力受工程擾動，在超前應力作用下 σv將逐漸增大，假設σH保持不變，此時應力圓的半徑先減小后增大，如圖5 中方向向右的箭頭所示，若剪應力達到底板抗剪強度，應力圓與強度包絡線相切，該處底板達到被動極限平衡狀態，如圖5 中應力圓②所示，此時σH=σ3、σv=σv1=σ1(σv1為應力圓②狀態下最大主應力)，根據摩爾?庫倫理論可知，當工作面前方底板達到或超過極限平衡狀態并隨應力的持續作用，Ⅰ區發生塑性變形，伴隨體積膨脹以壓力形式通過Ⅱ區向采空區方向(Ⅲ區)傳遞，導致底鼓同時形成連續的滑移面。當工作面推過之前超前應力作用的底板位置后，最大、最小主應力發生轉變，由于工作面的推進，σv減小，假設σH保持不變，此時應力圓的半徑先減小后增大，如圖5 中方向向左的箭頭所示，直至應力圓與強度包絡線相切，該處底板達到主動極限平衡狀態，如圖5中應力圓③所示，此時σH=σ1、σv=σv2=σ3(σv2為應力圓③狀態下最小主應力)，根據摩爾-庫倫理論可知，當工作面前方底板達到或超過極限平衡狀態并隨應力的持續作用，Ⅲ區發生塑性變形，導致采空區底板隆起，同時形成連續的滑移面。

圖5 煤層底板極限平衡狀態Fig.5 Limit equilibrium state of coal seam floor

為比較底板的被動極限平衡狀態和主動極限平衡狀態，根據極限應力圓與強度包絡線之間所得的關系式[14]：

式中：φ為巖石內摩擦角；c為巖石純剪切強度(黏聚力)。

將應力圓②和③的表達式代入式(1)可得底板被動和主動極限平衡狀態之間的關系：

需注意：①實際煤層開采工作面斜長有限，同時工作面斜長與底板采動破壞帶深度密切相關，導致對底板破壞進行定量分析存在困難，但采用半無限平面對底板破壞進行定性力學分析是可行的；② 對煤層下深度為z處單元體最大水平應力 σH保持不變的假設同實際不符，Ⅰ區的 σH應當隨頂板周期來壓而發生周期性變化，Ⅰ區在豎向應力作用下，應力必然向四周傳遞，由于采空區這一臨空面的存在，Ⅱ區、Ⅲ區的 σH隨Ⅰ區的 σH的變化而變化；③從回采過程中底板受力狀態的應力圓與底板強度包絡線之間的關系可知，達到應力圓②的狀態比應力圓③要困難，即Ⅰ區達到極限平衡狀態較難。

3 數值模擬論證

3.1 本構模型的選擇

底板破壞模擬一般采用摩爾-庫倫本構模型，此模型為理想狀態的彈塑性模型，不考慮黏聚力、內摩擦角、抗拉強度等材料參數隨塑性變形的變化情況[15]，同實際地質材料受力特征不符。此模型在底板塑性屈服前能較好地反映底板變形，但在底板發生塑性屈服階段開始后，同實際破壞有較大差異。底板塑性屈服后，在應力作用下呈應變軟化行為，底板在地應力的作用下產生微裂紋及巖體的相對滑動，底板強度將不斷降低并且越來越缺乏彈性，直到破壞以及剪切帶的產生。利用FLAC3D建立應變軟化摩爾-庫倫地層模型，能夠對圍巖塑性破壞后的力學狀態更準確表述。

實際巖石峰后的應變軟化過程的彈塑性剛度矩陣為一個不定矩陣[16]，導致應變軟化問題求解困難。為避免這一情況，可將巖石峰后應變軟化過程簡化為一系列的脆塑性過程[17-19]。FLAC3D內置的應變軟化模型為基于經典彈塑性理論將實際的應變軟化過程中黏聚力、內摩擦角、剪脹角與塑性剪切應變的函數近似為一組首尾相連的分段線性函數的模型。

基于巖石軟化相關文獻[20-21]，將其實驗數據擬合為巖石峰后黏聚力、峰后內摩擦角隨塑性剪切應變εs的指數函數：

為驗證擬合函數的合理性，將函數嵌入應變軟化本構模型，對參考文獻中的泥巖進行單軸壓縮數值模擬。模擬采用單軸壓縮試驗常用直徑(D)∶高(H)為1∶2 的圓柱體進行試驗，為較好地模擬真實試驗和呈現應力-應變全過程曲線，對試件端部施加恒定速度代替試件受壓情況。根據上述條件進行2 種本構模型的試驗，得到σ-ε曲線(圖6)。屈服前，2 種材料的σ-ε曲線一致且基本符合線彈性；屈服后，摩爾-庫倫材料與應變軟化材料的σ-ε曲線明顯不同。應變軟化材料的σ-ε曲線同真實試驗有較高的吻合度，可較好地描述巖石屈服軟化后的力學特征，說明根據式(3)、式(4)建立的應變軟化模型是合理可行的。

圖6 2 種本構模型應力-應變曲線Fig.6 Stress-strain curves of two constitutive models

3.2 模擬過程

3.2.1 初始模型建立

數值模擬的工況條件以林西礦12 煤2023 工作面為背景，由于12 煤層厚度0.41～8.48 m，煤層含夾矸0～1 層，夾矸厚度0.10～0.31 m，結構較簡單，埋深可至1 000 m 以下，故模擬工作面采高4 m，傾向與走向長度100 m×800 m，埋深1 000 m。地層信息參考林西礦深部ZK7 奧灰水位觀測孔鉆孔信息，對地層傾角、巖石力學參數相近及薄巖層進行適當簡化，最終模型尺寸長×寬×高為1 000 m×300 m×240 m，剖分網格數量804 000 個，節點數量833 748 個(圖7)，數值模型的巖石力學參數見表1。

圖7 煤層底板破壞數值模型Fig.7 Numerical model of coal seam floor failure

表1 數值模型巖石力學參數Table 1 Rock mechanical parameters of the numerical model

本構模型選用本文提出的應變軟化模型，對模型施加10 m/s2豎直向下的重力加速度，在模型頂部施加22.618 MPa 應力代替未建模的上覆巖層，底部采用豎直位移約束，四周采用水平位移約束，并在四周施加隨深度增加的側向應力。側向應力的大小參考沉積巖的應力分布規律[22]，根據模型煤層埋深1 000 m，得垂直應力(σv)∶最大水平主應力(σH)∶最小水平主應力(σh)為1∶1.133∶0.758，為研究原巖應力對結果的影響，在模擬中施加2 種相反地應力分布，設置初始地應力分布情況：σZ∶σX∶σY=1∶1.133∶0.758、σZ∶σX∶σY=1∶0.758∶1.133(X、Y、Z為模型坐標方位)。

3.2.2 考慮采空區接觸的方法

實際煤層開采中，基本頂隨周期來壓垮落，垮落后的頂板與底板接觸，發生頂底板之間的應力傳遞。但在FLAC3D模擬工作面回采中發現，頂底板在應力作用下接觸后，并不會進行應力接觸計算，在大變形模式下，甚至可清楚觀察到頂底板發生交叉的現象，這與實際情況不符，導致模擬與實際產生巨大偏差，所以在模擬工作面回采中需要對采空區進行接觸模擬。為體現真實回采中應力變化，開挖步距參考真實工作面周期來壓步距，模擬過程中工作面以20 m 為步距循環開挖。煤層回采后，回采區域由“應變軟化”模型轉為“空”模型。將實際采空區頂板垮落后堆積的碎石假設為彈性整體，利用“彈性(各向同性)”模型替換“空”模型，達到模擬采空區頂板垮落后頂底板應力接觸的目的。在彈性體參數的確定上，由于采空區的碎石是頂板垮落產生，將其視為裂隙發育的彈性整體，其彈性模量將大幅衰減，泊松比有所上升[23]。

數值計算時，工作面回采3 個循環后，模擬頂板垮落后的頂底板接觸，即“空”模型以20 m 為步距循環向“彈性”模型轉化，為避免彈性體對側向產生應力傳遞，在切眼、側幫及終采線附近不改變“空”模型。在工作面中心頂底板位置分別布置測點，記錄回采全過程豎直方向應力值和豎直位移量(圖8、圖9)。根據測點記錄發現，模擬采空區頂板垮落后的接觸與對采空區不做處理的結果具有顯著差異，考慮采空區接觸的結果更貼近實際。考慮采空區接觸的采空區頂底板應力得到一定程度的恢復，而不考慮接觸的采空區頂底板應力基本處于泄壓狀態；考慮采空區接觸的采空區頂底板豎直位移量明顯小于不考慮接觸的情況，并且考慮采空區接觸后的底板位移量表現出小幅的回落。

圖8 煤層頂底板豎直方向應力Fig.8 Stress in vertical direction of coal seam roof and floor

圖9 煤層頂底板豎直位移Fig.9 Vertical displacement of coal seam roof and floor

采用“應變軟化-空-彈性”模型轉變的方法，達到模擬采空區頂板垮落后應力傳遞的效果，彌補了以往煤層開采模擬中采空區頂底板不接觸的固有缺陷。

3.3 模擬結果及分析

運用自定義應變軟化本構關系和考慮采空區接觸的數值方法進行目標工作面的回采模擬，對地應力σZ∶σX∶σY=1∶0.758∶1.133 模擬過程中(回采40、100、800 m 平衡后)的塑性區(圖10)和累計塑性剪切應變率大于0.01 的區域(圖11)進行切片展示。

圖10 不同回采距離塑性區分布Fig.10 Plastic zone distribution in different mining distance

圖11 不同回采距離塑性剪切應變突出區Fig.11 Strain-shear-plastic outburst zone in different mining distance

1) 底板塑性區分析

模擬回采0～100 m 過程中，塑性區深度迅速增大，后隨回采的進行，塑性區深度基本穩定在23 m 左右，同實際底板導水裂隙帶深度的實測結果一致；底板塑性區上部的狀態為“shear-p、tension-p”，為剪切屈服和張拉屈服共存狀態，分布形態隨回采呈周期性分布；底板塑性區下部的狀態基本為“shear-p”，為剪切屈服狀態，將上部“shear-p、tension-p”狀態包圍。通過結合該區域應力及位移分布，底板塑性區上部的最小主應力呈拉應力，下部呈壓應力；在底板2 種塑性狀態交界附近的位移量有一定突變。結果表明，底板破壞區域上部為剪切、拉張交互破壞，下部為剪切破壞。

2) 塑性剪切區分析

底板塑性區全區包含塑性剪切狀態(圖10)，選取累計塑性剪切應變率大于0.01 的區域進行顯示(圖11)，其分布形態為斜向采空區的半包圍面狀結構。

截取2 種極限地應力條件下的3 個開采循環距離進行底板破壞分析，對塑性區、塑性剪切應變率大于0.01 進行整合處理(圖12)。發現原巖應力的改變對底板破壞規律幾乎沒有影響，塑性剪切應變突出區域將剪切、拉張交互破壞區(“shear-p、tension-p”)同剪切破壞區(“shear-p”)劃分開，即剪切帶內側為剪切、拉伸交互破壞，外側為剪切破壞。一般認為，巖石沿最大有效剪應力面形成剪切破裂面，而塑性剪切應變集中區域與最大有效剪應力集中區域基本一致，故認為滑移面即剪切破壞面是沿塑性剪切應變集中區域分布。通過利用優化后的數值方法得到的結果可知，煤層回采后底板破壞類型可分為剪切和拉張交互破壞、剪切破壞2 種類型，并可根據塑性剪切應變集中度，對底板滑移線進行三維可視化顯示，底板滑移面呈斜向采空區的半包圍面狀陣列分布。

圖12 2 種地應力下局部模擬結果Fig.12 Local simulation results under two ground stresses

將數值結果與本文提到的底板破壞力學分析結合，工作面回采前，底板巖體處于地應力平衡狀態，受采動影響底板初始應力狀態被打破，在采空區前方底板(Ⅰ區)產生超前支撐壓力，由于采空區這一臨空面的存在，Ⅰ區的大部分應力和位移通過Ⅱ區向Ⅲ區傳遞，致使采空區底板整體處于高圍壓、低軸壓狀態，底板破壞形式整體表現為塑性剪切破壞。對塑性剪切應變率較高區域進行顯示，該區域將底板塑性破壞區分為上下部分。朗肯土壓力理論中3 個破壞區是根據塑性剪切程度進行劃分，認為模擬結果中底板塑性區可根據塑性剪切應變率的大小進行劃分，劃分結果可與朗肯土壓力理論中3 個區進行對應。本文底板塑性區中塑性剪切應變率突出的區域(滑移面)為Ⅱ區和Ⅲ區之間的分界，Ⅰ區未顯現，同前文“Ⅰ區達到極限平衡狀態較難”對應。Ⅱ區的破壞形式為剪切破壞，Ⅲ區的破壞形式為拉張和剪切的交互破壞。

4 結 論

a.運用朗肯土壓力理論并結合考慮圍巖應變軟化和采空區接觸的數值方法研究了河北林西礦深部底板破壞特征，根據塑性剪切應變率的變化，對底板滑移面實現了三維顯示，并將底板塑性區與朗肯土壓力中的主動區、過渡區和被動區相對應，其中過渡區、被動區破壞形式分別為剪切破壞、拉張與剪切的交互破壞。

b.提出的考慮圍巖應變軟化和采空區接觸的FLAC3D數值方法，對煤層開采數值模擬實現了優化。在本構模型的選擇上，建立更貼合實際的應變軟化本構關系；對以往模擬中采空區頂底板不接觸的固有缺陷，采用“應變軟化-空-彈性”模型轉變的方法得到解決。該方法為煤層開采及需要考慮開挖后接觸的大變形工程的數值計算提供一種更貼合實際的模擬思路。

c.應變軟化本構模型能夠讓計算結果更貼合實際，但該模型需以大量巖石峰后黏聚力、內摩擦角衰減的數據為基礎進行建立，因此，還需對不同巖石峰后的力學現象進行深入研究。