煤層氣藏全流固耦合數(shù)學(xué)模型

未志杰 ，康曉東 ，劉玉洋 ，曾 楊

（1.海洋石油高效開發(fā)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100028；2.中海油研究總院有限責(zé)任公司，北京100028）

0 引言

煤層屬于天然裂縫性儲(chǔ)層，力學(xué)強(qiáng)度低，具有顯著的應(yīng)力敏感性。流固耦合作用，又稱地質(zhì)力學(xué)效應(yīng)。研究表明，在數(shù)值模擬過程中是否考慮流固耦合作用，對(duì)于準(zhǔn)確預(yù)測(cè)煤層氣產(chǎn)能至關(guān)重要［1-2］。相比其他裂縫性儲(chǔ)層，煤層的流固耦合作用更為復(fù)雜，不僅存在有效應(yīng)力效應(yīng)，還包括基質(zhì)膨脹或收縮作用?；|(zhì)因煤層氣吸附或解吸現(xiàn)象產(chǎn)生形變，吸附使基質(zhì)膨脹，導(dǎo)致有效滲流孔道和滲透率均減??；解吸使基質(zhì)收縮，導(dǎo)致有效滲流孔道和滲透率均增大。

為刻畫煤層流固耦合作用，研究人員［3-6］曾提出包括ARI模型、Palmer&Mansoori模型和Shi&Durucan模型等在內(nèi)的多種模型。ARI模型［3］為經(jīng)驗(yàn)公式，沒有地質(zhì)力學(xué)理論基礎(chǔ)，基質(zhì)膨脹或收縮作用引起的應(yīng)變量與煤層氣吸附量成正比；Palmer模型［4］基于地質(zhì)力學(xué)理論，認(rèn)為煤層是均質(zhì)各向同性線彈性孔隙介質(zhì)，并將基質(zhì)膨脹或收縮應(yīng)變等效類比為熱膨脹應(yīng)變；Shi模型［5］不同于 Palmer模型，滲透率與水平有效應(yīng)力呈對(duì)數(shù)關(guān)系。這些模型均屬于解析或經(jīng)驗(yàn)流固耦合模型，具有形式簡(jiǎn)潔直觀、便于與商業(yè)煤層氣模擬軟件結(jié)合的優(yōu)勢(shì)，但缺點(diǎn)是須引入較多強(qiáng)假設(shè)，如固定上覆應(yīng)力與單軸向應(yīng)變假設(shè)，導(dǎo)致滲透率計(jì)算結(jié)果失真，影響產(chǎn)能預(yù)測(cè)精度。由此通過引入煤巖形變本構(gòu)方程來準(zhǔn)確刻畫煤層地質(zhì)力學(xué)效應(yīng)，即建立全流固耦合數(shù)學(xué)模型，以期獲得更準(zhǔn)確的儲(chǔ)層物性參數(shù)及產(chǎn)量預(yù)測(cè)結(jié)果。

1 全流固耦合數(shù)學(xué)模型

煤層氣藏全流固耦合數(shù)學(xué)模型，包括多過程流動(dòng)方程和多孔介質(zhì)地質(zhì)力學(xué)方程。前者采用三孔雙滲模型來描述煤層孔隙結(jié)構(gòu)，可進(jìn)一步發(fā)展目前的雙孔單滲模型；后者通過計(jì)算煤層應(yīng)力、應(yīng)變等，可獲得更準(zhǔn)確的孔隙度和滲透率。

1.1 多過程流動(dòng)方程

煤層屬于非常規(guī)天然氣儲(chǔ)層，既是生氣巖又是儲(chǔ)氣巖，原始地層條件下煤層氣的主要賦存狀態(tài)為吸附態(tài)而非游離態(tài)。目前常用雙孔單滲模型來描述煤層［6］，即煤層由連續(xù)分布的裂縫系統(tǒng)以及離散分布的基質(zhì)系統(tǒng)構(gòu)成，相應(yīng)的流動(dòng)過程為地層孔隙壓力隨排水下降，煤層氣由基質(zhì)解吸并擴(kuò)散進(jìn)入裂縫，而后以滲流方式進(jìn)入井筒。在礦場(chǎng)應(yīng)用中，這一模型的計(jì)算結(jié)果與實(shí)際生產(chǎn)數(shù)據(jù)存在較大偏差（煤層氣產(chǎn)量被高估、產(chǎn)水量被低估、氣體突破時(shí)間比實(shí)際早）［7-9］，有必要對(duì)煤層多過程物質(zhì)運(yùn)移和煤層孔隙結(jié)構(gòu)進(jìn)行細(xì)分再認(rèn)識(shí)，以得到合理準(zhǔn)確的產(chǎn)量計(jì)算結(jié)果，為此進(jìn)一步發(fā)展出三孔雙滲模型。3套孔隙系統(tǒng)是指基質(zhì)、中孔隙系統(tǒng)和裂縫系統(tǒng)等?；|(zhì)離散分布，是煤層氣主要的吸附和存儲(chǔ)空間；中孔隙系統(tǒng)與裂縫系統(tǒng)雖分布連續(xù)，但物性差異顯著，中孔隙系統(tǒng)孔隙度大而滲透率小，裂縫系統(tǒng)孔隙度小而滲透率大，中孔隙系統(tǒng)提供地層水存儲(chǔ)空間，而裂縫系統(tǒng)提供主要的滲流通道。三孔雙滲模型的流動(dòng)過程為：首先進(jìn)行排水降壓，使煤層氣解吸，基質(zhì)內(nèi)煤層氣在濃度梯度驅(qū)動(dòng)下擴(kuò)散進(jìn)入中孔隙系統(tǒng)；然后以竄流方式進(jìn)入裂縫；最后以滲流方式進(jìn)入井筒并采出。

1.1 .1 質(zhì)量守恒方程

忽略煤層氣在水中的溶解和水分揮發(fā)，建立質(zhì)量守恒方程，由連續(xù)性方程和達(dá)西滲流方程構(gòu)成。

連續(xù)性方程為

達(dá)西滲流方程為

式中：ρ為流體密度，kg/m3；下標(biāo)l為流體相、η為孔隙系統(tǒng)、s為基質(zhì)；Φ為孔隙度，%；S為流體飽和度，%；V為運(yùn)動(dòng)速度，m/s；q為源匯項(xiàng)，kg/（m3·s）；t為時(shí)間，s；K為滲透率mD，Kr為相對(duì)滲透率；μ為黏度，mPa·s；P為流體壓力，Pa為梯度，·為散度。

將式（2）代入式（1），可得

將式（3）右邊第二項(xiàng)展開，可得

同時(shí)根據(jù)定義，固相運(yùn)動(dòng)速度Vs與固相位移u的關(guān)系滿足

煤層氣開采過程中固相位移與固相運(yùn)動(dòng)速度相對(duì)較小，整理式（4）與式（5）后得

式中：u為位移，m；e為煤巖體應(yīng)變。對(duì)于三孔雙滲模型，質(zhì)量守恒方程具體形式為

式中：下標(biāo)w為水相、g為氣相、f為裂縫、m為中孔隙系統(tǒng)；qgms為煤層氣解吸量，kg/（m3·s）；qgmf和qwmf分別為煤層氣和水在中孔隙系統(tǒng)與裂縫之間的竄流量，kg/（m3·s）；qgwell和qwwell分別為煤層氣和水的源匯項(xiàng)，一般指開發(fā)井的生產(chǎn)量或注入量，kg/（m3·s）。

1.1.2 吸附及擴(kuò)散方程

煤層氣吸附采用Langmuir等溫吸附曲線描述，基質(zhì)表面的平衡吸附量為

式中：C（Pgm）為煤層氣平衡吸附體積濃度，m3/m3；Pgm為中孔隙系統(tǒng)的煤層氣壓力，Pa；VL為基質(zhì)對(duì)煤層氣的Langmuir吸附量，m3/m3；PL為煤層氣Langmuir吸附壓力，Pa。

對(duì)于基質(zhì)中煤層氣的擴(kuò)散運(yùn)移，一般采用擬穩(wěn)態(tài)非平衡模型進(jìn)行描述，根據(jù)菲克第一定律，推導(dǎo)得出煤層氣擴(kuò)散方程

式中：C為基質(zhì)中煤層氣的平均濃度，m3/m3；τ是解吸附時(shí)間，s。

1.1.3 竄流方程

中孔隙系統(tǒng)與裂縫系統(tǒng)存在勢(shì)差，在此作用下產(chǎn)生竄流或物質(zhì)交換，中孔隙系統(tǒng)與裂縫系統(tǒng)之間的竄流量qgmf和qwmf的計(jì)算方程［9］分別為

式中：δ為形狀因子，m-2。

1.2 地質(zhì)力學(xué)方程

1.2.1 煤巖本構(gòu)方程

根據(jù)Biot線彈性理論，對(duì)于各向同性彈性的煤層，本構(gòu)方程為

式中：T 為應(yīng)力張量，Pa；Te為有效應(yīng)力張量，Pa；Γ為應(yīng)變張量；I為單位張量；G為剪切模量，Pa；λ為平均軸向模量，Pa；es是基質(zhì)體應(yīng)變。

應(yīng)變與位移的關(guān)系為

動(dòng)量距平衡方程為

聯(lián)立式（15）—（18），可得

對(duì)式（19）取散度，最終將轉(zhuǎn)化成如下簡(jiǎn)潔形式?

煤層形變受制于有效應(yīng)力與基質(zhì)膨脹或收縮效應(yīng)。實(shí)驗(yàn)表明，基質(zhì)體應(yīng)變量與氣體壓力的關(guān)系呈現(xiàn)Langmuir曲線特征［4］

式中：εL為基質(zhì)Langmuir體應(yīng)變量。

1.2.2 孔隙度與滲透率方程

煤層滲透率的影響因素包括有效應(yīng)力效應(yīng)、基質(zhì)膨脹/收縮作用和氣體滑脫效應(yīng)等。有效應(yīng)力升高導(dǎo)致滲透率降低［10-11］；煤層氣吸附或解吸導(dǎo)致基質(zhì)膨脹或收縮，同樣可使?jié)B透率發(fā)生變化［3］；此外，煤層滲透率低，氣體滑脫效應(yīng)同樣不可忽略［12］。

通過數(shù)學(xué)推導(dǎo)，裂縫與中孔隙系統(tǒng)的孔隙度為

式中：Vpη為孔隙系統(tǒng) η的孔隙體積，m3；Vb為煤巖總體積，m3。中孔隙系統(tǒng)與裂縫的孔隙體積變化方程［13］為

考慮到煤巖體應(yīng)變dVb/Vb=de和平均主應(yīng)力，孔隙度方程為

式中：cpn為孔隙體積壓縮系數(shù)，Pa-1；β為有效應(yīng)力系數(shù)；Kb為煤巖體積模量，Pa。

裂縫一般采用火柴棒束模型來描述，這樣滲透率與孔隙度呈三次方關(guān)系

式中：Pkf為裂縫Klinkenberg系數(shù)，Pa。上式右端第一項(xiàng)表征有效應(yīng)力和基質(zhì)收縮效應(yīng)，第二項(xiàng)表征氣體滑脫效應(yīng)。

對(duì)于中孔隙系統(tǒng)，火柴棒束模型已不再適用，采用Kozeny-Carman模型來表征滲透率與孔隙度的關(guān)系［14-16］，可表達(dá)為

式中：Pkm為中孔隙系統(tǒng)Klinkenberg系數(shù)，Pa。根據(jù)實(shí)驗(yàn)結(jié)果［17-20］，Pkf和 Pkm分別取值 0.28 MPa 和0.84 MPa。

2 數(shù)值求解

根據(jù)全流固耦合數(shù)學(xué)模型，控制方程包括質(zhì)量守恒方程［式（7）～（10）］、擴(kuò)散方程［式（12）］和地質(zhì)力學(xué)方程［式（20）］等，主要未知量包括水相壓力（Pwf和 Pwm）、水相飽和度（Swf和 Swm）、基質(zhì)煤層氣平均吸附濃度（C）和煤巖體應(yīng)變（e）等。采用全隱式有限體積法離散質(zhì)量守恒方程［21-24］，得

將式（27）展開，可得

對(duì)于擴(kuò)散方程式（12）和地質(zhì)力學(xué)方程式（20），有

采用 Newton-Raphson 迭代全隱式方程［25-28］求解上述離散方程，可獲得三孔雙滲全流固耦合煤層氣模擬器。主要未知量 Pwf，Pwm，Swf，Swm，C，e 被同時(shí)解出，據(jù)此同步更新煤層孔隙度/滲透率、流體PVT屬性、流度和傳導(dǎo)系數(shù)等，進(jìn)而獲得氣、油產(chǎn)量，然后進(jìn)入下一個(gè)時(shí)間步計(jì)算。

3 模型及算法驗(yàn)證

通過與商業(yè)模擬軟件Coalgas進(jìn)行對(duì)比分析，來驗(yàn)證三孔雙滲全流固耦合數(shù)學(xué)模型及模擬器的準(zhǔn)確性。因二者所采用的流固耦合模型不同，本次研究為全流固耦合數(shù)學(xué)模型，而Colagas使用的是ARI模型，所以為使計(jì)算結(jié)果具有可比性，可忽略煤層地質(zhì)力學(xué)作用，即使孔隙度和滲透率保持不變。

煤層概念模型的輸入?yún)?shù)為：煤層長(zhǎng)度為360 m、寬度為360 m、厚度為5 m；網(wǎng)格平面劃分為51×51，幾何結(jié)構(gòu)對(duì)稱，生產(chǎn)井位于中心并保持定壓生產(chǎn)；初始地層條件下裂縫系統(tǒng)孔隙度為0.5%，滲透率為5.0 mD，中孔隙系統(tǒng)孔隙度為5.0%，滲透率為0.05 mD；煤層氣吸附參數(shù)、流體PVT屬性參數(shù)、煤巖力學(xué)參數(shù)、裂縫系統(tǒng)與中孔隙系統(tǒng)物性參數(shù)如表1所列；煤層初始狀態(tài)是非飽和吸附，初始孔隙壓力為6.895 MPa，臨界解吸壓力為5.17 MPa。

煤層氣與水的產(chǎn)量模擬結(jié)果如圖1所示，模擬器計(jì)算結(jié)果與Coalgas的模擬結(jié)果一致，初步驗(yàn)證了全流固耦合數(shù)學(xué)模型的準(zhǔn)確性和可靠性。

4 流固耦合效應(yīng)的影響

流固耦合效應(yīng)包括有效應(yīng)力效應(yīng)和基質(zhì)收縮作用，下面分別剖析這2種作用對(duì)滲透率和煤層氣生產(chǎn)動(dòng)態(tài)的影響。

4.1 有效應(yīng)力效應(yīng)

為了評(píng)估有效應(yīng)力效應(yīng)的影響，可暫且忽略基質(zhì)收縮作用，這樣煤層在力學(xué)上等效于普通裂縫性儲(chǔ)層。有效應(yīng)力效應(yīng)作用程度依賴于煤巖力學(xué)強(qiáng)度，而楊氏模量是衡量煤巖力學(xué)強(qiáng)度的關(guān)鍵參數(shù)。本次研究選擇了3種不同力學(xué)強(qiáng)度的煤層，對(duì)應(yīng)楊氏模量E值分別為 1.25 GPa，2.50 GPa和 5.00 GPa，同時(shí)選取剛性儲(chǔ)層作為對(duì)照。不同力學(xué)強(qiáng)度儲(chǔ)層的煤層氣產(chǎn)量如圖2所示?？芍獜椥詢?chǔ)層的煤層氣產(chǎn)量與剛性儲(chǔ)層存在明顯差異，彈性儲(chǔ)層的煤層氣高峰日產(chǎn)量較低，且高峰出現(xiàn)時(shí)機(jī)較晚，這是有效應(yīng)力效應(yīng)作用的結(jié)果。楊氏模量越大，煤層氣高峰產(chǎn)量越大，且出現(xiàn)時(shí)機(jī)越早。

表1 煤層氣模擬輸入?yún)?shù)Table 1 Simulated input data of coalbed methane

圖1 全流固耦合模擬算法與Coalgas模擬結(jié)果對(duì)比Fig.1 Comparison of simulation results between the fully coupled fluid flow and geomechanics model and Coalgas

圖2 不同楊氏模量煤層的產(chǎn)氣量對(duì)比Fig.2 Comparison of gas rate for reservoirs with different Young's modulus

隨著煤層氣的開采，地層孔隙壓力減小，有效應(yīng)力（壓應(yīng)力）增大，導(dǎo)致孔隙收縮，孔隙度持續(xù)下降，而滲透率因氣體滑脫作用的影響而呈現(xiàn)不同的變化趨勢(shì)。以裂縫滲透率為例，不同楊氏模量下的滲透率隨孔隙壓力的變化情況如圖3所示，圖中對(duì)滲透率進(jìn)行了歸一化處理。對(duì)于剛性煤層，當(dāng)孔隙壓力由初始值（6.895 MPa）降低至臨界解吸壓力（5.170 MPa）之前，孔隙中無游離氣，滲透率保持不變；之后，煤層氣由基質(zhì)中解吸，滲透率受氣體滑脫作用影響而不斷升高。對(duì)于彈性煤層，滲透率的變化受有效應(yīng)力和氣體滑脫作用雙重作用影響，且二者作用方向相反，前者導(dǎo)致滲透率減小、后者導(dǎo)致滲透率增大：有效應(yīng)力效應(yīng)在開采初期占主導(dǎo)，滲透率減小，但后期氣體滑脫作用增大并逐漸占優(yōu)，滲透率可止降為升。

圖3 不同楊氏模量煤層井網(wǎng)格的裂縫滲透率變化情況對(duì)比Fig.3 Comparison of normalized fracture permeability of well block for reservoirs with different Young's modulus

圖4 不同Langmuir體應(yīng)變量煤層的產(chǎn)氣量對(duì)比Fig.4 Comparison of gas rate for reservoirs with different Langmuir volumetric strains

4.2 基質(zhì)收縮作用

為評(píng)估基質(zhì)收縮作用的影響，保持煤巖楊氏模量不變，通過改變基質(zhì)形變強(qiáng)度，相應(yīng)Langmuir體應(yīng)變量 εL分別取 0，0.006，0.012與 0.024等，煤巖楊氏模量為2.50 GPa。εL越大，單位煤層氣吸附量引起的基質(zhì)收縮量越大，對(duì)滲透率與孔隙度的影響也越大，εL=0即不考慮基質(zhì)應(yīng)變。

不同基質(zhì)形變強(qiáng)度下的煤層氣產(chǎn)量如圖4所示：Langmuir體應(yīng)變量εL對(duì)煤層氣產(chǎn)量的動(dòng)態(tài)影響顯著，εL越大，峰值產(chǎn)量越大，且出現(xiàn)時(shí)間越早。可見對(duì)煤層氣產(chǎn)量的影響而言，基質(zhì)收縮作用與有效應(yīng)力效應(yīng)作用方向相反。

不同基質(zhì)形變強(qiáng)度下裂縫系統(tǒng)與中孔隙系統(tǒng)的滲透率變化情況如圖5所示，當(dāng)εL≠0時(shí)，滲透率變化受有效應(yīng)力效應(yīng)、基質(zhì)收縮作用和氣體滑脫效應(yīng)等3種作用影響，其中有效應(yīng)力效應(yīng)會(huì)使?jié)B透率減小，但基質(zhì)收縮作用和氣體滑脫效應(yīng)卻使?jié)B透率增大。當(dāng)孔隙壓力高于臨界解吸壓力時(shí)，煤層氣無解吸，僅有效應(yīng)力效應(yīng)發(fā)揮作用，故滲透率不變；隨著煤層氣解吸，誘發(fā)基質(zhì)收縮作用和氣體滑脫效應(yīng)，在三者共同作用下，裂縫滲透率先減小后增大，εL=0.024時(shí)，最終滲透率甚至可達(dá)初始值的2.2倍。中孔隙系統(tǒng)滲透率僅在初期微弱下降，而后快速上升至初始值的2.3倍，這主要是氣體滑脫作用的影響。裂縫滲透率受基質(zhì)Langmuir體應(yīng)變量εL的影響顯著，但對(duì)中孔隙系統(tǒng)的影響幾乎忽略不計(jì)，這是因?yàn)橹锌紫断到y(tǒng)的力學(xué)強(qiáng)度大于裂縫系統(tǒng)（cpf＞cpm）。

圖5 不同Langmuir體應(yīng)變量煤層井網(wǎng)格的滲透率變化情況對(duì)比Fig.5 Comparison of normalized permeability of well block for reservoirs with different Langmuir volumetric strains

5 結(jié)論

（1）構(gòu)建了煤層氣藏全流固耦合數(shù)學(xué)模型，包括多孔介質(zhì)地質(zhì)力學(xué)方程和三孔雙滲流動(dòng)方程，前者定量刻畫了煤巖形變，獲得了更準(zhǔn)確的孔滲參數(shù)；后者進(jìn)一步完善了目前的雙孔單滲模型，實(shí)現(xiàn)了煤層氣多過程物質(zhì)運(yùn)移的準(zhǔn)確表征。

（2）開發(fā)了基于全隱式有限體積法的全流固耦合、三孔雙滲煤層氣數(shù)值模擬器，通過與商業(yè)模擬軟件Coalgas進(jìn)行對(duì)比，驗(yàn)證了其準(zhǔn)確性與有效性。

（3）地質(zhì)力學(xué)效應(yīng)對(duì)煤層氣生產(chǎn)動(dòng)態(tài)預(yù)測(cè)至關(guān)重要。有效應(yīng)力效應(yīng)與基質(zhì)收縮作用均可顯著影響孔滲參數(shù)及煤層氣產(chǎn)量，前者降低滲透率，后者增大滲透率，二者作用方向相反，但通常不可抵消。有效應(yīng)力效應(yīng)的作用強(qiáng)度在開發(fā)初期大于基質(zhì)收縮作用，后期逆轉(zhuǎn)，導(dǎo)致裂縫滲透率呈現(xiàn)獨(dú)特的“U”型變化，最終值甚至可達(dá)初始值的數(shù)倍。