考虑膜效应的岩土体三轴试验的离散元建模及数值模拟方法与流程

本发明涉及岩土体室内常规三轴试验的模拟方法，尤其是三维离散元建模和模拟。

背景技术：

常规三轴试验是土工试验中非常基本且重要的实验，它的设备简单，操作方便，可以很好地研究岩土体的受力变形特征，获得一些可靠的强度参数如抗剪强度。然而室内试验耗时长，数据往往具有很大离散性，对于复杂宏观现象，其微观机制很难揭示。数值模拟是一种很好地揭示三轴试验微观机制的手段，离散元法自提出以来，由于其对非连续介质较好的模拟能力，能清晰反映颗粒材料微观变形机制等优点，迅速被引入对岩土体材料的研究中。

三轴试验模拟时，在模拟装置底座、不透水石等组件时，常规离散元模型的边缘较为粗糙，使得其表面总体摩擦力较大，难以有效模拟。边界条件的施加尤其重要。在施加侧向围压时，室内试验通常采用液体施加恒定压力，即给样品外部套上一层乳胶膜。

三轴试验离散元模拟中，在模拟侧向恒定围压时，通常的做法是采用伺服控制系统控制的刚性板，这种方法操作复杂，虽然可以控制侧向围压恒定，然而会抑制侧向的大变形，使得模拟结果有一定的误差，且在设置加载板时，为了避免边界之间的相互作用，总要留出一点空隙。

技术实现要素：

为解决上述问题，本方法采用团簇模型来构建三轴试验模拟器的底板装置，并提出一种三维情况下弹性膜的模拟方法。

该方法主要是用单层方形网格状排列的球形颗粒模拟弹性膜。膜的两端固定在顶、底板上。采用团簇模型使膜单元更加紧密。在数值模拟过程中，忽略底板与膜的相互作用。由此建立起可以模拟弹性膜边界条件下三维三轴试验模拟装置。

本发明所公开的考虑膜效应的岩土体三轴试验的离散元建模方法，包括如下步骤：

s1、制样容器由顶板、底板及刚性圆筒组成，构建和制样容器直径相同的单层圆形顶板、底板，以及内径和制样容器直径相同的刚性圆筒，以完成制样容器的建模；

s2、在制样容器模型内生成岩土体颗粒样品，作为构建岩土体样品的样品单元不断沉积至指定高度，以形成目标尺寸的岩土体样品；

s3、将给定的岩土体样品单元的材料属性赋给所生成的岩土体样品，然后平衡岩土体样品单元间的受力，以完成岩土体样品的建模；

s4、构建弹性膜，并将步骤s1所形成的制样容器模型中的刚性圆筒替换为构建好的弹性膜，然后给弹性膜设置相应地材料属性，至此完成三轴试验离散元模型的建立。

作为一种优选方案，步骤s1中，采用团簇模型构建制样容器模型顶板、底板和刚性圆筒。

作为一种优选方案，步骤s3中，平衡岩土体样品单元间的受力后，根根据模拟的岩土体样品是否存在初始胶结确定是否胶结样品单元。

作为一种优选方案，步骤s2包括以下步骤：

21、在顶板、底板和刚性圆筒围成的空间内采用随机法生成粒径在指定粒径分布范围内的岩土体颗粒样品；

22、对颗粒样品施加随机方向的初速度，通过离散元迭代运算使颗粒样品运动和碰撞至随机位置；

23、对颗粒样品施加重力加速度g，通过离散元迭代运算，使颗粒样品自由堆积，制样容器的底板固定，顶板自由且在重力作用下与颗粒样品一同下落，以压平颗粒样品堆积体表面，直至压到指定高度后岩土体样品建模初步完成。

作为一种优选方案，步骤s4中，采用团簇模型构建弹性膜。

作为一种优选方案，采用团簇模型构建弹性膜包括以下步骤：

41、将一层呈正方形网格状排列的弹性膜单元紧密堆积形成弹性膜，弹性膜的两端固定在顶板和底板上，即限制顶板和底板处的弹性膜单元水平方向的运动。

42、设置弹性膜的材料属性；

43、通过迭代运算平衡弹性膜单元间的受力，至此完成弹性膜的构建。

作为一种优选方案，步骤s4中，限制顶板和底板处的弹性膜的沿水平方向运动。弹性膜的两端分别固定在顶板和底板上，以限制固定在顶板和底板上的弹性膜沿水平方向的运动，仅允许其沿竖直方向运动。

本发明还公开一种考虑膜效应的岩土体三轴试验的数值模拟方法，该方法基于具有上述特征的任意一种三轴试验离散元模型进行数值模拟，根据给定的加压条件，采用应力式加载对岩土体样品进行轴向施压，通过水压力对岩土体样品施加侧向围压。

作为一种优选方案，采用应力式加载是通过改变顶板的重力对每一级加载，每级位移小于最小样品单元的直径，且每级加载后均要再次平衡模型内部应力。

作为一种优选方案，通过水压力对岩土体样品施加侧向围压是指，将水压力均匀的施加到每个弹性膜单元上以施加给定的围压，即将围压等效成膜单元受力，且每次施加围压后需再次平衡模型内部应力。

作为一种优选方案，计算数值模拟时弹性膜单元的受力包括如下步骤：

s51、对弹性膜进行网格划分，将弹性膜的面积赋给每个网格所对应的弹性膜单元；

s52、将划分后的每个网格上所受的围压等效到各网格中心所对应的弹性膜单元上，并计算得到每个弹性膜单元所对应的压力。

本发明能为室内三轴试验模拟提供了一种更符合真实试验的数值模拟方法，具有以下有益效果:

(1)区别于传统数值模拟中的伺服墙，本发明采用膜单元来模拟三轴试验中的乳胶膜，其柔性边界更加符合试验情况。

(2)将建模所用刚性圆筒替换为弹性膜，在模拟侧向恒定围压时，弹性膜上的水压力等效到样品周围的弹性膜单元上，使得围压比较稳定。

(3)在顶板和底板建模时，采用团簇模型模拟平滑样品的顶板和底板，使表面更光滑，从而减小了摩擦力导致的边界效应。

(4)在进行数值模拟时，可通过算法使顶底板与弹性膜不发生相互作用，即忽略底板和弹性膜的相互作用，这样可以采用大面积的底板可以防止模拟时颗粒溢出泄露，解决了以往岩土体离散元三轴模拟时需要留出缝隙导致样品颗粒可能发生泄露这一技术难题。

附图说明

图1为实施例中考虑膜效应的岩土体三轴试验建模和模拟方法的实现流程图。

图2为实施例中样品单元生成流程图。

图3为实施例中团簇模型原理示意图。

图4为实施例中弹性膜构建流程图。

图5为实施例中弹性膜的受力计算流程图

具体实施方式

图1为考虑膜效应的岩土体三轴试验数值模拟方法的实现流程图，包括离散元建模和数值模拟实验两大部分，具体如下：

首先，设置样品材料的物理力学性质和模型的几何参数等，包括岩土体样品(简称样品)的尺寸、材料属性，作为岩土体样品的组成单元(即样品单元)的颗粒样品的粒径分布范围，弹性膜的材料属性等。然后开始进行考虑膜效应的岩土体三轴试验离散元模型的构建，即建模，包括如下步骤：

步骤s1、制样容器建模。

制样容器由顶、底板及刚性圆筒组成。根据输入的岩土体样品的尺寸，采用团簇模型构建和样品直径相同单层圆形顶板和底板，同样使用团簇模型建立单层空心圆筒模型，如图1中的1a所示。

如图3中的3a所示，在团簇模型中，两个直径为d的单元相互重叠，设重叠量为l0时两单元处于平衡状态，单元间相对位移xn由以下公式计算得到：

xn＝r-(d-l0)

r表示两单元的圆心之间的距离，两单元的相对位移xn为零，处于平衡状态，中间箭头线条所在位置可认为两单元的接触表面。图3中的3c和3d为团簇模型的法向正应力响应示意图。如图3c，xn增加，单元间产生拉力。如图3d，同样的，当相对位移达到断裂位移xb时，连接断裂，连接断裂后，当两单元相互挤压时，单元间仍然有压缩力。事实上，团簇模型将两单元间的平衡距离缩小，在计算相对位移时减去其初始的重叠量。通过这种方法，可以实现较为复杂模型的建立。

如图3中的3b所示，四个单元相互重叠接触，将此时的距离视为平衡距离，则四个球处于受力平衡状态。可通过此方式构建表面较为光滑的团簇模型，并用来建立有效的顶、底板模型。

步骤s2、样品单元生成。

离散元的模型有一系列颗粒样品构成，根据输入的颗粒样品的粒径分布，在刚性圆筒和顶、底板围成的空间内按行列生成颗粒样品,即岩土体样品的样品单元。采用这种随机速度法移动样品单元使样品颗粒更加均一，然后通过重力沉积至指定尺寸，即样品所需的高度或尺寸。

如图2所示，步骤s2中样品单元生成的具体过程如下:

s21、生成颗粒样品。在刚性圆筒和顶、底板围成的空间内按行列生成粒径在指定粒径分布的颗粒。

s22、随机改变颗粒位置。对颗粒施加随机方向的一定初速度，通过离散元迭代运算，颗粒运动和碰撞至随机位置，使样品颗粒更加均一。

s23、重力堆积。对颗粒施加重力加速度g，通过离散元迭代运算，使颗粒自由堆积。圆筒底板固定，顶板自由且在重力作用下与颗粒一同下落，以压平颗粒堆积体表面。直至压到指定体积，样品单元建模初步完成。

s24、平衡模型。通过迭代运算平衡样品单元间受力，直至样品颗粒间加速度近似为0。

步骤s3、材料设置。

根据输入的样品材料参数设置样品模型的力学性质，如弹性模量、抗压强度、抗拉强度、泊松比、摩擦系数等；然后平衡样品单元间受力，平衡结束后，根据需要确定是否胶结单元，如根据模拟的样品是否存在初始胶结的特点确定是否胶结单元。

步骤s1至s3为样品制作环节，模拟的是取样制样过程。s4的作用为将原本的刚性圆筒替换成弹性膜以优化用于施加围压的侧向边界。

步骤s4、弹性膜建模。

将建模所用的刚性圆筒替换为弹性膜，并设置相应地材料属性。弹性膜由一层方形网格状排列的弹性膜单元紧密堆积而成,同样使用团簇模型弹性膜的两端固定在顶板和底板上，这也与真实的室内试验采取的措施一致。膜单元设置完成后再次平衡模型内部应力。至此，完成用于数值模拟试验的三轴试验离散元模型构建。

图4所示为步骤s4中弹性膜建模的流程图，包括如下步骤：

s41、生成弹性膜单元。将刚性圆筒去掉，用弹性膜替代。弹性膜由一层方形网格状排列的弹性膜单元(简称膜单元)紧密堆积而成,同样使用团簇模型。弹性膜的两端与顶、底板相交处，限制水平方向的运动(x和y方向)，仅允许竖直方向(z方向)的运动。

s42、材料设置。设置弹性膜的材料属性，如某种橡胶膜的弹性模量为表一所示。

表一某种橡胶膜的材料参数

s43、平衡模型。材料设置完成后，由于替换了侧向边界，微观参数发生了变化，整个模型受力不再平衡，因此，需要再次平衡弹性膜单元间的受力，即通过迭代运算平衡样品单元间受力，直至样品颗粒间加速度近似为0。

至此，考虑膜效应的岩土体三轴试验的离散元模型建立完成，接下来对基于该模型所采用的数值模拟方法进行相应的说明。

步骤s5、数值模拟实验。

在数值模拟实验中迭代计算时，需要根据模型中样品的受力以计算颗粒的加速度等信息，然后以牛顿运动定律更新颗粒的位置。实施例中，样品颗粒受到的轴向力为加载压力和自重，侧向的围压则由弹性膜提供。

根据给定的加压条件进行试验，轴向力采用应力式施加，通过改变顶板重力的方法使加压更加均匀。在每一级加载之后，进行迭代计算，平衡模型。采用应力式加载时，每级位移要远小于最小样品单元直径，即小于生成的最小颗粒样品的粒径。每级加载后均要充分平衡，当加载至达到指定应变或加压完成时，试验结束，输出试验结果。

为了施加给定的围压，将水压力均匀的施加到每个膜单元上，将围压等效成膜单元受力。结合图5所示，数值模拟时每个膜单元受力的计算方式包括如下步骤：

s51、划分网格，以将弹性膜的面积赋给对应膜单元。

弹性膜由全体膜单元组成，为将周围水压力施加到样品上，将样品周围的弹性膜进行网格划分，将每个网格上所受的围压等效到网格所属的弹性膜单元。将每个弹性膜单元周围相邻的4个弹性膜单元的球心组成四边形，将该四边形各边的中点连线得到一新的四边形，即图中所示的阴影部分，由此得到的四边形即为中心弹性膜单元对应的一个网格，通过这种方式将弹性膜区域划分成一系列对应中心单元的四边形网格。上下边界上的弹性膜单元对应的网格，在生成弹性膜单元时让两端略微超出顶、底板，这样就可以采用同样的方式对底板处的弹性膜划分网格，超出顶、底板的弹性膜单元不需要划分网格。

s52、等效水压力，计算得到每个膜单元所对应的压力。

为了施加给定的围压，将水压力均匀的施加到每个膜单元上，将每个网格上的水压力之和fω等效到网格中心的膜单元上，fω的方向平行于四边形所在面的法线方向并指向膜的内侧。每个膜单元所对应的fω可以由下式确定：

fω＝σω·se

其中，σω为给定的围压，se为弹性膜单元所占的弹性膜的面积。

尽管以上结合附图对本发明的实施方案进行了描述，但本发明并不局限于上述的具体实施方案和应用领域，上述的具体实施方案仅仅是示意性的、指导性的，而不是限制性的。本领域的普通技术人员在本说明书的启示下，在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。