基于均匀化理论的页岩基质流固耦合尺度升级

严 侠, 王东英, 樊伟鹏, 刘丕养, 樊冬艳, 姚 军

(中国石油大学(华东)石油工程学院,山东青岛 266580)

页岩气资源分布广、储量大,但其基质渗透率极低,通常需要经过水力压裂才能进行商业化开采[1-5],压裂后的页岩储层具有多尺度孔缝特征[6-7],且处于复杂地应力场和渗流场的共同作用中,流固耦合效应明显[8-10]。其中页岩基质作为气体的主要储集空间,其流固耦合模型的建立是页岩气藏宏观流固耦合数值模拟的基础。由于页岩基质通常由微尺度的有机质和无机质组成[11-13],而这两种介质的力学性质差异较大[14-15],且气体在两者中的赋存方式和流动机制均不相同[12,16],要准确刻画气体在页岩基质中的流固耦合过程,需针对两种介质建立不同的微尺度模型,但由于计算效率问题,微尺度模型无法直接用于宏观模拟。因此有必要建立一种尺度升级方法,将微尺度上有机质和无机质各自特征有效地表征到页岩气藏宏观流固耦合模拟中。近年来,国内外学者分别采用流量等效法[17-20]和均匀化理论[21-24],对页岩基质中有机质分布对气体宏观渗流规律的影响进行了研究[25-26],但均未考虑流固耦合作用。对此,笔者基于均匀化理论,建立考虑微尺度上有机质和无机质各自特征的页岩基质流固耦合尺度升级方法,通过数值算例验证尺度升级方法的正确性,并分析有机质力学性质、体积分数和分布对页岩气藏宏观流固耦合数值模拟的影响。

1 微观尺度流固耦合模型建立

扫描电子显微镜(scanning electron microscope, SEM)扫描研究发现,在微观尺度上,有机质分散嵌入在无机质中[17]。假设有机质和无机质均满足连续介质假设[11-13],将页岩基质看成是由有机质和无机质组成的非均质多孔弹性介质[14],其中无机质中仅存在游离气,气体的运移机制主要为黏性流和Knudsen扩散,有机质中同时存在吸附气和游离气,气体的运移机制主要为黏性流、Knudsen扩散和表面扩散[16]。

为了便于书写,采用以下广义的流固耦合模型进行描述[27],表达式为

(1)

(2)

(3)

其中

式中,σ和σs分别为总应力张量和有效应力张量;p为孔隙压力;α为Biot系数;I为单位张量;C为弹性张量;us为骨架位移;e(us)为应变张量;β为综合压缩系数;t为时间;δ为克罗内克符号;φ为孔隙度;mg为吸附/解吸项;ρg为气体密度;v为气体流速;ka为气体的视渗透率;μ为气体黏度;Ks为基质骨架体积模量;Mg为气体摩尔质量;Z为气体压缩因子;R为通用气体常数;T为储层温度。

根据Langmuir等温吸附[28],气体吸附/解吸项mg的表达式[6]为

(4)

式中,ρr为岩石密度;ρgstd为标准条件下的气体密度;VL和pL分别为有机质的Langmuir体积和Langmuir压力。

无机质中气体的运移机制主要为黏性流和Knudsen扩散,其视渗透率kai[12,16]为

(5)

式中，Kn为knudsen数。

有机质中气体的运移机制主要为黏性流、Knudsen扩散和表面扩散,其视渗透率kak[12,16]为

(6)

其中

式中,rh和τ分别为孔隙半径和迂曲度;λ为气体分子的平均自由程;reff为孔隙有效半径;b=-1为滑移系数;Ds为吸附气的表面扩散系数;Cmax和θ分别为有机质的最大吸附量和气体覆盖率[12,16]。

以上页岩基质流固耦合模型的非均质系数在无机质中取值为

α=αi,C=Ci,φ=φi,δ=0,ka=kai,

rh=rhi,τ=τi,reff=rhi.

在有机质中取值为

α=αk,C=Ck,φ=φk,δ=1,ka=kak,rh=rhk,

式中,下标i和k分别表示无机质和有机质;dm为甲烷分子直径。

2 基于均匀化理论的尺度升级

在均匀化理论中,物理量特征函数依赖于两种不同尺度:一种是宏观尺度,用来表征介质特征函数在大尺度上的平缓变化;另一种是微观尺度,用来表征特征函数在小尺度上的快速振荡变化[29]。考虑一个特征长度为L的页岩基质,其由大量的周期性微尺度元胞组成,微尺度元胞的特征长度为l,引入不同特征长度之间的比值ε=l/L。上述页岩基质流固耦合模型中的变量σ、us、v和p可由以下渐进展开式近似:

φ=φ(0)(x,y,t)+εφ(1)(x,y,t)+ε2φ(2)(x,y,t)+….

(7)

式中,y=x/ε,y为微观尺度坐标;x为宏观尺度坐标;φ(i)为关于坐标x、y以及时间t的周期性变量。

此外,根据链式准则,可以得到:

(8)

把方程(7)代入方程(1)～(3)中,并采用式(8)进行变换,比较方程中不同阶数ε的系数。

通过比较方程(3)中ε-1的系数以及方程(1)中ε-2的系数可得

(9)

由此可知p(0)和us(0)仅在宏观尺度x上变化,与微观尺度y无关,即

p(0)=p(0)(x,t),us(0)=us(0)(x,t).

(10)

通过比较方程(2)中ε-1的系数以及方程(3)中ε0的系数可得

(11)

分别对v(0)和p(1)进行分离变量分析:

(12)

式中,ω和χ为关于y的周期性变量。

把式(12)代入式(11),可以得到元胞辅助方程为

(13)

式中,ei为笛卡尔坐标系中的i-方向的单位矢量。

对方程(12)中的v(0)在元胞上进行体积平均,可以得到宏观渗流方程为

(14)

式中,|Ω|为整个元胞的体积。

〈a〉表示变量a(a指本文中涉及到的变量)的体积平均,kequ为等效视渗透率张量,表达式分别为

(15)

(16)

比较方程(1)中ε-1的系数可得

(17)

由此可知,us(1)可以写成:

(18)

与坐标y相关的周期性变量ξpq满足:

δkqδlp).

(19)

比较方程(1)中ε0的系数可得

(20)

将方程(20)在元胞上进行体积平均,并采用散度定理和周期性边界条件,得到

(21)

根据体积平均定理和周期性边界条件,对方程(21)进行整理,可以得到宏观应力平衡方程为

(22)

其中

Cequijkl=〈Cijkl+Cijmneymn(ξkl)〉,αequ=〈α〉.

比较方程(2)中ε0的系数可得

(23)

将方程(23)在元胞上进行体积平均,并采用散度定理、体积平均定理和周期性边界条件,可以得到宏观质量守恒方程为

(24)

其中

式中,aTOC为元胞内有机质的体积分数。

最后,采用下标hm表示下标(0),将宏观等效连续介质流固耦合模型表示为

(25)

(26)

(27)

可以看出宏观等效连续介质流固耦合模型(25)～(27)的形式与方程(1)～(3)一致,但是由于同时考虑了有机质和无机质的存在,模型的系数发生了变化。元胞辅助方程(13)和(19)很难求得解析解,本文中采用有限单元法[10, 30]进行数值求解。

3 算例分析

3.1 模型验证

图1 有机质均匀分布模型示意图及不同网格划分结果Fig.1 Schematic of organic matter uniform distribution model and results of different grids

图2 粗网格单元等效视渗透率随压力的变化Fig.2 Relationship between equivalent apparent permeability and pressure for coarse grids

为了进一步验证提出的尺度升级方法的正确性,分别采用细网格和粗网格对有机质半径分别为0.02、0.04和0.05 mm(对应的有机质体积分数分别为0.031、0.126和0.196)的情况进行模拟,对比3种情况下两种方法计算得到的压力观测点(1.0 mm, 0 mm)处的压力和位移观测点(1.0 mm, 2.0 mm)处的y方向位移。压力和位移的对比结果如图4所示,可以看出两种方法计算得到的结果基本一致,验证了提出的尺度升级方法的正确性。此外,还可以看出随着有机质体积分数的增大,观测点压力下降速度减小,y方向位移增大,这是因为有机质的视渗透率和弹性模量均小于无机质,有机质体积分数增大导致基质整体的渗透性和弹性刚度降低,同时有机质体积分数增大还能增加气体吸附量。

图4 不同有机质体积分数下两种方法计算得到的压力和位移对比Fig.4 Comparisons of pressure and displacement of two different methods for various aTOC

3.2 实例计算及影响因素

3.2.1 有机质力学性质对气体生产的影响

如图5为某一100 m×100 m的页岩气藏模型及其内部有机质分布,假设其基质的表征单元体为1 mm×1 mm,有机质体积分数为0.138。气藏模型满足平面应变假设,力学边界条件如图5所示,流动边界为封闭边界,初始地层压力为20 MPa,生产井位于气藏模型左下角,井底流压为5 MPa,其余模型参数与3.1节算例一致。

图5 页岩气藏模型及其基质表征单元体Fig.5 Schematic of shale gas reservoir and its REV element

图6 基质表征单元体等效视渗透率随压力的变化Fig.6 Relationship between equivalent apparent permeability and pressure for REV element

根据计算得到的表征单元体等效参数,采用有限单元法对图5中的页岩气藏进行宏观流固耦合数值模拟,模拟过程中采用以下动态渗透率模型考虑岩石变形对渗透率的影响[31]:

(28)

式中,下角标0表示初始状态;Δεv和Δp分别为体积应变和气体压力变化量。

图7为有机质弹性模量取不同值时计算得到的累积产气量对比,可以看出有机质弹性模量越小累积产气量越高。这是由于有机质弹性模量降低,导致页岩基质的可压性增强,弹性开采过程中,相同压降导致的基质孔隙度减小越多,因此累积产气量越高。此外,从图7中还可以看出页岩基质的应力敏感程度较弱,即岩石变形导致渗透率降低进而降低累积产气量的效果不明显,这也是有机质弹性模量降低导致累积产气量增高的另一原因。

图7 不同情况下页岩气藏累积产气量对比Fig.7 Comparisons of cumulative production for different cases

3.2.2 有机质体积分数对气体生产的影响

根据计算得到的各个表征单元体的等效参数,考虑岩石变形对渗透率的影响,对图5中的页岩气藏进行宏观流固耦合数值模拟。图10为不同有机质体积分数情况下计算得到的累积产气量对比。可以看出,有机质体积分数越高,早期的累积产气量越低,而后期的累积产气量越高。这是由于有机质体积分数增大虽然增加了储层中气体储量,但也导致了基质等效视渗透率减小,以及基质的应力敏感程度增大。由于早期累积产气量受视渗透率的影响较大,因此有机质体积分数越高早期累积产气量越低,而后期累积产气量主要受气体储量的影响,因此有机质体积分数越高后期累积产气量越高。

图8 不同有机质体积分数的基质表征单元体示意图Fig.8 Schematic of REV elements with different aTOC

图9 不同有机质体积分数的基质表征单元体等效视渗透率随压力的变化Fig.9 Relationships between equivalent apparent permeability and pressure for REV elements with different aTOC

3.2.3 有机质分布对气体生产的影响

图10 不同有机质体积分数下页岩气藏累积产气量对比Fig.10 Comparisons of cumulative production for shale gas reservoir with different aTOC

图11 有机质分布不同的基质表征单元体示意图Fig.11 Schematic of REV elements with different organic matter distribution

同样,根据计算得到的3个表征单元体的等效参数,考虑岩石变形对渗透率的影响,对图5中的页岩气藏进行宏观流固耦合数值模拟。图13为不同有机质分布情况下计算得到的累积产气量对比,结合3种分布情况下的无机质连通情况,可以看出,无机质连通性越差,累积产气量越低。

图12 不同有机质分布的基质表征单元体等效视渗透率随压力的变化Fig.12 Relationships between equivalent apparent permeability and pressure for REV elements with different organic matter distribution

图13 不同有机质分布下页岩气藏累积产气量对比Fig.13 Comparisons of cumulative production for shale gas reservoir with different organic matter distribution

4 结 论

(1)采用均匀化理论可以对包含微尺度有机质和无机质的页岩基质进行流固耦合尺度升级。

(2)有机质弹性模量越小、无机质连通性越好,累积产气量越高,而有机质体积分数越高,早期的累积产气量越低,后期的累积产气量越高。