边坡失稳问题的非局部损伤本构建模

张 宇， 卢广达， 夏晓舟

（1.河海大学力学与材料学院工程力学系，南京 211100； 2.同济大学土木建筑工程学院，上海 200092）

边坡滑移破坏的数值模拟一直是工程界和学术界关注的焦点. 除有限元之外，发展了很多专门模拟破坏的数值方法，如界面元［1］、块体元［2］、流形元［2］、光滑粒子流体动力学（SPH）［3］、粒子有限元［4］、光滑粒子有限元［5］、物质点法［6］、无网格法［7］、扩展有限元法［8］以及近场动力学［9-10］. 离散元和块体元需提前定义滑移面，如果滑移面事先未知，则需要不断搜索，通过比较最速能量释放来确定最危险滑移面，在搜索过程中，网格需不断重构，非常繁琐. 流形元则是通过定义在覆盖上的函数进行加权平均来定义区域的位移模式，因而形式上比较灵活，能够构造出非连续的位移模式. 光滑粒子流体动力学（SPH）方法是针对每个粒子进行牛顿力学分析，以追踪粒子的运动轨迹. 粒子无网格法则是在局部范围引入径向基的核函数，因而能够反映应变局部化的破坏现象，但计算量较大，且本征边界条件不能直接施加. 扩展有限元是在传统的位移模式中，利用单位分解的思想，叠加一个反映裂纹处位移间断（或界面处应变间断）的附加不连续函数［11］以及裂尖处应力奇异的附加三角函数［12］，所以无须重构网格就能实现裂纹的追踪模拟，但裂纹的扩展机制受限于断裂力学理论，这对于裂纹分叉和动态裂纹发展仍然是一个挑战，且裂纹在发展过程中，尽管网格不需要重构，但附加自由度则需要随时进行更新以匹配当下的拓扑构型. 近场动力学某种意义上是一种超越了尺度概念的分子动力学模型，数值方法上类似于光滑粒子动力学，但物理内涵更为丰富，具体体现在结构上能够从连续状态自然地过渡到非连续状态. 在近场动力学模型中，每个点的力状态由其邻近点与该点的相对变形所决定，当该点与邻近所有的点发生关联，则表示连续，当该点与邻近所有点切断了联系，则表示在该点处不连续，处在这两种状态之间的为过渡态，即损伤介于0到1之间，因而能够突破连续介质力学中发生奇异或间断的微分运算，如涉及的非局部变形梯度. 因而近场动力学模型能够追踪裂纹的发展而无须更改初始构型的拓扑结构；尽管近场动力学模型是基于积分型的表达，某种程度上抹掉了某些奇异和间断的剧烈差异，但仍然保持着物理上的变形演变势特征，以驱动裂纹的发生. 实际中，理论上的应力奇异是不存在的，裂纹处的应变也是不存在的，因而非局部变形梯度的引入大大拓展了变形场的适用范围，使得拓扑结构不断变化的问题（如裂纹发展）能够在统一的框架下（不增加自由度）进行. 当然，近场动力学模型的主要问题是计算量大，存在边界效应，积分项中需要引入高阶微分算子才能达到精度上的要求，否则存在一定的数值振荡［13］，而与有限元混合建模是一个很不错的出路.

本构方面，主要分为局部型本构模型和非局部型的本构模型两类. 传统的塑性模型、损伤模型和弹塑性损伤模型都属于局部型本构模型，材料参数都是基于一定尺寸的试样测试的，因而参数和赋予的应力-应变关系曲线都是基于试样平均化获得的. 因此，在处理材料屈服后的力学行为时，局部化模型难以准确地反映应力释放和转移机制，也就很难捕捉应变局部化的效果. 应变梯度理论［14］（包括偶极理论，如Cosserat理论［15］）中物质点具有一定的尺寸大小（通常取材料的内禀尺度），能够反映物质点的微弯曲和微扭转机制，属于非局部化模型. 应变梯度理论模型对于裂尖场的变形能够精细地反映，但该理论模型最大的问题是物质点的内禀尺度不好确定，且在广义应力空间（包括偶应力）中构建屈服面时很难通过实验进行标定. 另一种非局部模型，就是通过引入非局部积分算子，对要关注的物理量（如塑性应变、损伤）进行加权平均化处理［16-17］，这种积分型的非局部化模型不仅能够克服网格依赖性问题，且某种程度上能够反映破坏过程中的应变局部化现象.

本文借鉴近场动力学中点对的思想，把微损伤定义在点对上，并在考察点的近场范围内通过加权平均化处理对该点的整体损伤进行评估，最后引入损伤导致的刚度退化模式，构建一个非局部化的损伤本构模型，并嵌入到有限元的计算框架下，实现边坡的滑移破坏模拟.

文章分为5个部分，在接下来的第2部分详细介绍非局部损伤本构模型的构建过程，包括键损伤模型、点的整体损伤评估和损伤导致的刚度退化模式；第3部分将展示多尺度有限元建模理念及通过四点剪切梁进行剪切荷载导致的破坏现象进行验证，以便应用到边坡的滑移破坏模拟，即第4部分的内容；第5部分则针对本文构建的模型进行深入的讨论和总结.

1 非局部损伤本构建模

模型任意一点的受力状态应该由其周边点与该点的相对位移所决定. 另外，实验给定的应力-应变关系曲线也是基于试样均匀化获得的，因此当试样内部出现微裂纹时，所发生应力释放和转移在局部模型中是反映不出来的，而这正是应变局部化的根源，所以在传统的损伤本构模型中呈现的损伤分布不够集中. 为此，本文受近场动力学非局部作用思想的启发，构建一个基于邻域点对（键）变形的非局部损伤本构模型.

1.1 键损伤模型

对于任一点x，考察其周边点x′与点x 形成的键的变形损伤情况. 定义点对x′x 的正伸长率ρ+为

图1 物质点对的相对运动Fig.1 The relative motion of a pair of points

在键的损伤演化方程中，损伤的发展速率受历史最大超越伸长率及控制参数γ影响，γ值取得越大，损伤发展就越快，图2 给出了不同γ 所对应的损伤随历史最大超越伸长率的演化曲线.

1.2 点的损伤状态评估

由点影响域内键损伤情况，对该点的整体损伤Ω 进行加权评估如下：

图2 键损伤演化曲线Fig.2 Damage evolution curves of bond

式中：V′为x点所占据的空间大小；φ( x′,x )为键损伤的权函数，且要求∫H(x)φ( x′,x )dV′=1.0. 由正伸长率的定义式（1）可知，对于均质材料，影响域内距离点x 最远的点x′形成的点对x′x 键更容易断裂，即影响域内远处点较近处点优先脱离影响范围，这与近场动力学的键伸长率定义不一样.

1.3 刚度退化模式

由胡克定律得：

式中：De为初始弹性矩阵；Ded为弹性损伤矩阵；εE为弹损伤应变；εe为弹性应变；εed为微裂纹闭合的那部分应变. 如果把介质区分为损伤部分和未损伤部分，则对于未损伤部分，有

由式（5）、（7）～（9），导出弹损伤矩阵与弹性矩阵之间的关系为

这就是二次刚度退化的模式，但实际情况是，随着损伤的发展，材料点的影响域范围在缩小，即影响半径l 随着影响域内部键的断裂而减小，可根据有效面积的减小定义影响半径随损失的演化关系为

式中：l0为点x 的初始影响域半径. 由公式（1）可知，影响半径的减小导致伸长率的增大，从而由键损伤判据（3）式和宏观损伤（4）式可获得这样的结论：随着损伤的发展，键更加容易断裂，即材料的抗力随着损伤的发展而逐步衰减，从而加速了刚度的退化，即宏观损伤Ω 与键损伤ω 是相互促进的.

2 模型验证

2.1 多尺度有限元建模

在大部分模型中，出现破坏的地方总是集中在很狭小的局部带中，因此只需在容易出现破坏的范围内划分足够密的网格，且采用提出的非局部损伤本构模型；其他部分的网格剖分只要保持工程精度即可，且采用弹性本构模型. 这样既能保持足够的精度，又能提高计算效率.

对于非局部损伤本构模型，为确保点的整体损伤评估精度，在给定近场范围，点对数量越多越好，至少保证每个点的近场范围内有40 个点对以上.另外，由于单元刚度矩阵是基于积分点进行加权求和的，因此，本文的点对直接通过单元积分点进行构建，为方便示意，以三角形单元为例，且每个单元布置一个积分点，如图3所示.

2.2 四点剪切梁

下面以四点剪切梁为例说明本文提出的非局部损伤本构模型是否适用剪切型的破坏模式. 梁模型、约束和加载如图4所示，梁的长度L=0.8 m，梁高H=0.2 m，开口裂缝长a=0.04 m，支座间距b=0.08 m.估算损伤破裂范围，并进行多尺度网格剖分，如图5所示. 对细网格单元，应用构建的非局部损伤本构模型，其中材料弹性模量E=40 GPa，泊松比ν=0.2，临界伸长率ρc=1.5e-4，控制参数γ=20.采用弧长法迭代，直到梁结构不能继续承载时结束计算，获得的损伤发展如图6所示，从图中可以看出，损伤从开口裂缝处发展，发展路径和试验破坏形态（图7）是一致的，且损伤分布比较集中，说明本文提出的非局部损伤本构模型是合理可靠的，能够模拟出剪切荷载的破坏形态.

图3 有限元空间离散以及物质点对示意Fig.3 Discretization of finite element space and point pairs of matter

图4 四点剪切梁Fig.4 Four-point shear beam

图5 梁的多尺度网格剖分Fig.5 Multi-scale meshing of beams

图6 四点剪切梁的损伤发展Fig.6 Damage development of four-point shear beam

图7 试验破坏后的形态Fig.7 Shape after test failure

3 边坡失稳模拟

某土质边坡，几何尺寸如图8所示，多尺度网格剖分如图9所示，弹性模量E=1 000.0 MPa，泊松比ν=0.45，临界伸长率ρc=1e-4，控制参数γ=10 . 文献［18］采用通过辅助刚块施加压力荷载P，本文采用位移加载模式，即对刚块接触面的节点施加相同竖向位移荷载，计算后获得的损伤分布如图10 所示，从图中可以看出一条明显的损伤滑移带，呈弧形状，与文献结果相同. 说明本文提出的方法完全适用于剪切滑移破坏问题，而且实现起来简单、高效.

图8 某土质边坡约束及受载Fig.8 Restraint and load of a soil slope

图9 土质边坡的多尺度网格剖分Fig.9 Multi-scale meshing of soil slope

图10 边坡的滑坡模拟Fig.10 Slope landslide simulation

4 讨论与结语

本文受近场动力学的思想启发，构建了一个非局部的损伤本构模型，该模型突破了连续介质力学中应力-应变关系均匀化标定的限制，这种均匀化标定对于材料出现破裂后是不太适合的，因为涉及应力释放和转移，这就是为什么传统损伤本构模型存在网格依赖性问题的原因. Bazant曾通过非局部化算子对损伤进行加权平均处理［17］，使得任意点的力学状态同其周边点联系起来，这样试样的均匀化应力-应变曲线赋予了该点的非局部范围，而不是该点自身，但在Bazant模型中容易出现应力锁死，即残余应力由于局部平均化的处理，即便损伤达到1.0时也消除不了. 本文提出的非局部损伤模型，其损伤本质上是不连续性程度的表征，因而损伤分布代表的就是缺陷的发展，而在传统唯象损伤模型［19］中，所谓损伤的分布实际上是退化刚度分布. 因而基于非局部损伤的严密定义，结合宏观损伤（非局部）与键损伤相互促进的影响机制，建立了随损伤发展的刚度退化模型和基于影响域范围逐步缩小的抗力衰减模式，很好地反映了结构在破坏过程中的应力释放和转移机制，使得一旦某个单元的损伤发生，其对节点力的贡献瞬间减弱很多，导致节点失衡速率加快，向损伤单元反向运动，形成对损伤单元的进一步拉扯或错动，达到应变局部化的效果.

相比于近场动力学模型，本文保留了连续介质力学框架的一些特性，比如弹性矩阵（包括弹损伤矩阵）天然地反映点对之间的相互作用，这是键基近场动力学模型所忽视的，尽管后面的近场动力学通过各种办法突破了泊松比的某些限制，但对于高泊松比材料如本文算例中的土质边坡，仍然不适用. 当然，如果采用态基的近场动力学模型，把连续介质力学的材料非线性本构模型转化到近场动力学中，也能达到边坡的滑移破坏效果［20］，但模型复杂，且不能如本文方法在统一的数值方法中（如有限元）实现多尺度模拟，达到高效计算的目的.

本文方法简单高效，容易拓展到多物理场的耦合计算；对于黏弹性损伤和蠕变问题，也容易拓展，只需对键损伤模型进行改造即可，因而本文方法具有重要的理论意义.