CN113779843A - 一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法 - Google Patents

Abstract

本发明公开了一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法，包括：根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的固体变形，获得固体变量；根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量；通过迭代，获得稳定状态下的固体变量和流体变量；稳定状态中固体变量和流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展；根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展，基于计算结果更新流体驱动压裂裂缝网络，并获得流体驱动压裂裂缝网络的形态。本发明的并行有限元‑离散元‑有限体积法，能有效地提高计算效率，相比其他数值方法可更高效地模拟岩体裂缝扩展过程。

Description

一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法

技术领域

本发明涉及岩石力学与工程技术领域，特别是涉及一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法。

背景技术

开采深部致密油气储层中的非常规油气藏主要的开采手段是大规模的水力压裂、超临界二氧化碳压裂等流体驱动压裂形成压裂缝网的工程技术，数值模拟可以较为有效地模拟流体驱动裂缝的起裂和扩展的全部动态过程。由于流体驱动多孔弹性岩体裂缝扩展过程涉及到多物理场、多尺度断裂等复杂过程，导致目前针对流体驱动压裂过程的数值模拟方法往往具有一定局限性，无法完全高效、可靠地模拟裂缝的起裂和扩展：(1)多物理场耦合分析方面，相对准确地模拟流体驱动压裂中的流-固耦合效应可以更加准确的模拟岩体裂缝扩展行为，流体和固体的物理场耦合以及动态演化过程分析使得分析过程复杂、计算过程耗时；(2)工程尺度模型分析方面，在传统的串行计算方法在求解岩体裂缝扩展问题时有一定的局限性，在进行工程尺度的数值模拟时往往需要数周乃至数月的时间，这样的计算效率无法满足实际工程需要。因此，开发有效的并行计算方法和方案成为实际工程中的需求。

发明内容

由于高效、可靠的数值方法和模拟技术的需求，及针对深部致密岩体的多物理场耦合、多尺度断裂行为、工程尺度储层体积压裂实施等复杂因素，对计算方法的精确性、高效性、适用性等均提出较高要求。本发明研究发展已经建立的岩体体积压裂有限元-离散元-有限体积方法与断裂准则，对固体变形、流体压力、单元连接与分离等计算过程构建基于OpenMP的并行计算方案。本发明的目的是提供一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法，以解决上述现有技术存在的问题。

为实现上述目的，本发明提供了如下方案：

本发明提供一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法，包括：

根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的固体变形，获得固体变量；

根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量；

通过迭代，获得稳定状态下的所述固体变量和所述流体变量；

基于所述稳定状态，所述固体变量和所述流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展；

根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展，基于计算结果更新水力压裂裂缝网络，并获得水力压裂裂缝网络的形态。

可选的，根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的岩体变形，获得固体变量之前，所述方法还包括：根据所述多孔弹性岩体建立数值模型，进行网格划分，基于流-固耦合方程将多孔弹性介质转换为固体有限元模型和流体有限体积模型。

可选的，根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的岩体变形，获得固体变量的过程中包括：通过所述有限元法计算固体变形的过程中，使用了OpenMP并行化来计算，获得固体变量，所述固体变量包括位移、速度和加速度。

可选的，所述使用了OpenMP并行化来计算的过程中，还包括：构建动力学惯量影响的平衡方程：

式中，u(x,y,z)＝(u(x,y,z),v(x,y,z),w(x,y,z))^T代表位移向量；

与

分别代表了速度向量与加速度向量；σ^e代表考虑孔弹性岩体介质的有效应力张量；f代表外力向量；ρ代表密度；c代表阻尼系数；

代表方程的微分算子；Ω代表求解域。

可选的，根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量的过程中，包括：将多孔弹性岩石基质破裂面上的裂缝和天然裂缝中流体流动的流体压力作为外荷载，通过公式(1)实现；

其中所述流体变量包括流体压力和速度。

可选的，所述流体压力的计算过程包括：利用有效应力原理，计算多孔弹性岩体基质中流体流动的流体压力，通过如下公式(2)实现：

σ^e＝σ-αpI (2)

式中，α为孔弹性岩体介质Biot常数，p为空隙水压力，I为单位矩阵。

可选的，所述断裂准则包括拉伸破坏准则

和剪切破坏准则

式中，σ为拉伸应力；τ为剪切应力；ε₀为拉伸应力达到最大值时的拉伸应变值；ε_f为发生拉伸断裂时的拉伸应变值；γ₀为剪切应力达到最大值时的剪切应变值；γ_f为发生剪切断裂时的剪切应变值；G_tf和G_sf分别是拉伸断裂能和剪切断裂能。

可选的，所述固体变量和所述流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展的过程中包括：

满足所述拉伸破坏准则和所述剪切破坏准则任何一个准则，就会发生裂缝起裂或扩展，形成所述流体驱动拉伸和所述剪切裂缝扩展。

可选的，根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展的过程中包括：

对固体变形的数值进行离散，再进行求解：

为求解固体变形的控制方程，有限元法将平衡方程转化为单元e上的矩阵形式：

M^eD(t)+C^eD(t)+K^eD(t)＝F^e,x,y,z∈Ω (5)

式中，D(t)表示由单元e的节点位移组合而成的位移向量；M^e，C^e和K^e分别表示质量矩阵、阻尼矩阵、刚度矩阵；

和

分别表示t时刻的节点速度和加速度的向量；F^e表示外荷载向量。其中F^e为：

式中，

表示体积力；

表示断裂表面的流体压力；

表示相邻单元节点力；

表示牵引边界上的力。

可选的，求解过程中，通过结合动态松弛法和显式迭代法进行求解，包括：

计算外部节点力、内部节点力和总节点力，其中所述外部节点力通过所述公式(6)计算获得；

计算节点加速度、节点速度和位移；

重复以上计算过程，直至收敛

本发明公开了以下技术效果：

本发明的并行有限元-离散元-有限体积法，可以有效地考虑流-固耦合效应，从而相比与其他数值方法可以更高校地模拟岩体裂缝扩展过程。通过并行计算可以提高计算效率，这使得基于OpenMP有限元法-离散元法-有限体积并行计算方法能更快速地完成数值模拟过程，更好地将该数值方法应用于实际工程情况。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例的流程示意图。

具体实施方式

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

本发明针对流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法，包括：根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的固体变形，获得固体变量；根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量；通过迭代，获得稳定状态下的固体变量和流体变量；基于稳定状态，固体变量和流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展；根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展，基于计算结果更新流体驱动压裂裂缝网络，并获得流体驱动压裂裂缝网络的形态。本发明的并行有限元-离散元-有限体积法，可以采用并行化方案有效地提高流-固耦合下流体驱动多孔弹性岩体裂缝动态扩展的计算效率，从而相比与其他数值方法可以更高效地模拟岩体裂缝扩展过程。流体驱动多孔弹性岩体介质裂缝中流-固耦合分析和并行计算的总体方案流程如上图所示。

首先，根据实际岩体得到物理、力学参数和几何信息。

然后，分析了多孔弹性岩体介质，并将其转换为有限元和有限体积模型。

在程序流程图中，左手边和右手边分别显示了用有限元法计算固体变形和用有限体积法计算流体压力均进行多线程并行化计算。固体的基本变量，如位移、速度和应力，由有限元计算，流体的基本变量，如流体压力和速度，由有限体积法计算。

通过迭代，固体变形和流体流动处于稳定状态，并获得稳定状态下的位移、应力、流体流动和流体压力。

最后，基于当前的应力状态，一旦离散元模型满足裂缝产生的断裂准则，流体驱动拉伸和剪切断裂扩展，离散元处理单元分离，本步骤亦进行多线程并行化计算。

随后，在下一轮的新裂缝网络中计算注入流体流量的增加，以流体压力驱动岩体压裂裂缝的扩展，并获得当前流体下的压裂裂缝网络形态。

流体驱动多孔弹性岩体断裂的控制方程，其中多孔弹性介质固体变形：

(1)考虑动力学惯量影响的平衡方程：

式中，u(x,y,z)＝(u(x,y,z),v(x,y,z),w(x,y,z))^T代表位移向量；

与

分别代表了速度与加速度向量，并且

与

均为位移向量对时间的导数；σ^e代表考虑孔弹性岩体介质的有效应力张量；f代表包括体积力，流体对裂隙表面的压力，相邻单元节点力与作用在交界面上的力在内的外力向量；ρ代表密度；c代表阻尼系数；

代表方程的微分算子；Ω代表求解域。有效应力张量σ^e由下式给出：

σ^e＝σ-αpI (2)

式中，α为孔弹性岩体介质Biot常数，p为空隙水压力，I为单位矩阵。

(2)几何方程

式中，u表示位移场，ε为小应变假设(Cauchy应变)下的固体应变张量。

(3)本构方程

σ^e＝D:ε (4)

式中，D是弹性系数张量。

(4)固体介质的边界条件

式中，

表示了边界Γ_u的位移场，

表示了作用在固体边界Γ_σ的外力场。

(5)固体介质的初始条件

a.本质边界条件：

位移边界Γ_u (5)

b.自然边界条件：

位移边界Γ_u

σ·n＝-pn,裂缝边界

式中，位移场

描述了边界Γ_u的位移，p描述了作用在裂缝边界上

的流体压力，

描述了作用在固体边界Γ_σ的外力场。

多孔弹性介质流体流动：

(1)基于达西定律，假设岩体(即含流体的多孔弹性岩体介质)中的流体流动为单相流动且忽略重力的作用，可将流动速度场进行简化得到速度场方程：

式中，v_m与v_f分别表示多孔弹性岩体介质中流体流动的速度场与水力裂缝中流体流动的速度场；p_m与p_f分别表示多孔弹性岩体介质中的压力变量与水力裂缝中的压力变量；k_m与k_f分别代表多孔弹性岩体介质的渗透性与水力裂缝的渗透性；μ表示粘滞系数。

(2)假设不考虑重力和毛细管力的影响，可得到多孔弹性岩体介质与水力裂缝中单相流动的压力场方程：

其中，S_m＝n/K_f,S_f＝n/K_f分别是多空介质的饱和度和水力裂缝的饱和度；n表示孔隙率，其中裂缝的孔隙度等于1；q表示外部流入流速；ε_V表示岩石基质的体积应变。压裂液流体通过以上方程，控制其在岩体的孔隙和裂缝中流动，可以自然实现裂缝面处压裂液相岩体基质的流动，实现压裂液的漏失效应。

(3)边界条件：

外边界Γ_q (9)

外边界上Γ_p (10)

其中，流体来源

表示在外边界Γ_q上的流体源泉，n_q表示在这个边界上的单位外法向向量；流体压力

表示在外边界上Γ_p的外部流体压力。

流-固耦合效应在本发明中使用的有限元-离散元-有限体积法中包括三个方面：

1.将多孔弹性岩石基质破裂面上的流体驱动裂缝和天然裂缝中流体流动的流体压力作为外荷载，通过公式(1)实现。

2.利用有效应力原理，计算多孔弹性岩体基质中流体流动的流体压力，通过公式(2)实现。

3.多孔弹性岩体基质变形对流体流动的影响采用体积应变方程，通过公式(8)实现。

断裂准则：

固体的破坏、滑动和断裂发生在单元之间的界面上。用最大拉应力准则和莫尔-库仑强度准则来判断单元的拉伸和剪切破坏：

式中，σ为拉伸应力；τ为剪切应力；ε₀为拉伸应力达到最大值时的拉伸应变值；ε_f为发生拉伸断裂时的拉伸应变值；γ₀为剪切应力达到最大值时的剪切应变值；γ_f为发生剪切断裂时的剪切应变值；G_tf和G_sf分别是拉伸断裂能和剪切断裂能。

本发明还进行了并行计算的分析：(1)在计算精度方面，使用的基于OpenMP的有限元-离散元-有限体积法与KGD模型和PKN模型得到的结果基本相同，具有较高的准确性，本发明所使用的并行算法可以准确的模拟多条件下流体驱动岩体裂缝扩展行为；(2)在计算效率方面，通过不同线程并行效率进行分析，发现八线程的加速比与七线程的加速比基本相同，加速比达到3.56，这说明了本发明所使用的并行算法在使用八线程进行并行计算时加速比达到最大值，继续增大线程数并不会有效地进一步提高计算效率；通过对比不同线程并行计算的结果，不同线程并行计算下的结果完全相同，说明了本发明所使用的并行计算方法的可靠性，可在不同计算环境下采用不同的线程数进行并行计算，从而更好地适应多种计算环境；(3)在单元类型方面，利用六面体单元和五面体建模和求解分析，验证了本专利所使用的基于OpenMP的有限元-离散元-有限体积并行计算方法在不同单元形势下解答有效性和不同网格密度下解答收敛性；(4)在工程尺度模型计算方面，本发明利用高效的并行计算方法研究了不同地应力、不同流体类型条件下工程尺度岩体裂缝扩展行为。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (10)

1.一种流体驱动多孔弹性岩体裂缝动态扩展的并行计算方法，其特征在于，包括以下步骤：

根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的固体变形，获得固体变量；

根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量；

通过迭代，获得稳定状态下的所述固体变量和所述流体变量；

基于所述稳定状态，所述固体变量和所述流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展；

根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展，基于计算结果更新水力压裂裂缝网络，并获得水力压裂裂缝网络的形态。

2.根据权利要求1所述的方法，其特征在于，根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的岩体变形，获得固体变量之前，所述方法还包括：根据所述多孔弹性岩体建立数值模型，进行网格划分，基于流-固耦合方程将多孔弹性介质转换为固体有限元模型和流体有限体积模型。

3.根据权利要求1所述的方法，其特征在于，根据基于OpenMP的有限元法多线程并行计算多孔弹性岩体的岩体变形，获得固体变量的过程中包括：通过所述有限元法计算固体变形的过程中，使用了OpenMP并行化来计算，获得固体变量，所述固体变量包括位移、速度和加速度。

4.根据权利要求3所述的方法，其特征在于，所述使用了OpenMP并行化来计算的过程中，还包括：构建动力学惯量影响的平衡方程：

式中，u(x,y,z)＝(u(x,y,z),v(x,y,z),w(x,y,z))^T代表位移向量；

与

分别代表了速度向量与加速度向量；σ^e代表考虑孔弹性岩体介质的有效应力张量；f代表外力向量；ρ代表密度；c代表阻尼系数；

代表方程的微分算子；Ω代表求解域。

5.根据权利要求4所述的方法，其特征在于，根据基于OpenMP的有限体积法多线程并行计算流体流动，获得流体变量的过程中，包括：将多孔弹性岩石基质破裂面上的裂缝和天然裂缝中流体流动的流体压力作为外荷载，通过公式(1)实现；

其中所述流体变量包括流体压力和速度。

6.根据权利要求5所述的方法，其特征在于，所述流体压力的计算过程包括：利用有效应力原理，计算多孔弹性岩体基质中流体流动的流体压力，通过如下公式(2)实现：

σ^e＝σ-αpI (2)

式中，α为孔弹性岩体介质Biot常数，p为空隙水压力，I为单位矩阵。

7.根据权利要求1所述的方法，其特征在于，所述断裂准则包括拉伸破坏准则和剪切破坏准则；

式中，σ为拉伸应力；τ为剪切应力；ε₀为拉伸应力达到最大值时的拉伸应变值；ε_f为发生拉伸断裂时的拉伸应变值；γ₀为剪切应力达到最大值时的剪切应变值；γ_f为发生剪切断裂时的剪切应变值；G_tf和G_sf分别是拉伸断裂能和剪切断裂能。

8.根据权利要求2或7所述的方法，其特征在于，所述固体变量和所述流体变量大于断裂准则时，多孔弹性岩体发生流体驱动拉伸和剪切裂缝扩展的过程中包括：

满足所述拉伸破坏准则和所述剪切破坏准则任何一个准则，就会发生裂缝起裂或扩展，形成所述流体驱动拉伸和所述剪切裂缝扩展。

9.根据权利要求2或4所述的方法，其特征在于，根据离散元法计算多孔弹性岩体的流体驱动拉伸和剪切裂缝扩展的过程中包括：

对固体变形的数值进行离散，再进行求解：

为求解固体变形的控制方程，有限元法将平衡方程转化为单元e上的矩阵形式：

M^eD(t)+C^eD(t)+K^eD(t)＝F^e,x,y,z∈Ω (5)

式中，D(t)表示由单元e的节点位移组合而成的位移向量；M^e，C^e和K^e分别表示质量矩阵、阻尼矩阵、刚度矩阵；

和

分别表示t时刻的节点速度和加速度的向量；F^e表示外荷载向量，所述F^e为：

式中，

表示体积力；

表示断裂表面的流体压力；

表示相邻单元节点力；

表示牵引边界上的力。

10.根据权利要求9所述的方法，其特征在于，求解过程中，通过结合动态松弛法和显式迭代法进行求解，包括：

计算外部节点力、内部节点力和总节点力，其中所述外部节点力通过所述公式(6)计算获得；

计算节点加速度、节点速度和位移；

重复以上计算过程，直至收敛。

Images