地球物理学报  2020, Vol. 63 Issue (6): 2386-2399   PDF    
基于数字岩心动态应力应变模拟的非均匀孔隙介质波致流固相对运动刻画
朱伟1, 赵峦啸2,3, 王晨晨4, 单蕊5     
1. 长江大学地球物理与石油资源学院, 武汉 430100;
2. 同济大学海洋地质国家重点实验室, 上海 200092;
3. 同济大学海洋与地球科学学院, 上海 200092;
4. 长江大学非常规油气湖北省协同创新中心, 武汉 430100;
5. 中煤科工集团西安研究院有限公司, 西安 710077
摘要:地震波传播激发的不同尺度的流固相对运动(宏观、中观和微观)是许多沉积岩地层中地震波频散和衰减的主要原因,然而野外观测和试验测量都难以对非均匀多孔介质孔隙压力弛豫物理过程进行精细刻画.通过数字岩石物理技术,本文建立了三个典型的数字岩心分别用于表征孔隙结构、岩石骨架和斑状饱和流体引起的非均质性,利用动态应力应变模拟技术计算数字岩心的位移和孔隙流体增量图像.通过分析和比较三个数字岩心的位移和孔隙压力增量图像,细致刻画了发生于非均匀含流体多孔介质内的宏观、中观和微观尺度的流固相对运动:1)宏观尺度的波致孔隙流体流动导致波长尺度上数字岩心不同区域的孔隙压力和位移差异;2)中观尺度的流体流动发生在软层与硬层之间、气层与液层之间;3)微观尺度的流体流动发生在孔隙内部或相邻孔隙之间.数值模拟试验也证明基于数字岩心的动态应力应变模拟技术可以从微观尺度上更好的理解波致孔隙流体流动发生的物理机理,从而为建立岩石骨架、孔隙流体、孔隙结构非均质性和弹性波频散-衰减特征的映射关系奠定基础.
关键词: 数字岩心      动态应力应变模拟      波致孔隙流体流动      非均匀性     
Characterization of wave-induced pore fluid flow based on dynamic stress strain simulation on digital rocks
ZHU Wei1, ZHAO LuanXiao2,3, WANG ChenChen4, SHAN Rui5     
1. College of Geophysics and Petroleum Resources, Yangtze University, Wuhan 430100, China;
2. State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China;
3. School of Ocean and Earth Science, Tongji University, Shanghai 200092, China;
4. Hubei Cooperative Innovation Center of Unconventional Oil and Gas, Yangtze University, Wuhan 430100, China;
5. Xi'an Research Institute, CCTEG, Xi'an 710077, China
Abstract: The wave-induced relative fluid-solid motion at different scales (macroscopic, mesoscopic and microscopic) is the main cause for wave dispersion and attenuation in many sedimentary rocks. Nevertheless, the detailed fluid pressure diffusion process in heterogeneous porous rocks is difficult to be observed and characterized by laboratory and filed measurements. To tackle this problem, on the basis of digital rock physics technique, three typical digital rocks are constructed to represent the heterogeneity caused by pore structure, rock skeleton and patchy saturation, respectively. Then dynamic stress-strain simulation is applied to calculate the displacement and pore pressure increment images of the digital rocks. By analyzing and comparing these images, details concerning the relative fluid-solid motion at macro-, meso- and micro-scale in fluid saturated heterogeneous porous media are described as follows. 1) Macroscopic flow yields pore pressure and displacement difference between different parts of digital rocks at wavelength scale. 2) Mesoscale flow mainly occurs between soft and stiff layers as well as gas and liquid saturated zones. 3) Microscopic fluid flow mainly takes place within or between adjacent pores. Numerical simulation also demonstrates that the dynamic stress-strain simulation based on digital rock can be used to effectively characterize the physics of wave-induced fluid flow at different scales, thus laying a foundation for establishing a link between the heterogeneities of rock skeleton, pore fluid, pore structure and wave dispersion and attenuation signatures.
Keywords: Digital rock    Dynamic stress strain simulation    Wave-induced fluid flow    Heterogeneity    
0 引言

地下的岩石尤其是沉积岩都是一般意义上的含流体非均匀多孔介质.当地震波在非均匀孔隙介质中传播时,由于非均质体的存在会引起孔隙压力的差异,流体就会在孔隙压力高的部分和孔隙压力低的部分之间交换,这样的物理过程通常被称为波致流体的流动(WIFF)(Müller et al., 2010).需要注意的是,虽然通常称之为流体的流动,但本质上这是孔隙压力弛豫的过程(Pride et al., 2004Batzle et al., 2006Zhao et al., 2015, 2017a).WIFF被认为是引起地震波频散和衰减的主要因素之一.波的频率以及岩石骨架、孔隙结构、流体状态的非均质性都会影响WIFF,从而影响地震波的频散和衰减特征.研究WIFF及其对弹性波传播和反射特征的影响对于综合利用多尺度的地球物理数据及其利用地震波的衰减特征来揭示地下介质的非均质性和流动性能都具有重要意义(Zhao et al., 2017b; Chen et al., 2018).

根据孔隙压力弛豫过程发生的尺度,通常将波动引起的孔隙流体的流动归纳为三大类:1)宏观尺度的WIFF一般发生在波长的尺度上,这是由于波峰和波谷引起的孔隙压力差异造成的(Biot, 1956a, b1962).2)在地震波挤压多孔岩石的过程中,由于较软的孔隙(比如微裂缝、孔隙间的喉道、破碎的颗粒接触)与较硬孔隙的压缩性不同,流体就会在较软的孔隙和较硬的孔隙间流动,这种发生在微观尺度的流体局部流动的物理过程通常被称为挤喷流(Mavko and Nur, 1975O′Connell and Budiansky, 1974Dvorkin et al., 1995Chapman et al., 2002Yao et al., 2015张广智等,2017).3)相比于发生在宏观波长尺度的Biot流和微观孔隙尺度的挤喷流,波场引起的孔隙压力弛豫还可以发生在比孔隙尺度要大但比波长尺度要小的中观尺度上,我们通常称之为中观流(White,1975Gurevich et al., 1997Pride et al., 2004Carcione and Picotti, 2006Ba et al., 2008, 2017巴晶等,2012Zhao et al., 2017b郭梦秋等,2018).

目前对WIFF物理机制和响应特征的理解一方面主要来自实验室测量的频散和衰减(Murphy Ⅲ,1984Batzle et al., 2006未晛等,2015Subramaniyan et al., 2015Spencer and Shine, 2016Yin et al., 2017), 这些实验室测量数据有利于标定影响弹性波频散特征的主控物理因素,但却难以对控制WIFF的物理过程和微观尺度的波致孔隙压力弛豫进行详细刻画.另外一方面,很多学者提出各种理论模型计算各种尺度的WIFF引起的弹性波速度频散和衰减特征,但这些理论模型大多数都是基于对实际复杂非均匀孔隙介质的简单抽象,往往与实际岩石还有一定差距.一定程度上,数字岩石物理技术既能够模拟实际岩心的复杂孔隙结构特征,又能够从微观尺度上观测WIFF引起的应力、孔隙压力和位移变化,这为更好地理解WIFF的物理机理和响应特征提供了非常好的技术手段.

数字岩石物理中的数字岩心相当于物理实验中的真实岩心.建立数字岩心的方法分为两大类.1)数学重建方法,主要包括过程法(田志等,2019)、模拟退火法(莫修文等,2016)、多点地质统计法(刘磊等,2018)和马尔科夫链-蒙特卡洛法(聂昕等,2016)等.数学重建方法建立的数字岩心在一定空间尺度上与真实岩心的孔隙结构具有某些统计相似性(朱伟和单蕊,2014).2)物理重建方法,主要包括CT成像—图像分割方法(刘向君等,2014Saenger et al., 2016Kong et al., 2019)和聚焦离子束切割—电镜扫描—图像分割方法(Kelly et al., 2016Zhu et al., 2019).物理重建方法建立的数字岩心与真实岩心的孔隙结构十分接近,但许多微裂缝仍有可能无法分辨,且成像分辨率与范围相互制约.

数字岩石物理中的弹性/黏弹性模拟技术主要包括四种.1)静力学模拟(Garboczi,1998Arns et al., 2002Zhang et al., 2016),假定数字岩心具有均匀的初始应变,利用共轭梯度迭代法计算最小弹性势能对应的应力应变场,由平均应变和平均应力计算等效弹性模量.由于假设骨架和流体均为线弹性体,模拟结果可认为是动态弹性模量的低频极限.共轭梯度迭代过程不能反映真正的应力应变演化过程,故静力学模拟无法展现数字岩心内部的WIFF.2)准静态应力应变模拟(Quintal et al., 2016, 2019),模拟数字岩心在加载阶跃应变之后的应力松弛过程,计算等效频散和衰减曲线.准静态模拟假设骨架为线弹性体,流体为牛顿流体,求解准静态方程.准静态模拟不考虑惯性力,与波传播过程有差异.3)透射波模拟(Saenger and Shapiro, 2002Saenger et al., 2011Sell et al., 2016Zhu and Shan, 2016印兴耀等,2016Li et al., 2019),模拟波在数字岩心中传播,利用传播时间和距离计算等效速度.透射波模拟的波场快照能够用于观测数字岩心内部的WIFF,但入射波的频带较宽,各频率的振幅不同,可能造成WIFF分析不易.4)动态应力应变模拟(Zhang and Toksöz,2012Zhu et al., 2017Wang et al., 2017Das et al., 2019),模拟数字岩心发生周期性受迫运动时内部的波场演化,利用时间迭代过程中输出的平均应变和应力曲线计算等效速度和衰减参数.动态应力应变模拟假设骨架为线弹性体,流体为牛顿流体或麦克斯韦体,求解动力学方程.Zhang和Toksöz(2012)首次利用动态应力应变模拟计算数字岩心在100 kHz的速度和衰减参数.Zhu等(2017)建立一种沿特定方向计算波传播速度和衰减参数的动态应力应变模拟方法,并展示了模拟过程中数字岩心的应力场和应变场.Wang等(2017)计算了二维含孤立孔隙模型的杨氏模量的频散和衰减.Das等(2019)利用COMSOL软件开展动态应力应变模拟,分析了二维含裂缝模型的频散曲线和裂缝扁率、流体黏滞系数的关系,研究了应变所致的孔隙压力变化.

在一定空间尺度上数字岩心表达了近似真实的微观结构和物质组成.因此,基于数字岩心的动态应力应变模拟不仅可以计算等效速度与衰减参数,也能用于观测数字岩心内部的WIFF.本文通过数字岩石物理技术设计了3个数字岩心.第一个模拟岩石骨架包含孔隙尺度非均质性且流体单一饱和的情形;第二个模拟岩石骨架含有中观尺度非均质性(比孔隙尺度要大)的情形,数字岩心的中间含有一层比较软的岩石骨架;第三个模拟岩石中含有斑状饱和流体的情形,数字岩心中间含有一层气层.通过三个数字岩心在动态应力应变模拟过程中波场图像细节的对比分析,刻画非均匀孔隙介质中的波致流固相对运动,这有利于从微观尺度上更好的理解WIFF发生的物理机理,从而为建立岩石骨架、孔隙流体、孔隙结构非均质性和弹性波频散衰减特征的映射关系奠定基础.

1 动态应力应变模拟的基本原理 1.1 边界条件

根据Zhu等(2017)的数值模拟方案,当计算纵波的传播性质时数字岩心的边界条件设置如图 1所示.在图 1中,长方体表示数字岩心,其上、下两个表面与坐标轴Z垂直,前、后两个侧面与坐标轴X垂直,左、右两个侧面与坐标轴Y垂直.

图 1 数字岩心的边界条件设置 Fig. 1 Setup of boundary condition on digital rock

数字岩心的4个侧面上应用了周期性边界条件.周期性边界条件是分子动力学数值模拟中较常用的一种边界条件(杨华,2016丁兢娜,2017).在这里相当于将数字岩心相同方向上的两个侧面的外表面连接在一起.当波从数字岩心内部入射至某个侧面时,波将从相同方向的另一个侧面进入数字岩心内部,且传播方向不发生变化.

数字岩心的下表面静止不动,上表面受到外部作用,发生受迫运动.数字岩心上表面质点的速度分量可以用下面的公式表示.

(1)

式中,t表示时间;vxvyvz分别是质点在XYZ方向上的速度分量;sin(2πft)是频率为f的正弦函数;A是控制数字岩心应变值的比例因子,一般取大于0的数值.

在公式(1)所示的外部作用下,数字岩心沿Z轴方向(即垂直方向)发生拉伸—恢复周期性变形运动.

1.2 波动方程数值模拟

假设数字岩心的骨架为均匀各向同性的线弹性介质,其本构关系为

(2)

式中,λμ是拉梅常数, σij表示应力,ekkeij表示应变,δij是克罗尼克函数.

假设孔隙流体为牛顿流体,其本构关系为

(3)

式中,p是流体压力,ηλημ表示黏滞系数,假设ηλ=ημ, 表示应变率.

数字岩心中质点的运动方程为

(4)

式中,ui表示位移,表示位移的二阶时间导数,ρ表示密度,σij, j表示应力的一阶空间导数,gi表示体力分量,在本文中gi=0.

方程(2)、(3)和(4)构成了二阶位移-应力格式波动方程.在边界条件作用下,求解波动方程可以得到随时间变化的位移场、速度场、应变场和应力(压力)场.

在数字岩心中求解波动方程的有效方法是旋转交错网格有限差分方法(Saenger et al., 2000, 2011).具有二阶精度的旋转交错网格的设置如图 2所示,位移和速度分量设置在网格的顶点,应变、应力、弹性模量、黏滞系数和密度设置在网格中心.物理量差分的计算分为两步:第一步计算网格对角线方向的差分,第二步计算坐标轴方向的差分.假设Ψi(i=1, 2, …, 8)表示网格顶点的位移或速度分量,则它们在网格中心处的差分ΔΨx,ΔΨy和ΔΨz

图 2 旋转交错网格示意图 Fig. 2 Rotated staggered grid

(5)

式中,Δx、Δy和Δz表示网格在XYZ方向的长度.对于立方体网格而言,Δxyz.

(6)

式中,di(i=1, 2, 3, 4)表示网格的四个对角线方向,Δr是网格对角线的长度.

当计算了位移和速度的有限差分后,可以计算网格中心处的应变和应变率,最后利用公式(2)和(3)计算网格中心处的应力.计算应力差分的方法与计算位移差分的方法相同.位移和速度的差分与弹性模量和黏滞系数所在位置相同,因此无需对弹性模量和黏滞系数进行插值.这保证了有限差分法在强烈非均匀的数字岩心中的数值稳定性,无需对流固边界做特殊处理.

1.3 稳定性条件和精度要求的限制

当采用动态应力应变模拟方法计算数字岩心的等效速度和衰减参数时,波长与非均质散射体直径的比值必须满足有效介质理论成立的范围(Mavko et al., 2009).因此,数值模拟的频率需足够低.当观测数字岩心内部的WIFF时,在满足稳定性条件和精度要求(Saenger et al., 2000, 2011Zhang et al., 2010张鲁新等,2010Zhang et al., 2014)的情况下,模拟频率可以适当提高.当数字岩心中硬孔隙较多(可以认为其孔隙纵横比较大)时,一般认为数字岩心微观尺度的WIFF只有在模拟频率比较高的条件下才会比较明显,这是因为其孔隙压力弛豫的时间往往比较短.

数字岩心像素网格的大小根据数值频散和流体趋肤深度的分析结果确定.对于骨架而言,限制网格大小的因素是数值频散.对于牛顿流体而言,需要综合数值频散和黏滞趋肤深度的空间采样.趋肤深度与波的频率、流体黏滞系数和密度有关.趋肤深度大于网格长度的三倍才能保证数值模拟的精度(Saenger et al., 2011).因此,网格尺度的上限由数值频散分析和趋肤深度空间采样分析所得的最小网格尺度决定.

稳定性条件决定数值模拟时间步长的上限.弹性固体对应的稳定性条件和牛顿流体对应的稳定性条件不同.对于牛顿型流体并采用1.2节所述的迭代格式,时间步长需满足Zhang等(2010)的稳定性条件.当牛顿流体的黏滞系数较大时,该稳定性条件对时间步长的限制较强.数字岩心中同时包含固体骨架和牛顿流体,时间步长的上限须同时满足上述两个稳定性条件.

2 波致流固相对运动的数字岩石物理刻画 2.1 数字岩心

研究所用的3个数字岩心如图 3所示,它们的大小均为200×200×200像素,每个像素的边长为5 μm.数字岩心A由砂岩的CT扫描图像经过阈值分割得到,消除孤立的孔隙和骨架像素点,孔隙度约为20%,平均等效孔隙半径约20 μm,最大等效孔隙半径约40 μm.数字岩心A的骨架由一种坚硬的介质构成,孔隙饱和单一液体.在数字岩心A的中部设置一厚度为80个像素的软骨架层,孔隙流体性质不变,得到数字岩心B.在数字岩心A的中部设置一厚度为80个像素的气饱和层,骨架性质不变,得到数字岩心C.数字岩心A、B和C中的硬骨架和软骨架不特指某一具体的矿物;液体和气体也不特指某一具体的流体.对数字岩心A进行改造得到数字岩心B和C增加了可观测的数字岩心的样本数量,又不至于使结构变化过度复杂造成对比分析困难.数字岩心各组分的弹性模量、黏滞系数和密度见表 1.数值模拟的频率为1 MHz,基本满足1.3节所述的稳定性条件和精度要求.

图 3 数字岩心(0:硬骨架,2:软骨架,3:气体,5:液体) Fig. 3 Digital rocks(0: stiff matrix, 2: soft matrix, 3: gas, 5: liquid)
表 1 数字岩心参数表 Table 1 Parameter of digital rocks

在本文中,数字岩心的最大等效孔隙直径约为8.0×10-5m,频率为1 MHz的简谐波的波长约为5.0×10-3 m.故针对本文中讨论的3个数字岩心,微观尺度是小于8.0×10-5 m的尺度,中观尺度是介于8.0×10-5m和5.0×10-3 m的尺度,宏观尺度是与5.0×10-3 m相当的尺度.需要指出的是,所有关于微观、中观和宏观尺度的界定都是相对的,虽然地震波的波长一般在几十米至几百米数量级,测井声波的波长一般在几十厘米的范围内,但数字岩心模拟的不同尺度的WIFF对地震波和声波尺度下的勘探仍然具有重要的启示和借鉴意义.

2.2 波场快照的静态观测

本文提取了数字岩心第一次达到最大拉伸应变时刻的波场快照,将数字岩心的三个正交方向的切片、位移场的三个正交方向的切片和孔隙压力增量的三个正交方向的切片组合成图,如图 4所示.垂直位移分量是指与坐标轴Z平行的分量,即与拉伸方向平行.水平位移分量是指与坐标轴XY平行的分量,即与拉伸方向垂直.位移的绝对值越大表示质点移动距离越大.孔隙压力增量是指数字岩心发生应变时流体压力的变化值,孔隙压力增量为负值表示流体压力相对初始值变小,孔隙压力增量为正值表示流体压力相对初始值增大.发生拉伸应变时,孔隙压力增量应为负值,孔隙压力增量的绝对值越大说明压力变化量越大,即变化幅度越大.

图 4 数字岩心和波场快照的组合排列图 Fig. 4 Snapshot arrangement of digital rock and wavefields

图 5是数字岩心A的模拟结果.总体上,垂直位移分量的空间变化不剧烈,骨架和相邻孔隙流体间的位移差异不明显,但也可以在局部区域内观测到较明显的位移差异.水平位移分量在骨架中的变化较小,在孔隙流体内的变化相对较大.流体与相邻骨架间的水平位移分量差异在某些区域十分明显.在水平切片中,孔隙压力增量存在空间变化,但没有明显的规律性.在垂直切片中,孔隙压力增量的数值(注意孔隙压力增量为负数,数值小但绝对值大,数值大绝对值小)表现为上部大、下部小的特征,即存在孔隙压力梯度.在一些狭窄孔隙内,存在孔隙压力增量的数值特别小的现象.

图 5 数字岩心A的模拟结果 各图的含义见图 4说明; 水平位移切片中的椭圆a、b和c用于讨论. Fig. 5 Simulation result of digital rock A Meaning of images is declared in Fig. 4. Ellipses named a, b and c are used for discussion.

图 6是数字岩心B的模拟结果,位移场和孔隙压力增量图像中黑色虚线框圈出的区域为软层的位置.与图 5相比,在软、硬层分界面附近,垂直位移分量在垂直方向的跃变较明显;而在软、硬层内部,垂直位移分量的变化不剧烈.在软、硬层分界面附近,由于骨架的垂直位移分量变化显著,骨架和相邻孔隙流体间的垂直位移差异增强.水平位移分量的变化特征与图 5相似.孔隙压力增量的空间结构与图 5明显不同,具有明显的三层性.软层的孔隙压力增量的数值小于硬层的孔隙压力增量的数值,顶部硬层的孔隙压力增量的数值大于底部硬层的孔隙压力增量的数值.孔隙压力增量垂向切片没有以软层中心为基点形成上、下对称的结构.软、硬层之间的孔隙压力增量的梯度较大,但并非突变.在软、硬层内部孔隙压力增量不但具有垂向变化,也具有横向变化.同样,在一些狭窄孔隙内存在孔隙压力增量的数值特别小的现象.

图 6 数字岩心B的模拟结果 各图的含义见图 4说明; 波场快照切片中的矩形线框表示软层的位置; 水平位移切片中的椭圆a、b和c用于讨论. Fig. 6 Simulation result of digital rock B Meaning of images is declared in Fig. 4. Rectangular wire frames in wave field snapshots represent the location of the soft layer. Ellipses named a, b and c are used for discussion.

图 7是数字岩心C的模拟结果,位移场和孔隙压力增量图像中黑色虚线框圈出的区域为气层的位置.在垂直切片中,骨架的垂直位移场与图 5中相应的垂直位移场相似;而在气层的上、下边界附近,骨架和邻近孔隙流体之间的垂直位移分量的差异非常明显.在上边界附近孔隙流体的垂直位移分量大于骨架;在下边界附近孔隙流体的垂直位移分量小于骨架,这表明液体和气体的分界面相对骨架发生运动.与图 5图 6相比,骨架的水平位移分量特征没有明显变化,但流体的水平位移分量的幅度明显减弱.在垂直切片中,孔隙压力增量也具有明显的三层性,中部气层的孔隙压力增量的数值最大,顶部液层的孔隙压力增量的数值居中,底部液层的孔隙压力增量的数值最小.孔隙压力增量没有以气层中心为基点形成上、下对称的结构.气、液层之间的孔隙压力增量的梯度较大.孔隙压力增量不但具有垂向变化,也具有横向变化.由于气层孔隙压力增量的绝对值以及空间变化均较小,故其内部细节不能体现.同样,在一些狭窄孔隙内存在孔隙压力增量的数值特别小的现象.

图 7 数字岩心C的模拟结果 各图的含义见图 4说明.波场快照切片中的矩形线框表示气层的位置.水平位移切片中的椭圆a、b和c用于讨论. Fig. 7 Simulation result of digital rock C Meaning of images is declared in Figure 4. Rectangular wire frames in wave field snapshots represent the location of the gas layer. Ellipses named a, b and c are used for discussion.
2.3 波场快照动态观测

本文提取了数字岩心C在第一个拉伸应变周期内20个等间隔时刻的垂直位移分量和孔隙压力增量的垂直切片.图 8为垂直位移分量的垂直切片,图 9为孔隙压力增量的垂直切片.图 8u表示图 8a图 8t所示的20张垂直位移分量快照所在的相对时间以及对应的拉伸应变的相对幅度.图 9u表示图 9a图 9t所示的20张孔隙压力增量快照所在的相对时间以及对应的拉伸应变的相对幅度.图 5图 7所示快照的时间对应图 8图 9中的“j”时刻.

图 8 数字岩心C在不同变形阶段的垂直位移分量图(a—t:位移;u:快照时间) Fig. 8 Snapshots of vertical displacement of digital rock C at different deforming stage(a—t: displacement; u: snapshot time)
图 9 数字岩心C在不同变形阶段的孔隙压力增量图(a—t:孔隙压力增量;u:快照时间) Fig. 9 Snapshots of pore pressure increment of digital rock C at different deforming stage(a—t: pore pressure; u: snapshot time)

图 8展示了数字岩心从发生拉伸应变,然后达到最大拉伸应变状态,最后恢复至原始长度的过程中,其内部垂直位移分量的变化过程.图 9展示了在这一过程中数字岩心内部孔隙压力增量的变化过程.随着数字岩心被拉伸,气层上、下边界附近的骨架和孔隙流体的垂直位移分量的差异先增大,然后逐渐减弱;孔隙压力增量的数值先逐渐减小,然后逐渐增大.在一个拉伸-恢复过程中,在相同的宏观应变状态下,数字岩心内部的垂直位移分量和孔隙压力增量分布是不对称的.通过这样一个周期内对垂直位移和孔隙压力变化的演化过程的观测,我们可以更好的刻画非均匀孔隙介质中孔隙压力弛豫的物理过程及其对位移和整体应变的影响,从而进一步有助于我们解剖弹性波频散和衰减特征背后的详细物理机制.

3 讨论

数字岩心在外力作用下发生垂向拉伸应变.孔隙压力增量的空间变化说明在孔隙流体中形成了压力梯度.孔隙压力梯度是流体运动的驱动力,是WIFF存在的间接证据.相邻的骨架与孔隙流体之间存在位移差异则直接表明数字岩心中存在WIFF.

在数字岩心的水平位移场快照中观测到了发生在孔隙内部或相邻孔隙之间的WIFF,往往表现为正负位移场相伴出现(如图 5图 6图 7的椭圆a和b内的孔隙),或者在狭窄孔隙的一端出现正或负位移场(如图 5图 6图 7的椭圆c内的孔隙),反映了孔隙拉伸过程中形成的挤喷流.挤喷流与局部孔隙结构有关,挤喷流在有些孔隙中比较明显(如图 5图 6图 7的椭圆a内的孔隙),而在另一些孔隙中相对较弱(图 5图 6图 7的椭圆b内的孔隙).挤喷流的结构和强度受中观尺度的物性变化的影响.在图 5图 6图 7的椭圆a内的孔隙中,当数字岩心含有气层时位移的幅度减弱.在图 5图 6图 7的椭圆b内的孔隙中,当数字岩心含有软层时位移的幅度减弱;当数字岩心含有气层时位移场的结构发生变化.

在含有软层(图 6)和气层(图 7)的数字岩心中观测到了中观尺度的WIFF.软层的孔隙压力变化幅度大于硬层,软、硬层之间的孔隙压力梯度较大,软层和硬层内部的压力梯度相对较小.中观尺度的WIFF发生在软、硬层的边界附近.气层的孔隙压力变化幅度小于液层,气、液层之间的孔隙压力梯度较大,气层和液层内部的压力梯度相对较小.中观尺度的WIFF发生在气、液层的边界附近.液体饱和软层与上、下层之间的压力梯度,气体饱和硬层与上、下层之间的压力梯度的方向是相反的.当数字岩心存在软层时,流体具有从硬层流向软层的趋势.当数字岩心存在气层时,流体具有从气层流向液层的趋势.液体饱和软层的孔隙压力变化幅度大于气体饱和的硬层,骨架和流体的压缩性共同决定了该变化特征,可以用Skempton系数B=dP/dσ=[1+Kϕ(Kfl-1-Ko-1)]-1(Mavko et al., 2009)说明它们对孔隙压力的影响.相比与液体饱和硬层,液体饱和软层的Kϕ减小,Kfl不变,Ko减小,故B变大.相比与液体饱和硬层,气体饱和硬层的Kϕ不变,Kfl减小,Ko不变,故B变小.因此,软层的孔隙压力变化幅度大,气层的孔隙压力变化幅度小.

图 5图 6图 7中数字岩心上部区域内的孔隙压力变化幅度大于下部区域内的孔隙压力变化幅度,上、下区域之间存在孔隙压力梯度,流体具有从数字岩心上部流向下部的趋势.这种现象可能是数字岩心上、下区域之间的应变不同造成的.如果从波在岩石中传播的角度来说,就是岩石中的某个区域与波形的相对位置关系.在与波峰或波谷对应的区域内岩石的应变较大,而其他区域内岩石的应变较小.频率为1 MHz的简谐波的波长约为数字岩心长度5倍,波长并非远远大于数字岩心的长度.由此,可以认为这种尺度的WIFF属于宏观尺度.

4 结论

非均匀孔隙介质中由于孔隙结构、岩石骨架和部分饱和流体引起的波致WIFF非常复杂,本文采用数字岩心动态应力应变模拟方法计算非均匀含流体多孔岩石发生周期性拉伸应变的位移场和孔隙压力增量场.通过比较三个表征不同类别非均匀孔隙介质的数字岩心的位移场和孔隙压力增量场的变化,本文证明通过数字岩心技术可以有效刻画并一定程度上解耦微观、中观和宏观尺度的波致孔隙压力弛豫物理过程:微观尺度的WIFF发生在孔隙内部或相邻孔隙之间;中观尺度的WIFF发生在软层与硬层之间、气层与液层之间;宏观尺度的WIFF发生在数字岩心的宏观区域之间.基于数字岩心动态应力应变模拟的波致孔隙流体流动刻画技术有利于提高非均匀孔隙介质发生在各种尺度上的流固相对运动物理机理的认识,从而为将来定量表征岩石孔隙结构、岩石骨架和部分饱和流体的非均质性与地震波频散和衰减特征提供支撑.

References
Arns C H, Knackstedt M A, Pinczewski W V, et al. 2002. Computation of linear elastic properties from microtomographic images:Methodology and agreement between theory and experiment. Geophysics, 67(5): 1396-1405. DOI:10.1190/1.1512785
Ba J, Nie J X, Cao H, et al. 2008. Mesoscopic fluid flow simulation in double-porosity rocks. Geophysical Research Letters, 35(4): L04303. DOI:10.1029/2007GL032429
Ba J, Carcione J M, Cao H, et al. 2012. Velocity dispersion and attenuation of P waves in partially-saturated rocks:Wave propagation equations in double-porosity medium. Chinese Journal of Geophysics (in Chinese), 55(1): 219-231. DOI:10.6038/j.issn.0001-5733.2012.01.021
Ba J, Xu W H, Fu L Y, et al. 2017. Rock anelasticity due to patchy saturation and fabric heterogeneity:A double double-porosity model of wave propagation. Journal of Geophysical Research:Solid Earth, 122(3): 1949-1976. DOI:10.1002/2016JB013882
Batzle M L, Han D H, Hofmann R. 2006. Fluid mobility and frequency-dependent seismic velocity-Direct measurements. Geophysics, 71(1): N1-N9. DOI:10.1190/1.2159053
Biot M A. 1956a. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅰ. Low-frequency range. The Journal of the Acoustical Society of America, 28(2): 168-178. DOI:10.1121/1.1908239
Biot M A. 1956b. Theory of propagation of elastic waves in a fluid-saturated porous solid. Ⅱ. Higher frequency range. The Journal of the Acoustical Society of America, 28(2): 179-191. DOI:10.1121/1.1908241
Biot M A. 1962. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 33(4): 1482-1498. DOI:10.1063/1.1728759
Carcione J M, Picotti S. 2006. P-wave seismic attenuation by slow-wave diffusion:effects of inhomogeneous rock properties. Geophysics, 71(3): O1-O8. DOI:10.1190/1.2194512
Chapman M, Zatsepin S V, Crampin S. 2002. Derivation of a microstructural poroelastic model. Geophysical Journal International, 151(2): 427-451. DOI:10.1046/j.1365-246X.2002.01769.x
Chen H Z, Innanen K A, Chen T S. 2018. Estimating P-and S-wave inverse quality factors from observed seismic data using an attenuative elastic impedance. Geophysics, 83(2): R173-R187. DOI:10.1190/geo2017-0183.1
Das V, Mukerji T, Mavko G. 2019. Numerical simulation of coupled fluid-solid interaction at the pore scale:A digital rock-physics technology. Geophysics, 84(4): WA71-WA81. DOI:10.1190/geo2018-0488.1
Ding J N. 2017. All-atom molecular dynamics studies of the protein-protein interactions of Ebolavirus[Ph.D.thesis] (in Chinese). Hefei: University of Science and Technology of China.
Dvorkin J, Mavko G, Nur A. 1995. Squirt flow in fully saturated rocks. Geophysics, 60(1): 97-107. DOI:10.1190/1.1443767
Garboczi E J. 1998. Finite element and finite difference programs for computing the linear electric and elastic properties of digital images of random materials. Gaithersburg, MD: National Institute of Standards and Technology.
Guo M Q, Ba J, Ma R P, et al. 2018. P-wave velocity dispersion and attenuation in fluid-saturated tight sandstones:Characteristics analysis based on a double double-porosity structure model description. Chinese Journal of Geophysics (in Chinese), 61(3): 1053-1068.
Gurevich B, Zyrianov V B, Lopatnikov S L. 1997. Seismic attenuation in finely layered porous rocks:Effects of fluid flow and scattering. Geophysics, 62(1): 319-324. DOI:10.1190/1.1444133
Kelly S, El-Sobky H, Torres-Verdín C, et al. 2016. Assessing the utility of FIB-SEM images for shale digital rock physics. Advances in Water Resources, 95: 302-316. DOI:10.1016/j.advwatres.2015.06.010
Kong L Y, Ostadhassan M, Hou X D, et al. 2019. Microstructure characteristics and fractal analysis of 3D-printed sandstone using micro-CT and SEM-EDS. Journal of Petroleum Science and Engineering, 175: 1039-1048. DOI:10.1016/j.petrol.2019.01.050
Li T Y, Wang R H, Wang Z Z. 2019. A method of rough pore surface model and application in elastic wave propagation. Applied Acoustics, 143: 100-111. DOI:10.1016/j.apacoust.2018.08.031
Liu L, Yao J, Sun H, et al. 2018. Reconstruction of digital rock considering micro-fracture based on multi-point statistics. Chinese Science Bulletin (in Chinese), 63(30): 3146-3157. DOI:10.1360/N972018-00221
Liu X J, Zhu H L, Liang L X. 2014. Digital rock physics of sandstone based on micro-CT technology. Chinese Journal of Geophysics (in Chinese), 57(4): 1133-1140.
Müller T M, Gurevich B, Lebedev M. 2010. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks-A review. Geophysics, 75(5): 75A147-75A164. DOI:10.1190/1.3463417
Mavko G, Nur A. 1975. Melt squirt in the asthenosphere. Journal of Geophysical Research, 80(11): 1444-1448. DOI:10.1029/jb080i011p01444
Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook:Tools for Seismic Analysis of Porous Media. New York, NY: Cambridge University Press.
Mo X W, Zhang Q, Lu J A. 2016. A complement optimization scheme to establish the digital core model based on the simulated annealing method. Chinese Journal of Geophysics (in Chinese), 59(5): 1831-1838.
Murphy Ⅲ W F. 1984. Acoustic measures of partial gas saturation in tight sandstones. Journal of Geophysical Research:Solid Earth, 89(B13): 11549-11559. DOI:10.1029/JB089iB13p1154
Nie X, Zou C C, Meng X H, et al. 2016. 3D digital core modeling of shale gas reservoir rocks:A case study of conductivity model. Natural Gas Geoscience (in Chinese), 27(4): 706-715.
O'Connell R J, Budiansky B. 1974. Seismic velocities in dry and saturated cracked solids. Journal of Geophysical Research, 79(35): 5412-5426. DOI:10.1029/JB079i035p05412
Pride S R, Berryman J G, Harris J M. 2004. Seismic attenuation due to wave-induced flow. Journal of Geophysical Research:Solid Earth, 109(B1): B01201. DOI:10.1029/2003JB002639
Quintal B, Rubino J G, Caspari E, et al. 2016. A simple hydromechanical approach for simulating squirt-type flow. Geophysics, 81(4): D335-D344. DOI:10.1190/GEO2015-0383.1
Quintal B, Caspari E, Holliger K, et al. 2019. Numerically quantifying energy loss caused by squirt flow. Geophysical Prospecting, 67(8): 2196-2212. DOI:10.1111/1365-2478.12832
Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2
Saenger E H, Shapiro S A. 2002. Effective velocities in fractured media:a numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting, 50(2): 183-194. DOI:10.1046/j.1365-2478.2002.00309.x
Saenger E H, Enzmann F, Keehm Y, et al. 2011. Digital rock physics:effect of fluid viscosity on effective elastic properties. Journal of Applied Geophysics, 74(4): 236-241. DOI:10.1016/j.jappgeo.2011.06.001
Saenger E H, Lebedev M, Uribe D, et al. 2016. Analysis of high-resolution X-ray computed tomography images of Bentheim sandstone under elevated confining pressures. Geophysical Prospecting, 64(4): 848-859. DOI:10.1111/1365-2478.12400
Sell K, Saenger E H, Falenty A, et al. 2016. On the path to the digital rock physics of gas hydrate-bearing sediments-processing of in situ synchrotron-tomography data. Solid Earth, 7(4): 1243-1258. DOI:10.5194/se-7-1243-2016
Spencer J W Jr, Shine J. 2016. Seismic wave attenuation and modulus dispersion in sandstones. Geophysics, 81(3): D211-D231. DOI:10.1190/geo2015-0342.1
Subramaniyan S, Quintal B, Madonna C, et al. 2015. Laboratory-based seismic attenuation in Fontainebleau sandstone:Evidence of squirt flow. Journal of Geophysical Research:Solid Earth, 120(11): 7526-7535. DOI:10.1002/2015JB012290
Tian Z, Xiao L Z, Liao G Z, et al. 2019. Study on digital rock reconstruction method based on sedimentological process. Chinese Journal of Geophysics (in Chinese), 62(1): 248-259.
Wang Z Z, Schmitt D R, Wang R H. 2017. Modeling of viscoelastic properties of nonpermeable porous rocks saturated with highly viscous fluid at seismic frequencies at the core scale. Journal of Geophysical Research:Solid Earth, 122(8): 6067-6086. DOI:10.1002/2017JB013979
Wei X, Wang S X, Zhao J G, et al. 2015. Laboratory study of velocity dispersion of the seismic wave in fluid-saturated sandstones. Chinese Journal of Geophysics (in Chinese), 58(9): 3380-3388.
White J E. 1975. Computed seismic speeds and attenuation in rocks with partial gas saturation. Geophysics, 40(2): 224-232. DOI:10.1190/1.1440520
Yang H. 2016. A molecular dynamics study on thermomechanical behaviors of active network polymers[Ph.D.thesis] (in Chinese). Beijing: Beijing JiaoTong University.
Yao Q L, Han D H, Yan F Y, et al. 2015. Modeling attenuation and dispersion in porous heterogeneous rocks with dynamic fluid modulus. Geophysics, 80(3): D183-D194. DOI:10.1190/geo2013-0410.1
Yin H J, Zhao J G, Tang G Y, et al. 2017. Pressure and fluid effect on frequency-dependent elastic moduli in fully saturated tight sandstone. Journal of Geophysical Research:Solid Earth, 122(11): 8925-8942. DOI:10.1002/2017JB01424
Yin X Y, Qin Q P, Zong Z Y. 2016. Simulation of elastic parameters based on the finite difference method in digital rock physics. Chinese Journal of Geophysics (in Chinese), 59(10): 3883-3890.
Zhang G Z, He F, Zhang J J, et al. 2017. Velocity dispersion and attenuation at microscopic and mesoscopic wave-induced fluid flow. Oil Geophysical Prospecting (in Chinese), 52(4): 743-751. DOI:10.13810/j.cnki.issn.1000-7210.2017.04.012
Zhang L X, Fu L Y, Pei Z L. 2010. Finite difference modeling of Biot's poroelastic equations with unsplit convolutional PML and rotated staggered grid. Chinese Journal of Geophysics (in Chinese), 53(10): 2470-2483. DOI:10.3969/j.issn.0001-5733.2010.10.021
Zhang W H, Fu L Y, Zhang Y, et al. 2016. Computation of elastic properties of 3D digital cores from the Longmaxi shale. Applied Geophysics, 13(2): 364-374. DOI:10.1007/s11770-016-0542-4
Zhang Y, Song LM, Deffenbaugh M, et al. 2010. A finite difference method for a coupled model of wave propagation in poroelastic materials. The Journal of the Acoustical Society of America, 127(5): 2847-2855. DOI:10.1121/1.3372640
Zhang Y, Toksöz M N. 2012. Computation of dynamic seismic responses to viscous fluid of digitized three-dimensional Berea sandstones with a coupled finite-difference method. The Journal of the Acoustical Society of America, 132(2): 630-640. DOI:10.1121/1.4733545
Zhang Y, Fu L Y, Zhang L X, et al. 2014. Finite difference modeling of ultrasonic propagation (coda waves) in digital porous cores with un-split convolutional PML and rotated staggered grid. Journal of Applied Geophysics, 104: 75-89. DOI:10.1016/j.jappgeo.2014.02.012
Zhao L X, Han D H, Yao Q L, et al. 2015. Seismic reflection dispersion due to wave-induced fluid flow in heterogeneous reservoir rocks. Geophysics, 80(3): D221-D235. DOI:10.1190/geo2014-0307.1
Zhao L X, Yao Q L, Han D H, et al. 2017a. Frequency-and angle-dependent poroelastic seismic analysis for highly attenuating reservoirs. Geophysical Prospecting, 65(6): 1630-1648. DOI:10.1111/1365-2478.12492
Zhao L X, Yuan H M, Yang J K, et al. 2017b. Mobility effect on poroelastic seismic signatures in partially saturated rocks with applications in time-lapse monitoring of a heavy oil reservoir. Journal of Geophysical Research:Solid Earth, 122(11): 8872-8891. DOI:10.1002/2017JB014303
Zhu L Q, Zhang C, Zhang C M, et al. 2019. Challenges and prospects of digital core-reconstruction research. Geofluids, 2019: 7814180. DOI:10.1155/2019/7814180
Zhu W, Shan R. 2014. Progress of digital rock physics. Oil Geophysical Prospecting (in Chinese), 49(6): 1138-1146.
Zhu W, Shan R. 2016. Digital core based transmitted ultrasonic wave simulation and velocity accuracy analysis. Applied Geophysics, 13(2): 375-381. DOI:10.1007/s11770-016-0550-4
Zhu W, Zhao L X, Shan R. 2017. Modeling effective elastic properties of digital rocks using a new dynamic stress strain simulation method. Geophysics, 82(6): MR163-MR174. DOI:10.1190/geo2016-0556.1
巴晶, Carcione J M, 曹宏, 等. 2012. 非饱和岩石中的纵波频散与衰减:双重孔隙介质波传播方程. 地球物理学报, 55(1): 219-231. DOI:10.6038/j.issn.0001-5733.2012.01.021
丁兢娜. 2017.基于分子动力学模拟方法对埃博拉病毒蛋白间相互作用的研究[博士论文].合肥: 中国科学技术大学.
郭梦秋, 巴晶, 马汝鹏, 等. 2018. 含流体致密砂岩的纵波频散及衰减:基于双重双重孔隙结构模型描述的特征分析. 地球物理学报, 61(3): 1053-1068. DOI:10.6038/cjg2018L0678
刘磊, 姚军, 孙海, 等. 2018. 考虑微裂缝的数字岩心多点统计学构建方法. 科学通报, 63(30): 3146-3157. DOI:10.1360/N972018-00221
刘向君, 朱洪林, 梁利喜. 2014. 基于微CT技术的砂岩数字岩石物理实验. 地球物理学报, 57(4): 1133-1140. DOI:10.6038/cjg20140411
莫修文, 张强, 陆敬安. 2016. 模拟退火法建立数字岩心的一种补充优化方案. 地球物理学报, 59(5): 1831-1838. DOI:10.6038/cjg20160526
聂昕, 邹长春, 孟小红, 等. 2016. 页岩气储层岩石三维数字岩心建模-以导电性模型为例. 天然气地球科学, 27(4): 706-715. DOI:10.11764/j.issn.1672-1926.2016.04.0706
田志, 肖立志, 廖广志, 等. 2019. 基于沉积过程的数字岩石建模方法研究. 地球物理学报, 62(1): 248-259. DOI:10.6038/cjg2019L0457
未晛, 王尚旭, 赵建国, 等. 2015. 含流体砂岩地震波频散实验研究. 地球物理学报, 58(9): 3380-3388. DOI:10.6038/cjg20150930
杨华. 2016.活性聚合物热力学行为的分子动力学研究[博士论文].北京: 北京交通大学.
印兴耀, 秦秋萍, 宗兆云. 2016. 数字岩石物理中弹性参数的有限差分计算方法. 地球物理学报, 59(10): 3883-3890. DOI:10.6038/cjg20161031
张广智, 何锋, 张佳佳, 等. 2017. 微观与介观波致流下的速度频散与衰减. 石油地球物理勘探, 52(4): 743-751. DOI:10.13810/j.cnki.issn.1000-7210.2017.04.012
张鲁新, 符力耘, 裴正林. 2010. 不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用. 地球物理学报, 53(10): 2470-2483. DOI:10.3969/j.issn.0001-5733.2010.10.021
朱伟, 单蕊. 2014. 虚拟岩石物理研究进展. 石油地球物理勘探, 49(6): 1138-1146.