一种用于精确计算离散元法模拟三轴试验柔性膜体积应变的方法

技术领域

本发明涉及岩土数值模拟技术领域，尤其是一种用于精确计算离散元法模拟三轴试验柔性膜体积应变的方法。

背景技术

随着计算机技术的普及，通过数值模拟方法研究具有离散特征的岩土体材料物理力学特性受到越来越多的科研工作者的关注。相较于实验室试验周期长，耗时久，人员多，开销大等弊端，数值模拟方法可以大大降低人工成本和试验周期。

离散元是一种常用于岩土领域的数值模拟方法，它通过考虑颗粒间的接触来研究岩土散体材料的物理力学特征，目前在岩土力学方面取得了越来越多的应用。特别是离散元法模拟三轴试验方面，取得了重要的进展。常见的通过离散元法模拟三轴试验的方法包括圆柱边界法，REV局部法，以及膜边界法等。这三种方法各有利弊。圆柱边界法将试验变形近似为圆柱形，计算简单，但难以模拟试样的不规则变形；REV局部法通过试样内部特征区域的局部行为表征试样的整体行为，可以大大减少计算量，但是局部和整体毕竟有差距；膜边界法可以显示剪切过程中样品的不规则变形，与其他边界条件相比，膜边界条件最能体现试样的实际变形情况，但仍有一个关键问题没有得到根本解决，即试样的不规则体积变形的精确计算。由于膜边界的优势在于反应试样在模拟过程中的体积变形，因此变形是不规则的，难以精确跟踪过程中的变化，见图1。目前，常用的计算膜边界体积应变的方法仍为近似法，代表性的膜边界体积应变计算方法包括第一种根据试样半径的变化来计算体积变化。通过确定边界粒子和样本中心线之间的水平距离，获得平均半径，再通过获取试样当前的高度，计算圆柱体的样本体积，该方法本质仍在体积计算时把试样不规则体积变形近似为圆柱形。第二种方法通过平均孔隙比和初始体积来评估体积变化，但是平均孔隙比是通过试样局部REV来获得的，无法如实地代表试样的整体情况。因此，这两种方法尤其是在试样在变形较大的情况下的误差不能忽略。

发明内容

本发明所要解决的技术问题在于，提供一种用于精确计算离散元法模拟三轴试验柔性膜体积应变的方法，可以精细计算试样在三轴剪切过程中任意不规则体积形变。

为解决上述技术问题，本发明提供一种用于精确计算离散元法模拟三轴试验柔性膜体积应变的方法，包括如下步骤：

步骤1、生成膜边界；

步骤2、在膜边界内生成互不接触的球形颗粒或不规则形状颗粒，直到空间内的整体颗粒达到试样预设的孔隙率；

步骤3、根据岩土体离散材料属性，针对试样中颗粒与颗粒之间、颗粒与立方体内壁之间的各个接触位置，分别设置接触模型；

步骤4、在试样的固结和剪切阶段，围压作用在膜边界上的力由围压与膜表面的乘积计算；

步骤5、将试样的体积分为三部分，即试样的上表面或下表面与中心点连接而形成的两个圆锥体，以及圆柱体的其余部分；

步骤6、在每个时间步长重复步骤4和步骤5即可时时记录膜边界的不规则体积变形。

优选的，步骤1中，膜边界条件是用规则排形的球颗粒组成，即膜颗粒，膜颗粒采用球颗粒。

优选的，步骤2中，采用颗粒膨胀法生成互不接触的球形颗粒或不规则形状颗粒。

优选的，步骤2中，针对立方体堆积体，采用夯实试样或振动试样的方式处理，制备满足预设孔隙率的试样。

优选的，步骤2中，针对立方体堆积体的顶部，采用上压板针对立方体堆积体的顶部进行击打，实现夯实试样方式处理，制备满足预设孔隙率的试样。

优选的，步骤2中，针对立方体堆积体的底部位置设置下底板，控制下底板以默认100Hz的频率振动，通过振动试样方式处理，制备满足预设孔隙率的试样。

优选的，步骤3中，接触模型包括线性接触模型，Hertz模型，平行粘接模型，抗滚动模型。根据不同的材料属性进行选择，进行简单估算可选用线性接触模型或Hertz模型，若有粘结的散体材料可选用平行粘结模型，若颗粒不规则形状较为明显，则选用抗滚动模型。

优选的，步骤4中，每三个相邻的膜颗粒的中心形成一个三角形，膜的表面积对所有三个相邻膜颗粒的中心形成的小三角形的求和；然后，通过膜颗粒和材料之间的接触，将膜上的约束力传递到试样上，在样品被大幅剪切后，相邻膜颗粒组成的三角形并不总是等边的。

优选的，计算三角形面积具体包括如下步骤：

步骤41、假设三角形三点坐标分别为P

步骤42、三角形面积计算公式为：

F＝p×A。

优选的，步骤5中，第三部分的体积是所有三维单纯形之和，三维单纯形为三个相邻膜颗粒的三角形和样品中心形成的体积，计算第三部分的体积具体包括如下步骤：

步骤51、圆锥体体积

其中A为上下底面积，H为上下底面分别到质心的距离；

步骤52、区域三内单个三维单纯形体积：

假设某三维单纯形四点坐标

P

其中J为雅可比行列式：

所有三维单纯形体积为

V

步骤53、总体积计算公式：

Vtotat＝Vcone1+Vcone2+Vrest。

本发明的有益效果为：本发明可以时时精确计算任意体积变形的试样体积应变，解决了基于离散元法的三轴试验膜边界试样模拟剪切形为的痛点问题，通过精确计算试样的体积形变有利于完善数值模型的精确程度，实现虚拟实验对室内实验的复制，模拟材料宏观属性，获取模拟材料的膨胀角、特征角等参数。

附图说明

图1为本发明膜边界不规则形变示意图。

图2为本发明的膜边界及膜颗粒分布示意图。

图3为本发明的体积区域划分及三维单纯形示意图。

具体实施方式

一种用于精确计算离散元法模拟三轴试验柔性膜体积应变的方法，包括如下步骤：

步骤A：生成膜边界。如图1所示，膜边界条件通常是用规则排形的球颗粒组成，简称膜颗粒。在初始阶段是圆柱形，与试样一致。散体材料填充在膜边界内。随着剪切的过程，膜边界上每个膜颗粒受力发生有限自由运动，由于膜颗粒不是统一行动，因此从整体看，膜变形即体现了试样的不规则形状。膜边界条件由相同的重叠膜颗粒组成，高度分布189个圆，每个圆201个颗粒，总共37989个膜颗粒，如图2(a)所示。膜颗粒的数量越大，越能精确表征试样的不规则变形细节，但是计算时间也会相应增加。因此，膜粒径的选择既考虑了计算效率，又考虑了模型精度。膜粒子的接触模型是线性接触键模型。然后进入步骤B；

步骤B：在膜边界内生成互不接触的球形颗粒或不规则形状颗粒,直到空间内的整体颗粒达到试样预设的孔隙率，然后进入步骤C；

采用颗粒膨胀法生成互不接触的球形颗粒或不规则形状颗粒。针对立方体堆积体，采用夯实试样或振动试样的方式处理，制备满足预设孔隙率的试样。针对立方体堆积体的顶部，采用上压板针对立方体堆积体的顶部进行击打，实现夯实试样方式处理，制备满足预设孔隙率的试样。针对立方体堆积体的底部位置设置下底板，控制下底板以默认100Hz的频率振动，通过振动试样方式处理，制备满足预设孔隙率的试样。接触模型包括线性接触模型，Hertz模型，平行粘接模型，抗滚动模型。

步骤C：根据岩土体离散材料属性，针对试样中颗粒与颗粒之间、颗粒与立方体内壁之间的各个接触位置，分别设置接触模型，然后进入步骤D；

步骤D：在试样的固结和剪切阶段，围压作用在膜边界上的力由围压与膜表面的乘积计算。每三个相邻的膜颗粒的中心形成一个三角形，如图2(b)所示。膜的表面积对所有三个相邻膜颗粒的中心形成的小三角形的求和。然后，通过膜颗粒和材料之间的接触，可以将膜上的约束力传递到试样上。在样品被大幅剪切后，相邻膜颗粒组成的三角形并不总是等边的，如图2(c)所示。

步骤D-1：

假设三角形三点坐标分别为P

步骤D-2：

三角形面积计算公式为：

步骤E：将试样的体积分为三部分，即试样的上表面或下表面与中心点连接而形成的两个圆锥体，以及圆柱体的其余部分，如图3所示。其中，第三部分的体积是所有三维单纯形(三个相邻膜颗粒的三角形和样品中心形成的体积)之和。

步骤E-1：圆锥体体积

其中A为上下底面积，H为上下底面分别到质心的距离。

步骤E-2：区域三内单个三维单纯形体积：

假设某三维单纯形四点坐标

P

其中J为雅可比行列式：

所有三维单纯形体积为

V

步骤E-3：

总体积计算公式：

V

步骤F：在每个时间步长重复步骤D和步骤E即可时时记录膜边界的不规则体积变形。