巷道圍巖不穩(wěn)定換熱系數(shù)計(jì)算方法探討

王紅兵 賈敏濤

(1.銅陵有色金屬集團(tuán)股份有限公司冬瓜山銅礦；2.中鋼集團(tuán)馬鞍山礦山研究院有限公司；3.金屬礦山安全與健康國(guó)家重點(diǎn)實(shí)驗(yàn)室；4.華唯金屬礦產(chǎn)資源高效循環(huán)利用國(guó)家工程研究中心有限公司)

?

巷道圍巖不穩(wěn)定換熱系數(shù)計(jì)算方法探討

王紅兵1賈敏濤2，3，4

(1.銅陵有色金屬集團(tuán)股份有限公司冬瓜山銅礦；2.中鋼集團(tuán)馬鞍山礦山研究院有限公司；3.金屬礦山安全與健康國(guó)家重點(diǎn)實(shí)驗(yàn)室；4.華唯金屬礦產(chǎn)資源高效循環(huán)利用國(guó)家工程研究中心有限公司)

為了更好地解決高溫礦井熱害，從圍巖導(dǎo)熱過(guò)程出發(fā)，在保證工程精度要求的前提下進(jìn)行必要假設(shè)，建立圍巖熱傳導(dǎo)微分方程，利用分離變量法求解出巷道圍巖內(nèi)部溫度分布函數(shù)；將其帶入圍巖不穩(wěn)定換熱系數(shù)定義式，求解出不穩(wěn)定傳熱系數(shù)理論計(jì)算式；最后利用Mathematica數(shù)學(xué)計(jì)算軟件進(jìn)行擬合計(jì)算，給出了不穩(wěn)定傳熱系數(shù)簡(jiǎn)化計(jì)算式，經(jīng)過(guò)擬合數(shù)據(jù)比較，簡(jiǎn)化計(jì)算式擬合誤差為0.005 990 04，證明該簡(jiǎn)化計(jì)算式可靠。

不穩(wěn)定傳熱系數(shù) Mathematica 圍巖散熱 地?zé)犷A(yù)測(cè)

研究表明，圍巖散熱是深井礦山的主要熱源，隨著礦井深度增加，所占比重還將進(jìn)一步升高。圍巖散熱量計(jì)算是確定采掘工作面降溫冷負(fù)荷、礦井降溫設(shè)備裝機(jī)冷負(fù)荷的依據(jù)，其準(zhǔn)確程度直接關(guān)系著工作面降溫冷負(fù)荷計(jì)算和礦井巷道風(fēng)流溫度的預(yù)測(cè)精度。

1 圍巖不穩(wěn)定換熱系數(shù)定義

前蘇聯(lián)學(xué)者А.Н.舍爾巴尼[1]在1953年提出，不穩(wěn)定換熱系數(shù)是指巷道圍巖深部未冷卻巖體與空氣間溫差1 ℃時(shí)，每小時(shí)從1 m2巷道壁面向(從)空氣放出(吸收)的熱量，是圍巖的熱物理性質(zhì)、井巷形狀尺寸、通風(fēng)強(qiáng)度及通風(fēng)時(shí)間等的函數(shù)。不穩(wěn)定換熱系數(shù)的解析式為:

(1)

1979年前后，我國(guó)學(xué)者岑衍強(qiáng)[2]和日本學(xué)者內(nèi)野健一等又提出無(wú)因次不穩(wěn)定換熱系數(shù)(經(jīng)時(shí)系數(shù))K的定義，可以通過(guò)無(wú)因次溫度來(lái)表達(dá)，即:

(2)

由式(1)、式(2)可知，無(wú)因次不穩(wěn)定換熱系數(shù)K與不穩(wěn)定傳熱系數(shù)Kτ存在如下關(guān)系:

(3)

由此可知，只要求得無(wú)因次不穩(wěn)定換熱系數(shù)K，便可由式(3)求得不穩(wěn)定傳熱系數(shù)Kτ。舍爾巴尼、平松良雄、岑衍強(qiáng)和內(nèi)野健一等的計(jì)算模型[3]，雖然定義方法不同，但實(shí)質(zhì)是相同的，均表明圍巖與風(fēng)流間的不穩(wěn)定換熱系數(shù)是畢歐準(zhǔn)數(shù)和傅里葉準(zhǔn)數(shù)的函數(shù)。

2 圍巖無(wú)因次不穩(wěn)定換熱系數(shù)的理論解

2.1 圍巖導(dǎo)熱微分方程

根據(jù)圍巖散熱特點(diǎn)，在保證工程精度要求前提下進(jìn)行一些假設(shè)[4]，建立圍巖的熱傳導(dǎo)微分方程式:

(4)

初始條件： 當(dāng)τ=0時(shí)，t=tr；

當(dāng)τ>0時(shí)，R→時(shí)，t→tr；

式中，R為空心圓柱體半徑，m；t0為風(fēng)溫，℃；ta為巷道壁溫度，℃。

2.2 導(dǎo)熱微分方程求解

利用分離變量法求解在上述邊界條件和初始條件下的熱傳導(dǎo)微分方程(4)，根據(jù)式(4)求圍巖的不穩(wěn)定換熱系數(shù)。

2.2.1 分離變量

設(shè)方程(4)具有兩個(gè)函數(shù)乘積形式的特解，其中一個(gè)函數(shù)X(R)只與坐標(biāo)有關(guān)，另一個(gè)函數(shù)θ(τ)只與時(shí)間有關(guān)，亦即

(5)

將式(5)帶入式(4)得

(6)

要使式(6)成立，相當(dāng)于對(duì)時(shí)間和坐標(biāo)定義域內(nèi)的任何一個(gè)τ和坐標(biāo)(R)均應(yīng)成立。即等式的左右兩邊都應(yīng)等于同一常數(shù)值D，即

(7)

(8)

2.2.2 求滿足邊界條件和常數(shù)微分方程的無(wú)窮多解

將式(7)積分得方程解為

(9)

式中，c1為積分常數(shù)。

常數(shù)D由熱傳導(dǎo)的物理意義確定：

對(duì)于溫度趨于平衡的熱過(guò)程，經(jīng)過(guò)較長(zhǎng)時(shí)間后應(yīng)達(dá)到一定的均勻溫度分布。導(dǎo)溫系數(shù)a已知為正數(shù)。因此，當(dāng)τ→時(shí)：①若D>0，則θ(τ)→，這與上述熱過(guò)程的實(shí)際物理意義不符，即溫度分布不應(yīng)趨于無(wú)窮大，所以，D不能大于0；②若D=0，則θ(τ)≡c，這與上述熱過(guò)程的物理意義也不相符，故D≠0；③若D<0，則θ(τ)隨時(shí)間的變化趨向于穩(wěn)定分布，這與研究的熱過(guò)程物理意義相符。

故D為小于0的任意常數(shù)，取D=-β2，則方程(7)和方程(8)變?yōu)?/p>

(10)

(11)

則特解(9)變?yōu)?/p>

(12)

于是熱傳導(dǎo)微分方程式(4)之特解可寫(xiě)為

(13)

解(13)符合熱傳導(dǎo)微分方程，若給予常系數(shù)c1和β以不同值，即可得無(wú)窮多的特解。依疊加原理，一般解等于符合微分方程式(4)的所有特解的總和，然后由邊界條件確定β值，由初始條件確定c1值。

式(11)可以寫(xiě)成兩個(gè)特解φ(R),ψ(R)之和，即

(14)

把式(14)帶入式(13)，得

(15)

式中，A=c1c2，B=c2c3.

要求得滿足初始條件的熱傳導(dǎo)微分方程式的一般解，可用上述特解之和，根據(jù)初始條件可以選擇A及B的適宜數(shù)值。

則熱傳導(dǎo)微分方程的一般特解為

(16)

(17)

2.3 不穩(wěn)定換熱系數(shù)理論解析式求解

根據(jù)式(2)、式(12)，按照K的定義需要對(duì)t求關(guān)于r偏導(dǎo)數(shù)，求得的K的理論公式：

(18)

(19)

式(18)為無(wú)量綱不穩(wěn)定換熱系數(shù)K與F0和Bi的理論表達(dá)式，可以看出隨著通風(fēng)時(shí)間的不斷延長(zhǎng)，K也不是恒定值。

3 圍巖無(wú)因次不穩(wěn)定換熱系數(shù)的簡(jiǎn)化計(jì)算式擬合

3.1 Mathematica數(shù)學(xué)軟件介紹

Mathematica是由美國(guó)Wolfram[6]研究公司開(kāi)發(fā)的一款廣泛使用的數(shù)學(xué)分析軟件,擁有強(qiáng)大的數(shù)值計(jì)算和符號(hào)運(yùn)算能力，也具有高精度的數(shù)值計(jì)算功能和強(qiáng)大的圖形處理功能。

Mathematica的Notebook界面下可以完成各種運(yùn)算，如函數(shù)作圖、求極限、解方程等，也可以編寫(xiě)結(jié)構(gòu)化程序。Mathematica系統(tǒng)中定義了許多功能強(qiáng)大的內(nèi)建函數(shù)(built-in-function)，這些函數(shù)可直接調(diào)用，使計(jì)算變得更加簡(jiǎn)單快捷。

3.2 擬合計(jì)算

無(wú)因次不穩(wěn)定傳熱系數(shù)的理論解是由特殊函數(shù)復(fù)合而成的，無(wú)法直接計(jì)算，可通過(guò)Mathematica數(shù)學(xué)軟件擬合不穩(wěn)定傳熱系數(shù)的簡(jiǎn)化計(jì)算式。

讓Bi以5為步長(zhǎng)從5增加到55，F(xiàn)0從0到20連續(xù)變化，利用Mathematica軟件畫(huà)出K(Bi，F(xiàn)0)的二維曲線列表圖2和三維曲面圖3。

圖1 Bi=55時(shí)方程(19)的左端曲線

圖2 K(Bi，F(xiàn)0)的二維曲線列表

圖3 K(Bi，F(xiàn)0)的三維曲面

以1為步長(zhǎng)對(duì)變量F0離散化，計(jì)算出11×20組K(Bi，F(xiàn)0)的三維數(shù)據(jù)。以此為基礎(chǔ)，利用Mathematica數(shù)學(xué)軟件中的二元擬合方法對(duì)K(Bi，F(xiàn)0)進(jìn)行擬合，得到如下擬合函數(shù)

(20)

式中，c1=-1.743 983,c2=0.011 916,c3=-0.000 323,c4=2.849 676×10-6,c5=0.013 902,c6=1.675 393,c7=-1.223 007。

將(Bi，F(xiàn)0)的離散取值代入式(20)，計(jì)算出相應(yīng)的K取值，再與K的擬合數(shù)據(jù)相比較，可以得出函數(shù)K(Bi，F(xiàn)0)的擬合誤差為0.005 990 04，該擬合精度較高。

圖2表明，在巷道掘成的通風(fēng)初期F0很小，K值隨著B(niǎo)i的增大而迅速增大，表明巖體與風(fēng)流活躍地進(jìn)行熱交換，這時(shí)對(duì)流換熱系數(shù)h值在熱交換中占主導(dǎo)地位。隨著通風(fēng)時(shí)間延長(zhǎng)(F0增大)，K值也遞減，F(xiàn)0增加到某一程度后，即使Bi增長(zhǎng)，K值的增長(zhǎng)幅度也很小，說(shuō)明這時(shí)圍巖與風(fēng)流已經(jīng)進(jìn)行了較充分的熱交換，冷卻帶的擴(kuò)大趨向變?nèi)酰锏辣诿娴臏囟戎饾u接近風(fēng)溫，這時(shí)巖石中的熱阻比壁面熱阻對(duì)K值影響更大。

4 結(jié) 語(yǔ)

(1)運(yùn)用分離變量法求解了導(dǎo)熱微分方程，給出了不穩(wěn)定換熱系數(shù)的解析式，具有重要的理論意義。

(2)采用Mathematica數(shù)學(xué)軟件中的二元擬合方法，給出了不穩(wěn)定換熱系數(shù)的簡(jiǎn)化計(jì)算公式，經(jīng)過(guò)擬合數(shù)據(jù)比較，簡(jiǎn)化計(jì)算公式擬合誤差為0.005 990 04，該擬合精度較高，為后續(xù)進(jìn)行準(zhǔn)確的工作面降溫冷負(fù)荷計(jì)算和礦井巷道風(fēng)流溫度預(yù)測(cè)提供了重要保障。

[1] А.Н.舍爾巴尼，О.А.克列姆涅夫，В.Я.茹拉夫連科.礦井降溫指南[M].黃翰文，譯.北京：煤炭工業(yè)出版社，1982.

[2] 岑衍強(qiáng)，胡春勝，候祺棕．井巷圍巖與風(fēng)流間不穩(wěn)定換熱系數(shù)的探討[J]．阜新礦業(yè)學(xué)院學(xué)報(bào)，1987，6(3):105-114．

[3] 平松良雄．通風(fēng)學(xué)[M]．劉運(yùn)洪譯．北京：冶金工業(yè)出版社，1981．

[4] 劉何清，何昌富，楊 威，等.巷道變溫圈內(nèi)溫度分布及不穩(wěn)定傳熱系數(shù)求解方法[J]．湖南科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2015(4)：7-13.

[5] 王義江．深部熱環(huán)境圍巖及風(fēng)流傳熱傳質(zhì)研究[D]．徐州：中國(guó)礦業(yè)大學(xué)，2010．

[6]Baumann,Gerd.MathematicaForTheoreticalPhysics[M]．2ndedition.Berlin：SpringerVerlag，2005.

Discussion on the Calculation Method of the Instability Heat Transfer Coefficient of the Surrounding Rock of Roadway

Wang Hongbing1Jia Mintao2,3,4

(1.Dongguashan Copper Mine,Tongling Nonferrous Metals Group Holding Co.,Ltd.;2.Sinosteel Maanshan Institute of Mining Research Co.,Ltd.;3.State Key Laboratory of Safety and Health for Metal Mines;4.Huawei National Engineering Research Center of High Efficient Cyclic and Utilization of Metallic Mineral Resources Co.,Ltd.)

In order to solve the thermal pollution problem of high temperature mine effectively,firstly,based on the heat conduction process of surrounding rock,the necessary assumptions is done on the premise of guarantee engineering accuracy requirement to establish the heat conduction differential equation,the temperature distribution function of the internal surrounding rock of roadway is solved by adopting method of separation variables;then,the temperature distribution function of the internal surrounding rock of roadway is taken into the definition equation of instability heat transfer coefficient,the theoretical calculation equation of instability heat transfer coefficient is obtained;finally,the fitting calculation is conducted by using Mathematica calculation software,the simplified calculation equation of instability heat transfer coefficient is given.The comparison results of the fitting data show that the fitting error of the simplified calculation equation is 0.005 990 04,it is indicated that the simplified calculation equation of this paper is reliable.

Instability heat transfer coefficient,Mathematica,Surrounding rock heat dissipation,Geothermal prediction

2016-09-01)

王紅兵(1969—)，男，工程師，244031 安徽省銅陵市獅子山區(qū)。