CN109063239B - 一种水热耦合的三维数值模拟方法 - Google Patents

Abstract

本发明公开了一种水热耦合的三维数值模拟方法，包括裂隙岩体以非连续面作为约束进行网格划分，将裂隙岩体划分为多个实体单元，且实体单元的边界与既有的非连续面重合；设置起粘结作用的节理单元或者界面单元；将非连续面处的节理单元或者界面单元标记为断裂；对裂隙岩体的岩石基质进行热传导计算；进行裂隙内流体的热力学计算；进行流体与岩石基质进行热量交换的计算；重复上述步骤至完成所有裂缝的计算。本发明能够模拟含任意复杂裂缝网络岩体中的流体流动和热量运移；且可以对水热耦合问题进行模拟。

Description

一种水热耦合的三维数值模拟方法

技术领域

本发明涉及能源和环境及岩土工程领域，尤其涉及一种水热耦合的三维数值模拟方法。

背景技术

随着地热资源的利用，干热岩地热资源的开发，大坝工程、深层采油作业、寒区隧道建设以及核废料储存库、环境工程中污染源的迁移和深部开采等工程实践的深入和发展，人们已越来越认识到裂隙岩体中流体流动和热量运移相互作用的重要性。然而现有的有限元-离散元耦合分析方法(FDEM)、非连续变形分析方法(DDA)、数值流形方法(NMM)上不能考虑水热耦合问题，因而也就难以用于地热、石油天然气等领域进行模拟分析。

发明内容

有鉴于此，本发明的实施例提供了一种在裂隙岩体流体流动和热量运移相互作用过程中能考虑水热耦合问题的水热耦合的三维数值模拟方法。

本发明的实施例提供一种水热耦合的三维数值模拟方法，包括以下步骤：

S1.裂隙岩体以非连续面作为约束进行网格划分，将裂隙岩体划分为多个实体单元，且实体单元的边界与既有的非连续面重合；

S2.在实体单元之间设置起粘结作用的节理单元或者界面单元；

S3.将非连续面处的节理单元或者界面单元标记为断裂；

S4.对裂隙岩体的岩石基质进行热传导计算；

S5.在步骤S3标记为断裂的节理单元或界面单元上进行裂隙内流体的热力学计算；

S6.在标记为断裂的节理单元或界面单元上进行流体与岩石基质进行热量交换的计算；

S7.重复上述步骤至完成所有裂缝的计算。

进一步，所述步骤S1中，裂隙岩体以节理、裂隙作为非连续面，采用四面体、六面体、二维/三维voronoi、三角形、四边形、其他任意多边形或者多面体中的一种或多种进行网格划分。

进一步，所述步骤S4中，岩石基质进行热传导的计算公式为：

式中：

为岩石基质中在前一时步和当前时步的温度，Δt为时间步长，Q_s为岩石基质的热流速，Q_e为岩石基质和流体的热交换量，C_p为比热容，M为岩石基质质量。

进一步，所述步骤S5中，裂隙内流体的热力学计算公式如下：

式中：

为裂隙内流体在前一时步和当前时步的温度，Δt为时间步长，Q_f为裂隙内流体的热流速，C_p为比热容，ρ_f为流体的密度，V为裂隙流体的体积。

进一步，所述步骤S6中，流体与岩石基质进行热量交换的计算公式如下：

式中：裂隙两侧边界上的岩石基质温度分别为

和

裂隙内流体的温度为T_f，裂缝的面积为A，岩石基质和流体的换热系数为h。

与现有技术相比，本发明具有以下有益效果：能够模拟含任意复杂裂缝网络岩体中的流体流动和热量运移；为有限元-离散元耦合分析方法(FDEM)、非连续变形分析(DDA)、数值流形分析方法(NMM)增添新的功能，使得这些方法可以对水热耦合问题进行模拟。

附图说明

图1是本发明一种水热耦合的三维数值模拟方法的一流程图。

图2是本发明一实施例裂隙岩体的一示意图。

图3是本发明一实施例对裂隙岩体进行网格划分后插入节理单元或者界面单元的一示意图。

图4是本发明一实施例中对非连续面处的节理单元或者界面单元进行标记的一示意图。

图5是本发明一实施例中流体与岩石基质进行热量交换的一示意图。

图6是本发明一实施例中含裂隙网络岩体的一示意图。

图7是本发明一实施例中模拟得到岩石温度随时间演化的一示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。

请参考图1，本发明的实施例提供了一种水热耦合的三维数值模拟方法，包括以下步骤：

S1.裂隙岩体以非连续面作为约束进行网格划分，将裂隙岩体划分为多个实体单元，且实体单元的边界与既有的非连续面重合；裂隙岩体以节理、裂隙作为非连续面，采用四面体、六面体、二维/三维voronoi、三角形、四边形、其他任意多边形或者多面体中的一种或多种进行网格划分。

S2.在实体单元之间设置起粘结作用的节理单元或者界面单元；

S3.将非连续面处的节理单元或者界面单元标记为断裂；

S4.对裂隙岩体的岩石基质进行热传导计算；

计算公式为：

式中：

为岩石基质中在前一时步和当前时步的温度，Δt为时间步长，Q_s为岩石基质的热流速，Q_e为岩石基质和流体的热交换量，C_p为比热容，M为岩石基质质量。

S5.在步骤S3标记为断裂的节理单元或界面单元上进行裂隙内流体的热力学计算；

计算公式如下：

式中：

为裂隙内流体在前一时步和当前时步的温度，Δt为时间步长，Q_f为裂隙内流体的热流速，C_p为比热容，ρ_f为流体的密度，V为裂隙流体的体积。

S6.在标记为断裂的节理单元或界面单元上进行流体与岩石基质进行热量交换的计算；

计算公式如下：

式中：裂隙两侧边界上的岩石基质温度分别为T_s ⁺和T_s ^-，裂隙内流体的温度为T_f，裂缝的面积为A，岩石基质和流体的换热系数为h。

S7.重复上述步骤至完成所有裂缝的计算。

本发明能够模拟含任意复杂裂缝网络岩体中的流体流动和热量运移；为有限元-离散元耦合分析方法(FDEM)、非连续变形分析(DDA)、数值流形分析方法(NMM)增添新的功能，使得这些方法可以对水热耦合问题进行模拟。

实施例1

对裂隙岩体，以节理、裂隙等非连续作为约束进行四面体、六面体、二维/三维voronoi、四边形、其他任意多边形或者多面体等其中一种或多种实体单元划分网格，生成的实体单元边界与既有的节理、裂隙等非连续面重合；这些实体单元之间设置起粘结作用的节理单元或者界面单元；将非连续面处的节理单元或者界面单元标记为断裂。在标记为断裂的节理单元或界面单元上进行裂隙渗流的计算；在标记为断裂的节理单元或界面单元上进行流体传导和对流传热计算；在标记为断裂的节理单元或界面单元上进行流体与岩石基质换热的计算；根据以上步骤循环，即完成含任意复杂裂缝网络岩体流体和岩石中水热耦合的分析。

以采用四面体单元对裂隙岩体进行网格剖分为例说明本发明。

设有如图2所示的裂隙岩体，将裂缝网络作为约束进行网格划分。

在四面体单元之间插入节理单元或者界面单元如图3所示。

将非连续面处的节理单元或者界面单元标记为断裂如图4所示；

在标记为断裂的节理单元或界面单元上进行热力学计算，即裂隙渗流、流体热传导和流体热对流计算；

流体在裂隙中流动，流体通过对流及热传导传热，流体的对流及热传导均在裂缝面上进行。因为在四面体单元的公共面上插入了起粘结作用的节理单元，节理单元断裂即形成裂缝面。因此岩体中的裂缝可由一系列的三角面组成，流体仅在这些三角面上流动。这样流体流动、流体热传导、流体热对流，本质上属于二维问题。如图3所示，为一裂缝面，由于已将天然裂缝作为约束条件划分四面体单元，并在四面体单元的公共边上插入起粘结作用的节理单元，因此划分网格后，该裂缝面自然变成由一系列空间三角面构成，如图4所示。相邻三角面之间的公共点构成计算的交汇结点，例如图4中的结点1-8。由于裂缝面各处的温度是不同的，因此可用交汇结点的温度来表示整个裂缝上的温度分布。任意一个裂缝三角面都构成了一个二维面状热流动路径。如图4中结点1,2,3所构成的裂缝三角面Δ123就构成了一个二维面状热流动路径。下面介绍，怎样基于图4所示的网格来计算整个裂缝面的流体温度分布。

以图4中的结点1为例，与节点1连接的三角面有Δ123，Δ134，Δ145，Δ156，Δ167，Δ78。以三角面Δ123为例，介绍一个三角面内的热流量求解步骤。设节点1,2,3的温度为T₁,T₂,T₃，在该三角面构成的局部坐标中(x,y位于三角面内，z轴垂直于三角面)，设这三个结点的局部坐标为(x_k,y_k,0)(k为节点号，k＝1,2,3)。假设三角面内的温度分布服从线性分布，那么三角面内任意一点的温度梯度为常量，可表示为：

由高斯散度定理，(1)式可写为：

其中，A是三角形单元的面积，n_i是外法向向量，

是m边的平均温度，

是m边的两顶点的坐标分量之差，∈_ij是二维置换张量。

根据立方定律，在i方向上的热流速可表示为：

其中q_i为i方向的热流速，

为流体的热传导系数张量，T为温度。

这样，将式(2)代入(10)式，即可得沿x,y方向的热流速，设为q_x,q_y.

于是，单位时间内，流入节点1的热流量可通过下式计算：

其中

是节点1所对的边的外法向单位向量，L⁽ⁿ⁾为节点1所对的边的边长。

这样我们便求得了三角形单元Δ123中流入节点1的热流量Q_Δ123。与上类似，我们可求得其他与节点1直接相连的三角形单元流入节点1的热流量。这样单位时间内，因为流体热传导流入节点1的总热流量可表示为：

裂隙内流体热对流模型

如图4所示，两个节点之间的热流量除了热传导产生的以外，还有通过热对流引起的流量，仍然以三角形单元Δ123为例，根据方程(2)，可求得三角形单元内任意一点的温度梯度为

这样，三角形单元Δ123内因为热对流而流入节点1热流量为

其中q_f是三角形单元Δ123内的流体流速矢量，采用如下方式求解，假设三角面内的流体压力分布服从线性分布，那么三角面内任意一点的流体压力梯度为常量，可表示为：

由高斯散度定理，(8)式可写为：

其中，A是三角形单元的面积，n_i是外法向向量，

是m边的平均流体压力，

是m边的两顶点的坐标分量之差，∈_ij是二维置换张量。

根据立方定律，在i方向上的流体流速可表示为：

其中q_i为i方向的流体流速，a为裂缝的水力张开度，μ为流体动力粘性，p为流体压力。

这样，将式(9)代入(10)式，即可得沿x,y方向的流速，便求得了最终流速矢量q_f.

由于节点1还与其他三角形单元直接相连，同理可求得其他三角形单元内通过流体对流流入节点1的热流量Q_Δ134，Q_Δ145，Q_Δ156，Q₁₆₇，Q₁₇₈，Q₁₈₂因此节点1因热对流而产生的总热流量为：

然后可据下式求得节点1的当前温度为：

其中，

为节点1在前一时步和当前时步的温度，C_f是流体的比热容，ρ_f是流体的密度，Q_f＝Q_f1+Q_f2为节点1的总热流量(流出为负，流入为正)，Δt为时间步长，其中V节点1的所有节理单元体积的一半。

在标记为断裂的节理单元或界面单元上进行流体与岩石基质进行热量交换的计算。

由于岩石基质的温度和裂隙中流体的温度通常会不一样，因此在裂隙处岩石基质和裂隙内流体会发生热量交换。这样如何处理流体在裂缝和岩石基质之间水岩热交换成为求解全域内温度分布的关键。为此，我们采用如下方式予以处理。设位于裂隙两侧边界上的岩石基质温度分别为T_s ⁺和T_s ^-，如图5所示，裂隙内流体的温度为T_f，裂缝的面积为A，岩石基质和流体的换热系数为h。于是单位时间内，岩石基质和流体之间的热量交换为：

这样进行裂隙内流体的热力学计算最终按下式进行更新：

而后续的岩石基质的热传导计算，最终采用下式进行更新：

根据以上步骤循环，即完成含任意复杂裂缝网络岩体流体和岩石中水热耦合的分析。

下面结合具体的实例对本方案的数值模拟方法进行说明：

如图6所示，为一长方体岩块，长为2m，宽为1m，含有J1-J9组成的裂缝网络。岩体的初始温度为100度，从左边界注入恒温度为20度的流体，流体流经裂缝网络并在裂缝处与岩石发生热量交换，然后从右侧边界流出。

如图7是采用本文模拟方法得到岩石温度随时间的演化，从这些图形可以看出裂缝处的岩石温度逐渐降低，可以对整个水岩换热过程进行模拟。

在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种水热耦合的三维数值模拟方法，其特征在于，包括以下步骤：

S1.裂隙岩体以非连续面作为约束进行网格划分，将裂隙岩体划分为多个实体单元，且实体单元的边界与既有的非连续面重合；

S2.在实体单元之间设置起粘结作用的节理单元或者界面单元；

S3.将非连续面处的节理单元或者界面单元标记为断裂；

S4.对裂隙岩体的岩石基质进行热传导计算；

S5.在步骤S3标记为断裂的节理单元或界面单元上进行裂隙内流体的热力学计算；

S6.在标记为断裂的节理单元或界面单元上进行流体与岩石基质进行热量交换的计算；

S7.重复上述步骤至完成所有裂缝的计算。

2.根据权利要求1所述的水热耦合的三维数值模拟方法，其特征在于，所述步骤S1中，裂隙岩体以节理、裂隙作为非连续面，采用四面体、六面体、二维/三维voronoi、三角形、四边形、其他任意多边形或者多面体中的一种或多种进行网格划分。

3.根据权利要求1所述的水热耦合的三维数值模拟方法，其特征在于，所述步骤S4中，岩石基质进行热传导的计算公式为：

式中：

为岩石基质中在前一时步和当前时步的温度，Δt为时间步长，Q_s为岩石基质的热流速，Q_e为岩石基质和流体的热交换量，C_p为比热容，M为岩石基质质量。

4.根据权利要求3所述的水热耦合的三维数值模拟方法，其特征在于，所述步骤S5中，裂隙内流体的热力学计算公式如下：

式中：

为裂隙内流体在前一时步和当前时步的温度，Δt为时间步长，Q_f为裂隙内流体的热流速，C_p为比热容，ρ_f为流体的密度，V为裂隙流体的体积。

5.根据权利要求4所述的水热耦合的三维数值模拟方法，其特征在于，所述步骤S6中，流体与岩石基质进行热量交换的计算公式如下：

式中：裂隙两侧边界上的岩石基质温度分别为

和

裂隙内流体的温度为T_f，裂缝的面积为A，岩石基质和流体的换热系数为h。

Images