本发明提供一种裂缝性地层压裂裂缝延伸轨迹预测方法，包括以下步骤：(1)收集地层参数；(2)收集施工参数；(3)建立多孔介质中压裂裂缝延伸控制方程组强形式；(4)建立多孔介质中压裂裂缝延伸控制方程组弱形式；(5)控制方程组有限元离散；(6)控制方程组线性化及迭代格式建立；(7)将步骤(1)和(2)获得的参数输入步骤(6)建立的迭代计算方程组，预测该参数下的裂缝延伸轨迹。本发明在裂缝延伸后不需要重新剖分网格，不需要引用额外准则来追踪裂缝轨迹。
1.一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，包括以下步骤：(1)收集地层参数；
(3)建立多孔介质中压裂裂缝延伸控制方程组强形式；
(4)建立多孔介质中压裂裂缝延伸控制方程组弱形式；
(7)将步骤(1)和(2)获得的参数输入步骤(6)建立的迭代计算方程组，预测该参数下的裂缝延伸轨迹。
2.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(1)中收集的地层参数包括：地应力参数、岩石弹性模量和泊松比、岩石本体临界粘聚力、天然裂缝长度、天然裂缝方向角、天然裂缝渗透率、天然裂缝临界粘聚力。
3.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(2)中收集的施工参数包括：注入排量、注入时间、压裂液粘度。
4.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(3)建立多孔介质中压裂裂缝延伸控制方程组强形式，包括以下内容：(3.1)多孔介质应力平衡方程建立
根据达西渗流定理以及Biot孔弹性理论可得到多孔介质中流体流动压力计算方程：式中：M为Biot模量，MPa，可通过公式(3)计算得到；α为Biot系数，可通过公式(4)计算得到；εii为体积应变；k为各向异性渗透率张量，m2，可通过公式(5)计算得到；μf为流体粘度，Pa.s；
式中：Kf为流体的体积模量，MPa；φ为多孔介质的孔隙度；K为多孔介质的体积模量，MPa；Ks为岩石骨架颗粒的体积模量，MPa；
式中：k0为岩石的初始渗透率，m2；k为岩石在受到外界力影响后的渗透率，m2，建立一种如公式(6)所示的计算方式；θ为最大主应变的方向角，如公式(7)所示：式中：b1，b2和b3为待定参数，需要通过实验数据拟合得到；ε1为最大主应变；εy为Y方向的应变；γxy为剪应变：(3.3)基于相场法的裂缝延伸控制方程建立
在水力压裂过程中，每单位体积完全饱和的多孔介质的总自由能可分解为三个部分，即：ψ(ε，ζ,c,▽c)＝ψeff(ε)+ψfluid(ε，ζ)+ψfrac(c,▽c) (8)式中：Ψeff为储存于岩石骨架中的弹性应变能密度，MPa；Ψfluid为流体能量密度，MPa；
岩石在损伤过程中会造成岩石骨架拉伸弹性应变能密度降低，而压缩弹性应变能密度保持不变，则将公式(8)中岩石骨架弹性应变能密度拆分为受损拉伸部分和压缩部分，即：式中：g(c)为衰减函数，定义衰减函数如公式(10)所示；ψ+eff和ψ-eff分别为拉伸和压缩弹性应变能密度，可通过公式(11)计算得到，MPa：g(c)＝(1-c)2 (10)
式中：λ和G为拉梅常数，MPa；εi(i＝1，2，3)为主应变；函数+＝(|x|+x)/2，-＝(|x|-x)/2；
式(8)中流体能量密度和裂缝能量密度分别可通过公式(12)和(13)计算得到：式中：l0为长度尺度参数，用于控制扩散裂缝区域宽度，m；Gc为临界能量释放率，MPa.m，可通过公式(14)计算得到；
裂缝相场演化由微力平衡方程描述；不考虑微惯性和外部微力的影响，多孔介质的微力平衡方程为：▽·Hm-Km＝0 (15)
式中：Hm为微牵引力，MPa.m；Km为内部微观力量，MPa；
式中：ρ为多孔介质密度，kg/m3；e为熵，J/(mol.K)；ζ为流体体积含量增量，如公式(17)所示：压裂过程可视为等温过程，则在该过程中总机械耗散是非负的，即可采用Clausius-Duhem不等式描述该过程：上式中总自由能密度对时间的导数可表示为：
将方程(16)、和(19)～(22)带入方程(18)可得到：所有热力学过程都必须满足上述不等式，也就是说 和 取任意值上述不等式都成立，则：将式(24)带入式(1)，应力平衡方程可改写为：式中：I为单位张量，二维情况下为[1 1 0]T；σ+eff为有效拉应力，MPa；σ-eff为有效压缩应力，MPa；
将方程(26)和(27)带入方程(15)可得到，压裂过程中裂缝相场演化方程：其中：
应力平衡方程、流体压力计算方程、裂缝相场演化方程所对于的边界条件分别如公式(31)～(33)所示；
方程(31)～(33)中， 和Γ分别为位移场、压力场和裂缝相场的Dirichlet边界， 和 分别为位移场、压力场和裂缝相场的Neumann边界。
5.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(4)建立多孔介质中压裂裂缝延伸控制方程组弱形式，包括以下内容：将偏微分方程(28)、(2)和(29)分别与乘以试探函数wu、wp和wc并在计算域积分，然后采用散度定理并结合边界条件(31)～(33)，可得到控制方程的弱形式：
6.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(5)控制方程组有限元离散，包括以下内容：采用4节点四边形单元对计算区域进行有限元离散，对于每个计算单元，位移、压力、裂缝相场和对应权函数以及对应梯度的插值形式分别如公式(37)和式(38)所示：式中：uh、ph和ch分别为计算节点处的位移、压力和裂缝相场值；whu、whp和whc分别为计算节点处的位移、压力和裂缝相场试探函数值；Nu、NP和Nc分别为位移场、压力场和裂缝相场的插值形函数；Bu和Buvol分别为应变矩阵和体积应变矩阵；Bp和Bc分别为压力和裂缝相场插值形函数的导数矩阵；
将方程(37)和方程(38)分别带入方程(34)～(36)可得到：式(40)中下标n代表第n个时间步的值；Δt为时间步长，s。
7.根据权利要求1所述的一种裂缝性地层压裂裂缝延伸轨迹预测方法，其特征在于，所述步骤(6)控制方程组线性化及迭代格式建立，包括以下内容：采用Newton–Raphson(NR)迭代法求解渗流-应力耦合方程组(39)和(40)，第i个迭代步的Newton–Raphson(NR)迭代格式为：其中：
通过方程(42)可求得第i个迭代步的位移增量δuh和压力增量δph后，进而可得到第i+1个迭代步位移和压力的试探解，即：获得第i+1个迭代步位移和压力的试探解后，然后通过求解方程(47)可得到裂缝相场值，在该过程中位移和压力值是固定的：其中
当位移、压力和裂缝相场都满足如公式(48)所示的收敛条件时，则迭代结束，进入下一时间步的计算，否则迭代继续进行；
||Ru||≤||Ru0||，||Rp||≤||Rp0||，||ci+1-ci||≤||Rc0|| (50)上式中，Ru0、Rp0和Rc0分别为位移场、压力场、裂缝相场收敛容差。
一种裂缝性地层压裂裂缝延伸轨迹预测方法
技术领域
[0001] 本发明涉及油气田增产改造领域，具体涉及一种裂缝性地层压裂裂缝延伸轨迹预测方法。
背景技术
[0002] 中国页岩储层、煤层、致密砂岩储层等非常规储层中含有丰富的天然气资源，但这些储层由于渗透率低，需要通过水力压裂技术在这类地层中建造人工裂缝，才能将这类地层中的天然气资源开采出来。在压裂施工前采用压裂裂缝模拟模型预测裂缝轨迹，有助于优化施工方案。而在压裂施工后，采用压裂裂缝模拟模型预测裂缝形态，有利于压后评估。因此压裂裂缝延伸轨迹预测技术成为水力压裂优化设计的核心技术。但由于页岩、煤岩、致密砂岩等储层中含有大量的天然裂缝，而这些天然裂缝的存在会影响水力裂缝延伸轨迹。
因此建立一种能够预测裂缝性地层中压裂裂缝延伸轨迹的方法对于裂缝性地层压裂优化设计具有重要意义。
发明内容
[0003] 本发明综合应用渗流力学、岩石力学、有限元理论、相场法理论等多学科知识，建立一种裂缝性地层压裂裂缝延伸轨迹预测方法。
[0004] 一种裂缝性地层压裂裂缝延伸轨迹预测方法，包括以下步骤：
[0005] (1)收集地层参数；
[0006] (2)收集施工参数；
[0007] (3)建立多孔介质中压裂裂缝延伸控制方程组强形式；
[0008] (4)建立多孔介质中压裂裂缝延伸控制方程组弱形式；
[0009] (5)控制方程组有限元离散；
[0010] (6)控制方程组线性化及迭代格式建立；
[0011] (7)将步骤(1)和(2)获得的参数输入步骤(6)建立的迭代计算方程组，预测该参数下的裂缝延伸轨迹。
[0012] 进一步的，所述步骤(1)中收集的地层参数包括：地应力参数，、岩石弹性模量和泊松比、岩石本体临界粘聚力、天然裂缝长度、天然裂缝方向角，、天然裂缝渗透率、天然裂缝临界粘聚力。
[0013] 所述步骤(2)中收集的施工参数包括：注入排量、注入时间、压裂液粘度。
[0014] 所述步骤(3)建立多孔介质中压裂裂缝延伸控制方程组强形式，包括以下内容：
[0015] (3.1)多孔介质应力平衡方程建立
[0016] 多孔弹性岩石的应力平衡方程为：
[0017]
[0018] 式中：σ为应力张量，MPa。
[0019] (3.2)流体压力计算方程建立
[0020] 根据达西渗流定理以及Biot孔弹性理论可得到多孔介质中流体流动压力计算方程：
[0021]
[0022] 式中：M为Biot模量，MPa，可通过公式(3)计算得到；α为Biot系数，可通过公式(4)计算得到；εii为体积应变；k为各向异性渗透率张量，m2，可通过公式(5)计算得到；μf为流体粘度，Pa.s；
[0023]
[0024]
[0025] 式中：Kf为流体的体积模量，MPa；φ为多孔介质的孔隙度；K为多孔介质的体积模量，MPa；Ks为岩石骨架颗粒的体积模量，MPa；
[0026]
[0027] 式中：k0为岩石的初始渗透率，m2；k为岩石在受到外界力影响后的渗透率，m2，建立一种如公式(6)所示的计算方式；θ为最大主应变的方向角，如公式(7)所示：
[0028]
[0029]
[0030] 式中：b1，b2和b3为待定参数，需要通过实验数据拟合得到；ε1为最大主应变；εy为Y方向的应变；γxy为剪应变：
[0031] (3.3)基于相场法的裂缝延伸控制方程建立
[0032] 在水力压裂过程中，每单位体积完全饱和的多孔介质的总自由能可分解为三个部分，即：
[0033]
[0034] 式中：Ψeff为储存于岩石骨架中的弹性应变能密度，MPa；Ψfluid为流体能量密度，MPa；Ψfrac为裂缝能量密度，MPa；
[0035] 岩石在损伤过程中会造成岩石骨架拉伸弹性应变能密度降低，而压缩弹性应变能密度保持不变，则将公式(8)中岩石骨架弹性应变能密度拆分为受损拉伸部分和压缩部分，即：
[0036]
[0037] 式中：g(c)为衰减函数，定义衰减函数如公式(10)所示；ψ+eff和ψ-eff分别为拉伸和压缩弹性应变能密度，可通过公式(11)计算得到，MPa：
[0038] g(c)＝(1-c)2 (10)
[0039]
[0040] 式中：λ和G为拉梅常数，MPa；εi(i＝1，2，3)为主应变；函数+＝(|x|+x)/2，-＝(|x|-x)/2；
[0041] 式(8)中流体能量密度和裂缝能量密度分别可通过公式(12)和(13)计算得到：
[0042]
[0043]
[0044] 式中：l0为长度尺度参数，用于控制扩散裂缝区域宽度，m；Gc为临界能量释放率，MPa.m，可通过公式(14)计算得到；
[0045]
[0046] 式中：σc为临界粘聚力，MPa；
[0047] 裂缝相场演化由微力平衡方程描述；不考虑微惯性和外部微力的影响，多孔介质的微力平衡方程为：
[0048]
[0049] 式中：Hm为微牵引力，MPa.m；Km为内部微观力量，MPa；
[0050] 多孔介质内部能量平衡方程可写为：
[0051]
[0052] 式中：ρ为多孔介质密度，kg/m3；e为熵，J/(mol.K)；ζ为流体体积含量增量，如公式(17)所示：
[0053]
[0054] 压裂过程可视为等温过程，则在该过程中总机械耗散是非负的，即可采用Clausius-Duhem不等式描述该过程：
[0055]
[0056] 上式中总自由能密度对时间的导数可表示为：
[0057]
[0058] 其中：
[0059]
[0060]
[0061]
[0062] 将方程(16)、和(19)～(22)带入方程(18)可得到：
[0063]
[0064] 所有热力学过程都必须满足上述不等式，也就是说 和 取任意值上述不等式都成立，则：
[0065]
[0066]
[0067]
[0068]
[0069] 将式(24)带入式(1)，应力平衡方程可改写为：
[0070]
[0071] 式中：I为单位张量，二维情况下为[1 1 0]T；σ+eff为有效拉应力，MPa；σ-eff为有效压缩应力，MPa；
[0072] 将方程(26)和(27)带入方程(15)可得到，压裂过程中裂缝相场演化方程：
[0073]
[0074] 其中：
[0075]
[0076] (3.4)控制方程对应边界条件
[0077] 应力平衡方程、流体压力计算方程、裂缝相场演化方程所对于的边界条件分别如公式(31)～(33)所示；
[0078]
[0079]
[0080]
[0081] 方程(31)～(33)中， 和Γ分别为位移场、压力场和裂缝相场的Dirichlet边界， 和 分别为位移场、压力场和裂缝相场的Neumann边界。
[0082] 所述步骤(4)建立多孔介质中压裂裂缝延伸控制方程组弱形式，包括以下内容：
[0083] 将偏微分方程(28)、(2)和(29)分别与乘以试探函数wu、wp和wc并在计算域积分，然后采用散度定理并结合边界条件(31)～(33)，可得到控制方程的弱形式：
[0084]
[0085]
[0086]
[0087] 所述步骤(5)控制方程组有限元离散，包括以下内容：
[0088] 采用4节点四边形单元对计算区域进行有限元离散，对于每个计算单元，位移、压力、裂缝相场和对应权函数以及对应梯度的插值形式分别如公式(37)和式(38)所示：
[0089]
[0090]
[0091] 式中：uh、ph和ch分别为计算节点处的位移、压力和裂缝相场值；whu、whp和whc分别为计算节点处的位移、压力和裂缝相场试探函数值；Nu、NP和Nc分别为位移场、压力场和裂缝相场的插值形函数；Bu和Buvol分别为应变矩阵和体积应变矩阵；Bp和Bc分别为压力和裂缝相场插值形函数的导数矩阵；
[0092] 将方程(37)和方程(38)分别带入方程(34)～(36)可得到：
[0093]
[0094]
[0095]
[0096] 式(40)中下标n代表第n个时间步的值；Δt为时间步长，s；
[0097] 所述步骤(6)控制方程组线性化及迭代格式建立，包括以下内容：
[0098] 本发明采用Newton–Raphson(NR)迭代法求解渗流-应力耦合方程组(39)和(40)，第i个迭代步的Newton–Raphson(NR)迭代格式为：
[0099]
[0100] 其中：
[0101]
[0102]
[0103]
[0104]
[0105]
[0106]
[0107] 通过方程(42)可求得第i个迭代步的位移增量δuh和压力增量δph后，进而可得到第i+1个迭代步位移和压力的试探解，即：
[0108]
[0109] 获得第i+1个迭代步位移和压力的试探解后，然后通过求解方程(47)可得到裂缝相场值，在该过程中位移和压力值是固定的：
[0110]
[0111] 其中
[0112]
[0113]
[0114] 当位移、压力和裂缝相场都满足如公式(48)所示的收敛条件时，则迭代结束，进入下一时间步的计算，否则迭代继续进行；
[0115] ||Ru||≤||Ru0||，||Rp||≤||Rp0||，||ci+1-ci||≤||Rc0|| (50)[0116] 上式中，Ru0、Rp0和Rc0分别为位移场、压力场、裂缝相场收敛容差。
[0117] 本发明提供的一种裂缝性地层压裂裂缝延伸轨迹预测方法，在裂缝延伸后不需要重新剖分网格，不需要引用额外准则来追踪裂缝轨迹。
附图说明
[0118] 图1实施例1计算区域及边界条件示意图；
[0119] 图2实施例1注入8s后裂缝轨迹图；
[0120] 图3实施例1注入16s后裂缝轨迹图；
[0121] 图4实施例1注入24s后裂缝轨迹图。
具体实施方式
[0122] 下面结合某井参数对本发明做进一步的详细说明，但不构成对发明的任何限制，其中计算所采用的基本参数见表1。
[0123] 表1实施例1计算所采用基本参数表
[0124]
[0125]
[0126] 第一步：建立计算物理模型，设定计算区域边界条件，如图1所示，该实施例中，计算区域为一20m×20m的正方形区域，计算区域中心有一初始长度为2m的裂缝，压裂液从初始裂缝中心注入。在距离注入点正上方3m处有一条与y方向夹角为30°长度为5m的天然裂缝。在计算区域右边界施加沿x方向的最小水平地应力20MPa，在计算区域上边界施加沿y方向的最大水平地应力25MPa，计算区域左边界沿x方向的位移固定为0，计算区域下边界沿y方向的位移固定为0。
[0127] 第二步：网格剖分，将计算区域均匀的剖分为80×80个四节点正方形单元。
[0128] 第三步：将表1中的参数带入本发明所建立的方程组进行模拟。模拟结果如图2～图4所示。由图2可知在水力裂缝与天然裂缝相交前，水力裂缝沿直线扩展，从图3可知，在该实施例中当水力裂缝与天然裂缝相交后，水力裂缝将沿天然裂缝扩展。由图4可知，水力裂缝沿天然裂缝扩展到天然裂缝尖端后，将转向沿最大水平地应力方向扩展。
