一种岩土三轴试验柔性薄膜边界的离散元模拟方法与流程

本发明属于岩土工程领域，尤其涉及一种岩土三轴试验柔性薄膜边界的离散元模拟方法。

背景技术：

三轴试验属于岩土工程领域最为重要的室内试验之一。在工程中，三轴试验广泛用于测试岩土体强度和变形参数，为设计和施工提供参数依据；在科研中，三轴试验则是研究土体各类力学特征的主要设备。

由于砂土，堆石料等土工材料具有天然的散粒特性，以颗粒为基本单元的的离散元算法广泛应用于各类岩土工程问题的研究中。然而，由于离散元算法采用的细观尺度参数不易通过物理试验测得，当前大多数针对岩土工程的离散元模拟，首先都要模拟常规岩土试验(如三轴试验)，通过不断调试离散元细观参数，直到模拟对象在宏观上具备和物理试验相似特征时，所用的离散元参数才被当作是一组可靠的参数。因而，基于离散元的三轴试验作为离散元模拟岩土工程问题的基础数值试验，其准确仿真是离散元模拟岩土工程问题的关键之一。

三轴试样在侧向通常由乳胶薄膜或者橡胶薄膜包裹在一定水压环境中，通过轴向加载，获得试样在一定围压下轴向加载时的应力应变关系。柔性的乳胶膜或橡胶膜，能在保证试样加载过程中自由变形的同时，稳定地传递水压作为试样的围压。有研究表明，薄膜的性质对于试样加载过程中的应力应变响应，以及剪切带的产生和发展具有重要影响。

中国发明专利(申请号：201710636809.7，专利名称：考虑膜效应的岩土体三轴试验的离散元建模及数值模拟方法与流程)提供了一种三轴试验薄膜模拟算法，该算法采用团簇单元模拟顶板，底板和薄膜，能在一定程度上模拟三轴试验加载后的变形。然而，该专利提供的模拟方法存在明显不足：1)无法模拟不同薄膜材料参数(如厚度，弹性模量和泊松比等)在三轴中的影响；发明人虽然提及要对薄膜输入弹性模量和泊松比等变形参数，然而参与离散元计算的参数均为颗粒刚度等细观参数，如何将宏观的弹性模量和泊松比等宏观参数赋与采用细观参数计算的团簇单元，存在现实的技术困难。2)合理捕捉三轴加载过程中的体积变形是成功分析岩土体变形行为的关键，然而该模型无法计算和合理表征三轴试验在加载过程中的体积变形；3)采用团簇单元作为上下加载顶板的模拟单元，由于上下加载顶板刚度大，团簇单元必须采用较高的刚度以合理模拟刚性加载板，然而，在显式离散元算法中，较大的刚度不得不采用较小的时间步长来保证模型计算的稳定性，该方案存在具体实施上的计算低效率。

技术实现要素：

本申请旨在至少解决现有技术中存在的技术问题之一。为此，本发明的目的之一在于提供一种岩土三轴试验柔性薄膜边界的离散元模拟方法。该方法能准确模拟散体材料在三轴试验过程中试样形态变化，应力应变关系和体积应变情况，能为散体材料的宏观参数测试和力学性能研究提供有力支撑。

为解决上述技术问题，本发明采用如下技术方案：

s1、制样容器建模；

制样容器模型包括上加载板、下加载板和侧向刚性墙体；

s2、岩土体样品的建模；

在制样容器内生成土体颗粒样品，基于离散元伺服原理，以制样容器为伺服墙体，将三轴试样加载到预定的应力状态；

s3、根据试验中所用薄膜的弹性模量、泊松比和薄膜厚度，计算构建薄膜模型所需的薄膜颗粒的细观参数：

其中：kn和ks分别为薄膜颗粒的法向刚度和切向刚度；e，ν和t分别为试验薄膜的弹性模量，泊松比和厚度；

s4、删除制样容器模型中侧向刚性墙体，将薄膜颗粒以六边形排列形式单层附着在岩土体样品的外围，薄膜颗粒间采用线性粘结模型粘连，进而完成六边形蜂窝结构的圆筒形薄膜模型的构建，将薄膜模型的两端分别固定在上加载板和下加载板上，至此完成三轴试验离散元模型的建立；

s5、施加静水围压；

s6、记录三轴加载过程中的体积应变；

s7、沿轴向实施应变式加载。

进一步的，所述步骤s5包括以下步骤：

1)将薄膜颗粒的速度固定为0；

2)计算整个离散元模型，并迭代至模型平衡，对所有岩土试样颗粒速度清零；

3)施加作用力到薄膜颗粒上，对于六边形颗粒排列来说，整个颗粒薄膜可以看作一系列三角形网络构成，每个三角形单元上承受的静水围压，都被构成此三角形单元上的三个颗粒均匀承担，则作用于每个颗粒上的等效静水围压作用力f⁰：

其中：σstatic为静水围压的大小；ni是第i个三角形单元的外法线方向，ni＝(n1，n2，n3)，si为第i个三角形单元的面积；每隔一定时步更新每颗薄膜颗粒应分配得到的等效薄膜力，可使得薄膜在自由变形过程中保持合理的静水围压；

4)释放位于上下加载板之间的薄膜颗粒速度；

5)将模型迭代计算至平衡状态。

进一步的，三轴试样在加载过程的体积，按如下公式计算：

其中：vs为试样体积，si为第i个三角形单元的面积，x1为第i个三角形单元中心在x轴方向上的坐标值，n1为第i个三角形单元外法线方向与x轴的方向余弦。

进一步的，施加静水围压时，每隔一定时步更新每颗薄膜颗粒应分配得到的等效静水作用力。

进一步的，薄膜模型的内半径为薄膜颗粒半径的35-100倍。

本发明与现有技术相比，本发明的技术效果在于：(1)能用细观尺度的颗粒刚度有效表征物理试验中的薄膜参数(弹性模量，泊松比，薄膜厚度)，合理描述不同类型薄膜对三轴试验模拟的影响；(2)薄膜在大变形过程中也能确保稳定和准确地加载静水围压；(3)能准确捕捉三轴试样在加载变形过程中体积的体积变形。

综上所述，本发明是一种优点突出、实施流程清晰且简单实用的薄膜模拟方法，可为基于离散元算法准确模拟岩土三轴试验提供技术支撑。

附图说明

图1为三轴试验初始试样制样容器；

图2为颗粒薄膜圆筒与上下加载板示意图；

图3为颗粒薄膜六边形蜂窝结构布置与代表性单元图；

图4为整体与局部坐标系示意图；

图5为三轴试验结束后的试样形态图(a：室内试验，b：刚性墙边界以及c：柔性薄膜边界)。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种岩土三轴试验柔性薄膜边界的离散元模拟方法，具体包括以下步骤：

步骤(1)：建立制样容器，制样容器模型包括上加载板、下加载板和侧向刚性墙体；

步骤(2)：岩土体样品的建模：在制样容器内生成土体颗粒样品，基于离散元伺服原理，以制样容器为伺服墙体，将三轴试样加载到预定的应力状态。

步骤(3)：计算薄膜颗粒的刚度参数。根据薄膜的弹性模量和泊松比，按如下公式计算薄膜颗粒的细观参数：

其中：kn和ks分别为薄膜颗粒的法向刚度和切向刚度；e，ν和t分别为试验薄膜的弹性模量，泊松比和厚度；具体推导过程见后文。

步骤(4)：安装颗粒薄膜。此步骤可分为以下4个子步骤：

1)删除侧向刚性墙体；

2)如图2所示，将薄膜颗粒以六边形排列形式单层附着在岩土体样品的外围，薄膜颗粒间采用线性粘结模型粘连，形成六边形蜂窝结构的圆筒形颗粒模型的构建，薄膜圆筒的半径应比侧向刚性墙体稍大，颗粒模型半径约为薄膜颗粒半径的35-100倍之间，颗粒间采用线性粘结模型，粘结强度需要足够大，使得颗粒薄膜的结构在加载过程中不被破坏。3)抑制上下刚性加载板与薄膜颗粒之间的相互作用，保留上下刚性加载板与三轴试样颗粒的相互作用。

步骤(5)：施加静水围压。此步骤可分为以下5个子步骤：

1)将薄膜颗粒的速度固定为0；

2)计算整个模型并迭代至平衡，对所有试样颗粒速度清零；

3)施加作用力到薄膜颗粒上(对应于试验中的施加静水围压。如图3a所示，对于六边形颗粒排列来说，整个颗粒薄膜可以看作一系列三角形网络构成，假设每个三角形单元上承受的静水围压，都被构成此三角形单元上的三个颗粒均匀承担，则作用于每个颗粒上的等效静水围压作用力f⁰，等于该颗粒与其周围六个颗粒分别构成的6个三角形区域上所受静水围压的矢量和。即：

其中：σstatic为静水围压的大小；ni和si分别是第i个三角形单元的外法线方向和面积；此公式可将围压的大小和方向准确加载到每颗薄膜颗粒上，每隔一定时步更新每颗薄膜颗粒应分配得到的等效薄膜力，可使得薄膜在自由变形过程中保持合理的静水围压。

4)释放位于上下加载板之间的薄膜颗粒速度，并保持固定位于上加载板以上位置和下加载板一下位置的颗粒速度；

5)将模型迭代计算至平衡状态。

步骤(6)：记录三轴加载过程中的体积变化。散体试样的剪切变形过程必然伴随着试样体积的变化。在本模型中，薄膜颗粒被设置为六边形布置，每相互相邻的三个颗粒总形成一个三角形，称此类三角形为一个单元。由于散体试样被柔性薄膜和上下加载板围成一个封闭的空间，因此，可考虑采用高斯散度定理进行试样体积的计算：

其中：vs为试样体积，s为封闭空间(颗粒薄膜与上下加载板)的表面，dv和ds分别为试样体积微分和薄膜表面的面积微分；si为第i个单元的面积，x1为第i个单元中心在x轴方向上的坐标值，n1为第i个单元外法线方向与x轴的方向余弦(即单元外法线方向与x轴夹角的余弦)

步骤(7)：通过赋予上下加载板稳定的移动速度，沿轴向对试样实施应变式加载。注意加载速度应尽量接近试验速度，以使得仿真模型处于准静态条件之下。当轴向应变达到预定应变值时，停止加载，三轴试验结束。

附：颗粒薄膜的细观参数公式推导过程

对于六边形型式的颗粒平面(如图3a),可选择其中的代表性单元(rve)为如图3b所示的六边形单元(单元0)，每个单元与周围6个单元(单元1-6)相邻。选择其中某两个颗粒为a和b，如图2所示，建立局部笛卡尔坐标系(x1′,x2′和x3′)，其中x2′方向连接颗粒a和b的中心；x1′方向与总体坐标系x1方向共轴。注：张量计算中，笛卡尔坐标系下的x轴，y轴和z轴，通常采用下标1,2,3表示；张量x的分量记法为xi(一阶)或者xij(二阶)，其中i和j可能是1,2,3中任意一项。

(1)应变能计算：

假设每对颗粒接触储存的能量由相互接触的两个颗粒均分，颗粒间接触力与接触位移为线性弹性关系，则储存在一个代表性单元中的能量为：

其中：kn和ks分别为颗粒的法向刚度和切向刚度；δu′n和δu′s分别为颗粒间的法向位移和切向位移。

(2)颗粒位移与应变的关系：

其中：和分别为颗粒a和b沿xj方向的坐标，为颗粒a与b之间的等效应变。

由于局部坐标系中x′2方向连接颗粒a和b的中心，所以：

其中：lab为颗粒a与b之间的距离，l2j为局部坐标系中x′2轴与整体坐标系中xj轴夹角的余弦。

为将整体坐标系中的相对位移δui和局部坐标系中的相对位移δu′i联系起来，需要采用如下坐标变化：

其中：lij局部坐标系中x′i与xj的余弦，即：

lij＝cos(e′i,ej)(8)

其中：e′i为局部坐标系中x′i轴的基向量，ej为整体坐标系中xj的基向量。

在图2所示坐标系下，a和b颗粒间的法向相对位移可以表示为：

a和b颗粒间的切向相对位移可以表示为：

将公式8和9进入公式3，则颗粒间的应变能可以表示为：

其中：lc,和u^c分别为接触距离，等效应变以及每个接触间的等效应应变能。

(3)应变能密度计算：

每个代表性单元(rve)的面积为：

假设薄膜的厚度恒为t,则每个代表性单元的体积为：

每两个相邻颗粒的距离，为颗粒半径的两倍(lc＝2r)：则颗粒薄膜的应变能密度为：

(4)应力张量与刚度张量计算：

根据弹性力学理论，应力张量可通过对应变能密度函数取相应应变张量的偏微分而获得，即：

弹性刚度矩阵，可通过对应力张量取相应应变张量的偏微分而获得，即：

(5)薄膜颗粒刚度参数的推导：

因为薄膜的半径远大于组成薄膜的颗粒半径，所以组成代表性单元的六个颗粒可看作在同一平面内。在附图4所示的坐标系下，局部坐标系可通过保持整体坐标系中x1轴不变，对x2和x3旋转相同的角度来获得。从整体坐标系到局部坐标系的方向余弦lij变为：

对于一般弹性体，其刚度矩阵可写为：

代入公式(16)和(17)进入(18)，可得：

假设薄膜在平面x2-x3内各向同性，则颗粒薄膜可看作是一种横观各向同性材料，根据广义胡克定理，其应力应变关系应满足:

其中：e1和e2为薄膜的弹性模量，ν1和ν2为薄膜的泊松比，g1和g2为薄膜的剪切模量。由于颗粒薄膜被假定其厚度在变形过程中保持不变，即e1＝∞：公式(20)变为：

其中，弹性常数e＝e1＝e2，ν＝ν1＝ν2，对应于各向同性平面内的弹性模量和泊松比。取出ε11和σ11项后，刚度矩阵变为：

通过对比(19)和(22)中的刚度矩阵，可得：

联立(23)和(24)可得出薄膜的颗粒刚度与弹性模量和泊松比的关系为：

以下将结合具体实施例对本申请作详细说明。

实施例

本实施例三轴试验参数如下：

参见图1，模拟步骤如下：

步骤(1)：生成初始试样。根据试验参数，在刚性圆筒墙体内生成试样，上下加载板为刚性墙体；使用伺服机制，将三轴试样加载到100kpa。

步骤(2)：计算薄膜颗粒的刚度参数。试验中所用薄膜的弹性模量，泊松比和颗粒厚度，计算薄膜颗粒的细观参数：

其中：kn和ks分别为薄膜颗粒的法向刚度和切向刚度；e、ν和t分别为试验薄膜的弹性模量、泊松比和厚度。

步骤(3)：安装颗粒薄膜。

1)删除侧向刚性墙体；

2)如图2所示，按照六边形排列型式生成颗粒圆筒，薄膜圆筒的半径应比侧向刚性墙体稍大，圆筒半径约为薄膜颗粒半径的40倍；颗粒间采用线性粘结模型，粘结强度为1e100kpa，使得颗粒薄膜的结构在加载过程中不被破坏。

3)抑制上下刚性加载板与薄膜颗粒之间的相互作用。

步骤(4)：施加围压。此步骤可分为以下5个

子步骤：

1)将薄膜颗粒的速度固定为0；

2)计算整个系统并迭代至平衡；

3)施加作用力到薄膜颗粒上(对应于试验中的施加静水围压。如图3a所示，对于六边形颗粒排列来说，整个颗粒薄膜可以看作一系列三角形网络构成，假设每个三角形单元上承受的静水围压，都被构成此三角形单元上的三个颗粒均匀承担，则作用于每个颗粒上的等效静水围压作用力，等于该颗粒与其周围六个颗粒分别构成的6个三角形区域上所受静水围压的矢量和。即：

其中：σstatic为静水围压的大小；ni和si分别是第i个单元的外法线方向和面积；此公式可将围压的大小和方向准确加载到每颗薄膜颗粒上，每隔一定时步更新每颗薄膜颗粒应分配得到的等效薄膜力，可使得薄膜在自由变形过程中保持合理的静水围压。

4)释放位于上下加载板之间的薄膜颗粒速度，并保持固定位于上加载板以上位置和下加载板一下位置的颗粒速度；

5)将系统迭代计算至平衡状态。

步骤(5)：记录三轴加载过程中的体积应变。采用高斯散度定理进行试样体积的计算：

其中：vs为试样体积，si为第i个单元的面积，x1为第i个单元中心在x轴方向上的坐标值，n1为第i个单元外法线方向与x轴的方向余弦。

步骤(6)：通过赋予上下加载板稳定的移动速度，对试样实施应变式加载。加载速度取5e-4m/s，以使得仿真模型处于准静态条件之下。当轴向应变达到16％时，停止加载。三轴试验结束时的室内试验试样变形情况、以刚性墙为边界的离散元试样变形情况和以柔性薄膜为边界的离散元试样变形情况分别如图5a﹑5b和5c所示，从5a﹑5b和5c可以看出本专利提出的柔性薄膜边界算法，相比传统的刚性墙边界算法，更能准确地表征三轴试验过程中的试样变形情况。

上述实施例仅仅是清楚地说明本发明所作的举例，而非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里也无需也无法对所有的实施例予以穷举。而由此所引申出的显而易见的变化或变动仍处于本发明的保护范围之中。