顺北油田某区块三维地应力场有限元精细建模与分析

陈修平， 沈新普， 刘景涛， 沈国晓

(1. 中国石化西北油田分公司 石油工程技术研究院， 乌鲁木齐 830001； 2. 中国石化缝洞型油藏提高采收率重点实验室， 乌鲁木齐 830001； 3. 中国石油大学(华东) 石油工程学院， 山东 青岛 266580)

顺北油田是塔里木盆地北部区域的主要碳酸盐储层油田.顺北油田所在区域有13条断裂带，地层裂缝较多.由于经历了多次地质构造运动，油田区块各处局部应力分布复杂，难以准确预测，从而导致无法准确计算安全泥浆窗口.实际钻井工程中井壁坍塌和泥浆漏失严重，发生钻井事故风险较高[1-2].本文的任务是结合顺北5号断裂带北部破碎性地层的单井实测信息及钻井工程中的阻卡情况，对目标区块三维地应力场进行精细化有限元建模计算，得到区块地应力的三维有限元空间分布数值结果，为钻井安全泥浆密度窗口进一步的选择提供理论依据和准确的参考数据.

区块地应力场的精准构建是近年来国内外工程力学的研究热点之一，若干研究者进行了不同的尝试[3-5]，采用不同的技术方法得到的地应力场各有特点.准确建立三维地应力场主要难点在于：1)主应力数值结果的大小与相应的实测值大小和实测值的方向不易同时满足，如果保证主应力分量大小的计算值和实测值一致，计算和实测结果的主方向则不容易保证一致；强调主应力计算结果和实测值的方向一致，则很难保证二者的大小一致.若干关于地应力场分析的文献中主应力结果与实测值的大小和方向均存在差别.2)碳酸盐储层的材料特性比较复杂，试验测量数据较少.

随着地质测量技术和计算技术的发展，近年来三维精细地应力场建模技术有了长足进步[6-8].目前解决上述问题的方法主要有：1)通过引入区块三维地震波数据，建立基于地震波的三维地质模型，在该三维地质模型的基础上建立三维区块有限元力学模型.2)进行单井地质力学分析，通过测井数据及钻井实测信息来计算地层材料参数，并确定所在井位上各个深度的地应力分量.3)根据单井地质力学分析结果和实测数据来建立区块地应力场的初始输入参数，并进行有限元数值计算，然后把三维有限元数值计算结果与单井实测结果相比较，根据两者误差精度情况改进输入参数取值，直到数值结果和实测结果两者精度位于合理范围.单井实测结果包括用地漏试验方法得到的最小水平主应力值，以及根据影像测井分析得到的最大水平主应力的方向.图1为三维精细地应力场的有限元建模与分析流程图，其中，Sh为最小水平主应力，SH为最大水平主应力.

图1 精细地应力场的技术路线流程图Fig.1 Flow chart of technical route for refined geostress field

1 顺北油田某区块三维有限元模型

图2为本文根据地震波数据并结合单井层位信息建立的顺北油田某区块的三维地质模型.模型中包括了11个地层，各层之间以地质分层界面为界，分别为：地表-T30界面；T30-T33界面；T33-T46界面；T46-T50界面；T50-T56界面；T56-T60界面；T60-T63界面；T63-T70界面；T70-T74界面；T74-T76界面；T76-模型底面，共11个地层.这里的T30至T76为根据地层沉积年代确定的层面.在该区块中，奥陶系碳酸盐油层的位置在T74-T76地层，深度约为7 000～7 600 m之间.

图2中区块地质模型长宽高尺寸分别为：长55 km，宽26 km，深度约9 980 m.地表海拔约为985 m，接近平坦.最下层是为了引入位移约束而加的一个层位.下底面设为平面.图2b为从地质模型转换而来的有限元网格.图2c为T74-T76奥陶系碳酸盐的油层网格.图中红点为A井目的层位置.有限元计算模型采用了127 765个C3D8R单元、137 088个节点对区块几何体进行离散.

2 单井地质力学分析

图3为A井的测井数据以及根据测井数据计算得到的地层参数.测井数据包括地层的伽马射线、压缩声波时长和剪切声速时长；地层参数包括杨氏模量、泊松比、密度.在计算杨氏模量和泊松比等地层参数时，采用了文献[9-11]给出的经验公式和计算原理.

图3 A井的测井数据和地层参数Fig.3 Logging data and stratum parameters of well-A

图4为A井单井的地应力分析结果.图4中各曲线分别为：紫线为垂向应力，蓝线为最小水平主应力，该最小水平主应力在5 400 m左右深度上的值与地漏试验点(红方块符号)相重合，说明该最小水平主应力曲线结果是合理的.细黑色曲线是最大水平主应力结果.粗黑色曲线为实际采用的泥浆密度曲线.在约4 800 m深度进行一次地漏试验，测得的漏失发生时的泥浆压力LOP用浅蓝色方块表示.在深度约7 500 m位置上有井漏发生，表明这两个位置上的地层中天然裂缝较多.在7 200 m深度上，发生过一次钻具阻卡，表明该位置上的钻井液/泥浆密度很接近坍塌压力，泥浆密度的安全系数不足.红色曲线代表计算得到的坍塌压力，即泥浆密度窗口下限.橄榄绿曲线代表地层孔隙压力分布曲线.

图4 A井的单井地应力分析结果Fig.4 Analytical results of single-well geostress for well-A

在距离地表1 000 m以内的范围内，由于种种原因，声波时长等测井数据缺少或精度不够，故1 000 m以上的单井分析结果不能使用.

3 区块地质力学有限元分析结果

使用图2中给出的三维有限元网格，将图3～4的单井分析结果和前述最大水平主应力方向等参数作为输入数据，建立了目标区块地应力场的有限元模型.该有限元模型的边界条件为四个侧面和底面的法向零位移约束.载荷为重力载荷，初始地应力参数的设置根据单井地应力分析结果确定.

为了体现弹性模量随深度变化的特性，把弹性模量的大小设为深度的函数，设计了用户子程序来实现这一模型功能特点.经过数值计算得到的三维有限元地应力分析数值结果如图5～7所示.

图5中最大水平主压应力方向接近正北方向.这与图6给出的T74-T76地层的最大水平主压应力方向不同，有超过20°的变化.该现象背后的主要原因是：不同的地层对应的造山运动方向不同，从而产生的构造应力方向也不同.

图5 上部T30-T33地层中最大水平主压应力的数值解Fig.5 Numerical solution of maximum horizontal principal stress in upper T30-T33 stratum

图6 T74-T76地层中最大水平主压应力方向Fig.6 Direction of maximum horizontal principal stress in T74-T76 stratum

图7 下部T74-T76地层中最小水平主压应力方向Fig.7 Direction of minimum horizontal principal stress in lower T74-T76 stratum

由于区块的应力模式属于“正断层应力”模式，即竖向应力为最大的主压应力分量，因此沿最小水平主压应力方向的钻井将具有最佳的井壁稳定性，也就是最佳的井轨迹方向.

图8为地应力场三维数值解的三个主应力分量在A井位置上随深度变化的分布曲线及其与A井的单井地应力分析结果的比较.图8中，Shminc为最小水平主压应力，Shmaxc为最大水平主压应力，Sv为竖向应力，3D为三维数值解，1D为单井解析解.

图9为图8结果在下部6 000～7 400 m深度段的地应力主分量解析解与数值解比较局部放大图.由图9可以看出，三维精细地应力场数值解在A井位置上与单井地应力分析数值十分接近，最大差别不到10 MPa，相对误差小于8%.图8中，由于上部网格划分得比较粗，对应的上部地应力数值很粗糙.

图8 地应力场三维数值解和单井解析解比较Fig.8 Comparison of 3D numerical solutions of geostress with their single-well analytical counterparts

图9 在下部深度上三维数值解和单井解析解比较Fig.9 Comparison of 3D numerical solutions at lower depth section with their single-well analytical counterparts

4 结 论

本文采用“综合三维地质体和单井地应力分析结果”的方法，为顺北油田某区块建立了三维精细地应力场.首先根据区块地震波数据信息建立精确的地质层位模型，并结合各个单井位置上的测井信息，校正/完善三维地震信息中得来的地层层位信息，得到完善的三维区块地质模型.该建模做法保证了结构变化引起的地应力方向变化能够得到准确体现.

由于区块模型进行地应力场有限元分析的主要目的是通过模型分析得到储层位置上精确的地应力场方向和主应力分量幅值的分布，结构模型最上面的上覆地层构造细节可以适当简化处理，以此来减小计算量.但是所对应的应力数值解在最上面地层中的分布比较粗糙.采用单井测井数据来计算地层材料和初始地应力相关参数，能够减少实测参数不足造成的输入参数不确定性，保证精细地应力场模型输入参数的准确性.

与实测结果相比，本文模型所得到的三维精细地应力场中最大主压应力的方向在A井位置上与实测值符合.最小水平主应力的数值解和A井位置上对应的解析解相比，误差小于8%，说明本文得到的精细地应力场数值解具有合理的精度.