2017, Vol. 32
2. 大庆钻探工程测井公司, 大庆 163412
2. Wireline Logging Company, Daqing Drilling & Exploration Corporation, Daqing 163412, China
随着人们对岩石储层认识的进一步加深，复杂储层的研究给传统的岩石物理实验带来了巨大的挑战，如低孔低渗致密砂岩储层、天然气水合物储层、裂隙发育的碳酸盐岩储层研究等.而以数值模拟技术为代表的新型的岩石物理研究方法，具有成本低、效率高、应用方便快捷等优点，并可弥补传统岩石物理实验的不足，受到国内外专家的广泛关注.
岩石物理的数值模拟技术主要包括：数字岩心模型的建立和岩石物理性质的数值模拟两部分.目前，比较典型的数字岩心重构方法是随机法和过程法.其中，随机法包括：高斯场法(Fatt, 1956)、模拟退火法(Hazlett, 1997；莫修文等，2016)、马尔科夫链蒙特卡洛法(Wu et al., 2004)、多点统计法(Okabe and Blunt, 2005)等.由随机法构建的三维数字岩心模型与实际的岩心切片具有相同的统计特性，但在三维孔隙空间结构上与真实岩心差异大，连通性差，导致所计算的地球物理性质与实际实验得到的数据不吻合.与随机法不同，过程法(Bryant and Blunt, 1992)是结合实际岩石的颗粒粒径分布，通过模拟沉积岩的沉积、压实、胶结过程来构建数字岩心.由于随机法无法真实展现实际岩心孔隙空间的传导特性，而过程法能够克服这一缺点，并具有高效、应用便捷等优点，所以本文选用过程法作为数字岩心的建模方法.
随着计算机技术的发展，数值模拟方法在研究复杂储层地球物理性质上具有巨大优势，并在近些年先后发展了随机行走算法(Brownstein and Tarr, 1979)、格子玻尔兹曼方法(Martys and Garboczi, 1992)、有限元法(Arns and Knackstedt, 2002)、有限差分法(马龙，2014)等.其中，随机游走法计算原理简单，比较容易实现，但由于没有考虑介质的各向异性，计算结果与实际岩心仍有一定的差距.格子玻尔兹曼方法的优点是易于处理复杂边界，并具有并行计算能力，但随着计算分量的增加，数值模拟耗时较长.有限差分法较为直观、简单，主要适用于规则的结构化网格，有限元法在处理各向异性介质和不规则的区域边界时则更为方便、有效，与有限差分法相比有限元法模拟结果更为精确.本文在前人研究的基础上，以过程法建模为基础，使用有限元法实现数字岩心的导电性模拟分析，探讨阿尔奇公式在微观领域的适用性，实现三维数字岩心模型的导电性各向异性分析，从而为新技术，新方法的开创奠定基础.
1 过程法建立三维数字岩心模型过程法是将二维岩心图像所获取岩石粒度分布、孔隙度等信息作为约束条件，通过模拟岩石沉积、压实和成岩过程来实现岩心模型重构.现将在数字岩心建模的过程中作如下假设：(1) 将沉积颗粒假设为球体；(2) 在沉积过程中只受重力的作用；(3) 颗粒之间不发生弹性碰撞；(4) 已经处于稳定状态的颗粒，不受后续沉降颗粒的影响.
1.1 沉积过程的数值模拟沉积算法模拟岩石沉积的过程要遵循重力势能最小原理.沉积步骤如下：(1) 确定数字岩心沉积区域大小及沉积厚度；(2) 在实际岩心粒径累积分布曲线上随机选取一个数值作为下落球的半径；(3) 沉积球在重力的作用下搜索下落过程中的最高落点球，如果没有找到，则球下落到底面；(4) 把沉积球在下落的过程中首次接触的落点球称为第一落点球，以第一落点球球心为原点建立球坐标系，按照重力势能最小原理，下落球沿第一落点球球面滚动，若滚动的过程中遇到其他沉积球，则称为第二落点球，此时以第一和第二落点球的连线为轴，并在重力的作用下继续转动，直到遇到第三落点球达到稳定状态；(5) 重复2~4过程，直到所有沉积颗粒均达到稳定状态，沉积过程模拟结束.
1.2 压实过程的数值模拟压实过程就是上覆岩层在重力作用下对颗粒挤压造成岩石总体积和孔隙度减小的过程.压实算法可设计为：在上覆地层压力的作用下，岩石颗粒通过向下移动一段距离来实现压实计算，并假定在沉积颗粒下移的过程中只改变颗粒的垂直坐标而并不改变成沉积颗粒的形状(Øren and Bakke, 2003)，公式为
|
(1) |
其中：z是新的z轴坐标，z0是原来的旧坐标，λ为压实因子，ξ模拟颗粒排列的变量.
1.3 胶结过程的数值模拟胶结过程是指未固结的疏松沉积物变成坚硬岩石的过程，这一过程常与压实过程同时进行，并伴随着胶结物的生长和含矿物质的溶解作用.目前在现有的技术条件下，仅仅可以模拟一些简单的成岩过程，如石英胶结物的生长和溶蚀作用可采用的模拟公式为(Øren and Bakke, 2002)：
|
(2) |
其中：R0(r)为沉积颗粒的初始半径；R(r)为从颗粒中心沿r方向的新半径；l(r)为沿沉积颗粒半径r方向，由颗粒表面到其沃龙诺伊多面晶胞表面的距离；ɑ为胶结系数；γ为反映胶结物在孔隙或喉道中的发育方向的常数.
1.4 建模实例及效果分析本文以实际的砂岩岩心为例，统计实际砂岩的颗粒粒径分布，并选取300个不同粒径的颗粒进行模拟计算.在模拟实际岩石形成的过程中，设定沉积区域的边长为2000 um，沉积颗粒个数为2500，在沉积模拟时，沉积颗粒的粒径是从粒度分布曲线上随机选取，如图 1所示，实际岩心的颗粒半径在100~200 um时颗粒数较多，模拟计算时被选中的概率也越大，而其余的颗粒编号被选中的概率较小，这符合岩石的实际沉积规律.压实过程模拟中压实系数设为0.05，模拟颗粒排列的变量取[-0.04, 0.04]之间的随机数，胶结过程的数值模拟只考虑石英胶结物的次生加大作用，可通过将沉积颗粒的半径增大1.1倍来实现.为使所构建的数字岩心模型与实际岩心的几何特征一致，可通过不断调节压实系数和胶结参数使模型孔隙度与实际岩石孔隙度误差最小，此时的模型即为最终所构建的数字岩心.
|
图 1 颗粒半径随机选择示意图 Figure 1 Schematic diagram of random selection of particle radius |
计算所得到的三维数字岩心模型如图 2所示，其中绿色部分表示骨架, 蓝色部分表示孔隙, 建模孔隙度分别为10%、15%和20%.过程法建模结果孔隙分布较为集中，连通性较好，孤立点较少，与实际岩心相比虽不能完全吻合，但从三维孔隙形态上能表征实际岩心的微观结构的分布特征.从而搭建起从微观到宏观过渡的桥梁，这将为实现在微观领域研究岩石导电规律奠定基础.
|
图 2 过程法建立的三维数字岩心 Figure 2 The 3D digital core of process-based method |
有限元方法计算岩心导电性的基础是变分原理.如图 3所示，在数字岩心模型两端施加一个固定电场，随着时间的变化其能量逐渐趋向于稳定，而数字岩心模型的总能量由其结点电压决定，可将问题简化为求取系统总能量的极值问题.
|
图 3 数值模拟供电示意图 Figure 3 The schematic diagram of numerical simulation |
采用有限元法计算三维数字岩心电阻率的步骤如下：(1) 数字岩心模型的离散化处理.对求解区域进行离散剖分，将连续的求解区域转化离散单元；(2) 结点单元分析.在单元结点上对求解函数进行插值，将微分方程中各变量用相应的插值函数组成线性表达式；(3) 整体分析.即借助变分原理，求解整体系统能量极值问题.
2.2 数字岩心导电性各向异性分析孔隙水电阻率设为1，骨架电阻率取1015 Ω·m，使用有限元法模拟含饱和水条件下砂岩数字岩心模型的导电特性.为实现数字岩心模型的不同方向下的导电性分析，分别计算了孔隙度在三个主轴方向随电阻率变化规律，并求出岩石指数a和胶结指数m.
如图 4所示，在X、Y、Z三个坐标轴方向上，电阻率随孔隙度都呈指数变化规律，这与阿尔奇公式的表现形式一致.当孔隙度小于0.15时，在三个方向上电阻率出现明显的差异，Z轴方向上电阻率最高，X轴其次，Y轴最低.孔隙度在0.15~0.2区间时，Z轴方向电阻率处于X、Y轴方向之间.当孔隙度大于0.2时，电阻率变随孔隙度的变化趋于平缓，三个方向的电阻率变化几乎没有明显差异.
|
图 4 数值模拟结果对比图 Figure 4 Comparison of numerical simulation results |
究其原因，可能是由于当孔隙度较大时，孔隙空间较大且基本处于完全连通状态，孔隙水中的离子可以自由运动，这时电阻率随孔隙度的变化相对较缓，在三个主轴方向的电阻率差异较小.当孔隙度较小时，孔隙结构较为密集，流体通道狭窄，导电性在三个主坐标方向分化明显.
2.3 阿尔奇公式适用性分析为了验证所建模型对阿尔奇公式的适用性，本文把数值模拟结果与理论计算结果进行了对比分析.根据阿尔奇公式可得地层因素F与岩石孔隙度ϕ的关系为
|
(3) |
其中：Rt为含饱和水岩石的电阻率，Rw为地层水的电阻率，ɑ为岩性系数，m为胶结指数.
当地层水电阻率设为1时，公式(3) 可变为
|
(4) |
为了研究所建模型的各向异性特性，现分别求取X、Y、Z三个主轴方向的电阻率与孔隙度的拟合关系式，并得出胶结指数m和岩心系数a的数值，见表 1.拟合相关系数R2均在0.98以上说明拟合公式具有较高的可信度.
|
表 1 阿尔奇系数计算结果表 Table 1 The calculation results of Archie coefficient |
由岩石物理实验结果分析可知，胶结指数m随孔隙度、渗透率的增大而增高，通常粉砂岩m值较小，粗岩性的m值较大；岩性系数a与m关系密切，a值大，m值就小，反之，a值小，m值就大.表 1可以看出，在三个主轴方向上，a值升高，m值则会相应的降低，这与实验结果相符.岩性系数在三个主轴方向的变化范围在0.5~0.8之间，而胶结指数在2.3~2.7之间变化，这与a为1，m为2的理论结果较为接近，说明所建立的数字岩心模型在三个坐标主轴方向上存在一定的导电性差异，而用有限元法计算数字岩心模型导电性能较好地展现岩石的导电规律.
以上探讨了当地层水电阻率为固定值时，电阻率数值模拟结果随孔隙度变化规律的分析.但孔隙度若为固定值的同一数字岩心模型，在不同方向上其导电性随着地层水电阻率将如何变化是值得探讨的问题.公式(5) 为
|
(5) |
由公式5可以看出，当数字岩心模型和孔隙度值确定时，在三个主轴方向上的岩性系数a和胶结指数m也为固定值，此时岩心电阻率随地层水电阻率的变化呈线性关系.分别计算了孔隙度为5%、10%和20%时的数字岩心电阻率，结果如图 5、图 6、图 7所示，在X、Y、Z三个主轴方向上电阻率数值模拟结果与地层水均呈近似线性变化规律，这与阿尔奇公式的变化规律一致.此外对于相同的数字岩心模型在三个主轴方向上数值计算的电阻率随地层水电阻率变化有明显差异，其中在孔隙度为5%时三个主轴方向上电阻率分化最大，孔隙度为10%的模型其次，孔隙度为20%的模型电阻率分化最小.在整体上看所建模型导电性变化规律上与理论结果完全相符，只是存在局部差异，这说明采用过程法建立的数字岩心模型能反映实际岩心孔隙空间结构特征，使用有限元法模拟数字岩心模型的导电性具有较高精度.
|
图 5 数值模拟电阻率与地层水电阻率交会图(孔隙度为5%) Figure 5 The intersection graph of simulation results and formation water resistivity (Porosity 5%) |
|
图 6 数值模拟电阻率与地层水电阻率交会图(孔隙度为10%) Figure 6 The intersection graph of simulation results and formation water resistivity (Porosity 10%) |
|
图 7 数值模拟电阻率与地层水电阻率交会图(孔隙度为20%) Figure 7 The intersection graph of simulation results and formation water resistivity (Porosity 20%) |
本文采用过程法建立了三维数字岩心模型，所建模型孔隙连通性好、孤立点较少，能反应真实岩石孔隙空间分布特征，能为后续研究岩石导电特性奠定基础.
3.2使用有限元法计算所建数字岩心模型电阻率，结果显示在三个主轴方向上表现出一定的导电性差异，说明过程法所建的模型表现出各向异性.
3.3研究表明：所建岩心模型的电阻率随孔隙度呈指数变化，计算得到的岩性系数a的值接近1、胶结指数m的值接近2，并且随着模型的改变，a值增大，m值就减小，反之，a值小，m值就大，这与理论结果一致.说明采用过程法建立的数字岩心模型能与实际岩心孔隙空间结构特征相吻合，使用有限元法模拟数字岩心模型的导电性具有较高可靠性.
致谢 在本文研究和写作过程中，得到了吉林大学李舟波教授和莫修文教授的关心和建议，在此表示诚挚谢意.| [] | 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[J]. Geophysics, 67(5): 1396–1405. DOI:10.1190/1.1512785 |
| [] | Brownstein K R, Tarr C E. 1979. Importance of classical diffusion in NMR studies of water in biological cells[J]. Physical Review A, 19(6): 2446–2453. DOI:10.1103/PhysRevA.19.2446 |
| [] | Bryant S, Blunt M. 1992. Prediction of relative permeability in simple porous media[J]. Physical Review A, 46(4): 2004–2011. DOI:10.1103/PhysRevA.46.2004 |
| [] | Fatt I. 1956. The network model of porous media. Ⅲ:Dynamic properties of networks with tube radius distribution[J]. Trans. AIME, 207: 164–181. |
| [] | Hazlett R D. 1997. Statistical characterization and stochastic modeling of pore networks in relation to fluid flow[J]. Mathematical Geology, 29(6): 801–822. DOI:10.1007/BF02768903 |
| [] | Ma L. 2014. Numerical simulation of saturation parameters in gas hydrate-bearing reservoirs (in Chinese)[Master's thesis]. Changchun:Jilin University. |
| [] | Martys N, Garboczi E J. 1992. Length scales relating the fluid permeability and electrical conductivity in random two-dimensional model porous media[J]. Physical Review B, 46(10): 6080–6090. DOI:10.1103/PhysRevB.46.6080 |
| [] | 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[J]. Chinese J. Geophys. (in Chinese), 59(5): 1831–1838. DOI:10.6038/cjg20160526 |
| [] | Okabe H, Blunt M J. 2005. Pore space reconstruction using multiple-point statistics[J]. Journal of Petroleum Science and Engineering, 46(1-2): 121–137. DOI:10.1016/j.petrol.2004.08.002 |
| [] | Øren P E, Bakke S. 2002. Process based reconstruction of sandstones and prediction of transport properties[J]. Transport in Porous Media, 46(2-3): 311–343. |
| [] | Øren P E, Bakke S. 2003. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects[J]. Journal of Petroleum Science and Engineering, 39(3-4): 177–199. DOI:10.1016/S0920-4105(03)00062-7 |
| [] | Wu K J, Nunan N, Crawford J W, et al. 2004. An efficient Markov chain model for the simulation of heterogeneous soil structure[J]. Soil Science Society of America Journal, 68(2): 346–351. DOI:10.2136/sssaj2004.3460 |
| [] | 马龙. 2014. 天然气水合物储层饱和度参数数值模拟[硕士论文]. 长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1014268418.htm |
| [] | 莫修文, 张强, 陆敬安. 2016. 模拟退火法建立数字岩心的一种补充优化方案[J]. 地球物理学报, 59(5): 1831–1838. DOI:10.6038/cjg20160526 |