一种基于PFC离散元的页岩水力压裂数值模拟方法

一种基于pfc离散元的页岩水力压裂数值模拟方法
技术领域
1.本发明涉及页岩气藏水力压裂技术领域，具体涉及一种基于pfc离散元的页岩水力压裂数值模拟方法。


背景技术：

2.我国页岩气储量巨大，探明可采页岩气资源量高达31.6万亿m3，资源潜力巨大。页岩气藏是一个“自生、自储”系统，页岩储层物性极差，多为纳米孔喉，孔隙度常为1％-7％，基质渗透率较小。国内外多采用水平井分段多簇压裂技术对页岩储层进行整体压裂，从而改善储层渗流能力，增大储层改造体积，增加页岩气产量。水力压裂作为有效开采页岩气的核心技术之一，广泛应用于石油天然气等行业，通常向页岩储层注入足够压力的压裂液，页岩断裂产生大量的水力裂缝，形成气液高速渗流通道，以改造页岩气储层，达到增产的目的。从水力压裂技术被开发至今，最受关注的是水力压裂作用下岩石的破裂机理和水力裂缝产生的方向和大小。尽管室内实验能够探索裂缝起裂规律，但难以精确表征岩石力学行为和水力压裂过程中的流动-位移-应力耦合问题。同时，数值模拟具有成本低、可重复及试验周期较短的优势，因而得到了更广泛的应用，大量工程实际的检验也验证了数值模拟的准确性，应用前景十分广阔。
3.由于页岩储层的复杂性，在解决水力压裂问题时，常规的有限元法、块体离散元法等数值模拟方法表现出一定的局限性。近些年来，基于离散元的数值模拟方法在处理岩土工程的问题上得到了广泛的应用。离散元方法与其他数值模拟方法对比，其最大的优势就是离散元数值模拟不受物体变形量的影响，针对多相介质耦合作用下的力学行为和岩体的破裂等问题更具有仿真性，能够直观的看出也岩体的变化过程和结果。pfc离散元法主要是通过圆盘与圆盘之间的相互作用来模拟页岩中颗粒单元之间的力学行为。pfc程序是在基于离散元理论的基础上通过显式差分的运算法则综合研发得到的，pfc离散元法与普通连续介质的数值模拟方法对比，最大的优点就是能够充分体现出岩体的实体结构及其内部颗粒单元之间的相互作用，这种以颗粒单元为基本研究对象，在针对页岩颗粒单元力学行为的观测与研究时更为直观，是一种运用非连续的计算方法处理复杂变形问题的基本理论。正是由于pfc离散元是以颗粒单元为基本的计算单元，以牛顿第二定律为基本的力学关系，因此在pfc离散元数值模拟中对颗粒单元性质的描述更能表现出研究对象的本质性。


技术实现要素：

4.为解决现有技术存在的不足，本发明提出了一种基于pfc离散元的页岩水力压裂数值模拟方法，基于pfc离散元方法探索水力压裂裂缝扩展机理，将渗流与离散元进行耦合，模拟页岩储层水力压裂裂缝起裂机制，进而指导页岩储层压裂设计；明确页岩水力压裂起裂机理和裂缝最终展布形态，揭示水力裂缝与页岩交互作用问题，实现页岩储层水力压裂最终裂缝展布形态的精准预测。
5.采用的技术方案为：
6.一种基于pfc离散元的页岩水力压裂数值模拟方法，包括以下步骤：
7.步骤一：建立外部墙体单元，在墙体内生成颗粒单元，在中心位置构建圆形区域，删除颗粒单元，建立注液井眼；
8.根据单轴压缩实验结果、巴西劈裂实验结果与单轴压缩模拟结果、巴西劈裂模拟结果对比分析输出的应力应变曲线，标定页岩细观参数；
9.赋予颗粒单元基本物性参数，通过平行粘结单元粘结相邻的颗粒单元，生成页岩储层离散元数值模型；
10.所述页岩细观参数包括所有颗粒单元基本物性参数及平行粘结单元参数；
11.步骤二：设定相邻颗粒单元围成的区域为流体域，设定颗粒单元接触垂直处为流体管道，在页岩储层离散元数值模型上建立网络流动模型；
12.步骤三：在页岩边界墙体施加围压，设定井眼注入流体的速度，进行渗流-位移-应力耦合的数值模拟计算；
13.步骤四：实时监测颗粒单元间粘结状态，判断粘结单元是否发生破坏，当粘结单元所承受应力大于粘结单元的法向强度或切向强度时，粘结单元为破坏状态，视为产生微裂缝；计算裂缝数量及水力裂缝扩展形态；计算结束后，导出并保存页岩储层水力裂缝扩展形态模拟结果。
14.所述流体域的变量计算如下：
15.在二维模型中，通过立方定律得到平行板通道的流量，如式(1)所示：
[0016][0017]
其中，w为通道的开度，单位m；δp为相邻流体域单元之间的压力差，单位pa；l为通道的长度，单位m；μ为流体的粘度，pa
·
s；
[0018]
流体域单元中每一个时步的压力变化，如式(2)所示：
[0019][0020]
其中，kf为流体体积模数，单位pa；vd为流体域单元的表观体积，单位m3；δt为时间步长，单位s。
[0021]
所述颗粒单元，粘结后的法向作用力、切向作用力和弯矩通过平行粘结模型计算；
[0022]
当法向接触力为拉力时，如式(3)所示：
[0023]
w＝w0 ng
ꢀꢀꢀ
(3)
[0024]
其中，w0为颗粒间接触力为0时对应的通道开度，称为残余开度，n为距离压缩因子，g为相邻颗粒表面的法向距离；
[0025]
当法向接触力为压力时,如式(4)所示：
[0026][0027]
其中，f0为通道开度是时颗粒单元间接触力的大小，f为颗粒单元之间的接触力，根据颗粒单元之间的叠合量，按力-位移定理得出。
[0028]
所述力-位移定理对于“颗粒单元-颗粒单元”接触，单位法向量ni的定义，如式(5)所示：
[0029][0030]
其中，和分别为颗粒单元1和颗粒单元2的位置坐标，d为两球球心之间的距离；在“颗粒单元-墙体”接触中，单位法向量ni则定义为球心距离墙体最大直线的方向；
[0031]“颗粒单元-颗粒单元”或“颗粒单元-墙体”之间的“重叠”量，如式(6)所示：
[0032][0033]
颗粒单元-颗粒单元或颗粒单元-墙体接触产生的接触力，分解为法向接触力和切向接触力如式(7)所示：
[0034][0035]
其中，法向接触力计算，如式(8)所示：
[0036][0037]
其中，kn为接触处的法向刚度，是一个割线模量，使总的法向位移与法向接触力发生联系；
[0038]
接触力矢量fi产生的切向接触力以增量的形式计算，当颗粒单元间的接触形成时，总的切向接触力先初始化为零，之后颗粒单元间的每个相对位移增量引起的弹性切向增量会累加到当前值上，从而产生切向接触力，通过更新每个时间步长δt内的单位法向量ni接触中心点xc的坐标来对接触运动，新的时间步长内产生的切向接触力等于当前时间步长开始时已经存在的切向接触力与切向弹性力增量之和，切向接触力增量如式(9)所示：
[0039][0040]
其中，ks为接触切向刚度，为切向速度；将法向接触力和切向接触力的值代入式中计算最终接触力。
[0041]
所述微裂缝的长度等于两颗粒单元的平均粒径，方向垂直于初始接触方向。
[0042]
本发明的有益技术效果：
[0043]
1、本发明提出的基于pfc离散元的页岩水力压裂模拟方法能够从微观尺度精确表征水力裂缝起裂和扩展机理，准确捕捉水力裂缝宏观展布形态。
[0044]
2、本发明提出的基于pfc离散元的页岩水力压裂模拟方法，能够精确再现水力裂缝与页岩交互作用，使得到的模拟结果与实际情况更加吻合。
附图说明
[0045]
图1为本发明实施例提供的pfc离散元进行页岩水力压裂数值模拟流程框图；
[0046]
图2为本发明实施例提供的页岩数值模型图；
[0047]
图3为本发明实施例提供的pfc软件2d模拟单轴压缩试验初始示意图；
[0048]
图4为本发明实施例提供的pfc软件2d模拟单轴压缩试件应力应变图；
[0049]
图5为本发明实施例提供的pfc软件2d模拟巴西劈裂试验初始示意图；
[0050]
图6为本发明实施例提供的pfc软件2d模拟巴西劈裂试件应力应变图；
[0051]
图7为本发明实施例提供的水力裂缝最终展布形态；
[0052]
图8为本发明实施例提供的水力裂缝最终展布形态速度云图。
具体实施方式
[0053]
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。
[0054]
本实施例提出了一种基于pfc离散元的页岩水力压裂数值模拟方法，基于pfc离散元方法探索水力压裂裂缝扩展机理，将渗流与离散元进行耦合，模拟页岩储层水力压裂裂缝起裂机制，进而指导页岩储层压裂设计；明确页岩水力压裂起裂机理和裂缝最终展布形态，揭示水力裂缝与页岩交互作用问题，实现页岩储层水力压裂最终裂缝展布形态的精准预测。
[0055]
采用的技术方案为：
[0056]
一种基于pfc离散元的页岩水力压裂数值模拟方法，如图1所示，包括以下步骤：
[0057]
步骤一：建立外部墙体单元，在墙体内生成颗粒单元，在中心位置构建圆形区域，删除颗粒单元，建立注液井眼；
[0058]
根据单轴压缩实验结果、巴西劈裂实验结果与单轴压缩模拟结果、巴西劈裂模拟结果对比分析输出的应力应变曲线，标定页岩细观参数；
[0059]
赋予颗粒单元基本物性参数，通过平行粘结单元粘结相邻的颗粒单元，生成页岩储层离散元数值模型，如图2所示；
[0060]
所述页岩细观参数包括所有颗粒单元基本物性参数及平行粘结单元参数；
[0061]
本实施例中，页岩水力压裂离散元数值模型的页岩细观参数由试验与模拟对比标定；建立标准试件100mm
×
50mm的二维页岩单轴压缩数值模型，如图3所示，进行单轴压缩数值模拟试验，输出单轴压缩应力-应变曲线，如图4所示，将页岩单轴压缩数值模拟结果与页岩单轴压缩试验结果对比；建立标准试件r＝50mm的二维页岩巴西劈裂数值模拟模型，如图5所示，进行巴西劈裂数值模拟试验，输出巴西劈裂应力-应变曲线，如图6所示，将页岩巴西劈裂数值模拟结果与页岩巴西劈裂试验结果对比；通过结果对比，标定页岩细观参数；
[0062]
步骤二：设定相邻颗粒单元围成的区域为流体域，设定颗粒单元接触垂直处为流体管道，在页岩储层离散元数值模型上建立网络流动模型；
[0063]
所述流体域的变量计算如下：
[0064]
在二维模型中，通过立方定律得到平行板通道的流量，如式(1)所示：
[0065][0066]
其中，w为通道的开度，单位m；δp为相邻流体域单元之间的压力差，单位pa；l为通道的长度，单位m；μ为流体的粘度，单位pa
·
s；
[0067]
流体域单元中每一个时步的压力变化，如式(2)所示：
[0068][0069]
其中，kf为流体体积模数，单位pa；vd为流体域单元的表观体积，单位m3；δt为时间步长，单位s；
[0070]
步骤三：在页岩边界墙体施加围压，设定井眼注入流体的速度，进行渗流-位移-应力耦合的数值模拟计算；
[0071]
步骤四：实时监测颗粒单元间粘结状态，判断粘结单元是否发生破坏，当粘结单元所承受应力大于粘结单元的法向强度或切向强度时，粘结单元为破坏状态，视为产生微裂缝；计算裂缝数量及水力裂缝扩展形态；计算结束后，导出并保存页岩储层水力裂缝扩展形态模拟结果；
[0072]
所述颗粒单元，粘结后的法向作用力、切向作用力和弯矩通过平行粘结模型计算；
[0073]
当法向接触力为拉力时，如式(3)所示：
[0074]
w＝w0 ng
ꢀꢀꢀ
(3)
[0075]
其中，w0为颗粒间接触力为0时对应的通道开度，称为残余开度，n为距离压缩因子，g为相邻颗粒表面的法向距离；
[0076]
当法向接触力为压力时,如式(4)所示：
[0077][0078]
其中，f0为通道开度是时颗粒单元间接触力的大小，f为颗粒单元之间的接触力，根据颗粒单元之间的叠合量，按力-位移定理得出；
[0079]
所述力-位移定理对于“颗粒单元-颗粒单元”接触，单位法向量ni的定义，如式(5)所示：
[0080][0081]
其中，和分别为颗粒单元1和颗粒单元2的位置坐标，d为两球球心之间的距离；在“颗粒单元-墙体”接触中，单位法向量ni则定义为球心距离墙体最大直线的方向；
[0082]“颗粒单元-颗粒单元”或“颗粒单元-墙体”之间的“重叠”量，如式(6)所示：
[0083][0084]
颗粒单元-颗粒单元或颗粒单元-墙体接触产生的接触力可分解为法向接触力和切向接触力如式(7)所示：
[0085][0086]
其中，法向接触力计算，如式(8)所示：
[0087][0088]
其中，kn为接触处的法向刚度，是一个割线模量，使总的法向位移与法向接触力发生联系；
[0089]
接触力矢量fi产生的切向接触力以增量的形式计算，当颗粒单元间的接触形成时，总的切向接触力先初始化为零，之后颗粒单元间的每个相对位移增量引起的弹性切向增量会累加到当前值上，从而产生切向接触力，通过更新每个时间步长δt内的单位法向量ni接触中心点xc的坐标来对接触运动，新的时间步长内产生的切向接触力等于当前时间步长开始时已经存在的切向接触力与切向弹性力增量之和，切向接触力增量如式(9)所
示：
[0090][0091]
其中，ks为接触切向刚度，为切向速度；将法向接触力和切向接触力的值代入式中计算最终接触力；
[0092]
所述微裂缝的长度等于两颗粒单元的平均粒径，方向垂直于初始接触方向；
[0093]
本实施例中，根据上述模拟即可得到页岩水力压裂裂缝最终展布形态，如图7所示；在该状态下，水力裂缝起裂和扩展受到应力状态影响，水力裂缝沿最大主应力方向扩展，最终形成双翼缝形态；并且，水力裂缝更有可能穿透页岩，因为它们倾向于向一个单一的方向扩展，且能量充足；
[0094]
图8为水力裂缝扩展速度场云图；颗粒单元运动速度云图呈现出沿水力裂缝扩展方向对称分布的规律；在高压水的作用下，水力压裂裂缝在尖端产生，向最大水平主应力的方向延伸扩展。