礦井風(fēng)網(wǎng)與采空區(qū)流場(chǎng)邊界耦合模型與求解技術(shù)

吳奉亮，高亞超，常心坦

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

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

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

1.1 整體模型

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

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

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

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

由于很難導(dǎo)出關(guān)于圖2中所有節(jié)點(diǎn)壓力或流量流速為未知數(shù)的控制方程組，采用分離計(jì)算的辦法求解整體模型。新風(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)的壓力來(lái)計(jì)算：

式中：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è)為來(lái)自采空區(qū)的流量。對(duì)于普通的通風(fēng)網(wǎng)絡(luò)，即圖2中的風(fēng)網(wǎng)部分，式(1)右側(cè)恒等于0。將式(2)代入式(1)得到的方程組以風(fēng)網(wǎng)節(jié)點(diǎn)壓力為未知數(shù)，因此，解此方程組來(lái)求解通風(fēng)網(wǎng)絡(luò)的方法稱(chēng)為節(jié)點(diǎn)壓力法。

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

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

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

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

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

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

式中：P為采空區(qū)中的壓力函數(shù)，Pa；K為滲流系數(shù)，m/s；G為整個(gè)采空區(qū)區(qū)域；p0為L(zhǎng)1上的已知壓力；nx和ny分別為邊界外法線上單位向量在x和y方向的分量。將圖2中三角形單元的3 個(gè)節(jié)點(diǎn)按逆時(shí)針?lè)较蛞来尉幪?hào)為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稱(chēng)為總體矩陣，可由每個(gè)單元的單元矩陣Ae合成。Ae是對(duì)稱(chēng)矩陣

AΔ為單元面積，m2；KΔ為單元滲透系數(shù)；系數(shù)b和c由三角形單元的i，j和k這3 個(gè)節(jié)點(diǎn)的坐標(biāo)x和y決定，對(duì)應(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向量中的序號(hào)，則Ae的1，2 和3 行或列分別與M的mi，mj和mk行或列對(duì)應(yīng)，合成是指單元矩陣元素累加到總體矩陣中。顯然，式(4)中Pg的元素只有前s個(gè)是未知數(shù)，因此，不用計(jì)算M的每一個(gè)元素。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時(shí)只需計(jì)算A和B分塊。

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

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

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

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

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

3 基于PARDISO的耦合模型求解

3.1 PARDISO特點(diǎn)

PARDISO 是針對(duì)大型稀疏線性方程組的求解器[20]。在算法上，采用LU 分解法；在性能上，采用共享內(nèi)存和分布式內(nèi)存的并行計(jì)算；在模塊化方面，PARDISO 以動(dòng)態(tài)鏈接庫(kù)向外提供封裝好的函數(shù)，可以被VC++等多種語(yǔ)言訪問(wèn)。此外，PARDISO 采用壓縮稀疏行格式(compressed sparse row，CSR)和壓縮稀疏列格式(compressed sparse column，CSC)存儲(chǔ)系數(shù)矩陣[21]，降低對(duì)矩陣元素的訪問(wèn)時(shí)間。這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 個(gè)單元，單元直角邊長(zhǎng)均為150 m，分支(1)設(shè)為固定風(fēng)量分支，初始風(fēng)量取30 m3/s。為便于對(duì)算例結(jié)果重現(xiàn)，對(duì)相關(guān)物理參數(shù)進(jìn)行簡(jiǎ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)動(dòng)黏性系數(shù)υ為14.6×10-6m2/s。

圖5 按圖2所示原理構(gòu)造的簡(jiǎn)單算例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ū)時(shí)，第1次迭代中得到的式(6)系數(shù)矩陣A如圖6所示，圖6中空格表示絕對(duì)的0 值，即沒(méi)有任何單元矩陣的元素合成到這些位置；非空格中的0值是由計(jì)算產(chǎn)生的，在條件變化時(shí)，其值可能不等于0。

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

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

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

將式(6)的右端向量b和以上3個(gè)一維數(shù)組傳入PARDISO的接口，即可實(shí)現(xiàn)基于PARDISO的線性方程組快速求解。圖5算例在風(fēng)網(wǎng)解算與采空區(qū)流場(chǎng)模擬之間迭代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 算法實(shí)現(xiàn)與復(fù)雜算例分析

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

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

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

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

4.2 模擬結(jié)果

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

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

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

圖9 風(fēng)網(wǎng)與采空區(qū)流場(chǎng)在邊界處耦合后的計(jì)算結(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 號(hào)采空區(qū)氧化升溫帶起始位置位于工作面后26.5 m，中心寬度為33.4 m，沿兩側(cè)邊界寬度為147.1 m；2號(hào)采空區(qū)氧化升溫帶位于工作面后11.5 m，中心寬度為28.5 m，沿兩側(cè)邊界的寬度為90.0 m。因此，當(dāng)確定用于控制采空區(qū)自然發(fā)火的最小推進(jìn)度時(shí)，應(yīng)采用氧化升溫帶在采空區(qū)邊界上的寬度。

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

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

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

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

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

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

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

作為對(duì)一個(gè)異常情況的分析，圖11模擬了1號(hào)調(diào)風(fēng)裝置意外破壞后的風(fēng)網(wǎng)流量與采空區(qū)速度場(chǎng)，其標(biāo)注信息與圖9中的相同。圖8中位于節(jié)點(diǎn)?和?之間的分支的風(fēng)阻為15.323 0 kg/m7，其中包含了1 號(hào)調(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ū)流場(chǎng)，即圖11所示結(jié)果。從圖11可知：1 號(hào)調(diào)風(fēng)裝置的損毀對(duì)2號(hào)工作面通風(fēng)系統(tǒng)的影響相對(duì)較小，但1號(hào)工作面通風(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，同時(shí)，1 號(hào)采空區(qū)漏風(fēng)量從2.30 m3/s降到0.79 m3/s。工作面風(fēng)量及采空區(qū)漏風(fēng)量的改變將影響這些區(qū)域瓦斯?jié)舛龋僭O(shè)采空區(qū)瓦斯涌出量不變，采空區(qū)漏風(fēng)量的減小將導(dǎo)致從1號(hào)采空區(qū)涌出氣體的平均瓦斯?jié)舛仍黾訛樵瓉?lái)的2.9 倍，加上工作面風(fēng)量減小，1 號(hào)工作面的危險(xiǎn)程度將明顯增加。

5 結(jié)論

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

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

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