王月紅，吳 怡，張九零，關雅潔

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

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

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

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

1 理論基礎

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

1.1 有限體積法

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

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

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

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

WC=W1+W2+W3

(1)

(2)

(3)

因為在采空區(qū)存在氣體的動力彌散現(xiàn)象，所以氧氣的擴散系數(shù)與風流速度大小成正比，故W3：

(4)

(5)

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

(6)

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

(7)

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

(8)

式中，n為有效孔隙度；ρO2為氧氣的密度，kg/m3；CO2為氧氣的濃度，mol/m3；kO2為氧氣的擴散系數(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ù)學模型，如(9)所示。

(9)

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

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

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

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

1.3.1 固體溫度場

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

Qsc+Qf=Qd+Qv

(10)

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

(11)

式中，k為滲透率，m2。

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

(12)

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

(13)

(14)

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

(15)

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

(16)

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

(17)

式中，λy為采空區(qū)冒落巖石的導熱系數(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為工作面推進速度，m/s。

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

1.3.2 氣體溫度場

設有溫度場Qg(Q,t)=Qg(x,y,t)，由能量守恒定律可知，單位時間內(nèi)，匯入封閉曲線范圍的氣體熱量Qgc，同控制體體內(nèi)固體對流散熱量Qd之和，應等于單元體內(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ū)氣體導熱系數(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)是移動坐標下采空區(qū)內(nèi)氣體溫度場能量方程。

2 軟件開發(fā)

2.1 軟件界面

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

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

2.2 軟件功能

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

3 數(shù)值模擬

3.1 網(wǎng)格劃分

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

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

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

然后在采空區(qū)中注入氮氣，選取一個注氮點，假設在該點注入氮氣，氮氣進入采空區(qū)后會在該點向四周擴散，劃分了如圖3所示的第二層網(wǎng)格。如圖所示坐標為(e，f)的點的氮氣含量系數(shù)為(e2+f2)1/2，用DI表示注入的氮氣量，則該點G的氮氣含量為DIG=DI×(e2+f2)1/2。同理，其余點的氮氣量也可據(jù)此求出，這樣就得出了采空區(qū)內(nèi)所有點的氮氣量。

圖3 注氮含量計算示意

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

圖4 網(wǎng)格劃分

3.2 參數(shù)設置

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

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

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

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

3.3 模擬結果

3.3.1 氧濃度場

對比未注氮(如圖5(a))和注氮(如圖6(a))的氧濃度圖可知，注氮使氧化帶變窄，窒息帶變寬，且氧化帶向工作面方向推進。因為隨著漏風氣流流入采空區(qū)中，氧氣持續(xù)向深部擴散，但氮氣的注入不斷稀釋氧濃度，使氧濃度下降。且溫度在注氮后逐漸降低，減輕了熱量積聚放熱量降低，使得耗氧速率減小，氧化帶寬度則變窄。

圖5 未注氮氧濃度模擬

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

圖6 注氮后氧濃度場模擬

3.3.2 溫度場

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

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

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

圖7(b)(注氮前)與8(b)(注氮后)三維圖對比發(fā)現(xiàn)，不僅注氮與未注氮最高點溫度不同，且在進風側附近，注氮后的溫度隨著采空區(qū)的深度的增加緩慢上升，而在回風側附近，上升速率較大。因為氮氣進入采空區(qū)后，由于氮氣的低溫降低了采空區(qū)內(nèi)遺煤溫度上升的速率，而在回風巷附近，由于逆風的作用，氮氣的溫度較低，所以大量的熱量積聚在上隅角的位置處，使得該處的溫度很高。

3.4 工程應用

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

4 結 論

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

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

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