煤與瓦斯突出多物理場分布特征的數值模擬研究

曹 偈，趙旭生，劉延保

(1.瓦斯災害監控與應急技術國家重點實驗室，重慶 400037； 2.中煤科工集團重慶研究院有限公司，重慶 400037)

2019年我國煤炭消費量占能源消費總量的57.7%[1]，煤炭在未來相當長一段時間內還將作為主要能源存在。我國是世界上煤與瓦斯突出災害最為嚴重的國家，近年來在大量關停突出煤礦,以及煤與瓦斯突出防治技術不斷創新的背景下，煤與瓦斯突出事故呈現下降趨勢。但由于煤礦井下條件的復雜性及煤與瓦斯突出防治技術的局限性，突出事故仍是影響煤礦安全生產的主要災害之一。

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

綜上所述，目前學者們針對其研究重點，采用不同數值分析軟件對煤與瓦斯突出過程的各物理場分布規律進行了初步探索。筆者在前期開展的煤與瓦斯突出物理模擬試驗基礎上，利用COMSOL模擬軟件對含瓦斯煤的氣固耦合方程進行解算，分析突出后煤體應力—變形—瓦斯場的分布規律，為突出機理的量化分析提供基礎數據。

1 煤與瓦斯突出數值模擬

1.1 含瓦斯煤氣固耦合方程

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

1)應力場控制方程

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

(1)

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

本研究中采用Drucker-Prager匹配Mohr-Coulomb準則來表征煤體的塑性失穩，其關系式為：

(2)

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

2)氣體流動場控制方程

煤中氣體流動場控制方程由氣體狀態方程、氣體含量方程、氣體運動方程及連續性方程構成。

①氣體狀態方程

(3)

式中：ρg為氣體密度；ρn為標準狀態下的氣體密度；pn為標準狀態下的氣體壓力；z、zn為氣體壓縮因子。

②氣體含量方程

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

(4)

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

③氣體運動方程

煤中由于裂隙發育的離散性導致氣體流動形態不同，可根據當地雷諾數決定，其運動可用廣義冪定律描述：

(5)

式中:v為氣體滲流速度；k為滲透率；μ為氣體動力黏度；m為氣體狀態指數，m=1～2，當m=1時，式(5)為達西滲流模型。

由于含瓦斯煤巖在應力環境中，應力的改變導致孔隙率及滲透率呈現動態變化?？紫堵视嬎愎饺缦拢?/p>

(6)

(7)

式中k0為初始滲透率。

④氣體運動的連續性方程

(8)

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

1.2 幾何模型及模型參數

研究團隊前期利用自主研制的大型煤與瓦斯突出模擬試驗裝置開展了垂直應力16 MPa、氣體壓力0.5 MPa條件下的突出模擬試驗[20]。試驗模型尺寸為1 500 mm×800 mm×800 mm，上覆200 mm厚的剛性壓板。突出煤粉15.86 kg，形成深10 cm、寬30 cm左右的突出孔洞?；诖舜瓮怀鑫锢砟M試驗構建了二維數值分析模型，模型煤層長1 500 mm、高800 mm，上覆剛性壓板長1 490 mm、高200 mm，將突出口孔洞簡化為半徑分別為150、100 mm的半橢圓形。突出幾何模型及網絡劃分如圖1所示。

圖1 突出幾何模型及網格劃分

根據物理試驗條件，幾何模型上部施加16 MPa垂直應力，左右兩邊及底部采用輥支承邊界，突出后突出孔洞內(圖1紅線所示)為自由邊界。突出前，煤層內部初始氣壓p0=0.5 MPa，煤層四周邊界均為零通量不透氣邊界(dp/dxi=0)；突出后，突出孔洞邊界為大氣壓。數值模擬以突出物理模擬試驗模型中取樣測試得到的煤層物性參數[20]為基礎，如表1所示。模擬分為2個步驟，首先利用固體力學模塊計算得到突出前煤層應力狀態；然后利用固體力學與PDE模塊計算突出后煤層應力、變形，以及氣體流動規律。

表1 模型物性參數值

2 模擬結果分析

2.1 數值模擬結果驗證

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

圖2 突出前后煤層應力分布的對比圖

由圖2可以看出，由于數值模型的簡化，數值模擬結果與試驗結果有一定差異。但數值模擬結果與物理試驗結果具有一致的變化規律，突出發生后在距突出口100～350 mm區域應力呈現降低趨勢，距突出口350～1 150 mm區域應力呈增大趨勢。證實了數值模型可用于對突出前后煤層物理場分布規律進行分析。

2.2 突出后各物理因素的分布規律

1)應力變形規律

突出后煤層應力、體積應變、塑性區分布云圖如圖3所示。突出發生前由于軸向應力施加在剛性壓板中心700 mm范圍內，形成了煤體上表面實際施加了非均勻分布力的特點，模擬了實際工作面煤體應力集中現象。

(a)垂直應力

(b)體積應變

(c)塑性區

由圖3可以得出，突出后煤體應力達到新的平衡狀態?？锥幢谏舷旅嬗捎谕怀鲞^程煤體被拋出而卸壓，應力降低；孔洞壁后方形成應力增高區，此處煤體體積應變增大，產生塑性變形，破壞較為嚴重。而到煤體深部受突出影響較小，煤體體積應變并不顯著。

2)氣體流動規律

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

(a)滲透率

(b)氣體壓力

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

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

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

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

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

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

由圖6可以看出，不同應力條件下煤體滲透率受突出影響程度不同。當應力較小時煤體塑性變形不顯著，滲透率變化較小；隨著應力增大，突出孔洞壁處煤體滲透率變化明顯，當應力大于 16 MPa 時，比煤體原始滲透率增加2～3個數量級。由此，對于高應力環境突出，較大的滲透率變化導致煤體氣體流動及壓力衰減較快，突出發生后在極短的時間內氣體壓力梯度減小，不具備突出持續發生的條件，突出煤量相對較小。

3 結論

利用Comsol Multiphysics中的固體力學和PDE自定義模塊構建了煤與瓦斯突出的數值分析模型，獲得了突出后煤體應力場、變形場、氣體流動場分布規律，并與前期開展的煤與瓦斯突出物理模擬試驗結果進行對比驗證，證實了數值分析方法的可行性。分析得到突出后孔洞壁附近煤體發生塑性變形，導致此處滲透率增大，且隨著應力增大，突出孔洞壁處煤體滲透率可比原始滲透率增加2～3個數量級。而氣體壓力從突出孔洞壁面到煤體深部呈現先急劇上升后降低至原始氣壓的變化趨勢，在孔洞壁處形成較高的壓力梯度，這是突出持續發生的必要條件。