基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法

摘要

本发明公开了一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，该方法包括以下步骤：S1、开展砂雨法与拟振动压实法的三维颗粒离散元模拟，制备不同相对密实度的虚拟压缩试样；S2、对虚拟压缩试样逐级施加轴向压力，进行侧限压缩试验的三维颗粒离散元数值模拟，得到各级压力下压缩稳定后的离散元数值模型；S3、对各级压力压缩后的离散元数值模型分别进行切片，并对切片图像进行处理，建立对应的孔隙结构模型；S4、针对建立的孔隙结构模型开展渗流的格子Boltzmann数值模拟。本发明克服了现有室内压缩渗透试验技术的不足，且原理简单、实施方便、计算效率高、应用效果好，验证了基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法的可行性和优越性。

权利要求

1.一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，该方法包括以下步骤：

S1、采用砂雨法与拟振动压实法进行离散元模拟，制备不同相对密实度的虚拟压缩试样；包括：

S1.1、在高于试样的模型空间中，按级配要求生成彼此不重叠的颗粒集合体，施加重力荷载，进行离散元法迭代计算，迭代中每隔一定计算循环将颗粒速度重新归零以控制颗粒下落后的最大速度不超过一定的速度阈值，下落完成并计算至平衡后获得颗粒堆积体的最松散状态；

S1.2、在颗粒堆积体的顶部施加设定的轴向压力，并设定颗粒间摩擦系数为0近似模拟最优振动作用，计算平衡后获得颗粒堆积体的最密实状态；

S1.3、维持颗粒堆积体顶部的轴向压力不变，从最松散状态开始逐步减小颗粒间的摩擦系数来模拟振动作用，获得颗粒间摩擦系数与颗粒堆积体孔隙率的关系，通过该关系并结合最松散状态和最密实状态，估算出获得给定相对密实度堆积体所需的颗粒间摩擦系数；

S1.4、在所需的颗粒间摩擦系数下对最松散堆积体进行设定的轴向压力的压实，即可获得给定相对密实度下的颗粒堆积体；

S2、对虚拟压缩试样逐级施加轴向压力，进行侧限压缩试验的三维颗粒离散元数值模拟，得到各级压力下压缩稳定后的离散元数值模型；

S3、对各级压力压缩后的离散元数值模型进行切片，并对切片图像进行处理，建立对应的孔隙结构模型；包括：

S3.1、将每级压力下变形稳定后的圆柱形试样沿轴向用一系列等间隔的切面进行切割，导出各切面处的模型切片图像，包括侧面圆柱形墙在切片上的圆形边界及球形颗粒在切片上对应的圆；

S3.2、对导出的切片图像进行批处理，将切片图像上每个像素的尺寸调整为与相邻切片间的竖向间隔相等，识别各切片上属于孔隙的像素，对应于格子Boltzmann模型中的流体单元，以此实现对离散元数值试样的体素化，并建立试样对应的孔隙结构模型；

S4、针对建立的孔隙结构模型开展渗流的格子Boltzmann数值模拟；包括：S4.1、将每级压力下变形稳定后的孔隙结构模型导入格子Boltzmann计算程序中，设定计算参数，包括：粒子碰撞模型、离散速度模型、流体边界条件，进行格子Boltzmann迁移与碰撞的迭代计算；

S4.2、待渗流场收敛到稳态后，停止迭代计算，根据所施加的水压差利用达西定律计算该级压力下数值试样的渗透率。

2.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S1中，制备的不同相对密实度的虚拟压缩试样的高径比不宜超过0.6，以保证后续压缩时试样受力均匀。

3.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S1.1中，控制颗粒下落后的最大速度不超过一定的速度阈值为：0.8m/s。

4.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S1.2、步骤S1.3、步骤S1.4中，设定的轴向压力为：14kPa。

5.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S1.2、步骤S1.3、步骤S1.4、步骤S2中，轴向压力的实现方法为：在试样顶部生成加载板，加载板由三层颗粒组成，每层颗粒按六边形或四边形规则排列，颗粒间的接触本构模型采用线性平行粘结模型，平行粘结的刚度取为试样颗粒间接触

100

刚度的50‑100倍，平行粘结的强度设置为1.0×10 Pa，以此模拟不发生变形的刚性加载板，将所需要的荷载转换为加载板顶层颗粒上的集中力进行施加，计算公式为：其中，Fplaten1为加载板顶层每个颗粒上需施加的轴向集中力，σ为所设定的轴向压应力，A为虚拟压缩试样的横截面面积，Nplaten1为加载板顶层颗粒的数量。

6.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S1.3中，颗粒间摩擦系数与颗粒堆积体孔隙率的关系由若干数值试验测点通过指数函数 拟合得到；式中，n为孔隙率，μfric为颗粒间摩擦系数，a、b、c为拟合参数。

7.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S2中，每级轴向压力按10kPa/5000steps的加载速度缓慢施加，以避免加载时动力作用的影响，且轴向加载过程中试样的侧面和底部为刚性边界所约束。

8.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S3.2中，切片图像处理的方法包括：切片图像处理过程中，将每张切片图像上的圆形墙边界向内缩进一与模型中最小颗粒半径相等的距离，使得边界与附近颗粒相交叠以免后续渗流模拟时在边界附近产生集中渗流。

9.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S3.2中，切片图像处理的方法包括：切片图像处理过程中，对圆形边界以外、图片正方形边界以内的区域标记为固体区域；

除了识别切片图像上的孔隙单元，还将固体区域进一步区分为流‑固边界单元和固体内部单元两类，具体作法是：若某固体单元周围的邻近单元均为固体单元，邻近单元包括相邻切片上对应的相邻单元，则该固体单元标记为固体内部单元；若有一个或多个邻近单元为孔隙单元，即流体单元，则该单元标记为流‑固边界单元。

10.根据权利要求1所述的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，其特征在于，所述步骤S4.1中，具体的计算参数包括：流体粒子碰撞模型，采用单松弛时间BGK碰撞模型；离散速度模型，采用D3Q19模型；边界条件为，模型四周壁面和固体颗粒表面为无滑移的流‑固边界，采用反弹法进行处理，流体流动方向沿试样轴向，由入口和出口的压力差驱动，压力边界采用Zou/He法进行处理；

所述步骤S4.1中，渗流计算过程中固体内部单元不执行迁移与碰撞，以减少计算耗时。

说明书

基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟

方法

技术领域

[0001] 本发明涉及岩土体压缩渗透试验的数值模拟仿真技术领域，尤其涉及一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法。

背景技术

[0002] 公路、铁路、市政、机场、水利等工程建设中均存在填土、开挖、建筑物施工等外部荷载变化引起地基、边坡、隧洞等地质体渗流量变化的问题，这一问题的研究对路基填筑和
开挖后的沉降计算、边坡隧洞等的稳定性分析及水环境评价等具有重要的理论与工程意
义。

[0003] 目前，有部分学者通过室内压缩渗透试验测定破碎岩体、砂土、黏土等在不同压力作用下的压缩变形及变形后的渗透系数。然而，这种室内压缩渗透试验存在诸多不足，主要
包括：(1)由于要设置测压管，试样的高径比偏大，轴向压力作用时下部土体的竖向压应力
将明显偏低，造成整个试样压缩并非均匀；(2)试验过程中的轴向压力是通过砝码施加的，
砝码的瞬时加载引起土体孔隙的瞬时减小会产生一个冲击动水压，影响了试验结果；(3)渗
透试验过程中，由于试样的侧壁为刚性边界，在边界附件易产生集中渗流，影响了试验结
果；(4)试验比较耗时且难以洞悉试样内部细观结构与渗流场的变化，阻碍了对岩土体压缩
变形及变形过程中渗透性变化的细观机理的进一步研究。

[0004] 随着计算机和数值计算方法的发展，虚拟试验仿真技术在岩土工程领域得到了广泛的应用。在诸多数值计算方法中，颗粒离散元法尤其适用于岩土材料在外荷载作用下的
变形模拟，格子Boltzmann方法擅长于模拟复杂多孔介质中的渗流。因此，有必要将这两种
先进的数值模拟方法相结合，提出一种基于颗粒离散元与格子Boltzmann的压缩渗透试验
数值模拟方法，从而克服上述室内压缩渗透试验的不足，为后续科学研究和工程应用提供
技术支撑。

发明内容

[0005] 本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，通过将固体力学领域的颗粒离散元方法
与流体力学领域的格子Boltzmann方法相结合，实现对虚拟岩土体试样在侧向压缩下的变
形及变形后的渗透性的有效模拟。

[0006] 本发明解决其技术问题所采用的技术方案是：

[0007] 本发明提供一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，该方法包括以下步骤：

[0008] S1、采用砂雨法与拟振动压实法进行离散元模拟，制备不同相对密实度的虚拟压缩试样；包括：

[0009] S1.1、在高于试样的模型空间中，按级配要求生成彼此不重叠的颗粒集合体，施加重力荷载，进行离散元法迭代计算，迭代中每隔一定计算循环将颗粒速度重新归零以控制
颗粒下落后的最大速度不超过一定的速度阈值，下落完成并计算至平衡后获得颗粒堆积体
的最松散状态；

[0010] S1.2、在颗粒堆积体的顶部施加设定的轴向压力，并设定颗粒间摩擦系数为0近似模拟最优振动作用，计算平衡后获得颗粒堆积体的最密实状态；

[0011] S1.3、维持颗粒堆积体顶部的轴向压力不变，从最松散状态开始逐步减小颗粒间的摩擦系数来模拟振动作用，获得颗粒间摩擦系数与颗粒堆积体孔隙率的关系，通过该关
系并结合最松散状态和最密实状态，估算出获得给定相对密实度堆积体所需的颗粒间摩擦
系数；

[0012] S1.4、在所需的颗粒间摩擦系数下对最松散堆积体进行设定的轴向压力的压实，即可获得给定相对密实度下的颗粒堆积体；

[0013] S2、对虚拟压缩试样逐级施加轴向压力，进行侧限压缩试验的三维颗粒离散元数值模拟，得到各级压力下压缩稳定后的离散元数值模型；

[0014] S3、对各级压力压缩后的离散元数值模型进行切片，并对切片图像进行处理，建立对应的孔隙结构模型；包括：

[0015] S3.1、将每级压力下变形稳定后的圆柱形试样沿轴向用一系列等间隔的切面进行切割，导出各切面处的模型切片图像，包括侧面圆柱形墙在切片上的圆形边界及球形颗粒
在切片上对应的圆；

[0016] S3.2、对导出的切片图像进行批处理，将切片图像上每个像素的尺寸调整为与相邻切片间的竖向间隔相等，识别各切片上属于孔隙的像素，对应于格子Boltzmann模型中的
流体单元，以此实现对离散元数值试样的体素化，并建立试样对应的孔隙结构模型；

[0017] S4、针对建立的孔隙结构模型开展渗流的格子Boltzmann数值模拟；包括：

[0018] S4.1、将每级压力下变形稳定后的孔隙结构模型导入格子Boltzmann计算程序中，设定计算参数，包括：粒子碰撞模型、离散速度模型、流体边界条件，进行格子Boltzmann迁
移与碰撞的迭代计算；

[0019] S4.2、待渗流场收敛到稳态后，停止迭代计算，根据所施加的水压差利用达西定律计算该级压力下数值试样的渗透率。

[0020] 进一步地，本发明的所述步骤S1中，制备的不同相对密实度的虚拟压缩试样的高径比不宜超过0.6，以保证后续压缩时试样受力均匀。

[0021] 进一步地，本发明的所述步骤S1.1中，控制颗粒下落后的最大速度不超过一定的速度阈值为：0.8m/s。

[0022] 进一步地，本发明的所述步骤S1.2、步骤S1.3、步骤S1.4中，设定的轴向压力为：14kPa。

[0023] 进一步地，本发明的所述步骤S1.2、步骤S1.3、步骤S1.4、步骤S2中，轴向压力的实现方法为：

[0024] 在试样顶部生成加载板，加载板由三层颗粒组成，每层颗粒按六边形或四边形规则排列，颗粒间的接触本构模型采用线性平行粘结模型，平行粘结的刚度取为试样颗粒间
100
接触刚度的50‑100倍，平行粘结的强度设置为1.0×10 Pa，以此模拟不发生变形的刚性加
载板，将所需要的荷载转换为加载板顶层颗粒上的集中力进行施加，转换公式为：

[0025]

[0026] 其中，Fplaten1为加载板顶层每个颗粒上需施加的轴向集中力，σ为所设定的轴向压应力，A为虚拟压缩试样的横截面面积，Nplaten1为加载板顶层颗粒的数量。

[0027] 进一步地，本发明的所述步骤S1.3中，颗粒间摩擦系数与颗粒堆积体孔隙率的关系由若干数值试验测点通过指数函数 拟合得到；式中，n为孔隙率，μfric为
颗粒间摩擦系数，a、b、c为拟合参数。

[0028] 进一步地，本发明的所述步骤S2中，每级轴向压力按10kPa/5000steps的加载速度缓慢施加，以避免加载时动力作用的影响，且加载过程中试样的侧面和底部为刚性边界所
约束。

[0029] 进一步地，本发明的所述步骤S3.2中，切片图像处理的方法包括：

[0030] 切片图像处理过程中，将每张切片图像上的圆形墙边界向内缩进一与模型中最小颗粒半径相等的距离，使得边界与附近颗粒相交叠以免后续渗流模拟时在边界附近产生集
中渗流。

[0031] 进一步地，本发明的所述步骤S3.2中，切片图像处理的方法包括：

[0032] 切片图像处理过程中，对圆形边界以外、图片正方形边界以内的区域标记为固体区域；除了识别切片图像上的孔隙单元，还将固体区域进一步区分为流‑固边界单元和固体
内部单元两类，具体作法是：若某固体单元周围的邻近单元均为固体单元，邻近单元包括相
邻切片上对应的相邻单元，则该固体单元标记为固体内部单元；若有一个或多个邻近单元
为孔隙单元，即流体单元，则该单元标记为流‑固边界单元。

[0033] 进一步地，本发明的所述步骤S4.1中，具体的计算参数包括：

[0034] 流体粒子碰撞模型，采用单松弛时间BGK碰撞模型；离散速度模型，采用D3Q19模型；边界条件为，模型四周壁面和固体颗粒表面为无滑移的流‑固边界，采用反弹法进行处
理，流体流动方向沿试样轴向，由入口和出口的压力差驱动，压力边界采用Zou/He法进行处
理；

[0035] 所述步骤S4.1中，渗流计算过程中固体内部单元不执行迁移与碰撞，以减少计算耗时。

[0036] 本发明产生的有益效果是：本发明的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，原理简单、实施方便、计算稳定、效果良好，可以克服常规室内压缩渗
透试验的不足，具体表现在：1)渗流由格子Boltzmann方法模拟，可方便地获取流场中任一
点的水压力，而无需像室内试验那样需在试样侧面布置测压管，因此试样高径比可设为较
小值，以保证试样的压缩比较均匀。2)侧限压缩变形由颗粒离散元法模拟，由3层粘结在一
起的颗粒组成刚性加载板，通过加载板颗粒逐级缓慢施加轴向压力，可有效减轻动力作用
的影响。3)通过对压缩变形后试样的图像处理，使得其边界与附近颗粒相交叠，进而避免了
渗流模拟时在边界附近产生集中渗流。4)颗粒离散元与格子Boltzmann方法都是细观力学
模拟方法，基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟能方便地获取压缩
变形过程中试样的细观结构和对应的渗流场，进而可加深人们对岩土体压缩变形及变形过
程中渗透特性变化的内在机制的理解。

附图说明

[0037] 下面将结合附图及实施例对本发明作进一步说明，附图中：

[0038] 图1是本发明实施例的基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法流程图；

[0039] 图2是本发明实施例的虚拟压缩试样模型示意图；

[0040] 图3是本发明实施例的侧限压缩模型示意图；

[0041] 图4是本发明实施例的切片图像处理示意图；

[0042] 图5是本发明实施例的轴向压力‑孔隙率‑渗透率关系曲线图；

[0043] 图6是本发明实施例的流速场云图。

具体实施方式

[0044] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不
用于限定本发明。

[0045] 如图1所示，本发明公开了一种基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法，包括以下步骤：

[0046] S1:开展砂雨法与拟振动压实法的离散元模拟，制备不同相对密实度的虚拟压缩试样。具体步骤为：

[0047] S1.1:在比试样高一些的模型空间中按级配要求生成彼此不重叠的颗粒集合体，施加重力荷载，进行离散元法迭代计算，迭代中每隔一定计算循环将颗粒速度重新归零以
控制颗粒下落后的最大速度不超过0.8m/s，下落完成并计算至平衡后获得颗粒堆积体的最
松散状态；

[0048] S1.2:在颗粒堆积体的顶部施加14kPa的轴向压力，并设定颗粒间摩擦系数为0近似模拟最优振动作用，计算平衡后获得颗粒堆积体的最密实状态；

[0049] S1.3:维持颗粒堆积体顶部的14kPa轴向压力，从最松散状态开始逐步减小颗粒间的摩擦系数来模拟振动作用，获得颗粒间摩擦系数与颗粒堆积体孔隙率的关系，通过该关
系并结合最松散状态和最密实状态，估算出获得给定相对密实度堆积体所需的颗粒间摩擦
系数；

[0050] S1.4:在所需的颗粒间摩擦系数下对最松散堆积体进行14kPa轴向压力的压实即可获得给定相对密实度下的颗粒堆积体。

[0051] 如图2所示为基于三维颗粒离散元软件PFC3D 5.0并经上述步骤建立的相对密度为0(即最松散状态)、0.5和1(即最密实状态)的颗粒堆积体，模型是直径为500mm的圆柱形，
侧边和底部由墙所限制，颗粒直径为20‑40mm，颗粒数量为2333，颗粒堆积体的高度均低于
300mm。

[0052] S2:施加轴向压力，开展侧限压缩试验的三维颗粒离散元数值模拟。如图3所示，轴向压力是通过在试样顶部生成加载板来实现的，加载板由三层颗粒组成，每层颗粒按六边
形规则排列，颗粒间的接触本构模型采用线性平行粘结模型，平行粘结的刚度取为试样颗
100
粒间接触刚度的100倍，平行粘结的强度设置为1.0×10 Pa，以此模拟不发生变形的刚性
加载板，将所需要的荷载转换为加载板顶层颗粒上的集中力进行施加。

[0053] 本实施例对图2中相对密度为0.5的数值试样开展了侧限压缩试验，施加的轴向压力分别为：50kPa、100kPa、200kPa、300kPa、400kPa，每级轴向压力按10kPa/5000steps的加
载速度缓慢施加，在加载到预设的压力后继续迭代计算至平衡。在模型中5个不同位置布设
测量球并测得5个孔隙率值，取平均值作为该级压力下数值试样变形稳定后的孔隙率。

[0054] S3:对压缩后的数值模型进行切片，并对切片图像进行处理，建立对应的孔隙结构模型。具体步骤为：

[0055] S3.1:将每级压力下变形稳定后的圆柱形试样沿轴向用一系列等间隔的切面进行切割，导出各切面处的模型切片图像，包括侧面圆柱形墙在切片上的圆形边界及球形颗粒
在切片上对应的圆；

[0056] S3.2:对导出的切片图像进行批处理，将切片图像上每个像素的尺寸调整为与相邻切片间的竖向间隔相等，识别各切片上属于孔隙的像素(对应于格子Boltzmann模型中的
流体单元)，以此实现对离散元数值试样的体素化，并建立试样对应的孔隙结构模型。

[0057] 如图4所示为虚拟压缩试样切片图像的处理过程及建立的孔隙结构模型，为统一起见，本实施例对各级压力下变形稳定后的试样从底部开始沿轴向长为200mm的部分进行
了切片，沿轴向切片共200张(格子Boltzmann模型的立方体单元边长即为1mm)。将每张切片
图像上的圆形墙边界向内缩进10mm，使得边界与附近颗粒相交叠以免后续渗流模拟时在边
界附近产生集中渗流。并且，将圆形边界以外、图片正方形边界(其边长对应的实际长度为
600mm)以内的区域作为固体区域，并将固体区域进一步区分为流‑固边界单元和固体内部
单元两类，具体作法是：若某固体单元周围的26个邻近单元(包括上下相邻切片上对应的相
邻单元)均为固体单元，则该固体单元标记为固体内部单元；若有一个或多个邻近单元为孔
隙单元(即流体单元)，则该固体单元标记为流‑固边界单元。在后续渗流格子Boltzmann模
拟过程中对固体内部单元不进行任何计算，即可大幅节省计算时间。

[0058] S4:针对建立的孔隙结构模型开展渗流的格子Boltzmann数值模拟。具体步骤包括：

[0059] S4.1:将每级压力下变形稳定后的孔隙结构模型导入格子Boltzmann计算程序中，设定粒子碰撞模型、离散速度模型、流体边界条件和相关的计算参数，进行格子Boltzmann
迁移与碰撞的迭代计算。

[0060] S4.2:待渗流场收敛到稳态后，停止迭代计算，根据所施加的水压差利用达西定律计算该级压力下数值试样的渗透率。

[0061] 本实施例的渗流格子Boltzmann模拟是在开源的分布式计算平台PALABOS上采用8个计算线程进行并行计算所实现的。流体粒子碰撞采用了单松弛时间BGK碰撞模型。离散速
度模型为D3Q19模型。模型四周壁面和固体颗粒表面为无滑移的流‑固边界并采用反弹法进
行处理。流体流动方向沿试样轴向，由入口和出口的压力差驱动，压力边界采用Zou/He法进
‑7
行处理。弛豫时间τ取为0.8，压力差为1×10 (格子单位)。如图5所示为模拟所获得的压力‑
孔隙率‑渗透率关系曲线，图6所示为格子Boltzmann模拟所获得的流速场云图。可以看到，
上述基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法可以较好地模拟岩土
体的侧限压缩变形及变形过程中渗透性的变化特征，这也进一步验证了本发明的优势。

[0062] 应当理解的是，以上仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本领域的技术人员在本发明所揭露的技术范围内，可轻易想到的变化或替
换，都应涵盖在本发明的保护范围之内。