一种含黏弹性流体孔隙介质纵波速度预测方法及其应用

1.本发明涉及地震岩石物理领域，具体涉及一种含黏弹性流体孔隙介质纵波速度预测方法及其应用。


背景技术：

2.孔隙介质的弹性波传播理论在油气勘探领域具有非常重要的理论和实际指导意义。弹性波传播特性与储层的物理性质密切相关。理解储层的弹性波波传播特性有利于更加明确的认识孔隙介质的基本特性。biot(biot,1956a,b,1962)提出了研究双相介质弹性波的动力学方程，描述了弹性波在饱和流体孔隙介质中传播的规律，并从理论上预测了慢纵波的存在。然而，biot模型假设岩石内部所有孔隙均一，且流体在波的激励下仅在波传播方向上发生管中层流，这些与实际非均质性岩石内部的复杂情况不相符的基本假设，使得该理论很难解释波的高频散和强衰减现象(dvorkin and nur,1993)。dvorkin与nur(1993)将biot流和squrit流放在一个力学模型中考虑，建立了bisq(biot/squirt)模型，有效的解释了波的高频散和强衰减现象。yang和zhang(2000,2002)将bisq模型推广至各向异性，随后学者们(yang et al.,2002；meng et al.,2003；li and tao,2009)基于bisq理论进行了波场模拟的研究，nie和yang(2008)将骨架黏弹性引入bisq模型研究弹性波传播特征。ba等(2011)考虑球形孔隙缩放振动，通过势能函数、动能函数、耗散函数重新推导了孔隙介质弹性波传播特征，形成了双重孔隙介质biot-rayleigh理论。zhang等(2021)基于固体复合材料的微分有效介质(dem)理论和双重孔隙介质的biot-rayleigh理论，建立了含包裹体的流体饱和无限孔隙介质的波传播理论。这些研究更加明确了饱和牛顿流体孔隙介质弹性波传播特性。
3.传统孔隙介质弹性波传播模型假设孔隙中饱和经典牛顿流体，然而大量理论和实验研究表明，现有牛顿流体孔隙介质波动方程已不适用于低孔、低渗致密油储层的复杂流动的情况。为了有效的利用地震技术提高剩余油采收率以及探寻非常规油藏，人们试图研究黏弹性流体对多孔介质动力学的影响。黏弹性是指在介质内部，之前时刻所发生的应变与变形，会对其后时刻的受力与变形产生影响。早在19世纪人们就意识到地球岩石中的黏弹性，1868年，maxwell在实验中首先发现材料的应力松弛现象，由此提出了黏弹性理论和经典的maxwell模型。黏弹性介质是指具有黏性和弹性两种特征得介质，其特征分布范围很广，部分介质的特征表现更接近于黏性流体介质，如maxwell黏弹性介质，部分则更接近于弹性固体介质，如kelvin黏弹性介质。黏弹性介质作为一种特殊的介质，在不同力学(地震学、材料力学、生物力学、连续介质力学)的研究中都具有重要的研究价值。
4.黏弹性流体是一种特殊的非牛顿流体，因此，针对黏弹性的研究还需要从非牛顿流体的研究中入手。非牛顿流体的特征表现在应力与应变的比值不是一个常数，具体表现可以通过流体黏度随频率的动态变化进行表征(陈文芳和蔡扶时，1983)。非牛顿流体的流动在工业中有着广泛的经验，如提高石油的采收率等(liu et al.,1998)。通过管道流对流体传输进行建模，可以实现对岩石孔隙介质中的流体流动进行有效的描述(bird et al.,
1987；liu et al.,1994)。tsiklauri等人在2001年的研究中，采用流体动力学的方法，考虑了在圆管中黏弹性流体的流动特征，并推广至多孔介质中，利用maxwell型黏弹性流体表征孔隙流体的黏弹性特征，建立了波传播模型，研究了弹性波响应特征的规律(tsiklauri and beresnev,2001)。
5.liu和kataube(1990)认为含稠油流体孔隙介质的液相在波激励下的同时存在弹性变形和黏性变形，并在理论上推导了慢横波的存在。carcione(1998)将标准线形体的本构关系引入biot模型中，使得快纵波及横波的频散和衰减显著增加。johnson等(1987)提出了动态渗透率的模型，表明波激励作用下的振动对应的流体渗流规律与频率有关。卢明辉等(2009)将biot模型和含黏性介质的流变学本构关系相结合，重新推导了波动方程，并定义了慢横波的特征频率。tsiklauri和beresnev(2001,2003)用maxwell流体替代牛顿型流体，研究流体饱和的多孔介质波传播特征。崔志文等(2010)通过引入狄波拉数(de)，实现将maxwell流体代替bisq论中的牛顿型流体，进行了喷射流机制的饱和非牛顿流体孔隙介质波传播模型研究，但该结果不影响弹性波的频散及衰减的大小。凌云等(2017)基于zener模型(zener and siegel,1945)与bisq模型，同时考虑孔隙流体的流动机制与固体骨架的黏弹性机制，分析了黏弹性孔隙介质的波传播特征。另外，一些学者通过将分数阶导数引入到渗流力学中，建立了孔隙流体介质非牛顿黏弹性流体的松弛模型，研究非牛顿型黏弹性流体在孔隙介质中的运动特性(friedrich,1991；shan et al.,2009；熊繁升,2016)。
6.常规油气探测中将孔隙流体视为牛顿流体。然而，在实际的油气勘探和开发中普遍存在非牛顿流体。稠油、沥青油、页岩油等均属于非牛顿型流体，它们是油藏开发工程领域的一个重要研究对象。沥青作为一种非常规油藏，其黏度值范围在10～1000pa
·
s(yuan et al.，2016)。页岩油作为国内重要的非常规油藏类型，成为目前国内油藏开发的热点。页岩油储集层具有低孔低渗的特点(邹才能等,2013；杨华等,2013)。黏度较高的页岩油在孔隙、裂隙中的运动特征属于非达西渗流(zeng et al.,2011)。因此，在对我国页岩油层系的致密储层的研究中，理应考虑流体的黏弹性特征。
7.目前的孔隙介质理论模型大多是基于牛顿流体的原理建立的，不能有效的将流体的黏弹性与复杂孔隙介质波传播理论相结合，因此，需要建立更准确的反应地层特征的黏弹性孔隙介质理论模型。


技术实现要素：

8.针对上述现有技术存在的问题，本发明提供了一种含黏弹性流体孔隙介质纵波速度预测方法及其应用，基于黏弹性流体的本构关系和bisq模型的应力应变关系，推导含黏弹性流体孔隙介质弹性波传播理论模型，形成了黏弹性bisq模型获取纵波速度及衰减的表达式，预测黏弹性流体孔隙介质的弹性波传播速度频散和衰减的特征，指导储层油气资源勘探开发。
9.第一方面，本发明提供了一种含黏弹性流体孔隙介质纵波速度预测方法，该方法包括如下步骤：
10.基于bisq模型应力应变方程及黏弹性流体的应力应变本构关系获取黏弹性流体bisq模型应力应变方程；
11.基于孔隙介质弹性波传播特征获取黏弹性流体bisq模型的波动方程；
12.基于所述黏弹性流体bisq模型的波动方程获取纵波的克里斯托弗尔方程；
13.基于所述纵波的克里斯托弗尔方程求解复波数；
14.基于所述复波数与波速和衰减之间的关系获取纵波速度和衰减预测值。
15.在一些实施方式中，所述bisq模型包括：
16.固相及液相微元的应力张量τ
ij
与s
ij
，分别为：
17.τ
ij
＝aeδ
ij
2ne
ij
qζδ
ij
；
18.s
ij
＝qeδ
ij
fζδ
ij
；
19.其中，τ
ij
与s
ij
分别为固相及液相微元的应力张量，e
ij
是应变张量，是应变张量，是应变张量，为求散度，u为系统中固相的位移，u为系统中液相的位移，a、n、q、f为弹性参数，具体的：
20.n＝μd；
[0021][0022]
其中，
[0023]
其中
[0024]
式中，μd为干燥岩石的剪切模量，kd为干燥岩石的体积模量，φ为孔隙度，k0为基质模量，k
fl
为流体的体积模量，j0，j1分别为0阶和1阶贝塞尔函数，r为特征喷射流长度，ρa为附加流体密度，ωc为特征频率，η为流体黏度，κ为渗透率，ρ
fl
为流体密度。
[0025]
在一些实施方式中，所述黏弹性流体bisq模型的液相的应力应变方程为：
[0026][0027]
式中，上标“·”表示对时间导数，η
*
(ω)为黏弹性流体的动态黏度，ω＝2πf为角频率，f为频率。
[0028]
在一些实施方式中，所述孔隙介质弹性波传播特征获取黏弹性流体bisq模型的波动方程，包括：
[0029][0030][0031]
式中，b＝ηφ2/κ，ρ1＝(1-φ)ρ
s-ρa，ρ2＝φρ
fl-ρa。
[0032]
在一些实施方式中，所述基于所述黏弹性流体bisq模型的波动方程获取纵波的克里斯托弗尔方程，包括：
[0033]
(1)获取固相位移及流相位移：
[0034]
u＝c1e
i(kx-ωt)
；
[0035]
u＝c2e
i(kx-ωt)
；
[0036]
(2)基于(1)中固相位移及流相位移代入所述黏弹性流体bisq模型的波动方程，获取纵波的克里斯托弗尔方程为：
[0037][0038][0039]
在一些实施方式中，所述基于所述纵波的克里斯托弗尔方程求解复波数k，包括：
[0040]
基于纵波的克里斯托弗尔方程的行列式为0求解复波数k。
[0041]
在一些实施方式中，所述基于所述复波数k与波速和衰减之间的关系获取纵波速度和衰减预测值，包括：
[0042][0043][0044]
第二方面，本发明提供了一种基于第一方面所述纵波速度预测方法在储层油气勘探开发中的应用，具体来说，基于第一方面所述纵波速度预测方法获得的纵波速度和衰减，预测含黏弹性流体孔隙介质的弹性波传播特征，指导储层油气勘探开发。
[0045]
本发明的一种含黏弹性流体孔隙介质纵波速度预测方法及其应用，具备如下有益效果：本发明基于bisq模型的应力应变方程，结合黏弹性流体的应力应变本构关系，推导含黏弹性流体的bisq模型的波动方程，预测孔裂隙介质黏弹性纵波速度和衰减。本发明的黏弹性流体bisq模型充分的考虑了黏弹性流体对孔隙介质弹性波传播特征的影响，有效的预测含黏弹流体的孔隙介质的纵波速度。
附图说明
[0046]
图1是本技术实施例中一种含黏弹性流体孔隙介质纵波速度预测方法的流程示意图；
[0047]
图2是黏弹性流体孔隙介质模型预测的纵波速度(图a)和衰减(图b)随频率变化关系图；
[0048]
图3为含黏弹性流体孔隙介质模型预测含沥青砂岩的纵波速度值与实验测量结果；
[0049]
图4为含黏弹性流体孔隙介质模型预测含沥青碳酸盐岩的纵波速度值与实验测量结果。
具体实施方式
[0050]
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅用以解释本发明,并不用于限定本发明。
[0051]
本技术实施例提供了一种含黏弹性流体孔隙介质纵波速度预测方法，基于bisq模
型的应力应变方程，结合黏弹性流体的应力应变本构关系，推导含黏弹性流体的bisq模型的波动方程，预测孔裂隙介质黏弹性纵波速度和衰减。分析含黏弹性流体孔隙介质中的弹性波传播速度频散和衰减的特征。
[0052]
本技术实施例的一种含黏弹性流体孔隙介质纵波速度预测方法，包括：
[0053]
步骤1，基于bisq模型应力应变方程及黏弹性流体的应力应变本构关系获取黏弹性流体bisq模型应力应变方程；
[0054]
在该步骤1中，具体包括如下步骤11、步骤12、步骤13：
[0055]
步骤11，获取黏弹性流体的应力应变本构关系，包括：
[0056]
黏弹性流体特征可以通过如下方程表示：
[0057][0058]
g(t-t
*
)为流体在发生弛豫时间段内的松弛模量，其表达式为：
[0059][0060]
其中，η0为初始黏度，β为应力松弛时间。
[0061]
对黏弹性流体特征方程(1)进行傅里叶变换，其表达式为：
[0062]
τ＝2η
*
(ω)d
1 (3)
[0063]
当孔隙中饱和黏弹性流体时，黏弹性流体在波激励作用下产生弛豫现象，而这种弛豫现象可以利用黏弹性流体中的maxwell本构方程进行描述，即通过引入动态黏度描述其对波传播特征的影响：
[0064][0065]
式中，η
*
(ω)为动态黏度，ω＝2πf为角频率，f为频率。
[0066]
步骤12，获取bisq模型的原有的固相及液相微元的应力张量τ
ij
与s
ij
，具体包括：
[0067]
基于bisq模型的基本假设，bisq模型固相和液相的应力应变方程为：
[0068]
τ
ij
＝aeδ
ij
2ne
ij
qζδ
ij (5)
[0069]sij
＝qeδ
ij
fζδ
ij (6)
[0070]
其中，τ
ij
与s
ij
分别为固相及液相微元的应力张量，分别为固相及液相微元的应力张量，为求散度，u为系统中固相的位移，u为系统中液相的位移，a、n、q、f为弹性参数，可以通过下式确定：
[0071]
n＝μd[0072][0073][0074][0075]
[0076][0077][0078]
式中，μd为干燥岩石的剪切模量，kd为干燥岩石的体积模量，φ为孔隙度，k0为基质模量，k
fl
为流体的体积模量，j0，j1分别为0阶和1阶贝塞尔函数，r为特征喷射流长度，ρa为附加流体密度，ωc为特征频率，η为流体黏度，κ为渗透率，ρ
fl
为流体密度。
[0079]
步骤13，在bisq模型应力应变方程的基础之上结合黏弹性流体的应力应变本构关系形成黏弹性流体bisq模型应力应变方程，包括：
[0080]
基于步骤11中的黏弹性流体本构方程的描述，对步骤12中bisq模型液相的应力应变方程进行改进，获取黏弹性流体bisq模型液相的应力应变方程：
[0081][0082]
其中，上标“·”表示对时间导数。
[0083]
当然，黏弹性流体bisq模型固相的应力应变方程仍然是上式(5)。
[0084]
步骤2，基于孔隙介质弹性波传播特征获取黏弹性流体bisq模型的波动方程；
[0085]
在该步骤2中，具体包括如下步骤21、步骤22：
[0086]
步骤21，获取岩石中固体、流体以及流固之间相互耦合的运动的方程：
[0087]
弹性波经过岩石时，流体与固体之前会产生相对运动，导致能量耗散。固体、流体以及流固之间相互耦合的运动的方程如下：
[0088][0089][0090]
其中，i,j＝1,2,3，表示空间上的三个方向，xj，为j方向上的空间坐标，ui与ui分别为i方向上的固相和流相位移，b＝ηφ2/κ，ρ1＝(1-φ)ρ
s-ρa，ρ2＝φρ
fl-ρa。
[0091]
步骤22，将式(5)与(8)代入式(9)与(10)，即得到含黏弹性流体饱和bisq模型的波动方程：
[0092][0093][0094]
步骤3，基于所述黏弹性流体bisq模型的波动方程获取纵波的克里斯托弗尔方程；
[0095]
该步骤3中具体包括：
[0096]
步骤31，在各项同性介质中，假定固相位移及流相位移表达式如下：
[0097]
u＝c1e
i(kx-ωt) (13)
[0098]
u＝c2e
i(kx-ωt) (14)
[0099]
步骤32，基于步骤31中固相位移及流相位移代入上述式(11)与(12)，获取纵波的克里斯托弗尔方程为：
[0100][0101][0102]
步骤4，基于所述纵波的克里斯托弗尔方程求解复波数k；
[0103]
为使c1与c2有非零解，式(15)及(16)的行列式必须等于零：
[0104][0105]
对式(17)的求解方法中，可以得到形式相同的一个关于复波数平方k2的方程：
[0106]
a1k4 a2k2 a3＝0 (18)
[0107]
其中，a1、a2、a3分别表达如下：
[0108]
a1＝(a 2n)(f-2η
*
(ω)iω)-q
2 (19)
[0109][0110][0111]
步骤5，基于所述复波数k与波速和衰减之间的关系获取纵波速度和衰减预测值。
[0112]
通过复波数可以求解得到纵波速度和衰减的表达式：
[0113][0114][0115]
本技术实施例中，还提供了一种基于上述纵波速度预测方法在储层油气勘探开发中的应用，具体的，该应用中，包括：
[0116]
通过数值模拟分析黏弹性bisq模型预测纵波速度与衰减的特征，在此基础之上结合超声波实验数据验证该模型的应用效果。
[0117]
基于上述一种含黏弹性流体孔隙介质纵波速度预测方法，下面以一具体实施例说明本发明的实施过程。
[0118]
在进行实施例之前，先对本发明的模型进行数值模拟，其岩石物理参数值见表1：
[0119]
表1黏弹性bisq模型的模拟参数
[0120]
[0121][0122]
本发明提出的模型与bisq模型数值模拟比较结果见图2。由图2a可知，频率较低时，两个模型计算的纵波速度一致，随着频率升高，纵波速度均增大，本发明提出的模型计算的纵波速度大于bisq模型计算结果。由图2b可知，本发明提出的模型计算的纵波衰减要大于bisq模型计算结果，且本发明的模型所获得的衰减峰与bisq模型相比会多一个黏弹性衰减峰。
[0123]
为了进一步说明本发明的效果，本实施案例是对加拿大阿尔伯塔省的沥青砂岩与沥青碳酸盐岩样本进行纵波速度的实验测量值与理论预测值的对比。
[0124]
实验测量频率为1mhz，在0～120℃范围内一共测量了11个温度值的沥青砂岩。孔隙度为38.19％，岩石干骨架密度为2680kg/m3，岩基质体积模量为32gpa，渗透率为100md，沥青密度为1080kg/m3，沥青的体积模量为2.85gpa，沥青黏度值为0.5pa
·
s。
[0125]
图3展示了沥青砂岩实际测量的纵波值与模型的预测值之间关系图，图中虚线表示实验测量纵波速度值，实线表示本发明所提出的黏弹性bisq模型预测结果。通过将模型预测的纵波速度与实验值的对比发现，本发明所提出的黏弹性bisq模型能有效的预测含黏弹流体的孔隙介质的纵波速度。
[0126]
实验测量频率为1mhz，在0～120℃范围内一共测量了12个温度值的沥青碳酸盐岩。孔隙度为15.5％，岩石干骨架密度为2700kg/m3，岩基质体积模量为64gpa，渗透率为100md，沥青密度为1080kg/m3，沥青的体积模量为2.85gpa，沥青黏度值为0.5pa
·
s。
[0127]
图4展示了沥青碳酸盐岩实际测量的纵波速度值与模型预测值之间关系图，图中虚线表示实验测量的纵波速度值，实线表示模型本发明所提出的黏弹性bisq模型预测结果。通过将模型预测的纵波速度与实验值的对比发现，本发明所提出的黏弹性bisq模型能有效的预测含黏弹流体的孔隙介质的纵波速度。
[0128]
综上所述，本发明充分的考虑了黏弹性流体对孔隙介质弹性波传播特征的影响，推导黏弹孔隙介质弹性波传播方程，形成了含黏弹性孔隙介质纵波速度及衰减的计算方法。
[0129]
以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。