赵兵 , 罗攀登, 张妹珠, 周志敏, 黄欣

(1.中国石油化工股份有限公司西北油田分公司, 乌鲁木齐 830011; 2.武汉大学科学技术发展研究院, 武汉 430072;3.武汉大学土木建筑工程学院, 武汉 430072)

随着油田勘探和开发的不断加深和扩展,碳酸盐岩储层逐渐进入了人们的视野,该种类型油藏具有储量大、分布广等优点,已成为了人们开采和研究的重点对象[1-2]。塔里木盆地的顺北油田是新开发的储量较大的碳酸盐岩断溶体储层之一,储层中的断裂破碎区域如天然洞穴、天然裂缝内储藏有充足的油气资源[3]。针对断溶体储层的特点(孔隙度低、渗透率低、非均质性强等)提出水力压裂开采技术,该技术通过增加储层中裂缝和溶洞之间的有效沟通路径,降低分隔性,从而实现油田增产[4]。研究表明,水力压裂过程中的水力裂缝的扩展轨迹与地应力具有很强的相关性[5]。顺北油田属于碳酸盐岩断控裂缝-洞穴型储层,该油田主要位于交汇断裂带区域,地质条件复杂。为了利用水力裂缝沟通尽可能多的缝洞储集体以实现增产,则必须掌握储层三维地应力场空间分布规律。因此,如何准确获得顺北油田三维地应力场的精细空间分布规律,成为利用水力压裂技术进行油藏开采增产的关键问题。

在地应力测量和分析方面,现场地应力测量方法根据测量原理可分为[6]：一是破裂法,如水压致裂法、声发射法和井壁崩落法;二是变形法,如孔径变形法和孔壁应变法等;三是地球物理法,包括微地震波法和地面电位法等;四是构造法,如宏观构造反推法、显微构造分析法。在实际情形中,地质体中的缝洞储集体以及破裂带会对地应力场分布情况产生较大影响,现场实测的地应力值存在一定的离散型,仅依靠现场测量无法全面准确地反映出地质体地应力场的空间分布。因此,目前许多学者采用有限元软件对地层进行模拟计算,以精细刻画地应力场空间分布。王爱国等[7]利用子模型法提取目标区域边界条件,在构造地应力场的尺度对黑山峡大柳树坝址区进行数值模拟分析计算,描述了地震断层处地应力场特性,对坝体破裂危险性进行了分析评价。白克新[8]针对含断层地层巷道覆岩应力分布分析问题,根据弹性力学理论,利用数值模拟软件建立了平面模型。他对存在断层及正常条件下地层的应力特征及位移规律进行了比较分析,并利用物理相似模拟方法获取开采工作面接近断层过程中的岩层破裂、垮落特征以及地层应力演化规律,提出了回采巷道的合理支护方案。李万才等[9]以某坝址岸坡作为研究对象,建立三维数值模型,分析了在河谷下切影响下地质体应力、位移场特征以及塑性区的演化过程。然而,对于地质结构复杂的油气藏储层来说,由于有限元等软件的建模能力较弱,难以建立精细的三维储层模型。

在石油工程中,常利用地质建模软件来构建精细的三维储层模型,辅助人员进行开采策略调整。在地质建模方面,Petrel软件作为一款的地质建模软件被广泛应用于油气开采行业。该软件采用角点网格构建三维地质模型,能够精细地刻画出地质体中的断层、侵入、褶皱、尖灭等地质构造情形,同时还整合地质体的地震监测数据、油藏属性等储层信息。在实际应用中,胡向阳等[10]为解决碳酸盐岩缝洞型油藏三维地质建模的难题,提出了多元约束三维地质建模方法。该方法能够表征了缝洞储集体的空间分布,并为该区的增储建产提供有效保障。孙先达等[11]为了实现营城组升平地区火山岩储层的三维可视化动态表述和展示,通过Petrel软件利用构造层面及断层数据建立了构造模型和断层模型。在此基础上,他们采用确定性建模和随机建模结合,构建了该区火山岩相三维地质模型,实现在三维空间上刻画了典型火山岩体的岩相特征,并研究其在三度空间上的变化规律。

为了满足地质-工程一体化的需求,目前许多学者利用地质建模软件完成了数值仿真模型的构建。郭颖星等[12]为解决包含断层的复杂地质构造建模问题,开发了Petrel与Adina等软件网格数据之间的接口程序。该程序可以实现地质模型网格数据与仿真软件网格数据的相互移植,从而实现有限元模型的精确建模。刘钰洋等[13]开发了一系列算法,并将其集成到Petrel2ANSYS中,实现了Petrel中角点网格的3D属性模型与ANSYS中有限元网格3D模型之间的双向转换。徐珂等[14]构建了高尚堡油田深层油藏南区三维岩石力学场,通过Petrel地质模型与ANSYS有限元模型联合建模,对该地区三维地应力场空间分布规律开展了研究,并为油藏开发过程中压裂优势区和优势井段的优选提供了理论依据。Park等[15]开发了PET2OGS一套3种算法来集成静态模型(Petrel)和动态模型(OpenGeoSys),转换有限差分法(FDM)网格为有限元法网格。

目前对油藏储层的应力场研究主要集中于油藏区域的整体应力场分布,没有针对缝洞型储集体油藏地应力场的研究,缺乏对断溶体储层区域及内部、不同断裂带交汇区三维地应力场分布规律的分析。AiFrac是适用于水力压裂、油气开采等领域的多物理场耦合作用岩体裂缝扩展软件[16]。现通过研发“地质-工程”耦合平台,将地质模型转化为AiFrac仿真软件使用的模型,完成地质模型和数值仿真一体化,以实现三维地应力场和水力压裂的精细化模拟,形成实时驱动的储层地质和压裂过程数字孪生平台。研究成果将为油气开采提供可靠的决策支持和优化方案。

1 “地质-工程”耦合平台

虽然Petrel等地质建模软件能够构建三维地质模型,对储层进行精细刻画并描述储层的空间分布和属性分布特征,但在计算地应力场方面效率较低,不能够胜任大型三维地应力场的模拟。相比之下,有限元等数值仿真软件相较于地质建模软件则具备很强的大型地应力场计算能力,但在三维精细化建模方面能力较弱。因此,将数值仿真软件与地质建模软件结合使用是十分有必要的,这种方法既可以精细表征储层特征,也可以准确地模拟的地应力场分布规律。利用“地质-工程”耦合平台,将地质建模软件与数值仿真软件AiFrac整合[17]。该平台将地质模型转换成AiFrac计算所需的几何数据,并将储层地质参数信息融入AiFrac网格数据中,从而实现三维地应力场的精细化仿真。

1.1 地质模型

三维地质建模软件能够创建三维构造框架模型,利用地球物理和地质知识库综合分析油藏开发的全过程,帮助人们理解复杂的地质问题。地质建模常常采用的网格为角点网格[18]。角点网格是由4个坐标线以及8个网格节点构成的网格体,该网格体内部又可以分成i×j×k个子单元(i、j、k分别为角点网格在X、Y、Z方向上的数量),且相邻子单元的相邻界面节点坐标可以不完全相同,以作为断层或者裂缝。对于单体结构单元,4条坐标线中的任何一条都可以减少到零长度,即网格单元存在着退化,一般退化的网格会出现在侵蚀等复杂地质构造处。其中,网格体在平面上可以分为i×j个子平面网格,并且单元网格的长、宽大小可变[18]。在平面外方向上将连接顶底网格点的网格面分成k层,并且该网格面可以是倾斜的。每个网格节点的坐标值可以根据顶底边界的节点坐标值进行内插计算得出。则A、E两点的X、Y坐标,可由顶部I点与底部J点的X、Y、Z坐标值通过线性内插的方式得到。

(1)

(2)

式中：A点坐标为(XA,YA,ZA);E点的坐标为(XE,YE,ZE);I点的坐标为(XI,YI,ZI);J点坐标为(XJ,YJ,ZJ)。

在建模时,首先,利用井分层点数据插值生成构造层面,再依据断层模型生成包含顶面、中面以及底面的空间格架网格模型,然后,将地层层面信息加入其中,最后,通过设置垂向单元的厚度、单元的个数进行三维网格的创建以及地层网格的精细化。

1.2 AiFrac数值仿真模型

采用的软件AiFrac基于有限-无网格法(finite element-meshfree method, FEMM)原理来模拟三维水力扩展和地应力场[19-20],该方法最早由Rajendran等[21]于2007年首次提出。有限-无网格法中,根据单元与裂隙的位置关系,将单元分为3种类型,即常规的普通单元、桥单元和裂隙单元,如图1所示。被裂隙穿过的单元定义为裂隙单元,与裂隙单元相邻的为桥单元,其余的单元为普通单元。同时,该方法定义两种节点,裂隙单元上所有节点为单位分解节点,其他节点为普通的有限元节点。

图1 有限-无网格法中的3种单元和两种节点

对于非裂缝的普通网格单元,采用单位分解法(partition of unity method, PUM)构建位移函数。对于一个四面体单元,其计算域Ω内位移函数由4个节点的权函数和4个节点的位移得

(3)

式(3)中：wi(x)为一系列非负的权函数;i为与节

对于裂隙穿过的裂隙单元,采用Shepard公式构造单元节点权函数[22],可表示为

(4)

式(4)中：φ′ι(x)为形函数构成的子权重函数。

(5)

式(5)的含义是一个被裂缝分割的单元只有节点都处在裂缝面的一侧才认为节点彼此可见。因此,对于考虑裂隙面穿过的任意单元中的计算点Pi,可以用式(6)计算φ′i。

(6)

式(6)中：φi可以由任意坐标为x∈ψΩ的计算点Pi的体积坐标计算得到,可表示为

(7)

式(7)中：v(P1P2P3P4)为顶点P1、P2、P3和P4的四面体体积;v[P(x)PiPjPk]为由单元中任意点P(x)和3个顶点Pi、Pj、Pk组成的四面体体积。

对于桥单元,权函数采用有限元中的形函数形式。

1.3 地质模型和AiFrac模型耦合

“地质-工程”耦合平台采用C++语言编写。C++语言同时拥有高级语言和汇编语言的优点,同时具有简洁灵活、数据结构丰富和程序执行效率高等特点,保证了平台能够快速高效地将地质模型中的角点网格单元转换为可以供AiFrac计算使用的四面体网格单元。主要使用步骤如下。

步骤1使用地质建模软件中的“Eclipse Keywords format”格式将储层地质体模型数据导出。

步骤2利用“地质-工程”耦合平台的读入功能,读取“.GRDECL”后缀名的数据文件,并获取SPECGRID、COORD和ZCORN等网格坐标参数,以及PERMX、PORO和ACTNUM网格属性信息,再对读入的数据进行处理,生成AiFrac仿真模型所需的“.inp”后缀名数据文件。

步骤3再将数据文件导入到AiFrac数值仿真软件中,实现地质模型和数值仿真的一体化。

“地质-工程”耦合平台运行逻辑如图2所示,具体如下：①检查输入模型的有效性并读取文件,解析角点网格模型数据,记录角点网格模型的节点、网格和属性信息;②判断构建的网格模型单元的有效性,将无效单元删除,同时移除无效节点,然后更新节点单元关系;③通过啮合重复节点,减少生成的网格中重复节点的数量,同时将角点网格转换成有限元网格;④将四面体网格的几何、属性信息输出,形成可供AiFrac仿真软件运行的输入文件。

图2 “地质-工程”耦合平台运行逻辑

2 顺北断溶体储层精细化建模

2.1 基本地质特征

选取塔里木盆地顺托果勒隆起构造带中的顺北5号断裂带部分区域作为研究对象,该区域的地质构造主要受断裂带的控制,储层类型为裂缝-洞穴型碳酸盐岩储集体。此外,顺北地区的主要目标层位为一间房组和鹰山组,在早期拉张环境下,该地区沉积了厚达3 000 m的寒武系-中奥陶统台地相碳酸盐岩地层。

2.2 储层力学参数的确定

通过现场采样进行室内岩石力学实验,以及处理分析测井资料,可以获取弹性模量、泊松比等岩石力学参数。以往的研究成果表明：通常情况下,断裂带的弹性模型等参数相比储层同性地层会有一定程度降低,而泊松比略大。目前还没有确定断裂带岩石力学参数的方法,一般可以将断裂带的弹性模量设置为储层同性地层杨氏模量的50%～70%[23-24]。断裂带的弹性模量和泊松比与其复杂程度有关,复杂程度越大,其弹性模量通常会越小,而泊松比则会越大。基于前期获得的岩石力学实验和测井数据,合理选取岩石力学参数[25]。具体岩石力学参数如表1所示。

表1 计算模型的岩石力学参数Table 1 Rock mechanics parameters of the calculation model

2.3 边界条件和载荷确定

水平最大主应力方向一般是利用电成像、偶极声波资料、古地磁及波速各向异性等方法获得。由于断裂构造、岩体结构和性质不同,造成地应力方向发生了显著变化,通过平均各单井的水平最大主应力方向数据,确定研究区域现今水平最大主应力方向为NE41.4°[25]。

现今地应力大小主要可以通过室内岩心实验、测井数据和水力压裂法等方法获得,通过对前期完成的单井实测数据取平均值,确定研究区的边界荷载分别为水平最大主应力177.78 MPa和水平最小主应力为137.01 MPa,垂向应力通常等于岩石的上覆岩层压力,可以利用密度测井估算获得[25]。

2.4 基于地质数据体的三维仿真模型建立

基于油田构造特征、三维地震解释和测井数据等工程勘测资料,利用地质建模软件建立断裂带储集体三维地质模型。地质模型采用平面分段、纵向分层的思路进行建立。具体而言,首先,基于地震解释等工程资料,采用地震解释和地质研究相结合的方式,对储集体的断层以及断裂组合特征进行描述;然后,基于地震数据,建立断裂带储集体空间模型;最后,结合试井、岩心、测井数据,建立储集体属性模型。由于断裂带岩体的孔隙度与储层原岩有较大差异,因此,可以通过孔隙度来表征断溶体,并将其与储层原岩进行区分。三维地质模型如图3所示。图3(a)和图3(b)分别为地质模型中的溶洞分布图和断裂带分布图,比例尺表征的是水平方向上的比例关系。

基于前期已构建的顺北5断裂带的地质模型,利用“地质-工程”耦合平台,构建基于地质数据体的三维仿真模型。由于研究区域的地质数据体三维模型尺寸大、网格数量多,使用精细数值仿真软件直接对原尺寸网格的模型进行计算会面临耗时较长、内存占用多等问题。因此,在对地质数据体三维模型进行网格转换前,会对模型网格进行适当的粗化。为了方便加载模拟真实的边界条件情形,将地质模型的周边未建模区域进行了适当填充,对其进行加框处理,从而生成一个长方体形状的数值仿真计算模型。模型边界方向分别对应实际研究区域水平最大主应力方向和水平最小主应力方向。采用1.3节中所介绍的方法,首先,将地质模型按照“Eclipse Keywords format”格式导出角点网格的几何和属性数据,然后,利用“地质-工程”耦合平台将数据转换生成后缀为.inp格式的AiFrac仿真输入文件,最后,开始AiFrac仿真计算。模型加载示意如图4所示,在水平方向上施加挤压应力。在垂向上施加重力。本次研究所建立的AiFrac三维仿真模型包含1.787×107个三维四面体单元、4.78×106个节点。

图4 模型加载示意图

3 三维地应力场分布特征

通过AiFrac模拟计算,可以得到顺北5号断裂带断溶体储层地应力场分布规律。图5(a)为最大主应力大小分布,可以发现,现今最大主应力主要以挤压应力为主,在远离断裂带的地层处,地应力场均匀分布,最大主应力的数值主要范围为170～180 MPa,而在断裂带的内部,最大主应力会偏小,其数值范围为小于170 MPa。

图5 储层区域主应力分布

图5(b)为中间主应力大小分布,可以发现,现今中间主应力同样以挤压应力为主,在远离断裂带的地层处,地应力很长均匀分布,中间主应力的数值主要范围为170～180 MPa,比最大主应力数值略微偏小,而在断裂带的内部,最小主应力会偏小,其数值范围为小于170 MPa。

图5(c)为最小主应力大小分布图,可以发现,现今最小主应力同样以挤压应力为主,在远离断裂带的地层处,地应力场均匀分布,最小主应力的数值主要范围为133～137 MPa,而在断裂带的内部,最小主应力会偏小,其数值范围为小于133 MPa。

图6为三轴主应力大小剖面图,图6(a)、图6(b)、图6(c)和图6(d)分别为剖面位置、最大主应力、中间主应力和最小主应力。可以发现,与前面的平面三轴主应力分布规律基本一致,在断裂带内部,三轴主应力大小偏小,而在断裂带端部及断裂带交汇,三轴主应力偏大。同时,由于垂向应力增加,三轴主应力大小整体上随着地层深度的增加而增加。

图6 储层区域应力场分布剖面

断裂带等非连续结构面对地应力的大小和方向都有着显著的影响,这是影响顺北地区地应力场分布的主要因素之一。由于岩石力学参数弱化,断裂带内部应力显著降低,产生明显应力释放现象。同时,断裂带端部附近则产生了明显的应力集中现象。另外,交叉等类型的断层交汇区也容易发生应力集中现象。此外,断裂带的规模越大,引起的地应力扰动带范围越大,同时,断裂带的倾角、走向,也会对地应力场的分布产生影响。综上所述,顺北断溶体储层区域及内部、不同断裂带交汇区三维地应力场分布差异明显,断裂带的规模、倾角、走向以及形态以及等因素在不同程度上均对地应力状态产生影响。

4 结论

针对顺北油田断溶体储层区域及内部、不同断层交汇区三维应力场精细空间分布刻画的难题,通过开发“地质-工程”耦合平台,提出了一种精细化三维地应力场仿真模拟方法,该方法实现了地质模型和数值仿真的一体化。此外,基于“地质-工程”耦合平台,开展了顺北碳酸盐岩储集体三维地应力场分布规律研究。得到以下主要结论。

(1)通过开发“地质-工程”耦合平台,可以将基于角点网格的地质模型信息转换到AiFrac数值仿真软件中,实现了地质模型和数值仿真一体化。该平台结合了两种方法的优势：一是克服了地质模型软件无法和力学模拟的软件衔接问题,实现了地质建模和地应力场反演的一体化;二是在三维地应力场数值模型中保留了地质模型的储层特征。可以快速和准确地将包含断层等复杂地质构造的地质模型转换成三维地应力场反演数值模型,保留地质模型中的网格和属性信息。

(2)以顺北油田为工程实例,基于开发的“地质-工程”耦合平台,将地质模型转换为AiFrac仿真模型,同时确定储层力学参数以及边界条件,实现了对碳酸盐岩储集体三维地应力场的精细研究。研究结果表明,断裂带等非连续结构面对三维地应力的大小和方向都有着显著的影响,这是影响顺北地区三维地应力场分布的主要因素之一。断裂带处材料参数的弱化,会导致在断裂带内部出现应力释放,从而引起断裂带内部应力降低,同时在断裂带端部附近则产生明显的应力集中现象。顺北断溶体储层区域及内部、不同断裂带交汇区三维地应力场分布差异也很明显,断裂带的规模、倾角、走向以及形态等因素在不同程度上都会对地应力状态产生影响。