CN104318032A - 流固耦合作用下的油田套损计算方法 - Google Patents

Abstract

本发明公开了一种流固耦合作用下的油田套损计算方法，该方法基于流固耦合理论，通过室内试验、数值计算和理论分析相结合的方法开展研究，根据流体力学理论及弹塑性理论，建立了油藏开发过程渗流数学模型及应力场计算数学模型，根据室内试验分析建立了储层岩石的弹塑性本构模型，借助有效应力原理，建立了开发过程流固耦合数学模型，根据流固耦合数学模型及求解思路，在已有商业软件基础上建立了储层流固耦合地应力预测及套损防控的数值模拟方法，以曲堤油田典型区块为例，给出了流固耦合作用下油田套损预测及防控的具体方法，并提出具体措施，对于改善油田套损状况，提高油井使用寿命具有重要的应用价值。

Description

流固耦合作用下的油田套损计算方法

技术领域

本发明属于油田设备技术领域，涉及一种流固耦合作用下的油田套损计算方法。

背景技术

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

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

油田现场大量应用的套管变形检测方法是铅膜法，利用专用管柱或钢丝绳下接铅膜，对套损几何形状等进行打印，然后对打印出来的印痕进行描绘、分析，提出套损点的形状、尺寸、深度位置。铅膜法的优点是迅速、方便和直观，但是铅膜直径大小的选择比较困难。直径过大，铅膜打印出来不在变形最明显处，印痕不清晰，直径过小，印痕不明显。

随着科学技术的进步，国内外专家学者研发和改进了许多工程测井仪器和测井技术从多方面检测套损状况。工程测井技术大多是通过检测到的物理信号来间接判断管内的缺陷情况，井径系列是油水井井身状况常规的检测手段，可以提供套管内径变化情况；声波测井系列中井壁超声成像测井可以提供全面、直观的套损状况，噪声测井则用于检测已形成的管漏、窜槽；方位系列用于确定套管变形的方位角度；磁测井系列可检查套管变形、错断、内外壁腐蚀及射孔质量；作为辅助的井温和注产系列可用于评价套管漏失和层间窜槽情况。

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

发明内容

本发明的目的在于提供流固耦合作用下的油田套损计算方法，解决了目前物理模型实验，很难模拟复杂地层套管受力和变形过程的问题。

本发明所采用的技术方案是按照以下步骤进行：

步骤1：根据流体力学理论及油藏开发涉及的运动状态，采用理论分析方法建立渗流场数学模型，渗流场数学模型包括单相渗流数学模型、油水两相渗流数学模型和渗流场边界条件；

步骤2：根据弹塑性力学理论，在对储层岩石岩心数据进行详细分类及分析的基础上，建立储层开发过程的应力场计算数学模型，以岩心的强度测试数据为基础，采用DP准则来描述储层岩石的弹塑性本构关系；

步骤3：以有效应力原理为桥梁，根据孔隙介质变形理论，建立储层岩石的孔隙度与孔隙压力的关系及渗透率与孔压、有效应力的关系，结合室内岩心应力测试结果，建立储层注水开发过程流固耦合数学模型；

步骤4：建立相应区块的地质模型：采用三维地质建模软件PETREL建立该模拟区块的地质模型，根据应力场计算的每个单元的节点坐标数据和和流场计算的块中心网格数据的几何对应关系，将PETREL建立的地质网格模型转化为应力场计算所需的有限元数值模型。

步骤5：将步骤4建立的地质模型导入eclipse，在eclipse软件中设置好流场边界及初始条件，根据步骤1中的流场单相流、两相流模型、边界类型及初始条件在eclipse中选择对应的流场模拟计算模型；结合每口井的生产情况，输入井的生产数据，在eclipse中设置好模型的边界及初始条件，计算模拟区域的孔隙压力；将步骤4转化出的有限元模型导入abaqus，将eclipse计算出的孔压导入abaqus作为初始应力条件，基于步骤2的应力场数学模型在abaqus中选择对应的应力计算模块，根据步骤3的渗透率随有效应力关系设计abaqus子程序来实现耦合计算过程中的有效应力及孔压对渗透率的影响计算，最终计算出储层的地应力及套管单元的受力，轴向载荷等大小，根据套管管材的强度判断套管的损坏情况，结合油田现场的油井生产制度给出典型生产工况下的套损预防措施。

本发明的有益效果是提供了流固耦合作用下的油田套损预测数学模型并通过仿真表现出来，能准确进行油田套损的防治。

说明书附图

图1是本发明实施例提供的控制体的质量守恒示意图；

图2是本发明实施例提供的深度坐标D示意图；

图3是本发明实施例提供的增量法示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图及具体实施例对本发明的应用原理作进一步描述。

本发明按照以下步骤进行：

步骤一：根据流体力学理论及油藏开发涉及的运动状态，采用理论分析方法建立渗流场数学模型，渗流场数学模型包括单相渗流数学模型、油水两相渗流数学模型和渗流场边界条件；

单相渗流数学模型：

通过研究直角坐标系中的一个小控制体，得到变形介质单相流体渗流的数学方程。假定围绕多孔介质区域的点P(x,y,z)有一个尺寸为δx、δy、δz(其边分别为平行于坐标轴x、y、z)的控制体。令表示密度为ρ的流体质量通量，即单位时间通过单位面积的质量。如图1所示控制体的质量守恒，设在x、y、z方向的分量为J_x,J_y,J_z。在短时间段δt内，通过该控制体的两个垂直于x方向的表面，流入超出流出的量可以用它们的差值表示：

[J_x|_{x-(δx/2),y,z}-J_x|_{x+(δx/2),y,z}]δyδzδt

近似地在P点周围按Taylor级数展开，并忽略二阶和高阶项，得：

$-({\∂J}_x/\∂x)\mathrm{\δx\δy\δz\δt}$

对另外两个方向重复上述过程，然后相加，可得到该控制体表面流入超出流出的数量为：

$-(\frac{{\∂J}_x}{\∂x}+\frac{{\∂J}_y}{\∂y}+\frac{{\∂J}_z}{\∂z})\mathrm{\δx\δy\δz\δt}$

按照质量守恒原理，该量必等于δt期间此控制体内质量的变化，即其中ΔV₀＝δxδyδz＝常数，为控制单元的体积。因此：

$-(\frac{{\∂J}_x}{\∂x}+\frac{{\∂J}_y}{\∂y}+\frac{{\∂J}_z}{\∂z})=\frac{\∂\left(\mathrm{\ρ\φ}\right)}{\∂t}---\left(1\right)$

$\mathrm{div}\stackrel{\→}{J}+\frac{\∂\left(\mathrm{\ρ\φ}\right)}{\∂t}=0$

质量通量可表示为：

$\stackrel{\→}{J}=\mathrm{\ρ}\stackrel{\→}{V}---\left(2\right)$

式中，ρ和表示平均值。

将式(2)代入到式(1)中，可得

$\frac{\∂\left(\mathrm{\ρ}{V}_x\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρ}{V}_y\right)}{\∂y}+\frac{\∂\left(\mathrm{\ρ}{V}_z\right)}{\∂z}+\frac{\∂\left(\mathrm{\ρ\φ}\right)}{\∂t}=0---\left(3\right)$

对于不变形的介质，φ＝常数，于是式(3)变为：

当考虑介质的变形时，方程(3)可写为：

$\frac{\∂\left(\mathrm{\ρ}{V}_x\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρ}{V}_y\right)}{\∂y}+\frac{\∂\left(\mathrm{\ρ}{V}_z\right)}{\∂z}+\mathrm{\φ}\frac{\∂\mathrm{\ρ}}{\∂t}+\mathrm{\ρ}\frac{\∂\mathrm{\ρ}}{\∂t}=0---\left(5\right)$

由于多孔介质和其中的流体是可压缩的，在等温情况下，油藏流体的密度ρ是压力p的函数，即：

式中，C_f为流体的压缩系数，

故 $\mathrm{\φ}\frac{\∂\mathrm{\ρ}}{\∂t}=\mathrm{\φ}{\mathrm{\ρ}}_0{C}_f\frac{\∂p}{\∂t}---\left(7\right)$

根据达西定律，流体在多孔介质中的渗流流速为：

$v_f=\frac{K}{\mathrm{\μ}}\mathrm{grad}p---\left(8\right)$

将式(8)代入到式(5)中，得：

$\frac{\∂}{\∂x}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂x}\right)+\frac{\∂}{\∂y}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂y}\right)+\frac{\∂}{\∂z}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂z}\right)+\mathrm{\φ}{\mathrm{\ρ}}_0{C}_f\frac{\∂p}{\∂t}+\mathrm{\ρ}\frac{\∂\mathrm{\φ}}{\∂t}=0---\left(9\right)$

上式简写为

$\▿\·(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\▿\·p)+\mathrm{\φ}{\mathrm{\ρ}}_0{C}_f\frac{\∂p}{\∂t}+\mathrm{\ρ}\frac{\∂\mathrm{\φ}}{\∂t}=0---\left(10\right)$

如果介质完全被流体所饱和，在流动期间，由于外部载荷的变化和(或)孔隙压力的变化，导致储层有效应力变化，使多孔介质发生变形。假定岩体只发生孔隙变形，颗粒是不可压缩的。则孔隙度与体积应变的关系为：

$\mathrm{\φ}=\frac{{\mathrm{\φ}}_0+{\mathrm{\ϵ}}_v}{1+{\mathrm{\ϵ}}_v}$

故 $\frac{\∂\mathrm{\φ}}{\∂t}=\frac{1-{\mathrm{\φ}}_0}{{(1+{\mathrm{\ϵ}}_v)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}---\left(11\right)$

式中，ε_V＝ε_x+ε_y+ε_z为体积应变，φ₀为初始孔隙度。

将式(11)代入到式(9)则得：

$\frac{\∂}{\∂x}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂x}\right)+\frac{\∂}{\∂y}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂y}\right)+\frac{\∂}{\∂z}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂z}\right)+\mathrm{\φ}{\mathrm{\ρ}}_0{C}_f\frac{\∂p}{\∂t}+\frac{\mathrm{\ρ}(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_v)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}=0---\left(12\right)$

如考虑源汇项，假设单位体积注入(采出)的质量流量为q，则变形介质中单相流体渗流的数学模型为：

$\frac{\∂}{\∂x}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂x}\right)+\frac{\∂}{\∂y}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂y}\right)+\frac{\∂}{\∂z}\left(\frac{\mathrm{\ρK}}{\mathrm{\μ}}\frac{\∂p}{\∂z}\right)+\mathrm{\φ}{\mathrm{\ρ}}_0{C}_f\frac{\∂p}{\∂t}+\frac{\mathrm{\ρ}(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_v)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+q=0---\left(13\right)$

油水两相渗流数学模型：

在油藏渗流中，通常存在油、水两相流动。油藏数值模拟的模型主要有黑油模型(Black Oil model)和多组分模型(Mutiple Composite Model)2种。黑油模型是指非挥发性原油的模型，是相对于油质极轻的挥发性油而言的，因其油质较重而色泽较深，故称之为黑油。随着凝析水田的开发以及各种提高采收率方法的试验和应用，油藏内不仅有多相流体渗流，而且还有相间传质现象(相变)，这是多组分模型研究的内容。本发明根据简化的黑油模型理论，将变形介质单相渗流的数学模型，拓展为油、水两相流体渗流的情况。

油、水两相渗流模型中：

(1)油藏中的渗流是等温渗流；

(2)油藏烃类只含油组分，且油组分是指将地层原油在地面标准状况下经分离后所残存的液体。储层中只有油、水两相参加渗流，每一相的渗流均符合达西定律；

(3)油、水之间不互溶。

在均质和各向同性的油藏中，油、水两相的运动方程分别为：

$v_o=-\frac{K{K}_{\mathrm{ro}}}{{\mathrm{\μ}}_o}(\▿p_o-{\mathrm{\ρ}}_{o}g\▿D)---\left(14\right)$

$v_w=-\frac{K{K}_{\mathrm{rw}}}{{\mathrm{\μ}}_w}(\▿p_w-{\mathrm{\ρ}}_{w}g\▿D)---\left(15\right)$

上两式中的D为深度坐标，如图2所示。

式中K为油藏岩石的绝对渗透率；K_ro、K_rw分别为油相、水相的相对渗透率；μ_o、μ_w分别为油、水的动力粘滞系数；ρ_o、ρ_w分别为油、水的密度；p_o、p_w分别为油相、水相的孔隙压力。

根据方程(13)可直接给出两相流体的渗流方程为：

油相：

$\▿\·\left[\frac{{\mathrm{\ρ}}_o{\mathrm{KK}}_{\mathrm{ro}}}{{\mathrm{\μ}}_o}(\▿p_o-{\mathrm{\ρ}}_{o}g\▿D)\right]+\mathrm{\φ}\frac{\∂\left({\mathrm{\ρ}}_o{S}_o\right)}{\∂t}+\frac{{\mathrm{\ρ}}_0{S}_o(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+q_o=0---\left(16a\right)$

水相：

$\▿\·\left[\frac{{\mathrm{\ρ}}_w{\mathrm{KK}}_{\mathrm{rw}}}{{\mathrm{\μ}}_w}(\▿p_w-{\mathrm{\ρ}}_{w}g\▿D)\right]+\mathrm{\φ}\frac{\∂\left({\mathrm{\ρ}}_w{S}_w\right)}{\∂t}+\frac{{\mathrm{\ρ}}_w{S}_w(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+q_w=0---\left(16b\right)$

式中q_o、q_w分别为油、水产量(产出或注入项，即源或汇项)；φ为储层岩石孔隙度；S_o、S_w分别为油、水相饱和度。

求解时，除了2个基本方程为，还需要辅助方程，它们是饱和度约束方程

S_o+S_w＝1    (17)

毛管压力是流体饱和度的函数，与各相压力之间的关系：

p_cow＝p_cow(S_w)＝p_o-p_w    (18)

状态方程(包括密度和粘度，即高压物性参数方程)：

ρ_o＝ρ_o(p_o)    (19)

ρ_w＝ρ_w(p_w)    (20)

μ_o＝μ_o(p_o)    (21)

μ_w＝μ_w(p_w)    (22)

相对渗透率与流体饱和度的关系：

K_ro＝K_ro(S_o)    (23)

K_rw＝K_rw(S_w)    (24)

式(19)、(20)中油、水密度ρ与压力p的函数关系，在不考虑温度影响的条件下，根据式(7)得：

$\frac{\∂{\mathrm{\ρ}}_o}{\∂t}={\mathrm{\ρ}}_o^0{C}_o\frac{\∂p_o}{\∂t}---\left(25a\right)$

$\frac{\∂{\mathrm{\ρ}}_w}{\∂t}={\mathrm{\ρ}}_w^0{C}_w\frac{\∂p_w}{\∂t}---\left(25b\right)$

又由于在式(16)中有：

$\mathrm{\φ}\frac{\∂\left({\mathrm{\ρ}}_o{S}_o\right)}{\∂t}=\mathrm{\φ}{S}_o\frac{\∂{\mathrm{\ρ}}_o}{\∂t}+\mathrm{\φ}{\mathrm{\ρ}}_o\frac{\∂S_o}{\∂t}---\left(26a\right)$

$\mathrm{\φ}\frac{\∂\left({\mathrm{\ρ}}_w{S}_w\right)}{\∂t}=\mathrm{\φ}{S}_w\frac{\∂{\mathrm{\ρ}}_w}{\∂t}+\mathrm{\φ}{\mathrm{\ρ}}_w\frac{\∂S_w}{\∂t}---\left(26b\right)$

即可得到考虑介质变形时的简化的油水两相黑油模型的基本微分方程：油组分方程：

$\▿\·\left[\frac{{\mathrm{KK}}_{\mathrm{ro}}}{{\mathrm{\μ}}_o}(\▿p_o-{\mathrm{\ρ}}_{o}g\▿D)\right]+\mathrm{\φ}{C}_o{S}_o\frac{{\mathrm{\ρ}}_o^0}{{\mathrm{\ρ}}_o}\frac{\∂p_o}{\∂t}+\mathrm{\φ}\frac{\∂S_o}{\∂t}+\frac{S_o(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+\frac{q_o}{{\mathrm{\ρ}}_o}=0---\left(27a\right)$

水组分方程：

$\▿\·\left[\frac{{\mathrm{KK}}_{\mathrm{rw}}}{{\mathrm{\μ}}_w}(\▿p_w-{\mathrm{\ρ}}_{w}g\▿D)\right]+\mathrm{\φ}{C}_w{S}_w\frac{{\mathrm{\ρ}}_w^0}{{\mathrm{\ρ}}_w}\frac{\∂p_w}{\∂t}+\mathrm{\φ}\frac{\∂S_w}{\∂t}+\frac{S_w(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+\frac{q_w}{{\mathrm{\ρ}}_w}=0---\left(27b\right)$

对于一般的非低渗透油田，油、水两相流动中，忽略毛管压力的作用，油相与水相压力相等，即：p₀＝p_w＝p，这样式(27a，27b)变为

油组分方程：

$\▿\·\left[\frac{{\mathrm{KK}}_{\mathrm{ro}}}{{\mathrm{\μ}}_o}(\▿p-{\mathrm{\ρ}}_{o}g\▿D)\right]+\mathrm{\φ}{C}_o{S}_o\frac{{\mathrm{\ρ}}_o^0}{{\mathrm{\ρ}}_o}\frac{\∂p}{\∂t}+\mathrm{\φ}\frac{\∂S_o}{\∂t}+\frac{S_o(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+\frac{q_o}{{\mathrm{\ρ}}_o}=0---\left(28a\right)$

水组分方程：

$\▿\·\left[\frac{{\mathrm{KK}}_{\mathrm{rw}}}{{\mathrm{\μ}}_w}(\▿p-{\mathrm{\ρ}}_{w}g\▿D)\right]+\mathrm{\φ}{C}_w{S}_w\frac{{\mathrm{\ρ}}_w^0}{{\mathrm{\ρ}}_w}\frac{\∂p}{\∂t}+\mathrm{\φ}\frac{\∂S_w}{\∂t}+\frac{S_w(1-{\mathrm{\φ}}_0)}{{(1+{\mathrm{\ϵ}}_V)}^2}\frac{\∂{\mathrm{\ϵ}}_V}{\∂t}+\frac{q_w}{{\mathrm{\ρ}}_w}=0---\left(28b\right)$

渗流场边界条件：对于实际问题，必须补充相应的定解条件，才能求解具体问题。定解条件包括初始条件和边界条件。

(1)边界条件：

油藏模拟中的边界条件主要包括两类：定压边界和定流量边界。所谓定压边界是指边界上每一点在每个时刻的压力是已知的或井的井底压力已知。这类边界条件可表示为：

p_G＝f₁(x,y,z,t)    (29)

上式表示边界上点(x、y、z)在t时刻的压力p_G，为给定的函数f₁(x,y,z)。

如在边界上流量已知，则这类边界条件为给定流量的边界，即第二类边界条件，表示为：

${\frac{\∂p}{\∂n}|}_G=f_2(x,y,z,t)---\left(30\right)$

式中：表示压力关于边界的外法线方向导数，f₂(x,y,z,t)为指定边界上的已知函数。

(2)初始条件：

求解非稳定的渗流问题，除了需要边界条件外，还需要初始条件。即初始时刻(t＝0)油层中的压力和饱和度分布。表示为：

p(x,y,z,0)＝p₀(x,y,z)    (31)

S(x,y,z,0)＝S₀(x,y,z)    (32)

步骤二：根据弹塑性力学理论，在对储层岩石岩心数据进行详细分类及分析的基础上，建立储层开发过程的应力场计算数学模型，以岩心的强度测试数据为基础，采用DP准则来描述储层岩石的弹塑性本构关系；

应力场数学模型：

多孔介质的应力场基本方程包括：平衡微分方程、几何方程。

平衡微分方程为：

$\frac{\∂{\mathrm{\σ}}_{\mathrm{ij}}}{\∂x_j}-f_{x_i}=0,i,j=\mathrm{1,2,3}---\left(33\right)$

式中，x_j为x，y，z三个坐标方向，为x_i方向的体积力。

考虑变形的几何条件，弹性体发生变形时，通过分析微元体的形变分量和位移分量之间的关系，建立起几何方程为：

${\mathrm{\ϵ}}_{\mathrm{ij}}=\frac{1}{2}(u_{i,j}+u_{j,i})---\left(34\right)$

对于弹性问题，由广义胡克定律，应力应变之间满足本构方程：

σ_ij＝D_ijklε_kl    (35)

式中，下标i，j，k，l表示坐标，u_i,j表示i方向的位移对坐标j的偏导，u_j,i表示j方向位移对i坐标偏导，ε_kl表示垂直于k坐标的平面上的l向应变。

弹塑性本构方程：

研究表明，岩土介质常常表现出非线性弹性及弹塑性力学行为，单纯的线弹性理论并不能完全反映耦合场中的介质的变形特点；用于求解非线性弹性及弹塑性问题的数值方法，通常采用图3所示的分段线性化处理方法，在每一增量步内，将总应变增量表示为弹性应变增量和塑性应变增量之和，也即：

{dε}＝{dε^e}+{dε^p}    (36)

式中：{dε}、{dε^e}和{dε^p}分别为总应变、弹性应变和塑性应变，其对应的应力增量分别为{dσ}、{dσ^e}和{dσ^p}。

弹性应变增量与应力增量之间的关系可由弹性模型确定，塑性应变增量与应力增量之间的关系则需要由相应的弹塑性本构方程确定。岩体的塑性本构理论包括屈服准则，加、卸准则及增量形式的应力-应变关系3个方面：

1)屈服准则

在复杂应力状态下，岩体内一点出现塑性变形时应力所应满足的条件称为屈服条件，屈服条件可用应力分量来表示，即可写成应力分量的函数：

F(σ_x,σ_y,σ_z,τ_xy,τ_yz,τ_zx)＝0    (37)

式(37)中6个应力分量与所选取的坐标系相关。采用不同的坐标系，就得到不同数值的应力分量，不便使用,为此，通常采用与坐标轴方向无关的量来表示，如主应力表示：

F(σ₁,σ₂,σ₃)＝0    (38)

或应力不变量表示：

F(I₁,J₂,,J₃)＝0    (39)

对于岩土，采用的屈服准则Mohr-Coulomb准则和Drucker-Prager准则，对于Mohr-Coulomb准则，当应力状态达到下述条件时，材料进入屈服状态：

式中：τ为最大剪应力，σ_n为同一平面内的正应力；C为岩体的介质的粘聚力，为介质的内摩擦角。

Mohr-Coulomb准则常写成用主应力表达的形式：

考虑产生屈服的各种可能的应力组合，Mohr-Coulomb准则的屈服面为一个角锥面；角锥的顶点在静水应力轴上。

由于Mohr-Coulomb准则的屈服面为角锥面，其角点在数值计算中常引起不便，为得到近似于Mohr-Coulomb曲面的光滑屈服面，Drucker－Prager采用如下公式代替：

$F=\mathrm{\α}{I}_1+\sqrt{J_2}-k=0---\left(42\right)$

式中：I₁为应力第一不变量，且有I₁＝σ₁+σ₂+σ₃；J₂为应力偏量第二不变量，且有： $J_2=\frac{1}{6}[{({\mathrm{\σ}}_1-{\mathrm{\σ}}_2)}^2+{({\mathrm{\σ}}_2-{\mathrm{\σ}}_3)}^2+{({\mathrm{\σ}}_3-{\mathrm{\σ}}_1)}^2].$

适当的选取常数α和k可以使Drucker-Prager屈服面接近于Mohr-Coulomb屈服面，例如，取：

则各截面上，Drucker－Prager屈服圆都与Mohr-Coulomb六边形的外顶点重合。如取：

则各截面上，Drucker－Prager屈服圆都与Mohr-Coulomb六边形的内顶点重合。

2)加、卸准则

当岩体介质达到屈服状态后，加载和卸载时的应力应变规律不同，对于理想弹塑性岩体，材料不发生强化，加载条件和卸载条件相同，其加卸载准则可表示为：

F(σ_ij)＜0    (弹性状态)

$F\left({\mathrm{\σ}}_{\mathrm{ij}}\right)=0,\mathrm{dF}=\frac{\∂F}{\∂{\mathrm{\σ}}_{\mathrm{ij}}}d{\mathrm{\σ}}_{\mathrm{ij}}=0$ (加载)    (45)

$F\left({\mathrm{\&sigma;}}_{\mathrm{ij}}\right)=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}<0$ (卸载)

对于强化材料，其加卸载准则可表示为：

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}>0$   (加载)

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}=0$ (中性加载)    (46)

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}<0$   (卸载)

增量形式的应力-应变关系

由塑性力学可知，塑性状态下，应力增量与应变增量的关系式为：

{dσ}＝([D]-[D]_p){dε}＝[D]_ep{dε}    (47)

式中：

[D]_ep＝[D]-[D]_p    (48)

$\left[D_p\right]=\frac{\left[D\right]\left\{\frac{\&PartialD;G}{\&PartialD;\mathrm{\&sigma;}}\right\}{\left\{\frac{\&PartialD;F}{\&PartialD;\mathrm{\&sigma;}}\right\}}^T\left[D\right]}{A+{\left\{\frac{\&PartialD;F}{\&PartialD;\mathrm{\&sigma;}}\right\}}^T\left[D\right]\left\{\frac{\&PartialD;G}{\&PartialD;\mathrm{\&sigma;}}\right\}}---\left(49\right)$

对于理想弹塑性模型，上式中A＝0；在Drucker公设成立的条件下，塑性位势能函数Q＝F，此时[D]_ep为对称矩阵。

式中，[D]表示弹性矩阵，[D]_p表示塑性矩阵：

此时，对于Mohr-Coulomb准则，有：

式中： $s_i=D_{i1}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x+D_{i2}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_y+D_{i3}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_z,(i=\mathrm{1,2,3})$

$s_i=G{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{kj}},(\mathrm{kj}=\mathrm{xy},\mathrm{yz},\mathrm{zx})$

$s_0=s_1{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x+s_2{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_y+s_3{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_z+s_4{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{xy}}+s_5{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{yz}}+s_6{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{zx}}$

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_x}=\mathrm{\&alpha;}+({\mathrm{\&sigma;}}_x-\mathrm{\&sigma;})/2\sqrt{J_2},(\mathrm{xy},\mathrm{yz},\mathrm{zx})$

步骤三：以有效应力原理为桥梁，根据孔隙介质变形理论，建立储层岩石的孔隙度与孔隙压力的关系及渗透率与孔压、有效应力的关系，结合室内岩心应力测试结果，建立储层注水开发过程流固耦合数学模型；

流固耦合数学模型：

根据Biot有效应力原理，作用在饱和岩土体上的外力由岩土介质内的骨架和孔隙水共同承担，外力所引起的岩土体内的总应力由介质骨架内的有效应力和孔隙水压力两部分组成，对于各向同性的岩土体，孔隙流体压力只能使介质发生体积改变而不能使介质产生形状改变，因此，介质的剪应力与孔隙压力无关，孔隙压力对介质正应力的影响在各个方向上相同，即：

σ_ij＝σ′_ij+αpδ_ij    (51)

式中：σ_ij为介质内总应力，σ′_ij为介质内的有效应力，p为孔隙水压力，δ_ij为Kroneker符号，α为Biot常数。

将有效应力公式(51)代入式(33)可得用有效应力和孔隙压力表示的平衡方程：

$\frac{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}^{\&prime;}}{\&PartialD;x_j}+\frac{\&PartialD;p}{\&PartialD;x}-f_{x_i}=0---\left(52\right)$

将岩土体变形场中的本构方程(37)，几何方程(36)代入由有效应力表示的平衡方程(52)可得用位移表示的平衡方程：

$G{\&dtri;}^2{u}_i-(\mathrm{\&lambda;}+G)\frac{\&PartialD;{\mathrm{\&epsiv;}}_v}{\&PartialD;x_i}-\frac{\&PartialD;p}{\&PartialD;x_i}+f_{x_i}=0,i=\mathrm{1,2,3}---\left(53\right)$

式中，λ、G为拉梅常数，ε_v为体积变形，且有其他各符号含义同前。

应力场的平衡方程式(53)中，p为流体孔隙压力，反映了流体渗流场对应力场的影响。

介质孔隙度与孔隙压力的关系

岩土体的孔隙度φ随孔隙压力p的变化可用压缩系数C_R来表征，其定义为

$C_R=\frac{1}{\mathrm{\&phi;}}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;p}---\left(54\right)$

由上式可有

$C_R\&PartialD;p=\frac{\&PartialD;\mathrm{\&phi;}}{\mathrm{\&phi;}}\&DoubleRightArrow;{\&Integral;}_{p_0}^p{C}_R\mathrm{dp}={\&Integral;}_{{\mathrm{\&phi;}}_0}^{\mathrm{\&phi;}}\frac{\mathrm{d\&phi;}}{\mathrm{\&phi;}}\&DoubleRightArrow;C_R(p-p_0)=\ln \frac{\mathrm{\&phi;}}{{\mathrm{\&phi;}}_0}$

式中φ₀为岩土初始孔隙度，从而：

$\mathrm{\&phi;}={\mathrm{\&phi;}}_0{e}^{C_R(p-p_0)}---\left(55\right)$

在地层压力变化不大高时，上式可按Taylor级数展开并取展开式的前两项

研究表明，岩土介质常常表现出非线性弹性及弹塑性力学行为，单纯的线弹性理论并不能完全反映耦合场中的介质的变形特点；用于求解非线性弹性及弹塑性问题的数值方法，通常采用图3所示的分段线性化处理方法，在每一增量步内，将总应变增量表示为弹性应变增量和塑性应变增量之和，也即：

{dε}＝{dε^e}+{dε^p}    (56)

式中：{dε}、{dε^e}和{dε^p}分别为总应变、弹性应变和塑性应变，其对应的应力增量分别为{dσ}、{dσ^e}和{dσ^p}。

弹性应变增量与应力增量之间的关系可由弹性模型确定，塑性应变增量与应力增量之间的关系则需要由相应的弹塑性本构方程确定。岩体的塑性本构理论包括屈服准则，加、卸准则及增量形式的应力-应变关系3个方面：

1)屈服准则

在复杂应力状态下，岩体内一点出现塑性变形时应力所应满足的条件称为屈服条件，屈服条件可用应力分量来表示，即可写成应力分量的函数：

F(σ_x,σ_y,σ_z,τ_xy,τ_yz,τ_zx)＝0    (57)

式(37)中6个应力分量与所选取的坐标系相关。采用不同的坐标系，就得到不同数值的应力分量，不便使用,为此，通常采用与坐标轴方向无关的量来表示，如主应力表示：

F(σ₁,σ₂,σ₃)＝0    (58)

或应力不变量表示：

F(I₁,J₂,,J₃)＝0    (59)

对于岩土，采用的屈服准则Mohr-Coulomb准则和Drucker-Prager准则，对于Mohr-Coulomb准则，当应力状态达到下述条件时，材料进入屈服状态：

式中：τ为最大剪应力，σ_n为同一平面内的正应力；c为岩体的介质的粘聚力，为介质的内摩擦角。

Mohr-Coulomb准则常写成用主应力表达的形式：

考虑产生屈服的各种可能的应力组合，Mohr-Coulomb准则的屈服面为一个角锥面；角锥的顶点在静水应力轴上。

由于Mohr-Coulomb准则的屈服面为角锥面，其角点在数值计算中常引起不便，为得到近似于Mohr-Coulomb曲面的光滑屈服面，Drucker－Prager建议采用如下公式代替：

$F=\mathrm{\&alpha;}{I}_1+\sqrt{J_2}-k=0---\left(62\right)$

式中：I₁为应力第一不变量，且有I₁＝σ₁+σ₂+σ₃；J₂为应力偏量第二不变量，且有： $J_2=\frac{1}{6}[{({\mathrm{\&sigma;}}_1-{\mathrm{\&sigma;}}_2)}^2+{({\mathrm{\&sigma;}}_2-{\mathrm{\&sigma;}}_3)}^2+{({\mathrm{\&sigma;}}_3-{\mathrm{\&sigma;}}_1)}^2].$

适当的选取常数α和k可以使Drucker-Prager屈服面接近于Mohr-Coulomb屈服面例如，取：

则各截面上，Drucker－Prager屈服圆都与Mohr-Coulomb六边形的外顶点重合。如取：

则各截面上，Drucker－Prager屈服圆都与Mohr-Coulomb六边形的内顶点重合。

2)加、卸准则：

当岩体介质达到屈服状态后，加载和卸载时的应力应变规律不同，对于理想弹塑性岩体，材料不发生强化，加载条件和卸载条件相同，其加卸载准则可表示为：

F(σ_ij)＜0    (弹性状态)

$F\left({\mathrm{\&sigma;}}_{\mathrm{ij}}\right)=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}=0$ (加载)    (65)

$F\left({\mathrm{\&sigma;}}_{\mathrm{ij}}\right)=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}<0$ (卸载)

对于强化材料，其加卸载准则可表示为：

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}>0$ (加载)

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}=0$ (中性加载)    (66)

$F=0,\mathrm{dF}=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_{\mathrm{ij}}}d{\mathrm{\&sigma;}}_{\mathrm{ij}}<0$    (卸载)

3)增量形式的应力-应变关系：

由塑性力学可知，塑性状态下，应力增量与应变增量的关系式为：

{dσ}＝([D]-[D]_p){dε}＝[D]_ep{dε}    (67)

式中：

[D]_ep＝[D]-[D]_p    (68)

$\left[D_p\right]=\frac{\left[D\right]\left\{\frac{\&PartialD;G}{\&PartialD;\mathrm{\&sigma;}}\right\}{\left\{\frac{\&PartialD;F}{\&PartialD;\mathrm{\&sigma;}}\right\}}^T\left[D\right]}{A+{\left\{\frac{\&PartialD;F}{\&PartialD;\mathrm{\&sigma;}}\right\}}^T\left[D\right]\left\{\frac{\&PartialD;G}{\&PartialD;\mathrm{\&sigma;}}\right\}}---\left(69\right)$

对于理想弹塑性模型，上式中A＝0；在Drucker公设成立的条件下，塑性位势能函数：

Q＝F，此时[D]_ep为对称矩阵。此时，对于Mohr-Coulomb准则，有：

式中： $s_i=D_{i1}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x+D_{i2}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_y+D_{i3}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_z,(i=\mathrm{1,2,3})$

$s_i=G{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{kj}},(\mathrm{kj}=\mathrm{xy},\mathrm{yz},\mathrm{zx})$

$s_0=s_1{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x+s_2{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_y+s_3{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_z+s_4{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{xy}}+s_5{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{yz}}+s_6{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{zx}}$

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_x=\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&sigma;}}_x}=\mathrm{\&alpha;}+({\mathrm{\&sigma;}}_x-\mathrm{\&sigma;})/2\sqrt{J_2}$

式(55)即为岩土体介质孔隙度与孔隙压力的耦合关系。

介质渗透率与孔隙压力、有效应力的关系

研究表明，油藏流固耦合问题中，储层变形后，随着孔隙度的变化的，其渗透张量不再是一个常数，而是孔隙度的函数，渗透系数张量与孔隙变形之间满足经验公式：

$K_{\mathrm{ij}}\left(\mathrm{\&Delta;\&phi;}\right)=K_{\mathrm{ij}}^0\exp [-\mathrm{\&alpha;\&Delta;\&phi;}]---\left(71\right)$

式中：K为渗透系数张量，Δφ为介质孔隙度的增量，α为耦合系数。

油藏流固耦合问题中，常常通过实验建立储层的渗透张量与孔隙压力和有效应力之间的关系，由于多孔介质中孔隙变化量与孔隙压力及单元体的应力张量之间满足关系式：

$\mathrm{\&Delta;\&phi;}=\frac{p}{R}-\frac{{\mathrm{\&sigma;}}_{\mathrm{ii}}}{3H}---\left(72\right)$

式中：σ为单元体的应力张量，Δφ为介质孔隙度的增量，p为孔隙压力，R、H为反映多孔介质变形性质的常数。

将式(56)代入式(55)可得渗透系数张量与孔隙压力及应力张量之间的关系式：

$K_{\mathrm{ij}}(p,{\mathrm{\&sigma;}}_{\mathrm{ij}})=K_{\mathrm{ij}}^0\exp [-\mathrm{\&alpha;}(\frac{p}{R}-\frac{{\mathrm{\&sigma;}}_{\mathrm{ii}}}{3H})]---\left(73\right)$

式(57)描述流固耦合过程中应力场(变形场)对渗透系数的影响。

步骤四：建立相应区块的地质模型：

采用三维地质建模软件PETREL建立该模拟区块的地质模型，根据应力场计算的每个单元的节点坐标数据和和流场计算的块中心网格数据的几何对应关系，将PETREL建立的地质网格模型转化为应力场计算所需的有限元数值模型。

步骤五：将步骤四建立的地质模型导入eclipse，在eclipse软件中设置好流场边界及初始条件，根据步骤一中的流场单相流、两相流模型、边界类型及初始条件在eclipse中选择对应的流场模拟计算模型；结合每口井的生产情况，输入井的生产数据，在eclipse中设置好模型的边界及初始条件，计算模拟区域的孔隙压力；将步骤四转化出的有限元模型导入abaqus，将eclipse计算出的孔压导入abaqus作为初始应力条件，基于步骤二的应力场数学模型在abaqus中选择对应的应力计算模块，根据步骤三的渗透率随有效应力关系设计abaqus子程序来实现耦合计算过程中的有效应力及孔压对渗透率的影响计算，最终计算出储层的地应力及套管单元的受力，轴向载荷等大小，根据套管管材的强度判断套管的损坏情况，结合油田现场的油井生产制度给出典型生产工况下的套损预防措施。

开发过程套管受力计算

开发过程中，由于注水或者产液而导致近井地层内应力场变化，引起的外部载荷综合作用使得套管内的应力过大造成套管损坏。

(1)热注过程套管受力理论分析

套管在热注生产过程中，受到如下力的作用：

内压力：套管内流体作用在套管内壁上的压强，可由地面套压和流体静压强公式计算。

地层和水泥环自重对套管的外挤力：固井水泥凝固后，水泥环与套管的作用力。

初始轴向力：由套管自重和水泥浆浮力共同作用的结果。形成于注入水泥过程中。

弯曲应力：由于井眼存在一定的曲率，当套管就位后，也要产生与井眼一致的弯曲，在套管内产生弯曲应力。

水泥浆对套管的外挤压力和套管内压力在套管内产生的应力

水泥浆对套管的外挤力:

$p_{\mathrm{co}}=g\underset{0}{\overset{z}{\&Integral;}}{\mathrm{\&rho;}}_c\cos \mathrm{\&alpha;dz}---\left(74\right)$

式中，p_co为水泥浆对套管的外挤力，ρ_c为水泥浆或钻井液密度，g为重力加速度，α为井角。

套管内压力:

$p_{\mathrm{ci}}=p_{\mathrm{cic}}+g\underset{0}{\overset{z}{\&Integral;}}{\mathrm{\&rho;}}_c\cos \mathrm{\&alpha;dz}---\left(75\right)$

式中，p_ci为套管内压力，p_cic为地面套压，ρ_o为套管内流体密度。

在水泥浆的外挤力和内压力的作用下套管的应力

根据厚壁筒理论:

${\mathrm{\&sigma;}}_{r1}=\frac{p_{\mathrm{ci}}{r}_{\mathrm{ci}}^2-p_{\mathrm{co}}{r}_{\mathrm{co}}^2}{r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2}-\frac{(p_{\mathrm{ci}}-p_{\mathrm{co}})r_{\mathrm{co}}^2{r}_{\mathrm{ci}}^2}{(r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2)r^2}---\left(76\right)$

${\mathrm{\&sigma;}}_{\mathrm{\&theta;}1}=\frac{p_{\mathrm{ci}}{r}_{\mathrm{ci}}^2-p_{\mathrm{co}}{r}_{\mathrm{co}}^2}{r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2}-\frac{(p_{\mathrm{ci}}-p_{\mathrm{co}})r_{\mathrm{co}}^2{r}_{\mathrm{ci}}^2}{(r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2)r^2}---\left(77\right)$

式中，σ_r1为水泥浆外挤力和内压力产生的径向应力，σ_θ1为水泥浆外挤力和内压力产生的周向应力，r_ci为套管内半径，r_co为套管外半径。

套管的初始轴向载荷和初始轴向应力

初始轴向载荷:

$F_a={\mathrm{\&rho;}}_s\mathrm{g\&pi;}\underset{z}{\overset{\mathrm{Lc}}{\&Integral;}}(r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2)\cos \mathrm{\&alpha;dz}+\mathrm{\&pi;}[p_{\mathrm{ci}}\left(\mathrm{Lc}\right)r_{\mathrm{ci}}^2-p_{\mathrm{co}}\left(\mathrm{Lc}\right)r_{\mathrm{co}}^2]---\left(78\right)$

式中，F_a为初始轴向载荷，L_c为套管柱下入深度，p_co(Lc)为套管下端的外挤力，p_ci(Lc)为套管下端的内压力，ρ_s为钢的密度。

初始轴向应力:

${\mathrm{\&sigma;}}_{z1}=\frac{F_a}{\mathrm{\&pi;}(r_{\mathrm{co}}^2-r_{\mathrm{ci}}^2)}---\left(79\right)$

式中，σ_z1为初始轴向应力。

套管的弯曲应力过计算点做井眼轴线的密切平面，则套管弯曲发生在平面内。若在主法线方向上的坐标用y表示，则有：

σ_b＝E_cK_by,(-r_co≤y≤r_co)    (80)

式中，y为主法线方向的坐标，σ_b为弯曲应力，E_c为套管钢材弹性模量，K_b为井眼曲率。

下面列举具体实施例对本发明进行说明：

以位于区堤油田曲九馆三和曲104-X3两个区块.为例，采用三维地质建模软件PETREL建立该模拟区块的地质模型，根据本发明方法得出曲九馆三区块地应力计算结果的仿真图形。

根据曲九馆三区块及曲104-X3区块流固耦合数值模拟结果分析，主要得到如下的结论：

1 随着开发的进行，油层水平x向位移逐渐增大，总体水平位移在增大，特别是在近井地带，油层位移变化最为明显，说明在开采过程中，井壁地带发生应力集中，尤其在QTQ9-61、QTQ9-49、QTQ9-45、QTQ9-40、QTQ9-X44、QTQ9-35这些井的周围，其位移明显高于整个油层其他区域，而根据实际生产动态历史资料查明，QTQ9-45、QTQ9-40、QTQ9-X44、QTQ9-35已发生了套损；

2 随着开发的进行，两个区块的油层Z向位移逐渐增大，且最终Z向位移向下，说明开发一段时间后，地层产生竖向向下的位移，储层岩石骨架将发生下沉；

3 随着开发的进行，油层位移的大小在变化，位移的方向不统一，具有区域性分布特点，主要是受各个井的流场变化来控制，油层开发还未达到平衡；

4 储层的相邻层，以及不相邻层之间的水平X、Y位移明显不一样，这种差异在X、Y两个方向均表现得比较明显，尤其在有生产井的地方，位移比一般储层要大，套管由于位移差异而产生的拉伸荷载有可能超过其所能承受的抗拉强度，将会引起拉伸套损；

5 储层在开发过程中，其相邻层之间的位移始终存在差异，其主要原因可能是不同层之间注采不平衡引起的，或者是由于不同层之间物性差异大，造成开采中的层间矛盾加剧，诱发不同层位采出程度不一样，孔压的变化也将产生差异；

6 整个开发过程中，不论是水平方向，还是在竖直方向，层和层之间都有位移差异，这种差异将是诱发拉伸套损的一大原因；

7 对于Q104-X3断块，断层发育，这种层间位移差异很可能会诱发某些层的顺层滑动而对套管产生剪切作用，从而诱发另一种套管破坏模式：剪切套损

8 随着开发的进行，油层水平X向应力在增大，且各区域应力大小主要受井流场的影响，在各生产井附近均有不同程度的应力集中，这将会导致这些井承受不均匀集中载荷作用；

9 开发过程中，油层三个方向地应力均在生产井附近集中，曲104-X316井的三向应力集中最为明显，而实际上该井已于2008年发生套变；

10 由于两个区块存在较多的泥岩层段，这些泥岩在开发中如果吸水就会发生膨胀，强度急剧下降，在同等地应力作用下，将率先进入塑性破坏阶段，将其受到的荷载转移到套管上，将可能导致挤压套损；

根据注采过程流固耦合变化规律及注采参数变化对套管受力影响敏感性分析，注水开发过程预防套损的主要措施有：

1 控制生产井井底压力，降低生产压差，在开采过程中生产井井底压力逐渐下降，在开发初期，可采取较高的生产井流压，当油层压力普遍降低后再进一步降低井底流压来获取更大的采油速度，这样可避免注采过程中地应力变化幅度过大而使套管承受较大挤压力；

2 控制注水速度，在注水初期阶段注水量可采取逐渐增加的小量变化方式，相应的单井日注水量采取逐渐增加模式，可避免单井日注水量太大而引起套管内压力过高；

3 控制注水压力，在保证注水量的前提下尽量采取较低的注水压力(不超过岩石的临界荷载)，可避免由于注水压力过高会引起近井地带油层孔隙压力高导致岩石发生塑性破坏，而使得套管承受非均匀荷载；

4 控制采注比或生产井日产液量，可避免由于采注比过大而引起地层孔隙压力下降太快从而导致套管承受的挤压力过高；

5 采用抗挤强度高的材料制作套管，可适应油藏情况下由于地质条件的变化而导致的所需套管强度提高的要求；

6 对于要采取酸化压裂等强化采油措施的注水井，建议采用内壁有防止酸化腐蚀薄膜材料的套管，可适应地下酸化腐蚀环境带来的套管强度下降；

7 在油层生产段和非生产段过渡的地方尽量采用强度高的套管，或者采用预应力套管可为生产过程中套管变形协调预留空间，从而减少套坏发生的可能性。

8 在开发过程中，尽量分层注采，或者选择合理的层间组合方式开采，可避免层间吸水及注采的矛盾，而诱发拉伸或者剪切套损。

以上所述仅是对本发明的较佳实施方式而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施方式所做的任何简单修改，等同变化与修饰，均属于本发明技术方案的范围内。

Claims (1)

1.一种流固耦合作用下的油田套损计算方法，其特征在于，该流固耦合作用下的油田套损计算方法包括以下步骤进行：

步骤以：根据流体力学理论及油藏开发涉及的运动状态，采用理论分析方法建立渗流场数学模型，渗流场数学模型包括单相渗流数学模型、油水两相渗流数学模型和渗流场边界条件；

步骤二：根据弹塑性力学理论，在对储层岩石岩心数据进行详细分类及分析的基础上，建立储层开发过程的应力场计算数学模型，以岩心的强度测试数据为基础，采用DP准则来描述储层岩石的弹塑性本构关系；

步骤三：以有效应力原理为桥梁，根据孔隙介质变形理论，建立储层岩石的孔隙度与孔隙压力的关系及渗透率与孔压、有效应力的关系，结合室内岩心应力测试结果，建立储层注水开发过程流固耦合数学模型；

步骤四：建立相应区块的地质模型：

采用三维地质建模软件PETREL建立该模拟区块的地质模型，根据应力场计算的每个单元的节点坐标数据和和流场计算的块中心网格数据的几何对应关系，将PETREL建立的地质网格模型转化为应力场计算所需的有限元数值模型；

步骤五：将步骤四建立的地质模型导入eclipse，在eclipse软件中设置好流场边界及初始条件，根据步骤一中的流场单相流、两相流模型、边界类型及初始条件在eclipse中选择对应的流场模拟计算模型；结合每口井的生产情况，输入井的生产数据，在eclipse中设置好模型的边界及初始条件，计算模拟区域的孔隙压力；将步骤四转化出的有限元模型导入abaqus，将eclipse计算出的孔压导入abaqus作为初始应力条件，基于步骤二的应力场数学模型在abaqus中选择对应的应力计算模块，根据步骤三的渗透率随有效应力关系设计abaqus子程序来实现耦合计算过程中的有效应力及孔压对渗透率的影响计算，最终计算出储层的地应力及套管单元的受力，轴向载荷大小，根据套管管材的强度判断套管的损坏情况。