雙系煤層采動導水裂隙演化規(guī)律的FDEM耦合模擬研究

李 浩，白海波，馬立強，康志勤，李志永，繆小成，武鵬飛，韋 婕

(1.太原理工大學 原位改性采礦教育部重點實驗室，山西 太原 030024；2.中國礦業(yè)大學 深部巖土力學與地下工程國家重點實驗室，江蘇 徐州 211116；3. 中國礦業(yè)大學 礦業(yè)工程學院，江蘇 徐州 211116；4. 徐州中礦巖土科技股份有限公司，江蘇 徐州 211116； 5. 晉能控股煤業(yè)集團，山西 大同 037003)

采動覆巖導水裂隙分布及演化規(guī)律的確定是開展水體下采煤、保水采煤、煤礦地下水庫設計、斷層防水煤巖柱留設、覆巖隔離注漿充填采煤設計等的重要基礎[1-3]。長期以來，覆巖導水裂隙高度往往根據“經驗公式”[4-5]來計算，一定程度上滿足了煤礦安全開采的需求。然而，經驗公式忽視了采動巖層連續(xù)→離散介質轉化力學特性，導致在某些特定開采條件下的預計結果與實際偏差較大，甚至發(fā)生突水事故[6-7]。

建立適用于采動巖體連續(xù)→離散介質轉化的耦合本構模型及相應的計算程序，是數值實現采動覆巖導水裂隙分布演化規(guī)律的重要手段。目前，學者廣泛采用“連續(xù)介質假設”，通過有限元法(FEM)、有限差分法(FDM)，分別基于理想彈塑性模型[8]、應變軟化塑性模型[9]、考慮圍壓對后繼屈服函數影響的塑性模型[10]、彈塑性損傷模型[11]，以Mohr-Coulomb屈服準則[8-9]、應變能[12]、損傷變量[11,13]為判據，得到塑性區(qū)、損傷區(qū)高度為采高的6.6～27.5倍。此外，眾多學者采用基于離散元方法(DEM)，將煤巖基質簡化為矩形塊、圓球形狀并賦予彈性體或剛體力學屬性，而基質間節(jié)理往往采用庫倫滑移模型或法向/切向接觸滑移模型[14-15]，以拉、剪應力為判據，實現采動覆巖破裂目的，根據牛頓第二定律計算得到覆巖裂隙帶高度為采高的11.2～21.3倍。

但是，上述成果大多基于單一(系)煤層地質條件得到。更重要的是，在復雜采動應力狀態(tài)下，完整覆巖往往由連續(xù)體轉化為隨機形狀的離散體，在此過程中準脆性巖石出現拉-剪混合型韌性斷裂、天然與新出現的結構面出現分離/壓縮/剪切摩擦耦合響應。連續(xù)-離散轉化過程深刻影響覆巖導水裂隙帶分布與演化，僅基于“彈塑性、損傷、線彈性斷裂力學、摩擦力學”的FEM或DEM等方法難以解決。

筆者緊密圍繞采動覆巖連續(xù)-離散轉化特性，依次建立巖石在復雜應力狀態(tài)下的混合型韌性斷裂本構、粗糙結構面壓縮-剪切摩擦本構，編制相應的有限元-離散元(FDEM)計算程序。在驗證本構模型合理性的基礎上，模擬研究雙系煤層采動導水裂隙演化規(guī)律。

1 巖石混合型韌性斷裂-壓剪摩擦本構

受荷巖石的連續(xù)→離散轉化過程包含3個顯著的力學響應特征：① 拉、剪斷裂為巖石的2種基本破壞類型，在復雜應力狀態(tài)下往往形成混合型斷裂模式；② 巖石是準脆性材料，在斷裂過程中表現出一定的韌性；③ 完全斷裂后，巖塊間粗糙結構面可出現分離、壓縮、剪切摩擦等復雜響應。以下建立巖石連續(xù)→離散介質轉化的混合型韌性斷裂-剪切摩擦(Mixed ductile fracture-shear friction，或MD-SF)本構模型，在此基礎上編制有限元-離散元FDEM數值計算程序。

1.1 混合型韌性斷裂本構

巖石韌性斷裂力學響應通過有限元FEM方法計算。不同斷裂模式下的斷裂過程分為彈性變形和韌性斷裂2個階段。彈性段本構關系可表示為

(1)

一旦達到如下應力條件，巖石開始進入韌性斷裂過程[16]：

(2)

為推導混合模式下材料完全斷裂的判據，首先根據荷載-位移曲線(圖1(c),(e))積分得到拉、剪及混合模式下的斷裂能Gc,n,Gc,S,Gc,m，表達式為

(3)

式中，Gc,j為單位面積上黏聚力單元的斷裂能，N/m;j為變量，表示n,S,m(即張拉型、剪切型、混合型)斷裂模式。

其次，在復雜應力狀態(tài)下，材料總是發(fā)生拉、剪混合型斷裂模式，為確定材料完全斷裂時拉、剪荷載的貢獻，根據剪切型斷裂能占總能量的比例，導出混合比ξ(無量綱)[16]表達式為

ξ=Gc,S/(Gc,n+Gc,S)

(4)

假設巖石在2個正交方向上的剪切性能相同，則材料在任意混合比ξ下完全斷裂的判據[17]為

(5)

圖2 拉、剪荷載下的荷載-位移曲線Fig.2 Load-displacement curves under tension and shear load

σc,n=(1-dc,n)Dc,nSc,n

(6)

(7)

(8)

(9)

由式(1)～(9)反映的巖石混合型韌性斷裂力學性質可由圖3表示,依次包括：① 彈性階段應力-應變關系；② 混合型韌性斷裂起始時刻應力表達式；③ 混合型斷裂過程中應力-應變關系；④ 混合模式下的完全斷裂準則。

圖3 混合型韌性斷裂力學響應示意Fig.3 Schematic diagram of mechanical response of mixed ductile fracture

一旦斷裂能滿足式(5)時，巖石完全斷裂。在數值計算中刪除該斷裂面上的黏聚力單元，從而消除其力學影響。

1.2 粗糙結構面壓縮-剪切摩擦本構

巖石完全斷裂后粗糙結構面的力學響應采用離散元方法(DEM)計算。在數值計算中，根據塊體接觸與相對運動趨勢，將塊體間關系劃分為分離、壓縮、剪切摩擦3類：① 當相鄰實體單元(模擬完整巖塊)的任意節(jié)點間距l(xiāng)>0，發(fā)生分離；② 若l=0且2者間存在壓應力時，塊體相互擠壓；③ 若l=0且2者間存在沿結構面的剪應力，則出現剪切摩擦。其中，法向、剪切應力判別方法根據筆者前期研究進行[16]。

若相鄰塊體分離(l>0)，2者間不存在相互作用力，其運動符合牛頓第二定律[18]。

若相鄰塊體擠壓，各要素本構關系為

(10)

式中，σn為壓縮應力，MPa；Nmax為結構面最大閉合量，mm，常參數，根據三維形貌掃描實驗確定；n為結構面閉合量，mm；Dn為結構面法向模量，GPa，本文取值為其相鄰巖塊參數。

若l=0，以SARGIN[19]的研究成果為基礎，進一步通過剪切實驗曲線擬合得到粗糙結構面壓剪摩擦本構關系為

(11)

式中，σS(σS,p)為(峰值)剪應力，MPa；S(sp)為(峰值)剪切位移，mm；DS為剪切剛度，GPa，3者均為常參數；參數p由σS,p、Sp、殘余剪應力σS,r和殘余剪位移Sr確定，σS,p,σS,r滿足如下強度準則[19]：

(12)

σS,r=σntanφr

(13)

其中，JRC為結構面粗糙度系數，無量綱；σUCS為結構面抗壓強度，MPa；φb為結構面基本摩擦角，無量綱；φr為結構面殘余變形階段的基本摩擦角，無量綱。需要強調的是，式(11)～(13)中的σp，DS，sp，p，σr，Sr等參數均為法向應力σn和結構面粗糙系數JRC的函數，其值根據圖4實驗獲取。

圖4 巖石斷裂與剪切摩擦實驗、數值模擬結果對比Fig.4 Comparison of experimental and numerical simulation results of rock fracture and shear friction

基于以上各式，建立混合型韌性斷裂-結構面剪切摩擦(mixed ductile fracture-shear friction，MD-SF)本構方程組，編制相應的有限元-離散元(FDEM)數值程序，進而模擬巖石從連續(xù)介質→離散介質的轉化過程(圖1)。

2 參數識別與本構方程驗證

巖石MD-SF本構方程的材料參數由斷裂力學實驗、結構面直剪實驗獲取。通過對比數值模擬結果與室內試驗結果，驗證本構方程合理性。

2.1 材料參數

巖石基質的力學參數根據文獻獲取[20]。張拉/剪切型斷裂參數根據三點彎曲實驗、加卸載條件下的貫穿剪切實驗獲取(表1)，基于各巖性試件的荷載-位移曲線，結合式(6)，(7)計算得到各自在拉、剪荷載下的損傷演化曲線(圖1(d),(f))；對于粗糙結構面剪切摩擦參數，通過如下方式獲?。簩r心斷口的三維形貌進行掃描，統計各巖性的JRC平均值，而后采用巴西劈裂方法制備粗糙結構面試件，選取與巖心斷口相同JRC值的試件開展直剪實驗。根據大同礦山4煤層埋深，確定法向應力σn=1，5，10 MPa。由此得到不同σn條件下粗糙結構面剪切摩擦參數(表1)。

表1 巖石力學參數Table 1 Rock mechanical parameters

2.2 本構方程驗證

MD-SF本構方程的適用性通過偏置預制裂縫的三點彎曲、貫穿剪切、粗糙結構面直剪試驗與模擬加以驗證。

通過對比不同混合比ξ條件下的斷裂力學實驗與模擬結果，驗證MD-SF本構方程中混合型韌性斷裂理論。根據斷裂力學實驗(圖1(a),(b))的試件尺寸、邊界條件，建立圖4(a)所示的數值計算模型。在模擬過程中，采用單調加載方式，位移加載速率為0.02 mm/min，加載時長10 min。模型的斷面上布置零厚度的黏聚力單元，而在其他位置處布置實體單元。單元力學屬性參見式(1)～(9)，力學參數由表1和圖1(d),(f)確定。

對于直剪試驗，以棱長100 mm立方體巖樣為基礎，采用巴西劈裂方法制備含粗糙結構面的試件，其JRC值為15(圖4(b))?？紤]到表1力學參數是根據法向應力σn=1,5,10 MPa得到的，為驗證本構模型合理性，設置σn=3,8 MPa，固定下盤，上盤的剪切加載速率為0.02 mm/min，當殘余應力穩(wěn)定后停止加載。根據直剪試驗條件建立相應的數值計算模型，結構面設置為平面，并在該處賦予法向、剪切摩擦接觸屬性(式(10)～(13))，其他位置處布置實體單元，力學參數由表1和圖1確定。通過實驗和FDEM計算程序(圖5)，得到實驗與模擬結果如圖4所示。分析圖4可得：

圖5 FDEM數值計算流程Fig.5 Block diagram of finite element-discrete element numerical calculation

(2)將σn=1,5,10 MPa條件下實驗所得力學參數(表1)代入MD-SF本構方程，通過圖5所示計算流程得到σn=3,8 MPa下的剪應力-剪位移曲線(圖4(d))。對比發(fā)現，實驗與模擬結果吻合程度較好。具體來說，當σn從3 MPa增加至8 MPa時，實驗(模擬)所得粗糙結構面峰值剪應力分別為4.82(4.77),10.68(10.94)MPa，峰值位移分別為0.89(0.93)，1.01(1.10)mm，殘余剪應力分別為2.04(2.42)，5.63(6.02)MPa，殘余剪位移分別為8.47(5)，7.13(5)mm。無論是模擬結果還是實驗結果，均反映峰值強度、位移隨σn的增加而增加，且模擬誤差介于4%～8%。但模擬所得殘余應力、殘余位移與實驗結果偏差較大，分別達到18%和69%。這是因為圖4所示數值模擬目的在于判斷MD-SF本構方程的合理性，考慮到計算成本，設置的計算時長只有250 min，殘余位移最大為5 mm。盡管如此，MD-SF本構方程也適用于粗糙結構面的剪切摩擦過程。

(3)需要注意的是，圖4(c)所示數值模擬所得曲線開頭有非線性段。這是因為在數值模型位移加載前期設置5 min的線性加載階段，加載速率逐漸從0增加至0.002 mm/min，避免形成沖擊荷載。

3 采動覆巖導水裂隙帶演化數值模擬

3.1 工程概況與數值計算模型

大同礦區(qū)西北部侏羅系、石炭—二疊系可采煤層共計5層。侏羅系4個煤層開采深度約142～248 m，現已開采完畢。目前主采煤層為石炭—二疊系山4號煤，采深約420 m，采厚6 m。經過長時間、大規(guī)模地多煤層開采，侏羅系煤系地層已形成相互連通的導水裂隙，并在煤層采空區(qū)形成高達4×104m3采空區(qū)積水(圖6)，嚴重威脅其下部山4號煤安全開采。

為研究大同礦區(qū)多煤層開采條件下采動覆巖運動及山4號煤的導水裂隙帶發(fā)育規(guī)律，根據工程地質情況建立如圖6所示的數值計算模型?？紤]MD-SF本構方程的高度非線性，模型尺寸控制為x=870 m(長),y=450 m(高),z=1 m(厚)，在各自方向上分別約束其法向位移，且在x,z方向上設置側壓力系數分別為0.52，0.74，y方向設置自重載荷。工作面推進速度為8 m/d。覆巖中主要結構面用泰森多面體表示[21-22]，其間距根據如下方法確定：依據山4號煤頂板粉砂巖周期垮落步距確定其主要結構面的間距，再根據各巖性巖石與粉砂巖抗拉強度之比，確定粉砂巖、細砂巖、中砂巖、泥巖、高嶺土砂巖、煤在x方向上的主要結構面平均間距(圖6紅色折線)分別為16.5，26.5，12.3，12.1，9.5，12.1，5.6 m，在y方向上的間距分別為10.3，20.0，10.6，5.0，4.7，6.0 m，若巖層厚度小于該值則取值巖層厚度，結構面初始寬度和JRC值見表1，而后賦予其“壓剪摩擦(SF)”力學屬性；在此基礎上，進一步在泰森多面體內部劃分六面體實體單元，單元尺寸根據各巖性巖芯平均長度之比與計算成本綜合確定，粉砂巖、細砂巖、中砂巖、泥巖、高嶺土砂巖、煤的實體單元平均尺寸分別為4.9，4.1，3.8，5.0，3.6，5.6 m(圖6黑色折線)，若巖層厚度小于該值則取值巖層厚度。實體單元賦予彈塑性本構模型[10]，而在單元邊界嵌入零厚度黏聚力單元，并賦予其“韌性斷裂-壓剪摩擦(MD-SF)”力學屬性，由此實現采動覆巖從連續(xù)體到離散體的轉變。材料參數如表1和圖1所示。

圖6 采空區(qū)積水分布及采動覆巖運動數值計算模型Fig.6 Water distribution in goaf and the numerical calculation model of mining overburden movement

3.2 模擬結果

通過上述簡化模型，得到采動覆巖導水裂隙帶數值計算結果(圖7)。由上述數值模擬結果可得：

圖7 不同工作面推進距離下巖石連續(xù)-離散轉化、覆巖裂隙分布及Mises應力場特征Fig.7 Characteristics of rock continuous-discrete transformation,overburden fracture distribution and Mises stress field under different mining face advancing distance

(1)在工作面推進方向上，覆巖導水裂隙介于工作面煤壁前方44.7 m與煤壁后方12.4 m之間。巖體脆-韌性力學特征、超前支承壓力、頂板天然結構面位置的疊加效應是產生上述現象的主要原因。山4號煤開采過程中，超前支承壓力范圍約為49 m，這會導致該范圍內脆性砂巖的結構面承受較大壓剪應力；隨著工作面向前推進，頂板巖體垮落，其附近巖體失去側向支撐而旋轉，導致在脆性明顯的砂巖內較大寬度的超前導水裂隙形成；反之，若工作面煤壁剛好經過前一條導水裂隙，且超前支承壓力不足以使后一條結構面出現法向位移時，則出現滯后型導水裂隙。與此同時，2條天然結構面之間的完整巖體內部出現復雜的拉、剪應力，導致完整巖體主要產生混合型斷裂裂縫，并最終由連續(xù)體轉化為離散體，使得覆巖導水裂隙進一步復雜化。

(2)采動覆巖導水裂隙高度。在山4號煤覆巖內且與煤層間距45，82，125 m處布置3條導水裂隙寬度測線(圖6)，統計裂隙寬度在1mm以上的導水裂隙總寬度隨工作面推進距離演化規(guī)律如圖8所示。

圖8 覆巖導水裂隙總寬度隨工作面推進距離曲線Fig.8 Curves of total width of water conducting fissure in overburden with mining advancing distance

圖8中上、中測線位于高嶺土砂巖層中，在工作面推進過程中，2個測線處導水裂隙總寬度呈現先增加后減小的趨勢，即從0增加至0.476 m(或1.053 m)再降低至0.143 m(或0.520 m)。以圖7(b),(c)為例解釋該現象的原因：當工作面推進至123 m時出現周期來壓，頂板砂巖層完全斷裂并發(fā)生約4°的轉動，裂隙寬度迅速增加至0.173 m；而隨著工作面推進至184 m，采動應力逐漸穩(wěn)定，該巖層回轉約2°，使得裂隙寬度降低至0.134 m。以上數據表明，煤層開采過程中更容易出現頂板突水事故。

由圖7,8可知，高嶺土砂巖中的采動覆巖裂隙寬度、密度遠遠小于其下部砂巖層。這不僅是因為砂巖層受采動應力較大，更重要的是，高嶺土砂巖具有更強的韌性破壞特征，在斷裂過程中大量能量轉化為塑性功，而高嶺土砂巖的裂縫表面能(即彈性能)是砂巖的63.57%，這意味著采動裂隙在脆性較為明顯的砂巖內更容易擴展。此外，砂巖較為顯著的脆性斷裂響應，導致采動裂隙快速發(fā)育至砂巖層頂界，巖塊偏轉并與周圍未破斷巖體形成最寬達3.347 m的寬大裂隙，而隨著采動應力降低巖塊回轉，裂隙寬度又會逐漸降低，并穩(wěn)定在1.726 m左右，形成波動升降曲線(圖8藍色線)。

綜上，山4號煤的采動覆巖導水裂隙帶高度不超過125 m，遠大于經驗公式所得74.1 m的結果。

為驗證覆巖導水裂隙帶模擬結果，進一步開展井下注水試驗。鉆孔距離開切眼280 m，仰角65°，鉆孔斜長135 m(垂高120 m)，待工作面開采完畢后采用分段注水方法得到水流漏失量。由圖9可知，當鉆孔斜長從7.6 m增加至23.4 m時，水流漏失量從2.2 L/min快速增加至26.5 L/min，意味著鉆孔進入導水裂隙帶；直至鉆孔斜長達到117.6～123.9 m，水流漏失量為4.2～2.7 L/min，達到穩(wěn)定。實測結果表明，山4號煤采動覆巖導水裂隙帶垂高約為112.3 m，從而驗證上述模擬結果的合理性。

圖9 注水試驗結果Fig.9 Water injection test results

4 結 論

(1)建立混合型韌性斷裂-剪切摩擦本構模型及相應的FDEM數值算法，采用實體單元、黏聚力單元和接觸對表示巖石基質、非貫通裂隙及貫通結構面，計算準脆性巖石彈塑性變形、混合型韌性斷裂及粗糙結構面分離/壓縮/剪切摩擦力學響應，數值實現巖體從連續(xù)介質到離散介質的轉化過程。

(2)相比較實驗結果，數值模擬所得不同拉、剪斷裂能混合比下巖石斷裂峰值荷載、峰值位移、斷裂能誤差為0.5%～3.0%，2.0%～11.0%，4.0%～16.0%；模擬所得不同法向應力下粗糙結構面剪切摩擦峰值荷載誤差為4%～8%，由此驗證混合型韌性斷裂-剪切摩擦本構模型的合理性。

(3)模擬條件下，開采過程中覆巖導水裂隙總寬度是覆巖運動穩(wěn)定后的2.26～7.11倍，表明煤層開采過程中更容易出現頂板突水事故。

(4)高嶺土砂巖顯著的韌性破壞特征有效控制其內導水裂隙發(fā)育高度，但在雙系煤層重復開采擾動下，覆巖導水裂隙發(fā)育高度約125 m，裂采比為20.8，遠超經驗公式所得結果。