孫曉元 王 川 孫英峰 邢云峰

（1.太原科技大學(xué)環(huán)境與安全學(xué)院，山西省太原市，030024；2.中國(guó)礦業(yè)大學(xué) （北京）資源與安全工程學(xué)院，北京市海淀區(qū)，100083；3.北京起重運(yùn)輸機(jī)械設(shè)計(jì)研究院，北京市東城區(qū)，100007）

受地質(zhì)構(gòu)造、煤層賦存及采掘方式多方面因素的影響，回采工作面及采空區(qū)的瓦斯流場(chǎng)表現(xiàn)為明顯的不均衡性和多變性。這種特點(diǎn)在以較快開采速度和較高開采強(qiáng)度為特征的現(xiàn)代綜合機(jī)械化采煤工藝中體現(xiàn)的尤為突出。作為典型的高瓦斯且具有煤與瓦斯突出危險(xiǎn)的礦井，華晉焦煤集團(tuán)沙曲礦開采的4＃煤層瓦斯壓力為1.52～1.57 MPa，瓦斯含量為7.3～17.82m3／t。根據(jù)現(xiàn)場(chǎng)測(cè)定，煤層透氣性系數(shù)為1.577～3.999 m3／MPa2·d，衰減系數(shù)0.01～0.016 （100d-1），抽放難度較大。而傳統(tǒng)的埋管抽放方式受空間非連續(xù)性的限制，尚未達(dá)到理想的抽放效果。鑒于此，亟需采用拖管抽采技術(shù)來(lái)改進(jìn)采空區(qū)抽放方法，并在對(duì)工作面和采空區(qū)瓦斯流分布特征進(jìn)行準(zhǔn)確的理論分析和數(shù)值模擬基礎(chǔ)上，確定最佳抽放點(diǎn)位置和合理抽放參數(shù)，為實(shí)際應(yīng)用提供技術(shù)支持。

1 風(fēng)流在采場(chǎng)中的運(yùn)動(dòng)模型

將風(fēng)流從進(jìn)入工作面進(jìn)風(fēng)巷，流經(jīng)工作面 （部分漏風(fēng)流經(jīng)采空區(qū)），然后從回風(fēng)巷排出的整個(gè)過(guò)程視為研究對(duì)象，該過(guò)程包括風(fēng)流在自由空間中的湍流過(guò)程和多孔介質(zhì) （采空區(qū)）的快速流動(dòng)過(guò)程，兩種過(guò)程所表現(xiàn)的流動(dòng)方式不同，遵循的控制方程也有所差異。

1.1 風(fēng)流在自由空間中的流動(dòng)方程

由于巷道和工作面正常通風(fēng)時(shí)其雷諾準(zhǔn)數(shù)遠(yuǎn)大于臨界雷諾數(shù)，因此國(guó)內(nèi)外學(xué)者往往從Navierstokes方程出發(fā)來(lái)描述其物理性質(zhì)。當(dāng)風(fēng)流為定常流動(dòng) （定常粘度、不可壓縮，即η，ρ=const）時(shí)，其平面二維動(dòng)量表達(dá)式 （N-S方程）為：

式中：u——風(fēng)流在自由空間中的速度矢量，m／s；

ρ——風(fēng)流密度，kg／m3；

p——風(fēng)流壓力，Pa；

ν——風(fēng)流的運(yùn)動(dòng)粘度，m2／s；

f——單位質(zhì)量風(fēng)流的體積力，N／kg。

式中：η——風(fēng)流的動(dòng)力粘度，Pa·s。

由式 （2）和 （3）可以看出，定常流動(dòng)N-S方程的因變量包括速度矢量u和壓力標(biāo)量p，為討論方便，將其用下標(biāo) “ns”來(lái)表示。

1.2 風(fēng)流在多孔介質(zhì)中的流動(dòng)方程

根據(jù)O 形圈理論，在距離工作面30m 范圍內(nèi)的采空區(qū)裂隙較為發(fā)育，且此范圍內(nèi)有較高速度的漏風(fēng)流經(jīng)過(guò)。實(shí)測(cè)表明，采空區(qū)內(nèi)的風(fēng)流基本上屬于介于N-S自由流動(dòng)和Darcy滲流之間的快速過(guò)渡流動(dòng)，因此，可以選用Brinkman方程對(duì)其進(jìn)行描述：

式中：ε——多孔介質(zhì)的孔隙率；

I——單位矢量；

F——流體阻力，N。

同樣，根據(jù)質(zhì)量守恒方程可以得出：

式（5）和 （6）組成了風(fēng)流在多孔介質(zhì)中運(yùn)移的控制方程，可以看出，Brinkman方程的因變量包括速度矢量u和壓力標(biāo)量p，為討論方便，將其用下標(biāo) “br”來(lái)表示。

1.3 風(fēng)流運(yùn)動(dòng)模型

通過(guò)上述分析可知，利用N-S方程和Brinkman方程可以較好描述風(fēng)流在自由空間和多孔介質(zhì)中的流動(dòng)過(guò)程，因此可以將兩者聯(lián)合起來(lái)構(gòu)建風(fēng)流從進(jìn)入進(jìn)風(fēng)巷道并流經(jīng)工作面 （部分經(jīng)過(guò)采空區(qū)）直至流出回風(fēng)巷道的運(yùn)動(dòng)模型，如圖1所示。實(shí)際上，將控制方程 （2）、（3）、（5）和 （6）聯(lián)立起來(lái)即可對(duì)工作面及采空區(qū)的流動(dòng)速度和壓力分布進(jìn)行求解，但要注意邊界的連續(xù)性，即在漏風(fēng)界面上存在uns=ubr和pns=pbr。

在進(jìn)行理論分析和數(shù)值模擬時(shí)，需考慮采空區(qū)滲透率和孔隙率的分布問(wèn)題，即受O 形圈影響不同的區(qū)域會(huì)有所差異。一般情況下，這些參數(shù)和碎脹性系數(shù)Kp（x，y）存在函數(shù)關(guān)系。而碎脹性系數(shù)的表達(dá)式為：

圖1 風(fēng)流運(yùn)動(dòng)模型

式中：Kp，0——初次冒落的碎脹系數(shù)；

Kp，1——壓實(shí)后的碎脹系數(shù)；

a0、a1、ξ1 ——調(diào)整系數(shù)；

L——工作面長(zhǎng)度，m。

2 瓦斯在采場(chǎng)中的運(yùn)移模型

上述討論研究了風(fēng)流在采場(chǎng)中的運(yùn)動(dòng)規(guī)律，然而往往更要關(guān)注的是瓦斯在采場(chǎng)中的運(yùn)移特征。為此，將瓦斯在采空區(qū)的運(yùn)移分為濃度擴(kuò)散和對(duì)流運(yùn)移兩個(gè)過(guò)程來(lái)進(jìn)行討論，這是因?yàn)槿魡渭兛紤]Fick擴(kuò)散第二定律來(lái)描述平面二維尺度空間內(nèi)的瓦斯非穩(wěn)定流動(dòng)場(chǎng)時(shí)，其方程表達(dá)式為：

式中：c——瓦斯?jié)舛龋琺ol／m3；

D——擴(kuò)散系數(shù)。

顯然，式（8）是基于純濃度梯度 （即濃度為唯一驅(qū)動(dòng)力）和質(zhì)量守恒定理得到的。但對(duì)于所要重點(diǎn)研究的采空區(qū)上隅角及其附近裂隙較為發(fā)育的區(qū)域則有所不妥，這是因?yàn)樵搮^(qū)域有較高速度的流體經(jīng)過(guò)，瓦斯運(yùn)移的驅(qū)動(dòng)力還應(yīng)包括通風(fēng)對(duì)流；該區(qū)域的瓦斯是存在有源補(bǔ)充的，匯源項(xiàng)包括遺煤及深部采空區(qū)不斷向外涌出的瓦斯。因此，需要將式（8）加以完善，改寫為：

式中：θs——瓦斯瞬時(shí)比例系數(shù)；Ri——匯源項(xiàng)；

u1——瓦斯平均流速，m／s。

關(guān)于平均流速的計(jì)算模型，文獻(xiàn)采用Darcy方程來(lái)進(jìn)行刻畫，鑒于本文研究區(qū)域的特殊性，故采用Brinkman方程來(lái)描述速度場(chǎng)分布，即聯(lián)立公式（5）、（7）和 （9）來(lái)探討采空區(qū)裂隙發(fā)育帶及上隅角附近的瓦斯運(yùn)移規(guī)律。

3 瓦斯在抽放滲流場(chǎng)中的運(yùn)移模型

在對(duì)瓦斯進(jìn)行抽放時(shí)，其有效抽放影響范圍內(nèi)瓦斯的流動(dòng)滿足質(zhì)量守恒方程：

式中：Q——單位體積煤巖體所含有的瓦斯量，

kg／t；

v1——鉆孔周邊的瓦斯?jié)B流速度，m／s；ρ1——瓦斯密度，kg／m3。

顯然，單位體積煤巖體中的瓦斯含量應(yīng)包括吸附瓦斯含量和游離瓦斯含量?jī)刹糠?。其中吸附瓦斯含量由Langmiur等溫吸附方程求得：

式中：Qa——吸附瓦斯含量，kg／t；

p1——瓦斯壓力，MPa；

a——吸附常數(shù)，單位體積煤體的極限瓦斯吸附量，m3／t；

b——吸附常數(shù)，MPa-1；

ρ0——標(biāo)準(zhǔn)大氣壓下的瓦斯密度，kg／m3。

與之相對(duì)應(yīng)的，若將瓦斯?jié)B流視為理想氣體在不考慮吸附解吸的等溫過(guò)程時(shí)，游離瓦斯含量為：

式中：Qf——游離瓦斯含量，kg／t；

ε——多孔介質(zhì)的孔隙率；

β——瓦斯氣體的壓縮因子，kg／m3·Pa。

顯然，存在公式：

假設(shè)瓦斯在煤層中的滲流模式為達(dá)西流動(dòng)，并考慮Klinkenberg效應(yīng)，則滲流的速度為：

式中：k——煤巖的滲透率，m2；

v1——滲流速度，m／s；

μ——瓦斯的動(dòng)力粘度系數(shù)，Pa；

?p——壓力梯度；

m——Klinkenberg等效因子。若考慮等溫過(guò)程多孔介質(zhì)孔隙率的變化：

式中：ζ——含瓦斯煤巖的體積應(yīng)變。

將式 （11）～ （15）聯(lián)立可得瓦斯在抽放滲流場(chǎng)中運(yùn)移的控制方程為：

4 數(shù)值模型及結(jié)果分析

通過(guò)上述分析，求得回采工作面風(fēng)流與瓦斯的流場(chǎng)分布及抽放滲流場(chǎng)控制方程的解析解是十分困難的，這是因?yàn)檫@些方程絕大部分都含有一階或二階偏微分項(xiàng)。因此，可以專門采用基于偏微分方程的有限元求解平臺(tái)COMSOL Multiphysics多物理場(chǎng)耦合分析工具。

4.1 幾何模型

算例選用沙曲礦14205綜采工作面，工作面長(zhǎng)度為158m，高度2.5m，進(jìn)風(fēng)巷和回風(fēng)巷凈斷面尺寸為4m×2.5m，采空區(qū)冒落帶高度為8m，為模擬O形圈內(nèi)的瓦斯運(yùn)移規(guī)律，走向長(zhǎng)度選為50m。

4.2 初始值和邊界條件

在14205綜采工作面正?；夭善陂g的生產(chǎn)班和非生產(chǎn)班多次測(cè)定其進(jìn)回風(fēng)巷和工作面的風(fēng)速、風(fēng)量和瓦斯?jié)舛鹊葏?shù)。通過(guò)統(tǒng)計(jì)分析，相關(guān)初始條件與其他邊界條件的取值：進(jìn)風(fēng)巷風(fēng)速172m／min，進(jìn)風(fēng)巷風(fēng)量1800m3／min，空氣密度1.29kg／m3，空氣動(dòng)力粘度5×10-5Pa·s，瓦斯動(dòng)力粘度1.34×10-5Pa·s，煤壁瓦斯涌出量34.21m3／min，落煤及遺煤瓦斯涌出量14.5m3／min，內(nèi)部采空區(qū)瓦斯涌出量49.85m3／min，工作面、采空區(qū)初始?xì)怏w壓力0.1MPa，采空區(qū)初始瓦斯?jié)舛?mol／m3，瓦斯分子擴(kuò)散系數(shù)0.0756 m2／h，采空區(qū)初始碎脹系數(shù)1.5，抽放管路負(fù)壓2kPa。

圖2 風(fēng)流在采空區(qū)和工作面運(yùn)動(dòng)模型的濃度分布云圖

4.3 模擬結(jié)果分析

單獨(dú)模擬風(fēng)流在采空區(qū)和工作面運(yùn)移時(shí)速度和濃度場(chǎng)分布，如圖2所示。其中在工作面運(yùn)動(dòng)方程采用N-S 方程，在采空區(qū) （多孔介質(zhì)）中采用Brinkman方程來(lái)刻畫速度場(chǎng)，將上述兩方程求得的速度場(chǎng)帶入空氣溶質(zhì)的對(duì)流傳遞中，以實(shí)現(xiàn)方程的耦合。圖2中橫坐標(biāo)代表模擬的不同時(shí)刻，從圖中可以看出，風(fēng)流大部分經(jīng)工作面由回風(fēng)巷排出，少部分風(fēng)流進(jìn)入采空區(qū)。由于采空區(qū)各處碎脹性系數(shù)設(shè)置的不同，其滲透率也有所差異。風(fēng)流在進(jìn)入采空區(qū)后，受風(fēng)流方向及滲透率的影響，采空區(qū)靠近進(jìn)風(fēng)巷及工作面的位置空氣濃度首先上升，原因是其區(qū)域位于O 形圈范圍內(nèi)。而中部區(qū)域由于滲透率較低，其濃度比O 形圈內(nèi)要低，見圖2 （c），但隨著時(shí)間步的推移，中部區(qū)域的濃度也逐漸上升，說(shuō)明瓦斯也逐步流入此區(qū)域，見圖2 （d），并最終達(dá)到平衡。

圖3模擬的是未采取抽放措施時(shí)采空區(qū)和工作面的瓦斯體積分?jǐn)?shù)分布云圖。其中瓦斯在采空區(qū)（多孔介質(zhì)）的平移擴(kuò)散時(shí)的速度場(chǎng)為Brinkman方程的模擬結(jié)果。由圖3可以看到，雖然在初始值設(shè)置的時(shí)候?qū)⒉煽諈^(qū)的瓦斯?jié)舛阮A(yù)定為3mol／m3，并且有內(nèi)部采空區(qū)和落煤 （遺煤）瓦斯的涌入，隨著時(shí)間的推移，采空區(qū)中高瓦斯?jié)舛人加械膮^(qū)域逐漸減小，并呈現(xiàn)出一定的規(guī)律性，即傾向方向上越接近進(jìn)風(fēng)巷瓦斯體積濃度越低，下降的速度越快，而越接近回風(fēng)巷道瓦斯體積濃度越高，且濃度下降的速率越慢；走向方向上瓦斯體積濃度隨深度的增大而逐漸增加；垂向方向上隨高度的增大瓦斯體積濃度也在增大；采空區(qū)中部滲透率較低區(qū)域的瓦斯?jié)舛认陆邓俾市∮贠 形圈范圍，說(shuō)明在裂隙發(fā)育帶的氣體流動(dòng)速率較快。

圖3 未采取抽放措施時(shí)采空區(qū)及工作面的瓦斯體積濃度分布云圖

圖4顯示的是抽放前后采空區(qū)及工作面瓦斯體積分?jǐn)?shù)與濃度等值面對(duì)比圖，這里采用了以縫代孔的方式來(lái)模擬抽放。通過(guò)分析可以發(fā)現(xiàn)，雖然單純利用通風(fēng)可以降低采空區(qū)及工作面的瓦斯，但一方面降低的速率較慢，另一方面在上隅角附近尤其是“三角區(qū)”范圍內(nèi)的瓦斯?jié)舛纫廊惠^大，如圖4（b）所示，仍存在發(fā)生瓦斯事故的危險(xiǎn)。因此，有必要采取抽放措施。圖4 （c）～4 （f）顯示了拖管抽放后的采空區(qū)及工作面瓦斯體積濃度等值面圖和分布云圖，對(duì)比圖3 （a）和圖3 （b），拖管抽放后上隅角瓦斯?jié)舛妊杆傧陆担榉?00000s后上隅角瓦斯?jié)舛冉禐?0.9%，其效果與未抽放1000000s時(shí)單純靠通風(fēng)降低瓦斯的效果 （9.27%）基本相似，當(dāng)拖管抽放1000000s時(shí)，上隅角附近的瓦斯?jié)舛纫呀档搅?.97%。經(jīng)過(guò)模擬發(fā)現(xiàn)，隨著時(shí)間的推移，當(dāng)拖管抽放時(shí)間達(dá)到2000000s時(shí)，上隅角的濃度降為0.9%，若能對(duì)上隅角進(jìn)行適當(dāng)?shù)姆舛?，減小該區(qū)域的滲透率，將能進(jìn)一步降低上隅角的瓦斯?jié)舛炔⑻岣叱榉判省?/p>

圖4 拖管抽放前后等值面與瓦斯體積分?jǐn)?shù)對(duì)比圖

圖5 拖管抽放深度及與回風(fēng)巷外幫距離參數(shù)的確定

為了確定14205綜采工作面拖管的合理位置，需要對(duì)不同高度、深度及距回風(fēng)巷外幫位置處的瓦斯抽放效果進(jìn)行模擬。顯然，抽放管內(nèi)的瓦斯?jié)舛入S高度的增加而增大，這是由瓦斯的升浮彌散特性所決定的。由此，只需要根據(jù)現(xiàn)場(chǎng)實(shí)際情況將瓦斯抽放拖管的位置抬到工程允許的適合高度即可，這里將其確定為2.3m。圖5顯示的是不同抽采深度和距離外幫不同位置時(shí)，上隅角及抽放管內(nèi)瓦斯?jié)舛龋ǔ榉艜r(shí)間均為2000000s），上隅角的取樣位置選擇在采空區(qū)縱深1m、距離底板2.3m、距回風(fēng)巷外幫0.3 m 處，抽放管的取樣位置選擇在抽放管內(nèi)距離瓦斯流出點(diǎn)1 m 處。由圖5 （a）可以看出，隨著抽放管進(jìn)入采空區(qū)深度的增加，拖管內(nèi)的瓦斯?jié)舛纫苍谥饾u上升，但當(dāng)拖管抽放口進(jìn)入采空區(qū)深度過(guò)大時(shí)，其距離工作面較遠(yuǎn)，上隅角位置的瓦斯?jié)舛瓤刂菩Ч^差，從安全考慮，選擇進(jìn)入采空區(qū)深度為15m 為宜；由圖5 （b）可以看出，由于采空區(qū)瓦斯流場(chǎng)分布的復(fù)雜性，抽放管并非越接近回風(fēng)巷外幫越好，確定準(zhǔn)則應(yīng)綜合上隅角瓦斯?jié)舛燃巴瞎軆?nèi)瓦斯?jié)舛染C合考量，經(jīng)過(guò)分析，選擇在距離外幫1.5m 處為宜。

綜上所述，沙曲礦14205綜采工作面拖管抽采的合理抽放參數(shù)為抽放高度2.3m，拖管進(jìn)入采空區(qū)深度15m，抽放位置距離回風(fēng)巷外幫1.5m。

5 結(jié)論

（1）采用COMSOL Multiphysics多物理場(chǎng)耦合分析工具可以合理準(zhǔn)確的研究回采工作面與采空區(qū)的速度場(chǎng)和濃度場(chǎng)特征，同時(shí)也可利用該工具來(lái)對(duì)拖管抽采措施進(jìn)行效果考察，從而為工程實(shí)踐提供了一定的參考依據(jù)。

（2）風(fēng)流在自由空間內(nèi)的運(yùn)動(dòng)規(guī)律遵循Navier-stokes方程，而在采空區(qū)O 形圈內(nèi)遵循Brinkman方程。同樣，瓦斯在采空區(qū)內(nèi)的速度場(chǎng)分布也符合Brinkman方程。風(fēng)流和瓦斯氣體的混合過(guò)程為對(duì)流擴(kuò)散和溶質(zhì)傳遞過(guò)程，運(yùn)用上述方程可以很好地刻畫工作面和采空區(qū)的流場(chǎng)分布規(guī)律。

（3）由于沙曲礦14205綜采工作面的特殊性，有必要采取拖管抽放措施。合理的抽放參數(shù)為抽放高度2.3m，拖管進(jìn)入采空區(qū)深度15m，抽放位置距離回風(fēng)巷外幫1.5m。

［1］ 吳兵，雷柏偉等.回采工作面上隅角瓦斯拖管抽采技術(shù)參數(shù)研究 ［J］.采礦與安全工程學(xué)報(bào)，2014（2）

［2］ 李東印，許燦榮等.采煤工作面瓦斯流動(dòng)模型及COMSOL數(shù)值解算 ［J］.煤炭學(xué)報(bào)，2012 （6）

［3］ 章爭(zhēng)榮，張湘?zhèn)?二維定常不可壓縮粘性流動(dòng)N-S方程的數(shù)值流形方法 ［J］.計(jì)算力學(xué)學(xué)報(bào)，2010（3）

［4］ 錢鳴高，許家林.覆巖采動(dòng)裂隙分布的 “O”形圈特征研究 ［J］.煤炭學(xué)報(bào)，1998 （5）

［5］ 張人偉，賀曉剛等.朱仙莊礦綜放面采空區(qū) “三帶”范圍的確定及應(yīng)用 ［J］ .采礦與安全工程學(xué)報(bào)，2008 （3）

［6］ 楊天鴻，陳仕闊等.采空垮落區(qū)瓦斯非線性滲流-擴(kuò)散模型及其求解 ［J］.煤炭學(xué)報(bào)，2009 （6）

［7］ Brinkman H C.A calculation of the viscous force exerted by flowing fluid on a dense swam of particles［J］.Applied Scientific Research，1947 （1）

［8］ 李宗翔，衣剛等.基于 “O”型冒落及耗氧非均勻采空區(qū)自燃分布特征 ［J］.煤炭學(xué)報(bào)，2012 （3）

［9］ 楊其鑾，王佑安.煤屑瓦斯擴(kuò)散理論及其應(yīng)用 ［J］.煤炭學(xué)報(bào)，1986 （3）

［10］ 李宗翔，紀(jì)書麗等.采空區(qū)瓦斯與大氣兩相混溶擴(kuò)散模型及其求解 ［J］.巖石力學(xué)與工程學(xué)報(bào)，2005（16）

［11］ 楊天鴻，徐濤等.應(yīng)力-損傷-滲流耦合模型及在深部煤層瓦斯卸壓實(shí)踐中的應(yīng)用 ［J］.巖石力學(xué)與工程學(xué)報(bào)，2005 （16）

［12］ 馬海峰，王磊等.掘進(jìn)煤巷瓦斯?jié)B流模型及運(yùn)移規(guī)律研究 ［J］.中國(guó)煤炭，2013 （8）

［13］ 李培超，孔祥言等.飽和多孔介質(zhì)流固耦合滲流的數(shù)學(xué)模型 ［J］ .水動(dòng)力學(xué)研究與進(jìn)展 （A 輯），2003 （4）

［14］ 尹光志，李銘輝等.基于含瓦斯煤巖固氣耦合模型的鉆孔抽采瓦斯三維數(shù)值模擬 ［J］ .煤炭學(xué)報(bào)，2013 （4）

［15］ 劉華鋒，陳輝等.順層鉆孔瓦斯抽采半徑確定及抽采效果考察研究 ［J］.中國(guó)煤炭，2013 （6）

［16］ 黃志安，張英華等.FLAC 在確定沙曲礦裂隙帶上下界中的應(yīng)用 ［J］.礦業(yè)研究與開發(fā)，2006 （1）

［17］ 王恩志，王洪濤等 .“以縫代井列”——排水孔幕模擬方法探討 ［J］.巖石力學(xué)與工程學(xué)報(bào)，2002（1）

［18］ 李宗翔，王繼仁等.采空區(qū)開區(qū)移動(dòng)瓦斯抽放的數(shù)值模似 ［J］.中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，2004 （1）

［19］ 郭玉森，林柏泉等.圍巖裂隙演化與采動(dòng)卸壓瓦斯儲(chǔ)運(yùn)的耦合關(guān)系 ［J］ .采礦與安全工程學(xué)報(bào)，2007 （4）

［20］ 邢云峰，孫曉元，孫英峰.合理鉆孔間距和瓦斯抽放時(shí)間的理論分析與數(shù)值模擬 ［J］ .中國(guó)煤炭，2014 （5）