覆巖破壞特征數(shù)值模擬研究現(xiàn)狀及進(jìn)展

王曉蕾, 王 晶

(1.呂梁學(xué)院礦業(yè)工程系， 呂梁 033000； 2. 煤礦機(jī)械裝備維護(hù)與檢測(cè)試驗(yàn)呂梁市重點(diǎn)實(shí)驗(yàn)室， 呂梁 033000)

中國(guó)是一個(gè)煤礦大國(guó)，煤炭在中國(guó)能源結(jié)構(gòu)中占有重要地位[1]，對(duì)于中國(guó)經(jīng)濟(jì)發(fā)展具有重要的促進(jìn)作用，是中國(guó)經(jīng)濟(jì)不可或缺的能源，雖然近幾年中國(guó)大力發(fā)展新能源，但是效果不明顯[2-4]?！赌茉窗l(fā)展戰(zhàn)略行動(dòng)計(jì)劃(2014—2020年)》指出，到2020年，一次能源消費(fèi)總量控制在48億t標(biāo)準(zhǔn)煤左右，煤炭消費(fèi)總量控制在42億t左右，其中煤炭占一次能源比重為62%。由此可見(jiàn)未來(lái)相當(dāng)長(zhǎng)的時(shí)間內(nèi)中國(guó)能源還是以煤炭為主[5-8]。

隨著中國(guó)煤炭資源的大量回采，優(yōu)勢(shì)地質(zhì)條件下煤炭資源已趨于殆盡，煤礦企業(yè)不得不開(kāi)采劣勢(shì)地質(zhì)條件下煤礦資源[9]，同時(shí)，隨著開(kāi)采深度的不斷增加，煤層瓦斯含量較大，煤層透氣性較低，嚴(yán)重威脅煤礦安全生產(chǎn)[10]。

煤層開(kāi)采前其應(yīng)力處于原始平衡狀態(tài)，當(dāng)煤炭資源回采后，受采動(dòng)和支撐應(yīng)力的作用[11]，煤層上覆巖層必然發(fā)生移動(dòng)，最終破斷，自上而下形成有規(guī)律的三帶[12]，煤層破壞高度對(duì)于煤礦水害治理和瓦斯防治具有重要意義[13]。煤層發(fā)育特征受地質(zhì)特征、開(kāi)采條件等多種因素影響，煤層破壞高度的有效確定對(duì)于煤礦安全生產(chǎn)具有重要意義[14]。對(duì)于煤層破壞高度的研究目前主要采用理論分析、模擬實(shí)驗(yàn)、現(xiàn)場(chǎng)測(cè)試的方法，其中數(shù)值模擬實(shí)驗(yàn)操作簡(jiǎn)便，模擬不同工況的覆巖破壞，是一種非常好的測(cè)試技術(shù)[15]。

現(xiàn)基于煤層開(kāi)采覆巖破壞特征數(shù)值模擬研究，總結(jié)煤層開(kāi)采覆巖破壞特征數(shù)值模擬研究現(xiàn)狀，分析煤層開(kāi)采覆巖破壞特征數(shù)值模擬研究的不足，針對(duì)目前煤層開(kāi)采覆巖破壞特征數(shù)值模擬研究存在的問(wèn)題，提出未來(lái)煤層開(kāi)采覆巖破壞特征數(shù)值模擬研究發(fā)展趨勢(shì)。

1 煤層開(kāi)采覆巖破壞

礦井中巷道開(kāi)掘之前不受任何人類(lèi)影響的巖體稱(chēng)為原巖體[16]。存在于原巖中的力為原巖應(yīng)力，巷道采掘前巖體處于彈性變形狀態(tài)[17]。假設(shè)上覆巖層為均質(zhì)連續(xù)介質(zhì)，巖體為半無(wú)限體，當(dāng)距地表一定高度時(shí)，其主要受上覆巖層的壓力作用[18-19]。

當(dāng)巷道開(kāi)挖后，原始應(yīng)力被破壞，應(yīng)力重新分布，巷道周?chē)l(fā)生應(yīng)力集中現(xiàn)象[20]。臨空煤巖體應(yīng)力升高，同時(shí)，其受力狀態(tài)由三向變?yōu)閮上蚴芰顟B(tài)，此時(shí)承重主要是采空區(qū)上部的煤巖體，其承受上部巖層應(yīng)力作用[21]。并將這種作用傳遞到工作面兩端的煤壁上以及工作面前后的煤巖體上[22]。在采空區(qū)區(qū)域內(nèi)為煤巖體卸壓區(qū)。而在工作面兩端煤壁上形成了應(yīng)力集中區(qū)[23]，如圖1所示。

圖1 工作面支撐壓力分布Fig.1 Working face support pressure distribution

采煤工作面前、兩端煤壁、采空區(qū)冒落的矸石分別承受的為前支承壓力[24]、側(cè)向的支承壓力以及承載支撐作用的后支撐壓力[25]。由于采空區(qū)的卸壓以及應(yīng)力集中區(qū)的存在，必然使得煤層頂板發(fā)生破壞[26]。最終自上而下形成有規(guī)律的三帶，如圖2所示。

Ⅰ為垮落帶； Ⅱ?yàn)榱严稁В?Ⅲ為彎曲下沉帶圖2 覆巖三帶Fig.2 Three zones of overburden

2 覆巖破壞數(shù)值模擬方法

隨著科學(xué)技術(shù)的不斷發(fā)展，中外學(xué)者根據(jù)覆巖破壞的特征及機(jī)理提出了眾多數(shù)值模擬軟件，主要包括UDEC、FLAC、RFPA、CDEM、PFC、ABQUES。這些模擬軟件對(duì)于煤礦安全生產(chǎn)起到了非常重要的作用。

2.1 UDEC

UDEC是一款基于離散單元法理論的一款計(jì)算分析程序[27]。它是處理不連續(xù)介質(zhì)的二維離散元程序，可以模擬非連續(xù)介質(zhì)承受靜載或者動(dòng)載作用下的響應(yīng)[28]。

該方法求解主要包括物理方程和運(yùn)動(dòng)方程兩部分，計(jì)算時(shí)采用迭代進(jìn)行，直到每塊物體不再出現(xiàn)不平衡力和力矩為止[29-30]。其物理方程為

F=-KU

(1)

式(1)中：F代表法向應(yīng)力；K代表節(jié)理的法向剛度系數(shù)；U代表塊體間的疊合量。

由于運(yùn)動(dòng)方程符合牛頓第二運(yùn)動(dòng)定律，采用動(dòng)態(tài)松弛法進(jìn)行求解，其方程為

(2)

式(2)中：u代表塊體形心的位移；c代表黏性阻尼系數(shù)；k代表剛度系數(shù)；m代表塊體單元質(zhì)量；f代表外載荷；t代表時(shí)間。

采用動(dòng)態(tài)松弛法[式(2)]進(jìn)行求解，是一種顯式求解。計(jì)算過(guò)程中對(duì)于時(shí)間差要求小，而且要求合理確定阻尼系數(shù)。

王國(guó)鋒[31]采用UDEC模擬軟件對(duì)趙莊煤礦覆巖三帶發(fā)育特征進(jìn)行了研究，以1307工作面作為模擬面，其模型長(zhǎng)度和高度分別是400 m和60 m，其模型如圖3所示。

圖3 UDEC數(shù)值模擬模型Fig.3 Simulation model of UDEC

以工作面煤以及覆巖巖層物理力學(xué)性質(zhì)作為模型參數(shù)，其力學(xué)參數(shù)如表1所示。

表1 UDEC模擬力學(xué)參數(shù)Table 1 Mechanical parameters of UDEC

模擬過(guò)程中選取左側(cè)邊界100 m進(jìn)行開(kāi)挖，得出開(kāi)挖過(guò)程中覆巖變化如圖4所示。

圖4 UDEC模擬覆巖破壞特征Fig.4 Fault characteristics of overburden of FLAC

由圖4可知，當(dāng)工作面回采20 m后，覆巖產(chǎn)生了新邊界，但未發(fā)生垮落，隨著工作面的繼續(xù)推進(jìn)，當(dāng)回采到30 m時(shí)，煤體上方發(fā)生了垮落，垮落高度可達(dá)13 m，當(dāng)回采到80 m時(shí)，采空區(qū)受到上覆巖體的作用，采空區(qū)被壓實(shí)，冒落帶達(dá)到最大值36 m，當(dāng)工作面回采到150 m時(shí)，裂隙帶發(fā)育距頂板48.5 m。該位置布置瓦斯抽放巷對(duì)于煤礦安全生產(chǎn)具有重要意義。

徐剛等[32]采用UDEC對(duì)崔家溝煤礦覆巖裂隙演化進(jìn)行了研究，分析了覆巖三帶裂隙演化規(guī)律，得出了垮落帶和斷裂西發(fā)育高度；黃清鋒等[33]采用該模擬軟件對(duì)某礦6303大傾角工作面進(jìn)行了模擬，得出了工作面覆巖破壞呈V形分布，采場(chǎng)后方破壞范圍最大，中部次之；許剛剛[34]采用該方法對(duì)礦井才懂裂隙演化規(guī)律進(jìn)行了研究，分析采動(dòng)裂隙隨開(kāi)采寬度變化特征，對(duì)于采動(dòng)裂隙定量化提供了實(shí)踐基礎(chǔ)；湯玉兵等[35]采用數(shù)值模擬對(duì)薛沙河下煤層覆巖裂縫發(fā)育特征和發(fā)育高度進(jìn)行了研究，建立了該礦的數(shù)值模擬模型，得出了導(dǎo)水裂縫帶發(fā)育高度。

2.2 FLAC

FLAC是1985年Cundall博士提出來(lái)的一種工程力學(xué)計(jì)算的有限差分程序，主要用于分析巖土工程[36-37]，該軟件具有多種結(jié)構(gòu)單元，能夠有效分析巖土及其他材料組成的復(fù)合結(jié)構(gòu)物，其內(nèi)置多種材料本構(gòu)模型[38]，能夠研究復(fù)雜工程的力學(xué)行為，得出各種云圖及曲線圖，具有很好的模擬效果[39]。

流固耦合機(jī)制是該模擬模型中比較復(fù)雜的內(nèi)容，其主要采用單相達(dá)西流進(jìn)行處理，其主要包括運(yùn)移定律、平衡定律、本構(gòu)法則、流體邊界條件4部分[40]。

對(duì)于流體的運(yùn)移主要采用達(dá)西定律進(jìn)行描述，即

(3)

流體的質(zhì)量遵循質(zhì)量平衡方程，即

?ξ/?t=-?qi/?xy+q0

(4)

式(4)中：ξ為流體質(zhì)量變化量；q0為初始滲流速度；xy為模型y軸方向長(zhǎng)度。

孔隙流體的反應(yīng)方程取決于飽和度的值，主要對(duì)完全飽和狀態(tài)下的本構(gòu)法則進(jìn)行分析。

?P/?t=M(?ξ/?t-α?ε/?t)

(5)

式(5)中：M代表比奧模量；ε代表體積應(yīng)變；α代表體積系數(shù)。

而流體的邊界初始條件，孔隙壓力為0，飽和度為1，因此透水邊界條件為

qn=h(p1-p2)

(6)

式(6)中：qn代表滲流速度分量；h代表巖土透水系數(shù)；p1代表邊界孔隙壓力；p2代表透水區(qū)域孔隙壓力。

張開(kāi)弦[41]采用FLAC對(duì)祁東煤礦開(kāi)采覆巖變形破壞規(guī)律進(jìn)行了研究。該煤礦生產(chǎn)能力2.4 Mt/a，存在較多含水層，最大厚度可達(dá)59.1 m，是最大額突水層，是礦井開(kāi)采主要威脅水源。

試驗(yàn)采用三維模擬模型，走向、傾向、高度分別是600、600、350 m。并在頂部施加一定載荷模擬覆巖自重。其力學(xué)參數(shù)如表2所示。

表2 FLAC數(shù)值模擬參數(shù)Table 2 Numerical simulation parameters of FLAC

根據(jù)力學(xué)參數(shù)和巖層厚度得出其計(jì)算模型如圖5所示。

圖5 FLAC數(shù)值模擬模型Fig.5 Simulation model of FLAC

共開(kāi)挖400 m，分40步并在模型兩邊設(shè)置保護(hù)煤柱。得出開(kāi)挖過(guò)程中覆巖位移特征如圖6所示。由圖6可知，當(dāng)工作面開(kāi)挖后，覆巖發(fā)生破壞，經(jīng)歷了垮落、開(kāi)裂、離層等幾個(gè)階段。煤層開(kāi)挖初始階段彎曲程度較小，隨著工作面的繼續(xù)開(kāi)挖，下沉繼續(xù)，隨著工作面距離的增加，覆巖下沉量越小，隨著工作天的繼續(xù)開(kāi)挖，周期性進(jìn)行垮落，當(dāng)推進(jìn)到200 m時(shí)，覆巖呈拱形，隨著繼續(xù)推移，當(dāng)工作面完全回采結(jié)束，覆巖云圖趨于完整。

圖6 FLAC模擬覆巖破壞特征Fig.6 Fault characteristics of overburden of FLAC

馮超等[42]采用FLAC對(duì)小保當(dāng)煤田煤層開(kāi)采引起的破壞進(jìn)行了研究，分析了導(dǎo)水裂縫帶發(fā)育高度；楊高峰等[43]采用該軟件對(duì)鎮(zhèn)城底煤礦8號(hào)煤層頂板兩帶發(fā)育高度進(jìn)行了研究，得出了覆巖導(dǎo)水裂縫帶發(fā)育高度，為礦井防治水提供了技術(shù)支撐；郝利生等[44]采用該軟件對(duì)斜溝煤礦18205工作面開(kāi)采覆巖兩帶破壞特征進(jìn)行了研究，得出了兩帶發(fā)育高度，并通過(guò)實(shí)測(cè)驗(yàn)證了模擬的準(zhǔn)確性；孫潤(rùn)等[45]采用該軟件對(duì)屯蘭煤礦工作面覆巖導(dǎo)水裂縫帶發(fā)育高度進(jìn)行了研究，確定了最大破壞高度。

2.3 RFPA

RFPA軟件主要應(yīng)用彈性力學(xué)來(lái)分析煤巖體破壞應(yīng)力特征[46-47]，在分析的過(guò)程中考慮了巖石本身的非均值性以及隨機(jī)性缺陷[48]。同時(shí)，模型中的單個(gè)單元對(duì)整體影響有限，當(dāng)模型中足夠多基元數(shù)量力學(xué)行為的疊加對(duì)模型整體產(chǎn)生影響[49]，就利用大量基元力學(xué)行為的集體效應(yīng)實(shí)現(xiàn)對(duì)巖石破壞失穩(wěn)過(guò)程的全程模擬[50]。

由于基元賦值符合Weibull分布，因此建立細(xì)觀與宏觀介質(zhì)力學(xué)聯(lián)系，其計(jì)算公式為

φ(α)=?/α0[(α/α0)?-1e-(α/α0)?]

(7)

式(7)中：α為煤巖體力學(xué)性質(zhì)；?為分布函數(shù)的性質(zhì)參數(shù)；α0為力學(xué)參數(shù)平均值；φ(α)為煤巖體介質(zhì)的力學(xué)統(tǒng)計(jì)分布密度。

假設(shè)模型的所有彈模平均值為E0，其分布函數(shù)積分為

(8)

式(8)中：φ(E)代表分布值；E代表彈性模量。由式(8)統(tǒng)計(jì)分布組成一個(gè)樣本空間，由于性質(zhì)參數(shù)的不同，得出的空間分布也不同。

段儉君[51]采用RFPA軟件對(duì)任家莊煤礦工作面開(kāi)采覆巖破壞高度進(jìn)行了研究。該煤礦主要開(kāi)采3、5號(hào)煤層，工作面頂板以泥巖、粉砂巖為主，采用走向長(zhǎng)壁一次采全高采煤方法。建立400 m×200 m的模型，其力學(xué)參數(shù)如表3所示。

表3 RFPA模擬力學(xué)參數(shù)Table 3 Mechanical parameters of RFPA

根據(jù)煤巖體力學(xué)性質(zhì)以及厚度進(jìn)行建模，得出其模型如圖7所示。

圖7 RFPA數(shù)值模擬模型Fig.7 Simulation model of RFPA

模擬煤層從左向右依次開(kāi)挖，開(kāi)挖長(zhǎng)度設(shè)置為120 m。共開(kāi)挖12步，每步10 m。得出工作面開(kāi)采過(guò)程中覆巖破壞特征如圖8所示。

由圖8可知，當(dāng)工作面推進(jìn)到30 m時(shí)，直接頂發(fā)生了垮落，隨著工作面的繼續(xù)推進(jìn)，當(dāng)推進(jìn)到50 m時(shí)，裂隙發(fā)育高度繼續(xù)向上部延展，并且老頂初次來(lái)壓，當(dāng)工作面推進(jìn)到80 m時(shí)，頂板周期性來(lái)壓，裂隙發(fā)育，裂隙帶發(fā)育高度達(dá)到35 m處。

圖8 RFPA模擬覆巖破壞特征Fig.8 Fault characteristics of overburden of RFPA

趙韶波[52]采用RFPA對(duì)寺河二號(hào)井覆巖破壞高度進(jìn)行了研究，確定了覆巖三帶發(fā)育高度，并用高位鉆孔對(duì)覆巖破壞高度進(jìn)行了驗(yàn)證；崔永青等[53]采用該軟件對(duì)馬堡煤礦工作面開(kāi)采覆巖三帶進(jìn)行了研究，得出垮落帶和裂隙帶破壞高度，與實(shí)測(cè)結(jié)果一直，驗(yàn)證了高技術(shù)的有效性；孫顯龍等[54]采用該軟件對(duì)煤礦頂板破壞特征進(jìn)行了研究，得出了工作面回采過(guò)程中覆巖破壞特征；顏林艷等[55]采用該技術(shù)對(duì)貴州煤礦開(kāi)采覆巖破壞特征進(jìn)行了分析，開(kāi)采承重覆巖重力形成應(yīng)力拱，構(gòu)造的存在是覆巖破壞的途徑。

2.4 PFC

PFC是顆粒流程序，它是由美國(guó)Itasca公司開(kāi)發(fā)的一種計(jì)算模擬軟件[56]，主要用于研究散粒體的分析，屬于離散元軟件范疇[57]。其理論基礎(chǔ)主要為牛頓第二定律和力與位移的關(guān)系，通過(guò)微觀角度來(lái)分析材料的宏觀行為[58]。當(dāng)顆粒所受到的不平衡力達(dá)到預(yù)設(shè)值時(shí)，模型則會(huì)停止計(jì)算，此時(shí)模擬結(jié)束[59]。

該顆粒流軟件主要包括物理方程和運(yùn)動(dòng)方程兩個(gè)方面[60]。

對(duì)于物理方程，主要分析顆粒之間的距離，其計(jì)算公式為

(9)

式(9)中：xA、xB代表兩顆粒的圓心位置。

兩顆粒之間的重疊量計(jì)算公式為

U=RA+RB-d

(10)

式(10)中：U代表顆粒之間重疊量；RA、RB代表兩顆粒的半徑。

對(duì)于運(yùn)動(dòng)方程，其加速度和角加速度計(jì)算公式為

ux(t0)=Fx/m

(11)

ωx(t0)=Mx/I

(12)

式中：ux代表計(jì)算度；Fx代表合力x方向的分量；m代表顆粒的質(zhì)量；I代表轉(zhuǎn)動(dòng)慣量；ωx代表角加速度；Mx代表顆粒受到力矩在x方向的分量；t0代表某一時(shí)刻。

陸敬鋒[61]采用PFC軟件對(duì)新疆大南湖礦區(qū)采動(dòng)覆巖運(yùn)移規(guī)律進(jìn)行了研究。建立模型長(zhǎng)、高分別是180、100 m。模型頂板施加2.5 MPa的應(yīng)力，其力學(xué)參數(shù)如表4所示。

表4 PFC模擬力學(xué)參數(shù)Table 4 Mechanical parameters of PFC

根據(jù)煤巖體的厚度以及力學(xué)參數(shù)建立數(shù)值模型如圖9所示。

圖9 PFC數(shù)值模擬模型Fig.9 Simulation model of PFC

將巖層分為20個(gè)巖層，煤層設(shè)置高度為3 m，每次開(kāi)挖10 m，共開(kāi)挖10次，開(kāi)采長(zhǎng)度為100 m。開(kāi)挖過(guò)程中覆巖破壞特征如圖10所示。

由圖10可知，當(dāng)工作面推進(jìn)到10 m時(shí)，工作面頂板發(fā)生下沉現(xiàn)象，直接頂出現(xiàn)了裂縫，直接頂沒(méi)有垮落；當(dāng)工作面推進(jìn)到50 m時(shí)，由于關(guān)鍵層的存在，能夠承受較大的載荷，抑制頂板下沉，垮落帶高度為6.56 m；當(dāng)工作面推進(jìn)到100 m時(shí)，裂隙向深部發(fā)展，含水層貫通，導(dǎo)水裂縫帶發(fā)育到下部含水層，上部隔水層沒(méi)有明顯破壞，未出現(xiàn)貫通裂隙。

圖10 PFC模擬覆巖破壞特征Fig.10 Fault characteristics of overburden of PFC

江成浩等[62]采用PFC軟件對(duì)滕東煤礦綜放工作面裂隙演化規(guī)律進(jìn)行了研究，得出了覆巖孔隙率；鐘江城等[63]采用該軟件對(duì)淺埋煤層覆巖裂隙演化進(jìn)行了研究，裂隙開(kāi)度與傾角是影響礦井突水突砂的主要影響因素；羅懷廷等[64]采用該軟件對(duì)固體充填采煤中的覆巖破壞特征進(jìn)行了研究，破碎矸石在側(cè)限壓縮條件下所產(chǎn)生的側(cè)壓力與其所受軸壓成正比，粒徑越小顆粒轉(zhuǎn)動(dòng)越活躍，大顆粒自身具備運(yùn)動(dòng)惰性。

2.5 CDEM

CDEM軟件是由中國(guó)科學(xué)院和北京極道成然科技有限公司聯(lián)合開(kāi)發(fā)的力學(xué)分析系列軟件，該軟件是基于連續(xù)介質(zhì)的力學(xué)離散方法[65-66]，主要用于模擬材料的變形和破壞過(guò)程，軟件單元之間通過(guò)虛擬彈簧傳遞相互作用力，在計(jì)算的過(guò)程中塊體單元可以實(shí)現(xiàn)斷裂[67]。模擬材料的破壞時(shí)通過(guò)界面的張開(kāi)、滑移和運(yùn)動(dòng)速度等表述，非常適合描述地質(zhì)體由連續(xù)向非連續(xù)的破壞過(guò)程[68]。

該軟件采用剛度矩陣法進(jìn)行求解，在每個(gè)迭代步通過(guò)式(13)計(jì)算單元自身的節(jié)點(diǎn)力，并將此節(jié)點(diǎn)力分配至單元對(duì)應(yīng)的節(jié)點(diǎn)上。

Fi=KiuiC

(13)

式(13)中：Fi代表單元i的節(jié)點(diǎn)力向量；Ki代表單元i的節(jié)點(diǎn)位移向量；ui代表單元i的單元?jiǎng)偠染仃?；C代表黏聚力。

其界面的彈簧力計(jì)算公式為

(14)

式(14)中：Fn、Fs分別代表彈簧的法向、切向應(yīng)力；Kn、Ks分別代表彈簧的法向、切向剛度；Δdn、Δds分別代表彈簧的法向、切向位移。

當(dāng)破壞計(jì)算時(shí)，采用莫爾庫(kù)倫準(zhǔn)則進(jìn)行修正，修正的計(jì)算式為

(15)

式(15)中：T代表抗拉強(qiáng)度；φ代表內(nèi)摩擦角。

焦振華[69]采用CDEM軟件對(duì)寺河二號(hào)井保護(hù)層開(kāi)采覆巖破壞特征進(jìn)行了研究。試驗(yàn)中采用模型長(zhǎng)和高分別是200 m×72.2 m，其力學(xué)參數(shù)如表5所示。

表5 CDEM模擬力學(xué)參數(shù)Table 5 Mechanical parameters of CDEM

根據(jù)煤巖體厚度及力學(xué)性質(zhì)建立相關(guān)模型如圖11所示。

圖11 CDEM模擬模型Fig.11 Simulation model of CDEM

為了消除邊界效應(yīng)，在模型的兩側(cè)留設(shè)50 m的煤柱，并在模型頂部施加7.5 MPa的巖層自重力，得出工作面回采中覆巖破壞特征如圖12所示。

圖12 CDEM數(shù)值模擬覆巖破壞特征Fig.12 Fault characteristics of overburden of CDEM

由圖12可知，當(dāng)工作面回采到20 m時(shí)，工作面頂板發(fā)生初次垮落，隨著工作面的繼續(xù)推進(jìn)，直接頂隨著采動(dòng)而發(fā)生冒落，當(dāng)工作面推進(jìn)到30 m時(shí)，基本頂發(fā)生垮落。最終工作面導(dǎo)水裂縫帶發(fā)育高度為30 m。

王振偉等[70]采用該方法對(duì)安家?guī)X煤礦覆巖破壞進(jìn)行了研究，得出了彎曲帶與垮落帶和裂隙帶的厚度比為1∶2，初次垮落步距為115～165 m；郝應(yīng)祥[71]采用該軟件對(duì)孫家溝煤礦覆巖開(kāi)采對(duì)地表的影響進(jìn)行了研究，得出了頂板破壞規(guī)律和位移的變化特征，確定了地表沉降范圍。

2.6 ABQUES

ABQUES是基于有限元的軟件，其主要研究的是相對(duì)簡(jiǎn)單的線性以及眾多復(fù)雜的非線性問(wèn)題?？梢院芎玫啬M各種形狀的單元庫(kù)以及眾多類(lèi)型材料的性能[72-73]。其適合沖擊載荷和爆炸、地震這種短暫的瞬時(shí)的動(dòng)態(tài)事件的分析與研究，對(duì)于高度非線性問(wèn)題也非常有效，包括非線性大變形等問(wèn)題的求解研究[74]。

其工作面開(kāi)采覆巖變化的黏聚力模型服從Drucker-Prager屈服條件，即

(16)

式(16)中：I1、J2分別代表應(yīng)力張量第一和第二不變量；α、k代表關(guān)于內(nèi)摩擦角和黏聚力的參數(shù)。

砌體界面單元的本構(gòu)關(guān)系為

t=Eε

(17)

式(17)中：t代表應(yīng)力向量；E代表彈性模量；ε代表名義應(yīng)變向量。

黏結(jié)單元指數(shù)計(jì)算公式為

(18)

式(18)中：D代表?yè)p傷變量；δ0代表初始的有效位移量；δmax代表加載過(guò)程中的最大位移量；δf代表失效時(shí)的有效位移量。

田森等[75-76]采用ABQUES對(duì)保護(hù)層開(kāi)采覆巖裂隙演化規(guī)律進(jìn)行了研究。該煤礦主采1號(hào)煤層，煤層埋深可達(dá)1 100 m，開(kāi)采區(qū)域長(zhǎng)度為90 m，其力學(xué)參數(shù)如表6所示。

表6 ABQUES數(shù)值模擬力學(xué)參數(shù)Table 6 Mechanical parameters of ABQUES

根據(jù)煤巖體的力學(xué)性質(zhì)和厚度建立煤層采動(dòng)有限元模型如圖13所示。

圖13 ABQUES數(shù)值模擬模型Fig.13 Simulation model of ABQUES

開(kāi)采過(guò)程中，通過(guò)殺死單元的方式進(jìn)行煤層的開(kāi)挖，煤層開(kāi)采共分9步進(jìn)行，每步開(kāi)挖90 m，共開(kāi)挖90 m，覆巖變化特征如圖14所示。

圖14 ABQUES模擬覆巖破壞特征Fig.14 Fault characteristics of overburden of ABQUES

由圖14可知，當(dāng)工作面推進(jìn)到30 m時(shí)，覆巖應(yīng)力發(fā)生了明顯變化，且靠近采空區(qū)中心應(yīng)力變化最為顯著，也是卸壓最快的區(qū)域；當(dāng)工作面回采到60 m時(shí)，覆巖應(yīng)力繼續(xù)發(fā)生著變化，同樣，采空區(qū)附近應(yīng)力變化最為明顯，形成了駝峰狀變化曲線；當(dāng)工作面推進(jìn)到90 m時(shí)，覆巖應(yīng)力重新分布，卸壓區(qū)范圍逐漸變大，采空區(qū)附近是最大，在支撐應(yīng)力區(qū)域內(nèi)卸壓效果較差。

王遷等[77]采用ABQUES軟件對(duì)司馬礦1112工作面煤層開(kāi)采覆巖變化特征進(jìn)行了研究，得出了下沉特征以及不同推進(jìn)距離下與下沉量的關(guān)系和工作面覆巖應(yīng)力變化特征。

3 存在問(wèn)題及發(fā)展趨勢(shì)

3.1 存在問(wèn)題

(1)對(duì)于工作面回采覆巖破壞特征的數(shù)值模擬研究，多以平面模型為主。而采礦工程是一個(gè)三維模型，其工作面的開(kāi)采覆巖變化特征是一個(gè)時(shí)空演變的復(fù)雜過(guò)程。單純的平面模型難以準(zhǔn)確地模擬覆巖破壞特征，難以反映開(kāi)采覆巖破壞的實(shí)質(zhì)問(wèn)題。

(2)與簡(jiǎn)單的巖土工程對(duì)比，工作面回采過(guò)程中覆巖移動(dòng)、損傷、破壞、失穩(wěn)是必然要發(fā)生的。目前的模擬軟件只能對(duì)理想的彈性、塑性以及損傷體進(jìn)行有效的變形和受力進(jìn)行分析。而采礦工程為更重要的是材料和結(jié)構(gòu)破壞后的力學(xué)特性、煤巖體是如何破壞以及破壞穩(wěn)定的變化過(guò)程。

(3)數(shù)值模擬軟件在模擬過(guò)程中會(huì)對(duì)邊界條件以及材料屬性進(jìn)行簡(jiǎn)化，多對(duì)模擬結(jié)果產(chǎn)生一定影響，并且模擬結(jié)構(gòu)離散形式不同，模擬的結(jié)果和精度也不相同，隨機(jī)性較大，造成模擬結(jié)果的可信度大大降低[78-80]。

3.2 發(fā)展趨勢(shì)

(1)隨著中國(guó)煤炭開(kāi)采技術(shù)的不斷發(fā)展以及環(huán)保和可持續(xù)發(fā)展理念的加強(qiáng)，提出了綠色開(kāi)采和智能化開(kāi)采，這些技術(shù)屬于較前沿課題，缺乏一定的技術(shù)支持。未來(lái)數(shù)值模擬應(yīng)該結(jié)合綠色開(kāi)采和智能化開(kāi)采技術(shù)研發(fā)新型的模擬模塊。適合煤礦技術(shù)的發(fā)展。

(2)采礦覆巖破壞是一個(gè)時(shí)空變化過(guò)程，平面模型不足以模擬開(kāi)采過(guò)程。未來(lái)應(yīng)研發(fā)可視化的三維模擬模型，高強(qiáng)度還原采礦過(guò)程中覆巖破壞特征，提高模擬準(zhǔn)確度。

(3)數(shù)值模擬僅僅考慮開(kāi)采深度和采厚因素，對(duì)于其他地質(zhì)、水文因素幾乎無(wú)涉及，造成模擬準(zhǔn)確度不高。未來(lái)數(shù)值模擬應(yīng)綜合覆巖破壞多影響因素，優(yōu)化煤巖體力學(xué)參數(shù)合理確定方法，引入新的思維方法和數(shù)學(xué)工具，提高模擬結(jié)果的準(zhǔn)確度。

(4)隨著煤礦開(kāi)采深度的不斷增加，對(duì)于軟巖在高地應(yīng)力蠕變特性突出。在現(xiàn)有數(shù)值模擬軟件中未考慮該問(wèn)題。未來(lái)數(shù)值模擬應(yīng)考慮高地應(yīng)力下軟巖蠕變特性，對(duì)于保證采礦工程的長(zhǎng)期穩(wěn)定性至關(guān)重要。