一种水力压裂裂缝延伸与多相流体流动模拟方法

1.本发明涉及油气田增产改造领域，具体涉及一种水力压裂裂缝延伸与多相流体流动模拟方法。


背景技术：

2.水力压裂技术是非常油气资源开采的重要手段，准确模拟压裂裂缝延伸轨迹是压裂施工设计优化的前提。
3.现有的水力压裂裂缝延伸轨迹模拟方法中将压裂液和地层流体考虑为同一种流体，但实际上，压裂液与地层流体通常为不混溶的两种流体，压裂液滤失到基质中并在高压下与地层流体形成两相流，流动过程需要克服毛管压力，各相流体的相对渗透率也会随其饱和度的变化发生变化，从而改变各相的流动性。研究两相流体流动的压裂裂缝延伸规律更贴近工程实际，因此，建立一种水力压裂裂缝延伸与多相流体流动模拟方法具有重要意义
4.目前常用的研究裂缝纵向扩展的数值模拟方法：边界元法、有限差分法、位移不连续法等在模拟裂缝延伸时需要预设裂缝扩展步长，且需要建立相交准则来判断裂缝扩展至界面后的延伸方向。基于相场法建立的裂缝延伸模型则不需要预先设定裂缝延伸方向和相交准则，该方法在处理裂缝延伸问题时不需要额外的规则来描述裂缝界面，裂缝延伸后不需要重新剖分网格。


技术实现要素：

5.本发明基于孔弹性理论、渗流力学、断裂力学、相场法等多学科知识，建立了一种新型裂缝扩展相场模型来描述固体变形与两相流体流动之间的耦合，考虑了裂缝扩展过程中的毛管压力。随着压裂液滤失到基质地层中，相场演化的驱动力和biot模量随之改变。引入各向异性绝对渗透率，且相渗透率与各相的饱和度有关。
6.发明内容
7.一种水力压裂裂缝延伸与多相流体流动模拟方法，包括以下步骤：
8.(1)收集工况和输入参数；
9.(2)建立多孔介质两相流应力平衡方程；
10.(3)建立两相流体流动控制方程；
11.(4)建立两相流相场演化方程；
12.(5)采用有限元数值离散方法和newton-raphson迭代法建立上述方程组的数值求解迭代格式；
13.(6)将步骤(1)中的工况和参数带入步骤(5)建立的求解迭代格式模拟不同工况下的水力裂缝延伸轨迹。
14.进一步的，所述步骤(1)中工况包括：裂缝性地层、多层叠置地层，输入参数包括：地应力参数、岩石弹性模量和泊松比、岩石抗拉强度、孔隙度、地层含油饱和度、地层压力、
地层流体和压裂液体积模量和粘度、岩石绝对渗透率、相渗曲线、压裂液注入时间、压裂液排量。
15.所述步骤(2)建立多孔介质两相流应力平衡方程，包括以下内容：
16.(2.1)多孔弹性岩石的两相流应力平衡方程为：
[0017][0018]
式中：σ
eff
为有效应力，可通过公式(2)表示；α(c)为biot系数，与相场变量有关，可通过公式(4)计算得到；i为单位张量；pw和po分别为水相压力和油相压力；sw和so分别为水相饱和度和油相饱和度。
[0019][0020]
式中，g(c)为衰减函数，可通过公式(3)计算得到；和分别表示拉应力和压应力。
[0021]
g(c)＝(1-c)2ꢀꢀꢀꢀꢀ
(3)
[0022][0023]
式中，k为多孔介质的岩石体积模量；ks为岩石骨架颗粒的体积模量。
[0024]
(2.2)应力边界条件
[0025][0026]
式中，为dirichlet边界上的位移；为作用在neumann边界上的应力；n为neumann边界上的方向向量。
[0027]
所述步骤(3)建立两相流体流动控制方程，包括以下内容：
[0028]
(3.1)多孔弹性岩石的两相流体流动控制方程为：
[0029][0030][0031]
式中，m(c)是biot模量，可通过公式(8)计算得到；θ是体积应变；k是各向异性绝对渗透率，可通过公式(7)计算得到；kw和ko为水和油的相渗透率，可由公式(9)计算得到；μw和μo为水和油的粘度。
[0032][0033]
式中，kw和ko为水和油的体积模量；φ为孔隙度。
[0034][0035]
式中，k
x0
和k
y0
为分别岩石基质的初始x方向和y方向的渗透率；k
fx
和k
fy
分别为裂缝的初始x方向和y方向的渗透率；w
x
和wy分别为x方向和y方向的渗透率加权系数，可通过公式(10)计算得到。
[0036][0037]
式中，w
x
和wy分别为x方向和y方向的缝宽，可由公式(11)和(12)计算得到；l
x
和ly分别为x和y方向的长度尺度因子。
[0038][0039][0040]
式中，ε
x
和εy分别为x和y方向应变；ε
cx
和ε
cy
分别为x和y方向的临界应变。
[0041]
(3.2)两相流体流动控制方程对应的边界条件
[0042][0043]
式中，和分别为dirichlet边界上的水相和油相压力；qw和qo为neumann边界上的水和油的注入量。
[0044]
所述步骤(4)建立的两相流相场演化方程，包括以下内容：
[0045]
(4.1)两相流相场演化方程为：
[0046][0047]
式中，为拉伸弹性应变能密度的最大值，可通过公式(14)得到。
[0048][0049]
式中，为拉伸应变能密度，可通过公式(15)计算得到。
[0050][0051]
式中，λ和g为拉梅常数；εi(i＝1，2，3)为主应变；macaulay括号算子《x》

＝(|x| x)/2和《x》-＝(|x|-x)/2用于区分拉伸和压缩状态。
[0052]
(4.2)相场演化方程对应的边界条件
[0053][0054]
进一步的，所述步骤(5)采用有限元数值离散方法和newton-raphson迭代法建立上述方程组的数值求解迭代格式，包括以下内容：
[0055]
将公式(6)和(7)相加，得到公式(17)
[0056][0057]
公式(1)、(7)、(14)、(17)组成非线性方程组，采用有限元数值离散方法对该非线性方程组进行离散，方程组(1)和(7)采用newton-raphson迭代法进行求解。在每个迭代步
内，裂缝相场不变，其迭代格式为：
[0058]
j(u，pw)iδ(u，pw)i r(u，pw)i＝0
ꢀꢀꢀꢀ
(18)
[0059]
式中，j(u，pw)i为雅可比矩阵，可通过公式(19)计算得到；δ(u，pw)i为第i个迭代步的位移和水相压力增量r(u，pw)为应力平衡方程和水相流体连续性方程的余量，可通过公式(20)和(21)计算得到
[0060][0061][0062][0063]
雅可比矩阵的各分量可通过公式(22)计算得到
[0064][0065][0066][0067][0068]
在公式(20)～(22)中，上标h代表节点的值，上标i-1代表第i-1个迭代步的值；上标n代表上一时间步的值。
[0069]
通过公式(18)计算得到第i个迭代步的位移增量和水相压力增量，计算过程中油相压力与水相饱和度为第i-1个迭代步的固定值。然后可通过公式(23)计算得到下一迭代步的位移和水相压力。
[0070]uh，i 1
＝u
h，i
δuh，
[0071]
计算得到第i 1个迭代步的位移和水相压力的试探解后，带入公式(24)求解第i 1个迭代步的油相饱和度。
[0072][0073]
式中，与可分别通过公式(25)和(26)计算得到。
[0074][0075][0076]
得到第i 1个迭代步的位移、水相压力和油相饱和度的试探解后，可通过公式(27)更新结点处的水相饱和度和油相压力，再带入公式(28)计算裂缝相场值。
[0077][0078]
式中，bm为常数。
[0079]
jcc
h，i 1
＝fcꢀꢀꢀꢀꢀ
(28)
[0080]
式中，jc和fc可分别由公式(28)和(29)计算得到。
[0081][0082][0083]
当位移、水相压力、油相饱和度和相场值都满足公式(30)所示的收敛条件时，迭代结束进入下一时间步，否则迭代继续进行。
[0084][0085]
上式中，和分别表示位移场、水相压力场、油相饱和度场和裂缝相场的收敛容差。
[0086]
本发明中的裂缝轨迹是自动确定的，裂缝延伸模型不需要预先设定裂缝延伸方向和相交准则，且可通过饱和度分布来描述压裂液的滤失情况。
附图说明
[0087]
图1为实施例1计算区域及边界条件示意图
[0088]
图2为实施例1第10秒的相场演化图
[0089]
图3为实施例1第20秒的相场演化图
[0090]
图4为实施例1第30秒的相场演化图
[0091]
图5为实施例1第10秒的含油饱和度分布图
[0092]
图6为实施例1第20秒的含油饱和度分布图
[0093]
图7为实施例1第30秒的含油饱和度分布图
[0094]
图8为实施例2计算区域及边界条件示意图
[0095]
图9为实施例2第10秒的相场演化图
[0096]
图10为实施例2第20秒的相场演化图
[0097]
图11为实施例2第30秒的相场演化图
具体实施方式
[0098]
下面结合不同工况对本发明做进一步的详细说明，但不构成对发明的任何限制。其中计算所用参数见表1、表2。
[0099]
表1裂缝性地层计算所采用基本参数表
[0100][0101]
表2多层叠置地层计算所采用基本参数表
[0102]
[0103][0104]
例1：如图1所示建立物理模型，计算区域为30m
×
16m的矩形区域，其中心处有一条沿最大水平地应力方向长度为4m的初始水力裂缝，天然裂缝的中心与注入点相距5m，天然裂缝与x方向夹角为45
°
。固定左边界在x方向上的位移和下边界在y方向上的位移。该区域由边长为0.25m的正方形单元剖分，共划分120
×
84个单元。注入压裂液排量为0.5m2/min，压裂液粘度为1mpa
·
s，注入时间30s，将表1中的参数带入本发明所建立的方程组进行模拟。模拟结果如图2～7所示(相场云图与含油饱和度云图具有一致性)，在该相交角下，水力裂缝两翼先沿着最大水平地应力方向延伸，上翼延伸至天然裂缝后转向，再沿天然裂缝右翼扩展，扩展至天然裂缝尖端后不再继续延伸，水力裂缝下翼则继续沿着最大水平地应力方向扩展。
[0105]
例2：如图8所示建立物理模型，计算区域为20m
×
20m的正方形区域，分为上、中、下三层，其中上下层为层厚7m岩性相同的阻隔层，中间层为层厚6m的产层。中间层的中心处有一条沿y方向长度为2m的初始水力裂缝。固定左边界在x方向上的位移和下边界在y方向上的位移。垂向应力为14mpa。该区域由边长为0.25m的正方形单元剖分，共划分80
×
80个单元。注入压裂液排量为0.3m2/min，压裂液粘度为1mpa
·
s，注入时间18s，将表2中的参数带入本发明所建立的方程组进行模拟。模拟结果如图9-11所示，在该参数条件下，产层中的水力裂缝先沿着垂向应力(y方向)扩展至交界面后转向，沿着交界面继续延伸，未进入上下阻隔层。