考虑土体参数空间变异性的边坡大变形破坏模式研究*

刘 鑫 王 宇 李典庆

( ①武汉大学水资源与水电工程科学国家重点实验室 武汉 430072)

( ②香港城市大学建筑学及土木工程学系 香港 999077)

0 引 言

岩土工程中不确定性不可避免，其中一个重要的不确定性是土体参数内在的空间变异性( Phoon et al.，1999; Wang et al.，2016b; 赵晶等，2018) 。由于边坡总是沿着土体中的最软弱区域滑动，土体参数的空间变异性与边坡失稳密切关联。大量研究已表明土体参数的空间变异性对边坡失效概率有较大的影响，这些研究大多采用随机极限平衡方法( Random Limit Equilibrium，RLEM) ( EI Ramly et al.，2002; Cho，2007; 杨 继 红 等，2007; Ji et al.，2012; 蒋水华等，2014; 杨智勇等，2019) 或随机有限元方法( Random Finite Element Method，RFEM)( Griffiths et al.，2004，2009; Hicks et al.，2010; 肖特等，2014; Li et al.，2016a) 。然而，由于极限平衡方法或有限元方法的限制，它们难以模拟土体的大变形破坏过程。实际上，边坡失稳是一个涉及启动、滑动和堆积3 个阶段的动态过程。但极限平衡方法仅满足静态平衡方程( 陈祖煜，2003) ，有限元方法在模拟大变形时可能遇到网格破坏问题导致计算不收敛。边坡失稳的全过程将决定边坡的失事后果，因而有必要研究考虑土体参数的空间变异性情形下边坡的大变形破坏过程及其特征。

近年来，物质点法( Material Point Method，MPM) ( Sulsky et al.，1994; 张雄等，2013) 被广泛应用于岩土领域的大变形问题中，如边坡渐进破坏( Zabala et al.，2011; Yerro et al.，2016) 、多相介质失稳( Bandara et al.，2015; Yerro et al.，2015; Wang et al.，2018) 、管涌( Yerro et al.，2017) 、静力触探试验( Ceccato et al.，2016) 等。物质点法与有限元法类似，是一种基于连续介质力学的数值模拟工具。物质点法将问题域同时划分为欧拉背景网格和拉格朗日物质点，背景网格用于求解动量方程，拉格朗日物质点携带材料和场变量( 如速度、加速度) 信息，使物质点可以自由变形，避免了有限元里重新划分网格带来的问题( Yerro et al.，2016) 。这使得物质点法适于分析边坡稳定性分析和大变形破坏全过程( 史卜涛等，2016; 王斌等，2017; 张巍等，2017) 。为了考虑土体参数的空间变异性，已有学者提出了结合直接蒙特卡洛模拟的随机物质点法( Random Material Point Method，RMPM) ( Wang et al.，2016a) ，但计算消耗较大，阻碍了它在实际工程中的应用。

为此，本文采用一种随机极限平衡－物质点法分析考虑土体参数的空间变异性情形下的边坡大变形破坏模式。该方法同时利用极限平衡方法简单、高效的优势和物质点方法模拟土体大变形破坏的能力，能够实现高效的边坡可靠度分析和风险评估。

1 随机极限平衡－物质点法

随机物质点法( Wang et al.，2016a) 一般采用直接蒙特卡洛模拟产生N 个随机样本进行不确定分析，边坡失效概率Pf的表达式为:

式中，Nf，RMPM为失效样本总数，即发生土体大变形破坏的随机样本数目; N 为随机样本总数，建议至少取N=10/Pf保证失效概率估计的变异系数小于0.3( Ang et al.，2007) 。

然而，边坡失稳通常是小概率事件，直接蒙特卡洛模拟需要产生大量的样本，考虑到每次物质点法分析本身计算耗时较长，直接蒙特卡洛模拟的计算总耗时可能过高。例如，对于边坡失效概率Pf=0.001，至少需要产生10 000 个随机样本。若采用随机物质点法完整分析所有样本，考虑每次分析耗时15 min，不考虑并行情况下计算总耗时将高达15万分钟，即2500 h。这对于工程实践是难以接受的。

事实上，蒙特卡洛模拟的N 个样本中仅有Nf，RMPM个失效样本，其他的未失效样本占大部分，对应在物质点法分析中没有发生大变形破坏。所以没有必要采用物质点法分析这些稳定的样本，这样可以大大提高计算效率。

本文提出的随机极限平衡－物质点法将结合随机极限平衡方法的计算效率优势和随机物质点方法模拟大变形破坏的能力，用于高效的边坡可靠度分析与风险评估( Liu et al.，2019) 。图1 绘制了本文方法的框架，包含4 个部分。第1 部分是采用传统极限平衡方法和物质点方法分别建立确定性分析的边坡模型，并保证两个模型的一致性，即拥有相同的几何形状、土体参数等。第2 部分采用随机场理论( Vanmarcke，2010) 作为不确定性模型模拟空间变异的土体，许多文献( 蒋水华等，2014; 肖特等，2014) 已有介绍，这里不再赘述。第3 部分是采用蒙特卡洛模拟进行不确定性分析。第4 部分为结果后处理，处理随机极限平衡方法和随机物质点方法得到的结果，如计算边坡失效概率，研究边坡大变形破坏模式等。

图1 随机极限平衡－物质点法框架Fig. 1 Illustration of the proposed framework

首先利用蒙特卡洛模拟产生N 个随机场样本，采用随机极限平衡方法进行不确定性分析，计算所有随机场样本的安全系数值。由于随机极限平衡方法计算简便、耗时短，可以轻易地完成N 个随机场样本的不确定性分析。接下来将所有随机场样本根据安全系数值从小到大排序，后续计算过程被划分为k 个迭代步。对于第k=1 次迭代，从排序后的样本序列中取第1 个随机场样本采用随机物质点分析，判断是否发生大变形破坏，若发生破坏则令失效样本总数目加1，采用式( 1) 计算失效概率。若不满足终止条件，则迭代步k=k+1。对于第k 次迭代( k＞1) ，从排序后的样本序列中取第k 个随机场样本采用随机物质点方法分析，判断是否发生大变形破坏，若发生破坏则更新失效样本数目，并采用式(1) 计算失效概率。迭代计算的终止条件定义为当失效概率收敛到一个确定的值，即迭代过程中计算的失效概率经过kt个迭代步保持几乎不变。根据经验，kt可取为1%～10%的样本总数目N。

2 边坡大变形破坏模式

为了验证本文提出的随机极限平衡－物质点法，采用一个双层土坡算例分析。图2 绘制了该土坡的几何形状。上层土高18 m，下层土高10 m，坡比为3︰4，水平距离被延长到112 m 用于堆积滑动体。已有研究采用随机有限元方法对该土坡进行可靠度分析( Li et al.，2016b) ，为了保持一致，土体参数取自该文献。两层土的重度均等于19 kN·m－3。两层土的不排水强度均服从对数正态分布。上层土不排水强度的均值为80 kPa，变异系数为0.30; 下层土不排水强度的均值为120 kPa，变异系数为0.30。采用指数型随机场模拟土体不排水强度的空间变异性，有关细节可参考文献( Li et al.，2015) 。水平相关长度为24 m，竖直相关长度为2.4 m。图2 绘制了随机场网格，网格大小为1 m×1 m，采用乔列斯基分解离散随机场。

图2 二层土坡几何形状Fig. 2 Geometry of a two-layer soil slope

对于边坡稳定分析，本文与前人一样假设滑动面为圆弧形并采用简化毕肖普法计算边坡安全系数。对于物质点方法分析，整个区域被离散为1 m×1 m 的网格，每个网格均匀分布4 个物质点。采用考虑线性软化的Drucker-Prager 模型作为本构模拟土体 的 弹 塑 性( Bandara et al.，2015; Wang et al.，2018) 。土体的残余强度取原始强度的50%，硬化模量为－50 kPa( 即软化) 。物质点法总模拟时长为20 s，每个时间步长约为0.000 75 s，重力加速度等于9.81 m·s－2。

首先采用直接蒙特卡洛模拟产生40 000 个随机场样本，采用文中的随机极限平衡－物质点方法分析每个随机场样本是否发生大变形破坏( 阈值为最大相对位移超过1 m，见Liu et al. ( 2019) ) 。最终得到一共有520 个样本发生破坏，对应随机物质点方法的失效概率等于Pf，RLEM=520/40000=1.30×10－2。该失效概率与随机极限平衡方法得到的1.12×10－2和随机有限元方法( Li et al.，2016b) 得到的1.06×10－2较为吻合。使用本文方法，无需采用随机物质点方法分析全部40 000 个随机场样本，实际仅耗费2512 个随机场样本即达到收敛，计算消耗仅约为原直接蒙特卡洛模拟的6%，大大提高了计算效率。

对于本算例，本文方法通过随机物质点方法分析可得到不同的破坏模式，根据滑动深度不同，可划分为浅层、深层、中层破坏模式，另外还有渐进破坏模式涉及到多重破坏模式同时或渐进发生。

2.1 浅层破坏模式

图3 边坡浅层破坏模式位移与不排水抗剪强度分布Fig. 3 Displacement and undrained shear strength distributions of a shallow slope failure mode

图3a、图3b 分别绘制了一种典型边坡浅层破坏模式的位移与不排水抗剪强度分布图。对于本算例，浅层破坏模式被定义为仅发生在上层土体。该破坏模式土体沿着上下层土的交界面滑动，顺着坡面下滑并堆积在坡脚处。图3 显示该破坏模式发生位移约为30 m。极限平衡方法得到的临界滑动面( 即对应安全系数最小的滑动面) 与物质点法中实际发生的滑动体形态较为一致。如图3b 所示，上层土靠近坡面部分的不排水抗剪强度较低，低于左侧土体，显著低于下层土。这种情况容易产生浅层破坏。

2.2 深层破坏模式

图4a 绘制了一种典型的深层破坏模式。深层破坏模式被定义为滑动带底部低于坡脚高度的破坏模式。该破坏模式滑动体积一般大于浅层破坏模式。该滑动体整体沿着某一个圆心转动发生失稳，滑动体形态保持较为完整，故被细分为径向深层破坏模式。图4b 绘制了对应的不排水强度分布。下层土底部及坡脚附近的不排水抗剪强度均较低，且下层土的软弱区与上层土的软弱区相连形成了弧形软弱带，可能促使径向深层滑动模式的发生。

图5a 绘制了另一种典型的深层破坏模式。与图4a 绘制的径向深层破坏模式不同，图5a 中滑动体主要沿着边坡底部水平切向滑动，转动较小。且滑动体被分为两个三角形的滑动体，两个滑动体之间存在明显的错动带。这可能是由于滑动体前后的滑动速度不一致造成滑动体沿着某个倾斜的软弱带发生错动。这种破坏模式被分类为切向深层破坏模式。该破坏模式的滑动距离约为24 m，大于图4a中的径向深层滑动模式。从图5b 中可以看出下层土的底部存在连续的水平软弱带，使得滑动体更易于沿着底部水平运动。值得注意的是，该破坏模式的形态不符合圆弧滑动面假设，说明圆弧滑动面的假设在考虑土体参数的空间变异性情形时可能不成立。

图4 边坡径向深层破坏模式位移与不排水抗剪强度分布Fig. 4 Displacement and undrained shear strength distributions of a radial deep slope failure mode

图5 边坡切向深层破坏模式位移与不排水抗剪强度分布Fig. 5 Displacement and undrained shear strength distributions of a tangential deep slope failure mode

2.3 中层破坏模式

图6 边坡中层破坏模式位移与不排水抗剪强度分布Fig. 6 Displacement and undrained shear strength distributions of an intermediate slope failure mode

中层破坏模式的滑动深度介于浅层和深层破坏模式之间，通常沿着坡脚切向滑动。图6a 绘制了一种典型的中层破坏模式的位移分布。与浅层破坏模式相比，中层破坏模式的滑动体积更大; 与深层破坏模式相比，滑动体的重力势能转换为动能后无需抵抗重力做功，故中层破坏模式的滑动距离通常大于浅层和深层破坏模式。图6b 绘制了该破坏模式的不排水抗剪强度分布。可以看出，坡脚高度之上的软弱区与上层土的软弱区相连，坡脚以下土体强度相对较高，从而导致了这种破坏模式。

2.4 渐进破坏模式

除此之外，本研究还发现了渐进破坏模式，即多种破坏模式同时或渐进发生。图7 绘制了一种典型的渐进破坏模式。图7a 和图7b 显示了时间为第4 s 和第8 s 时一个中层破坏模式在滑动，滑动形态与极限平衡方法得到的临界滑动面一致。图7c表明在第10 s 时一个新的浅层破坏模式逐渐形成，最终破坏形态见图7d，最终滑动距离达到40 m。由于渐进破坏模式涉及到多种破坏模式的发生，多个模式的滑动体积更大，滑动距离更远，造成的后果通常更为严重。需要指出的是，由于极限平衡方法或有限元方法难以模拟土体的大变形破坏全过程，它们难以得到后续的破坏模式，也就无法用于分析渐进破坏模式。说明了物质点法在模拟边坡大变形破坏全过程的优势。

图8 绘制了渐进破坏模式对应的不排水抗剪强度分布。它的土体不排水抗剪强度分布与中层破坏模式的类似，但上层土的强度较低，上层土在第1个中层破坏模式发生后失去支撑，导致了第2 个浅层破坏模式的发生。

图7 边坡渐进破坏模式位移演化Fig. 7 The displacement evolution of a progressive slope failure mode

3 结 论

本文采用随机极限平衡－物质点法分析了考虑土体参数空间变异性的边坡大变形破坏模式，得到了包括浅层、中层、深层和渐进一共4 种典型破坏模式，并分析了对应不排水抗剪强度与边坡破坏模式之间的关系。得到如下结论:

(1) 边坡的破坏模式与土体强度参数的空间分布密切相关。不同的土体强度参数空间分布可能引起不同的破坏模式，导致不同的失事后果。因此，在进行边坡风险评估与加固时，应重视勘察工作与勘察数据的利用，尽可能地量化土体参数的空间变异性。从而针对不同的可能破坏模式，做出更有效的加固措施。

( 2) 与极限平衡方法和有限元方法相比，物质点法具备土体大变形分析的能力，既可用于分析边坡是否失稳，又具备模拟边坡土体大变形破坏的全过程的能力。本研究提出的随机极限平衡－物质点法可以同时利用极限平衡方法的计算效率优势和物质点方法的土体大变形分析能力，可以实现高效地边坡风险评估。

( 3) 在考虑土体参数的空间变异性时，多数情况下极限平衡方法得到的圆弧滑动面与物质点法得到的滑动体形态较为一致，但在某些情况下圆弧滑动面假设可能失效( 如存在较大连通区域的软弱带时) 。此外，极限平衡方法难以用于考虑渐进破坏模式，即多个破坏模式同时或渐进发生情形。

( 4) 本文所采用的物质点方法仅分析单相土体，未来仍需拓展到多相用于考虑土－水的耦合作用，考虑孔隙水压力等更复杂工况。本文提出的随机极限平衡－物质点法仍然可用于边坡可靠度分析和风险评估。

( 5) 可以看出与传统有限元方法相比，物质点法适于边坡和其他岩土问题的大变形模拟和破坏分析，可以利用现有的土体本构模型，但物质点法在岩土领域的研究尚处于初期阶段，存在一些问题有待完善，如考虑土－水的多相耦合作用、隐式物质点法、应力波动和与有限元相比计算精度不高等问题。