CN110263362A - 基于孔隙密度流的岩土体离散元流固耦合数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，包括生成固体颗粒随机堆积模型、剖分识别孔隙网络、孔隙渗流方程建立、孔隙流体及相邻固体颗粒对固体的作用、固体位移对孔隙流体的作用、更新孔隙渗流参数，重复以上步骤直至固体颗粒平衡及孔隙流体渗流稳定，本发明极大地减小了计算量，并建立孔隙流体状态方程，通过密度天然地将温度场与渗流场耦合起来，类比宏观达西定律建立细观渗流方程实现孔隙渗流计算，基于孔隙尺度能高效地模拟较为复杂的宏观现象。

Description

基于孔隙密度流的岩土体离散元流固耦合数值模拟方法

技术领域

本发明涉及一种岩土体流固耦合的模拟方法，尤其是一种基于孔隙密度流的岩土体离散元法的数值模拟方法。

背景技术

岩土体是由固液气组成的三相复杂体系。地质工程中的灌浆、石油工程中的油气运移及开采、隧道的突水突泥、液化问题等都与流体与固体的相互作用有关。液气等流体的存在使得岩土体性质具有很大的复杂性与变异性，研究流体在岩土体中的运移与相互作用具有重要意义。长期以来基于宏观连续介质的渗流力学与固体力学理论对于复杂的宏观现象，难以揭示其微观机制，离散元法自提出以来，由于其对非连续介质较好的模拟能力，能清晰反映颗粒材料微观变形机制等优点，迅速被引入对岩土体材料的研究中。

基于离散元法的孔隙尺度的流固耦合方法主要包括有LBM-DEM法(格子玻尔兹曼-离散元法)和离散元平板流法。前一种方法将流体划分为不同密度的离散流体单元，格子步长通常比颗粒尺寸小一个量级，不同密度流体单元的压力不同，流体在压力差作用下发生运移，由于流体格子划分较小，这种方法可以模拟非常复杂的流体现象，但计算量巨大；后一种方法认为流体密度不变，将流体划分为孔隙流体和孔喉通道，流体运移导致体积改变，通过体积模量计算流体压力，计算量较小，但方法较为简化，无法模拟湍流等复杂流体现象。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，来模拟岩土体流固耦合问题。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，包括：(1)将岩土体离散为相互接触的固体颗粒，根据颗粒堆积骨架剖分识别孔隙，建立以孔隙流体单元为节点，孔喉通道相连接的孔隙流体拓扑差分网络；(2)根据流体密度与压力、温度的关系，建立孔隙流体单元的压力与密度、温度的函数关系式，以密度差异—压力差异驱动孔隙流体运动；(3) 相邻孔隙流体单元因密度、压力差异对流，固体颗粒通过接触力相互作用；孔隙流体对相邻固体颗粒产生压力作用，固体颗粒在相邻颗粒及流体合力作用下产生位移，一方面影响相邻孔隙体积及孔喉通道，另一方面影响孔隙流体拓扑网络。从而实现固体颗粒与流体间的相互作用；(4)基于MatDEM固体离散元-孔隙密度流的高效算法，通过时间步迭代计算，直至岩土体模型平衡，输出模拟结果；(5)考虑热效应。温度变化引起流体密度、压力变化，考虑孔隙流体间对流时亦传递热量，本发明亦能模拟温度场的耦合作用；(6)本发明从单个孔隙尺度入手，能模拟较复杂的宏观现象，同时不考虑单个孔隙内部具体流场，减小了计算量，具体包括以下步骤：

步骤1，将岩土体离散为相互接触的固体颗粒，根据颗粒尺寸、颗粒半径生成固体颗粒随机堆积模型。

步骤2，剖分识别孔隙网络：根据固体颗粒间距离设置阈值，将固体颗粒随机堆积模型孔隙空间剖分为一系列相互连通的较小孔隙，建立以孔隙流体单元为节点，孔喉通道相连接的孔隙流体拓扑差分网络，得到孔隙流体网络-固体颗粒骨架模型(图1)。

步骤3，孔隙渗流计算：

步骤31，根据孔隙流体温度和孔隙流体密度建立单个孔隙流体状态方程：

P＝f(ρ，T) (1)

其中，P为孔隙流体压力，ρ为孔隙流体密度，T为孔隙流体温度。

步骤32，相邻孔隙流体单元在压力差作用下沿相应孔喉通道渗流，类比达西定律，单个孔喉细观渗流方程：

对单个孔隙单元，流量流入为正，q为单个孔喉体积渗流量，K表示孔喉渗透系数，孔喉渗透系数K与孔喉宽度有关，A为孔喉通道面积，为平均压力梯度，其中P_j为相邻孔隙流体压力，P_i为中心孔隙流体压力，l为平均渗径。

步骤33，根据孔隙渗流量更新孔隙流体质量和温度：

更新孔隙流体单元质量：

其中，为第n+1次迭代时的第k个孔隙流体质量，为第n次迭代时的第k个孔隙流体质量，k表示第k个孔隙单元，Q^k为总质量渗流量。

更新孔隙流体单元温度：

其中，c为流体比热容，为第n+1次迭代时的第k个孔隙流体单元温度，为第n次迭代时的第k个孔隙流体单元温度，k表示第k个孔隙单元；T_i为对应渗流流体温度。

步骤4，孔隙流体及相邻固体颗粒对固体的作用：以固体颗粒为分析对象，计算相邻孔隙流体单元对固体颗粒的压力及相邻固体颗粒的作用力。根据牛顿力学及运动学方程计算固体颗粒合力、加速度、位移等。

步骤5，固体位移对孔隙流体的作用：根据固体颗粒位移更新孔隙体积及孔喉宽度，当固体颗粒相对位移较大时(图5b)，使得两颗粒距离小于距离阈值d_max，形成新的孔喉通道，孔隙间拓扑结构改变，重复步骤(2)重新剖分识别孔隙网络。

步骤6，更新孔隙渗流参数：根据步骤(2)、(5)更新孔隙流体质量、温度及体积，并计算孔隙流体密度，应用式(1)更新孔隙流体压力。同时根据步骤(5)的孔喉宽度更新孔喉渗透系数。

步骤7，重复步骤(3)-(6)直至固体颗粒平衡及孔隙流体渗流稳定。

优选的：步骤2得到孔隙流体网络-固体颗粒骨架模型中颗粒i和颗粒j之间颗粒距离：

d_ij＝D_ij-(R_i+R_j)

其中，d_ij表示颗粒i和颗粒j之间颗粒距离，D_ij为颗粒球心距离，R_i和R_j分别为颗粒i 和颗粒j半径，设置颗粒距离阈值d_max。

优选的：颗粒距离阈值d_max取0.1R，R表示两颗粒平均半径。

优选的：当颗粒距离d_ij<d_max时，识别为孔喉通道，并记录相应颗粒连接，通过闭环搜索算法将模型孔隙空间Ω剖分为n个较小的孔隙单元，记为孔隙单元Ω^k由N^k个固体颗粒围成，对应有N^k个孔喉通道与相邻孔隙连通。

优选的：根据中心孔隙压力、相邻孔隙压力建立平均压力梯度方程：

其中，J为平均压力梯度，ΔP＝P_j-P_i,P_j为相邻孔隙压力，P_i为中心孔隙压力，l为两孔隙平均渗径。

优选的：孔喉渗透系数K：

其中，K表示孔喉渗透系数，μ为流体粘滞系数，ω为等效孔喉宽度。

优选的：总质量渗流量Q^k：

其中，q_i表示第i个孔喉体积渗流量，ρ_i为流体密度。

优选的：等效孔喉宽度ω＝d+ω₀，d表示孔喉宽度，ω₀表示孔喉宽度d＝0时的等效孔喉宽度，与介质宏观渗透系数有关。

本发明相比现有技术，具有以下有益效果：

本发明提供的孔隙密度流流固耦合方法的优点在于：在单个孔隙尺度认为流体物理性质均一，极大地减小了计算量，并建立孔隙流体状态方程，通过密度天然地将温度场与渗流场耦合起来，类比宏观达西定律建立细观渗流方程实现孔隙渗流计算。本发明基于孔隙尺度能高效地模拟较为复杂的宏观现象。

附图说明

图1基于孔隙密度流的岩土体离散元流固耦合数值模拟实现流程图

图2孔隙渗流计算流程图

图3固体位移计算流程图

图4孔喉及孔隙渗流示意图，其中，图4a为孔喉示意图，图4b为孔隙渗流示意图。

图5颗粒位移示意图，其中，图5a为颗粒位移示意图1，图5b为颗粒位移示意图2。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

一种基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，如图1所示，包括以下步骤：

步骤10生成固体颗粒随机堆积模型。根据模型尺寸、颗粒半径等相关参数生成固体颗粒随机堆积模型；

步骤11剖分识别孔隙网络。定义相邻颗粒i、j距离(图4a)：

d_ij＝D_ij-(R_i+R_j)

其中D_ij为颗粒球心距离，R_i和R_j分别为颗粒i、j半径。设置颗粒距离阈值d_max(一般取0.1R)，当两颗粒距离d_ij<d_max时，识别为孔喉通道，并记录相应颗粒连接。通过闭环搜索算法将模型孔隙空间Ω剖分为n个较小的孔隙单元，记为孔隙单元Ω^k由N^k个固体颗粒围成，对应有N^k个孔喉通道与相邻孔隙连通；

步骤12孔隙渗流(图4b)。

图2为孔隙渗流计算流程图

步骤20计算单个孔隙流体单元压力。单个孔隙流体单元状态方程：

P＝f(ρ，T)

步骤21计算孔喉平均压力梯度。平均压力梯度：

其中ΔP＝P_j-P_i,P_j为相邻孔隙压力，P_i为中心孔隙压力，l为两孔隙平均渗径。

步骤22计算孔喉渗透系数。以二维渗流问题为例，对比泊肃叶定律，可取：

其中μ为流体粘滞系数，ω为等效孔喉宽度，一般取ω＝d+ω₀。

步骤23计算单个孔喉渗流量。定义流体流入时流量为正，单个孔喉由压力梯度引起体积渗流量：

q＝K(d)AJ

渗流量q正比于平均压力梯度J；A为孔喉通道截面积,对二维渗流问题，取A＝d；K为细观孔喉渗流系数，与孔喉宽度d有关。

步骤24计算单个孔隙总渗流量。对每个孔隙流体单元Ω^k，总质量渗流量：

q_i表示第i个孔喉体积渗流量，ρ_i为流体密度，取两渗流孔隙流体密度平均值。

步骤25更新孔隙流体单元质量。

其中为第n+1次迭代时的孔隙流体质量，为第n次迭代时的孔隙流体质量；

步骤26更新孔隙流体单元温度。由能量守恒：

得：

其中c为流体比热容，为第n+1次迭代时的孔隙流体质量，为第n次迭代时的孔隙流体质量；T_i为对应渗流流体温度；

步骤13孔隙流体及相邻固体颗粒对固体的作用。

图3为固体位移计算流程图。

步骤30计算单个孔隙流体压力；

步骤31计算单个颗粒所受相邻流体压力之和。

步骤32计算相邻固体颗粒作用力。

步骤33计算单个固体颗粒所受合力。

步骤34根据牛顿第二定律计算固体颗粒加速度。

步骤35根据运动学方程计算固体颗粒速度。

步骤36根据运动学方程计算固体颗粒位移。

步骤14固体位移对孔隙流体的作用。根据固体颗粒位移更新孔隙体积及孔喉宽度(图5a)；当相邻固体颗粒相对位移较大时(图5b)，孔隙拓扑结构改变，此时还需重新计算颗粒距离，重复步骤11重新剖分孔隙流体网络；

步骤15更新孔隙渗流参数。根据步骤25、步骤14更新孔隙流体质量和体积，并计算孔隙流体密度，结合步骤26更新流体温度并应用流体状态方程(式(1))更新孔隙流体压力；同时根据步骤14的孔喉宽度更新孔喉渗透系数；

步骤16判断固体颗粒是否平衡且流场稳定，否则重复步骤13-16直至迭代终止。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (8)

1.一种基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于，包括以下步骤：

步骤1，将岩土体离散为相互接触的固体颗粒，根据颗粒半径、颗粒级配生成固体颗粒随机堆积模型；

步骤2，剖分识别孔隙网络：根据固体颗粒间距离设置阈值，将固体颗粒随机堆积模型孔隙空间剖分为一系列相互连通的较小孔隙，建立以孔隙流体单元为节点，孔喉通道相连接的孔隙流体拓扑差分网络，得到孔隙流体网络-固体颗粒骨架模型；

步骤3，孔隙渗流方程建立：

步骤31，根据孔隙流体温度和孔隙流体密度建立单个孔隙流体状态方程：

P＝f(ρ，T) (1)

其中，P为孔隙流体压力，ρ为孔隙流体密度，T为孔隙流体温度；

步骤32，相邻孔隙流体单元在压力差作用下沿相应孔喉通道渗流，类比达西定律，单个孔喉细观渗流方程：

对单个孔隙单元，流量流入为正，q为单个孔喉体积渗流量，K表示孔喉渗透系数，A为孔喉通道面积，为平均压力梯度，其中P_j为相邻孔隙流体压力，P_i为中心孔隙流体压力，l为平均渗径；

步骤33，根据孔隙渗流量更新孔隙流体质量和温度：

更新孔隙流体单元质量：

其中，为第n+1次迭代时的第k个孔隙流体质量，为第n次迭代时的第k个孔隙流体质量，k表示第k个孔隙单元，Q^k为总质量渗流量；

更新孔隙流体单元温度：

其中，c为流体比热容，为第n+1次迭代时的第k个孔隙流体单元温度，为第n次迭代时的第k个孔隙流体单元温度，k表示第k个孔隙单元；T_i为对应渗流流体温度；

步骤4，孔隙流体及相邻固体颗粒对固体的作用：以固体颗粒为分析对象，计算相邻孔隙流体单元对固体颗粒的压力及相邻固体颗粒的作用力；根据牛顿力学及运动学方程计算固体颗粒合力、加速度、位移；

步骤5，固体位移对孔隙流体的作用：根据固体颗粒位移更新孔隙体积及孔喉宽度，当固体颗粒相对位移使得两颗粒距离小于距离阈值d_max，形成新的孔喉通道，孔隙间拓扑结构改变，重复步骤(2)重新剖分识别孔隙网络；

步骤6，更新孔隙渗流参数：根据步骤(2)、(5)更新孔隙流体质量、温度及体积，并计算孔隙流体密度，应用式(1)更新孔隙流体压力；同时根据步骤(5)的孔喉宽度更新孔喉渗透系数；

步骤7，重复步骤(3)-(6)直至固体颗粒平衡及孔隙流体渗流稳定。

2.根据权利要求1所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：

步骤2得到孔隙流体网络-固体颗粒骨架模型中颗粒i和颗粒j之间颗粒距离：

d_ij＝D_ij-(R_i+R_j)

其中，d_ij表示颗粒i和颗粒j之间颗粒距离，D_ij为颗粒球心距离，R_i和R_j分别为颗粒i和颗粒j半径，设置颗粒距离阈值d_max。

3.根据权利要求2所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：颗粒距离阈值d_max取0.1R，R表示两颗粒平均半径。

4.根据权利要求3所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：当颗粒距离d_ij＜d_max时，识别为孔喉通道，并记录相应颗粒连接，通过闭环搜索算法将模型孔隙空间Ω剖分为n个较小的孔隙单元，记为孔隙单元Ω^k由N^k个固体颗粒围成，对应有N^k个孔喉通道与相邻孔隙连通。

5.根据权利要求4所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：根据中心孔隙压力、相邻孔隙压力建立平均压力梯度方程：

其中，J为平均压力梯度，ΔP＝P_j-P_i，P_j为相邻孔隙压力，P_i为中心孔隙压力，l为两孔隙平均渗径。

6.根据权利要求5所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：孔喉渗透系数K：

其中，K表示孔喉渗透系数，μ为流体粘滞系数，ω为等效孔喉宽度。

7.根据权利要求6所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：总质量渗流量Q^k：

其中，q_i表示第i个孔喉体积渗流量，ρ_i为流体密度。

8.根据权利要求7所述基于孔隙密度流的岩土体离散元流固耦合数值模拟方法，其特征在于：等效孔喉宽度ω＝d+ω₀，d表示孔喉宽度，ω₀表示孔喉宽度d＝0时的等效孔喉宽度。