土石混合体边坡细观特征对滑面形成影响研究*

2020-09-10杨俊雄

工程地质学报 2020年4期

陈晓石崇杨俊雄

（①河海大学，岩石力学与堤坝工程教育部重点实验室，南京210098，中国）

（②河海大学岩土工程科学研究所，南京210098，中国）

0 引言

土石混合体是由粒径不同、强度较高的岩块和软弱的土体颗粒组成，一般存在于第四系松散堆积层中（Xu et al.，2009，2011；徐文杰等，2013）。其各种力学性质及特性和普通匀质土体差异较大，因此，由此类土体构成的岩土混合体边坡的稳定性影响因素较多，与普通匀质边坡的滑面发展和破坏机制有所不同（杨忠平等，2017；Wei et al.，2018；朱晨光等，2019）。这种类型边坡是大型土木工程的潜在威胁，因此研究该类介质滑动机理与稳定性具有重要意义。

土石混合体作为一种广泛分布在工程领域的重要地质材料，其力学性质及物理特性直接影响工程项目的安全稳定。Shakoor et al.（1990）研究了块石含量及尺寸对土石混合体抗压特性的影响。Wang et al.（2018）借助CT技术对三轴压缩试验中土石混合试样内部的微结构变化进行了观测。徐文杰等（2007a，2012）利用数值图像技术充分研究了土石混合体细观结构，将细观结构定量化，为土石混合体离散元数值模型建立提供基础；并进一步采用离散元及有限元方法模拟了土石混合体直剪试验，发现块石的存在使内部细观应力状态有结构性特征，含石量及试样尺寸均影响土石混合体的宏观抗剪强度和剪胀性（徐文杰等，2007a，2007b，2014）。李世海等（2004）等建立三维离散元随机计算模型，分析了土石混合体内部应力场分布与土石级配及岩石块度的关系。王环玲等（2019）通过颗粒流方法模拟了双轴压缩试验，研究了含石率、土体孔隙率、不均匀系数对土石混合体强度特性的影响。

土石混合体边坡广泛分布在中国西南地区，其稳定性对于水利水电等大型工程的建设尤为重要。由于土石混合体具有复杂的力学性质，Vessia et al.（2017）提出采用多维度理论模型计算土石边坡的安全系数；Medley et al.（2004）通过分析简单的土石混合地质模型，提出块石含量是影响边坡安全系数的重要因素；Napoli et al.（2018）建立了不同有限元边坡模型，并结合极限平衡方法，研究了不同块石含量边坡的稳定性。在国内的相关研究中，徐文杰（2009），徐文杰等（2008）以数字图像处理技术建立土石混合体边坡细观结构模型，并采用有限元强度折减法分析了块石对于边坡稳定性的影响；李亮等（2017）运用Flac3D模型，探讨了土石混合体边坡内部的塑性区扩展规律及土石界面强度参数对边坡安全系数的影响。刘顺青等（2019）提出了一种考虑不同含石率及块石随机分布的边坡稳定性分析方法。尽管目前在土石混合体边坡方面已有较多研究，但相关的研究更多集中于有限元方法（邵帅等，2014；张森等，2016；龚健等，2017），采用颗粒离散元的相关研究相对欠缺；且在以上研究中，均忽略了块石形态的影响，将块石形态简化为矩形或圆形。

颗粒离散元方法以颗粒（ball）和接触（contact）作为基本单元，可以直观地观测到边坡内部裂隙扩展及边坡失稳过程（程升等，2018；郑博宁等，2019）。基于此，本文以某一典型土石混合体边坡为例，建立土石混合体边坡颗粒流细观模型，依据室内试验数据标定土体材料参数。在合理细观模型的基础上，首先模拟均质边坡与土石混合体边坡的失稳过程，然后分析不同含石量边坡滑面发展过程的差异，探讨细观结构对滑面发展及边坡稳定性的影响。

1 土石混合体细观模型构造

1.1 土石混合体模型构建方法

土石混合体细观特征的构建，重在骨架颗粒的外轮廓和微观裂隙统计，在现场多采用统计窗方式。预先规划好地质统计窗，利用参照物拍照，然后借助AUTOCAD等工程软件绘制轮廓，进一步分析其细观特征。

数字图像处理方法是一种被广泛采用的土石混合体细观特征提取方法，但是数字图像方法的像素往往较高，如果将每一个像素都按照相应位置转化为数值模型，则单元、节点多，计算工作量大。

实际上，图像识别是一种有损识别方法，由于光照、阴影、拍照角度的差异，每一幅图片中土石区分均有差异，因此可以只采用数字图像识别的骨架颗粒轮廓线构造土石混合体的细观特征。故本文提出如下块石轮廓构造方法：

如图1a，将数字图像识别的骨架颗粒轮廓线作为边界线，对每一条边界进行读取（polyline），处于任一多段线内的像素属于块石、而不在任一多段线内的像素属于胶结物。这样即可将每一像素的性质（土或石）区分开，并借助这些多段线数据开展颗粒粒径、形状等信息的统计。最终，将这些信息导入PFC中作为模型生成的边界（图1b）。

图1 模型边界（徐文杰等，2008）Fig.1 Model boundaries（Xu et al.，2008）

1.2 土石混合体边坡模型的建立

首先，按照图1b所示的边界生成墙体（wall），再在边界内采用基本的颗粒生成命令形成完整的颗粒（ball）体系。生成的初始模型中颗粒是松散的，颗粒间的间隙较大，且应力分布极不均匀。此时，我们需要利用伺服机制（Cundall et al.，2008；Itasca.Consulting Group Inc.，2008）使初始模型达到一种密实且应力分布均匀的状态。

对边界wall施加一定的速度以表示恒定约束力，墙体移动产生的作用于墙体上的力可以表示为：

式中：kn为平均接触刚度；ΔUn为单位计算步内位移增量；Nc为与wall接触颗粒总数目；u·为墙体的运动速度；Δt为单位计算步时长。作用在墙体上的应力的变化值可以采用以下公式计算：

式中：A为与墙体接触作用的面积，二维情况下，通常取墙体的长度；于是可得wall的速度与单位计算步内应力增量的关系：

式中：G为伺服系数，不同计算步下，与墙体接触的总数目不同，每一计算步更新得到不同的伺服系数。伺服的过程就是通过调整墙体的运动使墙体与模型之间的接触力不断趋向于目标应力值的大小，墙体的接触力与目标应力值的差值可以表示为：

式中：σm为当前计算时步监测得到的墙体接触应力；σr为伺服结果需要实现的最终应力。由当前计算步的应力差值，根据式（5）和式（6）可以计算得到下一计算时步的墙体运动速度，不断迭代更新，直至Δσ′小于某一限定值。图2a为伺服前的模型图，图2b为伺服结束时的模型图。可以看出伺服后，模型的接触（contact）数目及分布均匀程度均更接近真实状态。

图2 接触分布形态Fig.2 The distribution of contact

在得到稳定的初始模型后，根据图1b所示的多边形边界，获取所有颗粒的位置信息，判断颗粒位置与多边形边界的关系，当判断位于多边形内时，则确定该颗粒属于岩石材料，位于同一多边形区域内的颗粒封装形成块石簇（clump）；并将边界内的其他颗粒定义为基质土材料。

通过以上模型构造方法，建立了如图3所示的颗粒离散元边坡数值模型。模型尺寸如图3a所示。模型共包含222i633个颗粒（ball），颗粒（ball）半径为0.5～1 cm。考虑到后续模拟滑坡过程的需要，在模型的左边界和下边界设置了边界条件，固定边界颗粒的运动，如图3中红色区域所示。

1.3 模型细观参数标定

图3 边坡细观结构模型Fig.3 Numerical model of slope mesoscopic structure

在颗粒离散元方法中，采用圆球（ball）模拟材料颗粒，并通过接触（contact）模拟颗粒之间的相互作用，通过为接触（contact）赋予不同的细观参数模拟材料的不同力学性质。故细观参数的标定是颗粒离散元数值模拟中的重要环节，合理细观参数的选取是数值模拟结果可信的基础（韩振华等，2019）。

在颗粒离散元中存在不同的接触本构模型（Itasca Consulting Group Inc.，2008），在本文中分别采用线性平行黏结模型（Linear Parallel Bond Model）和线性接触黏结模型（Linear Contact Bond Model）模拟岩石和土的力学性质。由于在滑坡过程中，块石不会发生破坏，故岩石的参数根据经验取值。对土的细观参数进行了标定。

在接触细观参数中，对宏观力学行为影响显著的主要包括黏结强度和有效模量。经相关研究发现（Potyondy，2004；Cho et al.，2007），接触有效模量直接影响材料的宏观杨氏模量；法向切向刚度比控制材料的泊松比；黏结的拉伸与剪切强度比与试样的破坏模式相关联。基于以上规律，并根据土体的室内试验曲线，如图4中虚线所示，首先确定强度参数和模量参数，保证不同围压下，试样的抗压强度与试验曲线基本一致。而后控制拉伸与剪切强度比保证试样表现出较强的塑性破坏特征。在标定过程中发现，当接触抗拉强度与剪切强度基本相等时，试样表现为明显的脆性破坏；随着抗拉强度与抗剪强度比值的增加，试样的塑性破坏特征增强；并且在低强度试样中，有效模量参数不仅影响试样的宏观杨氏模量，同时会对试样的强度产生影响。最终，标定得到的细观参数如表1所示，不同围压下的数值试验曲线如图4中实线部分所示。

除以上标定的强度特性参数外，在模拟滑坡过程中，阻尼参数的选取也尤为重要。在颗粒离散元中，共包含两种阻尼参数：局部阻尼和黏性阻尼。前者可以加快数值模拟的计算平衡，后者则反映颗粒碰撞过程中的能量消散。为了反映滑坡的真实过程，在模拟滑坡过程中不考虑局部阻尼作用，法向和切向黏性阻尼比分别取0.4和0.2。

图4 不同围压下土体应力－应变曲线Fig.4 Stress-strain curve of soil under different confining pressure

表1 土石混合体主要细观力学参数Table 1 The main mesomechanical parameters of soil-rock mixtures

1.4 边坡稳定性计算方法

强度折减法是由Dawson et al.（1999），Griffiths et al.（1999）提出的“包围方法”自动寻找安全系数。将宏观强度参数，如黏聚力和摩擦角，除以折减系数K作为新的材料参数，计算不同的折减参数并判断系统是否处于平衡状态，监测稳定与不稳定的临界状态，从而得到模型破坏时的安全系数。折减后ck、φk为：

在颗粒离散元中，对岩土体赋予的不同接触模型的细观力学参数并非宏观参数，因此需对该微观参数进行折减。其主要实现方法是通过对细观张拉强度、剪切强度以及摩擦系数进行参数折减，当模型的不平衡力比（ratio）能够达到1e－5时，则认为模型能够达到稳定状态，不平衡力比不能达到要求时，则认为模型为非稳定状态。具体折减参数如下：

2 模型计算结果

为了研究土石混合体细观结构对边坡滑面形成机制的影响，分别模拟纯土体边坡和土石混合体边坡的失稳过程，计算两种边坡的安全系数，记录边坡内部应力变化与颗粒间的传力机制，以及裂隙和滑面发展演变过程。对比两种边坡破坏机制的差异，分析土石混合介质对边坡滑面形成的影响。

2.1 初始边坡模型

分别将细观参数赋值于边坡细观模型，并在自重作用下平衡。图5是纯土体边坡自重平衡下的接触力链分布图。可以看出在重力作用下，坡体表层的力链分布稀疏，且坡体表层的接触力明显要小于坡体内部的接触力，坡脚位置力链密集区与稀疏区有较为明显的分界线，可以判定这是边坡的潜在滑面位置。

图6是土石混合体边坡在自重作用下平衡后的接触力链分布图，从图中可以看出边坡内部接触力链的分布情况与图5中纯土体边坡的力链分布明显不同，图6中力链的分布明显比纯土质边坡的力链分布复杂，尤其是在块石周围存在明显的剪切环，由于块石的存在，边坡内部接触力链的分布在遇到石块时会绕开石块，形成沿块体边缘的剪切闭环传力路径。且浅层的力链分布明显比纯土质边坡的接触力链分布密集。

2.2 边坡模型稳定性计算

图5 纯土体边坡接触力链分布Fig.5 The distribution of force chain of pure soil slope

图6 土石混合体边坡接触力链分布Fig.6 The distribution of force chain of soil-rock slope

为了准确描述该土石混合体边坡的稳定性，在边坡模型的浅层坡体内选取了5个监测点（图7），分别监测并记录不同强度折减系数下5个测点的位移值，进而求得土石混合体边坡的安全系数。监测结果如图8所示，其中图8a～图8d分别对应强度折减系数为1.0、1.05、1.10、1.15 4个工况下的监测点位移变化情况。从位移曲线中可以看出，强度折减系数为1.0时，5个监测点的位移曲线均在某一个水平上下浮动，且位移值均小于1 cm，故边坡处于稳定状态；强度折减系数为1.05时，测点1的位移明显较大，且其余测点位移也开始有上升的趋势，但最终位移均趋近于稳定，即边坡仍处于稳定状态；当强度折减系数达到1.10时，位移曲线形状发生突变，不再收敛于某一固定值。测点的位移均随计算时间增加而增加，且测点1的位移始终是最大的，故边坡失稳，发生滑坡；强度折减系数为1.15时，边坡的滑动更加明显，表层完全错动，即形成滑坡。从整个过程中可以得知：此边坡的安全系数在1.05～1.10之间。由于坡脚位移最大且最早发生，因此可推断此滑坡属于牵引式滑坡。采用同样的方法，经过计算发现，此纯土边坡的安全系数约为1.05。故此土石混合体边坡的安全系数要略高于相同条件下纯土体边坡的安全系数。

图7 边坡位移监测点Fig.7 Monitoring point of slope displacement

图8 不同强度折减系数下监测点位移Fig.8 Displacement of monitoring points under different strength reduction coefficient

2.3 边坡模型滑面特征

从以上计算结果可以看出，土石混合体边坡模型与纯土体边坡模型存在明显差异，为了分析土石混合体细观特征对滑面形成机制的影响，分别模拟纯土体边坡和土石混合体边坡的滑坡过程。在数值模型内部，微观裂隙的衍生过程即为边坡模型的失稳破坏过程。通过追踪模型内部的裂隙扩展路径，即可表征边坡内的滑面形成机制。图9、图10分别记录了纯土体边坡模型和土石混合体边坡模型计算25万步（图9a）、50万步（图9b，图10b）、75万步（图9c，图10c）、100万步（图9d，图10d）时的裂隙发育形态。

图9 纯土体边坡失稳过程Fig.9 The landslide process of pure soil slope

图10 土石混合体边坡失稳过程Fig.10 The landslide process of soil-rock mixture slope

从图9中可以看出，边坡表层的颗粒由于无墙体约束，在强度折减后黏结力降低，发生浅层的颗粒滑动。计算到50万步时边坡体内开始发生破坏，随着计算步数的增加，边坡内部滑面从底部尖端向上发展。与图10相比，纯土体边坡滑面的形成过程与土石混合体边坡模型滑面的形成过程存在显著不同。纯土体边坡滑面发展方向单一，微观破坏裂隙的分布沿滑面的发展方向相对分散。相较于纯土体边坡内微观破坏的离散性，土石混合体边坡内的微观破坏局部集中，且由于块石的存在，坡体内部的破坏形态明显发生改变，滑面的位置更趋近于表层，在破坏速度上，也明显低于纯土体边坡，相同计算时间内，纯土体边坡模型内的微裂隙数目约为土石混合体边坡模型微观裂隙数目的1.5倍，故块石的存在明显改变了边坡内部的细观结构，增强了坡体抵抗破坏的能力，提高了边坡的稳定性。同时，从图10b和图10c中可以发现，破坏滑面明显沿着两个方向发展，一条主滑面，一条次滑面，两条滑面均穿过块石之间的间隙，且裂隙的分布在其发展方向上相对集中。图10d中边坡体内发育形成一条贯通的可视滑面，最终的滑面位置存在于块石之间的间隙。

图11为与图1所示边坡相似的一处土石混合体边坡，地质情况相似，滑坡的形状也接近，此图片为边坡滑坡后的情形，从图中可以看出滑面的位置与滑坡体形态均与上述模拟结果基本吻合。

图11 现场滑面Fig.11 The sliding surface on site

2.4 不同含石率边坡滑面破坏特征

图12 不同含石率下边坡滑面扩展过程Fig.12 The spreading process of slope sliding surface under different rock content

为了进一步论证土石混合体细观特征对边坡破坏机理的影响，分别进行了土石混合体边坡中含石量为10%、20%、30%、50%的失稳过程模拟（图12）。含石率为10%时，滑坡体表现为典型的圆弧状破坏形式，滑面表现为自坡脚向坡顶贯穿的圆弧形态。含石率为20%时，在坡体底部，微裂隙沿块石间的薄弱位置发育。由于在坡体顶部，块石分布密集，滑面形成受阻，在块石周围散布少许微裂隙，且未完全贯通。故在边坡体内形成一条不完整的圆弧状滑动带。含石率为30%时，边坡破坏形式发生明显改变，裂隙在坡脚位置集中，且基本分布在边坡浅层附近，边坡的破坏形式由整体滑动变为坡脚局部滑动，在坡脚的表层形成一条局部的滑面。含石率为50%时，微裂隙的分布基本密布在坡脚位置，在坡脚的位置形成一条局部的圆弧滑面，其滑面半径约是含石率10%时的圆弧状滑面半径的1／4。就不同含石率下的滑坡体体积而言，含石率10%与20%时，滑坡体体积基本相同，无较大差异。含石率30%时，滑坡体体积明显减少，约为含石率10%时滑坡体的1／8。含石率继续增加至50%时，滑坡体体积继续减少，约为含石率30%时滑坡体体积的1／2。

故土石混合体细观特征对于边坡稳定性及滑面形成过程影响显著。不同含石率下，滑面的形态各异，尤其含石率超过30%时，边坡的失稳形式发生明显变化，滑面的形态也由整体向局部变化。