呂 閆 孫 亮

(1.山西華晉焦煤有限責(zé)任公司沙曲二礦;2.煤科集團沈陽研究院有限公司)

?

基于MATLAB的煤礦通風(fēng)與瓦斯流動一體化解算

呂 閆1孫 亮2

(1.山西華晉焦煤有限責(zé)任公司沙曲二礦;2.煤科集團沈陽研究院有限公司)

為快速地獲得某礦301綜放工作面通風(fēng)巷道內(nèi)通風(fēng)風(fēng)量與瓦斯流量網(wǎng)絡(luò)解算數(shù)據(jù)，在對工作面通風(fēng)系統(tǒng)參數(shù)測試基礎(chǔ)上，使用CAD軟件提取巷道坐標(biāo)數(shù)據(jù)，運用MATLAB軟件，結(jié)合圖論理論、流體力學(xué)理論，編制解算程序，程序可讀取巷道坐標(biāo)并且自動建立工作面網(wǎng)絡(luò)拓?fù)潢P(guān)系，并求解工作面通風(fēng)和瓦斯網(wǎng)絡(luò)解算各參數(shù)矩陣，利用Cross迭代算法解算工作面通風(fēng)系統(tǒng)內(nèi)的通風(fēng)巷道通風(fēng)風(fēng)量和瓦斯流量解算數(shù)據(jù)。通過與現(xiàn)場實測數(shù)據(jù)對比，程序解算數(shù)據(jù)基本正確。模擬結(jié)果表明，MATLAB軟件能快速便捷地處理工作面通風(fēng)風(fēng)量和瓦斯流量網(wǎng)絡(luò)解算，并具有較高的準(zhǔn)確性。

綜放工作面 通風(fēng) 瓦斯 MATLAB 網(wǎng)絡(luò)解算 Cross迭代算法

煤礦通風(fēng)風(fēng)量、瓦斯?jié)舛惹闆r通常由通風(fēng)部門的瓦斯檢測員每班進(jìn)行測定并填寫報表，但限于人力、物力因素以及時間、空間的限制，在同一時刻獲取礦井各巷道內(nèi)通風(fēng)風(fēng)量、瓦斯?jié)舛确植记闆r是不可能完成的。通過數(shù)值模擬手段獲取通風(fēng)系統(tǒng)中所有巷道的通風(fēng)風(fēng)量、瓦斯?jié)舛葦?shù)據(jù)，為礦井通風(fēng)系統(tǒng)優(yōu)化改造提供理論參數(shù)，為礦井制定相關(guān)安全技術(shù)措施提供理論指導(dǎo)，對保障礦井的安全生產(chǎn)具有重要的意義。為此，張少云[1]利用遠(yuǎn)程監(jiān)控系統(tǒng)采集實時數(shù)據(jù)，實現(xiàn)了通風(fēng)網(wǎng)絡(luò)的實時解算。魏引尚[2]提出了通風(fēng)網(wǎng)絡(luò)中瓦斯分布的計算模型，并使用Visual Basic 6.0編寫了專門的計算機程序，并給出了算例。鄧敢博[3]利用概率統(tǒng)計方法，以通風(fēng)網(wǎng)絡(luò)解算為基礎(chǔ)對通風(fēng)巷道瓦斯流動及分布規(guī)律進(jìn)行了研究，并使用Visual C++6.0編寫了相應(yīng)程序。以上軟件程序在解算過程中均需要大量處理矩陣運算，導(dǎo)致了其編制周期長，編寫過程過于復(fù)雜。

由于MATLAB軟件在通風(fēng)網(wǎng)絡(luò)矩陣運算處理、圖形顯示方面具有較為直觀的形式[4-7]，本文使用MATLAB軟件編寫通風(fēng)、瓦斯解算程序，并以某礦301綜放工作面通風(fēng)系統(tǒng)為算例，構(gòu)建通風(fēng)瓦斯網(wǎng)絡(luò)解算數(shù)學(xué)、物理模型，采用回路Cross迭代算法，通過圖形直觀顯示解算結(jié)果，為工作面安全生產(chǎn)提供可靠參數(shù)。

1 301綜放工作面概況

某礦301綜放工作面位于403采區(qū)南部，東為3條上山保護煤柱，西為井田邊界，南為井田邊界，北為403采區(qū)待采實體煤。工作面標(biāo)高在675～750 m，地面標(biāo)高在1 030～1 050 m，埋深約350 m。301綜放工作面為403采區(qū)首采工作面，其走向長2 845 m，傾斜長180 m，綜放開采，開采的4#煤層平均厚16.7 m，平均開采厚度為13.3 m，煤層傾角平均為22°。工作面通風(fēng)采用礦井負(fù)壓U型通風(fēng)方式，布置運輸順槽、回風(fēng)順槽以擔(dān)負(fù)通風(fēng)任務(wù)。距離工作面回風(fēng)順槽15 m(中對中)布置專用瓦斯抽放巷道以抽放采空區(qū)瓦斯。301綜放工作面通風(fēng)系統(tǒng)示意見圖1。

2 301工作面瓦斯情況分析

2.1 工作面瓦斯涌出量

工作面瓦斯涌出分為煤壁瓦斯涌出、采空區(qū)瓦斯涌出、落煤瓦斯涌出三大部分[8]。目前工作面絕對瓦斯涌出量約59.02 m3/min,其中地面低負(fù)壓抽采系統(tǒng)抽采混合量約500 m3/min,瓦斯?jié)舛燃s9.0%，純量為45 m3/min；井下高負(fù)壓抽采系統(tǒng)抽采混合量約117 m3/min，瓦斯?jié)舛燃s8.7%，純量為10.18 m3/min；回風(fēng)巷口風(fēng)量為1 745 m3/min，瓦斯?jié)舛燃s0.22%，風(fēng)排瓦斯為3.84 m3/min。301綜放工作面瓦斯涌出來源見圖2。

2.2 各巷道瓦斯涌出量

為統(tǒng)計301綜放工作面各巷道瓦斯涌出情況，同時測試巷道起、末節(jié)點瓦斯?jié)舛炔钪?，并且記錄巷道的風(fēng)量數(shù)據(jù)，即可以計算此條巷道的瓦斯涌出量。301工作面各條巷道瓦斯涌出情況見表1。

圖1 301工作面通風(fēng)系統(tǒng)示意

圖2 301綜放工作面瓦斯源示意

表1 301綜放工作面各巷道瓦斯涌出量

3 通風(fēng)瓦斯流動耦合仿真模擬數(shù)學(xué)模型

3.1 風(fēng)網(wǎng)拓?fù)潢P(guān)系

3.1.1 巷道坐標(biāo)提取

為節(jié)省時間，提高效率，使用CAD軟件作為繪圖系統(tǒng)。在CAD軟件系統(tǒng)中繪制“通風(fēng)網(wǎng)絡(luò)”圖層，在繪制通風(fēng)巷道中，使用單線表示巷道，并確保每一條巷道的通風(fēng)系統(tǒng)沒有大的變化。巷道連接處、通風(fēng)系統(tǒng)變化處使用不同的直線線條繪制。CAD軟件中，使用“工具”中“數(shù)據(jù)提取”，即可以提取圖層中巷道三維坐標(biāo)，并保存為Excel文件。

3.1.2 拓?fù)潢P(guān)系自動建立

拓?fù)潢P(guān)系的自動建立是為通風(fēng)、瓦斯網(wǎng)絡(luò)解算所需要的迭代解算矩陣提供必要條件的。使用MATLAB軟件讀取生成的“通風(fēng)網(wǎng)絡(luò)”坐標(biāo)Excel文件，當(dāng)讀取第一條巷道的三維坐標(biāo)(X1，Y1，Z1)和(X2，Y2，Z2)后，設(shè)置坐標(biāo)棧，并讓第一點坐標(biāo)入棧，并自動標(biāo)記為V1和V2節(jié)點。接下來程序會讀取第二點坐標(biāo)(X3，Y3，Z3)和(X4，Y4，Z4)，此時坐標(biāo)入棧策略是(X3，Y3，Z3)坐標(biāo)首先與棧中所有坐標(biāo)對比，如果有，則入棧并標(biāo)記為相應(yīng)節(jié)點，如果沒有，則入棧并標(biāo)記為棧中已經(jīng)存在的所有節(jié)點標(biāo)號最大值加1，即Vmax+1，其中Vmax表示棧中已經(jīng)存在的最大的節(jié)點標(biāo)號，(X4，Y4，Z4)節(jié)點入棧策略與此相同。當(dāng)所有巷道坐標(biāo)都已經(jīng)入棧并標(biāo)記節(jié)點編號，則循環(huán)結(jié)束，此時風(fēng)網(wǎng)拓?fù)潢P(guān)系自動建立成功。

3.2 通風(fēng)瓦斯網(wǎng)絡(luò)解算

通風(fēng)網(wǎng)絡(luò)解算中風(fēng)量與阻力需滿足節(jié)點風(fēng)量平衡定律、回路風(fēng)壓平衡定律及阻力定律三大定律。風(fēng)網(wǎng)解算的基本方程如下：

(1)

式中,bij為基本關(guān)聯(lián)矩陣元素；qj為分支流量；fi(q1,q2,…,qn-m+1)為回路阻力平衡方程；cij為基本回路矩陣元素；rj為分支風(fēng)阻,(N·s2)/m8；hj為回路附加阻力，Pa。

通過基本關(guān)聯(lián)矩陣可構(gòu)造m-1個線性無關(guān)的方程組，通過流量平衡方程與n個分支的阻力平衡方程聯(lián)系，可得到(m-1)+(n-m+1)個，即n個方程，解此方程，一般采用回路的Cross算法，其中基本回路迭代風(fēng)量修正值為[9]

(2)

式中,bij為基本關(guān)聯(lián)矩陣元素；Rj為與i鄰接的j分支的風(fēng)阻，(N·s2)/m8;Qj為j分支的風(fēng)量，m3/s；Pi為i回路的自然風(fēng)壓，Pa；Fi1和Fi2為風(fēng)機特性曲線方程表示的風(fēng)壓及風(fēng)機特性曲線的斜率。

對基本回路風(fēng)量進(jìn)行修正，直至風(fēng)量修正值滿足精度要求，即可求出風(fēng)網(wǎng)中所有巷道的風(fēng)量。

在解算通風(fēng)巷道內(nèi)所有風(fēng)量后，對巷道內(nèi)瓦斯流量進(jìn)行解算。通風(fēng)巷道內(nèi)的瓦斯流動解算與風(fēng)網(wǎng)解算原理基本相同，主要是滿足節(jié)點流量守恒定律，即文獻(xiàn)[2]中風(fēng)網(wǎng)內(nèi)任意節(jié)點相關(guān)分支的瓦斯流量代數(shù)和為零，即

(3)

式中,Qij為與i節(jié)點相關(guān)聯(lián)的分支j的有向風(fēng)量；n為節(jié)點相關(guān)聯(lián)的分支數(shù)。

在測定通風(fēng)巷道內(nèi)的瓦斯后，利用求解風(fēng)量即可求出網(wǎng)絡(luò)內(nèi)巷道瓦斯分布。

4 綜放工作面通風(fēng)瓦斯流動耦合仿真模擬

4.1 工作面通風(fēng)阻力系數(shù)測試

構(gòu)建301綜放工作面通風(fēng)瓦斯解算系統(tǒng)，首先要對工作面進(jìn)行通風(fēng)阻力系數(shù)測試，測試內(nèi)容主要包括名稱、支護、形狀，斷面積、周長、長度、風(fēng)量、巷道瓦斯涌出量等[10]。測試結(jié)果見表2。

表2 301綜放工作面通風(fēng)系統(tǒng)阻力測試結(jié)果

4.2 通風(fēng)瓦斯解算數(shù)據(jù)處理

預(yù)先根據(jù)拓?fù)潢P(guān)系處理并生成的矩陣如下

所示，使用回路Cross算法進(jìn)行迭代計算。

關(guān)聯(lián)矩陣：

參數(shù)矩陣：

樹：tree=e1，e8，e2，e3，e5。

余樹：cotree=e4，e6，e7。

基本回路矩陣：

4.3 解算結(jié)果與應(yīng)用效果分析

對301綜放工作面通風(fēng)瓦斯流動進(jìn)行解算，使用MATLAB軟件的figure函數(shù)[11-12]繪制通風(fēng)瓦斯實時解算結(jié)果(圖3)，詳細(xì)數(shù)據(jù)見表3。

圖3 301綜放工作面通風(fēng)瓦斯解算示意

由圖3、表3可知，通風(fēng)風(fēng)量模擬結(jié)果與現(xiàn)場實測基本吻合；瓦斯流量模擬結(jié)果除了專用瓦斯抽放巷(e4)的解算濃度為0.59%，與實測濃度8.7%誤差較大外，其余解算濃度基本吻合現(xiàn)場實測數(shù)據(jù)。專用瓦斯抽放巷濃度較大是因為采空區(qū)抽采流量較大。

5 結(jié) 論

(1)對某礦301綜放工作面通風(fēng)系統(tǒng)參數(shù)進(jìn)行測試，得到工作面通風(fēng)風(fēng)量、瓦斯涌出量等數(shù)據(jù)，計算了工作面各條巷道瓦斯涌出量。

(2)使用MATLAB軟件建立了301綜放工作面通風(fēng)瓦斯網(wǎng)絡(luò)解算拓?fù)潢P(guān)系，求解了Cross迭代計算所需要的關(guān)聯(lián)矩陣、回路矩陣和瓦斯網(wǎng)絡(luò)分布解算的參數(shù)矩陣等，得出通風(fēng)風(fēng)量、瓦斯流量的模擬解算結(jié)果，與實測數(shù)據(jù)基本吻合，驗證了其準(zhǔn)確性和可靠性。

表3 301綜放工作面通風(fēng)瓦斯解算結(jié)果

(3)將MATLAB軟件應(yīng)用到通風(fēng)瓦斯網(wǎng)絡(luò)解算中，具有編程效率高、解算結(jié)果圖形顯示方便等優(yōu)點，另外，可以將此程序編譯成exe可執(zhí)行文件，方便運行。

[1] 張少云.礦井風(fēng)網(wǎng)實時解算與安全預(yù)警系統(tǒng)的研究和實現(xiàn)[D].西安:西北工業(yè)大學(xué),2005.

[2] 魏引尚,常心坦.瓦斯在通風(fēng)巷道中流動分布情況研究[J].西安科技大學(xué)學(xué)報,2005(3):271-273.

[3] 鄧敢博.基于網(wǎng)絡(luò)計算礦井瓦斯分布規(guī)律研究[D].西安:西安科技大學(xué),2009.

[4] 李宗翔,王德民,溫永宇.礦井3D風(fēng)網(wǎng)圖及基于MATLAB仿真編程實現(xiàn)[J].安全與環(huán)境學(xué)報,2010,10(6):168-171.

[5] 袁 梅.MATLAB在礦井通風(fēng)網(wǎng)絡(luò)解算中的應(yīng)用[J].礦業(yè)工程,2009,7(5):63-65.

[6] 姜詩明,裴紹宇,郄雷敏.基于MATLAB的礦井通風(fēng)網(wǎng)絡(luò)解算程序編制[J].陜西煤炭,2010(6):28-29.

[7] 陳梅芳,丁德馨,張傳飛.基于MATLAB的礦井通風(fēng)網(wǎng)絡(luò)圖的矩陣表示及電算方法[J].礦業(yè)研究與開發(fā),2007,27(6):65-67.

[8] 煤炭科學(xué)研究總院撫順分院.AQ 1018—2006礦井瓦斯涌出量預(yù)測方法[S].北京：中國標(biāo)準(zhǔn)出版社，2006.

[9] 劉 劍,賈進(jìn)章,鄭 丹.流體網(wǎng)絡(luò)理論[M].北京:煤炭工業(yè)出版社,2002.

[10] 煤炭科學(xué)研究總院撫順分院，遼寧工程技術(shù)大學(xué).MT/T 440—2008礦井通風(fēng)阻力測定方法[S].北京：中國標(biāo)準(zhǔn)出版社，2009.

[11] 陳 澤,占海明.詳解MATLAB在科學(xué)計算中的應(yīng)用[M].北京:電子工業(yè)出版社,2011.

[12] 王海英,黃 強,李傳濤,等.圖論算法及其MATLAB實現(xiàn)[M].北京:北京航空航天大學(xué)出版,2010.

2016-08-06)

呂 閆(1983—)，男，副科長，助理工程師，033399 山西省呂梁市柳林縣穆村鎮(zhèn)沙曲村。