礦井采空區漏風問題的迎風有限元求解技術及其應用

吳奉亮，李智勝，常心坦

（1.西安科技大學安全科學與工程學院，陜西 西安710054；2.西安科技大學西部礦井開采與災害防治教育重點實驗室，陜西西安710054）

0 引 言

煤礦采空區漏風容易引發煤自燃，形成可爆瓦斯區域［1］甚至引發瓦斯爆炸。由于對采空區進行直接觀測十分困難，數值模擬成為國內外學者研究采空區漏風問題的有效手段，如判定采空區內的煤自燃三帶［2-3］，模擬煤自燃過程［4-5］及瓦斯運移規律［6-7］。目前模擬采空區漏風問題的途徑主要有采用商業CFD與自編程2種，相關的數值方法包括有限元、有限體積法等。盡管商業CFD（如Fluent，COMSOL等）已被大量學者使用，但它對理論知識要求高，目前還很難被現場人員所接受。自編程研究采空區漏風問題很可能是將CFD技術推向現場應用的有效途徑，因為通過自編程可將礦井通風問題的相關概念作為軟件的界面元素進行設計，將CFD算法封裝到后臺［8］，以此來降低CFD技術的使用門檻。有限元法相對其它數值方法具有對網格劃分要求低、易編程的特點，許多偏微分方程都可通過標準的伽遼金法拆分成代數方程組，這些特點使其成為自編程研究采空區漏風問題的主要數值方法：章夢濤、趙陽升等給出了采空區線性滲流模型的有限元求解方法［9-10］；丁廣嚷、李宗翔等研究了采空區非線性滲流的有限元求解模型［11-12］，基于Matlab開發了可以模擬采空區瓦斯涌出與耗氧情況［13］的G3程序。采空區漏風引起的瓦斯運移與遺煤耗氧驅動的氧濃度場變化都屬于“對流-擴散”過程。在對流起主導作用的采空區流體輸運過程中，基于標準伽遼金法導出的有限元控制方程給出的是一個不正確的振蕩解，Zienkiewicz，章本照等的研究表明這個問題可以通過迎風有限元法來解決［14-15］。盡管采空區中瓦斯與氧氣的運移過程都是“對流-擴散”問題，但在數學模型描述中具有不同的源項表現，采用普通對流擴散模型的控制方程進行求解迭代量大，并不高效。文中基于迎風有限元法思想，導出采空區“瓦斯-氧”對流擴散的統一模型，通過自編程將其用于煤自燃三帶與可爆瓦斯濃度覆蓋區域的判定。

1 采空區滲流特性參數分布規律

采空區流場的滲透系數E，孔隙率ε及冒落高度H等參數主要受采空區的冒落煤巖碎脹系數影響，研究表明碎脹系數呈“O”形圈分布，以圖1所示U形工作面通風系統為例，各參數可用以下公式計算［16］

圖1 U形通風工作面采空區布置圖Fig.1 Gob layout of U-shape working face ventilation system

式中 Kp（x，y）為給定點（x，y）的碎脹系數；Kp，max和Kp，min分別為初始和壓實后的碎脹系數；x0為開切眼到Y軸的距離；u為回采速度；t為回采時間；a0和a1為傾斜與走向方向冒落煤巖碎脹衰減系數；d0=y-｜y-0. 5W｜，其中 W為采空區寬度；ξ為“O”形圈的幾何調整系數；h為采高；υ為空氣運動粘性系數；g為重力加速度；d為粒子平均直徑。

2 采空區漏風問題的定解模型

2. 1 采空區流場非線性滲流定解模型

采空區中的風流在冒落煤巖體中的流動是各向同性非均勻介質滲流，由于采空區寬度與深度遠大于高度，因此三維采空區流場完全可簡化為如圖1所示的平面流場。若先假定其流動符合達西定律，則流速V，滲透系數E和壓力P之間的關系為

進一步考慮連續性方程，得到平面采空區流場的定解模型為（6）式微分方程［9-10］

式中 G為整個采空區區域，采空區與風網連接處的邊界L1為壓力邊界，p0為L1上的已知壓力值；L2為不漏風的邊界；nx，ny為邊界外法線上單位向量的分量。實際上采空區內，特別是鄰近工作面的范圍，一般為非線性的滲流，其流動符合Bachmat方程［17］

式中 β為粒子形狀系數。采用對滲透系數進行迭代修正的方法，可基于（6）式的線性模型求解以上非線性滲流，詳細的迭代原理及過程見文獻［12］。

2. 2 采空區中“瓦斯-氧”運移的定解模型

2. 2. 1 采空區瓦斯涌出與耗氧強度

采空區中單位厚度上的瓦斯涌出強度（mol/（m2h））計算式為［16］

式中 w′0和w0分別為采空區局部地點恒定的瓦斯涌出強度和下部煤層沿底板裂隙持續涌入瓦斯的強度，mol/（m2h）；w1和 w2是采空區遺煤和鄰近煤層瓦斯初始涌出強度，mol/（m2h）；λi是瓦斯涌出量的衰減系數，d-1.根據 Arrhenius公式［4］，煤的耗氧速率WO2可用下式描述

式中 A為指前系數；CO2為氧氣摩爾濃度，mol·m-3；Ea為表觀活化能，kJ/mol；R是氣體常數，J/（mol·Ｋ）；T絕對溫度，Ｋ.由于式中exp（-Ea／RT）項在煤自燃的早期階段變化很小，因此上式可簡化為

即煤的耗氧速度與氧濃度成線性關系（A′是線性系數），與文獻［2］中通過實驗得到的結論相同。考慮遺煤厚度H1的影響，則采空區中煤的耗氧速率應為

2. 2. 2 采空區中“瓦斯-氧”運移的數學模型

由于工作面回采速度遠小于采空區內氣體對流擴散的速度，因此采空區中的瓦斯與氧氣運移可視為穩態“對流-擴散”過程［18］，描述這一過程的偏微分方程及定解邊界條件為

Θ為 CH4或 O2中的一種；CΘ為氣體濃度；V′Θ為采空區氣體流速，V′Θ=V／ε；DΘ為氣體擴散系數；IΘ為源項中與CΘ無關的部分，sΘ為源項與CΘ存在線性關系的系數，IΘ與 sΘ的取值分別見表1；ГC為采空區的入風邊界，Гq表示回風邊界及不漏風的邊界；D為氣體擴散系數張量。D的計算式為［9］

式中 i，j取值 1或 2，分別為 x，y方向；D′為分子擴散系數，m2/s；τ為介質曲率；當 i=j時，δij=1，否則δij=0；aT，aL分別為縱向與橫向彌散度。

表1 公式（12）中IΘ和 sΘ的不同表現形式Table 1 Different forms of IΘand sΘ in Eq. （12）

3 “瓦斯-氧”對流擴散迎風有限元解算法

3. 1 用伽遼金法求采空區流場

若將待求采空區問題（如（6）式中的壓力P或（12）式中的濃度 CΘ）的真解記為 U（x，y），則它的有限元解是指將采空區離散成單元（文中采用三角形單元，設節點總數為n）后，用每個單元上的插值函數 UΔ（x，y）來近似 U（x，y），這個插值函數為

式中 ui，uj，uk為三角形單元的3個節點（按逆時針方向編號為 i，j，k）對應問題的待求量；Nη（η=i，j，k）是基函數，記 AΔ為單元面積，Ni=（ai+bix+ciy）／（2AΔ），Nj=（aj+bjx+cjy）／（2AΔ），Nk=（ak+bkx+cky）／（2AΔ），ai=xjyk-xkyj，bi=yj-yk，ci=xk-xj，aj=xkyi-xiyk，bj=yk-yi，cj=xi-xk，ak=xiyj-xjyi，bk=yi-yj，ck=xj-xi. 因此有限元方法的關鍵步驟是導出關于待求問題 ui（i=1，2，…，n）為未知數的控制方程組。加權余量法是最方便實施的一種方法，當在單元內為每個節點取定對應的權函數 ωη（η=i，j，k）時，對（6）或（12）式對應的問題導出的控制方程組的矩陣表示形式均可記為

伽遼金法是將權函數取為與基函數相同的一種加權余量法，對（16）式取ωη=Nη，積分得到［0，0，0］T，單元矩陣為一對稱陣，其計算式為

3. 2 采空區“瓦斯 -氧”運移的迎風有限元法

在（18）中取權函數ωη為基函數Nη（即采用伽遼金法）時，導出方程組的解是振蕩的［19］，解決辦法是對權函數進行修正以增大來流方向的效應，修正式為［14］

式中 hΔ為單元的特征長度，可取三角形單元的最長邊；α為由佩克萊特數Pe決定的系數，它們的計算式分別為

將（19）式代入（18）式積分后得到

式中 Sij項是源項中與濃度存在線型關系的項的積分結果，它的計算式為

綜上，在求得采空區流場后，將（23）式代入（15）式即可求得采空區氣體的濃度場。

4 應用實例

4. 1 模型實現與算例概況

文中基于 VC++與 ObjectARX［20］技術將以上模型及前后處理過程集成于AutoCAD中，方程組的求解采用PARDISO［21-22］并行求解器。應用實例為某礦U形工作面通風系統采空區，寬180 m，煤層開采厚度6m，工作面設計風量為25m3/s，入風隅角與回風流隅角的相對壓力分別為308.1和268.2 Pa；工作面推進度為3 m/d，圖2為回采90 d后形成的270 m長采空區（按滲透系數渲染）。模擬中涉及的其它參數取值見表2.

4. 2 計算結果

圖3采用流網顯示采空區流場的模擬結果，其中流線的流函數值間距為1.9×10-3，等壓線間距為0.8 Pa，流網顏色按流速大小渲染，綠、黃、藍分別表示按流速劃定的采空區自燃三帶：散熱帶（流速大于 0.004 m/s），氧化升溫帶（流速位于0.001 7至0.004 m/s）和窒息帶（流速小于0.001 7 m/s）。圖4（a）為模擬得到的采空區的氧濃度，采用氧濃度劃分三帶時，通常將濃度介于8%至18%的區域劃定了氧化升溫帶，如圖中黃色區域。由于風速過大將使熱量無法積聚，因此更合理的方法是用風速與氧濃度雙指標進行劃分，即氧化升溫帶應滿足

圖2 采空區滲透系數Fig.2 Distribution of hydraulic conductivity in gob area

表2 模擬參數Table 2 Parameters for the numerical simulations

式中 CO2-min=8%，Vmax=0.004 m/s.按（25）式將圖3與圖4（a）疊加得到圖4（b）所示的三帶（仍用綠、黃、藍區分）劃分結果。模擬結果中三帶的分布及自燃帶的形狀與文獻［23］采用FLUENT模擬的結果相似：采空區自燃帶位于工作面后一定位置，其寬度沿進風側向回風側逐漸減小。

圖3 采空區流網Fig.3 Flow net in gob area

圖4 采空區氧濃度場及自燃三帶Fig.4 Distribution of oxygen concentration and“three zones”for analyzing spontaneous combustion in gob

圖5為模擬得到的采空區瓦斯濃度場。藍色顯示采空區深部是瓦斯富積的區域，濃度大于15%，綠色是濃度低于5%的區域，黃色區域瓦斯濃度介于5%～15%、在正常空氣氧濃度下具有爆炸危險。從圖3，圖5可知漏風主要從2個隅角流入采空區與流回工作面，工作面與采空區邊界內側、近上隅角的49.5 m長范圍內（回風隅角）瓦斯濃度大于5%（5%～23%）。由于遺煤耗氧，使得采空區中氧濃度降低，不一定在任何地方都滿足瓦斯爆炸的需要。圖6左圖的瓦斯爆炸三角形給出了瓦斯與氧氣形成可爆氣體的濃度范圍［24］，當氧濃度低于12%時，將不滿足瓦斯爆炸的需要。據此將圖4（a）與圖5疊加，得到圖6右圖紅色區域所示的具有瓦斯爆炸危階性區域，其中綠色、黃色分別表示瓦斯含量過低和過高的區域。因此在圖6右圖所示的紅色區域內，任何點火源（包括煤自燃）均可引發瓦斯爆炸。因此針對模擬結果，對圖中標定的危險區域進行惰化，以及保證必要的推進度來防止煤自燃對避免采空區火災或爆炸事故的發生是十分必要的。

圖5 采空區瓦斯濃度場Fig.5 Field ofmethane concentration in gob area

圖6 瓦斯爆炸三角形（左）及采空區具有瓦斯爆炸危階的區域（右）Fig.6 Methane explosibility triangle and presumed location of the explosive gas zones in the gob

5 結 論

1）基于伽遼金法建立了采空區流場的求解模型。給出了工作面動態回采下采空區滲流參數的計算方法。針對采空區流場的定解模型，基于三角形單元離散的采空區，采用標準伽遼金法導出了以所有節點壓力為未知數的有限元控制方程組，為自編程研究采空區漏風問題提供了基礎；

2）基于迎風有限元法建立了采空區“瓦斯-氧”輸運過程的統一模型。分析了采空區瓦斯涌出與遺煤耗氧的特點，建立了采空區中“瓦斯-氧”輸運過程的對流擴散模型。針對標準伽遼金法求解對流擴散模型時解的振蕩問題，給出了迎風有限元法求解“瓦斯-氧”輸運過程的控制方程組；

3）開發了采空區漏風問題模擬軟件。基于以上模型，以及VC++，ObjectARX等技術開發了可視化模擬軟件，應用實例表明使用此軟件可方便地進行采空區模型的前后處理，并根據模擬結果給出采空區煤自燃三帶及具有瓦斯爆炸危險性的區域。