張 璐，王 剛,2，于俊紅，肖智勇，王 珂

(1.山東科技大學(xué) 山東省土木工程防災(zāi)減災(zāi)重點實驗室，山東 青島 266590；2.山東科技大學(xué) 礦山災(zāi)害預(yù)防控制省部共建國家重點實驗室培育基地，山東 青島 266590；3.北京市勘察設(shè)計研究院有限公司，北京 100038)

隨著全球氣候變暖問題的日益突出，二氧化碳封存技術(shù)作為減少大氣溫室氣體含量的有效方法引起了人們的注意。煤巖是一種典型的雙重介質(zhì)，將二氧化碳注入煤層不但可以減少大氣中溫室氣體的含量，也可以驅(qū)替孔隙中的瓦斯和石油，強化油、氣的采收率[1-5]。煤巖滲透性是決定瓦斯產(chǎn)量和二氧化碳注氣量的關(guān)鍵因素[6-8]，因此建立煤巖滲透性演化模型具有重要的工程實用價值。

煤層中吸附性氣體(二氧化碳和甲烷)吸附/解吸引起的煤巖膨脹/收縮變形對孔隙度和滲透率的影響已經(jīng)得到廣泛研究。Sawyer 等[9]提出了煤巖基質(zhì)孔隙度與氣體壓力和濃度的函數(shù)關(guān)系模型。隨后，很多學(xué)者在此基礎(chǔ)上進行了大量研究。Seidle等[10]在不考慮有效應(yīng)力變化引起的固體介質(zhì)彈性變形的情況下，建立煤巖解吸收縮變形時的滲透率變化模型。潘榮瑩[11]考慮頁巖儲層在固體變形和多重流動機制耦合作用下，建立頁巖儲層單孔介質(zhì)流固耦合模型。上述研究均把煤巖作為一個整體進行滲透率研究。為了建立更準(zhǔn)確的煤巖滲透率模型，Elsworth等[12]將煤巖視為由基質(zhì)-裂隙組成的雙重介質(zhì)，建立非吸附性物質(zhì)-水滲流情況下的煤巖介質(zhì)變形模型。之后Wu等[13]在其基礎(chǔ)上，考慮煤巖解吸收縮變形，忽略氣體運移方式對滲透率的影響，建立了煤巖雙孔雙滲模型。但實驗室條件下難以對基質(zhì)和裂隙滲透率區(qū)分測量，為了彌補這方面的不足，一些學(xué)者[14-15]引入有限元軟件COMSOL，采用數(shù)值模擬手段，通過理論模型組成的方程組建模，控制多場耦合作用來研究裂隙和基質(zhì)的滲透率演化。然而這些研究大多都是在不改變有效應(yīng)力的情況下進行的，且很少考慮基質(zhì)和裂隙之間氣體擴散作用的影響，使得數(shù)值計算獲得的滲透率結(jié)果與實際情況存在一定的誤差。

吸附氣體在多孔介質(zhì)中的質(zhì)量運移方式為Fick擴散、Knudsen擴散和表面擴散。本研究將煤巖視為雙重介質(zhì)，考慮有效應(yīng)力引起的煤巖彈性變形和吸附/解吸效應(yīng)引起的基質(zhì)膨脹/收縮變形的影響，以及基質(zhì)和裂隙之間的氣體擴散，建立煤巖雙孔雙滲模型；并利用有限元軟件COMSOL，探究原位條件下二氧化碳封存過程中的煤巖滲透性演化規(guī)律。

1 滲透率演化模型

將煤巖視為由基質(zhì)-裂隙組成的雙重介質(zhì)，其物理模型如圖1所示。二氧化碳?xì)怏w在煤巖中以游離態(tài)和吸附態(tài)兩種形式存在。煤巖基質(zhì)中的孔隙大多為微米級和納米級，因此，傳統(tǒng)的達(dá)西定律無法準(zhǔn)確地描述氣體在多孔介質(zhì)中的質(zhì)量運移規(guī)律?；|(zhì)孔隙中氣體運移方程仍遵循Darcy公式，引入視滲透率模型對基質(zhì)滲透率進行修正，形成修正Darcy流動方程[11]。綜合考慮有效應(yīng)力變化引起的煤巖彈性變形、吸附/解吸動態(tài)過程引起的基質(zhì)膨脹/收縮變形，氣體Fick擴散、Knudsen擴散、表面擴散對氣體質(zhì)量運移的影響，建立煤巖滲透性演化模型。假設(shè)：

圖1 煤巖雙重孔隙介質(zhì)模型

1)儲層是等溫的，氣體黏度在等溫條件下保持恒定；

2)雙重孔隙介質(zhì)是均質(zhì)、各向同性的；

3)吸附態(tài)氣體僅存在于基質(zhì)內(nèi)，且吸附-解吸規(guī)律遵循Langmuir等溫吸附方程；

4)煤巖中僅存在二氧化碳?xì)怏w，且視二氧化碳為理想氣體。

1.1 視滲透率模型

單組分吸附性氣體在孔隙介質(zhì)中的質(zhì)量運移方式為Fick擴散、Knudsen擴散和表面擴散。氣體分子之間相互碰撞產(chǎn)生Fick擴散，氣體分子與壁面碰撞產(chǎn)生Knudsen擴散，吸附態(tài)氣體分子沿壁面蠕動產(chǎn)生表面擴散[16]，如圖2所示。氣體在多孔介質(zhì)中的Fick擴散可采用Javadpour模型，Civan模型和DGM模型描述。三種模型在相同溫度、相同等效流動半徑條件下結(jié)果相差不大[17]?？紤]吸附氣體分子所占空間對氣體流動影響，可用“壁聯(lián)擴散”表征Knudsen擴散和表面擴散的總效應(yīng)[18]。由于壁聯(lián)擴散和滑脫效應(yīng)等同，在流動計算中可互換，因此本研究基于DGM模型，考慮Knudsen擴散和表面擴散效應(yīng)，用Klinkenberg系數(shù)建立修正滲透率模型，即：

圖2 單組分吸附氣體在基質(zhì)中的運移機制

(1)

式中：下標(biāo)m表示基質(zhì)，下標(biāo)f表示裂隙；km表示煤巖基質(zhì)的固有滲透率；pm為基質(zhì)氣體壓力；ρga為標(biāo)準(zhǔn)狀況下的氣相密度；Cadmax為氣體最大吸附量濃度；μ為氣體黏度；Dad為表面擴散系數(shù)；pL為Langmuir壓力；bk為Klinkenberg系數(shù)，與固有滲透率和孔隙度的變化有關(guān)[17]。

(2)

其中，φm和φm0分別為基質(zhì)當(dāng)前孔隙度和初始孔隙度。

當(dāng)煤巖中的氣體遵循Langmuir吸附規(guī)律時，表明擴散系數(shù)Dad與孔隙中氣體的覆蓋度有關(guān)，則表面擴散系數(shù)與基質(zhì)氣體壓力的關(guān)系[11]為：

(3)

式中Dad0為氣體覆蓋度為0時的表面擴散系數(shù)。假設(shè)裂隙中不發(fā)生吸附解吸行為，因此也不存在表面擴散現(xiàn)象。Knudsen擴散與分子運動的自由程有關(guān)，當(dāng)孔隙尺寸很小，小到與分子運動自由程一個級別時，分子與壁面碰撞的幾率遠(yuǎn)高于分子之間的碰撞幾率，Knudsen擴散越明顯；反之，當(dāng)孔隙尺寸越大，大到遠(yuǎn)遠(yuǎn)高于分子運動自由程時，分子之間的碰撞幾率占主導(dǎo)，遠(yuǎn)高于分子與壁面碰撞的幾率，Knudsen擴散可以忽略[19]。考慮到煤巖裂隙開度的數(shù)量級遠(yuǎn)高于分子運動自由程，因此不考慮裂隙中的Knudsen擴散，即：

kf,app=kf。

(4)

1.2 固有滲透率模型

煤巖基質(zhì)和裂隙所受有效應(yīng)力[15]：

(5)

式中，α、γ分別為基質(zhì)和裂隙的有效應(yīng)力系數(shù)。對于雙重介質(zhì)模型，精確考慮有效應(yīng)力變化引起的基質(zhì)和裂隙的彈性變形以及吸附/解吸效應(yīng)引起的基質(zhì)煤巖的膨脹/收縮變形，給出煤巖介質(zhì)體應(yīng)變表達(dá)式[20]。

(6)

式中：Δεv為煤巖介質(zhì)體積應(yīng)變；Km、Kf分別為基質(zhì)、裂隙的體積模量。設(shè)體積單元邊長為s，基質(zhì)寬度為a，裂隙寬度為b，則s=a+b。Δσem、Δσef分別為基質(zhì)和裂隙的有效應(yīng)力變化量，Δεs為煤巖介質(zhì)的吸附應(yīng)變[20]：

(7)

之后省略體應(yīng)變與線應(yīng)變方程的2次項，簡化為3倍關(guān)系，獲得裂隙開度變化量

(8)

從而，裂隙度φf可表示為[15]：

(9)

另外，由于煤巖基質(zhì)孔隙度等于基質(zhì)孔隙體積與基質(zhì)體積的比值，即φm=Vmp/Vm，對等式兩邊同時求微分，利用基質(zhì)體積模量和基質(zhì)孔隙體積模量之間的關(guān)系Kmp=Kmφm/α可以獲得基質(zhì)孔隙度變化方程。

最后，基于固有滲透率與孔隙度的三次方關(guān)系，得到基質(zhì)和裂隙的固有滲透率演化模型[20]：

(10)

(11)

式中，基質(zhì)有效應(yīng)力系數(shù)α=1-K/Km。K表示煤巖介質(zhì)的體積模量。

1.3 氣體滲流控制方程

氣體在基質(zhì)和裂隙中運移需要滿足氣體質(zhì)量守恒方程[17]，即：

(12)

(13)

其中：ρs表示煤巖密度，VL表示Langmuir體積常數(shù)。ρmg表示基質(zhì)中游離的二氧化碳?xì)怏w密度，狀態(tài)方程為[19]：

(14)

其中：Mg為氣體相對分子質(zhì)量，R為氣體常數(shù)，T為溫度。經(jīng)典Darcy方程在基質(zhì)孔隙中不再適用，因此采用視滲透率在Fick擴散、Knudsen擴散、表面擴散作用下對滲透率進行修正，因此：

(15)

基質(zhì)與裂隙之間的氣體交換量可表示為：

(16)

將式(13)～(16)代入式(12)，得到氣體在基質(zhì)和裂隙中滲流控制方程：

(17)

(18)

1.4 煤巖變形控制方程

以位移為基本未知量，考慮地層壓力改變對煤層變形的影響，煤巖變形控制方程可表示為[21]：

(19)

其中：G為剪切模量，u為位移量，ν為泊松比，f為體積力。式(17)～(19)構(gòu)成滲流-應(yīng)力耦合作用下煤巖雙重介質(zhì)滲流和變形控制方程組。給定邊界條件及初始條件，求解該方程組可得到煤層滲透性演化規(guī)律。

2 模型驗證

為驗證理論模型(17)～(19)的有效性，將模型結(jié)果與文獻[22]實驗結(jié)果對比。實驗采用的煤巖來自懷俄明州吉列煤礦，采取長334 mm，直徑69.5 mm的圓柱形煤巖作為試樣。主要實驗儀器為：巖心夾持器、加壓氣瓶、氣瓶調(diào)節(jié)閥、計量閥、背壓閥、壓力傳感器和流量計。將干燥后的煤巖試樣置于夾持器中，給定6.885 MPa(1 000 psi)的圍壓，二氧化碳?xì)怏w注入壓力范圍為0.5～5.1 MPa。采用同一個煤巖試樣，連續(xù)進行5次實驗，獲得不同注氣壓力作用下煤巖的滲透率比率數(shù)據(jù)，具體實驗步驟參見文獻[22]?；谠嚇訋缀纬叽绾褪芰Φ膶ΨQ性，將試樣縱向半截面簡化為2D平面模型，如圖3所示，二氧化碳從左側(cè)注入，右側(cè)流出。模型左側(cè)水平方向位移和底部豎直方向位移均為0，實驗的煤巖試樣相關(guān)參數(shù)如表1所示。

圖3 實驗?zāi)Ｐ?/p>

表1 實驗條件下的相關(guān)參數(shù)

模型結(jié)果與實驗數(shù)據(jù)的對比如圖4所示，其中滲透率比值的方程為[23]：

圖4 實驗數(shù)據(jù)與模型結(jié)果對比圖

(20)

從圖4可以看出，在較低孔隙壓力條件下,匹配結(jié)果的誤差很小。在孔隙壓力為4.08和5.12 MPa時誤差分別為6.01%和6.62%，位于正常范圍之內(nèi)。產(chǎn)生這一現(xiàn)象的原因是：為了控制不變量，采用一個試樣連續(xù)進行實驗，試樣在每一次實驗后都會產(chǎn)生一定量不可逆變形，如軸向位移。在已經(jīng)進行三次實驗的條件下，試樣的初始參數(shù)變化使模擬對環(huán)向和軸向長度等參數(shù)高估或低估，因此與實驗結(jié)果產(chǎn)生了誤差。綜上表明，模型的模擬結(jié)果與實驗數(shù)據(jù)吻合良好，證明建立的理論模型有較強的適用性。

3 數(shù)值模擬

利用建立的理論模型(17)～(19)和有限元軟件COMSOL，對原位狀態(tài)下二氧化碳封存過程中煤層滲透性的演化規(guī)律進行數(shù)值模擬。COMSOL Multiphysics 是可進行多物理場耦合的有限元模擬軟件，利用固體力學(xué)模塊并引入PDE模塊，基于理論模型(17)～(19)模擬控制物理場，實現(xiàn)模型的求解。

數(shù)值計算模型的模擬區(qū)位于鄂爾多斯盆地榆林區(qū)山西組山2段，其地層區(qū)構(gòu)造穩(wěn)定、地層平緩，斷裂活動不發(fā)育，且煤層較厚，為儲層內(nèi)部實現(xiàn)二氧化碳封存創(chuàng)造了條件，部分地層剖面圖如圖5(a)所示。山2段地層厚40～65 m，地層壓力22.94～28.87 MPa，平均26.7 MPa，地層平均溫度186.8 K，平均孔隙度6.2%，滲透率0.1～1 mD[24]。取儲層厚度為50 m，定義400 m×400 m的注入面積，采用豎井注入，井口半徑為0.1 m，由于結(jié)構(gòu)的對稱性，只模擬注入面積的四分之一(如圖5(b))。模型網(wǎng)格劃分如圖5(c)所示。模型左邊的水平方向和底部的豎直方向位移為0，上部和右側(cè)施加27 MPa的邊界壓力，邊界定義為無流動，注氣壓力3.45 MPa，初始儲層壓力1.5 MPa，模擬天數(shù)為10 000 d，二氧化碳?xì)怏w的相關(guān)參數(shù)如表2。取點A(100,100)作為分析點以分析在整個模擬期間基質(zhì)和裂隙的相關(guān)參數(shù)的演化，輸出孔隙壓力、滲透率和吸附應(yīng)變的結(jié)果如圖6～9。

圖5 原位狀態(tài)地質(zhì)背景及模型圖

表2 模擬原位狀態(tài)下相關(guān)參數(shù)[25-26]

圖6 氣體壓力隨時間的變化

4 討論

4.1 氣體壓力的變化

圖6為開始注氣到第10 000 d時基質(zhì)和裂隙內(nèi)的氣體壓力變化過程。從圖中可以看出，隨著注氣時間的增加，基質(zhì)和裂隙的氣體壓力均成多項式增加，并隨時間的增長逐漸接近注氣壓力。

煤巖被視為雙重孔隙介質(zhì)結(jié)構(gòu)，裂隙可分別與基質(zhì)和外界氣體進行氣體交換，而基質(zhì)只與裂隙氣體有交換作用，與外界氣體沒有直接的氣體交換作用，即氣體可直接流入裂隙而無法直接進入基質(zhì)內(nèi)。圖6顯示裂隙滲透率存在一個明顯的拐點，將該點之前的部分稱為第一階段，之后的部分稱為第二階段。在第一階段，裂隙氣壓在開始時刻即明顯增加，且之后增長速度明顯快于基質(zhì)氣壓。產(chǎn)生這一現(xiàn)象的原因是：開始注氣時，二氧化碳?xì)怏w立刻進入裂隙，裂隙氣體壓力迅速增加，但氣體無法直接進入基質(zhì)，只能通過基質(zhì)與裂隙之間的擴散作用，從裂隙進入到基質(zhì)，這是一個緩慢的過程，因此在這一階段基質(zhì)的氣體壓力只能隨著時間的推移緩慢增加。在第二階段，裂隙氣壓的增長速率變緩，出現(xiàn)這一現(xiàn)象的原因是：裂隙氣壓增大，與注氣壓力之間的壓力差降低，氣體進入裂隙的速度變緩，因此注氣量逐漸減少，裂隙氣體壓力增長幅度變緩。

4.2 滲透率的變化

隨著基質(zhì)和裂隙氣體壓力的變化，滲透率也發(fā)生相應(yīng)的改變。如圖7～8表示滲透率隨時間和氣體壓力的增加而降低的一個過程。

圖7 滲透率隨時間的變化

煤巖介質(zhì)中同時存在兩種機制：一種是在注氣過程中，隨著氣體壓力增加，有效應(yīng)力降低，使孔隙體積發(fā)生變化，導(dǎo)致滲透率改變；另一種是吸附作用，使得儲層固體組分發(fā)生膨脹，導(dǎo)致孔隙體積改變，滲透率發(fā)生變化。圖7所示基質(zhì)和裂隙的滲透率隨時間的推移而降低的過程，且裂隙滲透率降低的幅度大于基質(zhì)，滲透率下降到一定值時儲層內(nèi)形成一個動態(tài)平衡，使?jié)B透率維持在一個恒定的數(shù)值。為了探究壓力對滲透率的影響，氣體壓力的變化與基質(zhì)和裂隙滲透率之間的關(guān)系如圖8所示，基質(zhì)滲透率隨著氣體壓力的增加而線性下降，說明兩種機制對滲透率產(chǎn)生的凈影響與壓力成負(fù)相關(guān)。而裂隙滲透率在氣體壓力達(dá)到2.5 MPa前下降緩慢，之后才開始迅速下降，且下降速度大于基質(zhì)，結(jié)合圖6～7發(fā)現(xiàn)，造成這一現(xiàn)象的原因是注氣前期裂隙的滲透率較高，由于裂隙與注氣壓力之間的壓差較大，氣體進入裂隙的速度很快，因此該過程經(jīng)歷的時間很短。而吸附作用的發(fā)生是一個緩慢的過程，短時間內(nèi)基質(zhì)膨脹量很小，孔隙體積變化很小，因此在氣體壓力達(dá)到2.5 MPa之前，裂隙滲透率下降很少。之后由于裂隙內(nèi)氣體壓力與注氣壓力之間的壓差減小，氣體進入裂隙的速度減緩，裂隙壓力從2.5 MPa增加到3.45 MPa這一過程經(jīng)歷的時間較長，吸附作用發(fā)揮作用明顯，孔隙減小，因此這一階段滲透率下降的趨勢明顯。

圖8 滲透率隨壓力的變化

4.3 吸附及有效應(yīng)力作用

為了進一步探究導(dǎo)致滲透率下降的關(guān)鍵因素，將影響基質(zhì)和裂隙滲透率的吸附作用和有效應(yīng)力作用分離開。將僅考慮有效應(yīng)力作用引起的孔隙變形、僅考慮吸附作用引起的孔隙變形及同時考慮兩種作用引起的孔隙變形三種情況下的滲透率變化情況在圖9中進行對比。

圖9 僅有效應(yīng)力作用，僅吸附作用和綜合作用三種情況下的滲透率變化

從圖9(a)中可以看出，基質(zhì)中僅考慮吸附作用和同時考慮兩種作用時的滲透率曲線幾乎重合，而僅考慮有效應(yīng)力作用時的滲透率幾乎沒有變化，表明基質(zhì)滲透率由吸附作用控制，有效應(yīng)力變化產(chǎn)生的影響很小。由于吸附作用的產(chǎn)生，二氧化碳吸附于孔隙壁面，使得固體基質(zhì)發(fā)生膨脹，孔隙體積減少，基質(zhì)滲透率降低。同時，從圖9(b)中可以看出，裂隙中僅考慮吸附作用與同時考慮兩種作用下的滲透率比值隨著裂隙氣體壓力的增大，差異逐漸顯現(xiàn)。另外僅考慮有效應(yīng)力作用時的滲透率隨著裂隙氣體壓力的增大而增大，表明裂隙壓力越大，有效應(yīng)力作用越明顯。因此在這一過程中，裂隙滲透率由有效應(yīng)力和吸附作用共同影響，二者為競爭關(guān)系。結(jié)合圖8～9，可以發(fā)現(xiàn)：基質(zhì)中隨著吸附作用產(chǎn)生，固體組分膨脹，孔隙體積減小，滲透率降低，有效應(yīng)力作用對孔隙變化產(chǎn)生的影響很小，可以忽略，因此凈滲透率為下降趨勢。而裂隙中存在兩個相互競爭的機制：隨著吸附作用的產(chǎn)生，基質(zhì)固體組分發(fā)生膨脹，導(dǎo)致裂隙開度減小，滲透率減小；隨著氣體壓力的增加，有效應(yīng)力減小，裂隙體積增大，孔隙度增加，滲透率升高。裂隙滲透率在有效應(yīng)力和吸附機制的共同作用下，凈滲透率呈下降趨勢，可以認(rèn)為裂隙滲透率變化過程中吸附作用占主導(dǎo)。為了突出吸附作用的影響過程，吸附作用導(dǎo)致的應(yīng)變?nèi)鐖D10所示。

從圖10可以看出，隨著時間的推移，吸附產(chǎn)生的應(yīng)變逐漸增加，曲線增長與圖7對滲透率的影響一致，與圖9的結(jié)果相呼應(yīng)。最后系統(tǒng)到達(dá)動態(tài)平衡狀態(tài)時，吸附應(yīng)變也不再發(fā)生變化。圖10可以從側(cè)面證明滲透率演化趨勢的正確性。

圖10 注氣過程中的吸附應(yīng)變

4.4 注氣量預(yù)測

結(jié)合雙重介質(zhì)滲透率演化模型和數(shù)值模擬，可以獲得注氣10 000 d的注氣量預(yù)測如圖11，對其進行曲線擬合，獲得時間與注氣量的擬合表達(dá)式為：

圖11 注氣量的預(yù)測

y=7×10-19x5-2×10-14x4+3×10-10x3-2×10-6x2+0.007 4x+0.096 6，

R2=0.999 9。

式中，x表示注氣天數(shù)，y表示注氣量，單位為106m3，R2為相關(guān)系數(shù)。從圖11可以看出，現(xiàn)場原位二氧化碳封存在前3 000 d左右發(fā)生速度較快，之后速度開始降低，第6 000 d之后，注氣發(fā)生緩慢。因基質(zhì)滲透率比裂隙滲透率高兩個數(shù)量級，裂隙滲透率的變化對注氣量的影響相對較大。裂隙滲透率前期下降緩慢，注氣量增加較快；后期裂隙滲透率快速下降，注氣量增長速度變緩。

5 結(jié)論

將煤巖視為雙重孔隙吸附介質(zhì)，在有效應(yīng)力引起的彈性變形和吸附作用引起的膨脹變形作用下，綜合考慮Fick擴散、Knudsen擴散和表面擴散，建立改進的滲透率演化模型，并用實驗數(shù)據(jù)驗證該模型的有效性，最后將該模型用于二氧化碳封存技術(shù)，探究現(xiàn)場原位二氧化碳封存10 000 d過程中的煤巖滲透率演化規(guī)律及相關(guān)參數(shù)響應(yīng)。結(jié)合理論與數(shù)值模擬，得出以下結(jié)論：

1)注氣過程中，隨著氣體壓力的增大，裂隙滲透率的演化過程更為復(fù)雜，因此準(zhǔn)確預(yù)測裂隙滲透率的變化過程極為重要。

2)在較低的注氣壓力作用下，吸附作用比有效應(yīng)力作用對滲透率的影響更為明顯，但仍需考慮有效應(yīng)力作用，尤其是裂隙系統(tǒng)。若不考慮有效應(yīng)力的作用會使預(yù)測滲透率的結(jié)果產(chǎn)生較大誤差。

3)結(jié)合吸附作用、有效應(yīng)力作用和兩者同時作用的對比，可以得出基質(zhì)滲透率的下降是由吸附作用導(dǎo)致的固體組分膨脹、隙收縮引起的；而裂隙的凈滲透率變化是由有效應(yīng)力引起的彈性變形和吸附作用引起的基質(zhì)膨脹變形共同導(dǎo)致的，兩者為相互競爭的機制，其中吸附作用占主導(dǎo)。