適用于深層頁(yè)巖氣井生產(chǎn)模擬和預(yù)測(cè)的自擴(kuò)散流熱耦合模型

夏 陽(yáng) 韋世明 金 衍 陳康平

1. 油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室·中國(guó)石油大學(xué)（北京） 2. 美國(guó)亞利桑那州立大學(xué)

0 引言

四川盆地及其周緣埋深介于3 500 ～4 500 m 的深層頁(yè)巖氣資源豐富，占該區(qū)域頁(yè)巖氣總資源量的85%以上[1]。深層頁(yè)巖氣儲(chǔ)層具有高溫、高壓、高地應(yīng)力的特點(diǎn)[2]，準(zhǔn)確預(yù)測(cè)深層頁(yè)巖氣井產(chǎn)能可以為氣藏開(kāi)發(fā)方案的制訂提供重要的依據(jù)。

頁(yè)巖氣井生產(chǎn)數(shù)據(jù)統(tǒng)計(jì)結(jié)果表明，采用達(dá)西滲流模型計(jì)算得到的氣井產(chǎn)量往往低于實(shí)際氣產(chǎn)量[3]。現(xiàn)有的頁(yè)巖氣滲流模型則在達(dá)西滲流模型的基礎(chǔ)上，增加考慮滑脫、擴(kuò)散等因素的影響[4-7]，對(duì)達(dá)西滲流模型進(jìn)行修正。Klinkenberg[8]最早通過(guò)巖心滲透率測(cè)試實(shí)驗(yàn)，發(fā)現(xiàn)同一塊巖心的氣測(cè)滲透率大于其液測(cè)滲透率，并且氣體壓力越低，氣測(cè)滲透率越高；氣測(cè)滲透率與氣體壓力倒數(shù)近似成正比，此即為氣體滑脫效應(yīng)，又稱為Klinkenberg 效應(yīng)。2007 年，Javadpour等[9]將努森擴(kuò)散引入頁(yè)巖氣滲流模型中，并給出頁(yè)巖孔隙中努森擴(kuò)散系數(shù)的具體形式。2015 年，吳克柳等[10]采用努森數(shù)定義了頁(yè)巖氣流動(dòng)過(guò)程中努森擴(kuò)散和黏性流的權(quán)重系數(shù)，發(fā)現(xiàn)努森擴(kuò)散在低壓條件下對(duì)頁(yè)巖氣滲流的貢獻(xiàn)不可忽略。由于滑脫效應(yīng)僅發(fā)生在稀薄氣體中，而實(shí)測(cè)氣體滑脫因子介于0.1 ～1.0 MPa[11-13]，對(duì)于深層頁(yè)巖氣儲(chǔ)層而言，由于地層壓力較高，滑脫效應(yīng)對(duì)氣體流動(dòng)影響很小。Patzek[14]、Chen 等[15]指出，在深部致密氣儲(chǔ)層中氣體滑脫效應(yīng)和努森擴(kuò)散可以忽略。

頁(yè)巖儲(chǔ)層基質(zhì)滲透率極低，應(yīng)力敏感效應(yīng)使得生產(chǎn)過(guò)程中滲透率將進(jìn)一步降低[16-17]。由于深層頁(yè)巖儲(chǔ)層壓力、溫度高，在氣井生產(chǎn)過(guò)程中氣體流動(dòng)極大地依賴其自身膨脹，國(guó)內(nèi)外大多數(shù)學(xué)者將自擴(kuò)散效應(yīng)作為頁(yè)巖孔隙中氣體的主要流動(dòng)方式，并建立了初步的氣體自擴(kuò)散流動(dòng)模型[18-21]。Larry[22]研究了溶解氣驅(qū)油流動(dòng)機(jī)理，發(fā)現(xiàn)氣體膨脹是原油流動(dòng)的主要驅(qū)動(dòng)力。Kolesar 等[23]，Carlson 等[24]、Amann-Hildenbrand 等[25]隨后提出擴(kuò)散效應(yīng)是致密氣流動(dòng)的主要方式。Wan 等[20]基于氣體擴(kuò)散系數(shù)計(jì)算方法，模擬頁(yè)巖油吞吐生產(chǎn)過(guò)程中的油氣產(chǎn)量，對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了很好地?cái)M合。Cronin 等[21]基于Fick 定律，建立了不同氣體在頁(yè)巖油吞吐生產(chǎn)過(guò)程中的質(zhì)量擴(kuò)散方程，計(jì)算的油氣產(chǎn)量與現(xiàn)場(chǎng)生產(chǎn)數(shù)據(jù)吻合程度高，但模型中擴(kuò)散系數(shù)是基于對(duì)生產(chǎn)數(shù)據(jù)進(jìn)行擬合，尚未建立普適的自擴(kuò)散流動(dòng)模型。Hoteit[26]提出基于濃度梯度的Fick 定律不能描述致密儲(chǔ)層中的氣體擴(kuò)散行為。Jin等[27]從可壓縮流體的N-S 方程出發(fā)，在Klainerman等[28]提出的小馬赫數(shù)流動(dòng)理論基礎(chǔ)上，建立了致密氣一次開(kāi)采氣體自擴(kuò)散流動(dòng)控制方程，并基于現(xiàn)場(chǎng)實(shí)例對(duì)自擴(kuò)散流動(dòng)模型的特征和適用性進(jìn)行了研究。

深層頁(yè)巖儲(chǔ)層具有高溫、高壓的特征，CH4的體積黏度遠(yuǎn)大于其剪切黏度且受溫度的影響較大[29]。在頁(yè)巖氣井生產(chǎn)過(guò)程中，由于壓裂液冷卻、氣體膨脹使得近井溫度降低，溫度的變化將對(duì)氣體流動(dòng)產(chǎn)生影響，進(jìn)而影響深層頁(yè)巖氣井的產(chǎn)能[30]。截至目前，國(guó)內(nèi)外關(guān)于流體自擴(kuò)散流動(dòng)模型的研究都未考慮高溫、高壓的影響，尤其是溫度變化對(duì)氣體物理化學(xué)性質(zhì)及傳質(zhì)速度的影響規(guī)律尚未見(jiàn)到闡明。

為此，筆者基于Jin 等[27]提出的自擴(kuò)散流動(dòng)模型，建立了考慮滲流場(chǎng)和溫度場(chǎng)耦合的自擴(kuò)散流動(dòng)模型（以下簡(jiǎn)稱自擴(kuò)散流熱耦合模型），然后，將該模型與以達(dá)西定律和努森擴(kuò)散為基礎(chǔ)的滲流模型（以下簡(jiǎn)稱修正達(dá)西模型）計(jì)算結(jié)果進(jìn)行對(duì)比；基于自擴(kuò)散流熱耦合模型，分析考慮溫度場(chǎng)變化的深層頁(yè)巖氣自擴(kuò)散流動(dòng)行為，進(jìn)而探討了井底溫度變化對(duì)深層頁(yè)巖氣采出程度的影響；在此基礎(chǔ)上，將自擴(kuò)散流熱耦合模型應(yīng)用于四川盆地長(zhǎng)寧區(qū)塊上奧陶統(tǒng)五峰組—下志留統(tǒng)龍馬溪組一口頁(yè)巖氣水平井的生產(chǎn)模擬，以期為頁(yè)巖氣井產(chǎn)能的準(zhǔn)確預(yù)測(cè)提供支撐。

1 自擴(kuò)散流熱耦合模型的建立

在地層高溫高壓條件下，氣體狀態(tài)方程為：

式中ρG表示氣體密度，kg/m3；p 表示氣體壓力，Pa；M 表示氣體分子摩爾質(zhì)量，kg/mol；Z 表示氣體偏差因子；R 表示普適氣體常數(shù)，取值為8.314 J/(mol·K)；T 表示氣體溫度，K。

氣體偏差因子和剪切黏度隨氣體組成、壓力和溫度發(fā)生變化，筆者采用Heidaryan 等[31]提出的經(jīng)驗(yàn)公式，即式中ppr表示氣體對(duì)比壓力，取值為1；Tpr表示氣體對(duì)比溫度，取值為1；A1—A11表示由實(shí)驗(yàn)數(shù)據(jù)擬合得到的常系數(shù)，依次取值為3.252 838、-1.306 424× 10-1、-6.444 919 4×10-1、-1.518 028、-5.391 019、 -1.379 588×10-2、6.600 633×10-2、6.120 783×10-1、2.317 431、1.632 223×10-1、5.660 595×10-1；μ表示氣體剪切黏度，mPa·s；B1—B10表示由實(shí)驗(yàn)數(shù)據(jù)擬合得到的常系數(shù)，依次取值為1.022 872、 -1.651 432、5.757 386、-7.389 282×10-2、8.389 065× 10-2、2.977 476 10-1、-1.451 318、4.682 506、1.918 239、 -9.944 968×10-2。

當(dāng)考慮溫度場(chǎng)的變化時(shí)，氣體體積流量與熱通量可分別表示為[32]：

式中qG表示通過(guò)巖石單位面積的氣體體積流量，m/s；Dj表示自擴(kuò)散系數(shù)，m2/s；βT表示熱滲透系數(shù)，m2/(s·K)；qT表示通過(guò)巖石單位面積的熱通量，W/m2；CG表示氣體比熱容，kJ/(kg·K)；φ 表示基質(zhì)孔隙度；kG、ks分別表示氣體、巖石的熱傳導(dǎo)系數(shù)，W/(K·m)。

Dj計(jì)算式為[27]：

式中μb表示氣體體積黏度，mPa·s。

圖1 甲烷與T 關(guān)系曲線圖[29]

多孔介質(zhì)中流體的質(zhì)量守恒方程為：

含流體的多孔介質(zhì)熱平衡方程為[30-31]：

式中ρs表示巖石密度，kg/m3；Cs表示巖石比熱容，J/(kg·K)；αG表示氣體熱膨脹系數(shù)，K-1；KG表示氣體體積模量，Pa。

將式（1）、（4）代入式（7），則氣體自擴(kuò)散流動(dòng)控制方程為：

將式（4）、（5）代入式（8）中，得到溫度場(chǎng)控制方程[30-31]，即

其中

將式（9）、（10）結(jié)合，即為自擴(kuò)散流熱耦合模型，其中待求解變量為p 與T。該模型僅將游離氣的膨脹考慮為氣體流動(dòng)的驅(qū)動(dòng)力，暫未考慮吸附氣的影響。

2 模型分析

筆者針對(duì)壓裂井模擬區(qū)域建立了一維、二維模型（圖2），模型參數(shù)如表1 所示。一維模型的左端代表裂縫面，其截面積為1 m2，邊界條件為左端定壓定溫、右端封閉。筆者采用有限元軟件COMSOL中的數(shù)學(xué)模塊對(duì)式（9）、（10）的偏微分方程進(jìn)行求解。通過(guò)對(duì)比自擴(kuò)散流熱耦合模型（考慮溫度場(chǎng)的影響）與修正達(dá)西滲流模型的計(jì)算結(jié)果，明確自擴(kuò)散流熱耦合模型的正確性；在此基礎(chǔ)上，分析溫度場(chǎng)對(duì)氣體自擴(kuò)散流動(dòng)行為的影響，進(jìn)而研究溫度場(chǎng)變化對(duì)深層頁(yè)巖氣井產(chǎn)量的影響。

圖2 壓裂井模擬區(qū)域劃分示意圖

表1 模型參數(shù)表

2.1 計(jì)算結(jié)果對(duì)比

常規(guī)滲流模型認(rèn)為氣體流動(dòng)受到壓力梯度和濃度梯度的影響?；谶_(dá)西定律，增加考慮努森擴(kuò)散的影響，流量計(jì)算式為[6-7]：

式中J 表示通過(guò)單位巖石截面積的氣體質(zhì)量流量，kg/(m2·s)；K 表示巖石滲透率，mD；D 表示努森擴(kuò)散系數(shù)，m2/s；M 表示氣體摩爾質(zhì)量，g/mol；C 表示氣體濃度，mol/m3。

僅使用式（11）代替式（4）右側(cè)自擴(kuò)散流量時(shí)，就可以得到修正達(dá)西模型，之后可以采用有限元數(shù)值模擬軟件COMSOL 中的數(shù)學(xué)模塊進(jìn)行求解。本節(jié)基于一維模型（圖2）進(jìn)行相關(guān)計(jì)算。

圖3 自擴(kuò)散流熱耦合模型與修正達(dá)西模型計(jì)算孔隙壓力剖面對(duì)比圖

圖4 自擴(kuò)散流熱耦合模型與修正達(dá)西模型計(jì)算氣體密度剖面對(duì)比圖

通過(guò)對(duì)比自擴(kuò)散流熱耦合模型與修正達(dá)西模型計(jì)算的壓力剖面，發(fā)現(xiàn)在同一時(shí)間，自擴(kuò)散流熱耦合模型計(jì)算的壓力擴(kuò)散速度更快（圖3）。如圖3、4 所示，采用修正達(dá)西模型計(jì)算得到壓力剖面和氣體密度剖面，壓力、氣體密度均在裂縫面處最小。由于氣體密度、濃度與壓力正相關(guān)，因而圖3 中基于自擴(kuò)散流熱耦合模型計(jì)算的壓力剖面與修正達(dá)西模型計(jì)算得到的壓力剖面變化趨勢(shì)一致。然而，由于溫度場(chǎng)的影響，自擴(kuò)散流熱耦合模型計(jì)算得到的氣體密度剖面在基質(zhì)內(nèi)部出現(xiàn)最小值，且該最小值對(duì)應(yīng)的位置隨著生產(chǎn)時(shí)間延長(zhǎng)向基質(zhì)內(nèi)部移動(dòng)（圖4），可以看出，氣體自擴(kuò)散并非僅受氣體濃度梯度的影響[33]。如圖5 所示，采用自擴(kuò)散流熱耦合模型計(jì)算的氣體累計(jì)產(chǎn)出質(zhì)量遠(yuǎn)大于修正達(dá)西模型計(jì)算值，說(shuō)明基于氣體膨脹驅(qū)動(dòng)計(jì)算的產(chǎn)氣量大于基于氣體壓力和濃度梯度驅(qū)動(dòng)計(jì)算的產(chǎn)氣量，從而解釋了頁(yè)巖氣井生產(chǎn)初期氣產(chǎn)量通常大于修正達(dá)西模型計(jì)算值的原因[3]。

圖5 自擴(kuò)散流熱耦合模型與修正達(dá)西模型模擬 氣體累計(jì)產(chǎn)出質(zhì)量對(duì)比圖

2.2 溫度場(chǎng)變化對(duì)氣體自擴(kuò)散流動(dòng)行為的影響

下面分析溫度場(chǎng)變化對(duì)深層頁(yè)巖氣自擴(kuò)散流動(dòng)行為和產(chǎn)量的影響，本節(jié)計(jì)算采用的模型參數(shù)、邊界條件與2.1 部分相同，對(duì)于恒溫地層，則考慮地層溫度保持不變。

如圖6 所示，若不考慮溫度場(chǎng)變化的影響，自擴(kuò)散系數(shù)沿裂縫面向基質(zhì)內(nèi)部逐漸減??；考慮溫度場(chǎng)變化后，自擴(kuò)散系數(shù)沿裂縫面向基質(zhì)內(nèi)部呈先增大后減小的趨勢(shì)，且隨著生產(chǎn)進(jìn)行，自擴(kuò)散系數(shù)峰值所在位置逐漸向基質(zhì)內(nèi)部移動(dòng)。使用氣體狀態(tài)方程，將自擴(kuò)散系數(shù)計(jì)算式變形為：

式中bg表示氣體體積黏度與剪切黏度的比值。

由式（12）可以看出，Dj隨著T 升高而增大，隨著p 升高而減小。在裂縫面溫度低于儲(chǔ)層初始溫度時(shí)，從裂縫面向基質(zhì)內(nèi)部，地層溫度和孔隙壓力均不斷升高，自擴(kuò)散系數(shù)呈先增大后減小的趨勢(shì)。

不考慮溫度場(chǎng)變化時(shí)，自擴(kuò)散流熱耦合模型計(jì)算得到的氣體密度剖面（圖7）與修正達(dá)西模型計(jì)算得到的氣體密度剖面形狀相同（圖4）。如圖7 所示，隨著地層溫度降低，從裂縫面向基質(zhì)內(nèi)部，氣體密度先降低后增大，其原因是井底溫度的降低使得自擴(kuò)散系數(shù)降低，造成氣體流動(dòng)速度在裂縫面附近減小，引起氣體在近裂縫面區(qū)域發(fā)生質(zhì)量堆積，即氣體自擴(kuò)散的近井阻塞；隨著生產(chǎn)時(shí)間進(jìn)行，低溫區(qū)域范圍擴(kuò)大，發(fā)生阻塞的區(qū)域范圍增大，氣體密度谷值所在位置向基質(zhì)內(nèi)部移動(dòng)。

如圖8 所示，考慮井底溫度低于地層溫度時(shí)，預(yù)測(cè)的氣體累計(jì)產(chǎn)出質(zhì)量低于恒溫地層。由此可知，近井低溫將使頁(yè)巖氣井產(chǎn)量降低。

3 井底溫度變化對(duì)深層頁(yè)巖氣采出程度的影響探討

基于第2 節(jié)的模型參數(shù)，在裂縫面設(shè)定不同的溫度，分別低于初始地層溫度10 ～80 ℃，計(jì)算生產(chǎn)3 年末的氣采出程度（累計(jì)產(chǎn)氣量與模擬區(qū)域內(nèi)頁(yè)巖氣儲(chǔ)量的比值）。如圖9 所示，井底溫度越低，即溫差越大，氣采出程度越低；隨溫差增大，采出程度近直線下降，溫差從10 ℃增至80 ℃，采出程度降低2.3%。在壓裂過(guò)程中，壓裂液會(huì)冷卻裂縫面，并且隨壓裂液注入量增大和燜井時(shí)間延長(zhǎng)，裂縫面溫度還會(huì)進(jìn)一步降低。由于深層頁(yè)巖氣儲(chǔ)層具有高溫、高應(yīng)力特征，裂縫寬度更窄，在頁(yè)巖氣開(kāi)采過(guò)程中，上述兩方面的影響會(huì)更明顯，因而對(duì)于深層頁(yè)巖氣的開(kāi)采，應(yīng)該考慮溫度場(chǎng)變化所產(chǎn)生的影響。

圖9 溫度差對(duì)氣井生產(chǎn)3 年末氣采出程度影響曲線圖

4 實(shí)例井應(yīng)用

選取四川盆地長(zhǎng)寧區(qū)塊五峰組—龍馬溪組一口頁(yè)巖氣水平井進(jìn)行生產(chǎn)模擬，該井基礎(chǔ)參數(shù)如表2 所示，巖石和氣體熱力學(xué)參數(shù)如表1 所示。建立壓裂井二維模型（圖2），假設(shè)水力裂縫為無(wú)窮導(dǎo)流能力裂縫，則水力裂縫處邊界條件為定壓定溫，其中水力裂縫處的壓力則為由井口油壓計(jì)算得到的井底壓力[34]。

表2 實(shí)例井基礎(chǔ)參數(shù)表

如圖10 所示，考慮溫度場(chǎng)變化影響的模擬日產(chǎn)氣量與實(shí)際生產(chǎn)數(shù)據(jù)更加吻合。剛開(kāi)始生產(chǎn)時(shí)，若不考慮溫度場(chǎng)的變化，近裂縫處自擴(kuò)散系數(shù)計(jì)算值將大于實(shí)際值，從而使模擬日產(chǎn)氣量遠(yuǎn)高于實(shí)際結(jié)果。2014 年6 月18 日—7 月1 日，關(guān)井使得該井地層壓力得以恢復(fù)；由于井筒儲(chǔ)存效應(yīng)的影響，在開(kāi)井復(fù)產(chǎn)初期，實(shí)際日產(chǎn)氣量高于模擬結(jié)果。2014 年8 月12日以后，關(guān)井壓力恢復(fù)對(duì)產(chǎn)量模擬產(chǎn)生的影響消失，且考慮溫度場(chǎng)變化模擬的日產(chǎn)氣量與實(shí)際日產(chǎn)氣量的吻合度更好。在壓裂過(guò)程中，壓裂液將冷卻裂縫面，并且裂縫面溫度與壓裂液注入量、壓裂液浸泡時(shí)間均相關(guān)，而且裂縫面溫度與初始地層溫度相差越大，日氣產(chǎn)量越低。因此，對(duì)于采用大液量壓裂的深層頁(yè)巖氣井來(lái)說(shuō)，必須考慮溫度場(chǎng)變化的影響才能準(zhǔn)確預(yù)測(cè)單井氣產(chǎn)量。

5 結(jié)論

1）在相同參數(shù)下，采用自擴(kuò)散流熱耦合模型計(jì)算的頁(yè)巖氣產(chǎn)量大于修正達(dá)西模型計(jì)算值。

2）考慮溫度場(chǎng)變化時(shí)，自擴(kuò)散系數(shù)剖面出現(xiàn)峰值，氣體密度剖面出現(xiàn)谷值，并且峰值/谷值對(duì)應(yīng)數(shù)據(jù)點(diǎn)同步向地層內(nèi)部移動(dòng)，說(shuō)明自擴(kuò)散系數(shù)影響著氣體傳質(zhì)速率，且隨生產(chǎn)進(jìn)行，近井低溫對(duì)氣體自擴(kuò)散的影響范圍不斷增大。

圖10 長(zhǎng)寧區(qū)塊五峰組—龍馬溪組頁(yè)巖氣水平井生產(chǎn)模擬結(jié)果圖

3）當(dāng)井底溫度低于地層溫度時(shí)，近井氣體自擴(kuò)散系數(shù)減小，發(fā)生氣體近井阻塞現(xiàn)象，從而使氣井產(chǎn)量降低。

4）根據(jù)現(xiàn)場(chǎng)實(shí)例井的生產(chǎn)模擬結(jié)果，考慮溫度場(chǎng)的變化，自擴(kuò)散流熱耦合模型能夠更準(zhǔn)確地模擬深層頁(yè)巖氣井產(chǎn)量。