尹尚先,孟浩鵬,錢雙彬
(1.華北科技學院 安全工程學院,北京 101601;2.華北科技學院 理學院,北京 101601)
隨著時代發(fā)展,數(shù)值方法在巖土工程領域的應用逐漸普及,其中FLAC3D有限差分數(shù)值方法優(yōu)點眾多,在采礦領域得到廣泛認可與應用,眾多學者借助FLAC3D圍繞煤炭安全開采進行了大量研究。董書寧等[1]在改造奧灰頂部巖層段的判別準則研究中,利用數(shù)值計算分析了采深、采高和采寬等尺寸效應對底板破壞帶的影響;劉偉韜等[2]通過正交試驗設計,進行了7 個主控因素對底板破壞深度的影響研究,并對主控因素進行排序;劉新民[3]采用現(xiàn)場和模擬結合的方法進行沿空留巷對底板破壞深度的研究,認為無煤柱式的沿空留巷開采技術不會對底板破壞深度造成較大影響;朱斯陶[4]、朱廣安[5]、田雨桐[6]等的研究表明,數(shù)值模擬是研究采動影響下斷層活化規(guī)律的重要手段之一,可實現(xiàn)采動影響下對斷層多方位的定量分析;朱慶偉[7]、甘智慧[8]等運用數(shù)值方法對采動影響下覆巖結構演化和地面沉降進行了研究。上述研究成果在數(shù)值計算中一般采用摩爾-庫倫本構模型,該本構模型為理想彈塑性模型,對圍巖塑性屈服后的狀態(tài)無法準確描述,并且對拉格朗日法遵循連續(xù)介質假設而導致網格發(fā)生大變形但節(jié)點不接觸的固有缺陷[9]未進行深入研究,這2 個問題可能拉大數(shù)值模擬同實際情況之間的差距。大變形條件下不接觸的固有缺陷,導致無法模擬采空區(qū)頂板垮落后頂?shù)装褰佑|的應力傳遞現(xiàn)象,即“采空區(qū)不接觸”現(xiàn)象。觀察模擬工作面回采結果發(fā)現(xiàn):頂?shù)装褰佑|后并不進行接觸計算;頂?shù)装鍛咎幱谛箟籂顟B(tài)。真實情況是采空區(qū)在周期來壓過程中,基本頂發(fā)生周期性垮落與底板接觸,采空區(qū)的應力恢復隨時間變化呈指數(shù)函數(shù)關系[10]。因此,利用FLAC3D研究煤層回采后底板破壞時,應對圍巖應變軟化和采空區(qū)接觸進行考慮。
筆者將針對FLAC3D模擬工作面回采中本構模型的選擇和采空區(qū)不接觸的固有缺陷進行研究,以河北開平煤田林西礦2023 工作面底板導水裂隙帶實測為工程背景,建立考慮應變軟化和采空區(qū)接觸的工作面回采模型,并結合力學分析對模擬結果進行解釋,達到提高工作面回采數(shù)值計算準確性的目的。
林西礦位于河北開平煤田東南翼,地質構造以褶皺為主,開平主向斜穿過井田東側深部,西部有杜軍莊背斜和黑鴨子向斜;井田內大型斷層較少,小型斷層較發(fā)育;井田內揭露地層由老到新為:奧陶系、石炭系、二疊系和第四系。石炭系、二疊系為含煤地層,煤系基底為奧陶系灰?guī)r(簡稱奧灰)。林西礦目前主要生產水平為11 水平(?850 m)和12 水平(?1 000 m)。
12 煤位于二疊系下統(tǒng)趙各莊組,煤層底板至奧陶系自上而下為:石炭系上統(tǒng)開平組,上以K6 灰?guī)r頂界面與二疊系下統(tǒng)趙各莊組分界,下以唐山組K3 灰?guī)r頂界面與開平組分界;石炭系上統(tǒng)唐山組,上以K3 灰?guī)r頂界面與開平組分界,下以G 層鋁鐵質泥巖的底界面與奧陶系石灰?guī)r呈平行不整合接觸。趙各莊組為主要含煤組,厚度33.55~61.20 m,平均48.35 m,含11、12 兩層可采煤層;本組巖性頂部主要為黑色泥巖和灰色砂巖及褐灰色砂巖,其中砂巖向深部逐漸發(fā)展為黑色泥巖;中部及底部主要為灰色細?粗砂巖,淺部間有3~4 層砂礫巖,向深部礫巖直徑逐漸變小,砂巖粒度也逐漸變細且大部分砂巖為泥巖所替代。開平組層厚55.09~91.04 m,平均76.81 m,主要巖性為黑色粉砂質泥巖,砂巖次之,其中砂巖比例由淺向深逐漸減少。唐山組層厚64.44~80.55 m,平均70.18 m;地層巖性除K1、K2、K3 灰?guī)r和G 層鋁鐵質泥巖外,主要為黑色-深灰色泥巖和灰色砂巖。
根據(jù)林西礦深部ZK7 奧灰水位觀測孔鉆孔水位及12 煤層底板至奧灰頂界面間距計算得到:當12 煤層底板高程小于?936.3 m 時,工作面回采期間突水系數(shù)超過0.06 MPa/m。按照《煤礦防治水細則》,如果突水系數(shù)超限問題不解決,深部區(qū)域將無法進行安全帶壓回采,嚴重影響礦井采掘接替,威脅礦井生存。因此,有必要對12 煤層底板采動導水破壞帶深度進行實測及底板破壞分析等工作。
林西礦深部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。
為觀測林西礦2023 工作面底板采動導水破壞帶發(fā)育情況,在2023 工作面西側2023-2 工作面回風巷設計D01-1、D01-2 鉆孔,鉆孔鉆至2023 工作面法向向下30 m 左右位置,如圖1 所示,隨后進行分段壓水試驗,記錄不同位置流量穩(wěn)定后的漏失量,繪制鉆孔的漏失量曲線(圖2)。觀察鉆孔的漏失量曲線發(fā)現(xiàn),在垂距24 m 以后,漏失量下降明顯,注水壓力在漏失量穩(wěn)定后回彈明顯且接近初始注水壓力,說明該段裂隙不發(fā)育且貫通性差。在壓水試驗基礎上,為直觀了解底板破壞情況以及對壓水試驗結果進行驗證,對D01-2鉆孔進行鉆孔成像,獲取成像數(shù)據(jù)如圖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
根據(jù)D01-2 鉆孔成像結果顯示,圖3a 中14.2 m處巖性為灰色細砂巖,處于底板破壞帶邊緣,該處巖層完整性較好,但存細微層狀裂隙;圖3b、圖3c 中15.9~19.3 m 處,巖性為細砂巖,這2 處巖層及其之間的巖石破碎嚴重,縱橫裂隙發(fā)育明顯;圖3d 中26.6 m處,巖性為細砂巖,該處巖層及以下巖層完整,未見明顯裂隙。認為D01-2 鉆孔漏失量曲線在30 m 處的回彈是由于壓水試驗中的封堵裝置壓力過大致使原生裂隙張開導致。
綜合D01-1、D01-2 孔壓水試驗實測結果和鉆孔成像數(shù)據(jù),并參考鄰近趙各莊礦近似埋深的1237、2137 工作面正常底板實測采動導水破壞帶深度23、25 m[11],最終確定林西礦2023 工作面正常底板采動導水破壞帶深度為24 m。
運用朗肯土壓力理論對煤層底板破壞進行定性分析[12-13],圖4 表示具有半無限平面的煤層開采走向剖面圖,將底板塑性破壞劃分為Ⅰ區(qū)(主動區(qū))、Ⅱ區(qū)(過渡區(qū))、Ⅲ區(qū)(被動區(qū))3 個區(qū)域。圖4 表示回采過程中煤層底板受力狀態(tài)的應力圓與底板強度包線之間的關系。
圖4 煤層底板破壞分區(qū)Fig.4 Coal seam floor failure zones
在原巖應力狀態(tài)下,煤層下深度為z處單元體的應力為豎向應力σv(σz)、最大水平應力σH、最小水平應力 σh,取σH=σ1(最大主應力)、σv=σ3(最小主應力)的情況,用圖5 中的應力圓①表示,此狀態(tài)下應力圓距強度包絡線較遠,底板處于彈性平衡。隨工作面的推進,底板地應力受工程擾動,在超前應力作用下 σv將逐漸增大,假設σH保持不變,此時應力圓的半徑先減小后增大,如圖5 中方向向右的箭頭所示,若剪應力達到底板抗剪強度,應力圓與強度包絡線相切,該處底板達到被動極限平衡狀態(tài),如圖5 中應力圓②所示,此時σH=σ3、σv=σv1=σ1(σv1為應力圓②狀態(tài)下最大主應力),根據(jù)摩爾?庫倫理論可知,當工作面前方底板達到或超過極限平衡狀態(tài)并隨應力的持續(xù)作用,Ⅰ區(qū)發(fā)生塑性變形,伴隨體積膨脹以壓力形式通過Ⅱ區(qū)向采空區(qū)方向(Ⅲ區(qū))傳遞,導致底鼓同時形成連續(xù)的滑移面。當工作面推過之前超前應力作用的底板位置后,最大、最小主應力發(fā)生轉變,由于工作面的推進,σv減小,假設σH保持不變,此時應力圓的半徑先減小后增大,如圖5 中方向向左的箭頭所示,直至應力圓與強度包絡線相切,該處底板達到主動極限平衡狀態(tài),如圖5中應力圓③所示,此時σH=σ1、σv=σv2=σ3(σv2為應力圓③狀態(tài)下最小主應力),根據(jù)摩爾-庫倫理論可知,當工作面前方底板達到或超過極限平衡狀態(tài)并隨應力的持續(xù)作用,Ⅲ區(qū)發(fā)生塑性變形,導致采空區(qū)底板隆起,同時形成連續(xù)的滑移面。
圖5 煤層底板極限平衡狀態(tài)Fig.5 Limit equilibrium state of coal seam floor
為比較底板的被動極限平衡狀態(tài)和主動極限平衡狀態(tài),根據(jù)極限應力圓與強度包絡線之間所得的關系式[14]:
式中:φ為巖石內摩擦角;c為巖石純剪切強度(黏聚力)。
將應力圓②和③的表達式代入式(1)可得底板被動和主動極限平衡狀態(tài)之間的關系:
需注意:①實際煤層開采工作面斜長有限,同時工作面斜長與底板采動破壞帶深度密切相關,導致對底板破壞進行定量分析存在困難,但采用半無限平面對底板破壞進行定性力學分析是可行的;② 對煤層下深度為z處單元體最大水平應力 σH保持不變的假設同實際不符,Ⅰ區(qū)的 σH應當隨頂板周期來壓而發(fā)生周期性變化,Ⅰ區(qū)在豎向應力作用下,應力必然向四周傳遞,由于采空區(qū)這一臨空面的存在,Ⅱ區(qū)、Ⅲ區(qū)的 σH隨Ⅰ區(qū)的 σH的變化而變化;③從回采過程中底板受力狀態(tài)的應力圓與底板強度包絡線之間的關系可知,達到應力圓②的狀態(tài)比應力圓③要困難,即Ⅰ區(qū)達到極限平衡狀態(tài)較難。
底板破壞模擬一般采用摩爾-庫倫本構模型,此模型為理想狀態(tài)的彈塑性模型,不考慮黏聚力、內摩擦角、抗拉強度等材料參數(shù)隨塑性變形的變化情況[15],同實際地質材料受力特征不符。此模型在底板塑性屈服前能較好地反映底板變形,但在底板發(fā)生塑性屈服階段開始后,同實際破壞有較大差異。底板塑性屈服后,在應力作用下呈應變軟化行為,底板在地應力的作用下產生微裂紋及巖體的相對滑動,底板強度將不斷降低并且越來越缺乏彈性,直到破壞以及剪切帶的產生。利用FLAC3D建立應變軟化摩爾-庫倫地層模型,能夠對圍巖塑性破壞后的力學狀態(tài)更準確表述。
實際巖石峰后的應變軟化過程的彈塑性剛度矩陣為一個不定矩陣[16],導致應變軟化問題求解困難。為避免這一情況,可將巖石峰后應變軟化過程簡化為一系列的脆塑性過程[17-19]。FLAC3D內置的應變軟化模型為基于經典彈塑性理論將實際的應變軟化過程中黏聚力、內摩擦角、剪脹角與塑性剪切應變的函數(shù)近似為一組首尾相連的分段線性函數(shù)的模型。
基于巖石軟化相關文獻[20-21],將其實驗數(shù)據(jù)擬合為巖石峰后黏聚力、峰后內摩擦角隨塑性剪切應變εs的指數(shù)函數(shù):
為驗證擬合函數(shù)的合理性,將函數(shù)嵌入應變軟化本構模型,對參考文獻中的泥巖進行單軸壓縮數(shù)值模擬。模擬采用單軸壓縮試驗常用直徑(D)∶高(H)為1∶2 的圓柱體進行試驗,為較好地模擬真實試驗和呈現(xiàn)應力-應變全過程曲線,對試件端部施加恒定速度代替試件受壓情況。根據(jù)上述條件進行2 種本構模型的試驗,得到σ-ε曲線(圖6)。屈服前,2 種材料的σ-ε曲線一致且基本符合線彈性;屈服后,摩爾-庫倫材料與應變軟化材料的σ-ε曲線明顯不同。應變軟化材料的σ-ε曲線同真實試驗有較高的吻合度,可較好地描述巖石屈服軟化后的力學特征,說明根據(jù)式(3)、式(4)建立的應變軟化模型是合理可行的。
圖6 2 種本構模型應力-應變曲線Fig.6 Stress-strain curves of two constitutive models
3.2.1 初始模型建立
數(shù)值模擬的工況條件以林西礦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 奧灰水位觀測孔鉆孔信息,對地層傾角、巖石力學參數(shù)相近及薄巖層進行適當簡化,最終模型尺寸長×寬×高為1 000 m×300 m×240 m,剖分網格數(shù)量804 000 個,節(jié)點數(shù)量833 748 個(圖7),數(shù)值模型的巖石力學參數(shù)見表1。
圖7 煤層底板破壞數(shù)值模型Fig.7 Numerical model of coal seam floor failure
表1 數(shù)值模型巖石力學參數(shù)Table 1 Rock mechanical parameters of the numerical model
本構模型選用本文提出的應變軟化模型,對模型施加10 m/s2豎直向下的重力加速度,在模型頂部施加22.618 MPa 應力代替未建模的上覆巖層,底部采用豎直位移約束,四周采用水平位移約束,并在四周施加隨深度增加的側向應力。側向應力的大小參考沉積巖的應力分布規(guī)律[22],根據(jù)模型煤層埋深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 考慮采空區(qū)接觸的方法
實際煤層開采中,基本頂隨周期來壓垮落,垮落后的頂板與底板接觸,發(fā)生頂?shù)装逯g的應力傳遞。但在FLAC3D模擬工作面回采中發(fā)現(xiàn),頂?shù)装逶趹ψ饔孟陆佑|后,并不會進行應力接觸計算,在大變形模式下,甚至可清楚觀察到頂?shù)装灏l(fā)生交叉的現(xiàn)象,這與實際情況不符,導致模擬與實際產生巨大偏差,所以在模擬工作面回采中需要對采空區(qū)進行接觸模擬。為體現(xiàn)真實回采中應力變化,開挖步距參考真實工作面周期來壓步距,模擬過程中工作面以20 m 為步距循環(huán)開挖。煤層回采后,回采區(qū)域由“應變軟化”模型轉為“空”模型。將實際采空區(qū)頂板垮落后堆積的碎石假設為彈性整體,利用“彈性(各向同性)”模型替換“空”模型,達到模擬采空區(qū)頂板垮落后頂?shù)装鍛佑|的目的。在彈性體參數(shù)的確定上,由于采空區(qū)的碎石是頂板垮落產生,將其視為裂隙發(fā)育的彈性整體,其彈性模量將大幅衰減,泊松比有所上升[23]。
數(shù)值計算時,工作面回采3 個循環(huán)后,模擬頂板垮落后的頂?shù)装褰佑|,即“空”模型以20 m 為步距循環(huán)向“彈性”模型轉化,為避免彈性體對側向產生應力傳遞,在切眼、側幫及終采線附近不改變“空”模型。在工作面中心頂?shù)装逦恢梅謩e布置測點,記錄回采全過程豎直方向應力值和豎直位移量(圖8、圖9)。根據(jù)測點記錄發(fā)現(xiàn),模擬采空區(qū)頂板垮落后的接觸與對采空區(qū)不做處理的結果具有顯著差異,考慮采空區(qū)接觸的結果更貼近實際??紤]采空區(qū)接觸的采空區(qū)頂?shù)装鍛Φ玫揭欢ǔ潭鹊幕謴?,而不考慮接觸的采空區(qū)頂?shù)装鍛咎幱谛箟籂顟B(tài);考慮采空區(qū)接觸的采空區(qū)頂?shù)装遑Q直位移量明顯小于不考慮接觸的情況,并且考慮采空區(qū)接觸后的底板位移量表現(xiàn)出小幅的回落。
圖8 煤層頂?shù)装遑Q直方向應力Fig.8 Stress in vertical direction of coal seam roof and floor
圖9 煤層頂?shù)装遑Q直位移Fig.9 Vertical displacement of coal seam roof and floor
采用“應變軟化-空-彈性”模型轉變的方法,達到模擬采空區(qū)頂板垮落后應力傳遞的效果,彌補了以往煤層開采模擬中采空區(qū)頂?shù)装宀唤佑|的固有缺陷。
運用自定義應變軟化本構關系和考慮采空區(qū)接觸的數(shù)值方法進行目標工作面的回采模擬,對地應力σZ∶σX∶σY=1∶0.758∶1.133 模擬過程中(回采40、100、800 m 平衡后)的塑性區(qū)(圖10)和累計塑性剪切應變率大于0.01 的區(qū)域(圖11)進行切片展示。
圖10 不同回采距離塑性區(qū)分布Fig.10 Plastic zone distribution in different mining distance
圖11 不同回采距離塑性剪切應變突出區(qū)Fig.11 Strain-shear-plastic outburst zone in different mining distance
1) 底板塑性區(qū)分析
模擬回采0~100 m 過程中,塑性區(qū)深度迅速增大,后隨回采的進行,塑性區(qū)深度基本穩(wěn)定在23 m 左右,同實際底板導水裂隙帶深度的實測結果一致;底板塑性區(qū)上部的狀態(tài)為“shear-p、tension-p”,為剪切屈服和張拉屈服共存狀態(tài),分布形態(tài)隨回采呈周期性分布;底板塑性區(qū)下部的狀態(tài)基本為“shear-p”,為剪切屈服狀態(tài),將上部“shear-p、tension-p”狀態(tài)包圍。通過結合該區(qū)域應力及位移分布,底板塑性區(qū)上部的最小主應力呈拉應力,下部呈壓應力;在底板2 種塑性狀態(tài)交界附近的位移量有一定突變。結果表明,底板破壞區(qū)域上部為剪切、拉張交互破壞,下部為剪切破壞。
2) 塑性剪切區(qū)分析
底板塑性區(qū)全區(qū)包含塑性剪切狀態(tài)(圖10),選取累計塑性剪切應變率大于0.01 的區(qū)域進行顯示(圖11),其分布形態(tài)為斜向采空區(qū)的半包圍面狀結構。
截取2 種極限地應力條件下的3 個開采循環(huán)距離進行底板破壞分析,對塑性區(qū)、塑性剪切應變率大于0.01 進行整合處理(圖12)。發(fā)現(xiàn)原巖應力的改變對底板破壞規(guī)律幾乎沒有影響,塑性剪切應變突出區(qū)域將剪切、拉張交互破壞區(qū)(“shear-p、tension-p”)同剪切破壞區(qū)(“shear-p”)劃分開,即剪切帶內側為剪切、拉伸交互破壞,外側為剪切破壞。一般認為,巖石沿最大有效剪應力面形成剪切破裂面,而塑性剪切應變集中區(qū)域與最大有效剪應力集中區(qū)域基本一致,故認為滑移面即剪切破壞面是沿塑性剪切應變集中區(qū)域分布。通過利用優(yōu)化后的數(shù)值方法得到的結果可知,煤層回采后底板破壞類型可分為剪切和拉張交互破壞、剪切破壞2 種類型,并可根據(jù)塑性剪切應變集中度,對底板滑移線進行三維可視化顯示,底板滑移面呈斜向采空區(qū)的半包圍面狀陣列分布。
圖12 2 種地應力下局部模擬結果Fig.12 Local simulation results under two ground stresses
將數(shù)值結果與本文提到的底板破壞力學分析結合,工作面回采前,底板巖體處于地應力平衡狀態(tài),受采動影響底板初始應力狀態(tài)被打破,在采空區(qū)前方底板(Ⅰ區(qū))產生超前支撐壓力,由于采空區(qū)這一臨空面的存在,Ⅰ區(qū)的大部分應力和位移通過Ⅱ區(qū)向Ⅲ區(qū)傳遞,致使采空區(qū)底板整體處于高圍壓、低軸壓狀態(tài),底板破壞形式整體表現(xiàn)為塑性剪切破壞。對塑性剪切應變率較高區(qū)域進行顯示,該區(qū)域將底板塑性破壞區(qū)分為上下部分。朗肯土壓力理論中3 個破壞區(qū)是根據(jù)塑性剪切程度進行劃分,認為模擬結果中底板塑性區(qū)可根據(jù)塑性剪切應變率的大小進行劃分,劃分結果可與朗肯土壓力理論中3 個區(qū)進行對應。本文底板塑性區(qū)中塑性剪切應變率突出的區(qū)域(滑移面)為Ⅱ區(qū)和Ⅲ區(qū)之間的分界,Ⅰ區(qū)未顯現(xiàn),同前文“Ⅰ區(qū)達到極限平衡狀態(tài)較難”對應。Ⅱ區(qū)的破壞形式為剪切破壞,Ⅲ區(qū)的破壞形式為拉張和剪切的交互破壞。
a.運用朗肯土壓力理論并結合考慮圍巖應變軟化和采空區(qū)接觸的數(shù)值方法研究了河北林西礦深部底板破壞特征,根據(jù)塑性剪切應變率的變化,對底板滑移面實現(xiàn)了三維顯示,并將底板塑性區(qū)與朗肯土壓力中的主動區(qū)、過渡區(qū)和被動區(qū)相對應,其中過渡區(qū)、被動區(qū)破壞形式分別為剪切破壞、拉張與剪切的交互破壞。
b.提出的考慮圍巖應變軟化和采空區(qū)接觸的FLAC3D數(shù)值方法,對煤層開采數(shù)值模擬實現(xiàn)了優(yōu)化。在本構模型的選擇上,建立更貼合實際的應變軟化本構關系;對以往模擬中采空區(qū)頂?shù)装宀唤佑|的固有缺陷,采用“應變軟化-空-彈性”模型轉變的方法得到解決。該方法為煤層開采及需要考慮開挖后接觸的大變形工程的數(shù)值計算提供一種更貼合實際的模擬思路。
c.應變軟化本構模型能夠讓計算結果更貼合實際,但該模型需以大量巖石峰后黏聚力、內摩擦角衰減的數(shù)據(jù)為基礎進行建立,因此,還需對不同巖石峰后的力學現(xiàn)象進行深入研究。