有限體積法在采空區自然發火數學模型中的應用

王月紅，蔣冀萍

(華北理工大學 礦業工程學院，河北 唐山 063210)

0 引言

采空區自燃是采空區經常發生的一種現象，尤其是發生在距采空區有一定間隔的氧化升溫帶中，采空區分為冷卻帶、氧化升溫帶和窒息帶，氧化升溫帶的后面便是窒息帶，窒息帶內不發生自燃，因此只選取部分區域對采空區自燃進行研究即可。為了簡化計算，根據采空區不斷變化移動的特點，該項研究在模型中引入了移動坐標系的概念，使得采空區自然發火模型由靜坐標下移動邊界的非穩態的偏微分方程組變成了移動坐標下邊界固定的穩態的偏微分方程組，在很大程度上簡化了采空區自然發火三維數學模型的求解。影響采空區自然發火的主要因素有采空區內的風流狀況、氧氣濃度和溫度，而其中的溫度又分為遺煤溫度和氣體溫度，因此該項研究分別對采空區內流場、氧濃度場、遺煤溫度場和氣體溫度場進行了分析和相關數學模型的推導，分別建立了采空區流場模型、氧濃度場模型和溫度場模型，通過3個模型的耦合求解來建立采空區自然發火模型。

1 采空區流場模型

基于有限體積法“守恒性”的思想，建立采空區流場模型。有限體積法強調從物理角度建立離散方程，在進行推導的過程中物理概念清晰，離散方程的系數均具有一定的物理意義，每個離散方程均為有限體積法中的某個物理量守恒的表達式，并且從物理角度建立方程可以保證該方程具有“守恒性”。有限體積法的優勢在于其一方面不僅具有有限元方法的“靈敏性”，而且適用于處理復雜的范圍以及邊界問題；另一方面，有限體積法采用顯示算法，具有更加靈活的假設，可以克服泰勒展開離散的缺點，并且可以很好地解決復雜的工程問題。

1.1 達西定律

達西定律是基于有限體積法思想和質量守恒方程，在采空區內任取一個空間六面體微元，在空間六面體微元中任取一點來考察采空區漏風流的三維流態。有如下關系：

(1)

式中：ux、uy、uz分別為坐標軸3個方向上的速度分量；Kx、Ky、Kz分別為3個方向上的滲透系數；?p/?x、?p/?Y、?p/?Z分別為每個方向上的壓力梯度。

1.2 流場控制方程

假設有一流場，取流場內的某一點，該點被如圖1所示的六面微元體包圍，坐標軸正向為已知。

圖1 采空區內的某一六面微元體

單位時間內風流流進采空區、流出采空區導致的質量差如式(2)所示。

(2)

單位時間內風流流入控制體與流出控制體的差值和微元體內質量的增加量應該是相等的，如式(3)所示。

(3)

1.3 邊界條件

采空區的流場區域及邊界如圖2所示，Γ1、Γ2、Γ3、Γ4為采空區靠近進風側的面，Γ5和Γ6分別為底板和頂板。

圖2 采空區的流場區域及邊界

采空區靠近工作面邊界的壓力值可以測定，因此可表示為如式(4)所示。

P(x,y,z)|Γ4=p0(x,y,z)|(x,y,z)∈Γ4

(4)

式中：Γ4為已知邊界；P0(x，y，z)為已知滲流邊界上的風壓函數。

則流場模型可表示為：

(5)

2 采空區氧濃度場模型

2.1 Fick定律

Fick定律描述擴散作用，可以使用這條定律來求得擴散系數，該定律是由德國科學家阿道夫·菲克于1855年推導出來的。

(6)

式中：J為擴散通量，mol/(m2·s)；D為擴散系數或擴散度，m2/s；?為濃度，mol/m3；x為位置長度，m。

2.2 采空區氧濃度場模型

對采空區內任取的某一六面微元體進行質量守恒分析，在單位時間內，微元體內部的氧氣質量變化量(WC)主要由以下3個方面構成：(1)采空區風流流入微元體和流出微元體的氧氣質量差值(W1)；(2)氧氣消耗量(W2)；(3)濃度差異引起的流入微元體和流出微元體的氧氣質量之差(W3)，即：

WC=W1+W2+W3

(7)

其中：

(8)

(9)

(10)

(11)

將式(8)、(9)、(10)以及(11)代入(12)得：

(12)

為了便于計算，將式(12)中質量濃度改為摩爾濃度，即為式(13)。

(13)

2.3 邊界條件

采空區的氧濃度場區域及邊界如圖3所示，共有Γ1、Γ2、Γ3、Γ4、Γ5、Γ66個面的邊界。

圖3 采空區的氧濃度場區域及邊界

在實際的處理中，只有在工作面附近的Γ4邊界可以直接測量，因此對于進風側可以直接考慮空氣中的氧濃度；對于回風側，考慮到消耗氧氣的因素，需要對氧濃度進行實際的測量，但在該項研究中，回風側處的邊界根據第二類邊界條件處理，也就是垂直于工作面方向的氧濃度變化是零，如式(14)所示。

(14)

采空區內進風側和回風側，即Γ1邊界和Γ3邊界的耗氧量幾乎可以不計，所以在該方向上的氧濃度變化率為零，因此可根據第二類邊界條件進行處理，對于采空區的深部邊界(Γ2邊界)和上下邊界(Γ5、Γ6邊界)，同理根據第二類邊界條件處理，如式(15)所示。

(15)

即采空區氧濃度場模型如(16)所示。

(16)

3 采空區溫度場模型

由于難以直接、準確地判斷采空區的著火點位置和自燃環境的火情，因此在工程實踐中，采空區的防火、滅火工作不僅相對盲目，而且還浪費了大量的人力物力。近幾年，隨著綜采放頂煤技術的推廣應用，大大提高了生產效率，但同時也導致采空區內遺煤量更多、漏風更加嚴重的結果，因此對采空區內溫度場進行立體模型公式推導時，將溫度場分為遺煤溫度場和氣體溫度場，為及時預測以及研究采空區自燃問題提供了有效直觀的研究方法。

3.1 傅里葉定律

傅里葉定律是在1822年被法國著名科學家傅里葉提出的一條熱力學定律。表達式如式(17)所示。

(17)

式中：JT為熱流密度，w/m2；k為熱導率，w/(m·k)；dT/dx為該方向上的溫度梯度。

3.2 溫度場模型

3.2.1 采空區遺煤溫度場模型

對采空區放頂煤巖中小型控制體的熱平衡進行分析，得到如式(18)所示的守恒方程。

Qsc+Qf-Qd=Es

(18)

式中：Qsc表示由于導熱導致的熱量差，kJ/(mol·s)；Qf表示由內熱源導致的熱量，kJ/(mol·s)；Qd表示由對流換熱導致的熱量差，kJ/(mol·s)；Es表示由工作面移動導致的進出控制體的熱量差，kJ/(mol·s)。

其中：

(19)

(20)

(21)

(22)

將式(19)、(20)、(21)、(22)代入式(18)中即可得到如式(23)所示的采空區內固體顆粒的熱平衡方程。

(23)

式中：ts表示冒落煤巖溫度，K；tg表示氣體溫度，K；λs表示固體導熱系數，W/m·K；Ke表示對流換熱系數，J/(m2·s·K)；q(t)表示單位時間單位體積的煤氧化放熱強度，kJ/(m3·s)；Cs表示固體比熱容，kJ/(kg·K)；v0表示工作面推進速度，m/d；ρs表示固體密度，kg/m3；Sn表示單元體內固體顆粒與氣體對流換熱的表面積，m2。

3.2.2 采空區氣體溫度場模型

同理，對采空區中的任一氣體控制體進行熱平衡分析，可得如式(24)所示的能量守恒方程。

Qgc-Qh+Qd=0

(24)

式中：Qd表示單位時間內由固體與孔隙中氣體導致的對流換熱量，kJ/(mol·s)；Qh表示單位時間內由氣體流動帶進和帶出導致的熱量差，kJ/(mol·s)；Qgc表示單位時間內由導熱流進和流出控制體導致的熱量差，kJ/(mol·s)。

其中：

(25)

(26)

(27)

將式(25)、(26)、(27)代入(24)得：

(28)

式中：λg表示氣體導熱系數，W/m·K；cg表示氣體比熱容，J/(kg·K)；ρg表示氣體密度，kg/m3。

3.3 邊界條件

采空區溫度場模型的邊界條件和采空區氧濃度場的邊界條件相似，如圖3所示，由于采空區溫度場的溫度可以直接測定，因此對于臨近工作面進風側的邊界條件根據第一類邊界條件進行處理，表達式如式(29)所示。

(29)

對于采空區的其它邊界，可認為在該邊界坐標軸方向上的溫度沒有發生變化，表達式如式(30)所示。

(30)

將偏微分方程、采空區遺煤溫度場邊界條件和氣體溫度場邊界條件進行整合，可得到采空區溫度場數學模型，表達式如式(31)所示。

(31)

4 采空區自然發火三維數學模型

4.1 采空區自然發火三維數學模型

將式(5)、(16)、(31)進行聯立得移動坐標下采空區自然發火三維數學模型：

(32)

4.2 邊界條件

采空區自然發火三維數學模型的邊界條件如式(33)所示。

(33)

將偏微分方程與采空區流場、氧濃度場、以及溫度場的邊界條件相結合，可以得到完整的采空區自然發火三維數學模型，該數學模型的優勢不僅在于從能量守恒的角度建立，具有明確的物理意義，而且還能能準確地描述采空區自燃的過程，為采空區自燃的防治提供了理論依據。

5 結論

(1)根據工作面不斷移動的特點，建立了移動坐標系，使得采空區自然發火模型由靜坐標下移動邊界的非穩態的偏微分方程組變成了移動坐標下邊界固定的穩態的偏微分方程組，在很大程度上簡化了采空區自然發火三維數學模型的求解。

(2)基于有限體積法“守恒性”思想，利用氣體流動規律和達西定律，建立了采空區流場模型，描述了采空區內部壓力和風流速度的分布情況；根據氧氣擴散規律和質量守恒定律，建立了采空區氧濃度場模型，描述了采空區內部氧濃度的分布情況；根據熱傳導規律和傅里葉定律，建立了采空區氣體溫度場模型和采空區遺煤溫度場模型，描述了采空區內部溫度的分布情況，為采空區自然發火研究提供了理論依據。

(3)采空區自然發火模型是由采空區流場模型、氧濃度場模型、以及溫度場模型進行耦合求解來建立的，采空區的氣體壓力、氧氣濃度和溫度分布情況在不同程度上對采空區自然發火都具有一定的影響，采空區自然發火模型的建立為采空區內不同參數影響的判斷提供了理論依據。