1、DOI：10.19645/j.issn2095-0144.2022.10.004收稿日期：2022-09-08基金项目：国家自然科学基金（41771066）作者简介：李明（1998-），男，江苏连云港人，硕士研究生，研究方向：岩土工程，E-mail：。通信作者：刘恩龙（1976-），男，河南商丘人，教授，博士，主要从事土的本构关系与数值模拟研究，E-mail：。基于连续-离散模拟方法的冰川失稳过程与机理分析李明1，刘恩龙1,2，苏雨1，康建2，冯吉利3（1.四川大学 水利水电学院，四川 成都 610065；2.中国科学院 西北生态资源环境研究院，甘肃 兰州 730000；3.中国矿业大学（北京
2、）力学与建筑工程学院，北京 102206）摘要：在气候变暖的背景下，冰川运动逐渐活跃，发生冰崩灾害的可能性和危险性随之增加，研究冰川失稳过程与机理具有重要的现实意义。基于有关文献中的人工制备多晶冰室内静力三轴试验结果，采用FDEM方法建立多晶冰试样的数值模型标定参数、分析破坏机理，进行了不同条件下西藏阿里地区阿汝冰川失稳过程的二维数值模拟。通过FDEM方法计算冰体从连续到离散的变形特征，描述出阿汝冰川的从稳态开始、裂缝的产生及发展、局部滑移、崩解和碰撞反弹等破坏过程。同时，发现冰川和基岩的连接界面摩擦减小比冰川断裂增加对冰川失稳的影响更大。FDEM方法可以应用于冰川失稳破坏研究，在研究非连续问
3、题领域也有相当大的潜力。关键词：连续-离散（FDEM）；数值模拟；冰川失稳；阿汝冰川中图分类号：P343.6文献标志码：A文章编号：2095-0144（2022）10-0013-07第 58 卷 第 10 期2022 年 10 月GANSU WATER RESOURCES AND HYDROPOWER TECHNOLOGY甘 肃 水 利 水 电 技 术Vol.58，No.10Oct.，2022我国是世界上中、低纬度冰川最发育的国家，冰川资源储量4.31034.7103km3，位列世界第四1。青藏高原构造隆升强烈，气候敏感多变，孕灾环境较好，是地震、冰崩、冰湖溃决等灾害的多发区2，同时又是我国边
4、疆高寒地区和战略高地。因此，研究冰崩等高亚洲地区的灾害事件，对于“一带一路”建设和国民经济发展都具有重大意义。导致冰崩的原因有很多，如冰床坡度的陡增、陡坎的出现、冰中融水的分布等。EMMER等3研究了地震、强降雨等动态原因和深层冰的融化、静水压力等长期原因的影响及冰碛湖溃决的破坏机制。SCHNEIDER等4揭示了导致冰崩不同路径和流动特征的驱动因素。SCHAUB等5应用RAMMS和IBER软件进行了由冰崩引发的冲击波过程链的耦合数值模拟，量化了冰崩和冲击波之间的动量转化。NSSER等6对1856-2020年鲁帕尔河谷的5座冰川进行详细分析，强调了冰崩、浪涌的不稳定性和特定地点的地形特征对单个冰
5、川变化的重要性。童立强等7从引发冰崩的方式、冰崩运动的过程和影响因素等角度分析，将冰崩灾害类型分为冰崩直接灾害、冰崩-冰湖溃决灾害和冰崩-堵溃链式灾害。李晨毓等8以长江源各拉丹东地区的冰川为研究对象，发现该地区冰川变化趋势因分区与规模而异，部分区域有明显的冰川降级、分裂现象。近年来，使用连续-离散耦合分析方法（FDEM）研究科学和工程问题的学者越来越多，FDEM的应用范围迅速扩大，包括建筑物爆破拆除、混凝土重力坝的破坏模拟、隧道和边坡的失稳破坏、矿物加工模拟、岩体冒落破坏分析等9-13。在冰川监测、数值模拟方面，前人已经做了诸多研究工作，但对于冰川失稳过程和机理的研究还不够深入。在连续-离散耦
6、合分析方法中引入节理单元能有效地捕捉裂纹扩展路径，模拟渐进破坏过程。基于有关文献中的多晶冰静力三轴试验14结果，选取西藏阿里地区阿鲁错湖区发生冰崩的冰川体为研究对象，采用二维FDEM方法对静力作用下冰川的稳定性进行分析，探讨用连续-离散耦合分析方法13研究冰川堆积坡体起动机理和滑动过程的可行性。1模型试验数值模拟1.1FDEM原理FDEM的基本思路是初始阶段用有限元法描述连续域，在荷载作用下，依据一定断裂准则，连续域内如果产生断裂则不连续出现，从而形成离散体。实施方法为通过有限元离散化，在每对实体单元之间插入无厚度节理单元，节理单元失效之前，采用有限单元法模拟材料的连续变形。当拉伸或者剪切强度
7、达到材料峰值强度时，材料开始屈服并发生张开或剪切滑移。当材料断裂能释放后，节理单元破坏、裂纹形成，同时相应的裂纹单元从连续计算模型中被移除。接触检测算法一旦检测到单元接触对，就会执行接触相互作用算法来计算离散单元之间的作用力。需要注意的是，这些节理单元最开始并不是真实存在的节理，只有在断裂之后才成为真节理。由于微观结构缺陷和应力集中现象的存在，应力和应变场需要修改，基于应变的黏聚裂纹模型15已经在FDEM代码中得到实现。1.2数值试验方案采用二维平面模型进行双轴数值压缩试验来标定模拟参数，研究对象为-15 条件下的人造多晶冰样14，直径61.8 mm，长度125.0 mm，平均密度850 kg
8、/m3。利用LS-PrePost软件划分网格，如图1所示。计算域由5 040个三角形单元组成，节点数为2 620，其中最小单元尺寸为2.0 mm。FDEM法采用显式积分格式进行计算，为保证计算的收敛性和稳定性，计算时间增量步为510-5ms，总时间步数为2 500 000。静力模拟试验下加载板固定，上加载板加载速率为0.1 m/s。刚性加载板与数值试样之间的摩擦系数取为0.1。围压分别为0.5MPa、1.0MPa、2.0MPa、4.0MPa和6.0 MPa。计算过程中采集了试样的应力时程以及试样的变形破坏特征等数据，作为数值分析的内容。1.3计算参数单元的弹性模量和密度均由静力三轴试验14结果
9、得出，泊松比根据前人的经验选取。在FDEM数值模拟中，节理单元的参数包括抗拉强度ft、黏聚力c、内摩擦角、I型断裂能Gf、型断裂能Gf、断裂罚因子。黏滞阻尼=2h E（h为单元尺寸；E为弹性模量；为密度）。参考张红彪16对黄河冰力学性能研究中抗拉强度和应变速率的关系:t=AB（1）式中：t为冰的抗拉强度（MPa）；为应变速率；A和B为与温度有关的参数,通过试验模拟，当温度为-15 时，A=0.298 8，B=-0.070 0。故对多晶冰试样的抗拉强度取0.549 MPa。断裂韧度与温度的关系为：KIC=-29.792ln（-T）+64.945 7（2）式中：KIC为冰的断裂韧度（kPam0.5
10、）；T为温度（）。当温度为-15 时，KIC=145.62 kPam0.5，根据Gf=KIC2/E可得，Gf=25 N/m。对多晶冰试样的断裂能Gf取25 J/m2，Gf取50 J/m2。接触罚因子取10E，切向罚因子取E，断裂罚因子取5E。将3作为参数引入数值模拟过程中，黏聚力和内摩擦系数随3的变化而变化。加载块为刚性块，密度设置为7 850 kg/m3，弹性模量为200 GPa，泊松比为0.29。同时加载块不设置节理单元。最终确定的FDEM模拟参数见表1。1.4计算结果及分析绘出-15 时不同围压条件下静力三轴试验的数值模拟应力-应变曲线、体变-应变曲线，如图2、图3所示。由图2可知，模拟
11、试验得到的应力-应变曲线均表现为应变软化型。当轴向应变较小时，偏差应力随轴向应变增加而增大，达到峰值后，偏差应力开图1计算模型网格2022年第10期甘肃水利水电技术第58卷14始减小，表现出明显的峰后软化现象。整个过程表现出脆性材料的特点。随着侧向应力的增大，冰试样的最大偏差应力逐渐增大，在侧向应力达到4.0 MPa之后，最大偏差应力反而降低。残余强度随侧向应力的增大而不断增大。图3给出了体变-应变模拟结果。围压在0.5MPa、1.0 MPa 和 2.0 MPa 下，试样先体缩后体胀。高围压下，试样一直处于体积压缩状态。整体数值模拟结果与室内试验结果保持一致，为后续利用数值试验的参数进行冰川失
12、稳的数值模拟提供了支撑。2冰川失稳过程模拟与机理分析2.1研究区概况2016年7月17日至9月21日，位于西藏阿里地区，靠近阿鲁错湖的两座冰川坍塌。阿汝冰川的下部首先坍塌，高程在5 2506 150 m，坍塌的冰川流动了7.0 km，覆盖了89 km2的范围，并到达了阿鲁错湖18。仅两个月后，东南方向2.4 km处的第2座冰川又发生坍塌，之后碎裂成一条质量流，延伸5.0 km并覆盖了约7 km2的范围。结合GILBERT等19的实地考察，这两座冰川处于软且松散、高度易蚀的细粒沉积岩基础之上，该基础富含小摩擦角的黏土和粉砂。阿汝冰崩之前，在卫星图像上发现了冰川裂缝，推测可能是由分离区阿汝冰川内的
13、力平衡演化引起20。根据气象站观测的气温与降水数据分析21，该区域的快速升温很可能使冰的内部温度升高，热状态改变，再加上降水增加，形成了不利条件。青藏高原内部的其他诸多冰川可能正处于类似环境，冰川稳定性受到极大威胁。因此，选取阿汝1号冰川进行研究，以此了解冰川崩塌机理。2.2计算模型建立FDEM法最重要的优点是可以在坡体从一个稳态转化到另一个稳态的全生命周期内（稳定-微裂隙发育-贯通-拉裂-滑移-再次稳定），允许离散单元体产生位移和转动。因此，为了对阿汝1号冰川发生冰崩的现象进行数值模拟，观察破坏后破碎冰体沿基岩滑移的动态演变规律（图4），选取阿汝1号冰川发生冰崩前的二维断面进行研究，断面高程
14、为6 1234 910 m，长度为8 800 m，X方向大于8 800 m之后为阿鲁错湖。对阿汝 1 号冰川剖面进行建模，单元网格采用常应变平面三角形。图 5 为冰川计算模型几何、网格和边界条件，冰川上部为冰，最小单元尺寸为25 m，下部为软基岩，单元尺寸为100 m。该网表1FDEM模拟参数参数密度/（kg/m3）弹性模量/MPa泊松比黏滞阻尼/（kg/ms）接触罚因子/GPa切向罚因子/GPa抗拉强度/MPaI型断裂能/（J/m2）型断裂能/（J/m2）断裂罚因子/GPa多晶冰试样8508500.35174 310 8238.500.850.549162516504.25加载块7 8502
15、00 0000.293 256 643 487200.020.00/0.5 MPa1.0 MPa2.0 MPa4.0 MPa6.0 MPa图2静力三轴模拟试验应力-应变曲线10.05.00.01-3/MPa0.00.51.01.51/%-3-2-1012体变/%0.5 MPa1.0 MPa2.0 MPa4.0 MPa6.0 MPa0246轴向应变/%图3静力三轴模拟试验体变-应变曲线第10期李 明，等：基于连续-离散模拟方法的冰川失稳过程与机理分析第58卷15格图的计算域由2 755个三角形单元组成，节点数为1 691。底部边界为固定约束，侧向边界为法向约束。2.3计算参数选取由于在计算选取的
16、冰川剖面处，山顶和坡面上无人工建筑，因此不考虑一般构造应力，仅考虑冰川的自重应力（9.8 m/s2）。为了保证计算的稳定性和收敛性，计算时间增量步设置为110-5s。计算参数的选取和双轴试验类似，冰的参数按照多晶冰双轴数值试验选取。根据阿汝1号冰川以往的实地矿物调查，基岩岩性为松散、高度易蚀的细粒沉积岩，故基岩计算参数按照砂岩的特性进行选取。砂岩的密度取2 500 kg/m3，静态弹性特性取自 油层水力压裂22，对于多孔、松散至微胶结砂岩，弹性模量的平均值为 0.7104MPa，泊松比取0.2。另外，抗拉强度取 3.540 MPa23，黏聚力取10.000 MPa，内摩擦角取4524。冰和基岩之间的节理强度适当弱化，反复调整敏感性参数，综合判断数值计算结果的可靠性。综上所述，计算参数如表2和表3所列。2.4FDEM数值计算结果及分析阿汝双冰川的坍塌十分罕见。前人对冰崩发生的原因已有相关研究，KB等25认为，2016年夏季融雪和雨水引起冰川下水文系统的短期变化，基底阻力降低，又因为基岩岩性较软，粗糙度较低，摩擦力减小至极限剪切强度，致使冰川整体性失稳。ZHANG等20在2015年10月2
