吳奉亮，高亞超，常心坦

(1.西安科技大學(xué)安全科學(xué)與工程學(xué)院，陜西西安，710054；2.西安科技大學(xué)西部礦井開采與災(zāi)害防治教育部重點(diǎn)實驗室，陜西西安，710054)

通風(fēng)不僅為井下礦工提供足夠新鮮空氣，而且是防控井下各種污染、爆炸和窒息性氣體的基本措施。對于煤礦，供給采煤工作面的風(fēng)流會連續(xù)漏入采空區(qū)，這種難以避免的漏風(fēng)極易引發(fā)采空區(qū)遺煤自燃及瓦斯爆炸。我國采煤工作面40%～60%的瓦斯來自采空區(qū)漏風(fēng)，超過60%的煤礦火災(zāi)來自于采空區(qū)遺煤自燃[1-2]；美國煤礦廣泛采用大量風(fēng)流流經(jīng)采空區(qū)的Bleeder 形式通風(fēng)系統(tǒng)[3]?？梢?，掌握采空區(qū)中因漏風(fēng)形成的危險區(qū)域?qū)Ψ揽孛旱V爆炸與火災(zāi)事故具有重要的作用。由于礦井通風(fēng)系統(tǒng)的復(fù)雜性，難以直接觀測采空區(qū)流場，數(shù)值模擬成為國內(nèi)外學(xué)者研究這2個部分井下風(fēng)流的主要手段。礦井通風(fēng)系統(tǒng)(包括風(fēng)道、風(fēng)機(jī)和控風(fēng)設(shè)施等元素)常被描述成1 個一維的網(wǎng)絡(luò)模型，風(fēng)流在其中的流動遵循節(jié)點(diǎn)流量平衡和回路風(fēng)壓平衡定律，描述該定律的非線性方程組以井巷的風(fēng)量或節(jié)點(diǎn)的壓力作為未知數(shù)，未知數(shù)的數(shù)量常超過1 000個。長期以來，研究者圍繞此非線性方程組提出了牛頓法、斯拷特-恒斯雷法等多種不同解算方法[4-5]，許多由行業(yè)研究者自主開發(fā)的具有良好可視化[6]和并行計算[7]能力的解算軟件已被廣泛用于理解與控制這種復(fù)雜的通風(fēng)系統(tǒng)。與巷道中風(fēng)流的一維網(wǎng)絡(luò)模型不同，采空區(qū)中的風(fēng)流主要采用二維或三維場模型來描述，在計算流體力學(xué)CFD(computational fluid dynamics)技術(shù)支撐下，大量模擬研究很好地展現(xiàn)了風(fēng)流在采空區(qū)這個不可見區(qū)域中的流動過程及其可能引發(fā)問題的危險區(qū)域，如模擬采空區(qū)流場[8-9]、判定采空區(qū)中煤自燃危險區(qū)域[1-2,10]、模擬采空區(qū)煤自燃[11-13]和瓦斯運(yùn)移過程[1-2]。但現(xiàn)有模型多是將采空區(qū)與礦井通風(fēng)網(wǎng)絡(luò)分離，沒有考慮漏風(fēng)對采空區(qū)邊界條件的影響。采空區(qū)流場將隨著邊界條件的改變而改變，通風(fēng)系統(tǒng)的改變也將導(dǎo)致預(yù)先定義的邊界條件失效，因此，整體求解井下巷道與采空區(qū)中的風(fēng)流流動是一件十分有意義的工作。人們一直在嘗試建立礦井風(fēng)網(wǎng)與采空區(qū)流場的統(tǒng)一求解模型。文獻(xiàn)[14-15]將采空區(qū)劃分成許多縱橫交錯的分支，與礦井風(fēng)網(wǎng)形成一個整體。但與目前采空區(qū)模擬使用的主流技術(shù)CFD 相比，這種方法不嚴(yán)密。近年已有學(xué)者嘗試將網(wǎng)絡(luò)解算模型添加到商業(yè)CFD軟件中[16]，或?qū)FD 算法移置到網(wǎng)絡(luò)解算軟件中[17]，這些工作在一定程度上顯示了通風(fēng)網(wǎng)絡(luò)與采空區(qū)流場集成求解的優(yōu)越性，但還沒有很好地考慮兩者的耦合關(guān)系。本文利用礦井風(fēng)網(wǎng)與采空區(qū)流場的邊界流量來建立兩者的耦合模型，并研究模型的求解方法。

1 礦井風(fēng)網(wǎng)與采空區(qū)流場的邊界耦合模型

1.1 整體模型

采空區(qū)中的風(fēng)流在冒落煤巖體中的流動是各向同性非均勻介質(zhì)滲流，由于采空區(qū)寬度與深度遠(yuǎn)大于高度，三維采空區(qū)流場完全可簡化為如圖1所示的平面流場。礦井通風(fēng)網(wǎng)絡(luò)解算需要在井巷系統(tǒng)之上建立風(fēng)網(wǎng)結(jié)構(gòu)，采用有限元求解采空區(qū)流場則需要將求解區(qū)域劃分成單元。以圖1所示典型的U 形工作面通風(fēng)系統(tǒng)為例，高度簡化礦井風(fēng)網(wǎng)，可以建立圖2所示的幾何模型，該模型同時包含礦井風(fēng)網(wǎng)與采空區(qū)流場。

圖1 U形工作面通風(fēng)系統(tǒng)Fig.1 Layout of working face with U-ventilation system

圖2中節(jié)點(diǎn)③～⑥是風(fēng)網(wǎng)的邊界節(jié)點(diǎn)，節(jié)點(diǎn)⑦～⑩是采空區(qū)的邊界節(jié)點(diǎn)；邊界流量q′ij表示風(fēng)網(wǎng)邊界節(jié)點(diǎn)i與采空區(qū)邊界單元j之間的風(fēng)流交換量。對于采空區(qū)(老空區(qū))的漏風(fēng)密閉與風(fēng)網(wǎng)之間的邊界也可建立類似的虛擬分支，因此，當(dāng)引入足夠多的邊界節(jié)點(diǎn)與虛擬分支時，常規(guī)模型中分離的風(fēng)網(wǎng)與采空區(qū)流場在邊界上耦合成一個整體。采空區(qū)邊界節(jié)點(diǎn)壓力可通過風(fēng)網(wǎng)邊界節(jié)點(diǎn)計算，即圖2中的節(jié)點(diǎn)⑦可取節(jié)點(diǎn)③的壓力，節(jié)點(diǎn)⑧的壓力為節(jié)點(diǎn)③和④的壓力平均值。礦井風(fēng)網(wǎng)與采空區(qū)流場的邊界耦合模型將基于采空區(qū)流場及風(fēng)網(wǎng)的求解方法來研究如何求得q′ij，以及如何將q′ij反饋到風(fēng)網(wǎng)解算模型中。

1.2 邊界耦合模型及解耦流程

圖2 礦井風(fēng)網(wǎng)與采空區(qū)流場邊界耦合模型示意圖Fig.2 Schematic diagram of boundary coupling model between ventilation network and gob flow field

由于很難導(dǎo)出關(guān)于圖2中所有節(jié)點(diǎn)壓力或流量流速為未知數(shù)的控制方程組，采用分離計算的辦法求解整體模型。新風(fēng)網(wǎng)(圖2中風(fēng)網(wǎng)與邊界)的節(jié)點(diǎn)流量平衡方程將與常規(guī)通風(fēng)網(wǎng)絡(luò)不同，此處應(yīng)為

式中：N為風(fēng)網(wǎng)節(jié)點(diǎn)數(shù)；qij為連接結(jié)點(diǎn)i和j的風(fēng)網(wǎng)分支的流量，m3/s，由i和j節(jié)點(diǎn)的壓力來計算：

式中：Rij為連接結(jié)點(diǎn)i和j的分支的風(fēng)阻，kg/m7；pi和pj分別為節(jié)點(diǎn)i和j的壓力，Pa。式(1)左側(cè)是連接節(jié)點(diǎn)i的所有分支的流量之和，右側(cè)為來自采空區(qū)的流量。對于普通的通風(fēng)網(wǎng)絡(luò)，即圖2中的風(fēng)網(wǎng)部分，式(1)右側(cè)恒等于0。將式(2)代入式(1)得到的方程組以風(fēng)網(wǎng)節(jié)點(diǎn)壓力為未知數(shù)，因此，解此方程組來求解通風(fēng)網(wǎng)絡(luò)的方法稱為節(jié)點(diǎn)壓力法。

式(1)中加入邊界流量qik'，一方面，導(dǎo)致新風(fēng)網(wǎng)不再是封閉的，另一方面，在網(wǎng)絡(luò)解算與采空區(qū)流場模擬之間形成了相互依賴的耦合關(guān)系：采空區(qū)流場的模擬依賴于對新網(wǎng)絡(luò)的解算為其提供邊界上的節(jié)點(diǎn)壓力；同時新網(wǎng)絡(luò)解算，即式(1)中的邊界流量qik'又依賴于采空區(qū)流場模擬結(jié)果來確定。盡管2個部分看似僅在邊界上耦合，但求解過程不得不涉及完整的風(fēng)網(wǎng)及整個采空區(qū)。

圖3所示為風(fēng)網(wǎng)與采空區(qū)流場的迭代解耦流程。圖3中為第k次迭代中求得的所有邊界單元的邊界流量qik'構(gòu)成的向量；D為與之差的元素的絕對值的最大值，當(dāng)D小于給定精度ε時，求得風(fēng)網(wǎng)與采空區(qū)流場耦合模型的解。因此，從分離求解來看，圖2整體模型也是一個風(fēng)網(wǎng)與采空區(qū)流場之間的耦合模型。

需要指出的是：1)新網(wǎng)絡(luò)的不封閉性導(dǎo)致無法使用常規(guī)的回路風(fēng)量法進(jìn)行求解，對此本文采用風(fēng)網(wǎng)結(jié)點(diǎn)壓力法[18]來解決；2)當(dāng)多個采空區(qū)按圖2所示方式連接到風(fēng)網(wǎng)中時，圖3所示求解流程仍然有效，由于礦井風(fēng)網(wǎng)表示了所有井巷的風(fēng)路，以上用1個采空區(qū)為例的推導(dǎo)過程，實則建立了風(fēng)網(wǎng)與多采空區(qū)流場的邊界耦合模型及求解流程。自編程研究采空區(qū)漏風(fēng)流場，需要解算大型線型方程組，當(dāng)多個采空區(qū)連接到風(fēng)網(wǎng)時，解方程的計算量是普通采空區(qū)流場模擬計算量的若干倍。

圖3 風(fēng)網(wǎng)與采空區(qū)流場邊界耦合模型的迭代求解流程Fig.3 Iteration flow between ventilation network calculation and gob flow field simulation

2 有限元法求解采空區(qū)流場的控制方程組

采空區(qū)流場的控制方程組是指以采空區(qū)所有待求結(jié)點(diǎn)壓力為未知數(shù)的代數(shù)方程組。若將風(fēng)流在采空區(qū)中的流動先假設(shè)為線性滲流，則圖1采空區(qū)流場的定解數(shù)學(xué)模型為下列的微分方程：

式中：P為采空區(qū)中的壓力函數(shù)，Pa；K為滲流系數(shù)，m/s；G為整個采空區(qū)區(qū)域；p0為L1上的已知壓力；nx和ny分別為邊界外法線上單位向量在x和y方向的分量。將圖2中三角形單元的3 個節(jié)點(diǎn)按逆時針方向依次編號為i，j和k，相應(yīng)壓力分別記為pi，pj和pk。記L1邊界上的節(jié)點(diǎn)(壓力已知)數(shù)為w，不在L1邊界上的采空區(qū)節(jié)點(diǎn)(壓力待求)數(shù)為s；按待求在前、已知在后的順序?qū)⑺胁煽諈^(qū)節(jié)點(diǎn)壓力向量記為Pg=(Pu，Pb)，其中，待求壓力列向量Pu=(p1，p2，…，ps)T，已知壓力列向量Pb=(ps+1，…，ps+w)T?；谟邢拊碚撝械淖兎址ɑ蚣訖?quán)余量法均可導(dǎo)出關(guān)于Pg為未知數(shù)的控制方程組[19]：

式中：M稱為總體矩陣，可由每個單元的單元矩陣Ae合成。Ae是對稱矩陣

AΔ為單元面積，m2；KΔ為單元滲透系數(shù)；系數(shù)b和c由三角形單元的i，j和k這3 個節(jié)點(diǎn)的坐標(biāo)x和y決定，對應(yīng)關(guān)系分別是：bi=yj-yk，ci=xk-xj，bj=ykyi，cj=xi-xk，bk=yi-yj，ck=xj-xi。設(shè)mi，mj和mk分別表示單元的i，j和k節(jié)點(diǎn)的壓力在Pg向量中的序號，則Ae的1，2 和3 行或列分別與M的mi，mj和mk行或列對應(yīng)，合成是指單元矩陣元素累加到總體矩陣中。顯然，式(4)中Pg的元素只有前s個是未知數(shù)，因此，不用計算M的每一個元素。M具有圖4所示的分塊特點(diǎn)，即M的前s行、s列構(gòu)成的矩陣A是以采空區(qū)所有待求結(jié)點(diǎn)壓力為未知數(shù)的代數(shù)方程組(即控制方程組)的系數(shù)矩陣，因此，有

圖4 總體矩陣M的分塊表示Fig.4 The global matrix represented with block matrixes

式(6)右端列向量b可由分塊矩陣B與Pb求得，即b=BPb。因此，在合成M時只需計算A和B分塊。

事實上，采空區(qū)內(nèi)，特別是鄰近工作面的范圍，一般為非線性的滲流，可用式(7)所示的Bachmat方程描述：

式中：β為介質(zhì)顆粒形狀系數(shù)；d為平均調(diào)和粒徑，m；φ為孔隙度；υ為運(yùn)動黏性系數(shù)，m2/s。本文采用文獻(xiàn)[9]提出的對滲透系數(shù)進(jìn)行迭代修正的方法，基于式(3)線型模型求解式(7)的非線性模型。

通過式(6)求得節(jié)點(diǎn)壓力后，可通過達(dá)西定律公式V=-K?P求得邊界單元的流速，通過式(8)可求得圖2所示的邊界流量：

式中：lj為采空區(qū)邊界單元j在邊界上的邊長，m；h為采空區(qū)厚度，m；vjx和vjy分別為單元j在x和y方向上流速的分量，m/s；njx和njy分別為單元j所在邊界的外法線方向的單位矢量在x和y方向的分量。

綜上，求解耦合模型的計算量主要來自解(6)式大型稀疏線性方程組，它存在于非線性滲流迭代與圖3迭代的雙重循環(huán)之內(nèi)，再加上多采空區(qū)同時耦合于礦井風(fēng)網(wǎng)，其計算量是普通采空區(qū)流場模擬計算量的若干倍。

3 基于PARDISO的耦合模型求解

3.1 PARDISO特點(diǎn)

PARDISO 是針對大型稀疏線性方程組的求解器[20]。在算法上，采用LU 分解法；在性能上，采用共享內(nèi)存和分布式內(nèi)存的并行計算；在模塊化方面，PARDISO 以動態(tài)鏈接庫向外提供封裝好的函數(shù)，可以被VC++等多種語言訪問。此外，PARDISO 采用壓縮稀疏行格式(compressed sparse row，CSR)和壓縮稀疏列格式(compressed sparse column，CSC)存儲系數(shù)矩陣[21]，降低對矩陣元素的訪問時間。這2 種格式是利用PARDISO 進(jìn)行耦合模型求解的關(guān)鍵數(shù)據(jù)結(jié)構(gòu)。

3.2 基于PARDISO的耦合模型算例分析

按圖2所示原理，構(gòu)造圖5所示算例，風(fēng)網(wǎng)含有4條分支，分支(3)代表采煤工作面分支；采空區(qū)劃分為8 個單元，單元直角邊長均為150 m，分支(1)設(shè)為固定風(fēng)量分支，初始風(fēng)量取30 m3/s。為便于對算例結(jié)果重現(xiàn)，對相關(guān)物理參數(shù)進(jìn)行簡單設(shè)定，分支的風(fēng)阻均為0.05 kg/m7，單元的滲透系數(shù)KΔ均取0.1 m/s，采空區(qū)厚度h為6 m，平均粒徑d為0.04 m，介質(zhì)顆粒形狀系數(shù)β取0.07，空氣的運(yùn)動黏性系數(shù)υ為14.6×10-6m2/s。

圖5 按圖2所示原理構(gòu)造的簡單算例Fig.5 A simple example constructed according to the principle shown in Fig.2

取節(jié)點(diǎn)①為基點(diǎn)，壓力為0 Pa，按照耦合模型原理，第1次解算風(fēng)網(wǎng)得到節(jié)點(diǎn)②，③和④的壓力分別為135，90 和45 Pa，則采空區(qū)邊界節(jié)點(diǎn)⑤，⑧和?的壓力分別為90，67.5 和45 Pa。首次解算采空區(qū)時，第1次迭代中得到的式(6)系數(shù)矩陣A如圖6所示，圖6中空格表示絕對的0 值，即沒有任何單元矩陣的元素合成到這些位置；非空格中的0值是由計算產(chǎn)生的，在條件變化時，其值可能不等于0。

圖6 算例迭代中首次計算的系數(shù)矩陣Fig.6 Coefficient matrix for the first time in iterative calculation

采用PARDISO 解算式(6)，需將圖6所示矩陣的非0元素存入表1所示的3個一維數(shù)組a，Ai和Aj中。其中a存放元素值，Ai和Aj存放元素的位置信息，轉(zhuǎn)存規(guī)則是：在圖6中按從上到下逐行、從左到右逐列的順序依次掃描圖中的非0 元素，存入a中。Ai[k]表示圖6中第k行第1個非0元素在a中的序號，Aj[k]存放a[k]元素在圖中的列號。用表1所示的數(shù)據(jù)結(jié)構(gòu)來保存圖6中的非0元素，減少了大量的存儲空間。從圖6可知：式(6)實際是帶狀方程組，帶寬取決于各個單元結(jié)點(diǎn)序號差的最大值。為減少矩陣帶寬，在有限元網(wǎng)格劃分理論中常施加一定方法來降低帶寬，提高性能。采用PARDISO 引入的以上壓縮存儲格式后，不用再考慮帶寬對性能的影響，這在一定程度上降低了對采空區(qū)網(wǎng)格劃分的要求。

表1 圖6矩陣壓縮成的3個1維數(shù)組Table 1 Three one-dimensional arrays conforming to PARDISO rules

將式(6)的右端向量b和以上3個一維數(shù)組傳入PARDISO的接口，即可實現(xiàn)基于PARDISO的線性方程組快速求解。圖5算例在風(fēng)網(wǎng)解算與采空區(qū)流場模擬之間迭代8次后收斂，圖7所示為風(fēng)網(wǎng)分支流量、采空區(qū)單元流速及兩者之間的風(fēng)流交換量，節(jié)點(diǎn)(2)至(13)的壓力依次為137.86，92.86，47.78，92.86，75 .65，73.34，70.32，70.82，71.02，47.78，67.48和69.49 Pa。

4 算法實現(xiàn)與復(fù)雜算例分析

4.1 軟件開發(fā)和算例參數(shù)

本文基于VC++編制以上模型的軟件，簡稱i-MVS (integrated mine ventilation simulator)。i-MVS采用ObjectARX 技術(shù)[17]進(jìn)行風(fēng)網(wǎng)與采空區(qū)的可視化，通過動態(tài)鏈接庫使用PARDISO，按孔隙率呈“O”形圈分布設(shè)置采空區(qū)滲透系數(shù)[22-24]，可對任意風(fēng)網(wǎng)與平面采空區(qū)流場進(jìn)行統(tǒng)一計算。

圖7 圖5算例的模擬結(jié)果Fig.7 The simulation results of the example shown in Fig.5

圖8所示為由某礦2個采空區(qū)及高度簡化后的風(fēng)網(wǎng)形成的算例。工作面的微小分支按每100 m風(fēng)阻為0.020 8 kg/m7設(shè)定；1號和2號這2個采煤工作面設(shè)計風(fēng)量分別為30 m3/s和20 m3/s；1號工作面斜長300 m，已回采500 m；2 號工作面斜長200 m，已回采300 m；采空區(qū)厚度均取6 m。2個工作面共分成161 條分支，連接節(jié)點(diǎn)①和②的分支(風(fēng)機(jī)所在位置)設(shè)為初始風(fēng)量為100 m3/s的固定風(fēng)量分支，連接節(jié)點(diǎn)⑩與○13的分支代表風(fēng)網(wǎng)的其他部分，設(shè)計過風(fēng)量為45 m3/s。風(fēng)網(wǎng)共設(shè)有3 個調(diào)風(fēng)裝置，1 號和2 號位于每個工作面進(jìn)、回風(fēng)巷的聯(lián)絡(luò)巷，3 號用于控制2 號工作面的回風(fēng)量。模型中涉及的其他參數(shù)取值與3.2節(jié)的相同。

4.2 模擬結(jié)果

4.2.1 耦合模型計算結(jié)果與性能分析

圖9所示為風(fēng)網(wǎng)及采空區(qū)在邊界處耦合后的流量及流速計算結(jié)果。2個采空區(qū)分別用50條流線與等壓線形成流網(wǎng)。不漏風(fēng)的邊界為流函數(shù)值為0的流線，1 號和2 號采空區(qū)流線間的流量差距分別為0.046 m3/s 和0.025 m3/s，2 個采空區(qū)總的漏風(fēng)量分別為2.3 m3/s 和1.2 m3/s。1 號和2 號采空區(qū)等壓線差距分別為1 Pa和0.3 Pa，節(jié)點(diǎn)①和②之間的壓差表明系統(tǒng)總阻力838.6 Pa。由圖9可見：距離工作面越遠(yuǎn)，流線數(shù)量越少，采空區(qū)流量越??；大部分流線起于進(jìn)風(fēng)隅角，止于回風(fēng)隅角，表明采空區(qū)漏風(fēng)主要從工作面的2 個隅角流入與流出采空區(qū)。

圖8 模擬算例基本情況Fig.8 Simulated case condition

圖9 風(fēng)網(wǎng)與采空區(qū)流場在邊界處耦合后的計算結(jié)果Fig.9 Simulation results of the boundary coupling model between ventilation network and gob flow field

采空區(qū)中風(fēng)速為0.001 7～0.004 0 m/s 的區(qū)域代表煤自燃氧化升溫帶，從圖9可知，1 號采空區(qū)氧化升溫帶起始位置位于工作面后26.5 m，中心寬度為33.4 m，沿兩側(cè)邊界寬度為147.1 m；2號采空區(qū)氧化升溫帶位于工作面后11.5 m，中心寬度為28.5 m，沿兩側(cè)邊界的寬度為90.0 m。因此，當(dāng)確定用于控制采空區(qū)自然發(fā)火的最小推進(jìn)度時，應(yīng)采用氧化升溫帶在采空區(qū)邊界上的寬度。

PARDISO求解技術(shù)的性能分析。2個采空區(qū)共有24 660個單元，12 939個節(jié)點(diǎn)，其中，采空區(qū)節(jié)點(diǎn)有12 766個。風(fēng)網(wǎng)與采空區(qū)之間經(jīng)過6次迭代后收斂(ε=10-5)，共進(jìn)行了860次對式(6)的求解；在4核CPU(頻率2.2 GHz)的個人電腦上，完成耦合模型計算約2 min。

4.2.2 漏風(fēng)對采空區(qū)邊界條件的影響

作為對比，將圖3中不進(jìn)行迭代、只進(jìn)行1次網(wǎng)絡(luò)解算和采空區(qū)模擬結(jié)果稱為非耦合模型，進(jìn)行常規(guī)的風(fēng)網(wǎng)與采空區(qū)流場模擬。圖10顯示了2個工作面在2 種模型漏風(fēng)邊界沿線的壓力與漏風(fēng)量，耦合模型考慮采空區(qū)漏風(fēng)對邊界條件的影響，得到比非耦合模型更小的漏風(fēng)壓差與漏風(fēng)量。從圖10(a)可見：1號工作面在非耦合模型中兩端的壓力為56.2 Pa，耦合模型中為49.9 Pa，差值為6.3 Pa，是總壓力的12.6%。由圖10(b)可見：2 號工作面在耦合與非耦合模型下工作面壓差分別為14.5 Pa 與16.2 Pa，耦合模型比非耦合模型壓差小1.7 Pa，占工作面總壓差的12.0%。

對漏風(fēng)量的比較可以得到近似的結(jié)論：1號工作面在耦合與非耦合模型中的總漏風(fēng)量分別為138.0 m3/min 與152.8 m3/min，相 差14.8 m3/min，是耦合模型總漏風(fēng)量的10.7%；2 號工作面在耦合與非模型中的總漏風(fēng)量分別為72.0 m3/min 與80.7 m3/min，相差8.7 m3/min，是耦合模型總漏風(fēng)量的12.1%。因此，12%可以視為耦合模型對常規(guī)采空區(qū)模擬方法的改進(jìn)量。U形工作面通風(fēng)系統(tǒng)的采空區(qū)漏風(fēng)常小于其他形式通風(fēng)系統(tǒng)(如Y 形通風(fēng)系統(tǒng))。故考慮漏風(fēng)對采空區(qū)流場邊界條件的影響是必要的。

4.2.3 異常通風(fēng)模擬

圖10 工作面沿線壓力與采空區(qū)漏風(fēng)量分布Fig.10 Air pressures and air leakage along the working faces

圖11 1號調(diào)節(jié)裝置損壞后的風(fēng)網(wǎng)流量與采空區(qū)速度場Fig.11 Network air distributions and gob flow field colored after NO 1 regulator damaged

作為對一個異常情況的分析，圖11模擬了1號調(diào)風(fēng)裝置意外破壞后的風(fēng)網(wǎng)流量與采空區(qū)速度場，其標(biāo)注信息與圖9中的相同。圖8中位于節(jié)點(diǎn)?和?之間的分支的風(fēng)阻為15.323 0 kg/m7，其中包含了1 號調(diào)風(fēng)裝置風(fēng)阻(15.293 8 kg/m7)和巷道風(fēng)阻(0.029 2 kg/m7)。將此分支風(fēng)阻從15.323 0 kg/m7降到0.029 2 kg/m7得到此異常情況下的風(fēng)網(wǎng)與采空區(qū)流場，即圖11所示結(jié)果。從圖11可知：1 號調(diào)風(fēng)裝置的損毀對2號工作面通風(fēng)系統(tǒng)的影響相對較小，但1號工作面通風(fēng)系統(tǒng)幾乎被短路，47.00 m3/s的風(fēng)流從節(jié)點(diǎn)?直接流向節(jié)點(diǎn)?，工作面風(fēng)量從30.20 m3/s 減到12.40 m3/s，同時，1 號采空區(qū)漏風(fēng)量從2.30 m3/s降到0.79 m3/s。工作面風(fēng)量及采空區(qū)漏風(fēng)量的改變將影響這些區(qū)域瓦斯?jié)舛?，假設(shè)采空區(qū)瓦斯涌出量不變，采空區(qū)漏風(fēng)量的減小將導(dǎo)致從1號采空區(qū)涌出氣體的平均瓦斯?jié)舛仍黾訛樵瓉淼?.9 倍，加上工作面風(fēng)量減小，1 號工作面的危險程度將明顯增加。

5 結(jié)論

1)采煤工作面兩端的壓差是采空區(qū)漏風(fēng)的動力，沿線的壓力分布是采空區(qū)流場的定解條件、需要依靠對風(fēng)網(wǎng)的求解來確定；而采空區(qū)漏風(fēng)又會減小工作面的壓差、改變沿線的壓力分布。耦合模型將工作面細(xì)分成與采空區(qū)單元匹配的分支，通過將采空區(qū)邊界沿線的漏風(fēng)分段加入風(fēng)網(wǎng)邊界節(jié)點(diǎn)的流量平衡方程，導(dǎo)致風(fēng)網(wǎng)解算與采空區(qū)流場模擬形成耦合關(guān)系，通過對兩者的迭代計算求得更加準(zhǔn)確的解。

2)耦合模型實現(xiàn)了風(fēng)網(wǎng)與多采空區(qū)的統(tǒng)一求解。礦井通風(fēng)網(wǎng)絡(luò)包括全礦井的井巷系統(tǒng)，耦合模型實現(xiàn)了采空區(qū)流場與風(fēng)網(wǎng)在邊界上的無縫連接，使過去分離解算的采空區(qū)流場可以與風(fēng)網(wǎng)構(gòu)成一個整體，不用單獨(dú)為某個采空區(qū)流場設(shè)置邊界條件。

3)使用PARDISO 實現(xiàn)了對耦合模型的高效求解。采用自編程方法實現(xiàn)模型，將采空區(qū)流場模擬與礦井通風(fēng)網(wǎng)絡(luò)解算功能進(jìn)行了集成。采用PARDISO 求解耦合模型中的大型稀疏線性方程組保證了軟件的性能；使用PARDISO 引入的壓縮行、列存儲格式使得在采空區(qū)的網(wǎng)格劃分中不用考慮單元的結(jié)點(diǎn)編號對求解性能的影響。