采空區(qū)多參數(shù)流-固耦合溫度場模擬研究＊

張九零 陳慶亞 王月紅

（1.河北聯(lián)合大學礦業(yè)工程學院，河北省唐山市，063009；2.礦山熱動力災害與防治教育部重點實驗室，遼寧省阜新市，123000）

煤的自燃是環(huán)境散熱和煤體放熱相互作用的物理化學過程，煤自燃過程會產(chǎn)生大量的有毒有害氣體，嚴重威脅著工作面附近工作人員的安全。煤在自燃的過程中，煤體內(nèi)的各種氣體濃度和含量、溫度都在不斷變化，且溫度的分布不均勻，進而造成耗氧速率和放熱強度的變化，進一步造成采空區(qū)高溫區(qū)域的溫度、位置隨時間不斷變化。近年來，隨著科學技術(shù)的進步和煤礦管理水平的提高，煤礦自燃火災百萬噸發(fā)火率雖然每年有所下降，但是我國煤礦自燃火災仍較為嚴重，對煤礦的安全生產(chǎn)構(gòu)成嚴重威脅，所以研究影響采空區(qū)自然發(fā)火的因素是很有意義的。

1 基于有限體積法的采空區(qū)內(nèi)動態(tài)溫度場數(shù)學模型

根據(jù)采煤工作面隨著開采推進不斷移動的特點，依據(jù)有限體積法，利用傳熱學理論、能量守恒定律及多孔介質(zhì)理論對采空區(qū)的熱力分布進行分析，同時，由于采空區(qū)內(nèi)遺煤和氣體溫度之間存在差異，故將采空區(qū)內(nèi)的遺煤溫度和氣體溫度進行耦合關(guān)聯(lián)研究，建立移動坐標條件下的采空區(qū)溫度場多參數(shù)流-固耦合數(shù)學模型，如式 （1）所示。移動坐標的引進將采空區(qū)溫度場復雜的動態(tài)變化問題轉(zhuǎn)變?yōu)殪o態(tài)問題進行處理，避免了非穩(wěn)態(tài)方程組復雜、計算量較大的問題。

式中：vx——氣流沿x軸方向的分量，m／s；

vy——氣流沿y軸方向的分量，m／s；

λy——采空區(qū)內(nèi)遺煤導熱系數(shù)，W／m·℃；

λg——采空區(qū)內(nèi)氣體導熱系數(shù)，W／m·℃；

tg——采空區(qū)內(nèi)氣體溫度，K；

ts——采空區(qū)內(nèi)遺煤溫度，K；

ρg——采空區(qū)內(nèi)氣體密度，kg／m3；

ρs——采空區(qū)內(nèi)遺煤密度，kg／m3；

cs——采空區(qū)內(nèi)遺煤比熱容，kJ／ （kg·K）；

cg——采空區(qū)內(nèi)氣體比熱容，kJ／ （kg·K）；

Ke——對流換熱系數(shù)，J／ （m2·s·K）；

Sn——單元體內(nèi)遺煤與氣體對流換熱表面積，

m2；

q （t）——單位時間單位體積內(nèi)采空區(qū)遺煤的放熱量，kJ／ （m3·s）；

a——單元體的比表面積，1／m；

Δx——單元體的長，m；

Δy——單元體的寬，m；

Δz——單元體的高，m；

v0——工作面推進速度，m／s；

F——單元體的面積，m2；

?！獑卧w的邊界。

2 采空區(qū)動態(tài)多參數(shù)溫度場流-固耦合模擬

2.1 自燃模型數(shù)值計算范圍

為了確定邊界條件，需測定采空區(qū)內(nèi)上下兩巷位置的溫度，利用傳熱范圍邊界來替代采空區(qū)實際邊界，假定一個采空區(qū)的傳熱范圍，則這個假定的傳熱范圍邊界屬于絕熱邊界來處理，因為采空區(qū)的熱量傳遞不只限于在采空區(qū)的實際邊界內(nèi)進行，也與采空區(qū)相鄰的煤層頂?shù)装濉⒉煽諈^(qū)四周的煤壁等進行熱量交換，所以在采空區(qū)的邊界溫度或熱流通量不能明確的情況下，采空區(qū)溫度場的計算需要將采空區(qū)的事實邊界進行擴大，將采空區(qū)溫度的邊界外推到幾乎沒有熱流通量的地方，這樣就可以定義溫度場解算的邊界條件 （熱流通量為0）。根據(jù)經(jīng)驗將采空區(qū)自然發(fā)火模型的溫度場邊界條件分別向上和向下各推20m，以達到在此邊界屬于絕熱邊界，如圖1中虛線所確定的范圍為采空區(qū)溫度場模型的解算范圍。

圖1 采空區(qū)坐標系及數(shù)值求解范圍

2.2 網(wǎng)格劃分

對采空區(qū)進行三角形網(wǎng)格剖分，為在實際情況中減少計算量，在保證程序計算精度前提下，對采空區(qū)采用疏密不一的網(wǎng)格劃分方式，對采空區(qū)容易發(fā)生自燃的 “兩道兩線” （進風道、回風道、切眼線和停采線）處，因變化率較大 （本文未考慮停采的情況），進行網(wǎng)格加密處理，變化率較小的地方網(wǎng)格劃分較疏散，在不影響計算精度的情況下減小了計算量。采空區(qū)劃分后，共產(chǎn)生1000 個節(jié)點，劃分為4851個小單元。

2.3 參數(shù)選取

工作面長度100m，工作面傾角0°，工作面通風阻力為15Pa，工作面掘進速度為2.4m／d，采空區(qū)遺煤平均厚度為0.7m。采空區(qū)遺煤原始溫度30℃，工作面的進風溫度25℃，煤的密度1400kg／m3，臨界溫度下煤的平均比熱1200J／（kg·℃）。

2.4 主要采掘參數(shù)對采空區(qū)內(nèi)溫度場的影響

影響采空區(qū)溫度的因素有很多，除了本身煤體的性質(zhì)外，還和采煤方法及采煤工藝等因素有關(guān)。在采空區(qū)自然發(fā)火解算模型中提到了工作面長度、走向長度、推進速度、遺煤厚度等影響因素，由于篇幅有限，只對影響采空區(qū)溫度場變化的工作面長度、推進速度、通風阻力3個采掘參數(shù)進行分析。

2.4.1 工作面長度對采空區(qū)溫度的影響

為了研究采空區(qū)工作面長度對采空區(qū)溫度場的影響，其他參數(shù)不改變，對工作面長度50 m、80 m、120m、150m、200m 下采空區(qū)的溫度場進行解算，并用軟件Tecplot10.0 對解算的數(shù)據(jù)處理，以圖像的形式對采空區(qū)內(nèi)氣體溫度進行顯示，如圖2所示。

圖2 不同工作面長度時采空區(qū)溫度場分布示意圖

由圖2 （a）可知，采空區(qū)氣體溫度沿著x、y軸的方向逐漸升高，且沿著x 軸方向溫度變化劇烈。沿著x 軸方向在距離工作面100 m 處溫度由開始設(shè)定的25℃快速升高到41℃，在距離采空區(qū)130 m 處溫度逐漸降低到39℃，之后變化平穩(wěn)。沿著y 軸方向溫度逐漸升高，在距離回風口28m處溫度達到最大值41℃。圖2 （b）可知，采空區(qū)的溫度規(guī)律和圖2 （a）相似，都是沿著x、y 軸的方向逐漸升高，沿著x 軸方向溫度變化劇烈，工作面長度為200m 時采空區(qū)的溫度達到了45℃。

圖3 不同工作面長度下采空區(qū)內(nèi)溫度分布示意圖

圖3為不同工作面長度下采空區(qū)溫度分布示意圖，可以看出，隨著工作面長度的增加，溫度在不斷升高。由圖3 （a）可以看出，在工作面長度為50m 時采空區(qū)出現(xiàn)的最高溫度為41℃，隨著距離工作面由20 m 到200 m 時，在距離工作面20 m的條件下工作面的長度為0 m 時，采空區(qū)的溫度為25℃，到距離工作面200 m 的條件下工作面的長度為0m 時，采空區(qū)的溫度升高為37℃。由圖3（b）可以看出，在工作面長度為200 m 時采空區(qū)出現(xiàn)最高溫度46℃，隨著距離工作面由20 m 到200m 時，采空區(qū)溫度由30℃升高到46℃，出現(xiàn)最高溫度的位置向采空區(qū)深部移動。采空區(qū)最高溫度隨著工作面長度的增加而升高，雖然采空區(qū)出現(xiàn)最高溫度的地點距離工作面由100 m 延后到了150m，但是采空區(qū)煤自燃發(fā)生的可能性還是很大，所以選擇適當?shù)墓ぷ髅娴拈L度可以在一定程度上減小煤自燃的幾率。

工作面長度為50m 時，在采空區(qū)淺部130m處對溫度影響最大，溫度也達到最高，之后溫度下降，在采空區(qū)深部超過150m 時，溫度變化緩慢，溫度的變化不再受工作面長度的變化。工作面長度為200 m 時，在采空區(qū)淺部150 m 處影響最大，在采空區(qū)深部超過150 m 時，溫度變化緩慢，溫度的變化不再受工作面長度的變化。由此可以得出，工作面長度越長，采空區(qū)內(nèi)遺煤氧化的熱量就越不容易排出，為采空區(qū)遺煤升溫氧化自燃的發(fā)生提供了充足的時間。由此得出，工作面的長度越長，采空區(qū)遺煤的自燃速度就越快。

2.4.2 推進速度與采空區(qū)溫度場相互關(guān)系

為了研究推進速度對采空區(qū)溫度場的影響，在其他參數(shù)不變的情況下，增加采空區(qū)采煤的推進速度，速度由開始的1m／d分別增加到2m／d、3m／d、4m／d、5m／d。對采空區(qū)溫度場進行解算，并用軟件Tecplot 10.0對解算數(shù)據(jù)進行處理，對采空區(qū)內(nèi)氣體溫度以圖像的形式顯示，如圖4所示。

圖4 不同推進速度下采空區(qū)溫度場分布示意圖

由圖4 （a）可知，采空區(qū)氣體溫度沿著x、y軸的方向逐漸升高，且沿著x 軸方向溫度變化劇烈。沿著x 軸方向在距離工作面104 m 處溫度由開始設(shè)定的25℃快速升高到66℃，在距離采空區(qū)140 m 處溫度逐漸降低到50℃，之后變化平穩(wěn)。沿著y 軸方向溫度逐漸升高，在距離回風口60m處溫度達到最大值66℃。圖4 （b）中，采空區(qū)氣體溫度沿著x、y 軸的方向逐漸升高，沿著x 軸方向在距離工作面150 m 處溫度由開始設(shè)定的26℃升高到37℃，之后變化平穩(wěn)。

圖5 不同推進速度下采空區(qū)溫度分布

圖5 為不同推進速度下采空區(qū)溫度分布示意圖，在圖5 （a）中，隨著距離工作面越遠，采空區(qū)的溫度逐漸升高，溫度由開始的28℃升高到66℃。在圖5 （b）中，推進速度由1m／d增加到5 m／d時，采空區(qū)內(nèi)氣體出現(xiàn)的最高溫度明顯的降低，最高溫度由66.7℃降低到37.8℃，出現(xiàn)最高溫度的地點距離工作面的距離分別為104.5 m 和153m，速度越快則出現(xiàn)最高溫度的地點離工作面越來越遠。

推進速度為1 m／d 時，在采空區(qū)淺部104 m處對采空區(qū)內(nèi)溫度影響最大，溫度也達到最高，之后溫度下降，在采空區(qū)深部超過120 m 時，溫度變化緩慢，溫度的變化不再受工作面長度變化的影響。推進速度為5 m／d 時，在采空區(qū)淺部150 m處影響最大，在采空區(qū)深部超過160 m 時，溫度變化緩慢，溫度的變化不再受工作面長度變化的影響。

隨著工作面推進速度的加快，采空區(qū)內(nèi)最高溫度急劇降低，推進速度對采空區(qū)最高溫度影響較大，這是由于隨著工作面移動，流進采空區(qū)低溫氣流和采空區(qū)遺煤增多，促進與采空區(qū)之間的熱傳導，對采空區(qū)起到冷卻的作用，致使采空區(qū)溫度降低，有效地降低了采空區(qū)煤自然發(fā)火的發(fā)生幾率。

2.4.3 通風阻力對采空區(qū)溫度場的影響

為了研究通風阻力對采空區(qū)溫度場的影響，在其他參數(shù)不變的情況下，增加采空區(qū)的通風阻力，由開始的10Pa分別增加到20Pa、30Pa、40Pa、50Pa。利用采空區(qū)自然發(fā)火模型對采空區(qū)的溫度場進行解算，用軟件Tecplot10.0對解算數(shù)據(jù)進行處理，并對采空區(qū)內(nèi)氣體溫度以圖像的形式進行顯示，如圖6所示。

圖6 不同通風阻力時采空區(qū)溫度場分布示意圖

由圖6 （a）可見，隨著通風阻力的增大，采空區(qū)的最高溫度有所升高，采空區(qū)氣體溫度沿著x、y 軸的方向逐漸升高，且沿著x 軸方向溫度變化劇烈。沿著x軸方向在距離工作面122m 處溫度由開始設(shè)定的25℃快速升高到45℃，沿著y軸方向溫度逐漸升高，在距離回風口48m 處溫度達到最大值45℃。圖6 （b）中采空區(qū)的溫度沿著x、y軸的方向逐漸升高，在距離工作面160 m 處溫度由開始設(shè)定的28℃快速升高到49℃，之后變化平穩(wěn)。

圖7 為不同通風阻力下采空區(qū)溫度分布示意圖，在圖7 （a）顯示隨著距離工作面越遠，采空區(qū)的溫度逐漸升高，溫度由開始的28℃升高到45℃，圖7 （b）顯示溫度由開始的28℃升高到49℃。通風阻力由10Pa增加到50Pa時，采空區(qū)內(nèi)氣體出現(xiàn)的最高溫度有所升高，由45℃升高到49℃，出現(xiàn)最高溫度的地點距離工作面的距離由122m 延后到了158m。

圖7 不同通風阻力下采空區(qū)溫度分布示意圖

通風阻力為10Pa時，在采空區(qū)淺部120m 處對溫度影響最大，溫度也達到最高，之后溫度下降，在采空區(qū)深部超過130m 時，溫度變化緩慢，溫度的變化不再受工作面長度變化的影響。通風阻力為50Pa時，在采空區(qū)淺部160m 處影響最大，在采空區(qū)深部超過170 m 時，溫度變化緩慢，溫度的變化不再受工作面長度變化的影響。

雖然增大通風阻力后采空區(qū)內(nèi)的溫度增加的幅度較小，但通風阻力較大時，采空區(qū)的漏風量較大，增加了煤的氧化放熱幾率，使得采空區(qū)內(nèi)的溫度逐漸升高，且向離工作面遠的方向移動。盡管通風阻力會使采空區(qū)出現(xiàn)最高溫度的地點向遠離工作面的方向移動，但是溫度會逐漸升高，所以為了避免煤自燃的發(fā)生，可以選擇合適的通風阻力來保障安全生產(chǎn)。

3 結(jié)論

（1）通過有限體積法 （FVM）來建立采空區(qū)內(nèi)動態(tài)溫度場多參數(shù)耦合數(shù)學模型，將動態(tài)問題轉(zhuǎn)化為靜態(tài)問題，避免了計算量過大的問題，并且確定了采空區(qū)自然發(fā)火模型計算范圍，將采空區(qū)自然發(fā)火模型的溫度場邊界條件分別向上和向下各推20m來達到絕熱邊界。對采空區(qū)進行三角形劃分后選取參數(shù)，利用Tecplot10.0軟件對解算數(shù)據(jù)進行處理，并以圖像顯示采空區(qū)溫度場的變化情況，使得結(jié)果顯示簡單明了。

（2）當工作面長度增加時，采空區(qū)內(nèi)的最高溫度升高，雖然采空區(qū)出現(xiàn)最高溫度的地點逐漸遠離工作面，但是出現(xiàn)自燃的可能性還是很大，所以選擇適當?shù)墓ぷ髅骈L度可以有效降低采空區(qū)內(nèi)溫度，防止自燃事故的發(fā)生。

（3）當采空區(qū)采煤的推進速度增加時，采空區(qū)出現(xiàn)的最高溫度逐漸變低，由1 m／d到5 m／d溫度變化較為顯著，當推進速度超過5 m／d時，溫度變化較小，而且出現(xiàn)最高溫度的地點一直遠離工作面，所以要保證生產(chǎn)，降低采空區(qū)煤自燃的發(fā)生必須選擇合適的推進速度。

（4）隨著工作面通風阻力的增加，采空區(qū)出現(xiàn)的最高溫度有所升高，并且采空區(qū)出現(xiàn)高溫的地點向遠離工作面的方向移動，所以選擇適當?shù)耐L阻力也可以有效保證煤礦的安全生產(chǎn)。

［1］ 紀玉龍，時國慶等.綜放面采空區(qū)瓦斯抽放條件下氧濃度場的CFD 模擬 ［J］.中國煤炭，2012 （4）

［2］ 徐精彩，張辛亥等.煤層自燃膠體防滅火理論與技術(shù) ［M］.北京：煤炭工業(yè)出版社，2003

［3］ 羅海珠，梁運濤.煤自然發(fā)火預測預報技術(shù)的現(xiàn)狀與展望 ［J］.中國安全科學學報，2003 （3）

［4］ 鄧軍，徐精彩等.國內(nèi)外煤炭自然發(fā)火預測預報技術(shù)綜述 ［J］.西安礦業(yè)學院學報，1999 （4）

［5］ 李宗翔，韋涌清等.非均質(zhì)采空區(qū)氣-固耦合溫度場風有限元求解 ［J］ .昆明理工大學學報，2004（2）

［6］ 鄔劍明，遲克勇等.煤自然發(fā)火過程的模擬實驗研究 ［J］.中國煤炭，2008 （4）

［7］ 張國樞.通風安全學 ［M］.徐州：中國礦業(yè)大學出版社，2000

［8］ 張國樞，戴廣龍.煤炭自燃理論與防治實踐 ［M］.北京：煤炭工業(yè)出版社，2002

［9］ 牛立東，王海平 .水地莊煤礦井下自燃案例探討［J］.中國煤炭，2011 （7）

［10］ 王同科.一維二階橢圓和拋物線微分方程的高精度有限體積元方法 ［J］ .數(shù)值計算與計算機應用，2002 （4）

［11］ 尹哲 .多孔介質(zhì)滲流問題的對稱有限體積方法［D］.山東：山東大學，2007