基于随机裂隙的长距离调水工程渠道边坡稳定性模拟

王栋 方攀博 韩延成 黄兆虎 邓晓川 赵嘉诚

关键词：裂隙岩体；随机裂隙网络；流固耦合；强度折减法；边坡稳定性

很多长距离调水工程如南水北调中线和东线、引黄济青工程、胶东调水工程等采用明渠形式输水。明渠经过裂隙发育地区时，裂隙和孔隙会对明渠边坡内部的渗流场、位移场产生重要的影响，导致渠道产生渗透破坏，使边坡稳定性降低，严重影响明渠的输水能力和安全。

引起边坡失稳的原因有很多，主要原因有降雨入渗、地下水位变化、边坡地质构造等。降雨因素特别是历时长、强度大的暴雨入渗到边坡内部，影响土体骨架的力学性质，使基质吸力下降，抗剪强度降低，从而引起边坡失稳［１］。地下水位的变化会影响土体强度，使坡面产生不稳定变化［２］。马慧敏等［３］研究了工程特性较特殊的膨胀土（岩）渠段的边坡稳定性，并给出了分析边坡稳定性时参数取值的经验方法。很多调水工程边坡存在裂隙且现有的边坡稳定性分析软件不具备模拟岩体中随机裂隙网络的功能，使得模拟结果并不能完全反映实际情况下渠道的边坡稳定性情况，故模拟含裂隙边坡的稳定性较为困难。常峻岭［４］根据蒙特卡罗（Ｍｏｎｔｅ－Ｃａｒｌｏ）法模拟随机裂隙网络，研究了几何参数、力学参数对边坡稳定性的影响。邓东平等［５］研究了张裂缝位置及张裂缝水位对边坡稳定性的影响。目前边坡稳定性分析方法主要有极限平衡法和有限元强度折减法，因极限平衡法并不能求出变形情况，且不能反映实际条件下边坡所受内力，还需要假设滑动面位置及形状，故较多学者采用有限元强度折减法进行边坡稳定性分析［６－９］。孔不凡等［１０］和董金玉等［１１］将离散单元法和强度折减法结合模拟边坡稳定性，同样取得了较好的成果。

近年来因岩体裂隙而造成渠道边坡破坏的事故较多，裂隙渗透对渠道边坡稳定性的影响受到广泛关注。笔者以山东省胶东调水工程莱州段输水明渠为例，基于蒙特卡洛原理，生成符合实际岩体裂隙分布的二维随机裂隙网络，将达西定律模块和固体力学模块进行耦合，建立了渠道边坡二维随机裂隙模型，采用强度折减法对比研究了两种降雨条件（历时长而强度弱、历时短而强度强）和改变原始地下水位条件下两种介质模型（裂隙模型和孔隙模型）渠道的边坡稳定性。

１研究区概况

山东省胶东调水工程莱州段输水明渠位于莱州市虎头崖镇西北部的秀东村至后趴埠村。渠道边坡岩性主要为花岗岩，受岩体边坡中裂隙影响，强降雨期间渠道的地下水位变幅较大。

２随机裂隙网络生成

通过水文地质勘测资料及渠道两岸钻孔调查，对裂隙各参数分布规律进行计算分析，将研究区的裂隙分为两组，得到各组裂隙的初始几何参数，见表１。基于蒙特卡洛基本原理［１２］，利用Ｍａｔｌａｂ编制二维离散裂隙网络生成程序，生成的裂隙网络模型（如图１所示）尺寸为６０ｍ×３０ｍ（长×深）。山东省胶东调水工程莱州段（桩号６２＋８０９）输水明渠边坡单侧断面如图２所示。处理后最终生成随机裂隙分布图（见图３）。

３模型构建及原理

３．１渗流－应力耦合模型

３．１．１渗流微分方程

由于降雨对渗流及边坡稳定性具有很大的影响，因此渗流计算采用非饱和渗流Ｒｉｃｈａｒｄｓ方程［１３］：

３．１．２渗流－应力耦合原理

１）渗流场对应力场的影响。岩体中渗透压力对岩体的变形有重要影响，动水压力导致裂隙开度增大，岩体变形影响岩体的稳定性［１５］。Ｌｏｕｉｓ［１６］通过分析测试认为渗透压力对裂隙岩体的作用由３种体积力构成，即在一组平行节理中由水的黏性产生的切向力、静水压力或阿基米德压力、动水压力或渗透力。

２）应力场对渗流场的影响。当岩体受渗透体积力和其他外部荷载作用时，岩体中应力场发生改变，从而岩体裂隙宽度改变，进而影响渗透系数和岩体渗流场。根据达西定律，流体的渗透系数可表示为［１５］

可知，影响岩土体渗透性能的因素主要为流体性质和岩土体骨架性质［１７］。

３．２边界条件及相关参数

在渗流模型中，上边界和坡面上层边界为大气边界（在降雨条件下为降雨边界）；坡面下层和坡底位置即渠道衬砌部分实测渗漏量非常小，可以忽略不计；渠道内无水时渠道边坡内地下水水位与渠道内水位之差最大，扬压力最大，为最不利工况，因此渠道底部概化为无流动边界；模型左侧边界为水头边界即初始地下水位，右侧按年庚乾等［１８］研究的渗流边界处理，渗流边界方程如下：

在力学模型中，坡面上层边界、渠道坡面和渠底边界均为自由边界；模型左侧边界和右侧边界均设置为辊支撑边界（横向位移为０）；坡体底部边界为固定边界。对于渗流边界，当地下水位上升或者土体达到饱和状态时，水从边坡内不断渗出。边界设置如图４所示（黄色为渗流边界条件、绿色为力学边界条件）。

在渗流计算模块中饱和－非饱和渗流计算采用达西定律匹配理查兹方程作为基本方程，裂隙模型和孔隙模型除了饱和渗透系数不同外，其他采用相同的水力参数。非饱和水力参数采用ＶａｎＧｅｎｕｃｈｔｅｎ模型计算。固体力学模块计算模型采用线弹性本构模型搭配土壤塑性模型，采用Ｄｒｕｃｋｅｒ－Ｐｒａｇｅｒ准则匹配Ｍｏｈｒ－Ｃｏｕｌｏｍｂ准则作为破坏准则。相关计算参数由地勘结果、监测井抽水试验、室内力学试验得到，见表２、表３。

３．３强度折减法失稳判据

目前，用强度折减法来判断岩土边坡的破坏失稳判据主要有３种：１）有限元计算结果不收敛，说明此时边坡已经发生破坏，相应的应力应变无法确定，则前一个输入的折减系数对应的状态为此边坡的极限平衡状态［１９］；２）如果塑性变形区域贯通了整个土层，则表示该区域即将失稳，在塑性变形贯穿前的强度折减系数即为边坡的安全系数［２０］；３）边坡位移改变率发生突变時折减系数就是边坡的安全系数［２１］。

４结果分析

４．１模拟方案

共设置１２种模拟工况，见表４。

４．２渗流分析

本文用Ｃｏｍｓｏｌ有限元数值计算软件分别进行基于孔隙模型和裂隙模型的非饱和渗流模拟，得到任意时刻的渗流场。图５和图６为降雨强度为２０ｍｍ／ｈ和５ｍｍ／ｈ条件下７０ｈ和２８０ｈ的渗流场（云图代表孔隙水压力，面上流线代表达西速度场）。

渗流模拟结果表明，裂隙模型的入渗深度大于孔隙模型的，裂隙模型的体积含水率小于孔隙模型的。这是因为雨水进入裂隙后迅速下渗，入渗深度迅速增大，裂隙内的入渗水量下渗至地下水，所以入渗区内体积含水率较小，孔隙水压力减小。裂隙模型的入渗率大于孔隙模型的，随着降雨时长的增长，两种介质模型的入渗率分布趋近于稳定值。

４．３不同降雨条件边坡稳定性分析

根据上述渗流模拟结果，进行不同介质模型边坡塑性变形和位移数值模拟结果对比。

４．３．１塑性变形

通过数值模拟得到不同降雨条件边坡塑性变形的情况。可知，在２０ｍｍ／ｈ的降雨强度下，裂隙模型最大折减系数为１．４８，最大塑性变形为３．０１ｃｍ；孔隙模型最大折减系数为１．５０，最大塑性变形为１．７１ｃｍ。５ｍｍ／ｈ降雨强度下，裂隙模型最大折减系数为１．５０，最大塑性变形为１．６２ｃｍ；孔隙模型最大折减系数为１．５２，最大塑性变形为１．２４ｃｍ。由此可得：随着折减系数的增大，塑性变形区域从坡脚逐渐向坡顶移动，折减系数最大时从坡脚至坡顶形成了贯通区域，最大塑性应变发生在坡脚处。孔隙模型的安全系数大于裂隙模型的，且塑性变形小于裂隙模型的。

４．３．２位移

模拟得到在极限平衡状态不同降雨强度下裂隙和孔隙模型的位移分布云图，如图７、图８所示。由图７、图８可知，裂隙模型和孔隙模型得到的位移场大致相同，边坡整体的滑动趋势都是从坡顶逐渐向坡脚移动且坡脚处位移最大。２０ｍｍ／ｈ降雨强度下，裂隙模型坡脚处最大位移为７．４３ｃｍ，安全系数为１．４８；孔隙模型坡脚处最大位移为３．２３ｃｍ，安全系数为１．５０。５ｍｍ／ｈ降雨强度下，裂隙模型坡脚处最大位移为３ｃｍ，安全系数为１．５０；孔隙模型坡脚处最大位移为２．２９ｃｍ，安全系数为１．５２。纵向对比可知，雨强较小时边坡模型的安全系数较大，边坡更不易发生失稳；横向对比可得，相同降雨条件下，裂隙模型的边坡安全系数始终小于孔隙模型的，孔隙模型边坡更稳定。

４．３．３最大位移与安全系数的关系

基于孔隙模型和裂隙模型模拟得到不同降雨条件下边坡最大位移与安全系数的关系，如图９所示。由图９可知，边坡最大位移随安全系数的增大而增大，且增长趋势大致相同。当安全系数小于１．４６时，降雨强度和模型介质是否存在裂隙对边坡的最大位移影响很小，两种介质模型的最大位移几乎相同。安全系数为１．４６～１．４８时，最大位移增长率从大到小对应的降雨强度工况是裂隙模型２０ｍｍ／ｈ、裂隙模型５ｍｍ／ｈ、孔隙模型２０ｍｍ／ｈ、孔隙模型５ｍｍ／ｈ。安全系数为１．４８～１．５０时，５ｍｍ／ｈ降雨强度下，裂隙模型位移增长率大于孔隙模型的。综上，随着安全系数的增大，裂隙模型的位移增长率和最大位移均大于孔隙模型的，裂隙模型边坡较孔隙模型更易破坏且失稳时位移更大。

４．４不同地下水位边坡稳定性分析

４．４．１塑性变形

通过数值模拟得到不同地下水位下边坡塑性变形情况，数据见表５。由表５可知，两种介质模型边坡均形成了明显的贯通坡脚和坡顶的塑性变形区域，且裂隙模型塑性变形区域大于孔隙模型的。最大塑性变形发生在高地下水位条件下，裂隙模型和孔隙模型对应的最大值分别为５．１４、４．４１ｃｍ。随着地下水位的降低，折减系数逐渐增大，坡脚处的最大塑性变形区域减小。裂隙模型的边坡塑性变形均大于孔隙模型的。

４．４．２位移

模拟得到在极限平衡状态下不同地下水位裂隙和孔隙模型的位移等值线分布云图。可知，位移较大的部分从坡顶一直贯通至坡脚处。随着地下水位的降低，两种介质模型的安全系数逐渐增大，边坡稳定性逐渐增强，低地下水位条件下模型的边坡稳定性最佳。在２０ｍｍ／ｈ的降雨强度下，裂隙模型和孔隙模型的边坡安全系数都为１．５２，最大位移分别为７．４３、３．２３ｃｍ，裂隙模型位移较大的部分主要集中在边坡中下部，失稳时有可能发生牵引式滑动，其原因是降雨强度较大，降雨入渗时，雨水几乎全部通过裂隙向下渗透，入渗速度较快，入渗深度较大，坡脚处先达到饱和，从而使坡脚处的稳定性大幅减弱。在５ｍｍ／ｈ的降雨强度下，裂隙模型和孔隙模型边坡安全系数都为１．５２，最大位移分别为１４、１１ｃｍ，是１２种工况中最大的。原因是降雨历时太长且地下水位较高时，入渗水较快到达地下水位后上升至坡脚，致使坡脚处一定范围内岩体浸水全部达到饱和，坡脚长时间浸于水中承受较大孔隙水压力而产生较大的位移。裂隙模型安全系数始终小于等于孔隙模型的，在相同安全系数的前提下，裂隙模型的边坡最大位移大于孔隙模型的，孔隙模型的边坡稳定性更佳。相同地下水位、雨强较小时，边坡安全系数更大、更稳定。

４．４．３最大位移与安全系数的关系

基于孔隙模型和裂隙模型模拟得到改变初始地下水位条件下边坡最大位移与安全系数的关系，如图１０、图１１所示。可知，在２０ｍｍ／ｈ降雨强度下，位移增长率最大的是工况１，边坡稳定性最好的是工况６，工况６安全系數最大且最大位移最小；在５ｍｍ／ｈ的降雨强度下，位移增长率最大的是工况１１，边坡稳定性最好的是工况４。发生最大位移的分别是工况１１和工况１２（原因已在４．４．２中解释）。边坡最大位移随着安全系数的增大而增大，相同地下水位条件下裂隙模型边坡最大位移均大于孔隙模型的。

５结论

本文基于裂隙模型和孔隙模型，采用强度折减法对比研究了不同降雨条件、不同地下水位条件下两种模型对胶东调水工程莱州段输水明渠渠道边坡稳定性的影响，得到如下结论：

１）裂隙模型的安全系数小于等于孔隙模型的，且在相同安全系数的前提下，裂隙模型边坡的塑性变形和最大位移大于孔隙模型的。裂隙模型在位移突变时的位移增长率大于孔隙模型的，表明裂隙模型的边坡稳定性更差。

２）两种模型模拟结果均表明，在相同降雨强度条件下，高地下水位时边坡稳定性较差；在相同地下水位条件下，降雨强度较大时边坡稳定性较差。

３）裂隙模型在高雨强条件下位移较大的部分主要集中在边坡中下部，失稳时可能发生牵引式滑动。

因此，对裂隙发育地区的渠道边坡稳定性分析，按常规忽略裂隙的孔隙模型计算，结果可能与实际不符，结果偏于安全，会给工程运行带来潜在的安全风险。