CN110162934B - 随机颗粒垫层上土工膜局部变形三维仿真模型建立方法 - Google Patents

Abstract

本发明公开了一种随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，步骤1，绘制颗粒垫层区域；步骤2，确定垫层区域中不同粒径圆球颗粒的个数；步骤3，随机生成最大粒径为d₁的球；步骤4，随机生成d₂、d₃、…、d_n粒径的球；步骤5，生成随机颗粒垫层；步骤6，生成土工膜和步骤7，土工膜局部变形计算与破坏判定。本发明能使三维仿真模型建立速度快，计算工作量小，随机垫层颗粒的模拟效果更为逼真，颗粒垫层与土工膜之间的接触更接近真实状态，土工膜的局部变形计算结果更为准确，计算与判定结果对指导土工膜上覆荷载大小、土工膜厚度以及颗粒粒径级配的合理选择具有重要的科学研究意义和工程应用价值。

Description

随机颗粒垫层上土工膜局部变形三维仿真模型建立方法

技术领域

本发明涉及水利工程中土工膜变形数值模拟领域，特别是一种随机颗粒垫层上土工膜局部变形三维仿真模型建立方法。

背景技术

土工膜因具有防渗性能好、适应变形能力强、工程造价低及施工速度快等优点，已广泛应用于大坝、库盘、蓄水池、垃圾填埋场等防渗工程。土工膜防渗结构的设计通常需要土石颗粒材料作为膜下垫层支撑体。颗粒垫层上土工膜在上覆荷载作用下在颗粒空隙间容易发生局部变形，一旦超过允许值土工膜就会发生破坏，对工程防渗安全造成隐患。由于土工膜下垫层颗粒是不规则、不均匀的，且排布具有一定的随机性，采用真实颗粒材料作为土工膜下垫层进行试验时，土工膜在上覆荷载作用下产生局部变形的部位是无法预先判断的，给土工膜变形的测量带来了极大的困难。目前，尚无相关试验设备和方法可直接测试上覆荷载作用下真实颗粒垫层土工膜的工作性态，同时没有可直接用于计算随机颗粒垫层上土工膜局部变形的解析模型。因此，通过构造随机三维圆球颗粒垫层土工膜变形分析的有限元仿真模型，不失为研究土工膜局部变形及破坏机理的一种有效方法。

目前，土工膜下垫层颗粒的有限元数值模拟主要分为两类：一是生成三维均匀圆球颗粒垫层；二是在平面内生成二维圆形颗粒垫层，现有技术存在如下不足：

1.传统的三维有限元模拟方法通常将垫层颗粒假定为同一粒径的圆球，生成大小均匀的三维垫层颗粒。然而，实际工程中土工膜下垫层颗粒粒径是不均匀的，且排布是随机的，土工膜受到大小不一的颗粒约束后，在不同颗粒间隙内发生非均匀局部变形。因而，现有技术中的简化处理方式，将会导致土工膜在不同颗粒间隙的局部变形完全相同，这与实际情况不符。

2.现有技术在平面内生成二维圆形颗粒模型计算土工膜在颗粒间的局部变形时，只能反映相邻两个圆形颗粒的约束作用。而实际工程中，颗粒间隙周围多个颗粒的约束作用均会对土工膜的局部变形产生影响，因此将土工膜下垫层颗粒二维简化处理会导致土工膜局部变形计算失真。

3.现有技术生成的二维垫层颗粒模型时，往往是大小不同的各圆形颗粒的中心线在一个高度，使得垫层表面不在同一高度，而实际工程中颗粒垫层经过碾压平整后，其表面保持在同一高度，从而与土工膜紧密贴合。因此，现有二维垫层颗粒模型不能反映土工膜与颗粒垫层的真实接触关系。

发明内容

为了解决上述现有技术的不足，本发明提供了一种随机颗粒垫层上土工膜局部变形三维仿真模型建立方法。该方法能使三维仿真模型建立速度快，计算工作量小，随机垫层颗粒的模拟效果更为逼真，颗粒垫层与土工膜之间的接触更接近真实状态，土工膜的局部变形计算结果更为准确，计算结果对指导土工膜上覆荷载大小、土工膜厚度以及颗粒粒径级配的合理选择具有重要的科学研究意义和工程应用价值。

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

一种随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，包括如下步骤。

步骤1，绘制颗粒垫层区域：在xoy平面内绘制水平放置且边长为L的正方形，并以正方形的中心为原点，以过原点且平行于正方形相邻两边的直线为x轴和y轴，以过原点且垂直于x、y轴的直线为z轴，建立xyz坐标系。

步骤2，确定垫层区域中不同粒径圆球颗粒的个数，包括如下步骤：

步骤21，确定圆球颗粒粒径和每种粒径所占质量百分比：颗粒垫层中的随机颗粒通过圆球模拟，根据垫层颗粒材料的级配曲线，获得不同圆球颗粒粒径和每种粒径所占的质量百分比。

步骤22，确定不同粒径圆球颗粒的个数：根据步骤21中确定的每种粒径所占质量百分比可求得不同粒径圆球颗粒的个数；假设不同圆球颗粒粒径从大到小依次为d₁、d₂、…、d_n，每种粒径对应圆球的个数依次为N₁、N₂、…、N_n，则圆球个数的计算公式如下：

其中，N_i为第i种粒径颗粒对应圆球的个数；L为正方形的边长；p_i为第i种粒径颗粒所占的质量百分比；d_i为第i种粒径颗粒的粒径。

步骤3，随机生成最大粒径为d₁的球，包括如下步骤：

步骤31，随机生成第一个d₁粒径的球：在步骤1绘制的颗粒垫层区域内，先随机生成第一个d₁粒径的球心坐标[x₁,y₁,0]，并使x₁和y₁均在[(-L+d₁)/2,(L-d₁)/2]内，否则，重新生成球心坐标。

步骤32，划分边长为d₁的正方形区域网格：作步骤31生成的d₁粒径球的x向和y向四条切线，以d₁粒径为边长向四周划分网格，网格的总数为

将d₁粒径球的周围每一个网格编号并计算每个网格中心点的坐标。

步骤33，随机生成第二个d₁粒径的球：随机选取一个步骤32中的网格编号，在选取的网格内随机生成第二个d₁粒径球的球心坐标，并使第二个d₁粒径球的球心的x坐标或y坐标均在[(-L+d₁)/2,(L-d₁)/2]内，且任意两球之间无重叠；否则重复步骤33。

步骤34，随机生成剩余d₁粒径的球：重复步骤33，直至第N₁个d₁粒径的球全部随机生成。

步骤4，随机生成d₂、d₃、…、d_n粒径的球：将步骤3中粒径d₁的球依次更换为d₂、d₃、…、d_n粒径的球，重复步骤31至步骤34，即可随机生成d₂、d₃、…、d_n粒径的球。

步骤5，生成随机颗粒垫层：将N₂个d₂粒径的球沿z轴向上平移(d₁-d₂)/2，将N₃个d₃粒径的球沿z轴向上平移(d₁-d₃)/2，依此类推，将N_n个d_n粒径的球沿z轴向上平移(d₁-d_n)/2，使d₁、d₂、…、d_n粒径球的顶部均保持平齐。

步骤6，生成土工膜：在d₁粒径球心所在平面正上方d₁/2处，生成与边长为L的正方形区域大小相同的平面，将生成的平面沿z轴向上拖放h厚度，其中，h为土工膜厚度，即得土工膜；随机颗粒垫层和土工膜共同构成三维仿真模型。

还包括步骤7，土工膜局部变形计算与破坏判定，具体包括如下步骤：

步骤71，将三维仿真模型导入有限元计算软件：将步骤6形成的三维仿真模型导入有限元计算软件ADINA中，并对土工膜进行网格剖分。

步骤72，添加三维仿真模型的材料信息：将步骤6中的土工膜设置为超弹性材料，将垫层中随机圆球颗粒设置为刚体。

步骤73，三维仿真模型约束的施加：约束步骤6中土工膜三维仿真模型的四个侧面，在土工膜和垫层随机圆球颗粒之间设置接触关系。

步骤74，向土工膜施加上覆荷载：在土工膜上方施加竖直向下的均布荷载p。

步骤75，土工膜局部变形结果判定：在有限元软件ADINA中，完成三维仿真模型的计算，输出土工膜的变形和应变计算结果；土工膜的局部变形可通过应变大小衡量，若计算得到的最大应变大于设定允许值，表明土工膜已发生过大的局部变形，有出现破坏的隐患。

步骤75中，当计算得到的最大应变大于设定允许值时，通过调整土工膜上覆荷载大小、增加土工膜厚度或改良颗粒粒径级配，来确保土工膜运行安全。

步骤31中，生成第一个d₁粒径球的球心公式计算如下：

x₁＝(L-d₁)t₁+(-L+d₁)/2 (2)

y₁＝(L-d₁)t₂+(-L+d₁)/2 (3)

其中，t₁、t₂均为[0,1]间的均匀随机数。

步骤33中，生成的第二个d₁粒径球的球心坐标[x₂,y₂]的计算公式为：

x₂＝d₁t₃+x-d₁/2 (4)

y₂＝d₁t₄+y-d₁/2 (5)

其中t₃、t₄均为[0,1]间的均匀随机数；x、y分别为所选取网格的中心点的x坐标和y坐标。

步骤33中每两个球心之间的距离计算公式如下：

其中d为两个球心之间的距离；x_i、y_i分别为第i个球心的x坐标和y坐标；x_i+1、y_i+1分别为第i+1个球心的x坐标和y坐标。

本发明具有如下有益效果：

1.本发明根据实际工程垫层料的颗粒级配曲线确定垫层颗粒粒径和颗粒组成，进而生成随机三维圆球颗粒垫层，并添加土工膜和施加上覆荷载，建立了上覆荷载作用下土工膜局部变形的三维仿真模型，能较好地模拟实际工程中随机颗粒垫层上土工膜的受力变形情况。

2.本发明通过将生成的三维仿真模型进行有限元计算，可获得土工膜变形和应变等计算结果，据此可以判断土工膜发生局部变形时是否安全。判别方法为：若土工膜的最大应变值大于允许值，表明土工膜已发生过大的局部变形，会出现破坏；若土工膜的最大应变值小于允许值，表明土工膜安全运行。

3.根据实际工程中预期的土工膜上覆荷载大小及垫层颗粒级配，将本发明建立的随机颗粒垫层上土工膜局部变形的三维仿真模型用于有限元计算，若计算得到的土工膜最大应变值大于允许值，可通过增大土工膜厚度或减小上覆荷载大小或通过经验公式选择合理的垫层颗粒级配，以确保土工膜安全运行。

4.本发明生成的三维不同粒径大小的圆球垫层颗粒的顶部在同一高度，土工膜和垫层表面可以密贴合一致，合理反映了实际工程中碾压密实的垫层表面平整的初始状态。

附图说明

图1显示了本发明生成的随机颗粒垫层的平面效果图。

图2显示了本发明形成的三维仿真模型的效果图。

图3显示了垫层颗粒材料的级配曲线。

图4显示了25mm粒径时正方形区域网格剖分图。

图5显示了22mm粒径时正方形区域网格剖分图。

图6显示了14mm粒径时正方形区域网格剖分图。

具体实施方式

下面结合附图和具体较佳实施方式对本发明作进一步详细的说明。

本发明的描述中，需要理解的是，术语“左侧”、“右侧”、“上部”、“下部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，“第一”、“第二”等并不表示零部件的重要程度，因此不能理解为对本发明的限制。本实施例中采用的具体尺寸只是为了举例说明技术方案，并不限制本发明的保护范围。

一种随机颗粒垫层上土工膜局部变形的三维仿真模型建立方法，包括如下步骤。

步骤1，绘制颗粒垫层区域：在xoy平面内绘制水平放置且边长为L的正方形，并以正方形的中心为原点，以过原点且平行于正方形相邻两边的直线为x轴和y轴，以过原点且垂直于x、y轴的直线为z轴，建立xyz坐标系。本发明中，以L＝200mm为例进行说明。

步骤2，确定垫层区域中不同粒径圆球颗粒的个数，包括如下步骤。

步骤21，确定圆球颗粒粒径和每种粒径所占质量百分比：颗粒垫层中的随机颗粒通过圆球模拟，根据垫层颗粒材料的级配曲线，获得不同圆球颗粒粒径和每种粒径所占的质量百分比。假设垫层颗粒材料的级配曲线，如图3所示，则圆球颗粒粒径及质量百分比如下表1所示：

表1圆球颗粒粒径及质量百分比

步骤22，不同粒径圆球颗粒的个数确定：根据步骤21中确定的每种粒径所占质量百分比可求得不同粒径圆球颗粒的个数；假设不同圆球颗粒粒径从大到小依次为d₁、d₂、…、d_n，每种粒径对应圆球的个数依次为N₁、N₂、…、N_n，则圆球个数的计算公式如下：

其中，N_i为第i种粒径颗粒对应圆球的个数；L为正方形的边长；p_i为第i种粒径颗粒所占的质量百分比；d_i为第i种粒径颗粒的粒径。

如表1所示，粒径d₁、d₂、…、d_n从大到小依次为25mm、22mm、14mm、10mm、5mm和3mm，每种粒径对应圆球的个数N₁、N₂、…、N_n依次为8、9、18、29、91和210。

步骤3，随机生成最大粒径为d₁＝25mm的球，包括如下步骤。

步骤31，随机生成第一个d₁粒径的球：在步骤1绘制的颗粒垫层区域L＝200mm内，先随机生成第一个d₁粒径的球心坐标[x₁,y₁,0]，其球心公式计算如下：

x₁＝(L-d₁)t₁+(-L+d₁)/2 (2)

y₁＝(L-d₁)t₂+(-L+d₁)/2 (3)

其中，t₁、t₂均为[0,1]间的均匀随机数。

注意：需使x₁和y₁均在[(-L+d₁)/2,(L-d₁)/2](即[-87.5,87.5])内，否则，重新生成球心坐标。

步骤32，划分边长为d₁的正方形区域网格：如图4所示，作步骤31生成的d₁＝25mm粒径球的x向和y向四条切线，以d₁＝25mm粒径为边长向四周划分网格，网格的总数为

也即

将d₁＝25mm粒径球的周围每一个网格编号并计算每个网格中心点的坐标。

步骤33，随机生成第二个d₁＝25mm粒径的球：随机选取一个步骤32中的网格编号，在选取的网格内随机生成第二个d₁粒径球的球心坐标，并使第二个d₁粒径球的球心的x坐标或y坐标均在[(-L+d₁)/2,(L-d₁)/2](即[-87.5,87.5])内，且任意两球之间无重叠；否则重复步骤33。

步骤33中，生成的第二个d₁粒径球的球心坐标[x₂,y₂]的计算公式为：

x₂＝d₁t₃+x-d₁/2 (4)

y₂＝d₁t₄+y-d₁/2 (5)

其中t₃、t₄均为[0,1]间的均匀随机数；x、y为所选取网格的中心点的x坐标和y坐标。

步骤33中每两个球心之间的距离计算公式如下：

其中d为两个球心之间的距离；x_i、y_i分别为第i个球心的x坐标和y坐标；x_i+1、y_i+1分别为第i+1个球心的x坐标和y坐标。

步骤34，随机生成剩余d₁＝25mm粒径的球：重复步骤33，直至第N₁个d₁粒径的球全部随机生成。

步骤4，随机生成d₂＝22mm、d₃＝14mm、…、d_n＝3mm粒径的球：将步骤3中粒径d₁的球依次更换为d₂、d₃、…、d_n粒径的球，重复步骤31至步骤34，即可随机生成d₂、d₃、…、d_n粒径的球。

其中，图5显示了d₂＝22mm粒径时的正方形区域网格剖分图；图6显示了d₃＝14mm粒径时的正方形区域网格剖分图。

步骤5，生成随机颗粒垫层：将N₂个d₂粒径的球沿z轴向上平移(d₁-d₂)/2，将N₃个d₃粒径的球沿z轴向上平移(d₁-d₃)/2，依此类推，将N_n个d_n粒径的球沿z轴向上平移(d₁-d_n)/2，使d₁、d₂、…、d_n粒径球的顶部均保持平齐。

本实施例中，具体为：将9个22mm粒径的球沿垂直于正方形区域所在平面向上平移1.5mm，将18个14mm粒径的球向上平移5.5mm，将29个10mm粒径的球向上平移7.5mm，将91个5mm粒径的球向上平移10mm，将210个3mm粒径的球向上平移11mm，所有粒径的球平齐效果如图1和图2所示。

步骤6，生成土工膜：在d₁＝25mm粒径球心所在平面正上方d₁/2＝12.5mm处，生成与边长为L＝200mm的正方形区域大小相同的平面，将生成的平面沿z轴向上拖放h厚度，其中，h为土工膜厚度，优选为2mm，即得土工膜；随机颗粒垫层和土工膜共同构成三维仿真模型。

还包括步骤7，土工膜局部变形计算与破坏判定，具体包括如下步骤：

步骤71，将三维仿真模型导入有限元计算软件：将步骤6形成的三维仿真模型导入有限元计算软件ADINA中，并对土工膜进行网格剖分，优选x方向剖分份数为50、y方向份数为50和z方向份数为4。

步骤72，添加三维仿真模型的材料信息：将步骤6中的土工膜设置为超弹性材料，优选弹性模量E为300MPa，泊松比υ为0.46，将垫层中随机圆球颗粒设置为刚体。

步骤73，三维仿真模型约束的施加：约束步骤6中土工膜三维仿真模型的四个侧面，在土工膜和垫层随机圆球颗粒之间设置接触关系。

步骤74，向土工膜施加上覆荷载：在土工膜上方施加竖直向下的均布荷载p，优选为100kPa。

步骤75，土工膜局部变形结果判定：在有限元软件ADINA中，完成三维仿真模型的计算，输出土工膜的变形和应变计算结果；土工膜的局部变形可通过应变大小衡量，若计算得到的最大应变大于设定允许值，优选为12％，表明土工膜已发生过大的局部变形，有出现破坏的隐患。

步骤75中，当计算得到的最大应变大于设定允许值时，通过调整土工膜上覆荷载大小、增加土工膜厚度或通过经验公式选择合理的垫层颗粒级配，以确保土工膜运行安全。

其中，垫层颗粒级配合理选择的经验公式如下：

式中，d_i为第i种粒径颗粒的粒径；p_i为第i种粒径颗粒所占的质量百分比；d_max为垫层颗粒的最大粒径。

以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。

Claims (5)

1.一种随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，其特征在于：包括如下步骤：

步骤1，绘制颗粒垫层区域：在xoy平面内绘制水平放置且边长为L的正方形，并以正方形的中心为原点，以过原点且平行于正方形相邻两边的直线为x轴和y轴，以过原点且垂直于x、y轴的直线为z轴，建立xyz坐标系；

步骤2，确定垫层区域中不同粒径圆球颗粒的个数，包括如下步骤：

步骤21，确定圆球颗粒粒径和每种粒径所占质量百分比：颗粒垫层中的随机颗粒通过圆球模拟，根据垫层颗粒材料的级配曲线，获得不同圆球颗粒粒径和每种粒径所占的质量百分比；

步骤22，确定不同粒径圆球颗粒的个数：根据步骤21中确定的每种粒径所占质量百分比可求得不同粒径圆球颗粒的个数；假设不同圆球颗粒粒径从大到小依次为d₁、d₂、…、d_n，每种粒径对应圆球的个数依次为N₁、N₂、…、N_n，则圆球个数的计算公式如下：

其中，N_i为第i种粒径颗粒对应圆球的个数；L为正方形的边长；p_i为第i种粒径颗粒所占的质量百分比；d_i为第i种粒径颗粒的粒径；

步骤3，随机生成最大粒径为d₁的球，包括如下步骤：

步骤31，随机生成第一个d₁粒径的球：在步骤1绘制的颗粒垫层区域内，先随机生成第一个d₁粒径的球心坐标[x₁,y₁,0]，并使x₁和y₁均在[(-L+d₁)/2,(L-d₁)/2]内，否则，重新生成球心坐标；

步骤32，划分边长为d₁的正方形区域网格：作步骤31生成的d₁粒径球的x向和y向四条切线，以d₁粒径为边长向四周划分网格，网格的总数为

将d₁粒径球的周围每一个网格编号并计算每个网格中心点的坐标；

步骤33，随机生成第二个d₁粒径的球：随机选取一个步骤32中的网格编号，在选取的网格内随机生成第二个d₁粒径球的球心坐标，并使第二个d₁粒径球的球心的x坐标或y坐标均在[(-L+d₁)/2,(L-d₁)/2]内，且任意两球之间无重叠；否则重复步骤33；

步骤34，随机生成剩余d₁粒径的球：重复步骤33，直至第N₁个d₁粒径的球全部随机生成；

步骤4，随机生成d₂、d₃、…、d_n粒径的球：将步骤3中粒径d₁的球依次更换为d₂、d₃、…、d_n粒径的球，重复步骤31至步骤34，即可随机生成d₂、d₃、…、d_n粒径的球；

步骤5，生成随机颗粒垫层：将N₂个d₂粒径的球沿z轴向上平移(d₁-d₂)/2，将N₃个d₃粒径的球沿z轴向上平移(d₁-d₃)/2，依此类推，将N_n个d_n粒径的球沿z轴向上平移(d₁-d_n)/2，使d₁、d₂、…、d_n粒径球的顶部均保持平齐；

步骤6，生成土工膜：在d₁粒径球心所在平面正上方d₁/2处，生成与边长为L的正方形区域大小相同的平面，将生成的平面沿z轴向上拖放h厚度，其中，h为土工膜厚度，即得土工膜；随机颗粒垫层和土工膜共同构成三维仿真模型；

步骤7，土工膜局部变形计算与破坏判定，具体包括如下步骤：

步骤71，将三维仿真模型导入有限元计算软件：将步骤6形成的三维仿真模型导入有限元计算软件ADINA中，并对土工膜进行网格剖分；

步骤72，添加三维仿真模型的材料信息：将步骤6中的土工膜设置为超弹性材料，将垫层中随机圆球颗粒设置为刚体；

步骤73，三维仿真模型约束的施加：约束步骤6中土工膜三维仿真模型的四个侧面，在土工膜和垫层随机圆球颗粒之间设置接触关系；

步骤74，向土工膜施加上覆荷载：在土工膜上方施加竖直向下的均布荷载p；

步骤75，土工膜局部变形结果判定：在有限元软件ADINA中，完成三维仿真模型的计算，输出土工膜的变形和应变计算结果；土工膜的局部变形可通过应变大小衡量，若计算得到的最大应变大于设定允许值，表明土工膜已发生过大的局部变形，有出现破坏的隐患。

2.根据权利要求1所述的随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，其特征在于：步骤75中，当计算得到的土工膜最大应变大于设定允许值时，通过调整土工膜上覆荷载大小、增加土工膜厚度或改良颗粒粒径级配，来确保土工膜运行安全。

3.根据权利要求1所述的随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，其特征在于：步骤31中，生成第一个d₁粒径球的球心公式计算如下：

x₁＝(L-d₁)t₁+(-L+d₁)/2 (2)

y₁＝(L-d₁)t₂+(-L+d₁)/2 (3)

其中，t₁、t₂均为[0,1]间的均匀随机数。

4.根据权利要求1或3所述的随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，其特征在于：步骤33中，生成的第二个d₁粒径球的球心坐标[x₂,y₂]的计算公式为：

x₂＝d₁t₃+x-d₁/2 (4)

y₂＝d₁t₄+y-d₁/2 (5)

其中t₃、t₄均为[0,1]间的均匀随机数；x、y分别为所选取网格的中心点的x坐标和y坐标。

5.根据权利要求4所述的随机颗粒垫层上土工膜局部变形三维仿真模型建立方法，其特征在于：步骤33中每两个球心之间的距离计算公式如下：

其中d为两个球心之间的距离；x_i、y_i分别为第i个球心的x坐标和y坐标；x_i+1、y_i+1分别为第i+1个球心的x坐标和y坐标。

Images