基于FLUENT 動網格的煤礦采空區三帶動態特性模擬技術

王芳

（西安科技大學，陜西 西安710000）

采空區自燃三帶分布情況與工作面的推進有著密切聯系。工作面推進過程中，采空區邊界的移動不僅使采空區實體的幾何尺寸發生了連續的動態變化，同時也影響了采空區內部的漏風供氧條件、各組分氣體傳輸和三帶分布情況。工作面的推進首先造成采空區幾何邊界的移動，其次由于回采工作面的兩端壓差提供漏風動力[1]，采空區漏風邊界也隨之變化并直接改變了采空區氧氣濃度的傳輸邊界和分布范圍[2]，體現為采空區流場、氧濃度場等多場的動態變化[3]，因此對采空區自燃三帶的數值模擬應該在工作面動態推進情況下進行模擬[4]。

1 采空區自燃三帶數值模擬

由于采空區是一個封閉的區域，基本上無法進入，對采空區各場的變化規律進行數值模擬是研究采空區煤自燃危險的重要方式之一。結合采空區的實際情況，建立數學模型，并利用計算機對數學模型進行解算，分析研究采空區自燃“三帶”的分布規律是是數值模擬研究采空區的一般步驟。其中，CFD（Computational Fluid Dynamics，計算流體動力學）的方法是先建立一套描述采空區風流流場、溫度場及氧濃度場的數學模型，以現場實測數據為已知條件，選擇合理的離散方法和計算程序，通過模擬采空區漏風流場、溫度場和氧濃度場，量化分析采空區煤自燃在受到不同邊界條件影響下的分布規律，最后對計算結果進行相應的處理并顯示輸出[5]。CFD 中的FLUENT 作為求解器，采用了多種求解方法和多重網格加速收斂技術，有較高的收斂速度和求解精度。使用FLUENT 進行求解前，應先建立模型、劃分網格，根據模擬要求選擇穩態或非穩態模式，選擇需要解的方程（流動方程、能量方程、組分運輸方程等），確定需要附加的模型（多孔介質模型），在對邊界條件及材料物理性質進行設置后，對模型進行求解。對工作面動態回采情況下采空區進行數值模擬，是通過模擬開采進行一段時間內的采空區內各場變化情況實現的，各場的變化與時間相關，可以采用非穩態模式進行計算。工作面動態采情況下的采空區范圍隨著工作面的移動而不斷擴大，采用FLUENT 中的動網格模型來實現采空區工作面動態回采，并且建立孔隙率、煤氧反應速率等與時間相關的函數通過用戶自定義功能編入FLUENT 進行計算。

2 FLUENT 動網格技術

采空區范圍隨工作面推進而不斷擴大的過程是計算區域不斷擴大的過程，在Fluent 中解決該問題的方法通常為采用動網格模型。Fluent 中的動網格模型可用來模擬由于流體域邊界剛性運動或者邊界變形引起的流體域形狀隨時間變化的問題，每一個時間步上的體網格的更新是由解算器根據邊界的新的位置確定的，根據邊界的運動自動調節內部體網格的分布，通過指定運動區域的運動方式以及網格再生的方法可以使網格動態變化。

對于邊界移動的任意控制體積V 上的一般標量的 φ的守恒方程可寫為：

式中，V 為空間中大小和形狀都隨時間變化的控制體積，?V為控制體積的運動邊界，ug為運動網格的運動速度；ρ為流體密度；u 為流體速度；г 為耗散系數；S φ 是標量 φ的源項。

在FLUENT 中，每一個時間步上體網格的更新是由解算器根據邊界的新的位置來自動完成的，即解算器可以根據邊界的運動和變形自動地調節內部體網格節點的分布，使用起來非常地方便。在使用動態網格時，用戶僅需提供初始網格并在模型中指定任意區域的運動即可。關于運動的指定，FLUENT 允許用戶通過邊界profile 文件或者用戶自定義函數(UDF)或者六自由度( 6DOF )解算器來指定。當邊界運動函數較復雜時，一般利用UDF 宏指定運動，然而當需要指定的邊界為簡單的速度- 時間關系運動時，利用Profile 文件指定運動則更為方便高效。

3 采空區自燃三帶數值模擬

3.1 模擬過程

以某煤礦為例，對煤礦原型進行簡化，建立模型，如圖1 所示。配風量為36 m3/s，進回風巷寬4m，高4m，工作面斜長200 m，平均推進速度為4m/d。模擬研究時空范圍是工作面從已回采了200 m 的位置開始之后的45 天回采期，到時采空區總長度為380m。

圖1 采空區模型

所建模型含有兩個部分，工作面巷道和采空區。因此劃分網格時選擇Multizone 網格劃分方法進行劃分，Multizone 可以對目標區域進行自動分區，將目標區域自動分解成多個可以掃掠或是自由劃分的區域，只需要簡單的指定源面、設置網格控制參數等，即可對目標區域進行得到高質量的網格。相較于傳統的分割、掃掠生成六面體單元的網格劃分方法，Multizone 網格劃分方法省略了分割步驟，只需要進行適當參數設置即可生成高質量的六面體網格。

網格劃分質量的高低直接影響了FLUENT 結果的精度以及收斂的速度，判斷網格質量的高低主要依據兩個指標：Element quality 和Skewness。Element quality 基于給定單元的體積與邊長之間的比率，其值處于0 到1 之間，越接近于1 越好；Skewness是最基本的網格質量檢查項，與Element quality 相反的是，越接近于0 表示網格質量越高。網格完成后的模型如圖2 所示，共劃分100512 個六面體網格，Element quality 值為0.99，Skewness 值為0.00000004，由此可見網格質量較好，可以滿足計算精度的要求。

圖2 采空區模型網格劃分

工作面推進的過程即可視為簡單的速度- 時間運動，工作面的移架、放頂雖然不是連續進行的，但是從一個很長的時間周期角度來說，開采過程中每日進刀數是較為固定的，因此可以將工作面相對于采空區的運動看作為勻速運動，可以通過Profile 文件來描述。

如回采速度為4m/d，即0.000046m/s，Profile 文件如下：

((v_y transient 4 0)

(time 0 1 2 4)

(v_y -0.000046 -0.000046 -0.000046 -0.000046))

網格再生方法主要用來計算內部網格節點的調節，采用鋪層（Layering）與局部重構（Remeshing）相結合的方法進行計算。動態模型的模擬采用瞬態模式，根據所需要模擬的時長來設置時間步長以及時間步數。

時間步長大致由以下公式估計：

3.2 數值模擬結果

3.2.1 穩態結果

在不使用動網格方法的情況下進行穩態模擬，得到采空區氧氣濃度場的分布情況如圖3。按照氧氣體積濃度指標將模擬所得氧氣濃度場劃分成為的自燃“三帶”[6]，如圖4。氧氣濃度大于15%的范圍為散熱帶，漏風較大，不利于熱量累積，氧氣濃度為5%～15%的范圍為氧化升溫帶，氧氣濃度小于5%為窒息帶。可以看出進風側的氧化升溫帶寬度明顯比回風側的寬。

3.2.2 瞬態結果

使用動網格方法進行瞬態模擬，模擬采空區從200m 的初始長度推進至380m 的工作面動態推進過程，得到采空區氧氣濃度場的演變情況如圖5。初期采空區內氧氣濃度場變化較明顯，但隨著開采進行，言其濃度場逐漸趨于穩定，不再隨時間發生明顯變化。在氧氣濃度場不再隨時間發生明顯變化后，按照氧氣體積濃度指標劃分成為的自燃“三帶”，如圖6 所示。

圖3 采空區氧氣濃度場

圖4 采空區自燃“三帶”

圖5 工作面動態回采過程中的氧氣濃度場演變

圖6 采空區自燃“三帶”分布圖

3.2.3 結果分析

通過對穩態和瞬態結果進行對比可知，兩種模擬方法所得到的結果有一定差距。相較于穩態計算結果，瞬態計算出的氧化升溫帶的范圍更寬，而使用穩態方法進行計算的結果比較接近于使用瞬態方法的計算初期的結果。

由于瞬態計算即是在每個時間步內進行近似穩態計算，計算完整個計算域后再將所得的結果作為下一個時間步的初始值進行近似穩態計算，直至時間步的結束。穩態計算只需要最終迭代達到收斂即可，而瞬態計算則要求每個時間步內均達到收斂。若只考慮系統穩定后的狀態，那么就選擇使用穩態計算；若要考慮流場的演化情況，則需要使用瞬態。在實際研究中，選擇穩態方法還是瞬態方法取決于研究目標的側重點使用穩態計算的采空區自燃“三帶”的分布結果，表達了采空區各場處于穩定狀態時的情況，而使用FLUENT 動網格的瞬態計算，則更側重于工作面回采過程中氧氣濃度場的演變情況。

4 結論

4.1 通過對穩態和瞬態結果進行對比可知，兩種模擬方法所得到的結果有一定差距。使用穩態方法進行計算的結果比較接近于使用瞬態方法的計算初期的結果。

4.2 與穩態計算結果相比，使用FLUENT 動網格的瞬態計算，則更側重于工作面回采過程中氧氣濃度場的演變情況。