CN109505576A - 页岩水力压裂三维全耦合离散裂缝网络模拟方法及系统 - Google Patents

Abstract

本发明公开了一种页岩水力压裂三维全耦合离散裂缝网络模拟方法及系统，包括：生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；对离散裂缝网络进行网格划分，划分为矩形的裂缝单元；建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；利用迭代法对流体流动的控制和本构方程和力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。采用本发明提出的方法隐式的耦合了流体流动与裂缝诱导应力之间的关系，模拟多条水力裂缝的扩展以及天然裂缝的张开和剪切作用。

Description

页岩水力压裂三维全耦合离散裂缝网络模拟方法及系统

技术领域

本发明属于油气田开发领域，更具体地，涉及一种页岩水力压裂三维全耦合离散裂缝网络模拟方法及系统。

背景技术

经典的水力压裂模型是注入液体产生二维水力裂缝，并在地层持续扩展(Economides和Nolte，2000)。但是现场监测表明，水力压裂可能形成多裂缝、分支网络和其他的一些复杂形态。在很多实际应用中，经典的水力压裂模型是可以胜任工程模拟目的的，但是在页岩油气水力压裂过程中裂缝网络可能变得非常复杂，经典模型显得过于简化而不再适用。水力裂缝在页岩岩石中的扩展会受已有天然裂缝、层理面和其他弱面的影响这些影响可能终止裂缝的扩展，形成分支和网络的复杂情况。剪切滑移是水力压裂过程中另一项重要作用，其注入的流体会诱发已有裂缝的滑动，通过破坏矿物组成结构和形成裂缝面的分离提高裂缝导流能力。

用于描述水力压裂复杂过程的模型目前正在开发，其过程很有挑战性，因为需要同时模拟流体流动、裂缝变形引起的各种诱导应力、支撑剂运移、流体的复杂流变性和裂缝网络的几何形态。基于连续介质的双重孔隙模型或者假定理想几何形状的网状裂缝已经被用于模拟压裂井长期生产情况，但是用其模拟水力压裂进程可能会过于简化。常规水力压裂模型假定只有一个平面裂缝，因此需要改进这些方法以描述非常规页岩资源的压裂。离散裂缝网络(DFN)模型显式的描述了各个裂缝，使其可以更加实际地模拟压裂进程，但是DFN方法实施起来有难度，要面临较高的计算强度以及裂缝网络产生的复杂几何形态问题。非常规水力压裂增产过程中的裂缝网络复杂模型模拟已经成了一些油气研究团队近些年来的主要工作，开发了一些非常规油藏的增产模拟新方法。

然而，现有的模拟水力压裂的方法或者是忽略了裂缝变形产生的诱导应力、或者应用的是简化的二维模型、或者是使用理想的局限于小量裂缝的网络模型，或者是使用非常理想化的半解析近似解。因此，有必要提供一种页岩水力压裂三维全耦合离散裂缝网络模拟方法，在大型复杂三维裂缝网络模拟中，完全隐式的耦合流体流动与裂缝变形(张开或者滑动)形成的诱导压力。

发明内容

本发明模拟页岩水力压裂复杂的三维离散裂缝网络，完全耦合了流体流动与裂缝诱导应力之间的关系，分析裂缝性页岩地层多裂缝的扩展以及主裂缝与天然裂缝的相互作用，解释页岩压裂常常观察到的高净值压力、相对短的裂缝长度和广泛的微地震事件，这些现象在经典水力压裂模型中无法预测，根据地质条件和施工参数的不同，可以模拟各种页岩压裂机理，从单条长水力裂缝到包括数条天然裂缝和主裂缝的大范围的裂缝网络，从而指导页岩压裂设计，提高压裂改造效果。

根据本发明的一方面，提出了一种页岩水力压裂三维全耦合离散裂缝网络模拟方法，该方法包括：

生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；

对所述离散裂缝网络进行网格划分，划分为矩形的裂缝单元；

建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；

利用迭代法对所述流体流动的控制和本构方程和所述力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

优选地，所述建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用包括：

应用非线性经验公式显示正应力、裂缝张开、裂缝滑动与裂缝开度和裂缝导流能力之间的关系；

应用半解析的一维液体滤失模型计算压力随时间在裂缝中变化；

应用Forch-heimer方程模拟裂缝内的非达西压降；

应用垂直剖面非均匀应力场控制页岩各层间的水力裂缝高度。

优选地，应用所述半解析的一维液体滤失模型，每个裂缝单元的滤失只与当前的裂缝单元有关，此时压力传导距离表示为：

其中，d_inv表示为压力传导距离；c_t为地层综合压缩系数；μ为流体粘度；φ_init为初始孔隙度；k为渗透率；t为时间。

优选地，基于所述Forch-heimer方程，建立包含非达西传导系数的达西方程，所述达西方程表示为：

其中，P为流体压力；x为沿裂缝距离；T_adj为非达西传导系数；μ为流体粘度；T为传导系数；Q为达西流量；B为无量纲常数；ρ为流体密度。

优选地，利用边界单元法建立所述力学计算的控制和本构方程，所述力学计算的控制和本构方程包括：

正应力的边界条件方程；

剪应力的边界条件方程；

在最大剪切应力方向上滑动的条件方程。

优选地，所述正应力的边界条件方程表示为：

其中，为裂缝远场载荷应力；P为流体压力；Δσ_n为裂缝单元累计变形形成的诱导应力引起的正应力的变化。

优选地，所述剪应力的边界条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；η为辐射阻尼系数；v为滑移速度

优选地，所述在最大剪切应力方向上滑动的条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；v_d为倾斜方向的剪切应力；v_s为走向剪切应力。

优选地，所述交叉作用准则包括：提前指定可能形成裂缝的位置，假设裂缝在开始模拟之前与最小主应力方向垂直。

根据本发明的另一方面，提出了一种页岩水力压裂三维全耦合离散裂缝网络模拟系统，该系统包括：

存储器，其上存储有计算机可执行指令；

处理器，所述处理器执行计算机可执行指令时实现以下步骤：

生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；

对所述离散裂缝网络进行网格划分，划分为矩形的裂缝单元；

建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；

利用迭代法对所述流体流动的控制和本构方程和所述力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

本发明的有益效果在于：开发的页岩水力压裂模拟器及其方法可以描述复杂三维离散网络裂缝形态，考虑流体流动与裂缝诱导应力之间的全耦合关系，分析多条主裂缝同时扩展间及主裂缝与天然裂缝的相互力学作用以及局部应力场的变化；根据地质条件和施工参数的不同，该方法可针对天然裂缝发育程度不同的均质或裂缝性页岩地层，模拟水力压裂过程，可显示从单条长水力裂缝到包括天然裂缝和主裂缝的大范围的密集裂缝网络，可以解释页岩地层现场压裂观察到的高的净值压力、相对短的裂缝长度和广泛的微地震事件，而这些现象在经典水力压裂模型中无法预测。

本发明的其它特征和优点将在随后具体实施方式部分予以详细说明。

附图说明

通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出了根据本发明的页岩水力压裂三维全耦合离散裂缝网络模拟方法的步骤的流程图。

图2示出了根据本发明的一个实施例的模拟计算中预设的复杂裂缝网络的示意性结构图。

图3示出了根据本发明的一个实施例的模拟形成的最终压裂裂缝网络侧上方的俯视图。

具体实施方式

下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

实施例1

在该实施例中，根据本发明的页岩水力压裂三维全耦合离散裂缝网络模拟方法可以包括：生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；对离散裂缝网络进行网格划分，划分为矩形的裂缝单元；建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；利用迭代法对流体流动的控制和本构方程和力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

该实施例模拟页岩水力压裂复杂的三维离散裂缝网络，完全耦合了流体流动与裂缝诱导应力之间的关系，分析裂缝性页岩地层多裂缝的扩展以及主裂缝与天然裂缝的相互作用，解释页岩压裂常常观察到的高净值压力、相对短的裂缝长度和广泛的微地震事件，这些现象在经典水力压裂模型中无法预测，根据地质条件和施工参数的不同，可以模拟各种页岩压裂机理，从单条长水力裂缝到包括数条天然裂缝和主裂缝的大范围的裂缝网络，从而指导页岩压裂设计，提高压裂改造效果。

图1示出了根据本发明的页岩水力压裂三维全耦合离散裂缝网络模拟方法的步骤的流程图。下面参考图1详细说明根据本发明的页岩水力压裂三维全耦合离散裂缝网络模拟方法的具体步骤。

步骤101，生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络。

具体地，裂缝能够被用户特别指定或者随机生成，裂缝被随机生成是对长度分配、定向和空间密度使用未指明的统计。可能形成的水力压裂的位置和方向必须被提前指定，可能形成的水力压裂裂缝位置被用户明确的指定。在2D空间里，选择使用一个运算法则，能够自动寻找可能在整个储层中形成的裂缝，但是这些不能在3D空间里实施。

当随机定位裂缝的时候，规定要执行几个检测程序以避免有问题的几何：低角度的裂缝交叉，非常紧密连接的裂缝裂缝交叉和接近平行的不相交裂缝。对这些进行检测的一个原因是为了避免非常小的单元或者高的纵横比。第二个原因是边界单元法，如Okada(1992)方法在一定距离范围内是不准确的，这是单元大小的一小部分。伴随着足够的网格细化，这种不准确性能够被避免，但是不利的几何裂缝要求一个不能接受的高度的网格细化来避免数值问题。

通过改变裂缝形态，目前的技术已经发展到能够去处理裂缝交叉处裂缝网络的复杂几何问题。然而，这些方法仅仅被应用于处理流动问题，尚不清楚它们能否被应用于处理涉及地应力耦合的问题。

步骤102，对离散裂缝网络进行网格划分，划分为矩形的裂缝单元。

具体地，在一些复杂的三维离散裂缝网模拟中，网格划分是一项挑战性的任务，本发明实施了减少网格生成的复杂性和提高模拟效率的各种方法。首先，一种边界单元技术被用于计算裂缝变形产生的诱导应力大小。在离散的裂缝网络中，边界单元技术允许应力被一个仅包括裂缝的问题网格计算，而划分裂缝周围的网格并不是很重要；第二，半解析一维滤失模型用于计算流体滤失到基质中，这种方法假设滤失来自每个单独的一维的裂缝单元，且每个裂缝单元中的滤失不受邻近裂缝单元的影响。只要基质渗透率足够低，这种方法是非常精确的，且这种模拟只持续几个小时或者几天就可完成。

本发明将裂缝网格划分为矩形单元，对于压裂问题，这种形状的矩形单元对计算模拟的限制较小。因为裂缝单元必须是矩形的而且划分的网格具有相似性，所以在模拟中所有的裂缝必须是垂直的。

为了进一步简化均布网格划分的过程，裂缝高度被略微的调整以确保它们能够适应预先设定的网格。假如目标区域是h米高、N_elm个单元，则每个裂缝单元的高度等于为h/N_elm，通过调整，裂缝的顶部和底部是在z坐标位置，这样(h/2-z)/N_elm才是一个整数值(目标区域的中心位于z＝0处)。

裂缝单元长度在水平方向上是一个常量，l_elm。因为裂缝交叉，有时裂缝长度不能保证等于h_elm，所以裂缝长度允许稍长或稍短。

步骤103，建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用。

在一个示例中，建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用包括：应用非线性经验公式显示正应力、裂缝张开、裂缝滑动与裂缝开度和裂缝导流能力之间的关系；应用半解析的一维液体滤失模型计算压力随时间在裂缝中变化；应用Forch-heimer方程模拟裂缝内的非达西压降；应用垂直剖面非均匀应力场控制页岩各层间的水力裂缝高度。

具体地，应用有限体积法全隐式地求解流体流动方程，其时步长是自适应的选择过程。

流体流动是等温的单相流动(不包括支撑剂)，忽略重力的影响。流体具有恒定的粘度，是轻微的可压缩的，压缩系数是C_f，初始密度是ρ_init。

在裂缝内，应用达西定律求解非稳态质量平衡方程：

其中，E为孔隙开度(单位面积的裂缝存储的流体体积)；t为时间；q_flux为流体在裂缝间的质量流量，ρ为流体密度，k为渗透率，μ为流体粘度，P为流体压力；e为水力开度(裂缝允许流体流动的有效厚度)；q_leakoff为单位裂缝表面积滤失进周围基岩的质量流量；s为井的源项。

裂缝渗透率定义式为：

其中，e为水力开度。

当渗透率乘以水力裂缝开度，其结果就是传导系数T，经典的三次方定律：

其中，e为水力开度。

非线性关系是用来确定压力、应力和变形与空隙开度和裂缝开度之间的关系。发明中定义的闭合的裂缝表示其中的流体压力小于正应力且裂缝的两个壁面已相互接触。粗糙的裂缝壁允许闭合的裂缝保持一定的水力和空隙开度。定义了一个开启的裂缝，其中流体压力达到了正应力，裂缝壁面不再接触。对于天然裂缝，空隙开度和裂缝开度允许不同。这种方法的目的是包括微观和宏观尺度下的裂缝粗糙度的影响。

为了计算闭合裂缝单元的开度，使用下面的节点闭合的关系式：

其中，E₀为在有效应力为0时的空隙开度；δ'_n为正应力；δ_n,Eref为使空隙减小90％的有效正应力，有效的正应力δ_n,Eref等于正应力σ_n′减去流体压力P。

不同常数间的相似关系可用于计算裂缝开度e：

其中，e₀为在有效应力为0时的裂缝开度；δ'_n为正应力；δ_n,Eref为有效正应力；D_eff为累计滑移位移；φ_edil为剪切膨胀角，其可以关联剪切变形和裂缝传导性的增加。

D_eff的取值满足：

其中，D为某位置的累计滑移位移；D_eff,max为最大累计滑移位移。

如果裂缝微元是打开的，那么它的空隙开度和裂缝开度分别等于E₀和e₀，加上因为力学作用产生开度E_open。力学张开是在开裂的过程中发生在裂缝面之间的分离的量，空隙开度E表示为：

E＝E₀+E_open， (13)

其中，E₀为在有效应力为0时的空隙开度；E_open为力学作用产生开度。

裂缝开度e表示为：

e＝e₀+D_efftan(φ_edil)+E_open， (14)

其中，e₀为在有效应力为0时的裂缝开度；D_eff为累计滑移位移；φ_edil为剪切膨胀角；E_open为力学作用产生开度。

对于天然裂缝，E₀和e₀，是常数，但是对于水力裂缝单元，E₀和e₀值允许在模拟期间变化(每个单元的值都不同)。

算术平均用来计算两个不同传导率单元之间的平均传导率。在DNF模拟中，调和平均经常用于传导率。传导率是流体压力的一个高度非线性函数，函数平均可以创建一个强大的网格划分依赖性，只有高度的网格细化才能得到解决。计算过程中井筒储集效应可选择性的被包含，井筒储集系数可表示为C_w。

定义新形成的水力裂缝单元开度有很多的困难。在本方法中使用了不同的算法，它对于实现在测试压裂过程中压力瞬变行为模拟有实际的作用。

对于天然裂缝单元，E₀和e₀的值保持恒定。而对于水力裂缝单元E₀和e₀是不同的，它们被被允许增加直到达到一个极限值。被用于计算水力裂缝开度E₀的值被称为残余开度。

当水力裂缝单元起裂的时候，它被给予了一个非常小的残余开度，大约10微米。这种算法不能很好的保持质量平衡，因为一个新的单元起裂后。一个相应数量的流体不能在相邻单元被移除，但这种误差是非常小的，因为新的单元起裂时开度很小。本发明定义了一个残余开度参数，如果空隙和RAP的比残余开度大，那么残余开度会被增加到等同于空隙体积。当单元的残余开度增加，单元的全部空隙开度将保持不变。残余开度被允许增加到一个可能的最大值，这时可选择应用于天然裂缝单元的E₀的值。如果一个水力裂缝单元闭合，它的残余开度将保持恒定。在本方法中，我们使用50％的残余开度。

对于水力裂缝单元(不像天然裂缝)，裂缝开度e被设定为等于空孔开度E。这相当于设置e₀为E₀和设置为0。这种运算法则用于模仿自然形成的裂缝粗糙度，裂缝最初具有一个非常小的开度，然后断持续变化直到裂缝完全张开。

在一个示例中，应用半解析的一维液体滤失模型，每个裂缝单元的滤失只与当前的裂缝单元有关，此时压力传导距离表示为：

其中，d_inv表示为压力传导距离；c_t为地层综合压缩系数；μ为流体粘度；φ_init为初始孔隙度；k为渗透率；t为时间。

具体地，应用一维滤失模型，每个裂缝单元的滤失不被其他单元所影响，这种模型忽略了裂缝变形引起的地层流体压力的孔弹性作用。该方法的优点在于它给了一个高精度和有效的1D扩散方程求解方法，甚至包括裂缝中压力的任意变化。相反的，卡特滤失模型假设在裂缝中压力恒定，而这个简单地假定，严重的降低了模型的适用性。

相对于压力传导距离，如果裂缝之间间隔足够大、裂缝之间没有干扰，1D滤失的假设是合理的，压力传导距离可以通过公式(1)估算，通过公式(1)可知，在一个小时的模拟时间中，压力传导距离大约是3.2米，比模拟使用的典型的裂缝网络裂缝间距要小得多。

在一个示例中，基于Forch-heimer方程，建立包含非达西传导系数的达西方程，达西方程表示为：

其中，P为流体压力；x为沿裂缝距离；T_adj为非达西传导系数；μ为流体粘度；T为传导系数；Q为达西流量；B为无量纲常数；ρ为流体密度。

具体地，根据Forchheimer方程，一维裂缝中非达西流动可以表示为：

其中，P为流体压力；x为沿裂缝距离；μ为流体粘度；e为裂缝开度；Q为达西流量；B为无量纲常数；ρ为流体密度；T为传导系数。

根据Holditch和Morse(1976)，方程10可以用无量纲的非达西传导系数T_adj重写为达西方程的形式：

其中，P为流体压力；x为沿裂缝距离；T_adj为非达西传导系数；μ为流体粘度；T为传导系数；Q为达西流量；B为无量纲常数；ρ为流体密度。

在本方法中，非达西传导系数T_adj是唯一的能够显式计算的参数。在每次时间步的开始，通过单元密度、界面长度或两个单元的平均空隙开度分配质量流量以计算邻近单元间的达西流动速率，之后T_adj被计算而且在相应的时间步中作为一个常数。

除了求解流体流动方程，该方法还计算由于裂缝张开和滑动引起的应力变化。应力计算应用边界元法，假定弹性均匀和各向同性的地层、线弹性变形和小的应变。这种方法满足经典的连续介质方程：准静态应力平衡、相容方程、胡克定律。计算结果收敛于裂缝变形的解析解。这种计算是全数值化的，能够解决复杂的问题。边界元方法的一个优点是，它仅仅在裂缝周围进行网格划分。本发明假设断裂变形不影响基质流体压力而且基质流体压力的变化不影响裂缝面上的正应力。

每个单元的应力状态由三部分组成：n、s和d，分别代表正应力、走向剪切应力和倾斜方向的剪切应力。每种单元上的变形由三种不同的变化组成的，E_open，v_s和v_d。这些是力学的开启、走向滑动速度和倾斜的滑动速度。因为有三种类型的变形，所以对于每种裂缝单元来说有三种应力组成，九乘九的矩阵被运用。实际计算中，不是在每个时间步中每种都会变形，而且在每个时间步不是每种单元都会变形。所以在每个时间步中执行完全的增量是不必要的，与不相关的单元的相互作用系数矩阵列的变形可以忽略。

在一个示例中，利用边界单元法建立力学计算的控制和本构方程，力学计算的控制和本构方程包括：正应力的边界条件方程；剪应力的边界条件方程；在最大剪切应力方向上滑动的条件方程。

具体地，对于一个张开的裂缝微元，应力边界条件是有效正应力必须等于0。

在一个示例中，正应力的边界条件方程表示为：

其中，为裂缝远场载荷应力，是由原始条件计算得到，是系统的所有的裂缝微元变形累加的诱导应力产生的正应力的变化；P为流体压力；Δσ_n为裂缝单元累计变形形成的诱导应力引起的正应力的变化。

裂缝单元也允许剪切变形。因为裂缝是一个二维的面，有两个滑动速度分量，V_D和V_S，即滑动速度在倾滑方向(向下为正)，走滑方向(方向为正。，如果剪应力的大小超过了裂缝的摩擦阻力则裂缝产生滑动，使用滑动速度的大小变量v、从滑动的逆时针方向测量的滑动方向θ，可以使问题很方便的被参数化，走滑和倾角滑动速度分量可以由v和θ很容易地计算出。

对于每一个滑动单元，下面所描述的2个方程，必须满足于(1)应力边界条件和(2)滑动方向。对于那些不滑动的单元，这些方程从方程组省略。

对于张开的单元，应力边界条件是剪切应力的大小(包括应力的动态影响，剪切应力，表示与辐射阻尼项)必须等于零。

在一个示例中，剪应力的边界条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；η为辐射阻尼系数；v为滑移速度。

辐射阻尼系数η等于剪切模量除以2倍的剪切波速，这个项解释了在高滑动速度惯性的影响。辐射阻尼通常对模拟结果有一个可以忽略的影响，但它在数值计算中非常有用，因为在时间步接近零持续时，它保证了一个时间步的变形幅度接近零。

对于滑动闭合的单元，应用库伦法则平衡方程：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；η为辐射阻尼系数；v为滑移速度；μ_f为摩擦系数；σ_n为单元正应力；P为流体压力；S₀为内聚力。

剪切应力小于摩擦阻力的裂缝滑动(方程16的右边)不滑动。使S₀等于S_o，open可行，否则，由于从闭合到张开造成内聚力的突然消失，裂缝单元会经历快速滑动。

对于每个滑动单元，有一个滑动应力边界条件，如方程4或16，第二个方程需要控制滑动的方向。这个方程指定的滑动方向必须等于最大剪应力的方向。

在一个示例中，在最大剪切应力方向上滑动的条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；v_d为倾斜方向的剪切应力；v_s为走向剪切应力。

方程3、4、16和5都是隐式地求解。Δσ_n、Δτ_s和Δτ_d的值都是在每次更新非线性动力学迭代时间步中进行评估。

步骤104，利用迭代法对流体流动的控制和本构方程和力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

具体地，使用迭代法求解非线性系统方程，迭代耦合发生在流体流动和正应力边界条件代表的方程组(记为系统方程组1)及剪切滑移方程组(记为系统方程组2)之间。在迭代耦合方案中，首先对系统方程组1进行了设置得以求解，系统方程组2未知量依然不变，然后，系统方程组2求解，而系统方程组1的未知量不变，重复过程直到所有方程同时满足条件。

在每个方程体系中，都应用了迭代方案。形成相应迭代矩阵J，第i行和第j列由如下给出：

其中，R为剩余方程矢量(如果方程完全满足，R＝0)；X是未知矢量。

之后解出如下面公式(18)的线性方程组：

JdX＝-R， (18)

其中，dX为未知矢量的修正。

在dX计算之后，未知矢量X更改如下：

X_j+1+1＝X_j+dX_j， (19)

其中，下标j为迭代数字。

原始参数变更后，所有的次级参数，比如流体性质，也要更改。之后剩余矢量被重新计算。继续迭代直到R和dX的欧几里德和无穷大标准低于收敛标准。当评估剩余向量的范数时，对应于一个质量平衡方程的值(单位为kg/s)，按单位压力乘以时步时间除以流体密度、单元面积、及与压力有关的空隙开度导数的乘积。

流体流动和张开方程体系形成四种导数：和压力有关的质量平衡方程导数(方程6)，和张开位移有关的质量平衡方程导数(方程6)，和压力有关的实际正应力方程(方程3)的导数，和张开位移有关的实际正应力方程(方程3)的导数。

对微元裂缝的张开位移的导数全部非零，因为在边界单元法里，每个单元的变形诱导应力和其他单元的变形线性相关。这是有问题的，因为这意味着表示正向位移导数的迭代矩阵排列可能会紧密，会显著的降低计算效率。

为了避免这个问题，方法在在模拟的开始进行了预处理步骤，并且用量纲法分析交互作用因子得到了一个确定的最小限度。所有的剩余交互因子(大多数因子)在迭代矩阵里被设置为0，确保迭代矩阵是稀疏的。这个策略并没有降低方案的准确性，因为每次更改后，每个裂缝单元的应力会使用Hmmvp完全计算。因为裂缝单元间的交互作用因子随着距离迅速衰竭，迭代方法迅速收敛。这里有一个要求，E_open不能是负值，并且在每次的迭代步骤更改前，要进行检查，避免某个更改使得E_open变成负值。

方程的平移体系包含方程4、16、5。如果裂缝单元的剪切应力比滑动摩擦阻力小，就不发生滑动，不必在方程体系中描述该判断。方程体系的参数有变量v(滑动速度值)和θ(滑动方向角，从逆时针到冲击方向测得)。随着流体流动或者方程的张开体系，关联着有量纲交互因子的迭代矩阵的项低于某确定界限就被设置为0。迭代矩阵里的非零值必须要考虑到在方程14中滑动变形可以改变剪切应力和正应力。

参数θ被定义在-π和π之间。更正可能会使超出该范围，所以每次更改后，要进行检查，如果必要，要加上或者减去周期性的2π使得θ回到正常范围。参数v被定义为非负值，所以每次更改后都要进行检查，如果v比0小，要重新设置为0。这些检查需要在针对裂缝滑动变形做出修正之前。

裂缝单元状态被定义为是否张开或闭合以及滑动或者不滑动。在进行隐式计算时的一个难点是，裂缝单元状态可能随着时间步长变化。为了控制这个问题，每次迭代后要进行裂缝状态检查。如果一个裂缝单元的状态发生改变，就需要在一两个方程系统中添加或者移除。因此时间步长期间求解的方程体系会持续变化，并且在该时间步长开始时无法得知。正在变化的裂缝单元状态造成了导数解的不连续，但是很少会造成不收敛，因为：(1)使用相对小的时间步长，通常以秒为单位(2)每次时间步长内，迭代解都会使用。很少出现不收敛的情况，模拟器会放弃错误的时间步长，并用减少的时间步长重复计算。

在一个示例中，交叉作用准则包括：提前指定可能形成裂缝的位置，假设裂缝在开始模拟之前与最小主应力方向垂直。

具体地，设定一个交叉作用准则用来控制水力压裂裂缝的起裂、扩展和穿越。其限制条件是，可能形成裂缝的位置必须要提前指定，裂缝在开始模拟之前假设和最小主应力方向垂直。这是个简化的假设，因为裂缝变形产生的诱导应力叠加可能导致最小主应力局部旋转，使得水力压裂裂缝在不同方向形成。

和井眼连接的水力压裂单元被称为“井眼起裂”裂缝单元。在模拟的开始，它们是非激活状态(因为水力裂缝还没出现)。在每个时间步长，规定对每个裂缝单元进行下列检查：是否裂缝微元是激活状态的，它的压力等于井眼压力；是否该单元处于拉伸状态，有效应力正是负值；是否拉力会超过岩石的抗拉强度T_str。如果是，水力裂缝在该单元处形成。压裂开始后，连接井眼的裂缝单元被激活，为了避免对网格依赖性，所有离井眼确定距离(通常是5米)的裂缝单元也被激活。

连接天然裂缝的水力压裂裂缝单元是压裂起裂的单元。与井眼起裂单元相似，要进行如下检查：是否裂缝单元是活性的，是否有和链接的天然裂缝单元相同的压力，它是否处于拉伸状态，是否拉力会超过岩石的抗拉强度T_str，如果是，该单元和一定距离内的单元都要起裂。

裂缝扩展按线弹性断裂机理进行，所有和非活性水力裂缝单元毗邻的活性水力裂缝单元被当作裂缝终端，所有这些裂缝单元的应力强度因子K₁表示为：

其中，G为剪切模量；v_p为是泊松比；E_open为因力学作用产生的开度。

对于一个二维模拟，a被定义为裂缝单元的半长。对于我们的三维模拟，我们定义a为单元半长和半高中的较小值(通常它们非常接近，因为我们使用的裂缝单元接近正方形)。如果应力强度因子比裂缝韧性K_Ic大，裂缝从该单元处扩展。通过激活任何直接毗邻应力强度因子到达裂缝韧性的单元的水力压裂单元完成扩展。

为了判断水力裂缝是否穿过天然裂缝，提供了三种选择。一种选择是允许水力裂缝直接穿越天然裂缝。第二种选择是水力裂缝与天然裂缝相交时停止在原方向的扩展。如果发生了水力裂缝扩展终止，上述描述的裂缝扩展规则不允许在天然裂缝另一侧起裂水力压裂。第三种可用选择的标准基于剪切应力和正应力计算，既随着水力压裂靠近天然裂缝，判断天然裂缝是否在水力裂缝到达前发生滑动(会阻碍裂缝扩展)。当使用这个标准时，预测穿越仅仅基于远场载荷边界条件和裂缝方向，如果符合了终止条件，水力裂缝不会穿越天然裂缝。

本发明建立了一种更加有效的模拟方法，在大型复杂三维裂缝网络模拟中，完全隐式的耦合流体流动与裂缝变形(张开或者滑动)形成的诱导压力。使用这种方法，大型三维裂缝网格模拟可以在几个小时或内完成。这种水力压裂复杂三维裂缝模拟方法针对裂缝张开和滑动应用了合理的应力边界条件，在力学计算时使用了很有效的边界单元法，其加密网格后得到收敛解与经典文献中的解析解一致。

发明中两个关键点包括使用边界单元法进行力学计算以及应用半解析的一维滤失模型。这就产生了以下的限制条件：(1)力学计算假定均质线弹性变形；(2)滤失计算在低渗地层短时间范围内相对精确(几个小时到几天)。但是这些假定允许裂缝周围不进行网格加密，这极大的简化了计算过程，提高了计算效率。

本发明不涉及到水力裂缝在任意方向的起裂和扩展，可能形成的水力裂缝的位置必须预先给定。实际上，通过大量的潜在裂缝被定义，可以提高裂缝模拟的灵活性，确定裂缝形成的的具体位置和方向。

本实施例开发的页岩水力压裂模拟器及其方法可以对油田现场水力压裂范围内含有数千条裂缝的离散裂缝网络形成过程进行模拟，且模拟结果非常有效，其描述了复杂三维离散网络裂缝形态，考虑流体流动与裂缝诱导应力之间的全耦合关系，分析多条主裂缝同时扩展间及主裂缝与天然裂缝的相互力学作用以及局部应力场的变化；根据地质条件和施工参数的不同，该方法可针对天然裂缝发育程度不同的均质或裂缝性页岩地层，模拟水力压裂过程，可显示从单条长水力裂缝到包括天然裂缝和主裂缝的大范围的密集裂缝网络，可以解释页岩地层现场压裂观察到的高的净值压力、相对短的裂缝长度和广泛的微地震事件，而这些现象在经典水力压裂模型中无法预测。

应用示例

为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。

中石化涪陵页岩某水平井段分15段压裂，每段3簇，杨氏模量,36GPa，泊松比,0.23，水平主应力差12MPa，在某一段沿着水平井射孔簇中按14m³/min流量注入一个小时。

图2示出了在上述应用示例中模拟计算中预设的复杂裂缝网络的示意性结构图。图2中心的白色区域的出现是由观察角度的原因造成的，所有裂缝都是垂直的，主裂缝垂直与最小主应力方向，天然裂缝走向随机。

首先，设定裂缝性页岩地层复杂裂缝网络，包括主裂缝方向，天然裂缝的随机走向等，建立离散裂缝网络，并对离散裂缝网络进行网格划分，划分为矩形的裂缝单元。如图2所示，预设的裂缝性页岩地层复杂裂缝网含有大约3200个天然裂缝，其随机走向、90度倾角和随机定位，模拟的范围有800长(最大水平应力方向，X轴方向)，400米宽(最小水平应力方向，Y轴方向)和200米高，Z轴方向。天然裂缝高度从30米到50米，长度从30米到80米(所有裂缝都是长方形)，设定每间隔25米可能形成的水力裂缝，但每个水力裂缝初始高度仅有20米，这些有限高度的水力裂缝被垂直叠加，以致填满整个问题区域。页岩各层间的水力裂缝高度可以用垂直剖面非均匀应力场进行控制或模型直接假设水力压裂高度。

然后，建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；最后，利用迭代法对流体流动的控制和本构方程和力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越，模拟页岩水力压裂最终形成的裂缝网络。

模拟计算区域的边缘被设定为无流动边界，地层诱导应力的计算是假定裂缝被镶嵌在弹性的半空间(地球表面作为半空间界面)，在计算过程中，假定裂缝深度足够大，地球表面对计算结果没有影响。模拟的地层性质，流体压力，地应力状况基本上按照Barnett页岩所具有的参数进行。因为裂缝都是90度倾角，垂直主应力对裂缝应力没有影响，不必考虑。

复杂三维全耦合网络裂缝裂缝模型在每个的时间步长中，每个单元上都需要计算四个基本变量，包括压力、张开位移、滑动速度和滑动方向。在每个单元上相关的四个方程是：非稳态质量平衡方程(公式6)，包含正应力的边界条件方程(公式3)，包含剪应力的边界条件方程(公式4)和一个在最大切应力方向上滑动的条件方程(公式5)。

所有的方程都会被隐式进行求解，这意味着在每个时间步中，随着耦合方程的求解，依赖于基本变量变化(例如压力或者流体性质)的各种属性值会被重新确定，唯一的例外是一个非达西传导系数调整因子，它会被采用显示的方法进行计算。

图3示出了在上述应用示例中模拟形成的最终压裂裂缝网络侧上方的俯视图。为了更好的显示裂缝网络，压力小于20MPa的裂缝被忽略了。从图3中可以看出，网络裂缝中的净压力比单条裂缝要高。作为一个模型假设，水力裂缝高度限制在20米，所以当水力裂缝遇到天然裂缝时，它们不能在天然裂缝周围扩展，水力裂缝为了在地层中扩展更远，流体压力需提高到足够高从而张开阻碍裂缝扩展的天然裂缝。流体压力需要升高的原因是因为天然裂缝并不垂直于最小主应力。

在低天然裂缝开度和高应力各向异性情况下，天然裂缝中具有最小的导流能力。流体流动路径总是选择阻力最小的方向，张开相对与最小主应力方向垂直的天然裂缝，而避免张开垂直最大主应力方向的天然裂缝，所以，这种地层条件可能造成相对稀疏的裂缝网络。

对因为e₀更高既天然裂缝开度较大的情况下，压裂的裂缝网络会更密集，有更多的分支缝，并且较多液体滤失进入到水力裂缝周围的天然裂缝中，有着各种方向的天然裂缝能张开和接受流体，允许大量水力裂缝形成，组成一个复杂的裂缝网络。

本应用示例开发的页岩水力压裂模拟器及其方法可以描述复杂三维离散网络裂缝形态，考虑流体流动与裂缝诱导应力之间的全耦合关系，分析多条主裂缝同时扩展间及主裂缝与天然裂缝的相互力学作用以及局部应力场的变化；根据地质条件和施工参数的不同，该方法可针对天然裂缝发育程度不同的均质或裂缝性页岩地层，模拟水力压裂过程，可显示从单条长水力裂缝到包括天然裂缝和主裂缝的大范围的密集裂缝网络，可以解释页岩地层现场压裂观察到的高的净值压力、相对短的裂缝长度和广泛的微地震事件，而这些现象在经典水力压裂模型中无法预测。

本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。

实施例2

根据本发明的实施例，提供了一种页岩水力压裂三维全耦合离散裂缝网络模拟系统，该系统包括：

存储器，其上存储有计算机可执行指令；

处理器，处理器执行计算机可执行指令时实现以下步骤：

生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；

对离散裂缝网络进行网格划分，划分为矩形的裂缝单元；

建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；

利用迭代法对流体流动的控制和本构方程和力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

此系统的限制条件是所有裂缝必须是垂直的，力学计算中假设的是在线弹性和均匀介质中进行，不包括支撑剂运移，潜在的可能形成的水力裂缝的位置必须提前指定。为验证模型的可靠性，模拟首先对一条水力裂缝在考虑/不考虑滤失的情况下的扩展进行模拟，并与经典的水力压裂解析解进行了对比。现场水力压裂模拟是在一个天然裂缝密集发育的裂缝性地层进行，模拟过程显示了天然裂缝与主裂缝间相互作用，解释了常常在页岩地层压裂观察到的高的净值压力、相对短的裂缝长度和广泛的微地震事件，这些现象在经典水力压裂模型中无法预测。根据地质条件和施工参数的不同，模型可以预测各种压裂现象，从最小滤失的长水力裂缝到包括天然裂缝和主裂缝的大范围的密集裂缝网络。

本实施例开发的页岩水力压裂模拟器及其方法可以描述复杂三维离散网络裂缝形态，考虑流体流动与裂缝诱导应力之间的全耦合关系，分析多条主裂缝同时扩展间及主裂缝与天然裂缝的相互力学作用以及局部应力场的变化；根据地质条件和施工参数的不同，该方法可针对天然裂缝发育程度不同的均质或裂缝性页岩地层，模拟水力压裂过程，可显示从单条长水力裂缝到包括天然裂缝和主裂缝的大范围的密集裂缝网络，可以解释页岩地层现场压裂观察到的高的净值压力、相对短的裂缝长度和广泛的微地震事件，而这些现象在经典水力压裂模型中无法预测。

本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (10)

1.一种页岩水力压裂三维全耦合离散裂缝网络模拟方法，其特征在于，该方法包括：

生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；

对所述离散裂缝网络进行网格划分，划分为矩形的裂缝单元；

建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；

利用迭代法对所述流体流动的控制和本构方程和所述力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。

2.根据权利要求1所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，所述建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用包括：

应用非线性经验公式显示正应力、裂缝张开、裂缝滑动与裂缝开度和裂缝导流能力之间的关系；

应用半解析的一维液体滤失模型计算压力随时间在裂缝中变化；

应用Forch-heimer方程模拟裂缝内的非达西压降；

应用垂直剖面非均匀应力场控制页岩各层间的水力裂缝高度。

3.根据权利要求2所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，应用所述半解析的一维液体滤失模型，每个裂缝单元的滤失只与当前的裂缝单元有关，此时压力传导距离表示为：

其中，d_inv表示为压力传导距离；c_t为地层综合压缩系数；μ为流体粘度；φ_init为初始孔隙度；k为渗透率；t为时间。

4.根据权利要求2所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，基于所述Forch-heimer方程，建立包含非达西传导系数的达西方程，所述达西方程表示为：

其中，P为流体压力；x为沿裂缝距离；T_adj为非达西传导系数；μ为流体粘度；T为传导系数；Q为达西流量；B为无量纲常数；ρ为流体密度。

5.根据权利要求1所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，利用边界单元法建立所述力学计算的控制和本构方程，所述力学计算的控制和本构方程包括：

正应力的边界条件方程；

剪应力的边界条件方程；

在最大剪切应力方向上滑动的条件方程。

6.根据权利要求5所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，所述正应力的边界条件方程表示为：

其中，为裂缝远场载荷应力；P为流体压力；Δσ_n为裂缝单元累计变形形成的诱导应力引起的正应力的变化。

7.根据权利要求5所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，所述剪应力的边界条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；η为辐射阻尼系数；v为滑移速度。

8.根据权利要求5所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，所述在最大剪切应力方向上滑动的条件方程表示为：

其中，和为由于远场载荷的剪切应力分量；Δτ_d和Δτ_s为由于受力的累积变形而引起的剪应力的变化；v_d为倾斜方向的剪切应力；v_s为走向剪切应力。

9.根据权利要求1所述的页岩水力压裂三维全耦合离散裂缝网络模拟方法，其中，所述交叉作用准则包括：提前指定可能形成裂缝的位置，假设裂缝在开始模拟之前与最小主应力方向垂直。

10.一种页岩水力压裂三维全耦合离散裂缝网络模拟系统，其特征在于，该系统包括：

存储器，其上存储有计算机可执行指令；

处理器，所述处理器执行计算机可执行指令时实现以下步骤：

生成裂缝，提前指定可能形成的水力裂缝的位置，建立离散裂缝网络；

对所述离散裂缝网络进行网格划分，划分为矩形的裂缝单元；

建立流体流动的控制和本构方程和力学计算的控制和本构方程，模拟水力裂缝的扩展和天然裂缝的张开和剪切作用；

利用迭代法对所述流体流动的控制和本构方程和所述力学计算的控制和本构方程进行求解，并设定交叉作用准则来控制水力压裂裂缝的起裂、扩展和穿越。