基于有限元法的裂缝页岩计算机建模仿真研究

郭 晨，刘婉萍

(长安大学 信息工程学院，陕西 西安 710064)

0 引 言

随着世界能源需求量的不断增加，埋藏在浅层的常规石油资源已日渐枯竭[1]，寻找有开采价值的非常规资源成为重中之重的事情。页岩油页岩气是重要的非常规资源，其主要是指以页岩为主的页岩层系中所含的石油资源[2]。中国在四川盆地、苏里格等油田发现了储藏丰富的致密油气田，与此同时，随着勘探技术[3]的发展，页岩油气的开采成本也在不断降低。面对储藏丰富的页岩油气田，中国对于其开采和探测技术仍然处于摸索和依赖国外技术的阶段。

页岩作为一种致密沉积岩，是通过若干年的细微颗粒及薄层沉积形成的。大多数情况下，页岩在其水平方向可看作同一时期的沉淀物，在垂直方向可看作不同的沉淀物。因此可以将页岩看成不同介质组成的层状混合物。对于页岩混合物来说，在每层介质的形成过程中，由于各种地质结构以及氧化等原因导致在层与层之间形成很薄的薄层。此外，页岩在形成过程中，由于挤压等原因会形成裂缝[4]，这些裂缝会以不同形态存在于页岩的不同位置。这些薄层和裂缝的存在对页岩油气田的开采和探测具有一定的影响。因此国内外许多学者对此做了研究：美国Brace和加拿大Stesky等人对电测井响应进行了数值模拟，研究表明裂缝的存在会使地层的电特性出现显著变化[5-6]。2007年，中国石油大学的沈金松对水平裂缝进行了建模分析，分析了几种裂缝模型在不同围压条件下的电阻率变化规律[7]。2017年，唐巨鹏等人研究了预制裂缝角度对页岩水力压裂效果影响[8]。2018年，Matthew Josh和Han Tongcheng研究了定向裂缝对合成多孔砂岩介电特性的影响[9]。

在之前的研究中，对于电参数的求解问题，学者一方面通过目前已有的成熟的解析公式进行求解，另一方面利用实验室的测量手段进行求解。文中基于有限元法的计算机建模仿真，对不同几何结构的页岩模型进行仿真分析。对于裂缝页岩电特性的研究分析，学者只是对定向裂缝的电特性进行了研究，主要对层状结构、层状加薄膜结构以及含有不同形状、方向裂缝页岩结构进行了建模与仿真，分析不同结构对页岩电特性的影响并得出结论，结论为实际岩石的反演提供了理论依据，具有一定的应用价值。

1 有限元法

对于数值化分析方法，目前用于分析岩石电特性的主要有有限差分方法(FDM)[10-12]、有限元法(FEM)[13]、边界积分法(BIEM)[14]以及边界元法(BEM)[15]。文中主要用有限元法(finite element method，FEM)对页岩模型进行建模、仿真与计算。有限元法的核心内容是将一个连续三维几何体剖分成有限个面积不同、形状不同但是相连的几何体，即对求解的连续区域进行离散化，然后根据初始条件求解出每个剖分后区域电参数的近似解。最终，将每个离散区域进行总装，并利用随机法、迭代法、直接法等求解方法求解总装好的方程组，得出的解就是真解的近似解，若真解的近似解在误差范围内，则得出的近似解就是真解，若不满足，则需要重新计算，直至满足设计准则范围内的允许值为止。

有限元法的实现需要三个阶段：前期处理阶段、计算求解阶段和后期处理阶段。在前期处理阶段主要是对实际问题进行建模并剖分的工作，在后期处理阶段主要对结果的准确性进行评估分析，进而得到精确度高的解。

基于有限元的分析方法，国内外研究开发了许多功能强大的算法软件，文中主要利用COMSOL Multiphysics仿真软件进行建模与仿真。

2 仿真方法

有限元法的核心步骤分别为：确定物理场及几何区域；建立几何模型；有限元网格划分;确定初始条件;求解每个子域的近似解;联合离散域方程;求解等效电参数。

2.1 物理场选取

基于有限元法的COMSOL Multiphysics是一款高级多物理场的仿真软件，包含结构力学、声学、光学、电化学以及电磁分析等物理场。根据不同的实际问题选择合适的物理场，文中主要采用电磁方法对裂缝页岩模型的电学特性进行了建模及分析。因此，在仿真中，物理场都选取电磁波场。

2.2 仿真模型

由于页岩具有明显的层理结构，并且在每层形成过程中，介质表面会形成极薄的薄膜层或者由于挤压等原因会产生裂缝层。因此在文中主要有三种页岩结构模型：层状结构、层状加薄膜结构和椭球形内含物结构。图1是结构模型示意图。

图1 三种页岩结构模型

2.3 网格划分

图2是三种不同的网格划分示意图。

不同的网格划分影响最终求得的解的精度，网格划分越细，解的精度越高，所需的计算机内存越大，耗费的时间越长。为了保证计算解准确的情况下能快速计算页岩电参数，经过大量的仿真计算，利用标准网格进行剖分能快速有效地计算电参数。因此，文中在仿真中，都是用标准网格进行剖分。

图2 三种网格剖分图

2.4 等效电参数求解

最终仿真计算真解时，采用牛顿迭代法进行方程组计算，当相对容差小于10-3时终止计算。仿真计算结束后，通过COMSOL的后期处理功能，计算不同区域、不用方向的电场模量，根据以下公式求得不同页岩结构的等效介电常数。

(1)

其中，D是电通量密度，E是电场强度。

3 仿真结果

由于页岩在地形中的分布特点，主要对三种页岩结构进行建模仿真，其中包括层状结构、层状加薄膜结构和椭球形内含物结构。

3.1 层状结构

在地层中，由于页岩具有很明显的层理结构，在仿真建模中可以等效为层状模型，为了验证仿真结果的正确性，将COMSOL仿真结果与复折射率模型解析公式进行了对比，复折射率公式如下：

(2)

其中，ε1、ε2分别是上下两层介质的相对介电常数；d1、d2分别是上下两层介质的厚度，d=d1+d2。

模型尺寸及材料设置：两层立方体的长宽高分别为1 mm、1 mm、0.5 mm，下层介质的介电常数为20，改变上层介质的介电常数，结果如图3所示。

图3 COMSOL仿真结果与解析解对比

从图中可看出，COMSOL仿真结果与解析解的结果相吻合，由此证明了该模型与方法的正确性与有效性。

3.2 层状加薄膜结构

页岩具有明显的层理结构，但在每一层介质的交界面上会出现一层类似于薄膜的结构。对于含有薄膜结构的页岩来说，薄膜的厚度仅仅是整个模型尺寸的0.1%，或者更少。但是这种薄膜层对模型的等效介电常数影响极大。因此，文中对含有薄膜层结构进行了建模仿真。模型的尺寸及材料设置为：两层立方体的长宽高分别为1 mm、1 mm、0.5 mm，薄膜厚度为10-7mm，上层介质的介电常数为20，电导率为0.001 s/m，下层介质的介电常数为20，电导率为0.01 s/m，薄膜的介电常数为80，电导率为10-5s/m。结果对比如图4所示。

图4 含薄膜层与不含薄膜层仿真结果对比

由图中可看出，含有薄膜层结构模型的仿真结果在低频处出现了巨大增强现象，而在高频处和不含薄膜层结构的仿真结果基本一致。

3.3 椭球形内含物结构

页岩在实际形成过程中，由于受到挤压等原因产生一系列裂缝，这些裂缝以不同的形状、大小分布在页岩中，这对页岩的等效电参数带来了一定的影响，进而加大了对页岩的开采难度。在众多的裂缝结构中，椭球形裂缝结构最为常见，其中针状裂缝和圆盘状裂缝最具代表性。因此，文中对不同孔隙度下，不同角度分布的含裂缝页岩进行了建模仿真。模型如图5所示。

模型设材料置：立方体内包裹着椭球体，立方体的介电常数为4，电导率为0.01 s/m，椭球的介电常数为80，电导率为0.1 s/m。不同裂缝三方向上等效介电常数仿真图如图6所示(其中，图(a)和(b)的孔隙度为2.09%，图(c)和(d)的孔隙度为1.68%，图(e)的孔隙度为0.523%，图(f)的孔隙度为8.38%)。

图5 椭球不同角度(0°30°45°60°90°)分布示意图

图6 不同裂缝三方向上等效介电常数仿真图

由图6可以看出，含椭球形裂缝结构三方向上的等效介电常数值不相同，这主要是由于椭球在三个方向上的尺寸不同，从而对电场的阻碍程度不同所导致的。针状裂缝Y方向上的等效介电常数小于X、Z方向上的等效介电常数值，而圆盘状裂缝Y方向上的等效介电常数则大于X、Z方向上的等效介电常数值。针状裂缝和圆盘状裂缝在与水平方向的角度呈45°时，X、Z方向上的等效介电常数值相等。

综上所述，在宽频带内，若低频处的等效介电常数出现了巨大的增强现象，则在具有层理结构的页岩中夹杂着极薄的薄膜层；若三方向上的等效介电常数值不相等，则具有层理结构的页岩内部具有不规则的内含物，大多数情况下，等效其为椭球内含物；当Y方向上的等效介电常数小于X、Z方向上的等效介电常数值时，则在具有层理结构的页岩中包含着针状裂缝；当Y方向上的等效介电常数大于X、Z方向上的等效介电常数值时，则在具有层理结构的页岩中包含着圆盘状裂缝。研究结果可为实际岩石的反演提供理论依据。

4 结束语

主要运用有限元法对不同页岩内部几何结构进行建模与仿真，不仅解决了理论公式和实验测量无法进行内部几何结构的建模，而且形象准确地计算了不同几何结构页岩的电参数，根据其电特性规律总结出不同几何分布规律。仿真结果表明，含有薄膜层结构的页岩其电特性在低频处的等效介电常数出现了巨大的增强现象，含有椭球型裂缝结构的页岩其电特性在三方向上的等效介电常数值不相等。这些结论为实际页岩的反演提供了理论依据，具有一定的应用价值。