異常高壓氣藏流固耦合的數(shù)值模擬

鄭琴 陶自強 孟立新

中國石油大港油田公司勘探開發(fā)研究院

異常高壓氣藏在開發(fā)過程中，隨著氣藏壓力的下降，氣藏巖石骨架要承受比常規(guī)氣藏大得多的凈上覆壓力，儲層應(yīng)力敏感性極強，致使異常高壓氣藏具有比常規(guī)氣藏更強的流固耦合效應(yīng)。因此，研究異常高壓氣藏流固耦合效應(yīng)的理論及方法，對于真實模擬氣藏開采，指導(dǎo)氣田生產(chǎn)具有十分重要的現(xiàn)實意義。

縱觀國內(nèi)外流固耦合理論的研究現(xiàn)狀及發(fā)展動態(tài)發(fā)現(xiàn)，流固耦合已由簡單的單相孔隙介質(zhì)模型向復(fù)雜的兩相或三相連續(xù)介質(zhì)及擬連續(xù)或非連續(xù)介質(zhì)模型發(fā)展［1－4］。在石油工業(yè)領(lǐng)域，流固耦合理論在解決油藏及常規(guī)氣藏流體滲流及巖石小變形相互作用方面的研究較多，并且取得了較好的開發(fā)效果［5－8］，而對于異常高壓氣藏，多數(shù)參考文獻(xiàn)主要研究其開發(fā)動態(tài)及儲層物性變化，而對流固耦合理論在異常高壓氣藏開發(fā)機理上的研究很少報道。目前的流固耦合模型大多采用有限差分或有限差分與有限元法聯(lián)合求解［9－11］，在處理復(fù)雜形狀油氣藏及復(fù)雜邊界條件等問題上存在一定局限性，采用有限元法統(tǒng)一離散求解流固耦合模型的研究工作目前開展較少。

1 孔隙度、滲透率與體積應(yīng)變的關(guān)系

孔隙度與巖石變形的關(guān)系：

滲透率與孔隙度、體積應(yīng)變的關(guān)系為［12］：

2 異常高壓氣藏流固耦合模型的建立

2．1 滲流場方程的建立

2．1．1 運動方程

2．1．2 連續(xù)性方程

固體骨架的體積應(yīng)變用εV＝εx＋εy＋εz＝ui，i表示，經(jīng)整理可得最終整體連續(xù)性方程：

2．1．3 狀態(tài)方程

氣體密度的變化規(guī)律為：

2．1．4 滲流微分方程

假設(shè)是由于不可壓縮的固相巖土顆粒發(fā)生相對位移后重新組合造成了固體巖石的變形，故固體密度ρs為常數(shù)，ρs對時間的偏導(dǎo)數(shù)為零，可對式（4）進行簡化，再由式（3）可得異常高壓氣藏滲流微分方程：

2．2 應(yīng)力場方程的建立

2．2．1 平衡方程

2．2．2 應(yīng)變—位移方程（幾何方程）

可寫為矩陣形式ε＝Lu。

2．2．3 本構(gòu)方程

1）線彈性本構(gòu)方程

線彈性體的本構(gòu)關(guān)系是廣義胡克定律，其彈性關(guān)系（即應(yīng)力—應(yīng)變關(guān)系）為：

也可寫為矩陣形式σ′ij＝Deεij。

De在三維情況下，它是含有36個元素的對稱矩陣。

2）彈塑性本構(gòu)方程

彈塑性應(yīng)力—應(yīng)變本構(gòu)關(guān)系常用增量形式表示如下：

彈塑性矩陣［Dep］＝［De］－［Dp］，［Dp］表示塑性矩陣。

由彈塑性關(guān)聯(lián)流動法則可以推導(dǎo)出塑性矩陣的表達(dá)式：

2．2．4 修正的應(yīng)力—應(yīng)變關(guān)系

Biot提出含水多孔體的應(yīng)力—應(yīng)變關(guān)系應(yīng)該寫成：

其中α被稱為Biot耦合系數(shù)，δij為Kroneker符號，它定義為：δij＝0，i≠j；δij＝1，i＝j(luò)。

2．2．5 應(yīng)力場方程

由平衡方程（7），線彈性本構(gòu)方程（9）以及（12）給出的應(yīng)力—應(yīng)變關(guān)系，先消去應(yīng)力σij，再消去εij，可得用位移量ux、uy、uz表示的3個方程：

其中的變量為3個位移量和壓力。

2．3 定解條件的確定

2．3．1 滲流場控制方程的定解條件

初始條件：t＝0，p＝pi。

封閉外邊界條件：r＝re，p／n＝0。

定流量內(nèi)邊界條件：r＝rw，p／n＝fq（x，y，z，t）。

2．3．2 應(yīng)力場控制方程的定解條件

設(shè)巖石骨架所占空間區(qū)域為Ωd（d＝2或3），域邊界為Ω，如果設(shè)Γu為位移邊界，Γσ為應(yīng)力邊界，則ΩΓu∪Γσ，Γu∩Γσ＝0，邊界條件如下。

位移邊界條件：ui｜Γu＝珔ui。

應(yīng)力邊界條件：σijnj｜Γσ＝珡Ti。

初始位移條件：ui｜i＝0＝0。

3 流固耦合數(shù)學(xué)模型有限元離散

3．1 空間8節(jié)點等參元

采用空間8節(jié)點等參元對求解域進行離散。坐標(biāo)變換函數(shù)式和位移函數(shù)采用統(tǒng)一形式

將8個不同的形函數(shù)統(tǒng)一表示成：

利用雅可比矩陣建立整體坐標(biāo)（x，y，z）和局部坐標(biāo)（ξ，η，ζ）之間的一一對應(yīng)關(guān)系。即

在建立網(wǎng)格時，由于地層具有對稱性，只需要對四分之一個區(qū)域進行研究（圖1）。

圖1 空間8節(jié)點單元網(wǎng)格圖

3．2 流固耦合應(yīng)力場方程有限元離散

對8節(jié)點等參元，有u＝Na、p＝N珚p（珚p＝［p1，p2，…，p8］），由虛功原理，經(jīng)化簡可得：

其中，三維問題的m＝［1 1 1 0 0 0］T，二維問題的m＝［1 1 0］T。對于彈塑性形變和塑性形變，只需將上式中的彈性本構(gòu)矩陣D換為彈塑性本構(gòu)矩陣［Dep］和塑性矩陣［Dp］即可。單元形式的矩陣方程為：

增量形式的單元矩陣方程有：

3．3 流固耦合滲流場方程有限元離散

3．3．1 方程的空間離散

對方程（6）進行化簡可得：

取權(quán)函數(shù)為形函數(shù)N，且完全滿足邊界條件，通過Galerkin積分，經(jīng)變換得：

式（20）可簡寫為：

3．3．2 方程的時域離散［13］

空間離散后的流固耦合滲流微分方程含有函數(shù)、對時間的一次微分項，需進一步時間上離散。最終可得單元有限元方程：

把各單元的Ae、Be、Ce等都加以擴大到整個結(jié)構(gòu)的自由度的維數(shù)，然后疊加可得到整體系數(shù)矩陣，并將單元矩陣方程組裝成整體矩陣方程。

3．4 系數(shù)矩陣的確定

求解流固耦合有限元方程時，需確定有關(guān)單元系數(shù)矩陣和荷載矩陣的值，這類積分表達(dá)式中的被積函數(shù)通常包含了形函數(shù)及其偏導(dǎo)數(shù)的復(fù)雜組合，因此筆者采用高斯積分法來確定系數(shù)矩陣。

4 流固耦合模型算例

設(shè)原始地層壓力為60MPa，初始孔隙度為10%，初始滲透率0．03D，氣體原始壓縮系數(shù)為1．5 MPa－1，黏度為0．35mPa·s，氣藏巖石彈性模量為300MPa，泊松比為0．3，Biot系數(shù)為1．0，氣層有效厚度30m，產(chǎn)量為1×104m3／d，井筒半徑為0．1m，地層半徑為1 000m。應(yīng)力模型的邊界條件為：①氣藏外邊界的垂直位移設(shè)置為零；②氣藏內(nèi)邊界的水平位移設(shè)置為零；③氣藏頂部邊界設(shè)置為恒定的應(yīng)力邊界，其垂直應(yīng)力為65MPa，水平應(yīng)力為35MPa。

4．1 流體壓力比較

用非耦合模型和耦合模型兩種情況進行預(yù)測。求解有限元方程，得到兩個模型井底處某點的壓力隨生產(chǎn)時間變化的曲線圖（圖2）。

圖2 壓力隨時間變化趨勢圖

圖3是生產(chǎn)進行到300d時，井底所在平面上壓力沿徑向的三維分布圖。由圖中可看出近井地帶壓力的衰減比遠(yuǎn)井地帶要快，形成了 “壓降漏斗”。

圖3 生產(chǎn)300 d后壓力隨徑向距離變化趨勢圖

4．2 儲層位移比較

圖4、5表示不同時刻水平位移及垂直位移與徑向距離的關(guān)系。

圖4 不同時刻水平位移與徑向距離的關(guān)系圖

圖5 不同時刻垂直位移與徑向距離的關(guān)系圖

5 結(jié)論

1）如圖2所示，氣藏采用衰竭式開采，孔隙壓力逐漸減小，且非耦合模型的孔隙壓力下降比耦合模型要快。因此，氣藏滲流與應(yīng)力的耦合效應(yīng)對氣藏開發(fā)動態(tài)的預(yù)測有較大的影響，是不容忽視的因素。

2）由圖3可知氣藏應(yīng)力隨空間的變化規(guī)律是：井眼附近變化量較大，變化梯度也明顯較大，向油氣藏邊界很快減小。

3）從圖4、5中可以看出，只有在離井很近的地方，才會有明顯的位移，離井越遠(yuǎn)，位移的幅度也就越小。

3）筆者給出了異常高壓氣藏流固耦合模型建立及數(shù)值求解的方法，通過求解流固耦合模型可以得到流體壓力和巖石質(zhì)點位移的變化規(guī)律，實例分析結(jié)果與工程實際相符合。

符 號 說 明

φ為孔隙度，%；φ0為初始孔隙度，%；εV為巖石體積應(yīng)變，εV＝εx＋εy＋εz；K為滲透率，D；K0為初始滲透率，D；珗vg為氣體滲流速度，m／h；為氣體相對于固體的速度，m／h；μ為天然氣黏度，mPa·s；p為壓力，MPa；ρg為氣體密度，kg／m3；t為時間，h；ρs為固體密度，kg／m3；uii為巖石體積應(yīng)變，m；ρg0為初始?xì)怏w密度，kg／m3；cg為氣體壓縮系數(shù)，MPa－1；pi為初始壓力，MPa；σij為總應(yīng)力張量分量，MPa；fi為體力分量，MPa；εij為體積應(yīng)變張量分量，m；ui為巖石質(zhì)點的位移分量，m；σ′為固體骨架的有效應(yīng)力，MPa；λ為拉梅系數(shù)；G為剪切模量，MPa；γij為剪應(yīng)變張量分量，m；dσ為有效應(yīng)力增量，MPa；dε為應(yīng)變增量，m；［Dep］為彈塑性矩陣；a為Biot耦合系數(shù)，無量綱；F為屈服函數(shù)；A為硬化指數(shù)；Ni為形函數(shù)；ξ、η、ζ分別為局部坐標(biāo)系中的變量；ξi、ηi、ζi分別為局部坐標(biāo)系中節(jié)點i的坐標(biāo)分量；J－1為三階Jacobi矩陣的逆；L為應(yīng)變位移矩陣；D為彈性本構(gòu)矩陣；V巖石單元體積，m3；a為單元節(jié)點位移矩陣；f為單元體積力，MPa；T為單元表面力，MPa；Γ為巖石單元邊界；ae為單元節(jié)點位移矩陣；珚pe為單元節(jié)點應(yīng)力矩陣；Ke為單元剛度矩陣；Qe為孔隙壓力載荷矩陣；為單元等效節(jié)點荷載矩陣；Δae為單元節(jié)點位移增量；Δpe為單元等效孔壓增量；為n＋1時間點的單元節(jié)點應(yīng)力矩陣為n＋1時間點的單元應(yīng)變矩陣；為n時間點的單元應(yīng)變矩陣為n時間點的單元節(jié)點應(yīng)力矩陣；Δt為時間步長，h。

［1］TERZAGHI K．Theoretical soil mechanics［M］．New York：Tihn Willey，1943．

［2］陽仿勇．變形介質(zhì)氣藏流固耦合滲流理論及應(yīng)用研究［D］．成都：西南石油大學(xué)，2005：22－58．YANG Fangyong．Theory and application research of fluidsolid coupling seepage for deformation medium gas reservoir［D］．Chengdu：Southwest Petroleum University，2005：22－58．

［3］OSORIO J G，CHEN Heryuan．Numerical simulation of the impact of flow－induced geomechanical response on the production of stress－sensitive reservoirs［C］∥paper 51929－MS presented at the SPE Annual Technical Conference and Exhibition，14－17February 1999，Houston，Texas，USA．New York：SPE，1999．

［4］TORTILE W S，F(xiàn)AROUQ ALI S M．A framework for multiphase nonisothermal fluid flow in a deforming heavy oil reservoir［C］∥paper 16030－MS presented at the SPE Annual Technical Conference and Exhibition，1－4February 1987，San Antonio，Texas，USA．New York：SPE，1987．

［5］LEWIS R W，SUKIRMAN Y．Finite element modeling of three－phase flow in deforming saturated oil reservoirs［J］．International Journal for Numerical and Analytical Methods in Geomecharics，1993，17（8）：577－598．

［6］文成楊，杜志敏，唐石丹，等．裂縫—孔隙型碳酸鹽巖氣藏流固耦合新模型對比［J］．天然氣工業(yè)，2008，28（8）：89－91．WEN Chengyang，DU Zhimin，TANG Shidan，et al．A comparison of new fluid－solid coupling models for fracturepore type carbonate gas reservoirs［J］．Natural Gas Industry，2008，28（8）：89－91．

［7］雷霆，李治平．低滲氣藏變形效應(yīng)的處理方法和合理生產(chǎn)壓差的選擇［J］．天然氣工業(yè)，2007，27（2）：93－94．LEI Ting，LI Zhiping．Treatment of deformation effect and selection of reasonable production pressures differential for low－permeability reservoirs［J］．Natural Gas Industry，2007，27（2）：93－94．

［8］MINKOFF S E，STONE C M．Staggered in time coupling of reservoir flow simulation and geomechanical deformation：Step l － one－way coupling［C］∥paper 51920－MS presented at the SPE Annual Technical Conference and Exhibition，14－17February 1999，Houston，Texas，USA．New York：SPE，1999．

［9］FUNG L S K，BUCHANAN L．Coupled geomechanicalthermal simulation for deforming heavy－oil reservoirs［J］．Journal of Canadian Petroleum Technology，1994，33（4）：22－28．

［10］冉啟全，顧小蕓．彈塑性變形油藏中多相滲流的數(shù)值模擬［J］．計算力學(xué)學(xué)報，1999，16（1）：24－31．RAN Qiquan，GU Xiaoyun．Numerical simulation of multiphase flow in an elastoplastic deforming oil reservoir［J］．Chinese Journal of Computational Mechanics，1999，16（1）：24－31．

［11］范學(xué)平，李秀生，張士誠，等．低滲透氣藏整體壓裂流固耦合數(shù)學(xué)模擬［J］．石油勘探與開發(fā)，2000，27（1）：76－79．FAN Xueping，LI Xiusheng，ZHANG Shicheng，et al．Mathematical simulation of coupled fluid flow and geomechanical behavior for full low permeability gas reservoir fracturing［J］．Petroleum Exploration and Development，2000，27（1）：76－79．

［12］冉啟全，李士倫．流固耦合油藏數(shù)值模擬中物性參數(shù)動態(tài)模型研究［J］．石油勘探與開發(fā)，1997，24（3）：61－65．RAN Qiquan，LI Shilun．Study on dynamic models of reservoir parameters in the coupled simulation of multiphase flow and reservoir deformation［J］．Petroleum Exploration and Development，1997，24（3）：61－65．

［13］朱伯芳．有限單元法原理與應(yīng)用［M］．北京：中國水利水電出版社，1998 ZHU Bofang．Principle and application of finite element method［M］．Beijing：China Water Conservancy and Electricity Press，1998．