本发明涉及油气田增产改造领域，具体涉及一种煤砂互层穿层压裂射孔位置优化方法。
背景技术：
中国煤层气田含有丰富的天然气资源，该类气田除含煤层气外，在纵向上还含有丰富的致密砂岩气，形成煤砂岩互层，通过水力压裂技术在纵向上沟通煤岩和砂岩储层成为该类气田高效开发的一种方式。但由于层间界面以及不同层岩性和地应力差异的影响，若射孔位置不恰当则会使得水力裂缝只能在一个层里面延伸，因此通过人为优化射孔位置使得水力裂缝沟通多个储层具有重要意义。
而该项技术的核心是建立水力裂缝在多层中纵向延伸模型。常用的模拟裂缝在多层中纵向延伸的数值模拟方法包括:(1)palmer模型及其改进模型(2)位移不连续法(ddm)；(3)基于abaqus平台的有限元法；(4)扩展有限元法。但上述数值模拟方法在研究水力裂缝遭遇层间界面后的延伸轨迹时，需要建立相交准则来判断水力裂缝是穿过界面，还是开启界面。此外，现有裂缝延伸模型大多需要引入carter滤失系数来描述压裂液滤失现象。
技术实现要素：
针对上述技术问题，本发明通过综合应用biot多孔弹性理论、有限元理论、相场法理论、非线性方程数值求解方法等多学科知识，建立多孔弹性介质中水力裂缝纵向延伸数值计算模型，并基于此模型形成煤砂互层穿层压裂射孔位置优化方法。
一种煤砂互层穿层压裂射孔位置优化方法，包括以下步骤：
(1)收集输入参数；
(2)建立应力计算方程组；
(3)建立流体流动控制方程组；
(4)建立裂缝相场演化方程组；
(5)结合步骤(2)～(4)建立多孔弹性介质中水力裂缝纵向延伸数值计算模型；
(6)将步骤(1)获得的参数输入步骤(5)建立的模型，对比不同地层参数和射孔位置条件下裂缝轨迹，从而优化射孔位置。
进一步的，所述步骤(1)中收集的输入参数包括：地应力参数、不同层岩石弹性模量和泊松比、不同层岩石临界拉应力、不同层岩石渗透率、层间界面临界拉应力、层间界面渗透率、压裂排量、压裂液注入时间、压裂液粘度。
所述步骤(2)建立应力计算方程组，包括以下内容：
(2.1)应力平衡方程建立
多孔弹性岩石的应力平衡方程为
式中：σeff为有效应力，mpa，可通过公式(2)计算得到；α为biot系数；i为单位张量，二维情况下为[110]t；p为流体压力，mpa；
将方程(2)带入方程(1)中，应力平衡方程可重新写为：
式中：ψeff为储存于岩石骨架中的弹性应变能密度，mpa；ε为应变张量；g(c)为衰减函数，本发明定义衰减函数如公式(4)所示；c为裂缝相场，c＝1表示岩石完全破裂，c＝0表示岩石完好无损；ψ eff为拉伸弹性应变能密度，可通过公式(5)计算得到，mpa；ψ-eff为压缩弹性应变能密度，可通过公式(5)计算得到，mpa；σ eff为拉应力，mpa；σ-eff为压缩应力，mpa。
g(c)＝(1-c)2(4)
式中：λ和g为拉梅常数，mpa；εi(i＝1，2，3)为主应变；函数<x> ＝(|x| x)/2，<x>-＝(|x|-x)/2。
(2.2)应力平衡方程对应边界条件
上述应力平衡方程(3)需结合边界条件(6)才可以得到求解：
式中，
所述步骤(3)建立流体流动控制方程组，包括以下内容：
(3.1)多孔介质中流体流动连续性方程建立
多孔介质中流体流动连续性方程为：
式中：t是时间，s；ζ是流体体积含量的增量，可以通过方程(8)计算；v是流体流速，m/s，可通过方程(9)计算得到：
将公式(8)和(9)带入公式(7)，流体流动连续性方程重新写为：
式中：m为biot模量，mpa；εii岩石骨架体积应变；k为各向异性渗透率张量；μ为流体粘度，mpa.s。
(3.2)渗透率计算方程
对于二维情况，渗透张量可表示为：
其中
式中：kx和ky分别为x和y方向的渗透率，m2；k0为岩石基质初始渗透率，m2；wc是渗透率加权系数，代表了水力裂缝或层间界面对计算单元渗透率的贡献，本发明采用一种简单的渗透率加权系数计算公式，即wc＝w/he，w为裂缝宽度，由于在相场法中裂缝被转化分布于整个计算区域，因此需要计算所有单元的裂缝宽度，如公式(13)所示，he为有限元单元网格尺寸；kf为水力裂缝或层间界面渗透率，m2，可通过公式(14)计算得到；θ为裂缝面法向方向角或最大主应变方向角，可通过公式(15)计算得到。
w＝<ε1-εc> he(13)
式中：ε1为最大主应变；η为裂缝面形状参数；γxy为剪应变；εy为y方向应变；εc为岩石开始破裂时的临界拉应变，在相场法中，岩石开始破裂时的临界拉应变εc、临界拉应力σc和裂缝临界能量释放率gc三者之间可通过公式(16)联系起来：
式中：σc为临界拉应力，mpa；gc为裂缝临界能量释放率，mpa.m；l0为长度尺度参数，用于控制扩散裂缝区域宽度，m，如图1所示；e为岩石弹性模量，mpa.
(3.3)流体流动连续性方程对应的边界条件
上述流体流动连续性方程(10)需结合边界条件(17)才可以得到求解：
式中，
所述步骤(4)建立裂缝相场演化方程组，包括以下内容：
(4.1)尖锐裂缝的相场法近似
在相位场方法中，尖锐的裂缝γ通过一辅助相位场c(见图1)转化为扩散性裂缝γc(c)，而扩散性裂缝γc(c)可通过公式(18)描述：
其中，
(4.2)水力裂缝在多孔介质中延伸时自由能密度建立
水力裂缝在多孔介质中延伸时，多孔介质的总自由能密度ψ由储存于岩石骨架中的弹性应变能密度ψeff、储层于流体中的能量密度ψfluid、断裂能量密度ψfrac三部分组成，即
其中：
(4.3)基于变分原理建立裂缝相场演化方程
将方程(21)-(23)带入方程(20)可得到总自由能密度ψ的表达式。进而裂缝相场演化方程可通过ψ的变分导数确定，在与速率无关的条件下，演化准则可通过kuhn-tucker方程得到，即
则得裂缝相场演化方程为：
为了满足岩石损伤不可逆这一特性，相场演化方程(25)可改写为如下形式：
其中，h(ε，t)为整个过程中拉伸弹性应变能密度的最大值，可通过公式(27)计算得到：
(4.4)相场演化方程对应边界条件
相场演化方程(26)对应的边界条件如公式(28)所示：
式中：
进一步的，结合步骤(2)～(4)建立多孔弹性介质中水力裂缝纵向延伸数值计算模型；包括以下内容：
方程(3)、(10)、(26)组成了非线性方程组，本发明采用有限元法对该非线性方程组进行离散，采用向后欧拉法离散方程(10)中与时间有关的项。裂缝相场控制方程(26)采用picard迭代法求解，渗流-应力耦合方程组(3)和(10)采用newton–raphson(nr)迭代求解。在每一nr迭代步内，裂缝相场为固定值，则渗流-应力耦合方程组的nr迭代格式可写为：
式中：ru和rp分别为应力平衡方程和流体流动连续性方程的余量，如公式(30)和(31)所示；juu、jup、jpu和jpp为雅克比矩阵的四个分量，可通过公式(32)计算得到；δuh为第i个迭代步的位移增量，m；δph为第i个迭代步的压力增量，mpa。
公式(29)～(32)中，上标h代表计算网格节点处的值，下标n代表第n个时间步的值；δt为时间步长，s；nc、nu和np分别为裂缝相场、位移和压力有限元插值形函数；bu和buvol分别为应变矩阵和体积应变矩阵，bp为压力插值形函数的导数矩阵。
通过方程(29)求得第i个迭代步的位移增量δuh和压力增量δph后，第i 1个迭代步的位移和压力可表示为：
进而第i 1个迭代步的应变能历史函数h(ε，t)可求得，则第i 1个迭代步相场的试探解可通过方程(34)求得
其中
当位移、压力和裂缝相场都满足如公式(37)所示的收敛条件时，则迭代结束，进入下一时间步的计算，否则迭代继续进行。
||ru||≤tol||ru0||，||rp||≤tol||rp0||，||ci 1-ci||≤tol||rc0||(37)。
本发明采用相场法是基于变分原理的系统能量最小化建立起来，因此裂缝延伸轨迹是自动确定的。
本发明裂缝延伸轨迹和条件是自适应求解的，克服了现有技术中需要额外建立轨迹预测准则来判断裂缝延伸方向这一缺陷，且本发明不需要引入滤失系数来描述压裂液滤失现象。
附图说明
图1为实施例尖锐裂缝转化为扩散性裂缝以及边界条件示意图；
图2为实施例煤层射孔压裂计算区域及边界条件示意图；
图3为实施例煤层射孔压裂计算结束时裂缝相场分布图；
图4为实施例顶板砂岩射孔压裂计算区域及边界条件示意图；
图5为实施例顶板砂岩射孔压裂计算结束时裂缝相场分布图。
具体实施方式
下面结合我国山西某井对本发明做进一步的详细说明，但不构成对发明的任何限制，其中地层参数见表1。
表1实施例1计算所采用地层基本参数表
第一步：方案1：假设在煤层中心进行射孔压裂，计算区域如图2所示，计算区域被均匀的剖分为80×80个正方形单元，长度尺度参数l0为0.5m，压裂排量q为2.2×10-3m2/s，压裂液粘度为1mpa.s，注入时间为36s，模拟采用的时间步长为3s，并将表1中的参数带入本发明所建立的方程组进行模拟。模拟结果如图3所示，发现在这一地层条件下，在煤层中射孔压裂水力裂缝不能沟通上下砂岩储层。
第二步：方案2：由于在煤层中射孔压裂不能穿透上下砂岩层，因此改变射孔位置，在顶板砂岩中进行射孔压裂，计算区域如图4所示，将计算区域均匀的剖分为80×80个正方形单元，长度尺度参数l0为0.5m，压裂排量q为2.2×10-3m2/s，压裂液粘度为1mpa.s，注入时间为36s，模拟采用的时间步长为3s，并将表1中的参数带入本发明所建立的方程组进行模拟。模拟结果如图5所示，从图5可知，在顶板射孔压裂时，水力裂缝能够穿过上界面而进入煤层中继续延伸，但是当水力裂缝到达下界面时，水力裂缝将沿界面延伸，无法突破下界面而延伸到底板。
第三步：对比方案1和方案2可知，方案2开启的储层厚度大于方案1开启的储层厚度，因此选择方案2进行射孔。
1.一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，包括以下步骤：
(1)收集输入参数；
(2)建立应力计算方程组；
(3)建立流体流动控制方程组；
(4)建立裂缝相场演化方程组；
(5)结合步骤(2)～(4)建立多孔弹性介质中水力裂缝纵向延伸数值计算模型；
(6)将步骤(1)获得的参数输入步骤(5)建立的模型，对比不同地层参数和射孔位置条件下裂缝轨迹，从而优化射孔位置。
2.根据权利要求1所述的一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，所述步骤(1)中收集的输入参数包括：地应力参数、不同层岩石弹性模量和泊松比、不同层岩石临界拉应力、不同层岩石渗透率、层间界面临界拉应力、层间界面渗透率、压裂排量、压裂液注入时间、压裂液粘度。
3.根据权利要求1所述的一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，所述步骤(2)建立应力计算方程组，包括以下内容：
(2.1)应力平衡方程建立
多弹性岩石的应力平衡方程为：
式中：σeff为有效应力，mpa，通过公式(2)计算得到；α为biot系数；i为单位张量，二维情况下为[110]t；p为流体压力，mpa；
将方程(2)带入方程(1)中，应力平衡方程可重新写为：
式中：ψeff为储存于岩石骨架中的弹性应变能密度，mpa；ε为应变张量；g(c)为衰减函数，定义衰减函数如公式(4)所示；c为裂缝相场，c＝1表示岩石完全破裂，c＝0表示岩石完好无损；ψ eff为拉伸弹性应变能密度，可通过公式(5)计算得到，mpa；ψ-eff为压缩弹性应变能密度，通过公式(5)计算得到，mpa；σ eff为拉应力，mpa；σ-eff为压缩应力，mpa；
g(c)＝(1-c)2(4)
式中：λ和g为拉梅常数，mpa；εi(i＝1，2，3)为主应变；函数<x> ＝(|x| x)/2，<x>-＝(|x|-x)/2；
(2.2)应力平衡方程对应边界条件
上述应力平衡方程(3)需结合边界条件(6)才可以得到求解：
式中，
4.根据权利要求1所述的一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，所述步骤(3)建立流体流动控制方程组，包括以下内容：
(3.1)多孔介质中流体流动连续性方程建立
多孔介质中流体流动连续性方程为：
式中：t是时间，s；ζ是流体体积含量的增量，可以通过方程(8)计算；v是流体流速，m/s，可通过方程(9)计算得到：
将公式(8)和(9)带入公式(7)，流体流动连续性方程重新写为：
式中：m为biot模量，mpa；εii岩石骨架体积应变；k为各向异性渗透率张量；μ为流体粘度，mpa.s；
(3.2)渗透率计算方程
对于二维情况，渗透张量可表示为：
其中
式中：kx和ky分别为x和y方向的渗透率，m2；k0为岩石基质初始渗透率，m2；wc是渗透率加权系数，代表了水力裂缝或层间界面对计算单元渗透率的贡献，采用渗透率加权系数计算公式，即wc＝w/he，w为裂缝宽度，由于在相场法中裂缝被转化分布于整个计算区域，因此需要计算所有单元的裂缝宽度，如公式(13)所示，he为有限元单元网格尺寸；kf为水力裂缝或层间界面渗透率，m2，可通过公式(14)计算得到；θ为裂缝面法向方向角或最大主应变方向角，可通过公式(15)计算得到：
w＝<ε1-εc> he(13)
式中：ε1为最大主应变；η为裂缝面形状参数；γxy为剪应变；εy为y方向应变；εc为岩石开始破裂时的临界拉应变，在相场法中，岩石开始破裂时的临界拉应变εc、临界拉应力σc和裂缝临界能量释放率gc三者之间可通过公式(16)联系起来：
式中：σc为临界拉应力，mpa；gc为裂缝临界能量释放率，mpa.m；l0为长度尺度参数，用于控制扩散裂缝区域宽度，m；e为岩石弹性模量，mpa.
(3.3)流体流动连续性方程对应的边界条件
所述流体流动连续性方程(10)需结合边界条件(17)得到求解：
式中，
5.根据权利要求1所述的一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，所述步骤(4)建立裂缝相场演化方程组，包括以下内容：
(4.1)尖锐裂缝的相场法近似
在相位场方法中，尖锐的裂缝γ通过一辅助相位场c(转化为扩散性裂缝γc(c)，而扩散性裂缝γc(c)可通过公式(18)描述：
其中，
(4.2)水力裂缝在多孔介质中延伸时自由能密度建立
水力裂缝在多孔介质中延伸时，多孔介质的总自由能密度ψ由储存于岩石骨架中的弹性应变能密度ψeff、储层于流体中的能量密度ψfluid、断裂能量密度ψfrac三部分组成，即
其中：
(4.3)基于变分原理建立裂缝相场演化方程
将方程(21)-(23)带入方程(20)可得到总自由能密度ψ的表达式；进而裂缝相场演化方程可通过ψ的变分导数确定，在与速率无关的条件下，演化准则可通过kuhn-tucker方程得到，即：
则得裂缝相场演化方程为：
为了满足岩石损伤不可逆这一特性，相场演化方程(25)可改写为如下形式：
其中，h(ε，t)为整个过程中拉伸弹性应变能密度的最大值，可通过公式(27)计算得到：
(4.4)相场演化方程对应边界条件
相场演化方程(26)对应的边界条件如公式(28)所示：
式中：
6.根据权利要求1所述的一种煤砂互层穿层压裂射孔位置优化方法，其特征在于，结合步骤(2)～(4)建立多孔弹性介质中水力裂缝纵向延伸数值计算模型；包括以下内容：
方程(3)、(10)、(26)组成了非线性方程组，采用有限元法对该非线性方程组进行离散，采用向后欧拉法离散方程(10)中与时间有关的项；裂缝相场控制方程(26)采用picard迭代法求解，渗流-应力耦合方程组(3)和(10)采用newton–raphson(nr)迭代求解；在每一nr迭代步内，裂缝相场为固定值，则渗流-应力耦合方程组的nr迭代格式可写为：
式中：ru和rp分别为应力平衡方程和流体流动连续性方程的余量，如公式(30)和(31)所示；juu、jup、jpu和jpp为雅克比矩阵的四个分量，可通过公式(32)计算得到；δuh为第i个迭代步的位移增量，m；δph为第i个迭代步的压力增量，mpa；
公式(29)～(32)中，上标h代表计算网格节点处的值，下标n代表第n个时间步的值；δt为时间步长，s；nc、nu和np分别为裂缝相场、位移和压力有限元插值形函数；bu和buvol分别为应变矩阵和体积应变矩阵，bp为压力插值形函数的导数矩阵；
通过方程(29)求得第i个迭代步的位移增量δuh和压力增量δph后，第i 1个迭代步的位移和压力可表示为：
进而第i 1个迭代步的应变能历史函数h(ε，t)可求得，则第i 1个迭代步相场的试探解可通过方程(34)求得
其中
当位移、压力和裂缝相场都满足如公式(37)所示的收敛条件时，则迭代结束，进入下一时间步的计算，否则迭代继续进行；
||ru||≤tol||ru0||，||rp||≤tol||rp0||，||ci 1-ci||≤tol||rc0||(37)。
技术总结