有熱源巷道風(fēng)流溫度分布非穩(wěn)定模型及數(shù)值計(jì)算*

李宗翔，劉　宇，劉漢武，林　琳

(1.遼寧工程技術(shù)大學(xué) 安全科學(xué)與工程學(xué)院, 遼寧 阜新 123000；2.礦山熱動(dòng)力災(zāi)害與防治教育部重點(diǎn)實(shí)驗(yàn)室, 遼寧 阜新 123000 )

0　引言

有熱源巷道風(fēng)流溫度分布非穩(wěn)定計(jì)算，是災(zāi)變時(shí)期(熱害或火災(zāi))礦井通風(fēng)網(wǎng)絡(luò)溫度仿真計(jì)算的基礎(chǔ)和關(guān)鍵性問題[1-5]。有熱源巷道指的是巷道或工作面中設(shè)有大功率機(jī)電設(shè)備(放熱)或制冷設(shè)備(吸熱、負(fù)熱源)或巷道火災(zāi)火源放熱等，屬于典型流體分支中源的問題[6]。當(dāng)風(fēng)流從熱源經(jīng)過，溫度發(fā)生變化，熱源升溫，冷源降溫。溫度在井巷風(fēng)流中的傳播屬于有熱源的風(fēng)流對流-擴(kuò)散-換熱非穩(wěn)定過程。如果是單一條巷道或局部網(wǎng)巷道問題，可運(yùn)用Fluent軟件對3D巷道空間模型進(jìn)行精度更高的仿真計(jì)算[7]，多數(shù)火災(zāi)的仿真計(jì)算是按公式法計(jì)算溫度[8-9]。有熱源且非穩(wěn)定問題，初始溫度不是均勻分布就不能簡單套用解析式來計(jì)算。更復(fù)雜情況，如機(jī)電設(shè)備開停熱源時(shí)斷時(shí)續(xù)、火災(zāi)燃燒強(qiáng)度不斷變化或風(fēng)流紊亂等，都將形成非穩(wěn)定-非均勻的溫度分布，此時(shí)，必須用二階微分方程來描述，用差分?jǐn)?shù)值方法求解[9-12]。在礦井火災(zāi)計(jì)算中，火源都是按溫度值(1類條件)給出[11-13]。這里討論有熱源巷道非穩(wěn)定風(fēng)流溫度二階微分方程及有限元方法求解，解決熱源按放熱強(qiáng)度(2類熱流量條件)給出的熱源項(xiàng)處理算法，并將算法組構(gòu)到網(wǎng)域中，實(shí)現(xiàn)對全礦井風(fēng)流溫度分布計(jì)算等一系列問題做進(jìn)一步研究。

1　有熱源單一巷道風(fēng)流溫度非穩(wěn)定基本方程

被加熱的風(fēng)流(或高溫?zé)熈?一方面隨風(fēng)流整體移動(dòng)，另一方面由于溫度梯度和風(fēng)流脈動(dòng)作用，在風(fēng)流中逐漸彌散，產(chǎn)生溫度傳導(dǎo)和擴(kuò)散，同時(shí)熱風(fēng)流又與巷道圍巖產(chǎn)生熱交換。

圖1　通風(fēng)巷道熱源與風(fēng)流溫度分布Fig.1　Distribution of lane heat source and temperature in the ventilation tunnel

如果是在全礦復(fù)雜網(wǎng)域的背景下計(jì)算單一巷道問題，仍離不開一維井巷模型。通風(fēng)巷道熱源與風(fēng)流溫度分布情況如圖1所示，如果一維井巷通風(fēng)風(fēng)流的初始溫度分布是均勻、穩(wěn)定條件的計(jì)算，可根據(jù)風(fēng)流與圍巖熱交換，建立風(fēng)流溫度一階線性非齊次微分方程，求解方程導(dǎo)出穩(wěn)定溫度分布的析解式[1]。

從礦井網(wǎng)絡(luò)系統(tǒng)宏觀的角度看，單一巷道風(fēng)流運(yùn)動(dòng)可看做是一維的巷道通風(fēng)對流擴(kuò)散問題，礦井網(wǎng)絡(luò)系統(tǒng)火災(zāi)仿真計(jì)算仍然依據(jù)這一假設(shè)[10-12]，進(jìn)一步假設(shè)，熱源引起風(fēng)流溫度變化的時(shí)間，遠(yuǎn)小于長期的通風(fēng)與地層之間非穩(wěn)定熱交換時(shí)巷道圍巖溫度變化時(shí)間，計(jì)算時(shí)視巷道圍巖溫度不變化，并忽略熱濕換熱影響。那么，用溫度表示的熱交換對流擴(kuò)散微分方程為[11-12]：

(1)

式中：t為風(fēng)流溫度，℃；x為長度位置，m；kx為溫度擴(kuò)散系數(shù)，m2/s；CP為風(fēng)流定壓比熱容，J/(kg℃)；ρ為風(fēng)流密度，kg/m3；tr為巷道圍巖溫度，℃；qS為巷道中單位時(shí)間單位體積的放熱源項(xiàng)，J/(m3·s)；若qS以J/(ms)單位形式給出，除以巷道斷面積S；βf為單位時(shí)間單位體積風(fēng)流與巷道之間對流熱交換源項(xiàng)值，W/(℃· m3)，即平均到單位體積中換熱那部分值;τ為時(shí)間變量，s。

熱源qS分為穩(wěn)定熱源和非穩(wěn)定熱源，一般穩(wěn)定工作的機(jī)電設(shè)備散熱是穩(wěn)定熱源，礦井火災(zāi)的火源燃燒強(qiáng)度在不斷變化，一般屬于非穩(wěn)定熱源。熱源項(xiàng)大都是按第1類邊界條件，即熱源點(diǎn)按溫度給出的；還有一種形式是以熱源項(xiàng)給出，火源燃燒是能量釋放的過程，從能量釋放加以解釋時(shí)，應(yīng)按第 2類邊界條件，即放熱量給出。一般，穩(wěn)定熱源qS= const，非穩(wěn)定熱源用函數(shù)表示qS(x,τ)。

隨著含煙風(fēng)流的流動(dòng)，熱煙流與巷道圍巖存在熱交換。巷道風(fēng)流微元體如圖2所示，對一維巷道計(jì)算模型，有微分元體dx，微分元體積為 dV=Sdx，在dV中的熱交換的總熱流量取決于巷道壁的換熱，即:

βr=arUdx

式中:ar為換熱系數(shù)，W/(m2℃)；U為巷道周長，m。

圖2　巷道風(fēng)流微元體Fig.2　Lane heat source microelement

平均分?jǐn)偟絾挝惑w積中換熱那部分值定義為:

代入式(1):

(2)

或可寫成:

(2a)

式中:ar為巷道壁面換熱系數(shù)，W/(m2·℃)；Qt, s為井巷中單位長度上熱源升溫放熱強(qiáng)度，J/(m·s)。

方程(2a)各項(xiàng)的單位量綱統(tǒng)一于[℃·s-1]。

熱煙氣流的機(jī)械彌散伴隨著溫度的傳播，此外，即使不存在煙氣流彌散，也存在風(fēng)流空氣的熱傳導(dǎo)，顯然，熱煙氣流的機(jī)械彌散遠(yuǎn)大于空氣熱傳導(dǎo)。方程中的溫度擴(kuò)散系數(shù):

式中:λ為熱傳導(dǎo)系數(shù)，W/(m·℃);Ex為風(fēng)上流縱向熱煙氣流湍流彌散熱傳系數(shù)，W/(m·℃)。

2　巷道對流-擴(kuò)散問題的一維有限元分析

研究巷道一維對流-擴(kuò)散(彌散)問題可用的方法有差分方法和有限元法，這里只討論比較簡單的有限元法[13]。若將巷道按一維考慮，根據(jù)一維有限元單元的劃分，對單元e有插值函數(shù):

Φi(x) =ai+bix，(i= 1, 2 )

(3)

式中:ai,bi為插值系數(shù)，在每個(gè)單元中應(yīng)保證Φi(xi) =1，Φi(xj) =0，Φj(xj) =1，Φj(xi) =0，得到:

單元中濃度函數(shù)t=Φiti+Φjtj，其中ti,tj為待求的節(jié)點(diǎn)濃度(函數(shù))值。

由對流擴(kuò)散方程(1)，相應(yīng)的卡遼金積分表達(dá)式為 ：

根據(jù)Green-Gauss公式，對其分部積分 ：

單元近似函數(shù)t(τ) =ti(τ)Φi(x) 。對e單元，有限元特征式 ：

(4)

其中 ：

將式(3) 代入上式，積分得到 ：

(5)

(5a)

式中：[A]e，[B]e，[F]e和We分別為單元的矩陣、源向量及右端項(xiàng)向量；e代表單元。

以上單元分析通過單元合成將單元疊加到整體一維區(qū)域上，得到有限元求解的線性方程組(過程略)，即得到合成一維巷道中的總體方程為：

(6)

式中： [A]，[B]，[F]和W分別為由各單元合成的巷道總體的矩陣、源向量及右端項(xiàng)向量。

在MATLAB平臺上實(shí)現(xiàn)了上述過程一條巷道的一維有限元計(jì)算機(jī)程序編制工作(TF_GasTen1a.m)。

3　單一巷道中風(fēng)流非穩(wěn)定溫度分布模擬計(jì)算

3.1　計(jì)算的基本條件

巷道總長度為428 m，巷道參數(shù)：巷道斷面積S=11.2 m2，周長U=13 m，α=0.021 kg/m3，通風(fēng)風(fēng)量為15.1 m3/s，風(fēng)速為1.35 m/s。起點(diǎn)邊界條件為27 ℃，風(fēng)流密度1.2 kg/m3，風(fēng)流熱傳導(dǎo)系數(shù)0.027 W/(m·℃)；風(fēng)流定壓比熱1 005 J/kg℃，圍巖溫度27 ℃；廣義換熱系數(shù)7 W/(m2·℃)，Ex=1.35 m2/s。單元剖分尺度1.5 m(精度)；熱源點(diǎn)位置130 m，并假設(shè)：每件放熱機(jī)電設(shè)備(熱源)分配在一個(gè)單元的空間中，則We中的元素 ：

(7)

其中，h是單元長度，m。

3.2　巷道初始溫度恒定均勻分布的定解條件及求解

方程(2a)溫度分布t(x,τ) 求解的定解條件為 ：

邊界條件：t|x=0=t(0,τ) = 27℃；

初始條件：t|τ=0=t(x, 0) = 27,0

熱源項(xiàng)：Qs(130,τ) = 300 kW。

溫度的初始分布恒定，起點(diǎn)邊界條件是巷道入口風(fēng)流溫度始終保持27 ℃，有熱源巷道風(fēng)流溫度變化過程模擬結(jié)果如圖3所示。假設(shè)在流動(dòng)過程中井巷風(fēng)路中無外源氣涌出，熱源強(qiáng)度(機(jī)電設(shè)備功率)為300 kW。模擬時(shí)間步長取2 s，模擬結(jié)果如圖4所示，給出溫度的非穩(wěn)定模擬的動(dòng)態(tài)解。根據(jù)圖3的模擬結(jié)果，風(fēng)流在第175 s以后就基本趨于穩(wěn)定，為使之達(dá)到充分穩(wěn)定，這里以在第320 s時(shí)的溫度分布作為對比，可見已經(jīng)完全穩(wěn)定并重合。圖3(f)中的粗線是按解析法給出的解(對比)。

圖3　溫度均勻初始分布有熱源巷道風(fēng)流溫度變化過程的模擬結(jié)果Fig.3　The simulation results of temperature variation of lane with heat source

3.3　巷道初始溫度分布非恒定的定解條件及求解

在熱源穩(wěn)定的條件下，有熱源巷道風(fēng)流溫度分布最終趨于穩(wěn)定(見圖3)，但溫度不是均勻分布的。在礦井火災(zāi)的災(zāi)變時(shí)期，由于風(fēng)流紊亂，火災(zāi)煙流在復(fù)雜網(wǎng)絡(luò)中傳播，在某一時(shí)間點(diǎn)上井巷風(fēng)流的溫度分布也不是均勻的，存在較大的波動(dòng)狀態(tài)。將圖3中第320 s 趨于穩(wěn)定的溫度分布結(jié)果作為初始條件，將起點(diǎn)的邊界條件改為32℃，熱源強(qiáng)度改為100 kW，其他模擬條件不變，則方程(2a)求解的定解條件為 ：

邊界條件：t|x=0= 32℃；

初始條件 ：t|τ=0=t(x, 0)，0

熱源項(xiàng)：Qs(130,τ) = 100 kW。

風(fēng)流溫度計(jì)算結(jié)果如圖4所示。

在溫度邊界條件和源改變后，在非穩(wěn)定計(jì)算中，井巷溫度分布重新經(jīng)過一番調(diào)整變化過程，最后又達(dá)到新的穩(wěn)定狀態(tài)，圖4(f)中的粗線溫度分布是按解析法給出的一般解，顯然，解析法無法實(shí)現(xiàn)有熱源、非均勻分布的井巷風(fēng)流溫度分布的計(jì)算。

3. 4　有限元數(shù)值解的收斂性分析

在解的收斂性檢驗(yàn)方面，理論上，假設(shè)熱源前的風(fēng)流溫度是不受熱源影響的，圖3和圖4中的最后的一幅圖給出了穩(wěn)定解析解與非穩(wěn)定數(shù)值計(jì)算達(dá)到穩(wěn)定后的數(shù)值解的對比，熱源前一段的數(shù)值解與解析解兩者完全重合，說明有限元數(shù)值解的收斂性很好，計(jì)算精度滿足要求。

圖4　溫度非均勻初始分布有熱源巷道風(fēng)流溫度變化過程的模擬結(jié)果Fig.4　The simulation results of the temperature change process of the lane with heat source

4　結(jié)論

1) 數(shù)值模擬描述了有熱源巷道風(fēng)流溫度的平移和擴(kuò)散的流動(dòng)傳播過程。將井巷中設(shè)備熱源(或制冷負(fù)源)作為源項(xiàng)，以熱強(qiáng)度(熱流量)條件給出，構(gòu)建了巷道風(fēng)流溫度分布非穩(wěn)定變化求解方程的熱源表達(dá)方法。數(shù)值計(jì)算實(shí)現(xiàn)了有熱巷道溫度分布隨時(shí)間和位置的變化。強(qiáng)調(diào)指出，隨著熱源能量的傳播，經(jīng)過長時(shí)間或傳播距離的增長而逐漸趨于穩(wěn)定；礦井熱害主要決定因素還是地層溫度，設(shè)備熱源對風(fēng)流溫度影響與變化最終穩(wěn)定于地層溫度。

2) 提高計(jì)算精度是仿真計(jì)算的基本要求。數(shù)值計(jì)算的誤差原因是時(shí)間和空間剖分離散的誤差導(dǎo)致的，考慮到井巷各地點(diǎn)的分支風(fēng)量變化很大，宜于采用自動(dòng)化識別收斂性判定條件，確定剖分精度，算例中剖分單元長度1.5 m，取模擬時(shí)間步長 2 s 不發(fā)生振蕩。

3) 仿真計(jì)算模型(程序)可以規(guī)劃出井巷中的熱源的影響范圍，適用于井下熱害防治的工程設(shè)計(jì)與計(jì)算。以上算法同樣適用于非穩(wěn)定熱源，如設(shè)備開停以及火災(zāi)火源強(qiáng)度的變化屬于非穩(wěn)定熱源模型，在非穩(wěn)定計(jì)算過程中，隨時(shí)間調(diào)整熱源值，具體實(shí)現(xiàn)有待于進(jìn)一步研究。

[1]王省身，張國樞. 礦井火災(zāi)防治[M]. 北京：煤炭工業(yè)出版社, 1990.

[2]劉鐵民. 突發(fā)事件應(yīng)急預(yù)案體系概念設(shè)計(jì)研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2011, 7(8): 5-13.

LIU Tiemin. Design of the emergency plan system’s concept[J]. Journal of Safety Science and Technology, 2011, 7(8): 5-13.

[3] 盛偉，鄭海坤. 人體降溫服在礦井熱環(huán)境中的應(yīng)用綜述[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2013, 9(12): 95-101

SHENG Wei, ZHENG Hai kun. Literature review on application of body cooled suits in mine thermal environment[J]. Journal of Safety Science and Technology, 2013, 9(12): 95-101.

[4]宋遠(yuǎn)卓，代昌標(biāo)，陸愈實(shí). 煤礦火災(zāi)過程仿真模擬系統(tǒng)的研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2006,2(5):55-58.

SONG Yuanzhuo, DAI Changbiao, LU Yushi. Research on simulation system for coal mine fire process[J]. Journal of Safety Science and Technology, 2006,2(5):55-58.

[5]宋曉燕,謝中朋.三河尖煤礦熱力參數(shù)測定及熱濕源危害分析[J].中國安全生產(chǎn)科學(xué)技術(shù),2012,8(12):174-178.

SONG Xiaoyan, XIE Zhongpeng. Heat parameter determination of sanhejian mine and research of heat-harm cause[J]. Journal of Safety Science and Technology, 2012,8(12):174-178.

[6]李宗翔. 有源風(fēng)網(wǎng)模型及其應(yīng)用計(jì)算[J]. 煤炭學(xué)報(bào), 2010, 35(S1): 118-122.

LI Zongxiang. Containing the source ventilation network model and its application[J]. Journal of China Coal Society, 2010, 35(S1): 118-122.

[7]王凱，蔣曙光，張衛(wèi)青，等. 礦井火災(zāi)應(yīng)急救援系統(tǒng)的數(shù)值模擬及應(yīng)用研究[J]. 煤炭學(xué)報(bào), 2012, 37(5):857-862

WANG Kai, JIANG Shuguang, ZHANG Weiqing, et al. Numerical simulation and application research of mine fire emergency rescue system[J]. Journal of China Coal Society, 2012, 37(5): 857- 862.

[8]賈進(jìn)章，馬恒，劉劍. 礦井火災(zāi)時(shí)期溫度分布數(shù)值模擬[J]. 遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2003, 22(4): 460-462.

JIA Jinzhang，MA Heng，LIU Jian. Numerical simulation for temperature distribution during mine fire period[J]. Journal of Liaoning Technical University，2003,22(4):460-462.

[9]丁翠. 火災(zāi)時(shí)期礦山通風(fēng)巷道“特征環(huán)”與“關(guān)鍵環(huán)”的分布特征研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2017,13(8):55-60.

DING Cui.Study on distribution characteristics of “characteristic ring” and “key ring” in mine ventilation roadway during fire[J]. Journal of Safety Science and Technology,2017,13(8):55-60.

[10]張興凱. 礦井火災(zāi)時(shí)期風(fēng)流動(dòng)態(tài)過程的計(jì)算[J]. 煤礦安全, 1990, 20(6): 55- 60

ZHANG Xingkai. Mine fire during the dynamic process of calculation[J]. Safety in Coal Mines, 1990, 20(6): 55- 60.

[11]戚宜欣. 礦井火災(zāi)煙流溫度場及濃度場的數(shù)值模擬[J]. 西安礦業(yè)學(xué)院學(xué)報(bào), 1994(1):26-33

Qi Yixin. Numerical simulation of temperature and concentration fields of mine fire fume[J]. Journal of Xi’an Mining Institute, 1994(1): 26- 33.

[12]李宗翔，王雅迪，李林. 上行風(fēng)流火災(zāi)3D礦井通風(fēng)系統(tǒng)災(zāi)變過程仿真[J]. 煤炭學(xué)報(bào), 2015, 40(1): 115-121

LI Zongxiang, WANG Yadi, LI Lin. 3D simulation of disaster process in mine ventilation system during fire period[J]. Journal of China Coal Society, 2015, 40(1): 115-121.

[13]章本照. 流體力學(xué)中的有限元方法[M]. 北京：機(jī)械工業(yè)出版社, 1986.