基于FEM-SPH耦合方法的挡土墙墙后土体变形破损研究.pdf

引用格式：涂正楠，吴羿君，冯君，罗晗玲，杨涛 基于 耦合方法的挡土墙墙后土体变形破损研究 山地学报，（）：，（）：山 地 学 报 文章编号：（）：收稿日期（）：；改回日期（）：基金项目（）：四川省交通运输科技项目（）；四川省科技计划项目（，）。（）；（，）作者简介（）：涂正楠（），男，重庆江津人，硕士研究生，主要研究方向：岩土工程数值模拟。（），：通讯作者（）：冯君（）男，四川渠县人，博士，副教授，主要研究方向：岩土力学与边坡工程。（），：基于 耦合方法的挡土墙墙后土体变形破损研究涂正楠，吴羿君，冯 君，罗晗玲，杨 涛（西南交通大学 土木工程学院，成都 ）摘要：挡土墙发生较大位移时，墙后土体往往出现裂缝，应用传统有限元法（）计算墙后土体受力变形时存在困难。的耦合方法结合了 的高精度、高效率和 粒子处理大变形的能力，可以较好地模拟裂缝的衍生与发展演化。本文介绍了 方法的基本原理与 的耦合方法，使用 自适应耦合算法对平动位移模式下挡土墙墙后土体的裂缝衍生与墙后土压力进行了三维数值计算，并与试验结果进行了对比，结果表明：（）方法可以较好地模拟挡土墙墙后土体的裂缝衍生；（）墙后滑坡推力峰值在距离墙顶 处；（）挡土墙平动模式下，张拉裂缝与剪切裂缝交替出现，张拉裂缝的出现由远挡墙向近挡墙的顺序衍生，剪切裂缝的出现由土体下端向上端的顺序衍生。研究结果可以为解决岩土体大变形问题提供新的思路。关键词：光滑粒子流体动力学；挡土墙；模型试验；数值模拟；裂缝衍生中图分类号：文献标志码：挡土墙是支撑路基填土或防止填土或土体变形失稳的构造物。挡土墙上的土压力大小及分布规律，是挡土墙结构设计的基础。对于挡土墙应力应变、墙后土压力计算等研究在理论推导 、数值模拟 和模型试验 ，等方面已经取得了很多的成果。当挡土墙底部在位移较大时，墙后土体会产生裂缝，裂缝的产生会对模型计算网格造成较大的干扰，为了减少这种网格大变形带来的干扰，有学者 提出了使用插值无网格和弥散裂缝模型结合的方法模拟土体裂缝的产生与发展，但是该种无网格法需要较密的网格节点，存在精度不够高、计算效率低等问题。传统的 （）方法精度较高，但是其在处理大变形问题，尤其是土体这种抗拉强度较低的材料时，通常采用网格删除技术模拟开裂变形，其计算精度必然有所影响。为了解决大变形 对 计 算 结 果 产 生 的 干 扰，拉 格 朗 日 方 法（）和拉格朗日 欧拉法（）可以对网格进行重划分，但这会大大增加计算量并会导致部分数据的丢失 。无网格法是一种新兴的数值计算方法，因其在插值时不需要网格，仅需要计算在节点通过核函数构造的影响域，克服了依赖于网格的数值方法的一些缺点，在处理大变形及裂缝衍生等问题时具有明显优势。光滑粒子流体动力学方法（，）作为一种常用的无网格方法，前处理简单，不需要繁琐且冗杂的代码，同时兼备高精度与高效率而被广泛应用于大变形研究。岩土体中常用的本构关系 （）、（）、均能适用于 方法。由于 方法不需要网格划分和小变形假设的特性，其可以不受到网格畸变的影响，每个生成的粒子被单独赋予状态变量（质量、速度、位移、应力、应变等）和材料属性（重度、弹性模量、泊松比、内摩擦角、剪胀角、黏聚力等），方法 被 广 泛 应 用 于 滑 坡 模 拟、灾 害 防 治 预 测等 。目前已经有学者开始使用 法结合 方法解决工程中有关爆炸 和碰撞 的问题，但存在人为主观因素划分区域问题，不涉及 到 的转换。而土体为单一材料，人为划分大变形区域会有失精度。虽然 方法对于通常细微部件的材料变形和结构响应的计算分析精度不及同等尺度的 方法，但对于存在大变形的结构时，方法是通过解除粒子之间的约束来模拟开裂变形，不存在材料的删除缺失，故 模型的计算分析精度高于同等尺度的 模型。的耦合方法在有限元网格出现大变形的时候根据一定转化准则自适应转换为 粒子，结合了 的高精度、高效率和 粒子处理大变形的能力，可以模拟裂缝的衍生与发展演化，尝试解决传统方法无法准确模拟挡土墙受力变形的问题。本文介绍了 方法的基本原理与 的耦合方法，使用 的自适应耦合方法计算了墙后土压力的分布，再现了挡土墙墙后土体裂缝的衍生过程，与试验结果吻合较好，显示了 方法在大变形问题上计算精度的优势。理论基础与 耦合?方法基本理论光滑粒子流体动力学最早由 和 与 提出并应用于模拟星系动力学。随后，方法被广泛应用于液体、气体等流体的模拟，成为一种强大而灵活的工具。方法使用一组离散的粒子来表示流体，每个粒子具有独立的质量、位置、速度等属性。流体被视为由这些粒子组成的集合，通过对粒子的属性进行插值来模拟流体的性质。方法的构造方程主要由两个关键步骤组成 ：场函数核近似法、粒子近似法。在第一步中应用积分表示函数的核近似法对光滑核函数进行逐步积分，对于点 上的粒子的支持域内任意一点 ，其函数表达式 （）可以写为：（）（）（，）（）式中，为光滑核函数（），或简称为核（）（图 ），基点粒子的变量值可以通过其影响区域内相邻粒子的集合的贡献值和来近似；为光滑长度（），其决定了有多少粒子可以影响特定点的插值；为 的积分域。为粒子 的影响域；为比例因子，其决定了光滑核函数的作用范围。图 核函数 本文采用的核函数为 年 和 在三次样条函数基础上改进的 样条函数，其被广泛应用于 方法的计算，其表达式如下：（，）（，）（）（）式中，为满足光滑核函数正则化条件的系数，在山地学报 卷 第 期一维、二维、三维空间中分别有，；是在 和 处两粒子之间的相对距离，其中 是两点之间的距离。在第二步中粒子近似法将连续介质划分为一组具有特定属性的离散粒子，如质量、位置、速度等。每个粒子代表了一部分连续介质的特征。通过对粒子之间的相互作用建立数值模型，可以模拟和计算连续介质中的物理过程，每个粒子具有独立的质量、位置、速度、应变等属性。这些属性可以根据物理规律和数值模型进行更新和求解。然而，传统的 方法在边界区域存在粒子缺陷，积分精度受限导致求解不准确。为了改善粒子在边界处的不稳定性，国外学者 提出了不同的方法，其中，基于虚粒子的边界处理方法 是最常用的方法。当 粒子接近拉格朗日边界面时，通过平均相邻边界面形成一个虚拟平面。在虚拟平面的基础上，引入虚拟粒子，提高积分精度。改进边界条件下（图），表面两侧的 粒子不再相互作用。图 无滑移边界条件算法示意图（基于文献 修改）（）?耦合方法 法和 方法是计算力学中的两种重要的数值方法，由于其高精度的优点，法在工程应用中被广泛使用。然而，法难以解决网格畸变等大变形问题，而 方法存在拉伸不稳定性和应用边界条件施加时存在困难等限制。因此，一种常见的策略是在存在局部大变形的区域使用 方法，而其余区域使用计算代价较小精度更高的 法。通过施加适当的连续 非连续转换触发机图 与 接触示意图 制，可以实现 的耦合，从而完成对结构体力学变形行为全过程的精确模拟。与 之间的耦合主要依赖以下四个方面实现：（）接触面耦合：在 与 接触面处，通过接触力来耦合 和 ；（）边界条件耦合：在 边界处，通过对 的边界条件进行耦合；（）动力学耦合：通过耦合流体和固体的动力学方程，使得两者共同影响对方的运动；（）网格耦合：通过在网格将 和 耦合在一起，使得两者可以互相影响。为了实现 与 粒子之间的相互作用，在自适应 方法中对接触界面采用接触算法进行计算 ，有限元网格和 粒子之间的接触面（图 ），有限元节点上的 粒子为虚线表示的小圆形，虚拟粒子 的定义域在虚线圆形范围内。伴随着虚拟粒子运动到距离有限元节点两倍光滑长度，有限元节点和虚拟粒子间将产生一定大小的接触力，根据无网格粒子接触算法 ，各有限元节点可当成虚拟粒子，当粒子处在有限元节点支持域中，节点将受到 粒子引起的接触力，节点和粒子的接触力计算表达式为：（）（）（）（）（）式中，为支持域内的粒子（节点）数；为粒子（节点）的质量（）；为粒子（节点）的质量（）；为粒子（节点）的密度（）；为粒子（节点）的密度（）；代表接触 ，基于 耦合方法的挡土墙墙后土体变形破损研究图 粒子转换原理 算法的光滑核函数；是粒子间距（）；是粒子间光滑长度的平均值（）；是用户定义的标量接触刚度惩罚参数；是粒子 支持域内的粒子个数。本文采用 求解器探讨 自适应耦合方法，研究挡土墙墙后土体变形破损和裂缝发展情况，在挡土墙发生位移后，墙后土体塑性应变增大，基于塑性应变阈值的转化准则下，网格中单元塑性应变超过阈值后，损伤区单元开始由 网格转化为 粒子（图），其网格上原有的力学性质继承到 粒子上，破损土体与颗粒间的摩擦系数重生成，逐渐形成一条塑性贯通带，通常 方法在此时由于网格畸变严重导致计算终止，而 由于其无网格特性可以继续计算，墙后土体会产生第二条第三条剪切塑性带，直到计算结束，通过 自适应耦合方法可以深入探讨土体大变形后的裂纹发展及变形破损，方法转换流程如图 所示。模型试验?模型试验设计试验地点在西南交通大学岩土中心实验室，模型试验装置包括模型箱、数字摄影设备和数据采集系统（图），模型箱内部的净尺寸为 ，滑坡模型的滑体长为 ，宽为 ，模型试验的滑体土按砂土：石英砂：膨润土 的比例人工搅拌均匀，在搅拌过程中土体的含水率保持在 左右，根据直剪试验结果，模型滑体土的黏聚力为 ，内摩擦角为 ，基岩按砂图 方法转换流程图 图 试验模型设备：（）试验模型箱；（）摄影装置；（）数据采集装置 ：（）；（）；（）土：石英砂：膨润土 的比例人工搅拌均匀，含水率与滑体土保持一致，其中基岩土的黏聚力为 ，内摩擦角为 ，模型的滑体采用分层填筑的方法，每层滑体按 厚度进行压实直至指定高度。完成滑体土堆筑后，削坡至设计坡型，使用小山地学报 卷 第 期铁锹从外及里逐渐修整边坡，以减少对滑体土的扰动，其中滑坡模型尺寸及土压力盒的布测点埋置位置（图 ），在滑体与基岩的接触面位置处垫上一面光滑的聚乙烯薄膜，模拟模型的滑带。试验设计挡土墙向滑坡模型水平移动 ，试验中控制挡墙的水平移动速度，通过摄像技术及土压力传感器实时记录推动过程中墙后测点的应力变化，分析墙后土体的破坏形式及裂缝发展等。图 裂缝衍生过程及最终状态示意图：（）、（）、（）为裂缝 、的初始状态；（）、（）、（）为裂缝 、的最终状态 ：（），（）（），；（），（）（），?试验结果运用高清摄像技术，分析滑坡模型在整个试验过程中的位移和变形，通过滑坡模型试验架顶部机位拍摄的图像，研究挡墙后部滑体坡面的变形破损及裂缝发展与分布情况。裂缝衍生过程及最终状态结果如图 所示，滑坡坡面裂缝到墙体的距离及土体隆起高度如图 所示。由图像结果可知，在试验的初始阶段，距离挡墙 处的张拉裂缝首先发育，随着变形的发展，距离墙后 、处依次出现第二条、第三条裂缝，最终墙后大部分土体出现隆起现象，隆起高度约为 ，除了第一条张拉裂缝外未记录到明显贯通带。本次试验采用的微型土压力盒为丹东市虬龙传感器有限公司制作的 型土压力传感器，在试验过程中，可通过应变数据采集仪和土压力盒监测墙后土压力的变化情况，应变数据采集仪的采集间隔为 ，模型试验采集时长为 ，微型土压力盒图 模型箱尺寸与测点示意图：（）模型箱尺寸；（）测点布置 ：（）；（）埋置在挡墙不同位置处，用于监测挡墙后部应力变化情况，处理后得到各测点应力变化与时间的关系如图 所示。，基于 耦合方法的挡土墙墙后土体变形破损研究图 裂缝空间位置及最大隆起高度：（）裂缝空间位置；（）土体最大隆起高度 ：（）；（）图 土压力盒读数 图 计算模型示意图：（）模型尺寸图；（）模型网格图 ：（）；（）数值模拟建立与结果分析?数值模拟的建立参照模型试验建立数值计算模型如图 所示，网格均采用 ，共划分 个单元，本模型不考虑挡土墙的变形，因此赋予了挡土墙极大的弹性模量，挡土墙与土体接触面光滑无摩擦，土体与基岩间的摩擦系数按照 （）取值 ，为了使应力能在不同的部件中传递，在土体、基岩 与基岩 间采用绑定约束，粒子自适应最大应变转化阈值为 ，即土体 网格的最大应变达到 时视为变形过大，网格开始转化为 粒子，土体与基岩采用岩土工程中常用的 本构模型，模型参数如表 所示。数值模型系统采用幅值加载，重力在 内完山地学报 卷 第 期成线性加载，挡土墙在 内完成线性水平移动 ，换算速度为 。表 模型结构参数表 模型材料密度（）弹性模量（）泊松比黏聚力（）内摩擦角（）墙体 土体 基岩 基岩 图 裂缝衍生过程?裂缝衍生过程分析图 展示了挡土墙后土体变形及裂缝衍生的过程。在 时，重力加载完成，挡墙开始沿 方向移动。在第 时，墙后土体较远处出现了第一条张拉裂缝（平行于坡面）。这条张拉裂缝将土体切割为两个“土块”。同时，从墙体中部到上一个“土块”角点连线出现了一条与张拉裂缝相接的剪切裂缝（垂直于坡面）。剪切裂缝又将上一个“土块”切割为一大一小两个新的“土块”，第二条张拉裂缝开始发育。在第 ，第二条张拉裂缝发育完全，将上一个大“土块”又切割为两个新的“土块”。第一条剪切裂缝连通两条张拉裂缝。在第 时，第三条张拉裂缝发育完全。由于一条剪切裂缝的存在，第三条张拉裂缝没有贯穿土体形成贯通剪切带。第 时，第二条剪切裂缝从坡顶向张拉裂缝带蔓延，直到贯通第一条张拉裂缝，随着挡土墙的继续运动，在第 贯通第二条张拉裂缝。最终裂缝与土体隆起高度如图 所示。根据数值模拟结果，建立沿挡土墙后直线路径如图 所示的隆起高度数据点，突变处为 和 的分界线，即裂缝产生位置。最终土体隆起最大高度约为 ，三条剪切裂缝到挡土墙的直线距离依次为 、。数值模拟结果与试验结果的误差如图 ，结果表明除第二条裂缝外，数值模拟与模型试验的误差在 左右。试验中未明显观测到剪切裂缝的出现，这是由于剪切裂缝可能出现在边坡内部，还没有水平衍生蔓延到模型箱玻璃两侧。数值模拟能拟合出模型试验裂缝产生的位置与衍生趋势。但由于动力学计算导致网格在边界处存在奇异性，转化为 粒子后粒子间接触面积增大，进而导致了更大的位移，因此 计算结果会略大于试验值。?墙后滑坡推力分析挡土墙墙后土压力沿深度的分布情况如图 ，基于 耦合方法的挡土墙墙后土体变形破损研究图 最终裂缝状态及土体隆起高度：（）最终裂缝状态；（）土体隆起高度 ：（）；（）图 隆起高度与误差分析：（）土体隆起高度；（）模拟值与实验值误差分析 ：（）；（）所示。传统 法在计算深度 以下时与试验结果具有较好的吻合度，而深度超过 时由于下部土体网格扭曲变形过大，应力出现突变，结果已失真，不具有参考意义。方法计算结果与试验结果的趋势和数值吻合较好，展现出该算法的独特优势。墙后土压力峰值与试验结果相似，且均在距离挡土墙顶端约 处，说明对于存在大变形的岩土体问题，使用 方法能够更好地模拟挡土墙滑坡推力与墙后土体的变形破损问题。讨论传统 方法在处理大变形问题时往往由于网格畸变导致计算中止，难以得到可靠的结果，本文基于 基本理论提出一种 自适应耦合方法，通过对比室内模型试验与数值模拟裂缝衍生过程和滑坡推力的计算结果，验证了该方法的准确性，在保证精度的前提下为解决岩土体大变形问题提供了一条新的思路。同时通过该方法可以初步预图 墙后土压力沿深度方向分布 测挡土墙后土体裂缝出现位置，可为滑坡防治工程优化提供指导。本文虽然证明了 方法在计算岩土体山地学报 卷 第 期大变形问题的准确性，但是考虑到土体与基岩间在切向方向应当存在黏结相互作用，仅使用简化后的库伦摩擦 有失精度，在未来的研究中应充分考虑各土层间的相互作用，此外对于 自适应转化阈值的研究，即土体塑性应变达到多大时可以视为开裂状态 ，还需要更多研究加以完善。结论（）根据 方法的基本理论，引入 粒子与 节点间的接触理论，提出了 的自适应耦合方法，探究了从 网格到 粒子的转化原理与实现过程，通过岩土工程常用的 弹塑性本构模型，可以有效仿真挡土墙后土体的变形破损现象。（）模型试验与数值模拟在裂缝衍生现象及墙后滑坡推力的监测结果趋势吻合较好，滑坡推力沿挡墙方向先增后减，峰值在距墙顶 处，验证了数值模拟模型的可靠性与 方法进行三维挡土墙墙后土体变形模拟的可行性。（）挡土墙在水平位移模式下，张拉裂缝与剪切裂缝交替出现将土体切割为“土块”，新裂缝出现位置往往是连接上一个“土块”的一个角点与一条边界的中点，张拉裂缝的出现由远挡墙方向向近挡墙方向的顺序衍生，剪切裂缝的出现由土体下端向上端的顺序衍生。参考文献（）王鸿兴，孙大庆挡土墙后土体滑裂面及土压力变分法初探 岩土工程学报，（）：，（）：党发宁，张乐，王旭，等基于弹性理论的有限位移条件下挡土墙上土压力解析 岩石力学与工程学报，（）：，（）：牟太平，孙伟，邱志雄，等基于有限元近场动力学耦合方法的挡土墙变形破损分析 中山大学学报（自然科学版）（中英文），（）：，（）：周健，彭述权，樊玲刚性挡土墙主动土压力颗粒流模拟 岩土力学，（）：，（）：，（）：，（）：（）杨山奇，卢坤林，史克宝，等刚性挡土墙后三维被动滑裂面的模型试 验 岩 土 力 学，（）：，（）：顾慰慈，武全社，陈卫平挡土墙墙背填土中滑裂体形状的试验研究 岩土工程学报，（）：，（）：彭罛，袁会娜，张丙印无网格自动加密方法及其在土体裂缝分析中的应用 工程力学，（）：，（）：张琰高土石坝张拉裂缝开展机理研究与数值模拟 北京：清华大学，：，：，（）：，：，：，：，基于 耦合方法的挡土墙墙后土体变形破损研究 ，（）：，（），（）：，：，：，（）：，：胡?，谢谟文，王立伟基于弹塑性土体本构模型的滑坡运动过程 模拟 岩土工程学报，（）：，（）：，（）：陈诚，詹发民，周方毅，等 基于 算法的钢板接触爆破数值模拟研究 计算机仿真，（）：，（）：马上，王振清，陈叶青，等 基于 方法的水下爆炸重力坝数值模型构建 防护工程，（）：，（）：王翔宇，李志雨，孙树政，等 基于 方法的舷侧与冰山碰撞结构响应 船舶工程，（）：，（）：，（）：，：，：，光滑粒子流体动力学 一种无网格粒子法 韩旭，杨刚，强洪夫，译长沙：湖南大学出版社，：，：，：，：，：，（）：，：，（）：，（）：张志春，强洪夫，高巍然 接触算法在冲击动力学数值计算中的应用 固体力学学报，（）：，（）：，（）：费康，张建伟 在岩土工程中的应用 北京：中国水利水电出版社，：，：，：，：，山地学报 卷 第 期 ：，（）：，：，（，）：，（）（），：（）（）（），：；（责任编辑朱颖彦钟雨倩），基于 耦合方法的挡土墙墙后土体变形破损研究