采空區(qū)注氮軟件模擬及開發(fā)研究

王月紅，吳 怡，張九零，關(guān)雅潔

(1.華北理工大學(xué) 礦業(yè)工程學(xué)院，河北 唐山 063009；2.河北省礦山開發(fā)與安全技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北 唐山 063009)

由于采空區(qū)中存在的大量遺煤和適宜漏風(fēng)環(huán)境為煤氧化創(chuàng)造了條件，增大了自燃的危險(xiǎn)性[1]。而注氮作為防治遺煤自燃的有效措施之一，它具有稀釋采空區(qū)氧濃度、冷卻遺煤、增大采空區(qū)壓力、減少漏風(fēng)量、延長自然發(fā)火期等優(yōu)點(diǎn)，近年來被廣泛應(yīng)用[2-5]。

由于數(shù)值模擬成本低廉，只需在計(jì)算機(jī)上進(jìn)行模擬和數(shù)據(jù)處理，有時可達(dá)到理想化克服現(xiàn)場的困難。國內(nèi)外學(xué)者關(guān)于模擬采空區(qū)注氮技術(shù)的成果很多。在模擬漏風(fēng)流場和溫度場中，李宗祥[6-7]針對Y形通風(fēng)的工作面，綜合各種傳熱方程，并對模型參數(shù)及邊界設(shè)計(jì)后進(jìn)行注氮問題的數(shù)值模擬，得到了各種氣體(瓦斯、CO、氧氣等)、溫度和漏風(fēng)流的分布情況及變化規(guī)律，并且計(jì)算出了注氮與工作面推進(jìn)速度的配比關(guān)系。美國學(xué)者Yuan等[8]根據(jù)CFD模型研究采空區(qū)漏風(fēng)流場變化規(guī)律，并通過研究結(jié)果完善現(xiàn)場通風(fēng)技術(shù)，控制災(zāi)害的發(fā)生。羅新榮[9]采用Fluent軟件建立物理模型和數(shù)學(xué)模型，設(shè)置相應(yīng)的參數(shù)，并根據(jù)采空區(qū)滲透率和遺煤氧化范圍編制自定義函數(shù)，模擬不同抽采量和注氮形式對采空區(qū)溫度場的影響規(guī)律。在關(guān)于注氮前后氧濃度的變化規(guī)律上，劉玉良[10-11]將理論和數(shù)值模擬結(jié)合對淺埋深煤層注氮技術(shù)進(jìn)行研究，比較了注氮前后“氧化帶”的范圍。澳大利亞學(xué)者Ren[12]利用軟件模擬詳細(xì)研究了瓦斯的運(yùn)移規(guī)律和氧濃度分布情況，并結(jié)合理論知識控制自然發(fā)火問題。唐冠楚[13-15]在CFD模型下，模擬不同的抽采方式與注氮技術(shù)相結(jié)合自燃危險(xiǎn)區(qū)域的分布情況，在瓦斯抽采效率提高的同時防止采空區(qū)遺煤自燃的發(fā)生。

但是，國內(nèi)外均用流體動力學(xué)軟件進(jìn)行模擬，不能更針對性地研究采空區(qū)注氮防滅火技術(shù)。筆者基于有限體積法，利用VB語言編制注氮防滅火模擬系統(tǒng)，通過Tecplot軟件作為圖顯，進(jìn)行溫度場和氧濃度場的模擬。利用編制的軟件模擬采空區(qū)注氮，減少了現(xiàn)場的大量實(shí)驗(yàn)，應(yīng)用于各種類型下的采空區(qū)注氮問題。

1 理論基礎(chǔ)

在采煤工作持續(xù)進(jìn)行中，工作面向前推進(jìn)，相對區(qū)域也持續(xù)擴(kuò)大，因此選取一部分采空區(qū)作為研究對象，引入動態(tài)坐標(biāo)。在注入氮?dú)夂?，采空區(qū)氣體壓力、氧氣濃度、溫度均會產(chǎn)生影響?；谟邢摅w積法建立數(shù)學(xué)模型，利用數(shù)學(xué)語言對實(shí)際問題進(jìn)行定量化，使指標(biāo)體系可計(jì)算。

1.1 有限體積法

有限體積法(FVM)是將計(jì)算區(qū)域劃分為一系列不重復(fù)的控制體積，使每個網(wǎng)格點(diǎn)周圍有一個控制體積；將待解的微分方程對每一個控制體積積分，然后對積分式進(jìn)行離散化處理，再導(dǎo)出離散化方程。

有限體積法在流體流動和傳熱數(shù)值計(jì)算領(lǐng)域被廣泛應(yīng)用，與其他有限元和邊界方法比較，其基本思想簡潔，能夠在物理意義上表示控制體的通量平衡，建立的模型方程意義明確。因此，將有限體積法應(yīng)用到采空區(qū)注氮防滅火的模型構(gòu)建上，為解決礦井火災(zāi)易復(fù)燃、難撲滅等問題提供了有效支撐。

1.2 采空區(qū)注氮氧濃度場數(shù)學(xué)模型

根據(jù)質(zhì)量守恒定律：在單位時間內(nèi)，漏風(fēng)風(fēng)流通過任意封閉區(qū)域內(nèi)氧氣質(zhì)量的變化量(WC)主要由以下幾個方面引出：由于滲透風(fēng)流中的氧氣流出與流入任意封閉曲面之差(W1)；氧氣消耗量(W2)；濃度差異引起的彌散進(jìn)出任意封閉曲面的氧氣之差(W3)，如公式(1)，(2)，(3)所示。

WC=W1+W2+W3

(1)

(2)

(3)

因?yàn)樵诓煽諈^(qū)存在氣體的動力彌散現(xiàn)象，所以氧氣的擴(kuò)散系數(shù)與風(fēng)流速度大小成正比，故W3：

(4)

(5)

即將式(2)，(3)，(4)以及(5)代入(1)得：

(6)

在實(shí)際情況中，氧氣流動速率要大于采空區(qū)工作面的回采速度，即單位時間內(nèi)選取的微元ΔS內(nèi)的孔隙率和氧氣密度變化可以忽略不計(jì)，即式(6)可變成：

(7)

為了和Fick動律保持一致，將(7)式中質(zhì)量濃度改為摩爾濃度。

(8)

式中，n為有效孔隙度；ρO2為氧氣的密度，kg/m3；CO2為氧氣的濃度，mol/m3；kO2為氧氣的擴(kuò)散系數(shù)常數(shù)；vx為x方向氣體滲流速度，m/s；vy為y方向氣體滲流速度，m/s；v為氣體的滲流速度，m/s；u(t)為松散煤體內(nèi)的氧氣消耗速度，mol/(s·m3)。

確定邊界條件Γ1，Γ2，Γ3，Γ4后，得到氧濃度場的數(shù)學(xué)模型，如(9)所示。

(9)

式中，CO2(x，y)為回采工作面連接的邊界上的氧氣濃度函數(shù)，mol/m3。

根據(jù)模型可以分析出氧濃度在注氮整個過程中的變化，從而研究出注氮對采空區(qū)火災(zāi)的影響作用。

1.3 采空區(qū)注氮溫度場數(shù)學(xué)模型

在采空區(qū)中，除了巖石間導(dǎo)熱的過程外，還有氣體與氣體間互相導(dǎo)熱，它們導(dǎo)熱的方式是有區(qū)別的，除此之外，固體間的導(dǎo)熱還會受到氣體的影響，因此在采空區(qū)注氮模型中的溫度場分為固體溫度場和氣體溫度場。

1.3.1 固體溫度場

在一定的時間內(nèi)，導(dǎo)入和導(dǎo)出冒落巖石的熱量差Qsc、冒落巖石氧化放熱的熱量差Qf之和，也就是進(jìn)入控制體的熱量，應(yīng)該與導(dǎo)出控制體的熱量相等，也就是冒落巖石及孔隙中氣體的對流換熱量Qd和冒落巖石移動帶出控制體的熱量差Qv之和。即為公式(10)所示。

Qsc+Qf=Qd+Qv

(10)

由傅里葉定律得知，由邊界Γ流入曲面F全部熱量Qsc：

(11)

式中，k為滲透率，m2。

單位時間內(nèi)控制體的氧化放熱量為：

(12)

在松散煤體內(nèi)同時存在對流換熱和導(dǎo)熱兩種過程，這兩種熱量傳遞方式也相互影響，同時還受氣體運(yùn)動規(guī)律的影響，因此這是一個十分復(fù)雜的過程，利用牛頓冷卻公式進(jìn)行計(jì)算，如公式(13)所示。

(13)

(14)

將式(11)，(12)，(13)，(14)代入式(10)得：

(15)

根據(jù)傅里葉定律可知：

(16)

將式(16)代入(15)得：

(17)

式中，λy為采空區(qū)冒落巖石的導(dǎo)熱系數(shù)，W/(m·℃)；ts為冒落巖石溫度，K；qx為沿x軸方向的熱流通量；qy為沿y軸方向的熱流通量；Ke為對流換熱系數(shù)，J/(m2·s·K)；Sn為單元體內(nèi)冒落巖石與氣體對流換熱的表面積，m2；tg為氣體溫度，K；ts為冒落巖石溫度，K；q(t)為單位時間單位體積內(nèi)冒落巖石的放熱量，kJ/(m3·s)；ρs為固體的密度，kg/m3；cs為固體的比熱容kJ/(kg·K)；a為單元體的比表面積,a=Sn/ΔxΔyΔz，L/m；v0為工作面推進(jìn)速度，m/s。

即式(17)是移動坐標(biāo)下采空區(qū)內(nèi)冒落巖石能量方程。

1.3.2 氣體溫度場

設(shè)有溫度場Qg(Q,t)=Qg(x,y,t)，由能量守恒定律可知，單位時間內(nèi)，匯入封閉曲線范圍的氣體熱量Qgc，同控制體體內(nèi)固體對流散熱量Qd之和，應(yīng)等于單元體內(nèi)孔隙中氣體內(nèi)能的增加值Qh，即：

Qgc+Qd=Qh

(18)

由傅里葉定律得知，單位時間內(nèi)曲面邊界Γ流入曲面F全部熱量是Qgc：

(19)

熱量Qh可表示為：

(20)

對流換熱量Qd可表示為：

(21)

將式(19)，(20)，(21)代入(18)得：

(22)

根據(jù)傅里葉定律可知：

(23)

將(23)代入式(22)得：

(24)

式中，λg為采空區(qū)氣體導(dǎo)熱系數(shù)，W/(m·℃)；Ke為對流換熱系數(shù)，J/(m2·s·℃)；Sn為單元體內(nèi)冒落巖石與氣體對流換熱的表面積，m2；tg為氣體溫度，℃；ts為固體溫度，℃；a為單元體的比表面積，a=Sn/ΔxΔyΔz，L/m；vx，vy為氣流沿x軸、y軸方向的分量，m/s；ρg為氣體的密度，kg/m3；cg為氣體的比熱，J/(kg·℃)。

式(24)是移動坐標(biāo)下采空區(qū)內(nèi)氣體溫度場能量方程。

2 軟件開發(fā)

2.1 軟件界面

圖形顯示模塊的任務(wù)是把整個模擬過程處理存儲的數(shù)據(jù)用圖形顯示出來，具有一個功能齊全的圖形顯示界面是有限體積軟件中最基本的一個要求。對于大多數(shù)的用戶而言，圖形顯示功能是至關(guān)重要的，因?yàn)閳D形顯示能夠有效地減小大量數(shù)據(jù)結(jié)構(gòu)的堆積，進(jìn)而減小了結(jié)構(gòu)計(jì)算的錯誤。本設(shè)計(jì)軟件圖形顯示模塊分為3部分：登陸系統(tǒng)；參數(shù)輸入界面；結(jié)果匯總顯示。

顯示界面是Tecplot10.0作為繪圖以及數(shù)據(jù)分析情況下通用軟件，該軟件作為圖形顯示的優(yōu)點(diǎn)是易學(xué)易用，界面友好，能夠馬上讀入*.cas、*.dat文件等。因此基于Tecplot軟件的這些優(yōu)點(diǎn)，使其作為該軟件的圖形顯示模塊效果是非?？捎^的，而且該軟件較成熟，解決了不能準(zhǔn)確地處理數(shù)據(jù)的弊端。

2.2 軟件功能

該軟件的開發(fā)，不僅對采空區(qū)自然發(fā)火防治的研究起到了極大的作用，更是對注氮防滅火這一滅火技術(shù)有了更深入的研究。該軟件的開發(fā)可應(yīng)用到實(shí)際的生產(chǎn)中，通過各煤礦的實(shí)際情況，解算出注氮的最佳位置、最佳流量、最佳個數(shù)等實(shí)質(zhì)性問題；還可計(jì)算采空區(qū)注氮條件下“三場”(氧濃度場、流場、溫度場)的分布情況；避免了實(shí)際操作的困難，取而代之的是一次次的計(jì)算機(jī)上的操作，相對比人工操作精確而且可循環(huán)，試驗(yàn)性強(qiáng)，計(jì)算精度高，速度快，對煤礦解決火災(zāi)這個難題起到了極大的幫助。

3 數(shù)值模擬

3.1 網(wǎng)格劃分

首先通過計(jì)算一個控制體，達(dá)到計(jì)算出構(gòu)成采空區(qū)工作面所有單元的相關(guān)參數(shù)。為了清楚地表示出網(wǎng)格劃分的方法，用圖1所示MNPQ代表中心體積法中的矩形網(wǎng)格，四邊形abcd表示選取的控制體。沿ab方向?yàn)椴煽諈^(qū)深度，沿ad方向?yàn)楣ぷ髅娣较?，根?jù)采空區(qū)自然發(fā)火物理模型ab方向?yàn)椴煽諈^(qū)的進(jìn)風(fēng)巷，dc為回風(fēng)巷，在采空區(qū)中“兩道兩線”處為自燃易發(fā)部位，因此在對控制體abcd進(jìn)行網(wǎng)格剖分時在自燃易發(fā)部位處的網(wǎng)格加密。采空區(qū)其余工作面部位采取疏密不一的網(wǎng)格線進(jìn)行剖分，為了使得繪制的網(wǎng)格標(biāo)準(zhǔn)化和規(guī)范化，網(wǎng)格之間的距離引入數(shù)學(xué)理論中的等比數(shù)列且系數(shù)大于1，如圖2所示。

圖1 有限體積法矩形網(wǎng)格

圖2 控制體網(wǎng)格劃分

然后在采空區(qū)中注入氮?dú)?，選取一個注氮點(diǎn)，假設(shè)在該點(diǎn)注入氮?dú)?，氮?dú)膺M(jìn)入采空區(qū)后會在該點(diǎn)向四周擴(kuò)散，劃分了如圖3所示的第二層網(wǎng)格。如圖所示坐標(biāo)為(e，f)的點(diǎn)的氮?dú)夂肯禂?shù)為(e2+f2)1/2，用DI表示注入的氮?dú)饬?，則該點(diǎn)G的氮?dú)夂繛镈IG=DI×(e2+f2)1/2。同理，其余點(diǎn)的氮?dú)饬恳部蓳?jù)此求出，這樣就得出了采空區(qū)內(nèi)所有點(diǎn)的氮?dú)饬俊?/p>

圖3 注氮含量計(jì)算示意

通過將兩個網(wǎng)格疊加之后，從而進(jìn)行模型的數(shù)值計(jì)算。由Tecplot作為圖顯得到圖4的網(wǎng)格劃分形式。

圖4 網(wǎng)格劃分

3.2 參數(shù)設(shè)置

在林南倉煤礦井下采集具有代表性的煤樣進(jìn)行工業(yè)分析、元素分析等基礎(chǔ)實(shí)驗(yàn)，從而確定采空區(qū)注氮軟件開發(fā)基礎(chǔ)參數(shù)的設(shè)定，如表1所示。

表1 模型基礎(chǔ)參數(shù)

由采空區(qū)自然發(fā)火原理中的煤氧復(fù)合理論得出，空氣中的氧氣與煤的表面結(jié)合后會發(fā)生大量的化學(xué)反應(yīng)，并釋放大量的熱量，同時該過程更會加快煤與氧氣的反應(yīng)進(jìn)程，使得火災(zāi)或爆炸迅速發(fā)生。因此，針對這一理論知識，通過升溫氧化實(shí)驗(yàn)，確定煤氧化放熱參數(shù)，如表2所示。

表2 煤氧化放熱參數(shù)

3.3 模擬結(jié)果

3.3.1 氧濃度場

對比未注氮(如圖5(a))和注氮(如圖6(a))的氧濃度圖可知，注氮使氧化帶變窄，窒息帶變寬，且氧化帶向工作面方向推進(jìn)。因?yàn)殡S著漏風(fēng)氣流流入采空區(qū)中，氧氣持續(xù)向深部擴(kuò)散，但氮?dú)獾淖⑷氩粩嘞♂屟鯘舛龋寡鯘舛认陆?。且溫度在注氮后逐漸降低，減輕了熱量積聚放熱量降低，使得耗氧速率減小，氧化帶寬度則變窄。

圖5 未注氮氧濃度模擬

對比未注氮和注氮后的氧濃度場三維圖，發(fā)現(xiàn)兩者的氧濃度均是逐漸降低的，但氧濃度變化趨勢則由注氮前的凸形變?yōu)樽⒌蟮陌夹?。觀察圖5(b)和圖6(b)可知，三維圖中反曲線的凸率逐漸減小，說明隨著距回風(fēng)巷距離的增大，氧濃度改變的速率逐漸減小。證明了未注氮時在進(jìn)風(fēng)巷周圍由于漏風(fēng)流將氧氣帶入采空區(qū)深部，而到達(dá)回風(fēng)巷時，由于逆風(fēng)的作用，使得氧氣濃度的變化不再明顯；而注氮后，靠近進(jìn)風(fēng)巷的位置處氧濃度隨深度的改變呈凹形變化，且隨著與進(jìn)風(fēng)巷距離的加大，凹率逐漸減小。

圖6 注氮后氧濃度場模擬

3.3.2 溫度場

對比未注氮(如圖7(a))和注氮(如圖8(a))二維圖的溫度場可知，工作面附近的溫度較低，且注氮前后的溫度差只有2℃，說明注氮前后進(jìn)風(fēng)巷中與漏風(fēng)流的溫度相接近。但隨著風(fēng)流流向采空區(qū)深部，低溫氣體與采空區(qū)深部的冒落巖石及煤體進(jìn)行熱量交換，溫度急劇降低，因此采空區(qū)深部的溫度差較大。注氮前的最高溫度為70℃，注氮后的最高溫度僅為40℃，溫差達(dá)到了30℃左右。且注氮前后采空區(qū)溫度與深度呈正比例關(guān)系，采空區(qū)越深部溫度越高。注氮后進(jìn)風(fēng)巷附近的溫度有明顯的降低，而回風(fēng)巷周圍的溫度較總體有明顯的降低，但與注氮前相比較，高溫區(qū)域由注氮前的采空區(qū)深部位置變?yōu)榛仫L(fēng)巷部位。由于注入氮?dú)夂?，氮?dú)獗伙L(fēng)流帶入采空區(qū)深部，發(fā)生熱傳導(dǎo)降低原有遺煤溫度，熱量則積聚到回風(fēng)巷；又因回風(fēng)巷處的風(fēng)量較小，熱量不易被帶出采空區(qū)，導(dǎo)致回風(fēng)巷深部位置溫度最高。

圖7 采空區(qū)未注氮溫度場模擬

圖8 采空區(qū)注氮后溫度場模擬

圖7(b)(注氮前)與8(b)(注氮后)三維圖對比發(fā)現(xiàn)，不僅注氮與未注氮最高點(diǎn)溫度不同，且在進(jìn)風(fēng)側(cè)附近，注氮后的溫度隨著采空區(qū)的深度的增加緩慢上升，而在回風(fēng)側(cè)附近，上升速率較大。因?yàn)榈獨(dú)膺M(jìn)入采空區(qū)后，由于氮?dú)獾牡蜏亟档土瞬煽諈^(qū)內(nèi)遺煤溫度上升的速率，而在回風(fēng)巷附近，由于逆風(fēng)的作用，氮?dú)獾臏囟容^低，所以大量的熱量積聚在上隅角的位置處，使得該處的溫度很高。

3.4 工程應(yīng)用

根據(jù)林南倉煤礦的技術(shù)條件及生產(chǎn)能力確定相應(yīng)的注氮參數(shù)，礦井下通常運(yùn)用埋管方式進(jìn)行防滅火，沿回采工作面進(jìn)風(fēng)巷巷幫鋪設(shè)管路，伴隨著工作面推進(jìn)，管段注氮口會逐漸伸入深部，注氮口位置開始進(jìn)入氧化自燃帶時進(jìn)行注氮，每隔一段步距重新鋪設(shè)一條管路，在注氮的過程中采用封閉式注氮，即利用密閉墻上的注氮管向采空區(qū)中注入氮?dú)?。進(jìn)行現(xiàn)場注氮后發(fā)現(xiàn)溫度和氧氣濃度明顯降低，林南倉煤礦在正常生產(chǎn)的條件下是不易發(fā)生自燃的，因此注氮達(dá)到了相應(yīng)的效果，由此說明，該軟件的模擬結(jié)果切實(shí)可行，為日后采空區(qū)注氮技術(shù)研究起到了很大支撐作用。

4 結(jié) 論

(1)基于有限體積法，建立采空區(qū)注氮“氧濃度場”和“溫度場”數(shù)學(xué)模型。并利用VB語言編制注氮防滅火模擬系統(tǒng)，通過Tecplot軟件作為圖顯，開發(fā)了采空區(qū)注氮防滅火軟件?？蛇\(yùn)用到煤礦實(shí)際生產(chǎn)中，并有效解決煤礦火災(zāi)問題。

(2)根據(jù)模擬得到的氧濃度場，發(fā)現(xiàn)注氮前后氧化帶寬度變窄，并且向工作面方向推進(jìn)，窒息帶變寬。氧濃度變化曲線由凸形變?yōu)榘夹巍ｋS著回風(fēng)巷距離的增大，氧濃度變化速率越來越小。

(3)在溫度場中，注氮前后工作面附近的溫度差僅2℃，而采空區(qū)深部的溫度差為30℃，且采空區(qū)溫度與深度呈正比例關(guān)系。注氮后，進(jìn)風(fēng)巷附近溫度明顯降低，而回風(fēng)巷周圍溫度較總體明顯下降，但與注氮前比較高溫區(qū)域由深部位置變?yōu)榛仫L(fēng)巷位置。