基于FLUENT 動網(wǎng)格的煤礦采空區(qū)三帶動態(tài)特性模擬技術(shù)

王芳

（西安科技大學(xué)，陜西 西安710000）

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

1 采空區(qū)自燃三帶數(shù)值模擬

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

2 FLUENT 動網(wǎng)格技術(shù)

采空區(qū)范圍隨工作面推進(jìn)而不斷擴(kuò)大的過程是計(jì)算區(qū)域不斷擴(kuò)大的過程，在Fluent 中解決該問題的方法通常為采用動網(wǎng)格模型。Fluent 中的動網(wǎng)格模型可用來模擬由于流體域邊界剛性運(yùn)動或者邊界變形引起的流體域形狀隨時(shí)間變化的問題，每一個(gè)時(shí)間步上的體網(wǎng)格的更新是由解算器根據(jù)邊界的新的位置確定的，根據(jù)邊界的運(yùn)動自動調(diào)節(jié)內(nèi)部體網(wǎng)格的分布，通過指定運(yùn)動區(qū)域的運(yùn)動方式以及網(wǎng)格再生的方法可以使網(wǎng)格動態(tài)變化。

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

式中，V 為空間中大小和形狀都隨時(shí)間變化的控制體積，?V為控制體積的運(yùn)動邊界，ug為運(yùn)動網(wǎng)格的運(yùn)動速度；ρ為流體密度；u 為流體速度；г 為耗散系數(shù)；S φ 是標(biāo)量 φ的源項(xiàng)。

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

3 采空區(qū)自燃三帶數(shù)值模擬

3.1 模擬過程

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

圖1 采空區(qū)模型

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

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

圖2 采空區(qū)模型網(wǎng)格劃分

工作面推進(jìn)的過程即可視為簡單的速度- 時(shí)間運(yùn)動，工作面的移架、放頂雖然不是連續(xù)進(jìn)行的，但是從一個(gè)很長的時(shí)間周期角度來說，開采過程中每日進(jìn)刀數(shù)是較為固定的，因此可以將工作面相對于采空區(qū)的運(yùn)動看作為勻速運(yùn)動，可以通過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))

網(wǎng)格再生方法主要用來計(jì)算內(nèi)部網(wǎng)格節(jié)點(diǎn)的調(diào)節(jié)，采用鋪層（Layering）與局部重構(gòu)（Remeshing）相結(jié)合的方法進(jìn)行計(jì)算。動態(tài)模型的模擬采用瞬態(tài)模式，根據(jù)所需要模擬的時(shí)長來設(shè)置時(shí)間步長以及時(shí)間步數(shù)。

時(shí)間步長大致由以下公式估計(jì)：

3.2 數(shù)值模擬結(jié)果

3.2.1 穩(wěn)態(tài)結(jié)果

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

3.2.2 瞬態(tài)結(jié)果

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

圖3 采空區(qū)氧氣濃度場

圖4 采空區(qū)自燃“三帶”

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

圖6 采空區(qū)自燃“三帶”分布圖

3.2.3 結(jié)果分析

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

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

4 結(jié)論

4.1 通過對穩(wěn)態(tài)和瞬態(tài)結(jié)果進(jìn)行對比可知，兩種模擬方法所得到的結(jié)果有一定差距。使用穩(wěn)態(tài)方法進(jìn)行計(jì)算的結(jié)果比較接近于使用瞬態(tài)方法的計(jì)算初期的結(jié)果。

4.2 與穩(wěn)態(tài)計(jì)算結(jié)果相比，使用FLUENT 動網(wǎng)格的瞬態(tài)計(jì)算，則更側(cè)重于工作面回采過程中氧氣濃度場的演變情況。