礦井氣候多參數預測與通風網絡自動解算算法研究*

馬　恒，高　巍，周西華

(1. 遼寧工程技術大學 安全科學與工程學院，遼寧 阜新 123000； 2. 礦山熱動力災害與防治教育部重點實驗室，遼寧 阜新 123000)

0　引言

自通風網絡解算方法被提出以來，通過國內外學者不斷研究，該領域的技術日漸成熟，開發出大量的礦井通風網絡仿真軟件。但許多軟件對礦井氣候的考慮仍停留在靜態的層次上，即需要人事先測量出各網絡節點礦井氣候參數，然后再得出自然風壓。過程繁瑣耗時長，工作量大。如何進行快速、簡單且不失準確的計算是目前需要解決的問題。本文通過對風流的壓強和溫度變化規律的研究，在已知進風口風流溫度和壓強情況下，預測巷道各點風流溫度、壓強，進一步得出各巷道平均密度，得出礦井含自然風壓的通風網絡動態解算算法，并基于MATLAB軟件對算法進行程序設計，對算法進行驗證。

1　井巷風流密度計算

1.1　風流平均密度

井巷風流密度是影響自然風壓的一個重要參數。氣體的可壓縮性會對風流的密度造成變化，這個變化約為6%～8%[1]。由于影響較小，在一般情況下，將礦內風流視為非壓縮性流體來對待，即近似視為等容過程。由此平均密度計算公式[2]：

(1)

式中：p1，p2分別為始末節點空氣的壓強，Pa；ρ1，ρ2分別為始末節點空氣的密度，kg/m3；ρm為平均密度，kg/m3。

由公式可知，當空氣被設為不可壓縮流體時，密度的變化主要由壓強來決定。上述平均密度計算公式過于繁瑣，對于平均密度的計算，可以用兩端點的算術平均值[3]：ρm=(ρ1+ρ2)/2來代替，滿足計算要求。

1.2　風流密度計算

計算井下濕空氣密度[4]：

ρ=0.003 48(P-0.377 9φPs)/T

(2)

式中：P為濕空氣的壓強，Pa；T為濕空氣的絕對溫度，K；φ為濕空氣的相對濕度，%；Ps為飽和水蒸氣分壓力，Pa。

由公式可知，φ由測量的干濕溫度得來，Ps可根據T查表得出，因此，只要知道礦井巷道各點風流的溫度T和壓強P就可得出各個分支節點密度。這通常需要實地測量濕空氣壓力P，干球溫度td和濕球溫度tw才能得出。

2　礦井氣候多參數預測下的自然風壓算法的實現

現有自然風壓的計算均是通過實地測量濕空氣壓力，干球溫度和濕球溫度來進行的[5]，由于測點與數據較多，其過程不但繁瑣還無法達到動態計算的目的。對此，筆者考慮在已知井口大氣壓強p0和干濕溫度t0d和t0w的情況下，預測出井下各個節點的壓強和溫度，再通過各個節點的壓強和溫度能得出各個節點的風流密度。因此需要引入各節點風流溫度和各節點壓強計算公式。

2.1　風流溫度的計算

(1)井筒風流溫度的計算方程為[6]：

(σ-g/CPL)Qρ0 sCPL/λtK(α)

(3)

式中：T0m為井底(井筒末端)溫度，K；T0 s為地面(井筒始端)溫度，K；σ為地溫梯度，不同地區地溫梯度不同；Z為井筒始末兩端高差，m；L為巷道長度，m；λ為導熱系數，W/(m·K)。根據文獻[6]查表得出；K(α)為不穩定換熱系數，通過畢歐準數和傅立葉準數查表得出；Q為風流體積流量，m3/s；ρ0s為巷道入口的風流密度，kg/m3；CPL為干空氣的定壓比熱，一般取1 005 J/kg·K；g為重力加速度，m2/s，一般取9.81。

(2)普通巷道風流溫度的計算方程為[6]：

(4)

式中：Tgu為原巖絕對溫度，K；U為巷道周長，m；γ0為巷道當量半徑，即巷道有效截面和積濕周之比定義為水力半徑，當量半徑定義為二倍的水力半徑，m；根據文獻[6]，Δhp為由于各種熱源引起的風流熱焓增值，J/kg。

2.2　風溫預測的準確性分析

由公式(3)和(4)可知，風溫預測需考慮的參數過多，因此，保證每個節點計算的溫度可靠性是有必要考慮的。首先，干空氣的定壓比熱CPL在溫度250～300 K的范圍內，其數值變化在1 003～1 005之間，所以取1 005并不影響結果。其次，式中可以準確得到的參數：進風井口溫度T0 s(在井口容易測)，高差Z，巷道長度L，周長U，當量半徑γ0。在礦井建設的過程中這些參數在圖紙上都會標出，不需要再次測量。而原巖溫度Tgu，地溫梯度σ，導熱系數λ和不穩定換熱系數K(α)都與地質條件和巖石物理性質有關，因為建井前需要進行地質條件、水紋條件和巖層情況的考察，原巖溫度和地溫梯度都可以直接得出，而導熱系數λ和不穩定換熱系數K(α)需要知道巷道處在哪個巖層中，了解巖石性質，通過計算和查表也可以得出[6]。唯一不確定的參數就是風量Q和節點密度ρ0s，由于不能實測且在下文的網絡解算的過程中這2個參數數值不斷變化，因此會對風溫預測結果造成影響，由公式(3)和(4)可知，Q*ρ可作為一個參數來考慮，以圖1為例，假設其他的參數都已確定且不變，當進風口Q*ρ分別等于30，40，50時，各節點溫度值。分析結果見表1。

圖1　礦井通風系統簡Fig.1　Diagram of mine ventilation system

Q*ρT/KV1V2V3V4V530278.0141280.2022280.4145282.1756281.269840278.0144280.1572280.4067282.1422281.236650278.0145280.1296280.3975282.1202281.2146最大絕對差0.00040.07260.01700.05540.0552

結論：Q*ρ分別為30，40，50的情況下,溫度變化在0.000 4～0.072 6的范圍內，變化很小。因此風量Q和節點密度ρ對風溫的影響很小。而真正對風溫影響很大的高差，原巖溫度和始點溫度都能夠通過各種方式準確的得到，可以保證進回風井和普通巷道的溫度預測的準確性，至于擁有額外熱源的特殊巷道，文獻[6]、[7]介紹很清楚，主要考慮機器運作產生熱量Δhp，而產熱量可以根據機器的銘牌和實際運作功率等計算出。

2.3　各點壓強的計算

點壓強計算公式[7]：

(5)

式中：p1，p2分別為始末節點空氣的壓強，Pa；ρ1s，ρ1m為始末節點空氣的密度，kg /m3；g為重力加速度，m2/s，一般取9.81；Q1為風流體積流量，m3/s；Z1-2為始、末節點標高差，m。

2.4　自然風壓的計算

由前文可知，平均密度由ρm=(ρ1+ρ2)/2來代替。得到最終的自然風壓計算公式[8]：

(6)

式中：hn為礦井自然風壓，Pa；ρs，ρm為通路上各分支的始節點、末節點的密度，kg/m3；Zs-m為始、末節點標高差，m。

通過式(2)、(5)和(6)可以得出：

0.377 9φiPs)/Tis+ ((Pis+ 0.001 74(Pis-

(7)

3　通風網絡實時解算與程序實現

3.1　Cross迭代算法改進

Cross法，即斯考德-恒斯雷法[9]。其運算過程：按規定選擇獨立回路，運用高斯迭代法逐次求解回路內的風量修正量Δqi，直到滿足精度ε為止，見公式(8)和(9)。

(8)

(9)

在原有Cross迭代法的基礎上，加入上文自然風壓計算方法和原理[10-15]，實現含有自然風壓的通風網絡實時解算。在運算過程中，還需考慮以下幾個問題：

1)由式(7)、(8)和(9)可知，將各節點溫度和壓強同時進行迭代，存在迭代不收斂可能。為了避免此情況的發生，在程序編制中，對計算方法進行簡化。首先在不考慮自然風壓的情況下進行一次迭代計算，通過此計算得出各分支風量值，該分支風量作為計算各節點的風流溫度初始參數，在之后的迭代計算中將一直使用此溫度不變。

2)因為網絡各節點壓強和密度的計算均需已知彼此才能得出，在兩者均未知情況下，為了保證計算能進行，在初次迭代時需擬定各節點密度值，即假設各節點初始密度均等于井口空氣密度，這樣就可以計算出各點壓強，再利用得到的點壓強在下次迭代中得出各節點密度，依次循環。注意這個擬定的初始密度值只在初次計算應用。

3)需考慮濕度對自然風壓的影響，經驗可知[4]，濕度與巷道壁的潮濕度和蒸發率有關，這些變化都是復雜難測，因此，只能得出范圍值，現將通風系統分為進風段、用風段、回風段3個部分來簡化。相對濕度在井下變化趨勢：夏季，風流的相對濕度一直增加狀態，到回風段時會達到90%～100%；在冬季，風流的相對濕度剛開始會減小或不變，但到用風段與夏季變化趨勢一樣，到回風段時也會達到90%～100%。自然風壓預測的另一個主要因素：高差，高差主要在進風井和回風井，而用風段幾乎不存在高差，對自然風壓的預測中影響很小。進風井的平均相對濕度應該高于進風口風流相對濕度的0%～20%。除了青藏高原與新疆部分盆地，我國其他地區空氣的相對濕度都達到60%以上。因此，假如進風口處的相對濕度為60%，那么進風側的平均相對濕度應該在60%～80%之間，回風段的相對濕度在90%～100%之間，而用風段相對濕度就應該在80%～90%之間。取中間數就是進風段：70%；用風段：85%；回風段95%。應用取范圍中間數的相對濕度預測方法，可保證預測的相對濕度與真實的相對濕度誤差上下不超過10%。當然，這里的預測經驗主要應用于我國煤礦，煤礦絕大多數均為潮濕巷道。誤差分析：例如，溫度為0，5，15，20，25℃，巷道高差為100 m，壓強為100 000 Pa，假設預測相對濕度為70%而實際上為60%，密度與自然風壓計算誤差結果見表2。

表2　計算結果對比Table 2　Comparison of results

結論：由公式2和6，礦井溫度0～25℃的相對濕度每相差10%，風流密度就相差0.000 3～0.015，自然風壓每一百米的高差就產生0.29～1.45 Pa的誤差，自然風壓誤差在0.024%～0.116%，。因此采用分3段預測和取中間值的方法代替實測每個點的相對濕度，可以最大程度上減少誤差，誤差值對自然風壓的影響在可控的范圍內。受自然風壓影響越大的礦井預測越準確，誤差控制在1%以下。此方法為相對濕度對自然風壓造成誤差降到最低。

其算法程序流程為：

1)已知：G=(V,E)、總風量、分支風阻、初擬風量和各個分支始末兩端節點的高差。

2)對通風網絡G的分支和節點進行編號，擬定樹與余樹。

3)擬定基本回路矩陣：C。

4)在不考慮自然風壓和通風機風壓的情況下，對網絡各個分支進行初步分風(其目的是為初次擬定溫度和點壓強提供風量)。

5)已知：井口大氣壓p0與大氣溫度T0,求得井口密度ρ0。并根據公式(3)和(4)確定網絡各個節點的溫度。

6)先擬定網絡各節點密度ρi均為ρ0，并帶入公式(5)，通過公式(7)、(8)和(9)對其進行迭代運算。(考慮濕度影響)

7)檢查各參與迭代的回路風量修正量Δqi，若max(|Δqi|)<ε滿足精度，則進入步驟9，否則轉回6。

8)k=k+1。(k為迭代次數)

9)結束迭代，輸出結果

3.2　程序驗證

如圖2建立的一個簡單礦井模型。表3為程序是否含有自然風壓的2種情況下解算運行結果對比表，運算結果證明：雖然在不考慮自然風壓的情況下解算結果較準確，但與考慮自然風壓的情況下的結果相比還存在一定差距，因此在解算中對自然風壓的考慮是有必要的。

圖2　簡單礦井模型Fig.2　MVIS ventilation system network settlement

如圖2所示，已知井口壓強和溫度分別為p0=100 000 Pa，T0=275 K，各分支標高差Z1=300 m，Z2=20 m，Z3=0 m,Z4=20 m ,Z5=30 m,Z6=50 m,Z7=250 m。認為所有巷道長度L=500 m,R1=R2=R3=R4=R5=0.03N·s2·m-8，R6=R7=0.06 N·s2·m-8,導熱系數λ均為2 W/(m·K)，周長均為15 m，當量半徑γ0=3 m，設e1，e2，e3為進風段，e4為用風段，e5，e6，e7為回風段，由上文的相對濕度預估經驗，節點1到節點5相對濕度φ分別為70%，85%，85%，95%，95%，Δhp1-2=Δhp1-3=Δhp3-4=Δhp2-4=2 010 J，hf=30 Pa，Q1=Q7=30 m3/s。

表3　計算結果對比Table 3　Comparison of results

計算結果見表4，5和6，由圖3，4，5，6分析表明：經過7次迭代，回路風量修正量Δqi滿足精度且密度、壓強、自然風壓和風量在運算過程中均收斂，證明本程序滿足通風網絡解算的設計要求，也證明了礦井氣候多參數預測與通風網絡解算算法的結合是可實行的。

表4　迭代驗算Table 4　Iterative calculation

表5　迭代驗算Table 5　Iterative calculation

表6　迭代驗算Table 6　Iterative calculation

圖3　各節點壓強迭代數值變化曲線Fig.3　Each node pressure iteration value curve

圖4　各分支風量迭代數值變化曲線Fig.4　Each branch wind volume iteration curve

圖5　各節點密度迭代數值變化曲線Fig.5　Each node density iteration value curve

圖6　各回路自然風壓迭代數值變化曲線Fig.6　Natural wind pressure of each circuit volume iteration value curve

4　結論

1)通過對風流的壓強和溫度變化規律的研究，在已知進風井入口風流的溫度和壓強，預測計算出礦井巷道各節點風流的溫度和壓強，從而得出各個巷道平均密度，實現了快速、簡單且不失準確度的計算礦井自然風壓目的。

2)結合Cross網絡解算方法，在已知井口初值的情況下，實現了自然風壓和機械壓雙重作用下的復雜礦井含自然風壓的通風網絡自動解算的過程。

3)通過建立礦井模型對本程序進行驗證，證明了在通風網絡解算中對自然風壓的考慮是有必要的。又經過實驗驗證，本程序滿足通風網絡解算的設計要求，也證明了礦井氣候多參數預測與通風網絡解算算法的結合是可實行的。

[1]王海寧. 礦井風流流動與控制[M]. 北京：冶金工業出版社, 2007.

[2]林海飛，劉靜波，嚴敏. CO2/ CH4在煤儲層中擴散規律的分子動力學模擬[J]. 中國安全生產科學技術, 2017，13(1)：84-89.

LIN Haifei, LIU Jingbo, YAN Min.The molecular dynamics simulation of the diffusion of CO2/CH4in the coal reservoir[J]. Journal of Safety Science and Technology, 2017，13(1)：84-89.

[3]馬恒,張越,周西華. 含自然風壓通風網絡的自動解算及應用[J]. 金屬礦山, 2016，32(1) ：157-161.

MA Heng, ZHANG Yue, ZHOU Xihua,. Automatic calculation and application of ventilation network with natural wind pressure[J]. Metal Mine, 2016， 32(1) ： 157-161.

[4]周西華,王繼仁,洪林. 濕空氣密度的快速準確測算方法[J]. 礦業安全與環保, 2005，32(4) ：49-51.

ZHOU Xihua, WANG Jiren, HONG Lin. Fast accurate measurement of wet air density[J]. Mining Safety ﹠ Environmental Protection, 2005, 32(4): 49-51.

[5]王海寧,汪光鑫,劉紅芳. 礦井自然風壓適時計算與應用[J]. 礦業研究與開發, 2014， 34(2)： 77-80 .

WANG Haining, WANG Guangxin, LIU Hongfang. Timely calculation and application of mine natural wind pressure[J]. Mining Research and Development, 2014, 34(2): 77-80.

[6]馬恒,尹彬,劉劍. 礦井風流溫度預測分析研究[J]. 中國安全科學學報, 2010, 20(11): 91-95.

MA Heng, YIN Bin, LIU Jian. research on prediction airflow temperature in Mine[J]. Chinese Safety Science Journal, 2010, 20(11): 91-95.

[7]約阿希姆.福斯. 礦井氣候[M]. 北京：煤炭工業出版社,1989.

[8]李恕和,王義章. 礦井通風網絡計算的牛頓法[J]. 煤炭學報, 1982, 10(4): 52-62.

LI Shuhe，WANG Yizhang. The Newton method for the mine ventilation network[J]. Journl of China Coal Society, 1982, 10(4): 52-62.

[9]劉劍,賈進章,鄭丹. 流體網絡理論[M]. 北京: 煤炭工業出版社, 2002.

[10]宮良偉,何華,鄒德均. 通風網絡解算中自然風壓處理[J]. 礦業安全與環保, 2015, 42(6): 113-116.

GONG Liangwei,HE Hua,ZHOU Dejun.Treatment of natural wind pressure in ventilation network flowcalculation[J]. Mining Safety ﹠Environmental Protection, 2015, 42(6): 113-116.

[11]陳明生，陳伯輝，沈斐敏. 礦井通風優化評價指標體系權重確定[J]. 中國安全生產科學技術, 2011, 7(3): 22-25.

CHEN Mingsheng,CHEN Bohui, SHEN Feimin.Weight determination of evaluation index system in mine ventilation optimization[J].Journal of Safety Science and Technology, 2011, 7(3): 22-25.

[12]鐘德云，王李管，畢林. 復雜礦井通風網絡解算風網有效性分析[J]. 中國安全生產科學技術, 2014, 10(11): 10-14.

ZHONG Deyun, WANG Liguan, BI Lin.Analysis on ventilation network effectiveness in ventilation network calculation of complex mine[J]. Journal of Safety Science and Technology, 2014, 10(11): 10-14.

[13]鄧云峰,蓋文妹.基于動態網絡流的區域疏散分析模型與仿真[J].中國安全生產科學技術, 2016, 12(11): 24-30.

DENG Yunfeng,GAI Wenmei.Analysis model and simulation of regional evacuation based on dynamic network flow[J].Journal of Safety Science and Technology, 2016, 12(11): 24-30.

[14]周西華,王繼仁,盧國斌.回采工作面溫度場分布規律的數值模擬[J]. 煤炭學報,2002,27(1): 59-63．

ZHOU Xihua，WANG Jiren，LU Guobin.The numerical simulation of distribution law of temperature field in coal face[J]. Journal of China Coal Society, 2002,27(1) :59-63．

[15]王樹剛，徐哲，張騰飛．礦井熱環境人體熱舒適性研究[J]. 煤炭學報, 2010, 35(1): 97-100．

WANG Shugang, XU Zhe， ZHANG Tengfei. Human thermal com-fort for mine environment[J]. Journal of China Coal Society，2010， 35(1): 97-100．