强度折减条件下膜袋坝应力和变形响应特征离散元分析

刘磊， 尚静石， 李维朝， 王先鹏， 高连瑞

(1.中国地质大学(北京) 工程技术学院,北京 100083; 2.中国水利水电科学研究院，流域水循环模拟与调控国家重点实验室,北京 100038; 3.大唐双鸭山热电有限公司,黑龙江 双鸭山 155100)

煤炭是我国能源资源的重要组成部分。每年我国燃煤电厂所产出的粉煤灰量约占世界总产量的一半，我国粉煤灰综合利用率仅为70%[1]，过剩粉煤灰需要建设贮灰场或在原贮灰场上扩容贮放。目前贮灰场的扩容工程大多是在原贮灰场的基础上修建子坝。

膜袋坝可以像粉煤灰筑坝一样，采取水力冲填粉煤灰的方式填筑，施工简便快捷，可有效降低坝体土石料的用量，同时增加贮灰场的贮灰量[2]。膜袋坝还可以有效利用土工膜袋的力学特性[3]，增大膜袋坝的抗剪强度与稳定性，利用膜袋坝来修建子坝的工程数量不断增多。

关于土工膜袋的力学性质与工程应用，樊科伟等[4]研究了袋装石土工膜袋剪切力学特性，刘斯宏等[5]、李巧燕等[6]、潘洋等[7]研究了土工膜袋层间摩擦特性，刘斯宏等[8]研究了土工膜袋在加固地基中的应用，娄汉宇等[9]研究了土工膜袋与混凝土砂浆填充料相结合在堤防加固中的应用。这些研究也证明了土工膜袋具有良好的工程特性。

虽然膜袋坝具有较高稳定性，且使用量不断增加，但当前关于膜袋坝的变形特征、作用机理和计算方法等均尚不明确。当前计算膜袋坝稳定性及变形特征的计算方法有水平条分法、基于界面单元的有限元的强度折减法等。这些方法虽然可以很好地考虑土工膜袋的界面滑动特点，但都假设土工膜袋的界面是水平的，而实际上土工膜袋界面是柔性的，施工后的土工膜袋会受到不均匀固结、施工错缝等影响，导致土工膜袋界面并不是水平的而是会出现镶嵌现象。

鉴于此，本文利用离散元方法构建了实体土工膜袋，考虑土工膜袋界面非水平的镶嵌作用，通过强度折减的方法分析膜袋坝稳定性降低过程中的变形特征与受力特征，并对土工膜袋的加固机理进行分析。

1 工程概况

1.1 贮灰场概况

某贮灰场为山谷水灰场。山谷为南北走向，谷口朝北，贮灰场东西两侧为灌木覆盖的山坡。坝址处最低原地面标高为146.5 m、两侧最低凹处标高约205.0 m的山坡[10]。

1.2 膜袋坝设计概况

贮灰场坝体由主坝、子坝、灰场内隔坝组成，其中贮灰场内由北向南的两座隔坝把灰场划分出3个灰池，如图1所示。1号灰池面积12×104m2，2号灰池面积9×104m2，3号灰池面积8.5×104m2[10]。

图1 贮灰场俯视图

本次以2号灰池西侧A断面的两级水力充填膜袋坝为研究对象，分别如图2和图3所示。这两个子坝均是直接在贮灰场所贮粉煤灰上筑建。其中一级子坝宽约17.5 m、高约3.0 m；二级子坝宽约7.5 m、高约1.5 m；所有膜袋截面宽约2.5 m、高约0.5 m，材料为长丝机织土工膜，且材料质量满足GB/T 17640—1998的要求。

图2 膜袋坝断面图

图3 土工膜袋坝

2 膜袋坝离散元建模

2.1 土工膜袋及膜袋内粉煤灰建模

2.1.1 选取本构模型

离散元软件PFC3D中可以对颗粒之间的接触赋予不同的接触本构模型，从而模拟不同材料的特性。在同一模型中，可能包含多种材料类型，或者相同材料之间的力学性质差异较大。通过对模型中不同性质材料的分组胶结，赋予其不同的接触本构模型，从而解决单个模型中材料性质多元化的问题。PFC3D中包含10种本构模型，其中包括线性模型、接触黏结模型、线性平行黏结模型、Hertz接触模型、滞回接触模型、光滑节理接触模型、平直节理接触模型、滚动阻力线性模型、Burgess模型[11]。

土工膜袋及膜袋内粉煤灰分别选用接触黏结模型(Linear Contact Bond Model)和线性模型(Linear Model)。其本构特性分别如图4—5所示。图中：Fd为阻尼力；Fl为线性力；Fc为接触力；Mc为接触力矩；Dc为接触间距；kn为法向刚度；ks为切向刚度；gs为表面间隙；μ为摩擦系数；TF为抗拉强度；SF为抗剪强度。

图4 接触黏结模型接触特性

图5 线性模型接触特性

2.1.2 土工膜袋的离散元建模

土工膜袋建模过程分为两步：第一步，生成对应真实土工膜袋形状的离散元模型；第二步，对该模型进行胶结赋值。

完成第一步主要有两条途径：①PFC3D内置Fish语言结合“ball create”命令，首先把表达土工膜袋模型颗粒位置的数学函数写进Fish语言中，其次调用“ball create”命令用于生成颗粒。这种方法适用于外观形状比较规则的模型。②首先，通过AutoCAD、犀牛等辅助软件，生成PFC3D中可读文件，如“DXF”、“STL”等类型文件；其次，利用“geometry”命令将生成的文件导入PFC3D中，在PFC3D中生成表示模型形状的“geometry”模块；最后，把颗粒填进“geometry”区域内。文中采用第二种方法建立土工膜袋离散元模型，步骤如下。

1)使用三维AutoCAD软件生成两个闭合“STL”文件，用上述方法导入PFC3D中，得到一内一外两个“geometry”，如图6所示。然后，调用“wall import”命令在PFC3D中建立与“geometry”相同形状的两墙体。图6中绿色部分称之为内墙，灰色部分称之为外墙。2个三维墙体之间的距离处处相同，等于土工膜袋颗粒的直径(1.8 mm)。

2)墙体建立完成后，在墙体之间填充代表土工膜袋的颗粒。由于土工膜袋与其内部所填充的粉煤灰属于不同的材料，需要分别对其赋予不同的接触类型。所以，生成土工膜袋颗粒后，把此部分颗粒进行分组。分组后，对土工膜袋颗粒间的接触施加胶结，接触本构模型为接触黏结模型(Linear Contact Bond Model)。模型中的墙体暂时保留，以便向土工膜袋内填充粉煤灰。

图6 单个土工模袋建模形状

2.1.3 土工膜袋内粉煤灰的离散元建模

使用“ball distribute”命令，向内墙中(图6中绿色部分)填充直径为1.8 mm的粉煤灰离散元颗粒。粉煤灰颗粒接触本构模型为线性模型(Linear Model)。生成颗粒并运行一段时间后，删除与内墙重叠较多的粉煤灰颗粒，这么做的目的是减小模型应力集中现象，防止初始平衡过程中外部土工膜袋颗粒被内部具有较大速度的粉煤灰颗粒击破。

粉煤灰颗粒填充完成后，此时模型颗粒已经全部生成，但仍未删除墙体。需要注意的是，模型中墙体的存在会导致部分颗粒与墙体之间产生重叠，使这部分颗粒存在较大的应力。如果删除墙体后，不采取相应的措施使模型中的应力缓慢释放，就会导致模型颗粒之间的胶结会因不平衡力的突然增大而破坏。因此，删除墙体后，每运行一段时间就清零模型中颗粒的速度，使得应力缓慢释放。

离散元软件中通过1∶1的真实比例建立模型会产生巨量的颗粒，从而影响计算速度。鉴于此，通过缩尺的方式来建立模型，并利用离心机原理[12]在保证计算速度的基础上来实现对膜袋性能的真实分析。模型缩尺比例为:模型/原型=1/100，具体说明由“参数选取”一节给出。图7是根据上述建模过程所生成的土工膜袋，该膜袋长20 mm、宽25 mm、高5 mm。所列尺寸为模型中的数据。

图7 土工膜袋

2.2 土工膜袋界面

土工膜袋模型是由单独的颗粒组合胶结而成，粒子的相互作用通过内力和力矩实现。颗粒离散元法能对大变形物体进行准确的分析，且不受限于分析对象的变形量，这对于研究非连续介质力学有较大的优势。因此能在离散元软件中建立界面为柔性的土工膜袋模型。为直观表达土工膜袋受力的变形特性，模拟了对单个土工膜袋的重力加载试验。

单个土工膜袋重力加载数值模拟试验步骤如下：①建立数值试验模型，模型分为上下两部分，其上为土工膜袋，膜袋下方为位置固定的刚性墙体，墙体宽度为土工膜袋的一半，如图8所示。②对数值试验模型施加重力场，让土工膜袋自由下落。

图8 膜袋界面变形

从图8中可以明显看出，土工膜袋在受到弯矩时发生明显变形，左右两部分弯折成一定的角度。这使得在进行膜袋坝稳定性计算时，不必假设土工膜袋的界面都是水平的，在分析膜袋坝受力变形时能够考虑土工膜袋界面的变形、镶嵌等现象所带来的影响，使数值模拟结果更加接近实际情况。

2.3 膜袋坝整体离散元建模

2.3.1 结合施工顺序建模

膜袋坝如图9所示，一级子坝与二级子坝在图中用绿色部分表示。子坝均由单个土工膜袋堆砌而成。图9中其余颜色均表示粉煤灰，不同颜色表示不同施工顺序的粉煤灰。B1、B2、…、B10为参数折减过程中的位移监测点位置。

图9 膜袋坝

模型各部分的生成及固结时间顺序参考了施工总说明[13]。实际施工过程中，采用分层填筑的方法，即对每一层土工膜袋进行固结，并在具有代表性的位置处进行取样检测，满足要求后再堆砌新一层土工膜袋。且上下层膜袋之间错缝施工。鉴于此，膜袋坝模型建立顺序如下：首先，生成粉煤灰地基颗粒(图10中蓝色部分)并进行固结；其次，只导入一层土工膜袋并固结，最后循环该步骤直到一级坝体建造完成。二级坝体以及贮灰场内粉煤灰颗粒的建立方法和顺序与上述说明类似，示意图如图10所示。

图10 模型建立顺序

2.3.2 贮灰场内粉煤灰的离散元建模

粉煤灰属无黏性土料，可以忽略颗粒之间的黏结作用，选用线性模型(Linear Model)建模，其本构特性如图5所示。粉煤灰颗粒参数信息由“3 参数选取”一节给出。

2.3.3 高效建立模型的方法

因模型中包含许多独立的土工膜袋个体，这导致在单个计算范围内逐一生成各部分模型并固结的过程十分繁琐，且增多计算时间。文中采用生成“Brick”填充法，使单个膜袋与各层粉煤灰单独生成，然后保存为“Brick”文件，把这些文件看作组成模型的“块”。当组合模型用到这些“块”时就逐一导入。这种方法的优势在于能够短时间内生成大量颗粒，简化建模过程，节省建模时间。

3 参数选取

为更清楚地描述参数的选取，把模型分为4部分：地基与贮灰场内的粉煤灰、土工膜袋、墙体、边界条件，并依次介绍其参数及参数选取依据。

3.1 地基与贮灰场内的粉煤灰

3.1.1 地基与贮灰场内的粉煤灰参数的选取

该粉煤灰颗粒间接触本构模型为线性模型，颗粒直径为3 mm，法向刚度为1×106N/m,剪切刚度为1×106N/m，摩擦系数为0.5，局部阻尼常数为0.05，比重为2.12，颗粒数为23 889。

3.1.2 地基与贮灰场内的粉煤灰参数选取依据

参考该工程勘察报告[10]，粉煤灰比重取2.12，内摩擦角为28.4°，摩擦系数取其近似值，为0.5。通过查阅文献，粉煤灰的压缩模量为2.15～16.36 MPa[19-21]，取接触的法向刚度kn为1×106N/m，kn/ks=1，ks为剪切刚度。

3.2 土工膜袋

3.2.1 土工膜袋参数的选取

膜袋颗粒接触本构模型(接触粘结模型)、界面接触本构模型为线性模型，膜袋颗粒直径为1.8 mm，抗拉强度为20 N，法向刚度为1×103N/m，剪切刚度为1×103N/m，颗粒间摩擦系数为0.5，界面间摩擦系数为0.5，局部阻尼常数为0.05，比重为2.12，膜袋坝颗粒数为15 675，膜袋坝颗粒之间的接触总数为51 522。

3.2.2 土工膜袋参数选取的依据

参考膜袋坝施工说明[13]，土工膜袋质量满足GB/T 17640—1998的要求，因此，土工膜袋的抗拉强度参数同样满足此标准要求，模型中土工膜袋的抗拉强度(cb_tenf)设为20 N，值得注意的是此参数单位为N，这是因在PFC3D接触本构中是通过定义单个接触所能承受的最大拉力作为抗拉强度的。把此参数换算成GB/T 17640—1998中的单位后，其数值满足该标准中关于土工膜袋的抗拉强度的要求。

在离散元中模拟土工膜袋时，最重要的一点是突出土工膜袋界面的柔性，使土工膜袋界面的变形特征与实际相吻合。通过设置接触黏结本构模型中的法向刚度kn与剪切刚度ks来调节界面的变形特征。鉴于此，设置了单个土工膜袋重力加载数值模拟试验来调节参数值(图8)，试验具体步骤见第2.2节。通过试验验证，取kn=ks=1×103N/m时，土工膜袋界面的变形特性满足要求。

关于土工膜袋界面摩擦系数，刘斯宏等[5]分析了土工膜袋内材料粒径和土工膜袋的排列方式对土工膜袋界面摩擦系数的影响，指出在不同试验条件下，土工膜袋界面摩擦系数为0.442～0.619。李巧燕等[6]进行了土工膜袋的斜坡滑移试验，得到土工膜袋在不同固结天数下的摩擦系数在0.46～0.60之间。潘洋等[7]通过试验得到相邻土工膜袋之间的摩擦系数为0.564。参考上述试验结果并结合本文模拟的特点，土工膜袋界面摩擦系数取0.5。

3.3 墙体

3.3.1 墙体参数的选取

墙体与颗粒接触本构模型为线性模型，法向刚度为1×108N/m，剪切刚度为1×108N/m，摩擦系数为0.5。

3.3.2 墙体参数选取的依据

为了防止颗粒与墙体之间存在较大的重叠，墙体的法向刚度与剪切刚度一般设置为颗粒刚度的100倍。因此文中kn=ks=1×108N/m。

3.4 边界条件

3.4.1 边界条件的确定

模型X方向长度为0.46 m，模型Y方向宽度为0.02 m，模型坝体内粉煤灰高度为0.09 m，重力场为100g。

3.4.2 边界条件确定的依据

离散元中通过1∶1的真实比例建立模型，此模型中颗粒的数量将超出模拟计算能力范围，所以，一般采用缩尺的方法建立模型。但此方法导致模型中的应力水平同原型不一致，因此，许多学者[14-18]选用了离心机原理建模，即通过对缩尺后的模型施加高水平重力场，使模型与原型应力水平相似。鉴于此，本文建模采用了离心机原理，重力场取100g，模型尺寸/工程实际尺寸=1/100，模型尺寸如图9所示。

4 强度折减响应分析

4.1 强度折减法的应用

强度折减法是通过对研究对象强度参数的折减，直至其达到失稳破坏(失稳的判别方法主要包括位移突变法等)，从而得到模型的安全系数。此次模拟与常规计算边坡的安全系数不同，本文研究重点是通过对膜袋坝和粉煤灰强度参数的折减，研究膜袋坝在稳定性降低条件下的位移特征，分析稳定性降低过程中土工膜袋的力学响应，揭示土工膜袋加强坝体稳定性的机理。

4.2 强度折减位移分析

折减系数取值范围为1.2～3.0，增量为0.2。数值计算过程中以不平衡力小于1×10-5(这是一比值，无量纲)为终止条件，通常满足此条件的数值模型中的颗粒速度已经处于非常低的水平，可认为数值模型达到了稳定状态。

图11为膜袋坝在折减过程中水平位移变化云图，颜色区间为蓝色～红色，对应原型水平位移为-12～12 cm。图11中水平轴X的右向为正。由图11可看出，膜袋坝的水平位移随着折减系数的增大而增大，膜袋坝最大水平位移出现在坝体边坡表层土工膜袋中上部(表层土工膜袋位置如图11中的标示)。通过观测数值模型中的蓝色区域的发育过程(蓝色表示颗粒的水平位移相对较大)，发现滑动区域的发育过程呈现出由坝体边坡表层土工膜袋逐级向内延伸的特点。由于土工膜袋界面之间的镶嵌与摩擦作用，存在较大位移的土工膜袋影响了与其相接触膜袋的运动。这与连续材料边坡的失稳特征有较大的差异。

图11 膜袋坝水平位移随折减系数变化云图

为更直观地展示土工膜袋的水平位移随折减系数的变化，设置了B1、B2、…、B10等位移监测点，监测点位置如图9所示，图12中表示出各监测点水平位移随折减系数的变化情况。在图12中可以明显看出，在折减系数为2.0时曲线出现转折。这说明膜袋坝在折减系数由1.2增加至2.0时水平位移变化较大，而当折减系数从2.0增加到3.0时，土工膜袋的水平位移虽然也在增大，但增量却明显减小。均匀边坡的变形会随着强度折减量的增大而增加，当该增量达到临界值后，边坡变形量会出现骤增。这表明膜袋坝的变形特征与均匀边坡的截然不同。出现此现象的原因主要为土工膜袋的作用：土工膜袋不仅增强了单个膜袋内粉煤灰的抗剪强度，且膜袋与膜袋界面的镶嵌、摩擦作用，使得膜袋坝更稳定。为验证土工膜袋的作用，揭示土工膜袋特性在阻止坝体失稳过程中的内在机理，提取了数值模型中土工膜袋的拉力信息，用土工膜袋在抵抗失稳过程中的力学响应解释膜袋坝为何没有发生大的变形破坏。

图12 监测点水平位移随折减系数的变化

4.3 土工膜袋加强坝体稳定性的机理分析

在膜袋坝稳定性降低过程中，监测土工膜袋的力学响应是分析土工膜袋的受力特征及其加筋作用的基础。在离散元软件中利用Fish语言可以提取精确到每一个颗粒的受力信息，可实现对土工膜袋的细观力学分析。在数值模型中提取了第一层土工膜袋总拉力、第二层土工膜袋总拉力、膜袋坝(包括一级子坝和二级子坝)总拉力信息。膜袋编号如图13所示。

图13 膜袋编号

图14—16分别为第一层土工膜袋总拉力随折减系数的变化、第二层土工膜袋总拉力随折减系数的变化、膜袋坝总拉力随折减系数的变化情况。由图14—16可以明显看出，这3个图中拉力随折减系数的变化趋势相同，即土工膜袋在坝体失稳过程中拉力不断增加。这证明了土工膜袋在坝体失稳过程中产生了层间变形镶嵌、摩擦作用，从而使得膜袋的拉力增加，同时这种层间作用使膜袋坝整体更加稳定，具有更高的安全系数。

图14 第一层膜袋总拉力随折减系数的变化情况

图15 第二层膜袋总拉力随折减系数的变化情况

图16 膜袋坝总拉力随折减系数的变化情况

为展示单个膜袋拉力的变化趋势，提取了第二层6个膜袋的数据，图17为不同折减系数下第二层编号从1到6单个膜袋的总拉力，图18示出6个单膜袋总拉力随折减系数的变化曲线。从图17中可以看出，同一水平层次的膜袋，位置越靠右侧膜袋拉力越大。结合图17与图18可以看出，随着折减系数的增大，各个膜袋的拉力呈现增加的趋势。

图17 不同折减系数下第二层单个膜袋的总拉力

图18 第二层单个膜袋总拉力随折减系数的变化情况

上述表明，无论是单个膜袋还是整个膜袋坝，其拉力在抵抗坝体失稳过程中由于层间作用的产生，拉力都是不断增大的。与单一材料的碾压粉煤灰坝相比，膜袋拉力的增加与坝体的加筋作用使得膜袋坝成为一个更加稳定的整体，证明了土工膜袋特性在阻止坝体失稳过程中的作用。

5 结论

1)以实际膜袋筑坝工程为例，以PFC3D软件为基础，详细介绍了土工膜袋的离散元建模方法，利用“geometry”命令建立的内外两层墙体使建立土工膜袋数值模型的过程更加方便。同时，对土工膜袋的参数选取进行了详细介绍，并进行了土工膜袋界面数值试验，依据试验数值选取来保证土工膜袋界面的柔性特征。

2)结合实际膜袋坝的分阶段建立过程，提出了以单个土工膜袋为基础单元(称之为“块(Brick))”，当组合模型时需要用到这些“块”时就逐一导入的方法。这种方法的优势在于建立膜袋坝模型时，尺寸相同的膜袋只需生成一次，并把其作为模板使用，一次可导入多个膜袋，简化了建模过程，节省了建模时间。

3)使用了强度折减法，通过对膜袋坝强度的折减，不断降低膜袋坝的稳定性，展示膜袋坝在稳定性降低过程中的位移特征，分析了受柔性膜袋界面镶嵌、摩擦作用影响的膜袋坝位移特征与均匀边坡位移特征的不同。

4)利用离散元方法在微观力学分析方面的优势，提取了土工膜袋拉力。展示了膜袋坝在稳定性降低过程中土工膜袋的力学响应，突出了膜袋在加固坝体应用中的力学特性与加筋作用，结合其力学响应初步揭示了土工膜袋加强坝体稳定性的机理。

5)本文结合实际施工过程，利用离散元方法建立膜袋坝，对膜袋坝稳定性降低过程中的位移特征与受力进行了分析，这些工作为今后分析膜袋坝的工程特性等提供了有益的帮助。