煤與瓦斯突出多物理場(chǎng)分布特征的數(shù)值模擬研究

曹 偈，趙旭生，劉延保

(1.瓦斯災(zāi)害監(jiān)控與應(yīng)急技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400037； 2.中煤科工集團(tuán)重慶研究院有限公司，重慶 400037)

2019年我國(guó)煤炭消費(fèi)量占能源消費(fèi)總量的57.7%[1]，煤炭在未來(lái)相當(dāng)長(zhǎng)一段時(shí)間內(nèi)還將作為主要能源存在。我國(guó)是世界上煤與瓦斯突出災(zāi)害最為嚴(yán)重的國(guó)家，近年來(lái)在大量關(guān)停突出煤礦,以及煤與瓦斯突出防治技術(shù)不斷創(chuàng)新的背景下，煤與瓦斯突出事故呈現(xiàn)下降趨勢(shì)。但由于煤礦井下條件的復(fù)雜性及煤與瓦斯突出防治技術(shù)的局限性，突出事故仍是影響煤礦安全生產(chǎn)的主要災(zāi)害之一。

目前，學(xué)者們普遍認(rèn)可煤與瓦斯突出是地應(yīng)力、瓦斯、煤巖物理力學(xué)性質(zhì)綜合作用的結(jié)果，并利用自主研制的大型煤與瓦斯突出模擬裝置開(kāi)展物理相似模擬試驗(yàn)[2-8]，對(duì)突出機(jī)理進(jìn)行進(jìn)一步的量化分析，為災(zāi)害的預(yù)防控制提供了基礎(chǔ)數(shù)據(jù)。但由于開(kāi)展大型突出試驗(yàn)的難度較大且測(cè)試參數(shù)有限，在試驗(yàn)基礎(chǔ)上利用數(shù)值分析方法對(duì)突出過(guò)程煤與瓦斯等參數(shù)的變化規(guī)律進(jìn)行研究具有重要意義。徐濤[9]、王路軍[10]、顏愛(ài)華[11]、劉超[12]等科技工作者利用RFPA數(shù)值分析軟件研究了突出過(guò)程應(yīng)力場(chǎng)、滲流場(chǎng)等的變化規(guī)律；劉雪琴等[13]構(gòu)建了突出過(guò)程含瓦斯煤壓力場(chǎng)控制方程，以陽(yáng)煤五礦趙家分區(qū)煤礦為例，利用Fluent軟件分析了掘進(jìn)過(guò)程中瓦斯壓力、滲透率演變規(guī)律；崔聰[14]利用Abaqus軟件構(gòu)建三維有限元模型，分析了不同水平應(yīng)力條件下突出洞穴形體的孕育特征，并應(yīng)用平煤八礦突出現(xiàn)場(chǎng)洞穴特征與數(shù)值模擬結(jié)果進(jìn)行分析比較，證實(shí)了該方法的準(zhǔn)確性；劉永立[15]、王銳[16]等利用顆粒流理論對(duì)煤與瓦斯突出的臨界狀態(tài)、突出顆粒運(yùn)動(dòng)軌跡等進(jìn)行了模擬分析；藍(lán)航等[17]建立了煤巖動(dòng)力災(zāi)害統(tǒng)一能量方程，并利用FLAC模擬軟件計(jì)算分析了不同初始瓦斯壓力下圍巖能量及瓦斯涌出量變化，認(rèn)為瓦斯的不同參與程度導(dǎo)致了不同的動(dòng)力災(zāi)害顯現(xiàn)特征；AN F H等[18]建立了突出發(fā)生時(shí)考慮煤損傷的瓦斯膨脹能計(jì)算方法，利用COMSOL軟件模擬分析了不同氣壓、滲透率及加載條件下的突出瓦斯膨脹能。

綜上所述，目前學(xué)者們針對(duì)其研究重點(diǎn)，采用不同數(shù)值分析軟件對(duì)煤與瓦斯突出過(guò)程的各物理場(chǎng)分布規(guī)律進(jìn)行了初步探索。筆者在前期開(kāi)展的煤與瓦斯突出物理模擬試驗(yàn)基礎(chǔ)上，利用COMSOL模擬軟件對(duì)含瓦斯煤的氣固耦合方程進(jìn)行解算，分析突出后煤體應(yīng)力—變形—瓦斯場(chǎng)的分布規(guī)律，為突出機(jī)理的量化分析提供基礎(chǔ)數(shù)據(jù)。

1 煤與瓦斯突出數(shù)值模擬

1.1 含瓦斯煤氣固耦合方程

煤與瓦斯突出是一個(gè)復(fù)雜的力學(xué)作用過(guò)程[1]，在開(kāi)展數(shù)值分析過(guò)程中提出基本假設(shè)及簡(jiǎn)化：①煤體為各向同性的彈塑性多孔介質(zhì)；②不考慮煤巖變形破壞及氣體解吸滲流過(guò)程中的溫度因素；③煤與瓦斯突出后突出孔洞壁后方煤體變形為小變形；④氣體滲流服從廣義冪定律。

1)應(yīng)力場(chǎng)控制方程

由前人研究可知，考慮氣體壓力及吸附膨脹變形影響的煤體平衡微分方程[19]如下：

(1)

式中:G為切變模量；ui為煤體位移；ν為泊松比；α為有效應(yīng)力系數(shù)；p為氣體壓力；δij為Kronecker參數(shù)；E為煤骨架彈性模量；εs為氣體吸附引起的骨架變形；fi為體積力。

本研究中采用Drucker-Prager匹配Mohr-Coulomb準(zhǔn)則來(lái)表征煤體的塑性失穩(wěn)，其關(guān)系式為：

(2)

式中：I1、I2為應(yīng)力張量第一、第二不變量；J2為偏應(yīng)力張量第二不變量；C為黏聚力；φ為內(nèi)摩擦角；σi為主應(yīng)力分量，i=1，2，3。

2)氣體流動(dòng)場(chǎng)控制方程

煤中氣體流動(dòng)場(chǎng)控制方程由氣體狀態(tài)方程、氣體含量方程、氣體運(yùn)動(dòng)方程及連續(xù)性方程構(gòu)成。

①氣體狀態(tài)方程

(3)

式中：ρg為氣體密度；ρn為標(biāo)準(zhǔn)狀態(tài)下的氣體密度；pn為標(biāo)準(zhǔn)狀態(tài)下的氣體壓力；z、zn為氣體壓縮因子。

②氣體含量方程

煤中甲烷主要以游離態(tài)與吸附態(tài)存在，其中吸附瓦斯含量占90%以上。瓦斯吸附解吸滿足朗格繆爾方程，則煤中氣體含量計(jì)算公式為：

(4)

式中：M為煤體中的氣體含量；Ma為吸附氣體含量；Mb為游離氣體含量；a、b為吸附常數(shù)；η為煤質(zhì)校正系數(shù)；q為孔隙率。

③氣體運(yùn)動(dòng)方程

煤中由于裂隙發(fā)育的離散性導(dǎo)致氣體流動(dòng)形態(tài)不同，可根據(jù)當(dāng)?shù)乩字Z數(shù)決定，其運(yùn)動(dòng)可用廣義冪定律描述：

(5)

式中:v為氣體滲流速度；k為滲透率；μ為氣體動(dòng)力黏度；m為氣體狀態(tài)指數(shù)，m=1～2，當(dāng)m=1時(shí)，式(5)為達(dá)西滲流模型。

由于含瓦斯煤巖在應(yīng)力環(huán)境中，應(yīng)力的改變導(dǎo)致孔隙率及滲透率呈現(xiàn)動(dòng)態(tài)變化?？紫堵视?jì)算公式如下：

(6)

(7)

式中k0為初始滲透率。

④氣體運(yùn)動(dòng)的連續(xù)性方程

(8)

將公式(3)～(8)組合變換成偏微分形式，利用Comsol Multiphysics中的固體力學(xué)和PDE自定義模塊，對(duì)方程進(jìn)行數(shù)值求解，模擬突出煤層損傷與氣體滲流的耦合規(guī)律，并與前期開(kāi)展的物理模擬試驗(yàn)結(jié)果相互驗(yàn)證。

1.2 幾何模型及模型參數(shù)

研究團(tuán)隊(duì)前期利用自主研制的大型煤與瓦斯突出模擬試驗(yàn)裝置開(kāi)展了垂直應(yīng)力16 MPa、氣體壓力0.5 MPa條件下的突出模擬試驗(yàn)[20]。試驗(yàn)?zāi)Ｐ统叽鐬? 500 mm×800 mm×800 mm，上覆200 mm厚的剛性壓板。突出煤粉15.86 kg，形成深10 cm、寬30 cm左右的突出孔洞?；诖舜瓮怀鑫锢砟M試驗(yàn)構(gòu)建了二維數(shù)值分析模型，模型煤層長(zhǎng)1 500 mm、高800 mm，上覆剛性壓板長(zhǎng)1 490 mm、高200 mm，將突出口孔洞簡(jiǎn)化為半徑分別為150、100 mm的半橢圓形。突出幾何模型及網(wǎng)絡(luò)劃分如圖1所示。

圖1 突出幾何模型及網(wǎng)格劃分

根據(jù)物理試驗(yàn)條件，幾何模型上部施加16 MPa垂直應(yīng)力，左右兩邊及底部采用輥支承邊界，突出后突出孔洞內(nèi)(圖1紅線所示)為自由邊界。突出前，煤層內(nèi)部初始?xì)鈮簆0=0.5 MPa，煤層四周邊界均為零通量不透氣邊界(dp/dxi=0)；突出后，突出孔洞邊界為大氣壓。數(shù)值模擬以突出物理模擬試驗(yàn)?zāi)Ｐ椭腥訙y(cè)試得到的煤層物性參數(shù)[20]為基礎(chǔ)，如表1所示。模擬分為2個(gè)步驟，首先利用固體力學(xué)模塊計(jì)算得到突出前煤層應(yīng)力狀態(tài)；然后利用固體力學(xué)與PDE模塊計(jì)算突出后煤層應(yīng)力、變形，以及氣體流動(dòng)規(guī)律。

表1 模型物性參數(shù)值

2 模擬結(jié)果分析

2.1 數(shù)值模擬結(jié)果驗(yàn)證

前期開(kāi)展的物理模擬試驗(yàn)為了采集突出過(guò)程軸向應(yīng)力變化，在煤層垂直高度600 mm的水平面中心線上距突出口100、350、550、750 mm位置布置有應(yīng)力傳感器。由試驗(yàn)結(jié)果可知，在距突出口100 mm處的煤層應(yīng)力由于下方煤體破碎并拋出，形成突出孔洞，其應(yīng)力呈下降趨勢(shì)。而孔洞壁后方煤體由于應(yīng)力重新分布，垂直應(yīng)力有所升高，且越靠近孔洞壁處的應(yīng)力變化越大。為驗(yàn)證數(shù)值分析結(jié)果的可靠性，在模型相同位置設(shè)置觀測(cè)線。數(shù)值模擬與物理試驗(yàn)的突出前后煤層應(yīng)力分布對(duì)比如圖2所示。

圖2 突出前后煤層應(yīng)力分布的對(duì)比圖

由圖2可以看出，由于數(shù)值模型的簡(jiǎn)化，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果有一定差異。但數(shù)值模擬結(jié)果與物理試驗(yàn)結(jié)果具有一致的變化規(guī)律，突出發(fā)生后在距突出口100～350 mm區(qū)域應(yīng)力呈現(xiàn)降低趨勢(shì)，距突出口350～1 150 mm區(qū)域應(yīng)力呈增大趨勢(shì)。證實(shí)了數(shù)值模型可用于對(duì)突出前后煤層物理場(chǎng)分布規(guī)律進(jìn)行分析。

2.2 突出后各物理因素的分布規(guī)律

1)應(yīng)力變形規(guī)律

突出后煤層應(yīng)力、體積應(yīng)變、塑性區(qū)分布云圖如圖3所示。突出發(fā)生前由于軸向應(yīng)力施加在剛性壓板中心700 mm范圍內(nèi)，形成了煤體上表面實(shí)際施加了非均勻分布力的特點(diǎn)，模擬了實(shí)際工作面煤體應(yīng)力集中現(xiàn)象。

(a)垂直應(yīng)力

(b)體積應(yīng)變

(c)塑性區(qū)

由圖3可以得出，突出后煤體應(yīng)力達(dá)到新的平衡狀態(tài)?？锥幢谏舷旅嬗捎谕怀鲞^(guò)程煤體被拋出而卸壓，應(yīng)力降低；孔洞壁后方形成應(yīng)力增高區(qū)，此處煤體體積應(yīng)變?cè)龃?，產(chǎn)生塑性變形，破壞較為嚴(yán)重。而到煤體深部受突出影響較小，煤體體積應(yīng)變并不顯著。

2)氣體流動(dòng)規(guī)律

突出后煤體滲透率及氣體壓力分布云圖如圖4所示。

(a)滲透率

(b)氣體壓力

由圖4(a)可知，由于突出孔洞壁后方煤體破壞嚴(yán)重，其滲透率在此處達(dá)到最高，從孔洞壁面到煤體深部呈減小趨勢(shì)，在體積應(yīng)變最大值處達(dá)到最小，隨后恢復(fù)至原始滲透率。

由圖4(b)可知，煤體內(nèi)氣體壓力從孔洞壁面到煤體深部急劇上升，在應(yīng)力集中處達(dá)到最大值，且高于煤體原始?xì)鈮海S后恢復(fù)至原始?xì)鈮?，如圖5所示。這是由于應(yīng)力集中使煤體被進(jìn)一步壓實(shí)，滲透性降低，氣體難以向外流動(dòng)，使瓦斯在短時(shí)間內(nèi)迅速積聚，在突出孔洞壁面附近產(chǎn)生較高的壓力梯度。若此時(shí)突出壁面附近煤體破碎嚴(yán)重，且產(chǎn)生的氣體壓力梯度具有較高能量將煤體拋出，突出則會(huì)持續(xù)發(fā)生。

圖5 隨著時(shí)間推移突出口中軸線上氣體壓力分布曲線

由圖5還可以看出，隨著時(shí)間推移，突出孔洞壁面后方氣體壓力逐漸降低，氣體壓力梯度減小，5 s時(shí)間內(nèi)氣壓峰值消失，說(shuō)明突出的持續(xù)發(fā)生具有時(shí)效性。同一位置處，氣體壓力隨時(shí)間的衰減速率也呈現(xiàn)降低趨勢(shì)。

不同垂直應(yīng)力條件下突出口中軸線上煤體滲透率的分布曲線如圖6所示。

圖6 不同垂直應(yīng)力條件下的煤體滲透率分布曲線

由圖6可以看出，不同應(yīng)力條件下煤體滲透率受突出影響程度不同。當(dāng)應(yīng)力較小時(shí)煤體塑性變形不顯著，滲透率變化較??；隨著應(yīng)力增大，突出孔洞壁處煤體滲透率變化明顯，當(dāng)應(yīng)力大于 16 MPa 時(shí)，比煤體原始滲透率增加2～3個(gè)數(shù)量級(jí)。由此，對(duì)于高應(yīng)力環(huán)境突出，較大的滲透率變化導(dǎo)致煤體氣體流動(dòng)及壓力衰減較快，突出發(fā)生后在極短的時(shí)間內(nèi)氣體壓力梯度減小，不具備突出持續(xù)發(fā)生的條件，突出煤量相對(duì)較小。

3 結(jié)論

利用Comsol Multiphysics中的固體力學(xué)和PDE自定義模塊構(gòu)建了煤與瓦斯突出的數(shù)值分析模型，獲得了突出后煤體應(yīng)力場(chǎng)、變形場(chǎng)、氣體流動(dòng)場(chǎng)分布規(guī)律，并與前期開(kāi)展的煤與瓦斯突出物理模擬試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證，證實(shí)了數(shù)值分析方法的可行性。分析得到突出后孔洞壁附近煤體發(fā)生塑性變形，導(dǎo)致此處滲透率增大，且隨著應(yīng)力增大，突出孔洞壁處煤體滲透率可比原始滲透率增加2～3個(gè)數(shù)量級(jí)。而氣體壓力從突出孔洞壁面到煤體深部呈現(xiàn)先急劇上升后降低至原始?xì)鈮旱淖兓厔?shì)，在孔洞壁處形成較高的壓力梯度，這是突出持續(xù)發(fā)生的必要條件。