一种多孔弹性岩体热-流-固耦合压裂裂缝扩展模拟方法与流程

一种多孔弹性岩体热
‑
流
‑
固耦合压裂裂缝扩展模拟方法
技术领域
1.本发明涉及矿业新能源(煤层气、致密气、页岩气)的勘探与开发技术领域，特别是涉及一种多孔弹性岩体热
‑
流
‑
固耦合压裂裂缝扩展模拟方法。


背景技术：

2.致密油气储层中岩体的温度、应力及压裂液的流动是影响岩体水力压裂裂缝扩展行为的三个主要因素，三者之间相互耦合、相互作用、相互制约，形成裂隙岩体中温度场与渗流场、应力场(或应变场)之间复杂的耦合效应，因此进行岩体压裂的热
‑
流
‑
固(thermal
‑
hydro
‑
mechanical，简称thm)耦合机理分析非常必要。简单地说，在温度场、渗流场、应力场中，热作用可以在岩体中产生热应力并引起岩体的弹性模量、泊松比等固体特性的变化，同时由于温度变化引起压裂液密度的变化，影响了压裂液的流动。压裂液的流动可以通过传导或对流使热量更快地扩散开来，而岩体中压裂液的存在可以改变岩体的受力情况。岩体由于受到变形影响会在内部产生一定的热量耗散，而力学变形也会部分影响固体热学特性的变化，同时力学效应对于空隙度或流体的渗透性能也具有比较明显的影响。
3.目前关于岩体的热
‑
流
‑
固耦合研究主要有两种方法—离散法与等效连续法。离散法则是在确定每条裂隙的空间方位、隙宽等几何参数的情况下，以单条裂隙水流基本公式为基础，利用流入、流出各裂隙交叉点流量相等条件来求其水头值，而热学分析与此大致类似。离散法建立的这种模型假定岩块为刚性或可变形的，不连续面由节理单元模拟，更接近于研究实际。通过离散法的处理技术，发展出采用离散颗粒组装表征固体结构的离散单元法，该方法通过刚性颗粒体与颗粒体之间的粘结力(如弹簧力)形成结构体，该方法中离散颗粒体的数目、大小、颗粒粘结力等对解答影响显著，数值分析的工作量较大。等效连续法就是把裂隙中的渗流量、热流等平均到岩体中去，在得到等效连续介质模型之后利用经典的多孔介质分析方法来进行研究；但只要样本的代表性单元体积存在，并且相对于研究网格来说不是很大，则应该尽量采用等效连续介质模型进行分析。对于岩体热
‑
流
‑
固耦合，常用的数值模拟方法有三种：有限差分法、有限元法以及扩展有限元法。利用有限差分法或有限元法，将控制方程在时间、空间上离散，结合初始条件和边界条件即可计算出每个时间步中热
‑
流
‑
固耦合过程的变化情况。采用有限元模型，存在一定的局限性，网格在创建后不容易改变，需要进一步的网格细化。从数学角度看，有限差分法比有限元法近似程度高，但应用上后者远比前者简单灵活；有限体积法在有限差分法的基础上进行发展，目前多被应用于处理流体流动问题。经过开发的扩展有限元方法可以模拟多个水力裂缝，通过单元改善了裂缝扩展，但是计算结果的准确性仍然取决于网格生成质量。
4.将连续体中有限元应力解的精度与离散元中裂缝扩展的灵活性相结合，产生了一种耦合方法，为处理水力压裂问题提供了一种潜在的技术：将考虑多孔弹性岩体在水力压裂过程中的热(温度)、流体、固体等物理场耦合效应，利用有限元法和有限体积法计算各物理量，并通过离散单元法处理固体断裂，实现多孔弹性岩体在热
‑
流
‑
固耦合下水力压裂拉伸与剪切型裂缝的扩展模拟。


技术实现要素：

5.本发明的目的在于采用数值模拟方法考虑多孔弹性岩体水力压裂过程中的热
‑
流
‑
固耦合效应，从而得出更加精确的裂缝扩展行为和三个物理场模拟结果。
6.为实现上述目的，本发明提供了如下方案：
7.一种多孔弹性岩体热
‑
流
‑
固耦合压裂裂缝扩展模拟方法，包括：
8.s1、将多孔弹性岩体介质转化为有限元法模型和有限体积法模型，考虑流体驱动岩体水力压裂的固体和流体物理场耦合效应；
9.s2、基于所述有限元模型计算固体的变形，基于所述有限体积模型计算流体流动；
10.s3、基于所述多孔弹性岩体介质，建立压裂液与岩石基质间的热传递模型，计算所述压裂液与岩石基质间的热传递；
11.s4、将所述固体的变形、所述流体流动和所述压裂液与岩石基质间的热传递三个变量进行更新迭代，用于使固体变形和流体流动以及热扩散处于平衡状态；
12.s5、基于岩石材料的微观结构模型，在所述s4中的平衡状态下，若所述微观结构模型满足裂缝产生的断裂准则，采用离散元法处理固体断裂，产生水力拉伸和剪切裂缝的扩展，获得新的水力压裂裂缝网络；
13.s6、对所述水力拉伸和剪切裂缝进行评价，并分析所述水力压裂裂缝网络的形态。
14.优选的，所述s2中，通过有限元法计算所述固体的变形，通过有限体积法计算所述流体流动。
15.优选的，所述固体的变形包括固体的位移、速度和加速度，所述流体流动包括流体的孔隙压力和流速。
16.优选的，s3中，在所述压裂液与岩石基质间的热传递模型中，所述压裂液区域与所述岩石基质通过网格单元和多孔弹性岩体地层节点之间的接触单元进行传热。
17.优选的，所述网格单元向裂缝中流体传递的热量为：
18.q
c1
＝α
c
(t
n
)(t
n
‑
t
f1
)
19.q
c2
＝α
c
(t
n
)(t
n
‑
t
f2
)
20.其中，q
c1
和q
c2
分别为网络节点和多孔弹性岩体地层节点之间的接触热流，α
c
为网格节点t
n
内随温度变化的接触导热系数，t
f1
和t
f2
分别为对应的地层节点温度。
21.优选的，s5中，所述岩石材料的微观结构模型包括有限元单元和有限体积单元，模型配置为在所述多孔弹性岩体介质中固体使用有限元模型，多孔弹性岩体介质中流体采用有限体积模型。
22.优选的，所述岩石材料的微观结构模型中的每个块体需满足固体变形控制方程、流体流动控制方程以及断裂准则。
23.优选的，s5中，由所述水力拉伸引起的分离界面的法相位移w为：
[0024][0025]
其中，和分别是相邻两个有限元单元a和b的分离节点1和2的位移增量，n
12
是节点1和2之间接触的法向单位向量。
[0026]
本发明的有益效果为：
[0027]
本发明能够提供一种多孔弹性岩体在热
‑
流
‑
固耦合作用下压裂裂缝扩展的数值模拟方法，将连续体中有限元应力解的精度与离散元中裂缝扩展的灵活性相结合，并通过用于流体流动和热效应计算的有限体积法来完成未解决的热
‑
流
‑
固耦合效应分析，产生了一种高性能的有限元
‑
离散元
‑
有限体积耦合方法，实现多孔弹性岩体在热
‑
流
‑
固耦合下水力压裂拉伸与剪切型裂缝的扩展模拟。
附图说明
[0028]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0029]
图1为本发明方法流程图；
[0030]
图2为本发明实施例中通过网格和多孔弹性岩体地层节点之间的接触单元进行传热过程示意图；
[0031]
图3为本发明实施例中划分的几何区域和接触节点，其中a、b、c、d为单元，(a)为初始几何域，(b)为有限元域，(c)为离散有限元/离散元域，(d)为离散元域；
[0032]
图4为本发明实施例中a、b两个有限元单元接触节点结构示意图，其中为弹簧法向刚度，为弹簧切向刚度；
[0033]
图5为本发明实施例中流体在两个有限体积单元之间通过单元中心和界面流动示意图，其中ω
i
和ω
j
分别是有限体积单元i和j的域，其接口是q
ij
是从i到j的通量；
[0034]
图6为本发明热(温度)、流体、固体物理场的基本变量迭代过程示意图。
具体实施方式
[0035]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0036]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0037]
本发明提供一种基于热
‑
流
‑
固(thermal
‑
hydro
‑
mechanical，简称thm)耦合效应的压裂数值模拟方法，将连续体中有限元应力解的精度与离散元中裂缝扩展的灵活性相结合，产生了一种耦合方法，并通过用于流体流动和热效应计算的有限体积法来完成未解决的热
‑
流
‑
固耦合效应。为了计算多孔弹性岩体介质中固体变形和流体流动的控制方程，在算法中固体变形和断裂分析中采用了有限元和离散元方法，在流体流动和热效应的计算中采用了有限体积法法。利用上述数值方法，利用有限元
‑
有限体积
‑
离散元耦合方法，实现多孔弹性岩体在热
‑
流
‑
固耦合下水力压裂拉伸与剪切型裂缝的扩展模拟。
[0038]
首先，建立多孔弹性岩体储层的数值模型，引入几何和物理参数。然后，对以多孔弹性岩体介质进行分析，将其转化为有限元和有限体积模型。在程序流程图(如附图1所示)
的左侧和右侧分别显示了用有限元法计算的固体变形和用有限体积计算的流体流动。固体的基本变量(如位移、速度和应力)采用有限元法计算；流体的基本变量(如流体压力和速度)以及热扩散的基本变量(如温度梯度和热通量)采用有限体积法计算，将固体的变形、流体流动和压裂液与岩石基质间的热传递三个变量进行更新迭代，三个物理场场的基本变量迭代过程见附图6。岩石基质的固体变形与断裂采用有限元
‑
离散元算法(fe
‑
de算法)进行计算，岩体孔隙
‑
裂隙结构内的流体流动采用有限体积算法(fv算法)计算出流体压力场，岩体与流体的热传递采用有限体积算法(fv算法)计算出温度场。在迭代过程中，固体变形场传递到岩体孔隙
‑
裂隙结构中，影响结构的变形，从而影响流体压力场和温度场分布；流体压力场传递到岩体孔隙
‑
裂隙结构中，影响固体结构的受力和变形；温度场的变化，同时影响固体和流体的物理参数性质。上述迭代过程实现各物理场间的传递和相互迭代，形成热、流、固的耦合过程。通过迭代，固体变形和流体流动以及热扩散处于平衡状态，并获得稳定状态下的位移、应力、流体流动、流体压力以及温度。随后，基于当前的应力状态，一旦离散元模型满足裂缝产生的断裂准则，水力拉伸和剪切裂缝就会扩展。最后，在下一轮的新裂缝网络中计算注入水力压裂的压裂液流体，增加压裂液流量，驱动多孔弹性岩体裂缝扩展，并给出此过程中热、流、固等等物理场的变量解答。以下展开介绍上述过程的基本方程和求解方法。
[0039]
岩石材料的微观结构模型包括有限元单元和有限体积单元，模型配置为在所述多孔弹性岩体介质中固体使用有限元模型，多孔弹性岩体介质中流体采用有限体积模型。此外，每个单元的强度与其变形模量有关，根据材料的应力状态和本构关系得到每个单元的变形,单元运动的控制方程是考虑单元变形的平衡方程。数值模型必须满足固体变形、流体流动、热扩散的基本控制方程以及断裂准则，利用有限元法、有限体积法实现了多孔弹性岩体介质中流体流动、热扩散的全耦合离散化。
[0040]
(1)固体变形控制方程
[0041]
考虑动惯量影响的多孔弹性岩体的固体变形平衡方程为：
[0042][0043]
其中，u是位移矢量；和分别表示速度和加速度矢量，表示位移相对于时间t的导数；σ
e
代表考虑孔弹性岩体介质的有效应力张量；f指外部载荷矢量，包括体力、断裂面上的流体压力、弹簧力和牵引边界上的力；ρ是密度；c是阻尼系数；ω是求解域。为了完整性起见，式(1)中岩体动态变形过程考虑了惯性项与阻尼项。
[0044]
为了求解固体变形的控制方程，可以将平衡方程(1)转换为单元е上的以下矩阵形式：
[0045][0046]
其中，d(t)是通过组合单元e的节点位移组成的所谓位移矢量；m
e
、c
e
和k
e
分别是质量，阻尼和刚度矩阵；和分别表示时间t处的节点速度和加速度的矢量；和k
e
d(t)分别表示时间t处节点的合外力、阻力和弹簧弹力；f
e
是外部加载向量，f
e
可以表示为：
[0047][0048]
其中，是体积力；是断裂面上的流体压力；是弹簧力；是牵引边界上的力。上述方程可以采用常规的显示迭代方法进行求解。
[0049]
(2)压裂液和岩石基质间热传递计算
[0050]
压裂液区域与岩石基质通过接触和多孔弹性岩体地层节点的传热如图2所示。网格单元从裂缝中流体传递的热量：
[0051]
q
c1
＝α
c
(t
n
)(t
n
‑
t
f1
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0052]
q
c2
＝α
c
(t
n
)(t
n
‑
t
f2
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0053]
其中,q
c1
和q
c2
分别为网络节点和多孔弹性岩体地层节点之间的接触热流，α
c
为网格节点t
n
内随温度变化的接触导热系数，t
f1
和t
f2
分别为对应的地层节点温度。
[0054]
(3)流体速度场和压力计算
[0055]
水力压裂裂缝网络和原始多孔介质构成了压裂多孔介质，如图5所示。假设裂缝性多孔介质中为单相流，可采用简化形式的、忽略重力的多孔介质和水力裂缝中流体流动的达西定律来获得速度场：
[0056][0057][0058]
式中，ν
m
和ν
f
分别为孔隙介质和水力裂缝流体流动的速度场；pm和p
f
分别代表了多孔介质和水力裂缝的压力值；k
m
和k
f
分别表示孔隙介质和水力裂缝的渗透率；μ是粘度。
[0059]
在没有重力和毛细力的情况下，多孔介质和水力裂缝中单相流动的压力方程为：
[0060][0061][0062]
其中,分别为多孔介质和水力裂缝的蓄水系数；n为孔隙度，裂缝孔隙度等于1；q代表外部流体源；ε
ν
表示岩石基质的体积应变。
[0063]
以单元为中心的有限体积法适用于多孔介质和水力裂缝中的流体流动，为了在固体位移和流体流动的孔隙压力之间进行数据传递，在单元的节点上计算固体的变量(例如，位移，速度和加速度)，而流体的变量(例如，孔隙压力和流速)在单元的中心进行计算。
[0064]
(4)固体的断裂计算
[0065]
固体的破坏、滑动和断裂发生在单元之间的界面上，本方法采用最大拉应力准则和莫尔
‑
库仑强度准则来判断单元的拉、剪断裂：
[0066]
拉伸断裂的最大拉应力准则：
[0067]
剪切断裂的莫尔
‑
库仑准则：
[0068]
其中，σ和τ分别是法向应力和切向应力；抗拉强度；c是粘聚力；是内摩擦角。
[0069]
将单元的连续界面转换为断裂的不连续表面的策略描述如下。如果三院法向或切向应力满足方程(10)和(11)中的断裂准则，则单元的界面将分离，新分离的单元界面在水力压力的作用下继续变形，最终演变为典型的拉伸型或剪切型水力裂缝。水力裂缝的宽度由如下三院节点引起的分离界面法向位移决定：
[0070][0071]
其中，和分别是相邻两个有限元单元a和b的分离节点1和2的位移增量，n
12
是节点1和2之间接触的法向单位向量。
[0072]
图3给出了算法的断裂过程，其中初始几何域被离散为有限元域和离散元域，有限元/离散元单元之间的接触点是弹簧。如果满足断裂条件(即上文介绍的断裂准则)，则固体会在节点处分离而破裂，裂缝在单元边界处形成。如图4所示，为单元的刚度矩阵，其各分量表示弹簧的法向和切向刚度，为弹簧法向刚度，为弹簧切向刚度。
[0073]
本发明提供一种基于热
‑
流
‑
固(thermal
‑
hydro
‑
mechanical，简称thm)耦合效应的压裂数值模拟方法，将连续体中有限元应力解的精度与离散元中裂缝扩展的灵活性相结合，产生了一种耦合方法，并通过用于流体流动和热效应计算的有限体积法来完成未解决的热
‑
流
‑
固耦合效应。为了计算多孔弹性岩体介质中固体变形和流体流动的控制方程，在算法中固体变形和断裂分析中采用了有限元和离散元方法，在流体流动和热效应的计算中采用了有限体积法。利用上述数值方法，利用有限元
‑
有限体积
‑
离散元耦合方法，实现多孔弹性岩体在热
‑
流
‑
固耦合下水力压裂拉伸与剪切型裂缝的扩展模拟。
[0074]
以上所述的实施例仅是对本发明优选方式进行的描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。