孫曉元 王 川 孫英峰 邢云峰

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

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

1 風流在采場中的運動模型

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

1.1 風流在自由空間中的流動方程

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

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

ρ——風流密度，kg／m3；

p——風流壓力，Pa；

ν——風流的運動粘度，m2／s；

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

式中：η——風流的動力粘度，Pa·s。

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

1.2 風流在多孔介質(zhì)中的流動方程

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

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

I——單位矢量；

F——流體阻力，N。

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

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

1.3 風流運動模型

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

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

圖1 風流運動模型

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

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

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

L——工作面長度，m。

2 瓦斯在采場中的運移模型

上述討論研究了風流在采場中的運動規(guī)律，然而往往更要關注的是瓦斯在采場中的運移特征。為此，將瓦斯在采空區(qū)的運移分為濃度擴散和對流運移兩個過程來進行討論，這是因為若單純考慮Fick擴散第二定律來描述平面二維尺度空間內(nèi)的瓦斯非穩(wěn)定流動場時，其方程表達式為：

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

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

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

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

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

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

3 瓦斯在抽放滲流場中的運移模型

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

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

kg／t；

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

顯然，單位體積煤巖體中的瓦斯含量應包括吸附瓦斯含量和游離瓦斯含量兩部分。其中吸附瓦斯含量由Langmiur等溫吸附方程求得：

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

p1——瓦斯壓力，MPa；

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

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

ρ0——標準大氣壓下的瓦斯密度，kg／m3。

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

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

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

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

顯然，存在公式：

假設瓦斯在煤層中的滲流模式為達西流動，并考慮Klinkenberg效應，則滲流的速度為：

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

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

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

?p——壓力梯度；

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

式中：ζ——含瓦斯煤巖的體積應變。

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

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

通過上述分析，求得回采工作面風流與瓦斯的流場分布及抽放滲流場控制方程的解析解是十分困難的，這是因為這些方程絕大部分都含有一階或二階偏微分項。因此，可以專門采用基于偏微分方程的有限元求解平臺COMSOL Multiphysics多物理場耦合分析工具。

4.1 幾何模型

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

4.2 初始值和邊界條件

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

圖2 風流在采空區(qū)和工作面運動模型的濃度分布云圖

4.3 模擬結(jié)果分析

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

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

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

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

圖4 拖管抽放前后等值面與瓦斯體積分數(shù)對比圖

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

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

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

5 結(jié)論

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

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

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

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

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

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

［4］ 錢鳴高，許家林.覆巖采動裂隙分布的 “O”形圈特征研究 ［J］.煤炭學報，1998 （5）

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

［6］ 楊天鴻，陳仕闊等.采空垮落區(qū)瓦斯非線性滲流-擴散模型及其求解 ［J］.煤炭學報，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］.煤炭學報，2012 （3）

［9］ 楊其鑾，王佑安.煤屑瓦斯擴散理論及其應用 ［J］.煤炭學報，1986 （3）

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

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

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

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

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

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

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

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

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

［19］ 郭玉森，林柏泉等.圍巖裂隙演化與采動卸壓瓦斯儲運的耦合關系 ［J］ .采礦與安全工程學報，2007 （4）

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