基于DEM⁃CFD耦合方法的砂土渗透破坏强度劣化模型研究.pdf

第 卷第 期 年 月水利与建筑工程学报 ，：收稿日期：修稿日期：基金项目：国家自然科学基金重点项目（）；国家重点研发计划项目（）作者简介：邱松楠（），男，硕士研究生，研究方向为岩土力学等方面的研究。：通讯作者：黎晓冬（），女，博士研究生，研究方向为渗透防护措施等方面的研究。：基于 耦合方法的砂土渗透破坏强度劣化模型研究邱松楠，黎晓冬（国家开放大学，北京 ；北京工业大学 城市与工程安全减灾省部共建教育部重点实验室，北京 ；南京理工大学 机械工程学院，江苏 南京 ）摘要：渗透破坏对岩土类材料的冲蚀是造成土坝、桥梁等水利基础设施破坏的主要诱因之一。为了克服现有试验手段在还原渗透破坏发生发展过程中存在较大缺陷，利用离散元软件 中 颗粒流数值方法来模拟颗粒运移过程，并通过设计标准尺寸的砂柱模型进行模拟试验，结果表明，细颗粒含量对模拟砂柱的影响显著，渗透破坏导致砂柱的峰值强度值和临界强度值下降明显，强度劣化度不断增加。理论推导出 耦合劣化模型并加以验证，发现 耦合劣化模型可以在砂柱模型中进行渗透破坏过程计算，在整体趋势和计算精度方面也可得到良好的预测结果。研究成果可为渗透破坏引起的土体力学强度劣化分析及水利工程渗透破坏灾害的防控提供参考。关键词：渗透破坏；劣化度；流固耦合；离散元；数值模拟中图分类号：文献标识码：文章编号：（），（，；，；，）：，：；渗透破坏是水利工程防灾防护领域首要考虑的问题。在我国重大涉水项目工程中，时常出现因渗透破坏导致的管涌、地基沉降、边坡失稳等诸多问题，以雅鲁藏布江引水入疆工程为例，由于地理地质环境因素的限制，在砂土地基下进行隧道开挖、水渠建设需考虑渗透冲刷的影响。目前针对渗透破坏引发的土石坝溃坝的研究方法主要包括试验研究和数值研究。其中试验研究可根据场地条件划分为现场试验和室内试验，苗发盛等 基于渗透 环剪试验研究了滑坡土在不同渗透条件下强度弱化的性质。李军红 研究了非均质多孔介质在渗流耦合因素下细观结构变异特征及宏观强度劣化规律，并分析了细观结构变异对土渗透性和强度特性的影响。由于现场试验在场地选取、费用、时间成本及风险控制等方面存在较大困难，而基于相似理论的室内试验缺乏可靠的土工试验机，整体试验的操作难度高且供排水系统尚不完善。因此，数值模拟基于其成本低、灵活性大、可重复性强的特点，已被广泛应用于管涌破坏、溃坝致灾等水利工程领域。数值模拟方法中，离散元颗粒流 网格流体耦合模型方法（）可以较好的描述特定工况下研究区域内的颗粒分布情况，通过迭代计算速度场与压力场的变化，获得下一阶段颗粒的变化趋势。颗粒流方法很好的解决了固 固、流 固的多相耦合作用过程，相比于有限元程序颗粒流模拟渗透破坏工程具有很大的优势。诸多研究结果表明：采用离散元颗粒流 网格流体耦合模型方法能够模拟从低雷诺流到高雷诺流范围很大的流体运动，模型方法模拟多孔介质在不同雷诺数下的流体渗流和颗粒移动过程的结果真实度较高。倪小东等 利用颗粒流数值方法分析了管涌发展过程中模型内部细观变化规律，并于室内砂槽模型试验结果进行对比，证明了颗粒流程序研究管涌问题的合理性。周健等 对比了颗粒流模型试验与小尺寸细观模型试验结果，进一步验证了颗粒流方法的可行性和合理性。通过比对重塑粉砂试样在不同的时间间隔内涌砂量的变化，建立了涌砂量随时间和水力梯度的变化关系，并将渗透破坏细颗粒迁移模型与用于无侧限渗流分析的有限元程序和应力分析程序相耦合，分析了渗透破坏现象对土体和周围结构的影响。等 建立了二维离散元模型对高渗透过程的简化模拟，分析了渗透破坏过程中土体变形和颗粒重分布的相互作用，指出渗透破坏降低了土体强度，导致土体出现畸形应变。基于微观结构法分析了颗粒的应力 应变关系，发现在高应力条件下，土体产生较大变形，渗透破坏后的土体内摩擦角所获得的数值结果与现场堤坝上观察的破坏模式相一致。此外，离散元计算方法也已经成功应用于砂土的三轴剪切试验等分析中 ，可以有效的解释渗透破坏发生前后砂土力学性能和细观结构的变化。基于以上离散元颗粒流 网格流体耦合模型方法的研究成果，拟建立 渗透破坏数值模型，并通过分析累计涌砂量的变化，推导渗透破坏强度劣化模型，并结合试验结果与前人研究结果对本文劣化模型的适用性进行验证。基于 耦合的数值方法 求解方程采用颗粒流程序进行渗透过程的研究，首先要确保模拟渗流过程与多孔介质中的流动相适应。在多孔介质中，流体在雷诺数较低时应符合达西定律，当雷诺数超过一定值时，流体流动偏离达西定律，流速和水力梯度则成非线性关系。利用 内置 语言可定义流固两相的压力梯度方程和作用力方程，求解不可压缩流体中两相介质的 连续方程和运动方程。流体的作用力施加在砂土颗粒上，砂土颗粒的运动和孔隙率等也间接影响流体运动的参数及其状态，二者充分考虑流固耦合作用，具体的求解方法分为液体相方程和固体相方程 。接触模型颗粒间的接触刚度模型为线性接触模型（），服从广义胡克定律，为方便计算和模拟理想状态下纯砂土的渗流过程，将模拟过程中固体颗粒间粘结力设置为 ，而颗粒间的滑动模型则采用摩擦滑动模型（）。两种接触模型概念如图 所示。图 颗粒切向、法向接触模型水利与建筑工程学报第 卷 迭代程序具体迭代求解过程如图 所示。图 颗粒流数值模型迭代程序饱和砂土体固液耦合反应方程采用 中 （）模型，通过显式差分方式求解。在 时步中，颗粒位置决定颗粒间叠合量，从而确定粒间接触力的大小，将接触力与拖曳力相加，计算颗粒加速度，确定 时步颗粒的位移。流体计算基于交错网格法进行离散计算，流域内的颗粒移动引起土体孔隙率的变化，导致速度场、应力场分布发生改变，流体对颗粒的拖曳力发生改变，影响颗粒在 时步的位移，进而引起孔隙率的改变。渗透破坏的 数值模型 砂体模型此次模拟选定的颗粒粒径范围在 ，其中 的细颗粒含量（）为 ，的细颗粒含量为 ，的细颗粒含量为 ，的细颗粒含量为 ，的细颗粒含量为 。试样的尺寸高度 ，底面直径 。因为细颗粒数目较多，计算效率普遍较低，因此对流体网格划分单元。为了能真实模拟砂样的实际状态，砂样的生成按实际试验步骤进行。通过 函数在给定的模型空间内随机生成一定密实状态的基料，循环消除内力；经过足够长的时间后，砂样内部的不平衡力基本消散，砂样处于相对自然沉积状态。砂样内部共有 个颗粒，如图 所示。为了增加计算时步，提高计算精度，可放大颗粒半径，应用 等 提出的相似理论。在数值模拟分析中采用相似系数为 ，即颗粒放大 倍，各直径的颗粒数目由 等 的方法计算。在砂样饱和后，进行初始耦合计算，使颗粒体达到静水压力状态。施加边界条件：在砂样顶部施加 的竖向压力，保持底部压力为 。根据弹性梁在纯轴向荷载和纯切向荷载作用下的分析，得到颗粒法向和切向刚度数值，具体参数见表 。表 数值模型参数模块名称密度（）法向刚度（）切向刚度（）摩擦系数动力粘度（）颗粒 墙 流体 图 模拟砂样未发生渗透破坏的砂样（设为对照组）在 的围压下固结，固结后设置颗粒间的摩擦因数为 ，然后调整围压大小为试验设定数值。在剪切过程中，未发生渗透破坏的砂样和渗透破坏后的砂样分别在各自恒定围压下，按照恒定 的应变速率进行竖向压缩。墙和颗粒间的摩擦因数取值用于消除颗粒与墙的边界效应。试验方案三维视图下试样的示意图如图 所示，按照选定粒径范围的颗粒组成，四周设为刚性不透水边界，下砂面采用多孔底板模拟可蚀颗粒自由流出边界，上砂面模拟试验室的透水石。在上砂面施加大小为 的压力边界条件，在 的水力梯度条件下进行 耦合运算，使试样最终达到稳定渗流状态。渗流过程中在试样下砂面向下 位置处第 期 邱松楠，等：基于 耦合方法的砂土渗透破坏强度劣化模型研究设置承载板，用于计算流出颗粒数目，即累计涌砂量。计算完成后进行三轴剪切试验过程中，将下砂面的多孔底板替换为墙构件模拟的透水石，渗流过程边界和三轴剪切边界三维视图如图 所示。详细试验方案及参数设置见表 。图 模拟试样三维示意图图 渗流过程边界和三轴剪切边界三维视图表 试验方案试验名称初始细颗粒含量 渗流时间 围压 剪切速度（）结果分析 不同时刻细颗粒的运移情况渗流过程边界和三轴剪切边界三维图如图 所示，图 展示出土体内部细颗粒的运移情况，在试样内部多个位置同时出现细颗粒运移现象，渗流后期试样下部出现裸露的粗颗粒骨架。累计涌砂量随时间呈抛物线型变化，在渗流时间为 达到最大。该现象产生的原因是由于试样细颗粒含量较多且施加较低渗透坡降造成的。由图 可以看出，试样上部和底部颗粒的运移情况较中部颗粒更加明显，试样中部孔隙基本被完全充满，只有下游出口附近的细颗粒运移流失后，才能给上游侧的细颗粒运移提供空间，但由于大量细颗粒同时出现运移，可能造成局部淤堵，进而导致分层运移的情况。图 渗流过程边界和三轴剪切边界三维视图水利与建筑工程学报第 卷 渗透破坏前后的力链变化情况渗透破坏后力链的变化情况主要取决于细颗粒在土体应力传递结构中所起的作用。等 根据土体细颗粒含量大小，将土体应力传递方式分为 种情况。第 种情况，当细颗粒含量较少时，细颗粒松散分布在粗颗粒形成的孔隙中，几乎不参与土体内部应力传递；第 种情况，当细颗粒含量较多时，部分细颗粒开始承担应力，颗粒间形成弱力链，但仍以骨架颗粒为主传递应力；第 种情况，当细颗粒含量达临界值，颗粒足够充满孔隙时，粗细颗粒同时参与土体内部应力传递。上述模拟中 试样属于情况 ，其中部分细颗粒参与应力传递，形成了数量众多、错综复杂的弱力链，随着细颗粒不断流失，该情况下结构破坏方式为弱力链逐渐减少，强力链失去支撑，进而造成骨架应力传递结构失稳。渗透破坏前后的力链变化情况如图 所示，力链中线条的粗细代表接触力的大小，粗线条表示强力链，细线条表示弱力链。由图中 可看出，试样渗透破坏前的强弱力链差别不明显，而渗透破坏后，试样上部弱力链数量减少，同时强力链明显变粗，说明土体内部应力传递结构发生了明显改变。图 渗透破坏前后的力链变化情况 渗透破坏前后的应力 应变情况渗透破坏后模型应力应变特性与力链变化之间存在内在联系，模型内部的结构性损伤必然引起峰值强度的降低，三轴剪切模拟试验结果如图 所示，由图 可以看出，耦合运算方法下渗透破坏过程没有改变砂样的应力应变线的形状，不同渗流时间下的砂样应变软化性依然明显。与渗透破坏前模拟砂样的试验结果相比，软化后临界强度值下降明显，原因是由于模拟球状颗粒在三轴应力条件下快速达到临界状态时，大量力链断裂，承载能力快速下降，从而导致颗粒滑移，强度降低。与此同时，峰值强度也出现较为明显的变化，在 组试验中峰值强度 ，而在 组试验中峰值强度为 ，随着渗透破坏时间的增加，峰值强度劣化度逐渐增加。耦合劣化模型 基于累计涌砂量 的劣化模型引入劣化度 的概念，用于定量的反映土体力学性质变化规律，劣化度 为：前 后 前 （）式中：前为渗透破坏前的峰值强度；后为渗透破坏后的峰值强度。基于上述试验结果，模拟砂样的峰值强度劣化度随渗透破坏时间的增大而增大，且与累计涌砂量呈正比，即单位累计涌砂量损失过程中，强度劣化率为常数。设累计涌砂量产生 克后，砂样峰值强度为（），为可微函数，令未发生渗透破坏的砂样峰值强度为（），则从 克累计涌砂量到（）克的峰值强度劣化度为：（）（）（）（）（）式中：为单位累计涌砂量的峰值强度劣化度。变换得：（）（）（）（）即：第 期 邱松楠，等：基于 耦合方法的砂土渗透破坏强度劣化模型研究图 渗透破坏前后的土体应力 应变关系 （）（）（）对式（）积分可得：（）（）（）（）联立 （）（）（）（）并考虑实际状况下对峰值强度劣化度的影响差别，引入修正系数 对公式进行修正，则有：（）（）由式（）可看出，渗透破坏作用下模拟砂样峰值强度劣化度与累计涌砂量之间呈指数函数关系。但考虑到劣化度并非完全与累计涌砂量服从正比例函数关系，采用试验数据拟合的方式对上述 耦合劣化模型进行验证。劣化模型的验证 等 采用细颗粒控制试验方法，通过使用可溶盐颗粒代替细颗粒的方法，在控制的应力条件下将盐溶解在水中，以实现内部渗透过程期间细颗粒的指定损失。侵蚀过程中的径向和轴向变形均采用照相法测量，并进行了排水三轴压缩试验，研究结果表明损失不同数量细颗粒的土体应力应变行为中峰值摩擦角和临界摩擦角随着细颗粒的流失而减小。的试验结果得到了诸多试验证明和较为广泛的认可，代入其试验结果验证 耦合劣化模型的可靠性。图 计算结果与实测数据对比表 试验方案组别取值取值相对误差 拟合度 组 组 图 与表 反映了 耦合劣化模型对 试验结果的拟合程度，表明劣化度与累计涌砂量具有良好的相关性。结果证明，参数 的取值在附近，参数 在 附近。综上可知，水利与建筑工程学报第 卷耦合劣化模型公式可以在砂柱模型中进行渗透破坏过程计算，相比于双曲线模型，该模型在峰值流量的计算结果与实测值误差较小，在整体趋势和计算精度方面也可得到良好的预测结果。结论（）初始细颗粒含量对模拟砂土的影响显著，会改变土体内部结构和传力方式。具体表现为试样渗透破坏前的强弱力链差别不明显，而渗透破坏后，试样上部弱力链数量显著减少，渗透破坏后土体的峰值强度值和临界强度值下降明显，且随着渗透破坏时间的增加，峰值强度劣化度逐渐增加。（）前述试验结果和模拟结果充分证明，利用 中 颗粒流数值方法，是数值分析研究渗透破坏较直观、可靠的方法，模拟结果的相对误差也较低。（）渗透破坏作用下模拟砂样峰值强度劣化度与累计涌砂量之间呈指数函数关系，采用试验数据拟合的方式对上述 耦合劣化模型进行验证。结果表明：耦合劣化模型可以在砂柱模型中进行渗透破坏过程计算，在整体趋势和计算精度方面也可得到良好的预测结果。参考文献：苗发盛，赵帆程，吴益平，等 基于渗透 环剪试验的三峡库区童家坪滑坡滑带土强度特性研究 岩土工程学报，（）：李军红 细观结构变异引起的垃圾土渗透性及强度劣化研究 石家庄：石家庄铁道大学，倪小东，王媛，王飞 管涌的砂槽试验研究及颗粒流模拟 四川大学学报（工程科学版），（）：周健，周凯敏，姚志雄，等 砂土管涌 滤层防治的离散元数值模拟 水利学报，（）：，（）：，（）：，（）：徐小敏，凌道盛，陈云敏，等 基于线性接触模型的颗粒材料细 宏观弹性常数相关关系研究 岩土工程学报，（）：蒋明镜，胡海军 密实和松散颗粒材料等吸力三轴剪切试验离散元数值模拟 中南大学学报（自然科学版），（）：，（）：，：王恒通，王家全，唐毅，等 组合 颗粒加筋砂土三轴剪切试验离散元模拟分析 广西科技大学学报，（）：，：，（）：，（）：，（）：，（）：第 期 邱松楠，等：基于 耦合方法的砂土渗透破坏强度劣化模型研究