CN115455770B - 一种基于子结构的油气田井筒结构应力精细模拟方法 - Google Patents

Abstract

本发明提供了一种基于子结构的油气田井筒结构应力精细模拟方法，在粗网格模拟中，选取井筒所在的单元作为母单元，进行局部化的子结构建模并求解，推导了子结构的刚度矩阵和荷载列阵，并将子结构的刚度矩阵根据单元连接性进行了子块划分，使子结构既能与母单元完全衔接，又能体现套管、水泥环等精细结构。本发明能够在粗网格模拟的基础上精细的反应套管、水泥环处的应力，用较低的计算量实现了较高的计算精度。

Description

一种基于子结构的油气田井筒结构应力精细模拟方法

技术领域

本发明属于油气藏数值模拟领域，涉及一种基于子结构的油气田井筒结构应力精细模拟方法。

背景技术

我国各油、气田井筒损坏十分严重，国外同样也存在井筒破坏的情况，根据近30年来国内外的文献研究和近10年来的动态跟踪研究，发现井筒损坏是国内外长期存在的问题，而且一直没有得到解决。对于长期开采的油气田而言，长期的注采改变了油气藏初始地质环境，造成地层压力的不平衡，同时注入流体于原有油气在油气藏孔隙中渗流动态过程中，对地层地应力场和变形场会产生影响，注入流体对储层力学特性产生影响；注入流体进入断层，对断层滑移的诱发作用等，都使得油气藏地质条件变的更加复杂，油气井井筒的受力情况更为复杂。

井筒破裂的形态可以分为错断、裂开、腐蚀穿孔等。错断是最严重的一种井筒变形形态，井筒在水平方向错断开，断开处附近伴随弯曲，这种形式的变形主要由于井筒受强大的剪应力造成；井筒裂开是由于射孔或井筒钢材本身缺陷造成，油田浅层水酸性较高时也会造成井筒腐蚀破裂。井筒密封性破坏主要表现在井筒连接处，由于拉伸造成脱扣及井筒丝扣质量原因导致套外返油、气、水。井筒变形的检测技术与方法为掌握井下井筒变形的形态，更好的研究井筒变形损坏机理，并且为套损防控、气井耐压力评价提供了可靠的技术资料。

鉴于用物理模型实验，很难模拟复杂地层井筒受力和变形过程，因此，用计算机数值模拟技术为主要手段及现代地应力场的研究是今后复杂地层井筒损坏机理与预防措施研究发展方向。

发明内容

本发明的目的在于提供一种基于子结构的油气田井筒结构应力精细模拟方法，解决了目前物理模型实验，很难模拟复杂地层井筒受力和变形过程的问题。

本发明采用的技术方案是：一种基于子结构的油气田井筒结构应力精细模拟方法，包括以下步骤：

步骤1：在油藏数值模拟的粗网格中，提取油气井筒所在的单元作为母体单元；

步骤2：对母体单元进行精细划分，生成由套管、水泥环逐步过渡到母单元的精细网格，形成含套管-水泥环系统的子结构；

步骤3：建立子结构系统的单元刚度矩阵和荷载列向量；

步骤4：提取粗网中数值模拟的结果，采用形函数插值法对子结构精细网格内的各单元进行初始化；

步骤5：对子结构开展有限元计算，获取精细网格内的应力、变形及结点位移数据；

步骤1包含以下步骤：

步骤101：根据油气井的中心坐标(x_w,y_w),在粗网格中查找平面坐标包含点(x_w,y_w)的单元作为母单元；

步骤102：提取母单元的8个节点坐标(x_i,y_i,z_i)i＝1,...,8；

步骤103：计算母单元各个边的边长l_i i＝1,...,8，取单元的特征长度

步骤2包含以下步骤：

步骤201：设置套管的壁厚h_c和套管层的单元数n_c，计算套管层的单元厚度t_c＝h_c/n_c。设置水泥环的壁厚h_f和水泥环层t_f的单元数n_f，计算水泥环层的单元厚度t_f＝h_f/n_f；

步骤202：设定子结构内径向单元尺寸的增长系数k，设定环向单元划分数m；

步骤203：计算子结构径向划分的层数N：

步骤204：计算自套管向外第i层单元的径向尺寸t_i和环向尺寸r_i:

步骤205：根据各层单元的径向尺寸t_i和环向尺寸r_i形成子结构平面各层的精细网格；

子结构内的单元为自套管向外的辐射状网格，内层网格较为精细可以反应套管、水泥环系统；外层网格较大，可以和粗网格衔接。

步骤3包含以下步骤：

步骤301：计算子结构内各单元的刚度矩阵[h_i]；

步骤302：根据子结构内各单元的连接关系，将各单元的刚度矩阵[h_i]组装，集成为子结构的总体刚度矩阵[H]；

步骤303：分析子结构内节点的连接关系，将节点划分为与母单元共享节点、与母单元不直接相连的内部节点以及连接共享节点与内部节点的过渡节点3部分；

步骤304：采用高斯-约当消去法，对总体刚度矩阵[H]进行矩阵运算，得到刚度矩阵[H]分块表达式：

式中第1列的子块与共享节点相对应，称为出口刚度矩阵；第2列与内部节点相对应，称为内部刚度矩阵；第3列与连接节点相对应你，称为连接刚度矩阵。

步骤4包含以下步骤：

步骤401：对子结构中的单元E_i，计算出其第j个积分点的坐标J(i,j)；

步骤402：在母单元中，根据坐标J(i,j)，插值计算该点处的应力σ(i,j)，压强p(i,j)，温度T(i,j)；

步骤403：将插值得到的应力σ(i,j)，压强p(i,j)，温度T(i,j)赋给单元E_i的第j个积分点。

步骤5包含以下步骤：

步骤501：结合步骤304中的刚度矩阵[H]，分析单元荷载，形成有限元方程组：

[H]{φ}＝{R}

其中{φ}为待求的位移组成的向量；

步骤502：求解有限元方程组，得到各节点的位移；

步骤503：计算子结构内各单元的应力。

本发明的有益效果在于：本发明提供了一种基于子结构的油气田井筒结构应力精细模拟方法。传统的粗网格模拟中，由于单元尺寸过大，导致套管、水泥环等细微结构被忽略而难以精确的模拟的套管、水泥环等井身结构的应力；传统的细网格模拟中，由于要体现套管、水泥环等井身结构，单元尺寸较小导致单元数量巨大而难以计算。本发明提供的基于子结构的油气田井筒结构应力精细模拟方法，在粗网格模拟的中，选取井筒所在的单元作为母单元，进行局部化的子结构建模并求解，能够在粗网格的基础上精细的反应套管、水泥环处的应力，用较低的计算量实现了较高的计算精度。有效的解决了传统粗网格模拟中套管、水泥环等井身结构被忽略的问题和传统细网格模拟中单元数量巨大难以计算的问题。

附图说明：

图1示出了根据本发明的一个实施例的基于子结构的油气田井筒结构应力精细模拟方法的流程图。

图2示出了根据本发明的一个实施例的子结构生成过程示意图。

图3示出了根据本发明的一个实施例的从粗网格中提取含井筒结构的母单元的示意图。

图4示出了根据本发明的一个实施例的子结构的平面划分图。

图5示出了根据本发明的一个实施例的子结构三维图。

图6示出了根据本发明的一个实施例的子结构计算的应力分布图。

图7示出了根据本发明的一个实施例的子结构内套管和水泥环的应力分布图。

具体实施方式

下面结合实例和附图对本发明作进一步说明。

下面对本发明的建模方法具体介绍如下：

一种基于子结构的油气田井筒结构应力精细模拟方法如图1所示，所述方法具体包括如下步骤：

步骤1：在油藏数值模拟的粗网格中，提取油气井筒所在的单元作为母体单元。

对于实施例，粗网格如图2所示，具有300×200×15个单元格，待提取的注水井的平面坐标为(1320.5，2764.5)。粗网格中每个单元格均为6面体，具有8个节点，单元格的形状和节点编号规则如图2所示。

其具体的实现步骤为：

步骤101：根据油气井的中心坐标(1320.5,2764.5),在粗网格中查找平面坐标包含点(1320.5,2764.5)的单元作为母单元。提取过程如附图3所示，具体实施方法为：对粗网格中每一个单元格k,如果其8个节点坐标

满足以下条件，则选中该单元格作为母单元：

式中，x_w为井点平面位置的横坐标，对于本实施例有x_w＝1320.5；y_w为井点平面位置的纵坐标，对于本实施例有y_w＝2764.5。

步骤102：提取母单元的8个节点坐标(x_i,y_i,z_i)i＝1,...,8。

经过步骤102，提取到的13个母单元，以母单元E1为例，各的节点坐如表1所示，其余各母单元的节点坐标分布于此相似。

表1母单元E1各节点坐标

步骤103：计算母单元各个边的边长l_i i＝1,...,8，取单元的特征长度

采用103的方法，计算得到的母单元E1的特征长度为2.05m。

步骤2：对提取的母单元进行精细划分，生成包含套管-水泥环系统的子结构，本实施例中，所采用的套管水泥环系统的结构如图4所示。对步骤1中提取的各母单元，具体实施如下：

步骤201：换分套管和水泥环单元。本实施例中，套管壁厚为h_c＝18.3mm，沿径向套管单元数n_c＝5，则套管层的单元厚度t_c＝h_c/n_c＝18.3/5＝3.66mmmm。水泥环层的壁厚为h_f＝76.2mm，沿径向水泥环单元数n_f＝5，泥环层的单元厚度t_f＝h_f/n_f＝76.2/5＝15.24mm

步骤202：设定子结构内径向单元尺寸的增长系数k，设定环向单元划分数m。本实施例中，子结构内单元尺寸的增长系数取为k＝1.45，环向单元划分数取为m＝40。

步骤203：采用下式计算子围岩径向划分的层数N：

本实施例中，母单元的特征长度为L＝2.05m，水泥环的外层单元厚度t_f＝15.24mm，代入上式计算可得围岩部分径向单元层数N＝10。

步骤204：计算自水泥环向外第i层单元的径向尺寸t_i和环向尺寸r_i:

本实施例中，采用前述各数据，计算得到的各层单元的径向尺寸t_i和环向尺寸r_i如表2所示。

层号	径向尺寸(mm)	环向尺寸(mm)
			0	15.24	2.38
1	22.25	5.88
			2	32.48	10.99
3	47.42	18.44
			4	69.24	29.31
5	101.09	45.19
			6	147.6	68.38
7	215.5	102.23
			8	314.63	151.65
9	459.36	223.81
			10	670.67	329.16

采用计算的得到径向尺寸和环向尺寸，构建子结构的平面形式如图5所示，划分后的子结构的三维结构如图6所示。

步骤3：建立子结构系统的刚度矩阵。具体实施如下：

步骤301：建立子结构系统内各单元的刚度矩阵。

为建立子结构系统的单元刚度矩阵，首先需要建立各子单元的插值形函数，对于单元中各节点，采用局部坐标，则各节点插值形函数可表示为：

式中，ξ_i,η_i和ζ_i分别为第i个节点的3个局部坐标。

对各插值形函数求导，并组成偏导数矩阵B可得：

式中：N_i(i＝1,2,...,8)分别为8个节点插值形函数，计算如式(4)所示。

对于每个单元，根据其对应的弹性模型E和泊松比μ，计算其弹性矩阵[D]为：

式中：E和μ分别为对应的套管、水泥环和岩石材料的弹性模型和泊松比。

进而，可构建立子结构内部各单元的刚度矩阵，其计算公式为：

式中：[B]为形函数偏导数矩阵，计算如式(5)所示；[D]为弹性矩阵，计算如式(6)所示。

对于上式中的积分，可以采用高斯积分算法进行。

步骤302：根据子结构内各单元的连接关系，将各单元的刚度矩阵

组装形成子结构的总体刚度矩阵[H]。

采用单元叠加的形式，将单元刚度矩阵[h_i]_24×24叠加到子结构总体刚度矩阵中，本实例中，子结构共有1600个单元，800个节点，因此总体刚度矩阵[H]的维数为2400×2400。可采用自由度映射的方法将各单元的刚度矩阵叠加到总体刚度矩阵[H]中，其具体实施方法为：对于单元刚度矩阵

其局部自由度i对应的总体自由度为I，局部自由度j对应的总体自由度为J。则将元素

累加到总体刚度矩阵的[H_IJ]位置。

步骤303：将分析子结构内节点的连接关系，将节点划分为与母单元共享的出口节点、与母单元不直接相连的内部节点以及连接共享节点与内部节点的过渡节点3部分。本实施例中，出口节点、连接节点和过渡节点的的划分如图4所示。

步骤304：将采用高斯-约当消去法，对总体刚度矩阵[H]进行矩阵运算，得到刚度矩阵[H]分块表达式：

其中，第1列代表的子块为出口刚度矩阵，用于子结构和母体单元的连接，第2列代表的子块为连接刚度矩阵，用子结构和套管内壁单元的连接；第3列代表的子块为内部刚度矩阵，代表套管内壁部分。

步骤4：提取粗网中数值模拟的结果，采用形函数插值法对子结构精细网格内的各单元进行初始化。

步骤401：对子结构中的单元E_i，计算出其第j个积分点的坐标J(i,j)。

本实施例中，子结构中各单元均为8节点6面体单元，每个方向取2个积分点，8个积分点的局部坐标系的坐标如表3所示。

表3各积分点的坐标

积分点编号	ξ坐标	ζ坐标	ψ坐标
				J1	-0.577	-0.577	-0.577
J2	+0.577	-0.577	-0.577
				J3	+0.577	+0.577	-0.577
J4	-0.577	+0.577	-0.577
				J5	-0.577	-0.577	+0.577
J6	+0.577	-0.577	+0.577
				J7	+0.577	+0.577	+0.577
J8	-0.577	+0.577	+0.577

对于每个积分点，根据该单元8个节点的坐标的，采用如下公式计算其在全局坐标系的坐标。

式中：N_i为各个节点插值形函数，计算如式(4)所示。

步骤402：在母单元中，根据坐标J(i,j)，插值计算该点处的应力σ(i,j)，压强p(i,j)，温度T(i,j)。

式中：σ_ij,p_ij,T_ij分别为待求的点处的应力、压强和温度，

为插值形函数，

T_i ^j为各个节点处的应力、压强和温度。

步骤403：将插值得到的应力σ(i,j)，压强p(i,j)，温度T(i,j)赋给单元E_i的第j个积分点

步骤5：对子结构开展有限元计算，获取精细网格内的应力、变形等单元数据或结点位移数据。

步骤501：结合步骤304中的刚度矩阵[H]，形成有限元方程组，方程组可简写为：

式中，H为经过高斯-约当消去后形成的分块刚度矩阵，列向量φ为待求的位移组成的向量，列向量F为等效节点力组成的向量。

步骤502：求解有限元方程组，得到各节点的位移。计算中，可以采用QR分解、LU分解，Cholesky分解等多种方法，求解有限元方程组。

步骤503：计算子结构内各单元的应变和应力。

其中单元应变的计算方法公式为：

{ε}^e＝[B]{φ}^e

式中：{ε}^e为单元应变列向量；[B]为形函数偏导数矩阵，计算如式(5)所示；{φ}^e为单元各节点列向量。

单元应力的计算公式为：

{σ}^e＝[D]{ε}^e

式中：{σ}^e为单元应力列向量，{ε}^e为单元应变列向量；[D]为弹性矩阵计算如式(6)所示。

本实施例，取套管弹性模量为210GPa，水泥环弹性模量为3.8GPa,岩石弹性模量为18GPa，套管内压为27MPa，母单元内的水平大主应力为85.7MPa，水平小主应力为60MPa。采用本发明的子结构方法，计算了套管和水泥环的应力，表4以环向应力为例，给出应力的计算结果：

表4套管外层应力计算结果

根据本实施例计算得出的单元应力和应变，绘制了子结构的应力云图如附图6所示。精细子结构能够反应地层中套管、水泥环等细微结构的应力和变形，套管和水泥环处的应力云图如附图7所示。

Claims (3)

1.一种基于子结构的油气田井筒结构应力精细模拟方法，其特征是：包括以下步骤：

步骤1：在油藏数值模拟的粗网格中，提取油气井筒所在的单元作为母体单元；

步骤2：对母体单元进行精细划分，生成由套管、水泥环逐步过渡到母单元的精细网格，形成含套管-水泥环系统的子结构；

步骤3：建立子结构系统的单元刚度矩阵和荷载列向量；

步骤4：提取粗网中数值模拟的结果，采用形函数插值法对子结构精细网格内的各单元进行初始化；

步骤5：对子结构开展有限元计算，获取精细网格内的应力、变形及结点位移数据；

所述的步骤1包含以下步骤：

步骤101：根据油气井的中心坐标(x_w,y_w),在粗网格中查找平面坐标包含点(x_w,y_w)的单元作为母单元；

步骤102：提取母单元的8个节点坐标(x_i,y_i,z_i)i＝1,...,8；

步骤103：计算母单元各个边的边长l_ii＝1,...,8，取单元的特征长度

所述的步骤2包含以下步骤：

步骤201：设置套管的壁厚h_c和套管层的单元数n_c，计算套管层的单元厚度t_c＝h_c/n_c,设置水泥环的壁厚h_f和水泥环层t_f的单元数n_f，计算水泥环层的单元厚度t_f＝h_f/n_f；

步骤202：设定子结构内径向单元尺寸的增长系数k，设定环向单元划分数m；

步骤203：计算子结构径向划分的层数N：

式中：L为母单元的特征长度；

步骤204：计算自套管向外第i层单元的径向尺寸t_i和环向尺寸r_i:

步骤205：根据各层单元的径向尺寸t_i和环向尺寸r_i形成子结构平面各层的精细网格；

子结构内的单元为自套管向外的辐射状网格，内层网格较为精细可以反应套管、水泥环系统；外层网格较大，可以和粗网格衔接；

所述的步骤3包含以下步骤：

步骤301：计算子结构内各单元的刚度矩阵[h_i]；

步骤302：根据子结构内各单元的连接关系，将各单元的刚度矩阵[h_i]组装，集成为子结构的总体刚度矩阵[H]；

步骤303：分析子结构内节点的连接关系，将节点划分为与母单元共享节点、与母单元不直接相连的内部节点以及连接共享节点与内部节点的过渡节点3部分；

步骤304：采用高斯-约当消去法，对总体刚度矩阵[H]进行矩阵运算，得到刚度矩阵[H]分块表达式：

式中第1列的子块与共享节点相对应，称为出口刚度矩阵；第2列与内部节点相对应，称为内部刚度矩阵；第3列与连接节点相对应，称为连接刚度矩阵。

2.根据权利要求1所述的基于子结构的油气田井筒结构应力精细模拟方法，其特征是：步骤4包含以下步骤：

步骤401：对子结构中的单元E_i，计算出其第j个积分点的坐标J(i,j)；

步骤402：在母单元中，根据坐标J(i,j)，插值计算该点处的应力σ(i,j)，压强p(i,j)，温度T(i,j)；

步骤403：将插值得到的应力σ(i,j)，压强p(i,j)，温度T(i,j)赋给单元E_i的第j个积分点。

3.根据权利要求1所述的基于子结构的油气田井筒结构应力精细模拟方法，其特征是：步骤5包含以下步骤：

步骤501：结合步骤304中的刚度矩阵[H]，分析单元荷载，形成有限元方程组：

[H]{φ}＝{R}

其中{φ}为待求的位移组成的向量；

步骤502：求解有限元方程组，得到各节点的位移；

步骤503：计算子结构内各单元的应力。

Images