一种暂堵裂缝动态压力与扩展轨迹的预测方法与流程

1.本发明涉及石油天然气工程领域，尤其涉及一种暂堵裂缝动态压力与扩展轨迹的预测方法。


背景技术：

2.储层进行缝内暂堵压裂施工后会产生扩展裂缝，主流的裂缝扩展模型主要由流体连续性方程、裂缝缝宽方程、裂缝压降方程组成。裂缝扩展模型的建立经历了由简单到复杂，由二维到三维的转变，三维裂缝扩展模型中加入了裂缝缝高方程。发展过程中较为经典的模型是khristianovich和perkins等人建立的kgd和pkn模型，分别适用于缝高远大于缝长和缝长远大于缝高的情况，现常被用来进行对比实验。但在缝内暂堵重复压裂过程中，形成的初始裂缝会诱发井筒周围应力分布二次改变，产生诱导应力场，其会对后续裂缝起裂产生重要影响。
3.为了准确描述诱导应力场的方向与大小，sneddon和elliott最早提出用解析方法计算裂缝诱导应力场，并进一步推导得到了诱导应力场计算公式；后续学者又研究出了包括离散元方法、有限元方法和扩展有限元方法等的数值方法，计算了裂缝的诱导应力分布和相互干扰。解析法计算简单，但没有充分考虑地层性质、裂缝形状、流体作用等因素；数值方法求解精度高，但对参数输入要求也高，计算量相对较大。
4.为了解决这个难题，crouch等人创立了名为位移不连续方法的半解析法。与其他计算诱导应力场的方法相比，其只需进行一维网格划分，运算速度得到极大提升，同时能够考虑多种因素影响，计算精度较高。虽然建立了计算诱导应力场的新方法，但在使用位移不连续方法进行研究时没有考虑到流体渗流、压力扩散和孔隙压力变化带来的多孔基体变形和应力重分布，其会对诱导应力场产生影响从而影响裂缝扩展。
5.同时，对于所有的裂缝扩展模拟研究来说，精确计算多个裂缝中的流量分布是非常重要的。现在的流量分布计算方案基于流体的稳态分布理论，使用牛顿—拉夫逊法求解流量分布的非线性方程组，再通过展开近似解的泰勒级数、展开形成矩阵形式和高斯消去法等数学方法进行求解。但此方法中对于每个时间步长的每次迭代计算都需要求解雅可比矩阵，这不仅会导致多次求解使裂缝扩展模型运算复杂，矩阵的计算不准确还将导致计算结果不稳定甚至发散。


技术实现要素：

6.有鉴于此，本发明实施例的目的在于提供一种暂堵裂缝动态压力与扩展轨迹的预测方法。
7.为达到以上技术目的，本发明提供以下技术方案。
8.一种暂堵裂缝动态压力与扩展轨迹的预测方法，其特征在于，包括以下步骤：
9.通过孔隙介质位移不连续法，结合压力扩散方程建立水力裂缝诱导应力场模型，计算出地层任意点诱导应力分量的大小与孔隙压力；
10.通过对裂缝受力情况和薄弱点位置进行分析得到裂缝起裂模型，将水力裂缝诱导应力场、地应力场、封堵压力与摩阻压力进行叠加得到裂缝扩展应力场；
11.结合岩石变形模型，基于基尔霍夫定律和水电相似理论，使用阻力法得到转向新裂缝与老裂缝之间流量与压力分布情况；
12.结合裂缝扩展准则、缝宽缝长变化方程、裂缝扩展应力场与流量压力分配情况进行模拟得到裂缝扩展轨迹和其他参数变化情况。
13.进一步地，计算所需基本参数包括地层参数、压裂施工参数和暂堵参数。
14.进一步地，所述地层参数包括原始地层压力参数、地应力参数、岩石力学参数，所述压裂施工参数包括裂缝尺寸参数、压裂液参数，所述暂堵参数包括暂堵段长度、暂堵位置、扩展步长。
15.进一步地，所述通过孔隙介质位移不连续法与压力扩散方程结合还包括：
16.将裂缝基质界面流速与应力影响系数结合建立诱导应力场，将断裂段的应力变化和位移之间对于断裂变形随时间的局部关系离散化，与压力扩散方程中的孔压进行耦合求解求得与流速有关的应力边界影响系数。
17.进一步地，所述结合裂缝扩展准则、缝宽缝长变化方程、裂缝扩展应力场与流量压力分配情况进行模拟得到裂缝扩展轨迹和其他参数变化情况还包括：
18.耦合过程先从流量分布的第一个猜测开始，计算各段的压降，再分别对老裂缝压降、转向新裂缝压降、总压降和根段压力进行评估，然后使用耦合式得到流量分布新的猜测，重复该过程直到结果收敛。
19.本发明提供了一种暂堵裂缝动态压力与扩展轨迹的预测方法，该方法结合压力扩散方程建立水力裂缝诱导应力场模型，可以计算出地层任意点诱导应力分量的大小与孔隙压力。通过对裂缝受力情况和薄弱点位置进行分析得到裂缝起裂模型，将水力裂缝诱导应力场、地应力场、封堵压力与摩阻压力进行叠加得到裂缝扩展应力场；结合岩石变形模型，基于基尔霍夫定律和水电相似理论，使用阻力法得到转向新裂缝与老裂缝之间流量与压力分布情况；结合裂缝扩展准则、缝宽缝长变化方程、裂缝扩展应力场与流量压力分配情况进行模拟得到裂缝扩展轨迹和其他参数变化情况。本技术方案开发了一种新方法用于研究缝内暂堵裂缝转向的技术原理和计算新老裂缝间的流量与压力分配情况，可用于模拟典型压裂井组的缝内暂堵裂缝扩展轨迹，得到裂缝参数变化情况的同时使计算效率得到提高。
附图说明
20.图1为本发明裂缝边界离散与坐标变换示意图。
21.图2为本发明缝内暂堵裂缝起裂示意图。
22.图3为本发明暂堵裂缝流量压力分配示意图。
23.图4为本发明阻力法耦合示意图。
24.图5
‑
图8为本发明计算实例的结果。
具体实施方式
25.结合附图和本发明具体实施方式的描述，能够更加清楚地了解本发明的细节。但是，在此描述的本发明的具体实施方式，仅用于解释本发明的目的，而不能以任何方式理成
是对本发明的限制。在本发明的教导下，技术人员可以构想基于本发明的任意可能的变形，这些都应被视为属于本发明的范围。
26.在本发明中提出了一种暂堵裂缝动态压力与扩展轨迹的预测方法，该方法包括以下步骤：
27.1、收集计算所需基本参数，计算所需基本参数包括地层参数、压裂施工参数、暂堵参数，地层参数包括原始地层压力参数、地应力参数、岩石力学参数，压裂施工参数包括裂缝尺寸参数、压裂液参数，暂堵参数包括暂堵段长度、暂堵位置、扩展步长。
28.2、通过孔隙介质位移不连续法，结合压力扩散方程建立水力裂缝诱导应力场模型，计算出地层任意点诱导应力分量的大小与孔隙压力；
29.根据crouch等人创立的位移不连续方法，可得裂缝上下两个面的位移量可分别用u(x,0

)和u(x,0
‑
)表示，当|x|<裂缝半长时，裂缝法向位移不连续量d
n
和切向位移不连续量d
s
与位移量之间还存在如下关系：
[0030][0031]
式中：d
s
，d
n
——裂缝切向(即x轴方向)和法向(即y轴方向)的位移不连续量，m；u
x
(x,0
‑
)，u
x
(x,0

)——裂缝下表面和上表面在x轴方向错动量， m；u
y
(x,0
‑
)，u
y
(x,0

)——裂缝下表面和上表面在y轴方向错动量，m。
[0032]
选取上下表面利用格林函数进行位移量计算，然后规定y＝0时的连续性条件，再依据弹性力学理论对位移场进行亥姆霍兹分解，将应力表达式与位移表达式进行组合，最后联立上下两表面表达式可得到裂缝单元体的d
s
、d
n
以及地层任意处的位移场与地层任意处的应力分量为：
[0033][0034]
式中：u
x
，u
y
——水平和垂直位移，m；ν——泊松比，无量纲；σ
xx
，σ
yy
——x和y方向上正应力，mpa；τ
xy
——剪切应力，mpa；g——剪切模量， mpa；f(x,y)——只与单元坐标有关的模拟函数，f下标表示其对x，y的各阶导数，a表示裂缝单元体半长。
[0035]
裂缝渗流主要分为沿裂缝面切向和法向的流动。切向流的计算过程可用立方定律表示，但实际使用过程与该定律假设情况矛盾，故认为缝内压强均匀，忽略切向流影响。
[0036]
法向流方面，流体向岩体的渗透行为可以认为是扩散过程。点源密度基本方程通过积分变换法求解，再结合质量守恒定律求出待定系数，可得瞬时点源解的密度函数：
[0037][0038]
式中：ρ(x,y,t)——注入流体后平面内流体密度，kg/m3；ρ
i
——初始平面内流体密度，kg/m3；m——注入流体质量，kg；χ——扩散系数，无因次；x
′
， y
′
——注入流体点坐标，无因次；t——注入时刻，无因次；τ——注入停止时刻，无因次。
[0039]
微可压缩流体的密度与压力之间呈线性关系，可由关联式得到平面场点任意时刻的孔隙压力。若对注入点持续注入流体，即从0时刻到τ1时刻持续注入，则需对其进行时间维度的积分，假设注入速率恒定，称此方程称为压力扩散方程：
[0040][0041]
式中：p(x,y,t)——流体扩散后的孔隙压力，mpa；p
i
——初始孔隙压力， mpa；c
f
——压缩系数，mpa
‑1；v——注入流体体积，m3；q
int
——每单位裂缝长度每单位地层厚度的界面流速，m/s；
[0042]
在具有多个裂缝段的裂缝网络中，裂缝之间存在相互作用。将裂缝离散成 n个边界单元，第i个断裂段的应力变化受系统中所有断裂段变形的影响。第i 个断裂段的正应力和剪应力变化可从法向位移不连续值和切向位移不连续值的依赖关系中导出。通过结合恒定流体注入与产出体积率和裂缝段孔隙弹性位移不连续性的解，单元i中点的剪应力和法应力可以由j单元的位移不连续量与其裂缝内流体流速计算得到：
[0043][0044]
式中：——第j个断裂单元对第i个断裂单元的影响系数，物理意义为单元j上的单位切(法)向常位移不连续量相应的单元i中点的剪切力(法向力)，mpa。
[0045]
在本技术方案中位移不连续解和界面流速是相互依赖的。位移不连续性值会随时间变化，则从初始齐次状态开始每个边界积分避免了体积计量积分的需要。因此，在对每个时间步的增量进行数值积分时，必须包括该次增量的所有先前增量。根据位移不连续增量计算了第i个裂缝段上的诱导应力：
[0046][0047]
式中：——第ξ时间步第j个断裂段的对应物理量增量；——第ξ时间步第j个断裂段的对应物理量增量；——时间步长h处第j个断裂段的先前对应物理量增量，从1求和到ξ
‑
1；和——第j个断裂单元在时间步h对第i个断裂单元的影响系数。
[0048]
实际情况下单元坐标系和整体坐标系之间有一定夹角，如图1所示需在两坐标系之间进行坐标转换才可列出正确的方程。最后结合局部和全局关系将应力分量表达式向单元i的局部坐标系变换，即可得到应力边界影响系数的表达式：
[0049][0050]
对于裂缝网络中的任何裂缝，变形必须符合裂缝变形模型。如图1所示将裂缝网络离散为小裂缝段，并通过结合岩石变形模型将第i个断裂段的有效法向应力变化和法向位移之间对于断裂变形随时间的局部关系离散化，与压力扩散方程中的孔压进行耦合求解可求得与流速有关的应力边界影响系数。
[0051][0052]
式中：k
n
——法向刚度，n/m，；k
m
——系统体积模量，pa；k
s
——固体体积模量，pa；φ
d
为扩张角，
°
。
[0053]
3、通过对裂缝受力情况和薄弱点位置进行分析得到裂缝起裂模型，将水力裂缝诱导应力场、地应力场、封堵压力与摩阻压力进行叠加得到裂缝扩展应力场；
[0054]
如图2所示在已失效的裂缝内泵入暂堵剂，将产生两个薄弱区：薄弱区1 和薄弱区2，薄弱区壁面的受力实际上是由水力裂缝延伸产生的诱导应力和地层流体压力变化产生的诱导应力与原地应力的叠加。
[0055]
规定缝内暂堵起裂过程中转向新裂缝转向角度为：
[0056][0057]
式中：γ
β
——转向角度，(
°
)；τ
xy
——在xoy坐标下剪应力分量，mpa；σ
x
——在x坐标下正应力分量，mpa；σ
y
——在y坐标下正应力分量，mpa。
[0058]
裂缝扩展应力场主要由封堵压力、周向应力、摩擦阻力组成。
[0059]
封堵产生的压力主要由暂堵剂堆积产生，由封堵渗透率结合达西公式得出，本技术方案根据封堵情况直接给出封堵压力具体数值。
[0060]
周向应力是指物体沿着圆周方向产生的应力，在裂缝延伸过程中其指代的是让裂缝闭合的力，即原地应力与诱导应力的叠加。
[0061]
根据之前建立的直角坐标系，则人工裂缝周围一点i受到的x方向和y方向上的应力分量可以转化为：
[0062][0063]
式中：g
ij
——三维修正系数，无量纲；d
ij
——裂缝单元i与j的距离，mm； h
f，ji
——裂缝高度，mm；ε，η——经验系数，通常取ε＝1，η＝2.3；σ
h
、σ
h
——表示最大最小水平主应力，mpa；ν——储层岩石泊松比，无因次；α—— biot多孔弹性系数，无因次；p
p
——当前地层压力，mpa；p
e
——原始地层压力，mpa。
[0064]
在x和y坐标下正应力和剪应力分量转化到裂缝γ
β
角方向的应力，可以得到在γ
β
角方向坐标下的正应力和剪应力分量。
[0065]
扩展宽度方程为：
[0066][0067]
式中：w
f,i
——扩展缝宽，m；v——岩石泊松比，无因次；e——岩石弹性模量，mpa。
[0068]
根据裂缝流动方程可以计算出转向新裂缝与老裂缝扩展所产生的压降为：
[0069][0070]
式中：μ
w
——流体粘度，mpa
·
s。
[0071]
运用体积守恒求取扩展所需要的总时间为：
[0072][0073]
式中：c——滤失系数，m/min
0.5
，一般取0.0004；δl
f，i
——扩展缝长， m。
[0074]
转向新裂缝扩展到第i段时第j段缝长等于步长，老裂缝扩展缝长为：
[0075][0076]
扩展缝高方程为：
[0077]
[0078]
式中：k
ic
——岩石断裂韧性，mpa
·
m
0.5
。
[0079]
裂缝发生非平面延伸时，裂缝扩展由岩体ⅰ型破坏和ⅱ型破坏共同影响。 erdogan和sih最早提出了
ⅰ‑ⅱ
型复合裂缝扩展理论。
[0080]
最大周向应力准则为：
[0081][0082]
式中：kⅰ——ⅰ型应力强度因子，mpa
·
m
0.5
；kⅱ——ⅱ型应力强度因子， mpa
·
m
0.5
。
[0083]
式(17)左侧等于或大于右侧时岩体破坏，破坏方向与周向应力极值方向相同，这也是裂缝进一步扩展的方向。
[0084]
基于此，olson提出kⅰ和kⅱ可以根据裂缝末端法向和切向位移不连续量得到：
[0085][0086]
同时根据弹性力学理论，在周向应力取最大值方向，裂缝末端剪切应力一定等于0：
[0087][0088]
岩石发生破裂，裂缝开始延伸后，裂缝的几何参数和裂缝内的流体流量与缝内压力都会发生变化。在之前研究内容的基础上分别用到了基尔霍夫第一和第二定律，得到了具体表达式后，再通过阻力法对数值方程式进行耦合来方便进行计算。
[0089]
暂堵裂缝单元离散模型如图3所示，示意图展示了某一井筒主裂缝加入暂堵剂形成转向新裂缝(hf
12
)和老裂缝(hf
14
)的情形。根据物理模型和基本假设，将转向新裂缝离散为j1个单元，将老裂缝离散为j1 ns个单元。
[0090]
基于kirchoff第一定律，在进行暂堵压裂时，裂缝扩展到第i段时，水力裂缝总流量为q1(i)，总流量被分到转向新裂缝和老裂缝，转向新裂缝流量为 q
12
(i)，则老裂缝流量q
14
(i)为q1(i)
‑
q
12
(i)，即：
[0091]
q1(i)＝q
12
(i) q
14
(i),(i<j1)
ꢀꢀꢀꢀꢀꢀꢀ
(20)
[0092]
式中：q1(i)——压裂液注入到转向新裂缝的第i段时主裂缝流量，m3/min。
[0093]
又基于kirchoff第二定律，将转向新裂缝根端o2作为参考点，建立主裂缝中流体压力平衡准则。即转向新裂缝根端o2的压力等于转向新裂缝缝内的压力损失、转向新裂缝根部周向应力之和，也等于暂堵段暂堵压差、主裂缝hf
13
段压力损失、老裂缝hf
14
段摩阻损失以及老裂缝hf
14
端部周向应力之和。则扩展到第i段转向新裂缝时压力平衡方程为：
[0094][0095]
式中：——主裂缝o2端流体压力，mpa；δp
ft,ji
——转向新裂缝扩展到第i段时
根据式(13)计算出的第j段裂缝压降，mpa；σ
βy,ft,i
——转向新裂缝扩展到第i段时根据式(11)计算的周向应力，mpa；δp
t,i
——老裂缝扩展到第i段时的暂堵压差，mpa；——老裂缝扩展到第i段时的主裂缝hf
13
部分的裂缝压降，mpa；δp
fz,ji
——老裂缝扩展到第i段时的老裂缝hf
14
部分的第j段裂缝压降，mpa；σ
βy,fz,i
——老裂缝扩展到第i段时的周向应力，mpa。
[0096]
结合扩展应力场与扩展准则，即可根据压力、裂缝缝长缝宽和流量分布计算结果判断水力裂缝是否发生扩展并确定延伸方向。
[0097]
4、结合岩石变形模型，基于基尔霍夫定律和水电相似理论，使用阻力法得到新裂缝与老裂缝之间流量与压力分布情况；
[0098]
如图4所示，n型断裂系统类似于具有n个分支的电路网络，使用类似的电阻概念来计算流量在多个裂缝中的分布。假设每个裂缝都有一定的阻力，那么裂缝接收的流量与其阻力成反比。
[0099]
通过将式(20)和式(21)改写为式(22)来定义转向新裂缝和老裂缝的阻力：
[0100][0101]
根据式(22)，可得裂缝的阻力是每单位流量的总压降。则每个裂缝的流量速率可以写成：
[0102][0103]
式(23)代入式(21)得：
[0104][0105]
定义r
tot
为总阻力，表达式与并联电路中的总电阻相仿：
[0106][0107]
把式(25)代入式(24)得：
[0108]
[0109]
式(22)代入式(20)得到：
[0110][0111]
裂缝断裂阻力r
i
由式(22)定义，总阻力r
tot
由式(25)定义，都是流量速率q
12(i)
、q
14(i)
的函数；式(27)代表未知数与方程个数相同的不动点迭代公式，可以证明这个非线性方程组至少线性地局部收敛。
[0112]
耦合的最后结果可以表示为：
[0113][0114]
在与模型耦合过程中，先从流量分布的第一个猜测开始，计算各段的压降，再分别对r
12(i)
、r
14(i)
、r
tot
、进行评估，然后使用式(28)得到流量分布新的猜测，最后重复该过程直到结果收敛。
[0115]
5.计算实例与分析
[0116]
结合裂缝扩展准则、缝宽缝长变化方程、裂缝扩展应力场与流量压力分配情况进行模拟得到裂缝扩展轨迹和其他参数变化情况。
[0117]
下面以一个实际算例作为示例，基于本技术中公开的方法，进行了预测裂缝轨迹、压力、流量动态等参数预测，分析了参数变化的原因。
[0118]
(1)基本参数
[0119]
表1基础参数
[0120][0121]
(2)工况计算结果
[0122]
根据表1基础参数进行了运算，计算结果如图5
‑
图8所示。
[0123]
如图5所示，暂堵剂泵入后在原有水力裂缝的基础上形成了转向新裂缝和老裂缝，新裂缝起裂点符合缝内暂堵起裂物理模型中的薄弱点位置，并且老裂缝延伸距离远远小于转向新裂缝延伸距离，形成了稳定的封堵。
[0124]
如图6所示，刚开始由于暂堵剂的泵入，使得暂堵压力增加，从而使得转向新裂缝以高于老裂缝的排量开始扩展，但随着转向新裂缝有向着水平裂缝转向的趋势，下一时刻，转向新裂缝周向应力逐渐减小，转向新裂缝排量下降，老裂缝排量上升，最后两者稳定下来。值得说明的是，转向新裂缝的排量与老裂缝排量之和为9m3/min。
[0125]
如图7所示，刚开始由于暂堵剂的泵入，使得暂堵压力增加，但因为转向新裂缝有向着水平裂缝转向的趋势，下一时刻，转向新裂缝周向应力突然减小，缝内压力有一个减小的过程；随着时间的增加，由于缝内摩阻的增加大于了周向应力的减小，使得压力开始逐渐增加；在10min左右由于转向新裂缝方向愈加趋于最大水平主应力方向，裂缝周向应力有增加的过程，后续同样由于缝内摩阻的增加大于了周向应力的减小，使得压力开始逐渐增加。
[0126]
如图8所示，老裂缝最大缝宽0.57mm，转向新裂缝最大缝宽2.64mm；并且随着扩展长度的增加裂缝宽度逐渐减小，模拟情况符合实际压裂情况中裂缝宽度分布情况。
[0127]
以上通过实施例对本发明进行具体描述，有必要在此指出的是，本实施例仅是本发明的优选实施例，并非对本发明作任何限制，也并非局限于本文所披露的形式，不应看作是对其他实施例的排除。而本领域人员所进行的改动和简单变化不脱离本发明技术思想和范围，则均属于本发明技术方案的保护范围内。