李 坤 林 琦 楊紅陽

（安徽工業(yè)大學，安徽 馬鞍山 243002）

自燃是自然界中較為常見的災害現(xiàn)象，給人類的生產(chǎn)、生活、健康，乃至大氣環(huán)境等方面均帶來了極大的危害[1-2]。僅以煤炭行業(yè)為例，我國近十年發(fā)生的195 起煤礦火災事故中九成以上是由自燃引起的，造成了共計200 余人死亡，損毀的煤炭資源高達千萬噸[3-4]。

自燃具有潛伏期長、隱蔽性強、熄滅困難等特點，在潛伏期內(nèi)進行人工干預的效果明顯優(yōu)于自燃發(fā)生后的滅火治理，故對潛伏期內(nèi)自熱演化過程的準確掌握尤為重要[5]，特別是谷倉、煤堆等大型堆積態(tài)系統(tǒng)?，F(xiàn)階段，大型系統(tǒng)自燃規(guī)律預測的相關研究工作仍停留在宏觀層面，多對系統(tǒng)內(nèi)部結構采用均質(zhì)假設。自燃系統(tǒng)往往是典型的非均質(zhì)多孔結構，其內(nèi)部的空隙變化會直接影響氣流組織分布，進而改變了系統(tǒng)內(nèi)部各區(qū)域的生熱和蓄熱能力。因此，探究堆積態(tài)系統(tǒng)內(nèi)部的空隙率分布規(guī)律是自燃現(xiàn)象準確預測的重要前提條件。

當前，工程尺度下的堆體內(nèi)部非均質(zhì)空隙率分布測試仍是一個難題。雖然CT 掃描等可以在不損壞堆體的前提下獲得三維堆體的真實結構[6-8]，但其樣本范圍多限于厘米量級。而對于大型的堆體，顆粒的松散滑動使得空隙率數(shù)據(jù)較難完好無誤地提取。堆積態(tài)系統(tǒng)是通過顆粒從高空傾倒、下落形成，若能通過模擬仿真方法來重現(xiàn)物料顆粒的真實下落堆積過程，就可獲得堆體內(nèi)部結構，實現(xiàn)堆體幾何結構的數(shù)值重構[9]，進而為堆體內(nèi)部空隙率分布特性分析提供數(shù)據(jù)支撐。離散元法恰好可以模擬大量顆粒碰撞、遷移，追蹤顆粒運動，確定下落顆粒實時位置，實現(xiàn)堆體的數(shù)值重構。

該文基于離散元方法模擬真實煤粒下落堆積過程，研究不同當量粒徑和下落高度工況下的堆體非均質(zhì)空隙分布，并將堆體空隙分布變化規(guī)律和與空隙率相關的變量編成自定義函數(shù)進行非均質(zhì)煤自熱模擬，使得堆體自燃研究從“黑箱”均質(zhì)假設擴展到“灰箱”非均質(zhì)。

1 堆積計算模型

1.1 物理模型

工程中的物料堆積多通過傳送帶送至特定高度，從落料口自由下落并堆積。本研究為簡化過程，采用顆粒從漏斗落下過程來近似物料送出落料口的過程。在傾倒初期，總質(zhì)量為10 kg 顆粒從漏斗頸部小口流出，進行自由落體運動，直至接觸到底部圓盤，發(fā)生碰撞、反彈，進而緩慢堆積，最后積疊成一個小堆體。為保證模擬結果更貼近現(xiàn)實，通過查詢Rocky 幫助手冊確定煤粒的形狀參數(shù)，其垂直縱橫比為1.25，水平縱橫比為0.7，平滑度為1，轉(zhuǎn)角數(shù)為11，超二次曲面角為6.2。煤體與平臺的物性參數(shù)[10]中漏斗與地面所用的鋼材剪切模量為8.00×1010Pa，煤的剪切模擬量為4.70×108Pa，鋼與煤的泊松比均為0.3，煤與煤的恢復系數(shù)為0.64，煤與鋼的恢復系數(shù)為0.65。

1.2 數(shù)學模型

1.2.1 假設

為使問題簡化以便數(shù)學建模，作如下假設：

1）忽略顆粒之間的黏性作用；

2）忽略顆粒的破碎作用；

3）煤粒非剛體，接觸過程中允許出現(xiàn)微小變形，以彈簧模型來近似反映顆粒的擠壓。

1.2.2 運動控制方程

式中：mi為顆粒i的質(zhì)量，kg；Vi為顆粒i的速度，m/s；t為時間，s；ki為與顆粒i接觸的顆粒數(shù)；Fij為顆粒i與j的接觸力，N；g為重力加速度，m/s2。Ii為顆粒i的轉(zhuǎn)動慣量，kg·m2；wi為顆粒i的角速度，rad/s；Tij為顆粒j作用于顆粒i的力矩，N·m。

1.2.3 碰撞模型

顆粒相互碰撞，其間的接觸力可分解為法向與切向接觸力。

式中：Fn為顆粒i與j的接觸力的法向分量，即法向接觸力；Fτ為顆粒i與j的接觸力的切向分量，即切向接觸力。顆粒間的法向接觸力可用滯回線性彈簧模型描述，而顆粒間的切向接觸力可用線性彈簧庫侖極限模型描述。

1.3 模型驗證

為了驗證模型的準確性，本文模擬分析了實驗臺尺度下的煤粒下落堆積過程，并與實驗結果[10]進行對比分析，發(fā)現(xiàn)煤體顆粒積疊形態(tài)的模擬和實驗堆積角分別為28.27°和28.51°，僅為0.84%的相對誤差，模擬與實驗結果的堆積輪廓線基本吻合，呈現(xiàn)出較好的一致性。

2 空隙分布

收集分析1.5 cm 當量粒徑顆粒從0.55、0.65、0.75、0.85、0.95 m 的高度下落以及1.1、1.3、1.5、1.7、1.9 cm 當量粒徑顆粒從0.65 m 高度下落所形成的堆體空隙率數(shù)據(jù)。結果發(fā)現(xiàn)，堆體的垂直方向上的空隙率隨垂直高度的增大而先減小后增大，堆體水平方向上的空隙率基本保持不變，而下落高度與顆粒粒徑對堆體的空隙分布影響較小。

3 煤自熱模擬

煤堆水平方向空隙基本恒定，隨垂直高度的增大而先減小后增大，則煤堆整個區(qū)域可用垂直方向先減小后增大的變化趨勢規(guī)律描述，應用此規(guī)律研究非均質(zhì)煤堆自熱過程。

3.1 幾何模型及簡化假設

考慮實際煤場中煤體堆積要防止滑坡以及再運輸方便，多將其堆積為橫截面為等腰梯形的大壩體。如圖1（a）所示，該模型為二維，由兩個域組成：

圖1 氣流作用下的煤堆（m）

1）多孔域煤堆為等腰梯形，底邊長30 m，高10 m，底角為45°；

2）剩余區(qū)域代表煤堆周圍開放空間的流場，長為90 m，寬為30 m。

將該模型進行網(wǎng)格劃分，如圖1（b）所示，采用全局尺寸為1 m 的四邊形網(wǎng)格填充外部流場，局部尺寸為0.5 m 的三角形網(wǎng)格填充多孔域。煤堆與流場交界區(qū)域質(zhì)量、熱量交換最為頻繁，是自熱首先發(fā)生的地方，故在此處建立5 層厚度均為0.1 m的附面層薄層。

煤的自燃過程涉及多種反應，因此很難建立精確的數(shù)值模型，故必須進行一些簡化[11]：

1）煤堆物性參數(shù)不隨自熱過程發(fā)生改變；

2）采用一步全局氧化反應來研究煤自燃過程；

3）忽略煤堆中影響較小熱輻射效應和水分的作用；

4）流場與煤堆是一個整體，質(zhì)量、熱量交換均在其內(nèi)部。

3.2 控制方程

煤一步全局氧化反應可簡單表示為[12]

其 中：vO2、vCO2、vCO和（2vO2-2vCO2-vCO）代 表O2、CO2、CO 和煤氧復雜固體產(chǎn)物的化學計量系數(shù)。在30～170 ℃范圍內(nèi)，根據(jù)不同的反應機理將煤氧化反應分為三個階段：30～60 ℃為第一階段，60～100 ℃為第二階段，100～170 ℃為第三階段。具體參數(shù)設置見表1。

表1 數(shù)值模擬中煤氧化的設置

3.2.1 能量方程

煤在低溫下的氧化速率極為緩慢，將其產(chǎn)生的熱量定為源項，寫成用戶定義函數(shù)（UDF）。固體煤顆粒和氣流之間的溫差對煤堆內(nèi)部的熱相互作用影響頗大，需要將各自儲存的能量以及二者的熱能交換量分別表示出來。

1）固體能量守恒

式中：n為煤堆內(nèi)部的空隙率；ρc為煤粒的密度，kg/m3；Cpc為煤粒的定壓比熱容，J/（kg·K）；Tc為煤體的溫度，K；t為時間，s；λc為煤的導熱系數(shù)，W/（m·K）；d為煤粒的直徑，m；h為煤堆的對流傳熱系數(shù)，W/（m2·K）；Tg為煤堆內(nèi)部氣體的溫度，K；r為耗氧速率，kmol/（m3·s），ΔH為消耗單位摩爾氧氣產(chǎn)生的氧化熱，kJ/mol。方程左邊表示煤顆粒的內(nèi)能增長，右邊第一項為固體煤中的熱擴散，第二項為內(nèi)部氣體相互作用的熱對流，第三項為煤氧化產(chǎn)生的熱量源項。

2）氣體能量守恒

式中：ρg為煤堆內(nèi)部氣體的密度，kg/m3；Cpg為內(nèi)部氣體的定壓比熱容，J/（kg·K）；v為氣流速度矢量；P為壓力，Pa；λg為內(nèi)部氣體的導熱系數(shù)，W/（m·K）。方程左邊第一項表示氣流的瞬態(tài)能量上升，第二項表示氣流的熱對流；方程右邊第一項表示氣流中的熱擴散，第二項表示與固體煤相互作用的熱對流。

3.2.2 組分守恒方程

在整個自熱過程中，氮氣既不消耗也不產(chǎn)生，因此主要研究反應物氧氣以及生成物二氧化碳、一氧化碳的組分守恒。

氧氣：

式中：Co、Ci、Cj分別表示氧氣、二氧化碳、一氧化碳的氣體濃度，kmol/m3；Do、Di、Dj分別表示氣體混合物中氧氣、二氧化碳、一氧化碳的質(zhì)量擴散系數(shù)，kg·m2/s；DT,o、DT,i、DT,j分別表示氣體混合物氧氣、二氧化碳、一氧化碳的熱擴散系數(shù)，kg·m2/s。

3.2.3 多孔介質(zhì)中的動量平衡

加入粘性項對N-S 方程進行修正，求解多孔介質(zhì)煤堆的動量平衡，得到：

其中多孔介質(zhì)中兩個較為重要的力為粘性力與慣性力，其表達式：

式中：τ為黏滯應力張量，N/m2；Fb為體積力矢量，N/m3；μ為氣體粘性，kg·m/s2；k為煤基質(zhì)的滲透率，m2；C1為粘性阻力，m-2；C2為慣性阻力，m-1。

綜上可知，空隙率是眾多自熱因素的因變量，空隙率變化將會使得這些參數(shù)相應地發(fā)生改變。因此，必須將空隙率變化規(guī)律以及與空隙率相關的多孔介質(zhì)煤堆內(nèi)的粘性阻力、慣性阻力、氧氣的消耗速率、二氧化碳和一氧化碳的生成速率等參量編譯成用戶自定義函數(shù)并進行代碼調(diào)試，直至確定無誤后，方可進行非均質(zhì)煤堆的自熱模擬計算。

3.3 煤自熱模擬結果

將參數(shù)和求解的邊界條件設置完畢，采用耦合算法、二階迎風格式進行瞬態(tài)計算。在保持計算穩(wěn)定性與收斂性的前提下盡可能地減少計算時間，將計算時間步長定為1 h，監(jiān)測煤堆在熱失控前的溫度場、壓力場以及組分場的變化。

如圖2 所示，煤體10 d 升高3.5℃， 30 d 升高11℃， 60 d 升高26.5℃，90 d 升高100℃，隨時間的推移，升溫的速度極速加快，逐漸接近熱失控。溫度場的整體溫度隨時間的推移逐漸升高，溫度較為明顯地分成三個區(qū)域：煤堆的邊界處由于與氧氣接觸最為頻繁、反應產(chǎn)生熱量最多而溫度最高；再者是煤堆的中部區(qū)域；核心半圓處由于氧氣最少、反應產(chǎn)熱不足而溫度最低。隨著時長的增加，高溫區(qū)域逐漸從與風流接觸的外邊界縮小至迎風側(cè)、背風側(cè)，再慢慢收縮成靠近兩側(cè)與底部的兩個高溫點而成為煤堆安全的最大威脅。

如圖3 所示，煤堆內(nèi)部的壓力場主要可分為兩個區(qū)域，靠近邊界的壓力加大，內(nèi)部的壓力較小。隨著時間的推移，低壓區(qū)域逐漸減小，由近似梯形變?yōu)榘霗E圓形，半橢圓的高度也是逐漸減小，到最后幾乎貼著底部邊界。

圖3 煤堆流場隨時間的壓力變化圖

關于煤堆內(nèi)的各類別氣體濃度變化，其氧氣場與溫度場的變化相近，中心有低氧區(qū)域較低溫區(qū)偏大，二氧化碳與一氧化碳的濃度變化則與氧氣濃度場變化相反。這是因為氧氣要想通過重重煤顆粒的阻擋進入到中心區(qū)域是頗為困難的，進入到中心區(qū)域的氧氣本就很少并快速與周邊的煤粒反應生成二氧化碳、一氧化碳以及煤氧復合固體；相對而言，熱量傳導就容易得多，中心外部發(fā)生反應會將產(chǎn)生的熱量往中心傳導，致使低氧區(qū)中的低溫區(qū)面積減小。

4 結論

為了消除現(xiàn)有堆體自燃預測中均質(zhì)假設對結果準確性的影響，該文基于離散元方法模擬研究了大量顆粒從高空下落、堆積形成堆體的過程，獲取了空隙分布的基本信息，并借此研究非均質(zhì)煤的自熱過程。分析概括系統(tǒng)內(nèi)非均質(zhì)性分布及其煤自熱情況，得如下結論：

1）堆體內(nèi)的空隙隨垂直高度增加而先減小后增大，在水平方向上近似均勻。

2）研究范圍內(nèi)的傾倒高度大小變化以及粒徑粗細變化對堆體內(nèi)部空隙率影響較小。

3）非均質(zhì)煤自燃中高溫區(qū)域逐漸從與風流接觸的整個外邊界縮小至煤堆的迎風側(cè)、背風側(cè)，再慢慢收縮成兩個靠近迎風側(cè)、背風側(cè)與底部接觸的較其他區(qū)域溫度高出很多的點，這兩個“高溫點”將成為煤堆安全的最大威脅。