采空區多參數流-固耦合溫度場模擬研究＊

張九零 陳慶亞 王月紅

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

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

1 基于有限體積法的采空區內動態溫度場數學模型

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

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

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

λy——采空區內遺煤導熱系數，W／m·℃；

λg——采空區內氣體導熱系數，W／m·℃；

tg——采空區內氣體溫度，K；

ts——采空區內遺煤溫度，K；

ρg——采空區內氣體密度，kg／m3；

ρs——采空區內遺煤密度，kg／m3；

cs——采空區內遺煤比熱容，kJ／ （kg·K）；

cg——采空區內氣體比熱容，kJ／ （kg·K）；

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

Sn——單元體內遺煤與氣體對流換熱表面積，

m2；

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

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

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

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

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

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

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

Γ——單元體的邊界。

2 采空區動態多參數溫度場流-固耦合模擬

2.1 自燃模型數值計算范圍

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

圖1 采空區坐標系及數值求解范圍

2.2 網格劃分

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

2.3 參數選取

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

2.4 主要采掘參數對采空區內溫度場的影響

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

2.4.1 工作面長度對采空區溫度的影響

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

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

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

圖3 不同工作面長度下采空區內溫度分布示意圖

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

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

2.4.2 推進速度與采空區溫度場相互關系

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

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

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

圖5 不同推進速度下采空區溫度分布

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

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

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

2.4.3 通風阻力對采空區溫度場的影響

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

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

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

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

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

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

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

3 結論

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

（2）當工作面長度增加時，采空區內的最高溫度升高，雖然采空區出現最高溫度的地點逐漸遠離工作面，但是出現自燃的可能性還是很大，所以選擇適當的工作面長度可以有效降低采空區內溫度，防止自燃事故的發生。

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

（4）隨著工作面通風阻力的增加，采空區出現的最高溫度有所升高，并且采空區出現高溫的地點向遠離工作面的方向移動，所以選擇適當的通風阻力也可以有效保證煤礦的安全生產。

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

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

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

［4］ 鄧軍，徐精彩等.國內外煤炭自然發火預測預報技術綜述 ［J］.西安礦業學院學報，1999 （4）

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

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

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

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

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

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

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