结构性土坡变形破坏数值模拟方法的对比分析

汪儒鸿 ，周海清，彭国园

(1.陆军勤务学院军事设施系, 重庆 401331；2.西北核技术研究所, 陕西 西安 710024)

0 引 言

自然界中的土大都具有结构性[1]。结构性的存在对土的力学特性主要有两方面的影响：(1)结构性土较重塑土具有更大的强度和模量；(2)结构性的破损衰化使土体产生较大的体积变形。结构性土坡失稳破坏往往具有突发、快速等特点，如花岗岩残积土斜坡、黄土斜坡等[2-4]。这类土坡的破坏特点与一般非结构性土坡以滑坡为主的渐进破坏[5]区别较大(图1)。

图1 云南省保山市一处结构性土滑坡Fig.1 A granite residual soil landslide in Baoshan City, Yunnan Province

土体结构性的破坏是结构性土坡突然失稳的主要原因。但是传统的极限平衡分析法难以反映结构性土强度变化的非线性特点[6]，不少学者尝试借助数值模拟的方法来对这类特殊土边坡的稳定性加以分析：陈昌禄等[7]将结构性参数引入到黄土边坡的稳定分析中，对传统莫尔-库伦公式进行了修正，结合有限元强度折减法来求解结构性土边坡的安全系数，在一定程度上反映了人工开挖条件下黄土结构性对边坡稳定性产生的影响；李荣建等[8-9]以此为基础，先后考虑了初始结构性参数以及黄土的抗拉性对边坡稳定性的影响，对陈昌禄提出的公式进行了修正，进一步发展了引入结构性参数的边坡数值分析方法；此外，周家伍等[10]研究同样遵循了这一模式，即在提出新颖本构模型的基础上，借助强度折减法来对结构性边坡的稳定性加以研究。以上都是借助有限元或有限差分法的数值模拟研究。除此以外，目前分析结构性土坡稳定性的数值模拟方法还包括颗粒离散元方法。颗粒离散元方法(PFC) 属于定量评价方法中的一种[11]，我国学者已经对 PFC 在边坡工程中的应用开展了大量的研究：周健[12]对砂性土坡和黏性土坡分别进行数值模拟，分析了不同微观参数对土坡破坏型式的影响; 崔溦[13]等以大型复杂堆积边坡为依托，通过室内试验与模拟实验结果对比确定边坡土体的细观力学参数，进而通过建立堆积边坡离散元模型分析了边坡的稳定性，确定了边坡的破坏方式。

总的来说，数值模拟方法确实有助于结构性土坡失稳问题的研究，但由于相关数值分析方法尚处发展阶段，基于不同软件开展的研究还存在不少问题。因此，研究选取最适合分析结构性土坡失稳问题的数值模拟软件，并找出最有效的分析方法，对研究土体结构性对边坡稳定性的影响及合理评估不同工况下结构性土坡的安全系数有着重要经济意义与应用价值。

1 结构性土坡数值模拟研究的关键点

花岗岩残积土边坡及黄土边坡存在的特殊破坏模式与普通滑坡最大的不同在于它的突发性，即失稳前缺乏明显的前期变形预兆，因此可将这类破坏称之为边坡的突变失稳。而土体的结构性对突变失稳有什么影响、结构性土坡是否都存在突变特点则应该是数值模拟研究中的重点。结合结构性土的强度变化特性，在展开结构性土坡稳定性方面的研究时必须突出以下两个关键点：

(1)土体结构性的反映。即如何在数值模拟过程中合理反映结构性土的强度特点随微观结构的破坏而发生的变化。

(2)边坡失稳突发性的反映。即如何在数值分析过程中模拟出边坡突变失稳的破坏效果，并将其与结构性土坡的土体参数量化结合起来。

2 连续和非连续介质数值模拟方法在结构性土体中的应用

2.1 连续介质分析方法

以PLAXIS、FLAC3D为代表的有限元、有限差分软件在对对象进行计算时必须先由使用者指定材料的宏观本构模型。在不同的宏观本构模型中，土体强度的变化可以是线性的，也可以是非线性的，但都只能遵循系统内置的计算模式。目前主要的内置模型包括： 1.空单元模型(开挖模型) 2.三种弹性模型(各向同性、正交各向异性和横向各向同性)3.七种塑性模型(Drucker-Prager模型、摩尔-库伦模型、应变硬化/软化模型、多节理模型、双线性应变硬化/软化多节理模型、D-Y模型和修正的剑桥模型)。主流有限元、有限差分软件的内置本构模型基本涵盖了大多数土体的强度变化特点，但都没有针对结构性土所开发的模型，这与结构性土复杂的强度变化特性有关。针对这一问题，以陈昌禄等[7]、李荣建等[8-9]为代表的学者专门就结构性土的数值计算模型进行了开发，并取得了一定的成果，但存在问题如下：

(1)涉及结构性土岩土本构关系的理论尚处发展阶段，更未达成共识，基于某一类观点所开发的数值计算模型不一定准确。

(2)有限元、有限差分软件所定义的材料本构模型只能反映包括黏聚力在内的宏观土体参数的变化，并不能直观反映土体微观结构的改变对土体强度变化的影响[13]。

2.2 非连续介质分析方法

以PFC为代表的颗粒离散元法的基本原理来源于分子动力学，是从微观结构角度研究介质力学特性和行为的工具。它的基本构成为圆盘和圆球颗粒。然后利用边界墙(WALL)约束。计算时不需要材料参数定义宏观本构关系和对应的参数，而是采用局部接触来反映宏观问题，因此只需要定义颗粒和黏结的几何和力学参数，从微观角度来说模拟结构性土具有很大优势。而在反映材料强度特性方面，PFC2D通过设置微观粒子的接触本构模型来模拟材料的宏观本构特性。接触本构模型分为：接触刚度模型；接触滑动模型；黏结模型。以黏结模型中的平行黏结模型为例，平行黏结模型是通过在颗粒间注入一定截面形状和尺寸的胶结材料的黏结，除了能传递力还能传递弯矩，这类似于结构性土在二元胶结模型[14]中的定义。因此，尝试以平行黏结模型作为接触本构模型来对结构性土进行模拟，微观参数的选取如表1所示。微观颗粒组合完毕之后，采用PFC2D内置的双轴压缩功能[15]来对材料的宏观特性进行检测，以验证PFC2D模拟结构性土的可行性。

表1 平行黏结微观参数设置Table 1 Parallel adhesion micro parameter setting

图2 双轴压缩模拟结果Fig.2 Biaxial compression simulation results

在围压恒定为100 kPa的条件下，双轴压缩的模拟结果如图2所示。材料在宏观上有明显的应变软化现象，这是因为随着应变的不断积累，微观颗粒之间的胶结物质被不断破坏，造成土体强度迅速下降。这与结构性土的强度变化特征[16-17]是相符的，在一定程度上验证了离散元软件模拟结构性土的可行性。

2.3 基于两类数值模拟方法的算例试验

分别就有限元法和离散元法在模拟结构性土坡失稳破坏的准确性方面进行对比。为更加直观地反映两类分析方法的区别，以某处真实花岗岩残积土边坡为例进行边坡失稳分析。

2.3.1算例说明

图3所示，为边坡简化分析模型。从A点所在的坡脚算起，边坡垂直高度为20 m，坡角约为60°，边坡土体的土性参数如表2所示。

图3 算例简化模型Fig.3 Simplified model of the landslide

天然容重/γunsat(kN·m-3)渗透系数Kx=Ky(m/ d)刚度E/kPa泊松比v黏聚力C/kPa内摩擦角φ/(°)200.0012 0000.352925

2.3.2试验方法说明

离散元方法：通过PFC2D颗粒离散元的可视化技术，可实现对边坡破坏过程的全程监测。首先，以实际算例为基础生成边坡的二维离散元模型，每隔一定距离给颗粒上色以便于观察。颗粒的微观接触模型选用平行黏结模型以反映边坡土体的结构性；然后，通过重力增加法使边坡破坏直至表现出明显的破坏特点；最后，与有限元方法的数值模拟结果进行分析比对。

有限元方法：选用有限元软件PLAXIS2D进行结构性土坡失稳特点的分析。在土体本构模型的选择上，以谢定义等[18]提出的综合结构势为基础，对传统莫尔-库伦模型进行改进，实现以结构性参数的大小来表征土体的结构性强度，进而通过多次折减土体结构性参数来使得边坡刚好发生极限破坏。最后，对不同折减系数下边坡渐进破坏的情况进行监测，生成剪应变增量云图及塑性点云图。除此以外，为突出不同工况对边坡稳定性的影响，对算例边坡进行两组有限元分析，一组单纯折减边坡土体的结构性参数，另一组在此基础上附加坡顶荷载。

2.3.3试验结果及分析

三组边坡的试验结果如图4～图6所示，结合其破坏特点的不同可分别称为有限元推移式滑坡、有限元牵引式滑坡以及离散元边坡失稳破坏。从整体破坏特征来看，两组有限元分析方法模拟出的结构性土坡破坏均为产生完整滑面的渐进破坏；而离散元分析方法模拟的结构性土坡失稳则观察不到完整的滑面，而是一种介于崩塌和滑坡之间，同时具备多种破坏特点的失稳破坏。

对比图4和图5可以看出：在不同工况影响下有限元边坡的破坏情况是不同的。在单纯折减边坡土体的结构性参数时，有限元边坡首先从坡脚产生塑性区域，而后塑性区不断向顶部发展并贯通，最终形成完整滑面；而在坡顶施加竖向荷载以后，有限元边坡的失稳变成了从坡顶首先产生塑性区域，而后再向坡脚处发展。对比图4c及图5c的土体拉伸点(黑点)及土体塑性点(白点)分布情况也可以看出两组有限元边坡破坏特点的区别。但总的来说，通过有限元方法模拟出的结构土坡破坏在不同工况下仅仅表现出了典型推移式滑坡或牵引式滑坡的特点，并没有反映出以黄土边坡、花岗岩残积土边坡为代表的部分结构性土坡突发、快速的破坏特征。相比之下，离散元分析方法模拟出的边坡破坏在失稳初期无明显位移变形，坡体内部产生了巨大的张拉裂隙，边坡失稳几乎与裂隙出现同时发生；到失稳中期，随着边坡失稳继续发展，坡顶裂隙不断发育并伴随着土体崩塌滑落，坡脚也产生了后缘张拉裂隙并发生局部崩塌；到失稳后期，边坡基本完全破坏，整个破坏过程更加类似于逐级倾覆式的岩坡崩塌，有明显的脆性破坏特征，并没有形成完整的滑面，较好的表现出了部分结构性土坡破坏时区别于一般土质滑坡的特点。

图4 有限元牵引式滑坡Fig.4 Finite element traction landslide

图5 有限元推移式滑坡Fig.5 Finite element push landslide

图6 离散元边坡失稳破坏Fig.6 Discrete element slope damage

3 基于两类数值模拟方法的结果对比

在反映结构性土强度变化本质方面：以有限元法为代表的连续介质分析方法受限于宏观本构模型的搭建，难以反映结构性土复杂的强度变化特性。同时，有限元法无法从微观角度对结构性土的强度变化加以分析；相比之下，以PFC2D为代表的颗粒离散元方法能从微观角度对模拟土体的特性进行定义，接触本构模型中的平行黏结模型被证明能较好的表现结构性土的应变软化现象。

在表现结构性土坡的破坏特征方面：有限元法模拟出的边坡破坏效果依旧是以形成完整滑面为主的渐进破坏，边坡失稳前位移变形较大，发展速度缓慢；而以PFC2D为代表的颗粒流软件则能模拟出类似于岩坡脆性破坏的特点。边坡破坏过程中内部裂隙发育且不断伴随着崩塌，更符合诸如花岗岩残积土边坡或黄土边坡破坏所具有的突发、快速的破坏特点，因此更加贴合实际。

4 结论

(1)以PFC2D为代表的颗粒离散元软件更适合分析结构性土坡稳定性方面的问题。花岗岩残疾土边坡及黄土边坡的失稳破坏与一般素土边坡区别较大，常表现出突发、快速、灾前位移变形较小的失稳特点。考虑到这两类边坡都属于土体结构性比较强的边坡，在对其进行稳定性分析时，土体微观结构的破坏对边坡土体强度的影响必须被考虑进去。通过目前常见的数值模拟方法的比较可以看出，以PFC2D为代表的颗粒离散元软件更适合分析结构性土坡稳定性方面的问题。颗粒离散元软件不仅能合理反映结构性土复杂的强度变化特点，同时也能较好的模拟出部分边坡存在的突发、快速的破坏模式。

(2)在利用颗粒离散元软件对结构性土坡的失稳问题进行分析时，合理的反映边坡土体的结构性是分析问题的关键。在这一过程中，应该采用数值模拟与真实试验相结合的方法，除了利用PFC2D自带的双轴压缩功能对微观颗粒材料进行宏观特性测试以外，还应该就真实材料进行室内双轴或者三轴压缩实验。数值模拟的结果应该以真实试验的结果为参照，不断修正才能合理反映结构性土的强度变化特性。