在油气开采过程中,孔隙流体压力的变化会引起岩石骨架应力的变化,从而导致油藏物性参数的改变,而这些改变又反过来影响孔隙流体的流动和压力分布^[1-5]。对裂缝性油藏或低渗油藏等对应力较为敏感的油藏而言,这一效应更为明显^[6-11],因此需要考虑流固耦合进行计算。Terzaghi^[12]引入有效应力的概念来研究流固耦合问题,提出了土的一维固结模型。Biot^[13-14]建立了三维固结理论,并将其理论推广到了各向异性的多孔介质。Geertsma^[15]与Verruijt^[16]对Biot理论进行了推广,进一步发展了多相渗流与孔隙介质耦合作用的理论模型。Savage等^[17]将Biot理论应用于横观各向同性的孔隙介质中。Zienkiewicz等^[18]考虑了固相材料非线性的问题,在Biot的三维固结方程基础上提出了广义Biot公式。

针对流固耦合问题的数值模拟研究大多采用有限元方法。然而采用有限元方法离散时,位移和孔隙压力不宜采用相同的插值函数,从而导致编程实现较为复杂。Akai等^[19]提出了一种有限元-有限差分混合方法来求解饱和多孔介质的Biot方程。Oka等^[20]则将这一方法应用到了饱和土液化的数值分析。在此基础上,本研究建立了一种流固耦合分析的有限元-有限体积混合方法。这一方法对平衡方程与渗流方程分别采用不同的数值离散方法,并可以对两个方程使用不同密度的计算网格。通过对一维固结问题和二维单井定压生产问题的数值模拟,验证了本文方法的有效性。

将多孔介质视为固相与液相的混合体,用角标s代表固相,用角标f代表液相。用n代表孔隙率,设均匀化后固相材料密度为(-overρ)^s=(1-n)ρ^s,液相密度为(-overρ)^f=n ρ^f。

对固相和液相分别建立力的平衡方程:

(-overρ)^s(¨overu)^s_i-R_i=(σ^s_ij)/(x_j)+(-overρ)^sb_i,(1)

(-overρ)^f(¨overu)^f_i+R_i=(σ^f_ij)/(x_j)+(-overρ)^fb_i。(2)

式(1)～(2)中:(¨overu)^s与(¨overu)^f分别为固相与液相加速度; σ^s_ij与σ^f_ij分别为固相与液相的应力张量; x_j为空间坐标; b_i为体积力; R_i为液相对固相的作用力。

引入Zienkiewicz的u -p假定,即认为液相对于固相的加速度很小,代入各相的平衡方程中,并引入各相应力与孔隙压力的关系,可得:

(·overω)_i=-k/μ((p)/(x_i)+ρ^f(¨overu)^s_i-ρ^fb_i)。(3)

式(3)中:(·overω)_i为液相对固相的相对速度; k为渗透率; μ为液相黏度; p为孔隙压力; x_i为空间坐标。

式(3)即为流固耦合问题中考虑了固相位移的达西定律表达式,可见在考虑了固相位移后达西定律中出现了固相位移的相关项,这也导致了离散后的平衡方程和连续性方程中会出现耦合项。

将固相与液相的平衡方程相加,可以得到总体的平衡方程:

ρ(¨overu)^s_i=(σ_ij)/(x_j)+ρ b_i。(4)

式(4)中:ρ=(-overρ)^s+(-overρ)^f=(1-n)ρ^s+nρ^f为混合体的平均密度。

根据雷诺输运定理,可以得到固相和液相的质量守恒方程,分别有:

((-overρ)^s)/(t)+(((-overρ)^s(·overu)^s_i))/(x_i)=0,(5)

((-overρ)^f)/(t)+(((-overρ)^f(·overu)^f_i))/(x_i)=0。(6)

将两方程相加,并进行化简后可以得到:

((·overω)_i)/(x_i)+(·overε)^s_ii+n((·overρ)^f)/(ρ^f)+(1-n)((·overρ)^s)/(ρ^s)=0。(7)

式(7)中:(·overε)^s_ii为固相的体应变变化率。

代入考虑固体骨架位移后的达西定律表达式式(3),假定固体介质不可压缩,并引入液相的压缩系数C,整理后可得渗流连续性方程:

k/μ(-ρ^f(¨overε)^s_ii-(²p)/(x_ix_i))+(·overε)^s_ii+nCp^·=0。(8)

分别对上述平衡方程和连续性方程采用有限元法和有限体积法进行离散。

对力学平衡方程进行有限元法离散,可以得到在空间域内离散后的平衡方程:

[M]{u^.._n}+[K]{Δu_n}-{K_V}p_d={F_d}-{R_d}_|t。(9)

式(9)中:[M]为质量矩阵; {u_n}为节点位移; [K]为刚度矩阵; {K_V}为体刚度向量; {F_d}为荷载相对增量; {R_d}_|t为有效应力相对增量。

对式(9)在时间域内采用Newmark β方法离散,从而可得到有限元法离散后的力平衡方程。

对上述流固耦合问题的连续性方程采用有限体积法进行离散。在单元内对连续性方程取积分弱形式,并引入高斯定理,可得:

∫_V(-k/μρ^f(¨overε)^s_ii+(·overε)^s_ii+n Cp^·_d)dV-∫_Sk/μ(p_d)/(x_i)n_idS=0。(10)

图1 有限体积法中相邻的2个单元
Fig.1 Two adjacent elements in FVM

式(10)中:p_d为孔隙压力相对初始状态的增量; n_i为垂直于边界的方向向量; V为积分区域体积; S为积分区域表面积。

对式(10)等号左侧的第一项仍采用有限元法离散,对最后一项采用有限体积法离散,对于图1的平面四边形网格情形,该项可写为:

∫_Sk/μ(p_d)/(x_i)n_idS=∑⁴_i=1∫_bi(v_{x i}n_{x i}+v_{y i}n_{y i})db。(11)

式(11)中:v_{x i}、v_{y i}分别为2个方向的流入速度; n_{x i}、n_{y i}分别为方向向量的分量; b_i为单元的第i条边。

采用Newmark β方法在时间域离散,并对变量进行简化:

φ=(ρ^fk)/μ,(12)

α'_i=k/((φ-γΔt)μ)((b_{y i})/(s _{x i})+(b_{x i})/(s_{y i})),(13)

α'=∑⁴_i=1α'_i,(14)

A'=1/(Δt(φ-γΔt))∫_VnCdV。(15)

式(12)～(15)中:t为时间坐标,Δt为时间步长。

离散后可得连续性方程为:

-{K_V}^T{u^.._n}_{|t +Δt}+(A'_|t-α')p_{d|t +Δt}+∑⁴_i=1α'_ip_{di|t +Δt}=

-1/(φ-γΔt){K_V}^T({u^._n}_|t+(1-γ)Δt{u^.._n}_|t)+A'_|tp_d|t。(16)

式(16)中:t与t+Δt代表变量的时间步取值。

流固耦合问题的数值分析中,力的平衡方程与渗流方程的计算区域往往是不同的,渗流方程的计算区域对平衡方程的计算区域而言只是其一个子区域^[21]。以油藏数值模拟为例,真正需要进行渗流方程求解的只有储油层的区域; 而力平衡方程的计算应当在包含整个储层的更大区域内进行,这样才能保证对应力状态进行正确的求解。

图2 有限元法的网格和有限体积法的单元
Fig.2 Grids in FEM and elements in FVM

有限元-有限体积混合方法对采用有限元法求解力平衡方程的区域和有限体积法求解渗流方程的区域,使用不同精细程度的网格单元。具体实现方法为:在整个计算区域进行有限元的网格划分; 而在渗流场涉及的储层区域,将有限元的网格进一步细分,即将每个有限元中的单元离散为多个更小的单元作为有限体积法中使用的单元。

下面以储层区域平面四边形网格为例进行说明。如图2所示,实线表示有限元法中所使用的网格,而虚线表示的是有限体积法中使用的细分的子单元。图2中节点5～8围成的子单元,其节点位移向量可以通过其所在有限元单元的节点位移插值函数得到:

{u}={u₅,u₆,u₇,u₈}^T=[N']{u_n}。(17)

式(13)中:{u}为子单元位移向量; {u_n}为单元节点位移向量; [N']为插值函数组成的转换矩阵。

于是,可以分别得到有限元单元的平衡方程与子单元的连续性方程如下:

([M]+β(Δt)²[K]_{|t+Δ t}){u^.._n}_{|t+Δ t}-{K_V}p_avg=

{F_d}_{|t +Δ t}-{R_d}_|t-[K]_{|t +Δ t}(Δt{u^._n}_|t+(1/2-β)(Δt)²{u^.._n}_|t),(18)

-{K_V}^T[N']{u^.._n}_{|t +Δ t}+(A'_|t-α')p_{di|t +Δ t}+∑⁴_j=1α'_ip_{dij|t +Δ t}=

-1/(φ-γΔt){K_V}^T[N']({u^._n}_|t+(1-γ)Δt{u^.._n}_|t)+A'_|tp_di|t。(19)

式(18)中:p_avg为单元的平均孔隙压力。可见,对于每个有限元单元,方程由1个力平衡方程和9个连续性方程组成,对拼装后的总体方程组求解,就可得到整个问题的计算结果。

流固耦合问题的计算中,需要考虑应力或应变对固体介质孔隙度和渗透率的影响。本文采用如下的经验公式^[22]表示孔隙度及渗透率与应力的关系:

n=n_r+(n₀-n_r)exp(a×σ'_M),(20)

k=k₀exp[c×(n/(n₀)-1)]。(21)

式(20)～(21)中:n₀为初始孔隙度; n_r为残余孔隙度; k₀为初始渗透率; σ'_M为有效体应力; a与c为实验测得的常数。

首先采用一维弹性固结问题的算例来验证本文方法的准确性。计算一块1 m×1 m的岩体,设定弹性模量E为1 000 kPa,泊松比ν为0,孔隙比e为0.428 6,载荷为顶部均布载荷q为200 kPa,设定顶部为排水边界,其他3个边界不排水。

图 3给出了该岩体竖向固结位移随固结时间变化的数值模拟结果,与Terzaghi的一维固结问题理论解对比,二者吻合很好,说明了本文有限元-有限体积混合方法在求解流固耦合问题的准确性。

对于一维固结问题,采用非结构平面四边形网格划分的算例进行计算,如图4所示。在此算例中水平与竖直尺寸分别为5、10 m,2个侧边及底边均为不排水边界,顶部为排水边界,施加荷载为顶部均布荷载。岩体材料弹性模量E为20 GPa,泊松比ν为0.2,孔隙比e为0.4,渗透率k为1 Darcy,水的黏度μ为1×10^-3 Pa·s,静止土压力系数K₀为0.8,材料密度ρ为1.71×10³ kg/m³。

图3 一维固结问题数值解与理论解
Fig.3 Analytical and numerical solutions for 1D consolidation problem

以计算区域顶端的中点位移为固结沉降的代表值,改变顶部荷载的大小,考察土体最终竖向固结沉降位移与施加荷载大小的关系,并与理论解进行对比。如图5所示,数值计算结果与理论解吻合较好。

图4 一维固结问题算例的网格
Fig.4 Grids used in 1D consolidation problem

图5 最终固结沉降数值解与理论解
Fig.5 Analytical and numerical solutions for final consolidation displacements

对二维单相流体储层的单井定压生产问题进行模拟,取计算区域为一块长5 000 m、深4 000 m向上延伸至地表的区域,储层位于地下3 000 m处,储层厚度约为50 m。采用非结构网格对整个计算区域进行网格划分,如图6所示。计算时约束了两侧的竖向、底边的竖向及水平的位移。

仅设有一口生产井,生产井位于储层的正中部,进行定压生产。算例中对储层内与储层外的区域采用不同的材料参数。对于储层区域,弹性模量E₁为20 GPa,泊松比ν₁为0.2,孔隙比e₁为0.4,渗透率k₁为1 Darcy,密度ρ₁为1.71×10³ kg/m³; 对于储层以外的区域,弹性模量E₂为20 GPa,泊松比ν₂为0.2,孔隙比e₂为0.1,渗透率k₂为1×10³ Darcy,密度ρ₂为1.91×10³ kg/m³。储层中液相的密度ρ^f为1.00×10³ kg/m³,液相的黏度μ为1×10^-3 Pa·s,静止土压力系数K₀为0.8。

随着井定压生产的进行,孔隙压力的降低起初只在生产井的周围出现,最终逐渐发展至整个储层区域,图7给出了生产井定压生产1个月之后储层孔隙压力变化(P-P₀)的分布。

图6 单井生产问题的有限元网格
Fig.6 Grids used for FEM in single well problem

图7 生产1个月后孔隙压力变化的分布
Fig.7 Pore pressure change distribution after one month of production

图 8给出了生产1个月后,计算区域内竖向位移的分布情况。可见随着生产井定压生产,整个计算区域出现了明显的变形,而且井单元处竖向位移显著,这反映了实际生产中的井沉降现象。孔隙度和渗透率对应力的敏感性对流固耦合问题的计算结果有着很大的影响。图9给出了考虑应力对孔隙度与渗透率的影响时的生产井产量,可见产量出现了明显的减少。

图8 生产1个月后的竖向位移分布
Fig.8 Vertical displacement after one month of production

图9 考虑应力对孔渗系数影响时的产量
Fig.9 Productions calculated using fixed and stress-dependent porosity and permeability

基于多孔介质力学、渗流力学及u -p假定,建立了渗流与岩体变形耦合问题的基本方程,并采用一种有限元-有限体积混合方法对方程进行离散和求解,此方法对平衡方程与连续性方程分别采用不同的数值离散方法和不同精细程度的网格进行计算,既能保证渗流区域的计算精度,又能够提高计算效率。

一维固结问题的数值计算结果与理论解吻合较好,验证了该方法的准确性。对单井定压生产问题的算例进行了模拟,给出了生产过程中孔隙压力与竖向变形的分布,计算结果表明,当考虑孔隙度和渗透率随应力变化时,定压生产井的产量有明显的减少。