適用于深層頁巖氣井生產模擬和預測的自擴散流熱耦合模型

夏 陽 韋世明 金 衍 陳康平

1. 油氣資源與探測國家重點實驗室·中國石油大學（北京） 2. 美國亞利桑那州立大學

0 引言

四川盆地及其周緣埋深介于3 500 ～4 500 m 的深層頁巖氣資源豐富，占該區域頁巖氣總資源量的85%以上[1]。深層頁巖氣儲層具有高溫、高壓、高地應力的特點[2]，準確預測深層頁巖氣井產能可以為氣藏開發方案的制訂提供重要的依據。

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

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

深層頁巖儲層具有高溫、高壓的特征，CH4的體積黏度遠大于其剪切黏度且受溫度的影響較大[29]。在頁巖氣井生產過程中，由于壓裂液冷卻、氣體膨脹使得近井溫度降低，溫度的變化將對氣體流動產生影響，進而影響深層頁巖氣井的產能[30]。截至目前，國內外關于流體自擴散流動模型的研究都未考慮高溫、高壓的影響，尤其是溫度變化對氣體物理化學性質及傳質速度的影響規律尚未見到闡明。

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

1 自擴散流熱耦合模型的建立

在地層高溫高壓條件下，氣體狀態方程為：

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

氣體偏差因子和剪切黏度隨氣體組成、壓力和溫度發生變化，筆者采用Heidaryan 等[31]提出的經驗公式，即式中ppr表示氣體對比壓力，取值為1；Tpr表示氣體對比溫度，取值為1；A1—A11表示由實驗數據擬合得到的常系數，依次取值為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表示由實驗數據擬合得到的常系數，依次取值為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。

當考慮溫度場的變化時，氣體體積流量與熱通量可分別表示為[32]：

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

Dj計算式為[27]：

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

圖1 甲烷與T 關系曲線圖[29]

多孔介質中流體的質量守恒方程為：

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

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

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

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

其中

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

2 模型分析

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

圖2 壓裂井模擬區域劃分示意圖

表1 模型參數表

2.1 計算結果對比

常規滲流模型認為氣體流動受到壓力梯度和濃度梯度的影響。基于達西定律，增加考慮努森擴散的影響，流量計算式為[6-7]：

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

僅使用式（11）代替式（4）右側自擴散流量時，就可以得到修正達西模型，之后可以采用有限元數值模擬軟件COMSOL 中的數學模塊進行求解。本節基于一維模型（圖2）進行相關計算。

圖3 自擴散流熱耦合模型與修正達西模型計算孔隙壓力剖面對比圖

圖4 自擴散流熱耦合模型與修正達西模型計算氣體密度剖面對比圖

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

圖5 自擴散流熱耦合模型與修正達西模型模擬 氣體累計產出質量對比圖

2.2 溫度場變化對氣體自擴散流動行為的影響

下面分析溫度場變化對深層頁巖氣自擴散流動行為和產量的影響，本節計算采用的模型參數、邊界條件與2.1 部分相同，對于恒溫地層，則考慮地層溫度保持不變。

如圖6 所示，若不考慮溫度場變化的影響，自擴散系數沿裂縫面向基質內部逐漸減小；考慮溫度場變化后，自擴散系數沿裂縫面向基質內部呈先增大后減小的趨勢，且隨著生產進行，自擴散系數峰值所在位置逐漸向基質內部移動。使用氣體狀態方程，將自擴散系數計算式變形為：

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

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

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

如圖8 所示，考慮井底溫度低于地層溫度時，預測的氣體累計產出質量低于恒溫地層。由此可知，近井低溫將使頁巖氣井產量降低。

3 井底溫度變化對深層頁巖氣采出程度的影響探討

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

圖9 溫度差對氣井生產3 年末氣采出程度影響曲線圖

4 實例井應用

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

表2 實例井基礎參數表

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

5 結論

1）在相同參數下，采用自擴散流熱耦合模型計算的頁巖氣產量大于修正達西模型計算值。

2）考慮溫度場變化時，自擴散系數剖面出現峰值，氣體密度剖面出現谷值，并且峰值/谷值對應數據點同步向地層內部移動，說明自擴散系數影響著氣體傳質速率，且隨生產進行，近井低溫對氣體自擴散的影響范圍不斷增大。

圖10 長寧區塊五峰組—龍馬溪組頁巖氣水平井生產模擬結果圖

3）當井底溫度低于地層溫度時，近井氣體自擴散系數減小，發生氣體近井阻塞現象，從而使氣井產量降低。

4）根據現場實例井的生產模擬結果，考慮溫度場的變化，自擴散流熱耦合模型能夠更準確地模擬深層頁巖氣井產量。