吳大衛(wèi)，邸 元，孫建芳

(1.北京大學(xué) 工學(xué)院，北京 100871；2.中國石油化工股份有限公司石油勘探開發(fā)研究院，北京 100083)

在油氣開采過程中，孔隙流體壓力的變化會(huì)引起巖石骨架應(yīng)力的變化，從而導(dǎo)致油藏物性參數(shù)的改變，而這些改變又反過來影響孔隙流體的流動(dòng)和壓力分布[1-5]。對(duì)裂縫性油藏或低滲油藏等對(duì)應(yīng)力較為敏感的油藏而言，這一效應(yīng)更為明顯[6-11]，因此需要考慮流固耦合進(jìn)行計(jì)算。Terzaghi[12]引入有效應(yīng)力的概念來研究流固耦合問題，提出了土的一維固結(jié)模型。Biot[13-14]建立了三維固結(jié)理論，并將其理論推廣到了各向異性的多孔介質(zhì)。Geertsma[15]與Verruijt[16]對(duì)Biot理論進(jìn)行了推廣，進(jìn)一步發(fā)展了多相滲流與孔隙介質(zhì)耦合作用的理論模型。Savage等[17]將Biot理論應(yīng)用于橫觀各向同性的孔隙介質(zhì)中。Zienkiewicz等[18]考慮了固相材料非線性的問題，在Biot的三維固結(jié)方程基礎(chǔ)上提出了廣義Biot公式。

針對(duì)流固耦合問題的數(shù)值模擬研究大多采用有限元方法。然而采用有限元方法離散時(shí)，位移和孔隙壓力不宜采用相同的插值函數(shù)，從而導(dǎo)致編程實(shí)現(xiàn)較為復(fù)雜。Akai等[19]提出了一種有限元-有限差分混合方法來求解飽和多孔介質(zhì)的Biot方程。Oka等[20]則將這一方法應(yīng)用到了飽和土液化的數(shù)值分析。在此基礎(chǔ)上，本研究建立了一種流固耦合分析的有限元-有限體積混合方法。這一方法對(duì)平衡方程與滲流方程分別采用不同的數(shù)值離散方法，并可以對(duì)兩個(gè)方程使用不同密度的計(jì)算網(wǎng)格。通過對(duì)一維固結(jié)問題和二維單井定壓生產(chǎn)問題的數(shù)值模擬，驗(yàn)證了本文方法的有效性。

1 數(shù)學(xué)模型的建立

1.1 力學(xué)平衡方程

對(duì)固相和液相分別建立力的平衡方程：

(1)

(2)

引入Zienkiewicz的u-p假定，即認(rèn)為液相對(duì)于固相的加速度很小，代入各相的平衡方程中,并引入各相應(yīng)力與孔隙壓力的關(guān)系，可得：

(3)

式(3)即為流固耦合問題中考慮了固相位移的達(dá)西定律表達(dá)式，可見在考慮了固相位移后達(dá)西定律中出現(xiàn)了固相位移的相關(guān)項(xiàng)，這也導(dǎo)致了離散后的平衡方程和連續(xù)性方程中會(huì)出現(xiàn)耦合項(xiàng)。

將固相與液相的平衡方程相加，可以得到總體的平衡方程：

(4)

1.2 連續(xù)性方程

根據(jù)雷諾輸運(yùn)定理，可以得到固相和液相的質(zhì)量守恒方程，分別有：

(5)

(6)

將兩方程相加，并進(jìn)行化簡(jiǎn)后可以得到：

(7)

代入考慮固體骨架位移后的達(dá)西定律表達(dá)式式(3)，假定固體介質(zhì)不可壓縮，并引入液相的壓縮系數(shù)C，整理后可得滲流連續(xù)性方程：

(8)

2 離散方程的建立

分別對(duì)上述平衡方程和連續(xù)性方程采用有限元法和有限體積法進(jìn)行離散。

對(duì)力學(xué)平衡方程進(jìn)行有限元法離散，可以得到在空間域內(nèi)離散后的平衡方程：

(9)

式(9)中：[M]為質(zhì)量矩陣;{un}為節(jié)點(diǎn)位移;[K]為剛度矩陣;{KV}為體剛度向量;{Fd}為荷載相對(duì)增量;{Rd}|t為有效應(yīng)力相對(duì)增量。

對(duì)式(9)在時(shí)間域內(nèi)采用Newmarkβ方法離散，從而可得到有限元法離散后的力平衡方程。

對(duì)上述流固耦合問題的連續(xù)性方程采用有限體積法進(jìn)行離散。在單元內(nèi)對(duì)連續(xù)性方程取積分弱形式，并引入高斯定理，可得：

(10)

圖1 有限體積法中相鄰的2個(gè)單元Fig.1 Two adjacent elements in FVM

式(10)中：pd為孔隙壓力相對(duì)初始狀態(tài)的增量；ni為垂直于邊界的方向向量；V為積分區(qū)域體積；S為積分區(qū)域表面積。

對(duì)式(10)等號(hào)左側(cè)的第一項(xiàng)仍采用有限元法離散，對(duì)最后一項(xiàng)采用有限體積法離散，對(duì)于圖1的平面四邊形網(wǎng)格情形，該項(xiàng)可寫為：

(11)

式(11)中：vx i、vy i分別為2個(gè)方向的流入速度；nx i、ny i分別為方向向量的分量；bi為單元的第i條邊。

采用Newmarkβ方法在時(shí)間域離散，并對(duì)變量進(jìn)行簡(jiǎn)化：

(12)

(13)

(14)

(15)

式(12)～(15)中：t為時(shí)間坐標(biāo)，Δt為時(shí)間步長。

離散后可得連續(xù)性方程為：

(16)

式(16)中：t與t+Δt代表變量的時(shí)間步取值。

3 儲(chǔ)層區(qū)域采用精細(xì)網(wǎng)格

流固耦合問題的數(shù)值分析中，力的平衡方程與滲流方程的計(jì)算區(qū)域往往是不同的，滲流方程的計(jì)算區(qū)域?qū)ζ胶夥匠痰挠?jì)算區(qū)域而言只是其一個(gè)子區(qū)域[21]。以油藏?cái)?shù)值模擬為例，真正需要進(jìn)行滲流方程求解的只有儲(chǔ)油層的區(qū)域；而力平衡方程的計(jì)算應(yīng)當(dāng)在包含整個(gè)儲(chǔ)層的更大區(qū)域內(nèi)進(jìn)行，這樣才能保證對(duì)應(yīng)力狀態(tài)進(jìn)行正確的求解。

圖2 有限元法的網(wǎng)格和有限體積法的單元Fig.2 Grids in FEM and elements in FVM

有限元-有限體積混合方法對(duì)采用有限元法求解力平衡方程的區(qū)域和有限體積法求解滲流方程的區(qū)域，使用不同精細(xì)程度的網(wǎng)格單元。具體實(shí)現(xiàn)方法為：在整個(gè)計(jì)算區(qū)域進(jìn)行有限元的網(wǎng)格劃分；而在滲流場(chǎng)涉及的儲(chǔ)層區(qū)域，將有限元的網(wǎng)格進(jìn)一步細(xì)分，即將每個(gè)有限元中的單元離散為多個(gè)更小的單元作為有限體積法中使用的單元。

下面以儲(chǔ)層區(qū)域平面四邊形網(wǎng)格為例進(jìn)行說明。如圖2所示，實(shí)線表示有限元法中所使用的網(wǎng)格，而虛線表示的是有限體積法中使用的細(xì)分的子單元。圖2中節(jié)點(diǎn)5～8圍成的子單元，其節(jié)點(diǎn)位移向量可以通過其所在有限元單元的節(jié)點(diǎn)位移插值函數(shù)得到：

{u}={u5,u6,u7,u8}T=[N′]{un}。

(17)

式(13)中：{u}為子單元位移向量；{un}為單元節(jié)點(diǎn)位移向量；[N′]為插值函數(shù)組成的轉(zhuǎn)換矩陣。

于是，可以分別得到有限元單元的平衡方程與子單元的連續(xù)性方程如下：

(18)

(19)

式(18)中：pavg為單元的平均孔隙壓力。可見，對(duì)于每個(gè)有限元單元，方程由1個(gè)力平衡方程和9個(gè)連續(xù)性方程組成，對(duì)拼裝后的總體方程組求解，就可得到整個(gè)問題的計(jì)算結(jié)果。

流固耦合問題的計(jì)算中，需要考慮應(yīng)力或應(yīng)變對(duì)固體介質(zhì)孔隙度和滲透率的影響。本文采用如下的經(jīng)驗(yàn)公式[22]表示孔隙度及滲透率與應(yīng)力的關(guān)系：

(20)

(21)

4 算例分析

4.1 一維固結(jié)問題的計(jì)算

圖3 一維固結(jié)問題數(shù)值解與理論解Fig.3 Analytical and numerical solutions for 1D consolidation problem

首先采用一維彈性固結(jié)問題的算例來驗(yàn)證本文方法的準(zhǔn)確性。計(jì)算一塊1 m×1 m的巖體，設(shè)定彈性模量E為1 000 kPa，泊松比ν為0，孔隙比e為0.428 6，載荷為頂部均布載荷q為200 kPa，設(shè)定頂部為排水邊界，其他3個(gè)邊界不排水。

圖3給出了該巖體豎向固結(jié)位移隨固結(jié)時(shí)間變化的數(shù)值模擬結(jié)果，與Terzaghi的一維固結(jié)問題理論解對(duì)比，二者吻合很好，說明了本文有限元-有限體積混合方法在求解流固耦合問題的準(zhǔn)確性。

對(duì)于一維固結(jié)問題，采用非結(jié)構(gòu)平面四邊形網(wǎng)格劃分的算例進(jìn)行計(jì)算，如圖4所示。在此算例中水平與豎直尺寸分別為5、10 m，2個(gè)側(cè)邊及底邊均為不排水邊界，頂部為排水邊界，施加荷載為頂部均布荷載。巖體材料彈性模量E為20 GPa，泊松比ν為0.2，孔隙比e為0.4，滲透率k為1 Darcy，水的黏度μ為1×10-3Pa·s，靜止土壓力系數(shù)K0為0.8，材料密度ρ為1.71×103kg/m3。

以計(jì)算區(qū)域頂端的中點(diǎn)位移為固結(jié)沉降的代表值，改變頂部荷載的大小，考察土體最終豎向固結(jié)沉降位移與施加荷載大小的關(guān)系，并與理論解進(jìn)行對(duì)比。如圖5所示，數(shù)值計(jì)算結(jié)果與理論解吻合較好。

圖4 一維固結(jié)問題算例的網(wǎng)格Fig.4 Grids used in 1D consolidation problem

圖5 最終固結(jié)沉降數(shù)值解與理論解Fig.5 Analytical and numerical solutions for final consolidation displacements

4.2 單井生產(chǎn)問題的模擬

對(duì)二維單相流體儲(chǔ)層的單井定壓生產(chǎn)問題進(jìn)行模擬，取計(jì)算區(qū)域?yàn)橐粔K長5 000 m、深4 000 m向上延伸至地表的區(qū)域，儲(chǔ)層位于地下3 000 m處，儲(chǔ)層厚度約為50 m。采用非結(jié)構(gòu)網(wǎng)格對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行網(wǎng)格劃分，如圖6所示。計(jì)算時(shí)約束了兩側(cè)的豎向、底邊的豎向及水平的位移。

僅設(shè)有一口生產(chǎn)井，生產(chǎn)井位于儲(chǔ)層的正中部，進(jìn)行定壓生產(chǎn)。算例中對(duì)儲(chǔ)層內(nèi)與儲(chǔ)層外的區(qū)域采用不同的材料參數(shù)。對(duì)于儲(chǔ)層區(qū)域，彈性模量E1為20 GPa，泊松比ν1為0.2，孔隙比e1為0.4，滲透率k1為1 Darcy，密度ρ1為1.71×103kg/m3；對(duì)于儲(chǔ)層以外的區(qū)域，彈性模量E2為20 GPa，泊松比ν2為0.2，孔隙比e2為0.1，滲透率k2為1×103Darcy，密度ρ2為1.91×103kg/m3。儲(chǔ)層中液相的密度ρf為1.00×103kg/m3，液相的黏度μ為1×10-3Pa·s，靜止土壓力系數(shù)K0為0.8。

隨著井定壓生產(chǎn)的進(jìn)行，孔隙壓力的降低起初只在生產(chǎn)井的周圍出現(xiàn)，最終逐漸發(fā)展至整個(gè)儲(chǔ)層區(qū)域，圖7給出了生產(chǎn)井定壓生產(chǎn)1個(gè)月之后儲(chǔ)層孔隙壓力變化(P-P0)的分布。

圖6 單井生產(chǎn)問題的有限元網(wǎng)格Fig.6 Grids used for FEM in single well problem

圖7 生產(chǎn)1個(gè)月后孔隙壓力變化的分布Fig.7 Pore pressure change distribution after one month of production

圖8給出了生產(chǎn)1個(gè)月后，計(jì)算區(qū)域內(nèi)豎向位移的分布情況?？梢婋S著生產(chǎn)井定壓生產(chǎn)，整個(gè)計(jì)算區(qū)域出現(xiàn)了明顯的變形，而且井單元處豎向位移顯著，這反映了實(shí)際生產(chǎn)中的井沉降現(xiàn)象?？紫抖群蜐B透率對(duì)應(yīng)力的敏感性對(duì)流固耦合問題的計(jì)算結(jié)果有著很大的影響。圖9給出了考慮應(yīng)力對(duì)孔隙度與滲透率的影響時(shí)的生產(chǎn)井產(chǎn)量，可見產(chǎn)量出現(xiàn)了明顯的減少。

圖8 生產(chǎn)1個(gè)月后的豎向位移分布Fig.8 Vertical displacement after one month of production

圖9 考慮應(yīng)力對(duì)孔滲系數(shù)影響時(shí)的產(chǎn)量Fig.9 Productions calculated using fixed and stress-dependent porosity and permeability

5 結(jié) 論

基于多孔介質(zhì)力學(xué)、滲流力學(xué)及u-p假定，建立了滲流與巖體變形耦合問題的基本方程，并采用一種有限元-有限體積混合方法對(duì)方程進(jìn)行離散和求解，此方法對(duì)平衡方程與連續(xù)性方程分別采用不同的數(shù)值離散方法和不同精細(xì)程度的網(wǎng)格進(jìn)行計(jì)算，既能保證滲流區(qū)域的計(jì)算精度，又能夠提高計(jì)算效率。

一維固結(jié)問題的數(shù)值計(jì)算結(jié)果與理論解吻合較好，驗(yàn)證了該方法的準(zhǔn)確性。對(duì)單井定壓生產(chǎn)問題的算例進(jìn)行了模擬，給出了生產(chǎn)過程中孔隙壓力與豎向變形的分布，計(jì)算結(jié)果表明，當(dāng)考慮孔隙度和滲透率隨應(yīng)力變化時(shí)，定壓生產(chǎn)井的產(chǎn)量有明顯的減少。

[1] 周志軍.低滲透儲(chǔ)層流固耦合滲流理論及應(yīng)用研究[D].大慶：東北石油大學(xué),2003.

[2] 田杰.巖體滲流的流固耦合問題及其工程應(yīng)用[D].北京：中國科學(xué)院研究生院(滲流流體力學(xué)研究所),2005.

[3] 董平川,郎兆新,徐小荷,等.油井開采過程中油層變形的流固耦合分析[J].地質(zhì)力學(xué)學(xué)報(bào),2000,6(2):6.

[4] EL-AMIN M F, NEGARA A, SALAMA A, et al. Simulation of coupled flow and mechanical deformation using implicit pressure-displacement explicit saturation (IMPDES) scheme[C]//SPE Middle East Unconventional Gas Conference and Exhibition. Abu Dhabi: SPE,2012.

[5] ALPAK F O. Robust fully-implicit coupled multiphase-flow and geomechanics simulation[J]. SPE Journal,2015,20(6).

[6] 劉建軍,裴桂紅.裂縫性低滲透油藏流固耦合滲流分析[J].應(yīng)用力學(xué)學(xué)報(bào),2004,21(1):36.

[7] LI S, LI X, ZHANGhang D. A fully coupled thermo-hydro-mechanical, three-dimensional model for hydraulic stimulation treatments[J]. Journal of Natural Gas Science and Engineering,2016,34:64.

[8] BLANCO-MARTIN L, RUTQVIST J, DOUGHTY C, et al. Coupled geomechanics and flow modeling of thermally induced compaction in heavy oil diatomite reservoirs under cyclic steaming[J]. Journal of Petroleum Science and Engineering,2016,147:474.

[9] REN G T, JIANG J M, YOUNIS R. Fully-coupled XFEM-EDFM hybrid model for geomechanics and flow in fractured reservoirs[C]//SPE Reservoir Simulation Conference. Houston, Texas: SPE,2017.

[10] ANDERSON O, NILSON H M, RAYNAUD X. Coupled geomechanics and flow simulation on corner-point and polyhedral grids[C]//SPE Reservoir Simulation Conference. Houston, Texas: SPE,2017.

[11] EKER E, UZUN I, KAZEMI H. Numerical simulation of dual-porosity multiphase flow using poroelasticity theory[C]//SPE Reservoir Simulation Conference. Houston, Texas: SPE,2017.

[12] TERZAGHI K. Theoretical soil mechanics[M]. London: Chapman and Hall,1959.

[13] BIOT M A. General theory of three-dimensional consolidation[J]. Journal of Applied Physics,1941,12(2):155.

[14] BIOT M A. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range[J]. The Journal of the Acoustical Society of America,1956,28(2),168.

[15] GEERTZMA J. The effect of fluid pressure decline on volumetric changes of porous rocks[J]. Transactions of the American Institute of Mining Engineers,1957,210(12):331.

[16] VERRUIJT A. Elastic storage of aquifers[M]//Flow through porous media. New York: Academic Press,1969.

[17] SAVAGE W Z, BRADDOCK W A. A model for hydrostatic consolidation of Pierre shale[J]. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts,1991,28(5):345.

[18] ZIENKIEWICZ O C, SHIOMI T. Dynamic behaviour of saturated porous media: the generalized Biot formulation and its numerical solution [J]. International Journal for Numerical and Analytical Methods in Geomechanics,1984,8:71.

[19] AKAI K and TAMURA T. Study of two-dimensional consolidation accompanied by an elasto-plastic constitutive equation [J]. Proceedings of Japan Society of Civil Engineers,1978,269:98.

[20] OKA F, YASHIMA A, SHIBATA T, et al. FEM-FDM coupled liquefaction analysis of a porous soil using an elasto-plastic model[J]. Applied Scientific Research,1994,52(3):209.

[21] MINKOFF S E, STONE C M, BRYANT S, et al. Coupled fluid flow and geomechanical deformation modeling[J]. Journal of Petroleum Science and Engineering,2003,38(1/2):37.

[22] RUTQVIST J, WU Y S, TSANG C F, et al. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock[J]. International Journal of Rock Mechanics and Mining Sciences,2002,39(4):435,437.