煤粒瓦斯定壓吸附數(shù)學(xué)模型及數(shù)值解算

徐浩，秦躍平，毋凡，劉佳，褚翔宇，劉曉薇

中國(guó)礦業(yè)大學(xué)(北京) 應(yīng)急管理與安全工程學(xué)院，北京 100083

掌握瓦斯在煤體中的流動(dòng)規(guī)律對(duì)預(yù)測(cè)煤層氣抽采產(chǎn)量、分析儲(chǔ)層成藏機(jī)理和優(yōu)化生產(chǎn)策略具有重要意義[1-2]。目前，針對(duì)煤基質(zhì)中瓦斯吸附解吸的流動(dòng)機(jī)理并沒有統(tǒng)一和清楚的認(rèn)識(shí)，學(xué)術(shù)界也一直對(duì)此展開了深入的研究和探討[3-5]。

有學(xué)者認(rèn)為瓦斯在煤基質(zhì)中吸附解吸是受濃度差驅(qū)動(dòng)的，可以用菲克擴(kuò)散理論來描述。但是菲克定律中常擴(kuò)散系數(shù)只能描述剛開始較短時(shí)間甲烷的擴(kuò)散過程，其理論預(yù)測(cè)結(jié)果與整個(gè)解吸時(shí)間尺度上的實(shí)驗(yàn)值存在顯著的偏差[6-9]。后續(xù)有研究人員提出了基于菲克定律的動(dòng)擴(kuò)散系數(shù)數(shù)學(xué)模型[10-15]，并驗(yàn)證了該模型可以表征氣體在煤粒中的運(yùn)移全過程。雖然將擴(kuò)散系數(shù)變成與時(shí)間有關(guān)的函數(shù)表達(dá)式能夠保持模擬和實(shí)驗(yàn)結(jié)果的一致性，但是它卻違背了菲克定律最初設(shè)定的常擴(kuò)散系數(shù)的假設(shè)條件，這可能會(huì)造成所建模型的物理意義與理論基礎(chǔ)不明晰。另外，動(dòng)擴(kuò)散系數(shù)存在一些待定參數(shù)，會(huì)導(dǎo)致計(jì)算過程復(fù)雜化。文獻(xiàn)[7-9，16]基于菲克和達(dá)西定律建立了在變壓、定壓條件下的煤粒瓦斯吸附和解吸數(shù)學(xué)模型，并通過大量的實(shí)驗(yàn)結(jié)果驗(yàn)證了數(shù)值解算結(jié)果，發(fā)現(xiàn)煤粒中瓦斯的流動(dòng)是符合達(dá)西定律而不是菲克定律。達(dá)西定律認(rèn)為瓦斯在煤粒中的流速(體積流量)與壓力梯度成正比，它是基于較大尺寸孔隙內(nèi)的流動(dòng)。然而煤粒中存在大量的微孔隙，學(xué)者們大都認(rèn)為達(dá)西定律在微孔隙中是不成立的。因此，將達(dá)西定律用來描述瓦斯在煤粒中的流動(dòng)過程也受到一些爭(zhēng)議。文獻(xiàn)[8-9]指出煤粒的微孔隙中存在大量的吸附態(tài)瓦斯，煤粒中的瓦斯流動(dòng)主要是游離瓦斯的運(yùn)動(dòng)。煤粒中的瓦斯無論是何種狀態(tài)，其本質(zhì)均可以理解為傳質(zhì)。

基于此，本文首先設(shè)計(jì)了定壓情況下的煤粒瓦斯吸附實(shí)驗(yàn)。提出了瓦斯在煤粒中的流動(dòng)是受游離瓦斯密度梯度驅(qū)動(dòng)的理論模型，即瓦斯質(zhì)量流量與游離瓦斯密度梯度成正比。之后建立了以新理論模型為基礎(chǔ)的煤粒瓦斯定壓吸附數(shù)學(xué)模型，并重新構(gòu)建了基于達(dá)西定律的瓦斯吸附數(shù)學(xué)模型。通過將兩種模擬解算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，直觀判斷所提出的理論模型是否與實(shí)驗(yàn)相符。詳細(xì)討論了理論模型中的關(guān)鍵比例系數(shù)與時(shí)間和壓力的關(guān)系，旨在進(jìn)一步了解煤粒中瓦斯流動(dòng)的基本規(guī)律。

1 煤粒瓦斯定壓吸附實(shí)驗(yàn)

1.1 煤樣制備

煤樣取自山西水峪煤礦。將從煤礦井下現(xiàn)場(chǎng)采集的煤樣密封保存，運(yùn)回實(shí)驗(yàn)室后將其置于真空烘箱中，在105 ℃下烘干2 h。之后取出煤樣放入干燥器中進(jìn)行冷卻達(dá)到室溫，再將煤樣放入破碎機(jī)中破碎，并用實(shí)驗(yàn)篩進(jìn)行篩分。根據(jù)本實(shí)驗(yàn)要求，最終得到的煤粒為60～80目，粒徑范圍為180～250 μm。實(shí)驗(yàn)前還需要把煤粒放入真空干燥箱中，在105 ℃的溫度下干燥8 h以除去煤粒中的水分。

每組實(shí)驗(yàn)選取6.5 g的煤粒來測(cè)定干煤樣對(duì)瓦斯氣體的吸附情況。

1.2 實(shí)驗(yàn)系統(tǒng)

實(shí)驗(yàn)儀器采用高溫高壓氣體吸附分析儀H-Sorb 2006，實(shí)驗(yàn)系統(tǒng)主要包括溫度控制系統(tǒng)、等溫吸附系統(tǒng)和數(shù)據(jù)采集系統(tǒng)等三部分，如圖1所示。每次實(shí)驗(yàn)均將溫度控制在35 ℃，以保證是在恒溫條件下開展的煤粒瓦斯吸附實(shí)驗(yàn)。

圖1 實(shí)驗(yàn)系統(tǒng)

1.3 實(shí)驗(yàn)步驟

實(shí)驗(yàn)之前檢測(cè)系統(tǒng)的氣密性是否良好，必須在氣密性良好的情況下開展實(shí)驗(yàn)，而且每組實(shí)驗(yàn)均需重新檢測(cè)氣密性以保證結(jié)果的準(zhǔn)確性。利用真空充氦的方法測(cè)定樣品罐和參考罐的實(shí)際體積，并計(jì)算出樣品罐的自由空間體積。用真空泵對(duì)系統(tǒng)進(jìn)行抽真空處理，之后向參考罐中通入瓦斯，待壓力平衡后連接參考罐和樣品罐。使樣品罐中的煤粒在初始?jí)毫?.5 MPa、1 MPa、2 MPa、4 MPa的4種情況下吸附瓦斯。在這過程中樣品罐的壓力會(huì)不斷降低。為保證瓦斯在定壓環(huán)境下進(jìn)行吸附實(shí)驗(yàn)，需要不斷地貫通樣品罐和參考罐，使樣品罐始終保持在設(shè)置的恒定壓力值上，當(dāng)吸附解吸平衡時(shí)實(shí)驗(yàn)停止；工控機(jī)的數(shù)據(jù)采集系統(tǒng)可以記錄每秒樣品罐壓力的具體數(shù)值，并通過前后時(shí)間點(diǎn)的壓差來計(jì)算樣品罐中的煤粒累計(jì)瓦斯吸附量[16]。

2 達(dá)西數(shù)值模型

2.1 數(shù)學(xué)模型

文獻(xiàn)[16]的研究結(jié)果表明，定壓邊界條件下煤粒中的瓦斯解吸吸附過程符合達(dá)西定律(壓力梯度驅(qū)動(dòng))。本文在此基礎(chǔ)上建立了球坐標(biāo)下的煤粒瓦斯定壓吸附的流動(dòng)方程：

(1)

式中，a，b為朗格繆爾吸附常數(shù)，m3/t，1/MPa；ρs為煤視密度，t/m3；P為瓦斯壓力p的平方，MPa2；t為吸附時(shí)間，s；A為與游離瓦斯含量有關(guān)的系數(shù)，m3/(t·MPa)；f為孔隙率；λ為瓦斯透氣性系數(shù)，m2/(MPa2·d)；r為煤粒中心到其他位置的距離，m。

系數(shù)A的計(jì)算公式[17]為

(2)

式中，Tn為標(biāo)準(zhǔn)情況下的溫度，K；pn為標(biāo)準(zhǔn)情況下的壓力，取0.101 325 MPa；T為實(shí)驗(yàn)溫度，K。

流動(dòng)方程的初始和邊界條件為

(3)

式中，p0為吸附前的煤粒中的初始?jí)毫?，MPa；pw為煤粒外表面的壓力，MPa；P0為初始?jí)毫Φ钠椒?，MPa2；Pw為煤粒外表面壓力的平方，MPa2；R為煤粒半徑，m。

煤粒吸附實(shí)驗(yàn)前對(duì)煤粒脫氣處理，所以初始時(shí)刻煤粒中的瓦斯壓力約為0。與煤粒瓦斯變壓吸附過程不同的是，在定壓吸附實(shí)驗(yàn)過程中，煤粒外部空間的壓力始終保持設(shè)定值。因此，煤粒外表面壓力的平方Pw是不變的。

2.2 有限差分?jǐn)?shù)值解算

相對(duì)于解析解，數(shù)值解可以不用設(shè)定諸多的假設(shè)，準(zhǔn)確性也更高[18]。本文使用有限差分的數(shù)值方法求解2.1節(jié)中的數(shù)學(xué)方程[7-9]。將球形煤粒沿球的半徑從球心到球表面劃分為N個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)間距等比變小，編號(hào)為0，1，2，3，…，N，如圖2所示。

圖2 球形煤粒節(jié)點(diǎn)劃分

節(jié)點(diǎn)所在的球面用圖2中的洋紅色實(shí)線表示。以兩個(gè)相鄰節(jié)點(diǎn)間的中心作同心球面，得到綠色虛線所示的球面。相鄰虛線球面之間形成球殼，而在中心處形成一個(gè)實(shí)心球體，每一個(gè)球殼或小球包含一個(gè)節(jié)點(diǎn)。這樣可以得到以0點(diǎn)為中心的實(shí)心球和包含各節(jié)點(diǎn)的N個(gè)球殼。

根據(jù)質(zhì)量守恒定律，節(jié)點(diǎn)1到N-1所形成的各自的球殼的非穩(wěn)態(tài)瓦斯流動(dòng)差分方程為

(4)

(5)

式中，i，i-1，i+1為劃分節(jié)點(diǎn)編號(hào)；j為時(shí)間節(jié)點(diǎn)編號(hào)；Δtj為第j個(gè)時(shí)間步長(zhǎng)。

節(jié)點(diǎn)0處的非穩(wěn)態(tài)瓦斯流動(dòng)的差分方程為

(6)

(7)

節(jié)點(diǎn)N處的非穩(wěn)態(tài)瓦斯流動(dòng)的差分方程為

(8)

由式(4)至式(8)構(gòu)成了第j時(shí)刻以N個(gè)節(jié)點(diǎn)瓦斯壓力為未知量的完備方程組。已知一個(gè)時(shí)間步長(zhǎng)以及上一時(shí)刻的壓力，可以計(jì)算下一時(shí)刻的壓力。在此基礎(chǔ)上，可以得到任意時(shí)間煤粒內(nèi)的瓦斯壓力。計(jì)算的步長(zhǎng)采用等比步長(zhǎng)，這樣可以在保障精確度的情況下盡量節(jié)省計(jì)算時(shí)間。式(5)和式(7)的右邊為第j時(shí)刻節(jié)點(diǎn)瓦斯壓力的非線性表達(dá)式，需采用高斯迭代的方法進(jìn)行求解。具體參照文獻(xiàn)[7-9]中的解算步驟。

依據(jù)N節(jié)點(diǎn)和N-1節(jié)點(diǎn)的壓力值，計(jì)算j時(shí)刻單位質(zhì)量煤粒累計(jì)吸附的瓦斯體積含量：

(9)

式中，QV為單位質(zhì)量煤粒累計(jì)吸附瓦斯體積，cm3/g。

3 游離瓦斯密度梯度驅(qū)動(dòng)的擴(kuò)散模型

3.1 數(shù)學(xué)模型

煤粒中吸附態(tài)的瓦斯數(shù)量與煤粒表面積有關(guān)；游離態(tài)的瓦斯分子能自由移動(dòng)，其數(shù)量與空間有關(guān)。原始煤體中瓦斯的游離態(tài)與吸附態(tài)是處在動(dòng)態(tài)平衡的環(huán)境中，一旦外部卸壓，游離態(tài)的瓦斯便向外運(yùn)移，而氣態(tài)瓦斯的擴(kuò)散與其密度梯度成正比[8]。在此基礎(chǔ)上，提出了煤基質(zhì)中瓦斯質(zhì)量流量與游離瓦斯密度梯度成正比的新模型，其比例系數(shù)為瓦斯微孔道擴(kuò)散系數(shù)。數(shù)學(xué)表達(dá)式如下：

(10)

式中，Jm為煤基質(zhì)中瓦斯質(zhì)量流量，即單位時(shí)間內(nèi)通過單位面積的瓦斯質(zhì)量，g/(m2·s)；Dm為微孔道游離瓦斯擴(kuò)散系數(shù)，m2/s；ρg為游離態(tài)瓦斯密度，g/m3；l為等密度線外法線方向的長(zhǎng)度，m。

游離瓦斯可視為理想氣體，根據(jù)理想氣體狀態(tài)方程，則有

(11)

式中，p指瓦斯氣體壓力，MPa；Rs為通用氣體常數(shù)，8.314 J/(mol·K)；M為瓦斯的摩爾質(zhì)量，16 g/mol；Km為煤粒的微孔道擴(kuò)散系數(shù)，g/(MPa·m·s)。

達(dá)西模型中的流動(dòng)方程表示的是體積流量，而游離瓦斯密度梯度模型描述的是質(zhì)量流量，則新模型下的煤粒瓦斯沿徑向流動(dòng)的連續(xù)方程為

(12)

式中，ρc為標(biāo)準(zhǔn)情況下的瓦斯密度，取7.17×102g/m3。

初始和邊界條件為

(13)

與定壓吸附實(shí)驗(yàn)條件保持一致，煤粒內(nèi)部的初始?jí)毫0為0，外表面的壓力pw始終不變。

3.2 有限差分?jǐn)?shù)值解算

節(jié)點(diǎn)還是按照?qǐng)D2劃分，節(jié)點(diǎn)1到N-1所形成的各自球殼的非穩(wěn)態(tài)瓦斯流動(dòng)的差分方程為

(14)

節(jié)點(diǎn)0所形成的各自球殼的非穩(wěn)態(tài)瓦斯流動(dòng)差分方程為

(15)

節(jié)點(diǎn)N所形成的各自球殼的非穩(wěn)態(tài)瓦斯流動(dòng)差分方程為

(16)

式(14)至式(16)構(gòu)成了第j時(shí)刻以N個(gè)節(jié)點(diǎn)瓦斯壓力為未知量的完備方程組。同樣需要采用高斯迭代的方法，求解式(14)和式(15)中關(guān)于瓦斯壓力的非線性表達(dá)式。基于Visual Basic程序編寫相應(yīng)的解算代碼，程序結(jié)構(gòu)流程如圖3所示。

圖3 程序結(jié)構(gòu)流程

依據(jù)N節(jié)點(diǎn)和N-1節(jié)點(diǎn)的壓力值，計(jì)算j時(shí)刻單位質(zhì)量煤粒累計(jì)吸附瓦斯的質(zhì)量：

(17)

式中，Qm為單位質(zhì)量煤粒累計(jì)吸附瓦斯質(zhì)量，g/g。

4 實(shí)驗(yàn)與模擬結(jié)果的對(duì)比

為了便于比較兩種理論模型解算的結(jié)果，這里將基于游離瓦斯密度梯度理論模型得到的煤粒累計(jì)吸附瓦斯質(zhì)量含量Qm轉(zhuǎn)換為吸附體積含量QV。二者的轉(zhuǎn)換關(guān)系為

(18)

數(shù)值模擬中參數(shù)的選取見表1。其中，a、b、ρs、f等均由實(shí)驗(yàn)手段測(cè)取，λ和Km是在解算程序中不斷調(diào)試確定的。實(shí)驗(yàn)和模擬得到的是累計(jì)瓦斯吸附體積含量QV與時(shí)間t的變化關(guān)系。由于程序解算中的時(shí)間步長(zhǎng)是采用等比方式，導(dǎo)致時(shí)間數(shù)據(jù)的波動(dòng)情況太大，因此取時(shí)間(h)的對(duì)數(shù)值為橫坐標(biāo)，吸附量為縱坐標(biāo)，分別繪制了不同初始?jí)毫η闆r下的曲線趨勢(shì)，如圖4所示。

表1 模擬參數(shù)取值

由圖4可以看出，在不同初始?jí)毫η闆r下，基于達(dá)西理論和游離瓦斯密度梯度理論模型得到的煤粒瓦斯吸附模擬結(jié)果的變化趨勢(shì)，均與實(shí)驗(yàn)數(shù)據(jù)保持一致。由于數(shù)值模擬理想化的假設(shè)條件，導(dǎo)致一些模擬數(shù)值點(diǎn)和實(shí)驗(yàn)點(diǎn)存在些許誤差，但從整體上來看也是可以接受的。這在一定程度上說明新提出來的理論模型能夠用來描述煤粒瓦斯的流動(dòng)過程。另外，通過模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行匹配，可反演出兩種關(guān)鍵系數(shù)(λ和Km)的具體值，它們隨壓力的變化關(guān)系如圖5所示。

圖4 模擬與實(shí)驗(yàn)結(jié)果對(duì)比

圖5 關(guān)鍵系數(shù)隨壓力變化趨勢(shì)

由圖4可以看出，λ和Km均不隨時(shí)間的變化而改變；圖5則顯示出λ隨初始?jí)毫Φ脑龃蠖鴾p小，且變化幅度較大，從0.5 MPa到4 MPa透氣性系數(shù)下降了將近8倍；而Km在不同壓力情況下幾乎保持不變，微孔道擴(kuò)散系數(shù)受壓力影響很小。眾所周知，從理論上計(jì)算瓦斯抽采量時(shí)，需要清楚地認(rèn)識(shí)到瓦斯在煤基質(zhì)中的流動(dòng)規(guī)律。要保證計(jì)算結(jié)果的準(zhǔn)確性，達(dá)西模型中的透氣性系數(shù)必須隨壓力的增大而減小，而本文新提出的游離瓦斯密度梯度理論模型中的微孔道擴(kuò)散系數(shù)則不受此限制。通常，開采煤層中各處的瓦斯壓力分布不均勻，且受擾動(dòng)的影響不斷變化，此種情況下新假說的優(yōu)勢(shì)顯而易見。

5 討 論

學(xué)術(shù)界針對(duì)瓦斯在煤粒中的運(yùn)輸機(jī)制還沒有統(tǒng)一的認(rèn)識(shí)，大都認(rèn)為煤粒中的瓦斯流動(dòng)符合菲克擴(kuò)散定律(濃度梯度驅(qū)動(dòng))[19-21]。本課題組也對(duì)此作了深入研究，發(fā)現(xiàn)基于菲克定律的瓦斯流動(dòng)模型所預(yù)測(cè)的結(jié)果與實(shí)驗(yàn)值存在明顯的偏差，只有當(dāng)菲克定律的擴(kuò)散系數(shù)隨時(shí)間不斷減小時(shí)才能與實(shí)驗(yàn)結(jié)果相吻合[7-9，22-23]。之后經(jīng)過變壓、定壓下的煤粒瓦斯吸附和解吸等多種條件下的大量實(shí)驗(yàn)證明：基于達(dá)西定律(壓力梯度驅(qū)動(dòng))建立數(shù)值模型的模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相一致，其中的透氣性系數(shù)是不隨時(shí)間變化的，因此認(rèn)為瓦斯在煤粒中的流動(dòng)機(jī)理服從達(dá)西定律而不是菲克定律。

然而后續(xù)的研究中卻發(fā)現(xiàn)，在不同初始?jí)毫l件下，根據(jù)實(shí)驗(yàn)數(shù)據(jù)反演出的煤粒透氣性系數(shù)卻出現(xiàn)很大的差別，透氣性系數(shù)與試驗(yàn)初始?jí)毫瞥煞幢汝P(guān)系[23]。但實(shí)際上，如果菲克和達(dá)西理論模型是正確的，則其中的關(guān)鍵比例系數(shù)不應(yīng)隨時(shí)間或壓力的變化而變化。故而無論是菲克定律還是達(dá)西定律都不能很好地描述煤粒中的瓦斯流動(dòng)過程。根據(jù)達(dá)西定律只能推出瓦斯的體積流量與壓力梯度成正比，而不是質(zhì)量流量(單位時(shí)間單位面積上通過的瓦斯質(zhì)量)與壓力梯度成正比，這可能是在不同壓力條件下透氣性系數(shù)出現(xiàn)巨大差異的根本原因。

本文提出的煤粒中瓦斯質(zhì)量流量與游離瓦斯密度梯度成正比的新理論模型得到了較好的驗(yàn)證。通過將模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)匹配可知，其中的關(guān)鍵比例系數(shù)(微孔道擴(kuò)散系數(shù))與壓力、時(shí)間等參數(shù)無關(guān)。

6 結(jié)論與展望

本文開展了定壓條件下的煤粒瓦斯吸附實(shí)驗(yàn)，提出了游離瓦斯密度梯度驅(qū)動(dòng)的煤粒中瓦斯流動(dòng)模型，并與達(dá)西模型進(jìn)行對(duì)比。主要結(jié)論如下：

(1) 采用有限差分法解算兩種數(shù)學(xué)模型的模擬結(jié)果均與實(shí)驗(yàn)數(shù)據(jù)相一致。另外微孔道擴(kuò)散系數(shù)與壓力無關(guān)，透氣性系數(shù)卻隨壓力的增大而減小。

(2) 瓦斯在煤粒中的吸附運(yùn)移規(guī)律可以用游離瓦斯密度梯度驅(qū)動(dòng)的擴(kuò)散模型來表述，而且比壓力梯度驅(qū)動(dòng)的達(dá)西定律以及濃度梯度驅(qū)動(dòng)的菲克定律更加合理。

作為一種初步探討，本文僅實(shí)施了定壓條件下的煤粒瓦斯吸附實(shí)驗(yàn)來驗(yàn)證新提出的游離瓦斯密度梯度驅(qū)動(dòng)瓦斯擴(kuò)散模型。實(shí)驗(yàn)樣本相對(duì)較少，今后需要進(jìn)行全面系統(tǒng)研究，對(duì)該模型進(jìn)行廣泛驗(yàn)證。