一种基于离散元流固耦合的多簇压裂裂缝形态预测方法与流程

本发明涉及页岩压裂技术领域，具体涉及一种基于离散元流固耦合的多簇压裂裂缝形态预测方法。

背景技术：

页岩气存于以富有机质页岩为主的储集岩系中的非常规天然气，是连续生成的生物化学成因气、热成因气或二者的混合，可以游离态存在于天然裂缝和孔隙中，以吸附态存在于干酪根、黏土颗粒表面，还有极少量以溶解状态储存于干酪根和沥青质中；

随着勘探技术的突破和商业性开发程度的加深，目前页岩气储层开发已逐步向深层进军。据测算，全球深层页岩气资源量巨大，勘探开发前景十分广阔。深层页岩气开发是世界级难题，深层页岩压裂面临诸多难题：随井深增加，压裂液沿程摩阻增加，施工排量受限，缝内净压力低，造缝能力减弱，砂液比难以提升；温度及压力增加，岩石塑性特征增强，裂缝起裂延伸困难；闭合压力高，支撑剂嵌入严重，各种尺度裂缝窄，裂缝导流能力递减快，长期导流能力难以维持。因此，在深层页岩中形成高效支撑的裂缝网络，提高改造体积，关键是要认识深层页岩压裂缝网形成机理和导流能力的优化，从而提出针对性工艺技术方法。

技术实现要素：

针对上述问题，本发明提供一种基于离散元流固耦合的多簇压裂裂缝形态预测方法,基于离散元方法，考虑层理和天然裂缝的影响，分别通过沿裂缝高度方向和沿裂缝长度方向的裂缝扩展模型，模拟分析裂缝在缝高方向和缝长方向的扩展规律。

本发明采用下述的技术方案：

一种基于离散元流固耦合的多簇压裂裂缝形态预测方法，包括以下步骤：

s1、建立流动方程，根据二维圆盘颗粒模型建立流体压力的增量方程：

式中，kf为流体的体积模量，vd为流体空间的体积，δt为时步，∑qδt-δvd为流体变化引起流体空间内部体积的变化；

s2、建立流固耦合方程，根据所要达到的流体与固体的耦合作用建立流体的变化引起的流体空间内压力的变化方程：

式中，kf为流体的体积模量，vd为流体空间的体积，δt为时步；

s3、模拟垂向上层理对裂缝在高度上的扩展，进而明确层理对缝高的影响程度的主控因素，最后对压裂工艺和施工参数进行优化；

s4、模拟裂缝沿裂缝长度方向上的扩展，进而明确裂缝延伸的主控因素，最后对压裂工艺和施工参数进行优化。

优选的，所述二维圆盘颗粒模型为：

1)、将颗粒围成的域看作为一个流体空间，并假定空间内的所有参数保持不变；

2)、颗粒间的胶结处存在一个流体通道，负责与其他流体空间相连，可以让液体自由流入流出。

优选的，所述所要达到的流体与固体的耦合作用的条件为：

1)、接触力的变化应引起颗粒间流体通道的变化，或者控制流体通道的开启或关闭；

2)、流体空间内流体体积的变化造成区域压力的改变；

3)、流体空间内的压力作用在包围该空间的颗粒上。

优选的，所述层理对缝高的影响程度的主控因素包括垂向应力差、杨氏模量、抗张强度、层理密度、断裂韧性、排量。

优选的，所述裂缝延伸的主控因素包括天然裂缝密度、水平应力差、杨氏模量、抗张强度、断裂韧性、排量。

本发明的有益效果是：

本发明基于离散元方法，考虑层理和天然裂缝的影响，分别建立起沿裂缝高度方向和沿裂缝长度方向的裂缝扩展模型，模拟分析了裂缝在缝高方向和缝长方向的扩展规律，明确了裂缝延伸的主控因素；

本发明分析了水力裂缝扩展规律，层理密度、垂向应力差是裂缝沿高度方向上扩展的主要控制因素；杨氏模量、天然裂缝密度、水平应力差是裂缝沿长度方向上扩展的主要控制因素。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍，显而易见地，下面描述中的附图仅仅涉及本发明的一些实施例，而非对本发明的限制。

图1为本发明二维圆盘颗粒domain(域)的范围示意图；

图2为本发明domain(域)内部压力按比例分布在周围颗粒上的分布示意图；

图3为本发明裂缝纵向扩展模型示意图；

图4为本发明层理条数与缝高关系示意图；

图5为本发明垂向应力差与缝高关系示意图；

图6为本发明断裂韧性与缝高关系示意图；

图7为本发明杨氏模量与缝高关系示意图；

图8为本发明抗张强度与缝高关系示意图；

图9为本发明裂缝横向扩展模型示意图；

图10为本发明杨氏模量与缝宽关系示意图；

图11为本发明杨氏模量与半缝长关系示意图；

图12为本发明杨氏模量与改造体积关系示意图；

图13为本发明断裂韧性与缝宽关系示意图；

图14为本发明断裂韧性与半缝长关系示意图；

图15为本发明断裂韧性与改造体积关系示意图；

图16为本发明抗张强度与缝宽关系示意图；

图17为本发明抗张强度与半缝长关系示意图；

图18为本发明抗张强度与改造体积关系示意图；

图19为本发明裂缝密度与缝宽关系示意图；

图20为本发明裂缝密度与半缝长关系示意图；

图21为本发明裂缝密度与改造体积关系示意图；

图22为本发明应力差与半缝长关系示意图；

图23为本发明应力差与改造体积关系示意图；

图24为本发明排量与缝宽关系示意图；

图25为本发明排量与半缝长关系示意图；

图26为本发明排量与改造体积关系示意图；

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例的附图，对本发明实施例的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于所描述的本发明的实施例，本领域普通技术人员在无需创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

除非另外定义，本公开使用的技术术语或者科学术语应当为本公开所属领域内具有一般技能的人士所理解的通常意义。本公开中使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。“上”、“下”、“左”、“右”等仅用于表示相对位置关系，当被描述对象的绝对位置改变后，则该相对位置关系也可能相应地改变。

下面结合附图和实施例对本发明进一步说明。

如图1至图26所示，一种基于离散元流固耦合的多簇压裂裂缝形态预测方法，包括以下步骤：

s1、建立流动方程：如图1所示，图中的domain(域)看作为一个流体空间，假定空间内的所有参数保持不变。颗粒间的胶结处存在一个流体通道，负责与其他流体空间相连，可以让液体自由流入流出。假设孔径宽度为a，通道长度l为该胶结键的长度，根据立方定律，有：

式中，k代表渗透率，p2、p1为通道两侧流体空间的压力(空间2和空间1)，并假定从空间2流入空间1的方向为正方向。

当整个介质没有裂缝时，依然可以发生液体的流动，因此假定孔径在没有载荷作用时依然存在一个残余孔径a0当孔径两侧出现压应力时，该孔径会缩小，但如果让孔径完全闭合，需要的压应力为无限大；相反地，如果通道两侧为拉应力，则孔径会增大。基于该假设，适用于压应力情况下，则有：

式中，f0是给定值，相当于孔径处的初始压应力，f为孔径处新增的压应力值。

可以发现，当f＝f0时，a＝a0/2。

如果胶结键处是拉应力，则孔径等于残余孔径与胶结键两侧颗粒真实距离б之和，

a＝a0+mδ(3)

式中，m代表间距增大系数，常规取1。

如图2所示，每一个流体空间都会从与其他空间有接触的流体通道交换流体，总量记为∑q，在每个时步内，流体压力的增量可以通过下式计算，假定流体流入为正，

式中，kf为流体的体积模量，vd为流体空间的体积，δt为时步，∑qδt-δvd为流体变化引起流体空间内部体积的变化；

s2、建立流固耦合方程，

所要达到的流体与固体的耦合作用的条件为：

1、接触力的变化应引起颗粒间流体通道的变化，或者控制流体通道的开启或关闭；

2、流体空间内流体体积的变化造成区域压力的改变；

3、流体空间内的压力作用在包围该空间的颗粒上。

步骤s1已经完成前条件1、2，针对条件3，将作用颗粒上的力看作为施加在胶结键处的力，因为流体空间内部一定是多边形，有：

fi＝pnis(5)

式中，ni即代表多边形第i条边的方向矢量，s代表第i条边的边长，p代表流体空间内部的压力。为了保证算法的稳定性，需要对时步的大小进行限制。假设流体空间内所有流体通道两侧发生的压力改变均为δpp，可以计算出该情况下流体空间内流量的总变化值为：

式中，r为包围流体空间周围颗粒的平均粒径，2r即可表示通道长度l。n为流体通道的个数，根据该方程，我们可以计算流体的变化引起的流体空间内压力的变化，有

式中，kf为流体的体积模量，vd为流体空间的体积，δt为时步；

在上述的压力改变δpp下，产生的新压力变化δpr应小于前者，因此可以求得时步的上限为

s3、通过模拟垂向上层理对裂缝在高度上的扩展，进而明确层理对缝高的影响程度，最后对压裂工艺和施工参数进行优化；

如图3所示，根据裂缝纵向扩展模型，输入已统计的目标井区压裂井参数，包括岩石弹性模量，泊松比，抗张强度，断裂韧性，地应力，天然裂缝参数(密度，倾角，倾向，尺寸等)：

表1纵向扩展模拟输入参数表

1)、层理密度：分别模拟了层理条数为2、4、6、10和20条情况下缝高扩展，结果如图4所示；

模拟结果表明随着层理密度增大，层理数量增加，扩展过程中突破的层理数量增加。由于层理的开启和压裂液向层理的滤失，裂缝高度方向的扩展受到限制，缝高随着层理密度增加不断减小，但降幅也逐渐减小。

2)、垂向应力差：分别模拟了垂向应力差为5、10、15、20、25mpa下缝高扩展，结果如图5所示；

模拟结果表明随着垂向应力差增加，作用于层理面上的应力增加，层理开启难度增加，压裂液滤失减小。因此，裂缝高度随着垂向应力差的增加而增加，但是增幅较小。

3)、断裂韧性：分别模拟了断裂韧性为0.2、0.4、0.6、0.8、1mpa·m^0.5下缝高扩展，结果如图6所示；

模拟结果表明随着断裂韧性增加，裂缝扩展需要的能量增加，因此裂缝缝高随之减小，但影响较小。

4)、杨氏模量：分别模拟了杨氏模量为10、15、20、30、40gpa下缝高扩展的，结果如图7所示；

模拟结果表明随着杨氏模量增加，裂缝高度会略微上涨，但不会突破新的层理，表明杨氏模量对裂缝高度的影响较小。

抗张强度：分别模拟了抗张强度为4、6、8、10、12mpa下缝高扩展，结果如图8所示；

模拟结果表明随着排量增加，裂缝不断突破层理，缝高增加。

s4、模拟裂缝沿裂缝长度方向上的扩展，进而明确裂缝延伸的主控因素，最后对压裂工艺和施工参数进行优化。

如图9所示，根据裂缝横向扩展模型，输入已统计的目标井区压裂井参数，包括岩石弹性模量，泊松比，抗张强度，断裂韧性，地应力，天然裂缝参数(密度，倾角，倾向，尺寸等)，如表2所示：

表2横向模拟输入参数表

1)、杨氏模量：分别模拟了杨氏模量为10、20、30、40gpa下多裂缝扩展，结果如图10-图12所示：

模拟结果表明弹性模量增加，水力裂缝缝宽减小，缝长增加，改造体积减小；其中杨氏模量增加对缝宽的影响较为明显，达到30％。

2)、断裂韧性：分别模拟了断裂韧性为0.5、0.75、1、1.25、1.5mpa·m^0.5下多裂缝扩展，结果如图13-图15所示：

模拟结果表明断裂韧性增加水力裂缝缝宽增加、缝长和改造体积减小，但变化幅度较小；在区块断裂韧性的变化范围内，该参数对水力裂缝扩展的影响较小。

3)、抗张强度：分别模拟了抗张强度为2、6、10、14、18mpa下多裂缝扩展，结果如图16-图18所示：

模拟结果表明抗张强度增加，水力裂缝缝宽增加、缝长和改造体积减小，但变化幅度均较小；在目标井区，抗张强度的变化范围对水力裂缝扩展的影响很小。

4)、天然裂缝密度：分别模拟了天然裂缝密度为1、2、4、6、8、10条/m下多裂缝扩展，结果如图19-图21所示：

模拟结果表明天然裂缝密度增加，水力裂缝缝宽和缝长减小，改造体积增加。天然裂缝密度增加水力裂缝的复杂程度增加，但是裂缝控制的区域将会减小。

5)、水平应力差，分别模拟了应力差为5、10、15、20mpa下多裂缝扩展，结果如图22-图23所示：

模拟结果表明应力差增加，水力裂缝缝宽和缝长有一定增加，裂缝改造体积减小较为明显；应力差增加裂缝间的干扰相对减弱。

同时结合天然裂缝与水力裂缝的交互角度可知，增加应力差，天然裂缝面上的正应力增加，天然裂缝的开启和剪切破坏的难度均会增加。

6)、排量：分别模拟排量为8、10、12、14、16、18m³/min下多裂缝扩展，结果如图24-图26所示：

模拟结果表明当排量小于12m³/min以主缝为主，排量大于14m³/min后以主缝+分支缝为主；当排量大于14m³/min后，改造体积的增加幅度有限。

模拟应力差为10mpa和20mpa下的排量为8、10、12、14、16、18m³/min的多裂缝扩展，低排量下天然裂缝开启数量有限，增加排量有利于天然裂缝的开启，高应力差天然裂缝延伸受限制。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。