某复杂古滑坡堆积体在蓄水条件下的稳定性评价

阚露 张发明

摘要：文章以位于某水电站坝址区右岸古滑坡Ⅱ区堆积体为研究对象，建立了一个较为复杂的三维模型，并从主应力、塑性区、剪应力和位移四个方面对古滑坡堆积体在蓄水条件下的变形和稳定性情况进行了详细的模拟和分析。计算得到的四个方面的结果都相互吻合相互印证。文章以FLAC3D中的单元界面能自动依附于指定范围内模型表面生成的特性为思路来模拟地下水位。研究结果表明，堆积体虽然未出现整体性滑动破坏，但要注意坡脚局部地区坍塌和土层内部局部滑弧滑动，具体直观地在模型上指出堆积体上的不稳定区域，为后期针对性的加固治理方案提供参考，为工程优化设计提供借鉴意义。

关键词：堆积体；三维数值模拟； 稳定性

中图分类号：文献标志码：A

文章编号：1672-1098（2013）03-0000-00

崩塌、滑坡发生后形成的堆积体变形及稳定性问题一直是工程地质领域研究的重点问题之一，滑坡形成的堆积体结构松散，强度低，在受到扰动后容易失稳，对水电站的水库正常运营会造成很大的影响，尤其是大型堆积体如果发生破坏，很可能会堵塞河道，后果不堪设想。迄今为止，研究者们已提出了如极限平衡分析法、块体理论法[1]、数值方法[2-3]等计算条件相对简单的传统方法，也提出了人工神经网络理论[4]、模糊数学理论[5]、灰色系统理论[6]、可靠性理论[7]等基于未确知性理论的计算方法。

本文采用FLAC3D数值模拟的方法对某水电站坝址区右岸古滑坡Ⅱ区堆积体进行稳定性分析， 由于FLAC3D采用的是显式差分法求解微分方程，不形成刚度矩阵，每一步计算仅需要很小的内存。在求解过程中通过叠加每一时步的小变形获得大变形求解，从而仅占用计算机很少的内存就可以模拟大变形，并且在模拟材料的塑性破坏和塑性流动方面具有优势。但是由于FLAC前处理功能较弱，建立复杂的三维模型非常困难，故本文选用ANSYS进行堆积体前期建模，并划分计算网格，再导入FLAC3D程序中计算，充分利用了两种软件的优点，优化分析过程。

2计算原理

21有限差分原理

有限差分法就是在利用数值计算方法求解偏微分方程时，用有限差分近似公式代替每一处导数，从而将求解偏微分方程的问题转化成求解代数方程的问题。

FLAC有限差分公式：由高斯散度定理有

∫snids=∫Afxi （1）

式中：∫s为封闭曲面上沿边界的积分； ni为曲面s的单位法向量；f为标量、向量或张量； xi为位置向量； ds为增量弧长；∫A为对表面积A积分。

定义f在面A上的梯度平均值为：

=1A∫AfxidA （2）

式中：<>表示求平均值。

将（2）代入（1）后得到

=1A∫sni fds （3）

对一个三角形子单元

=1A∑sniΔs （4）

式中：s为三角形某一条边的边长，等式右边的求和在三角形三条边上加总。取为相应边上的平均值[8]。

应力应变：用每一边速度矢量均值ui代替（4）式中的f，ui取各条边两端点的结点（即差分网格的角点）a和b的速度平均值，则

=12A∑N[（uai+ubi）niΔs]≈uixj （5）

同理可以求出的值。由几何方程就可以求出其应变率

eij=12uixj+ujxi （6）

由材料的本构方程和相应的边界条件，就可以求得应力增量。对各向同性的材料，有

σij=λδijθ+3μeij （7）

式中：λ，μ为拉梅常数；θ为体积应变，当时i=j，δij为1，否则，δij为零。

这样通过上述各式的迭代求解，就可以得出每一迭代时步对应的单元的应力应变值[9]。

22强度折减法

本文通过强度折减法来计算堆积体的稳定系数。强度折减法的原理是通过对岩土体的强度指标C和Φ值不断地折减，反复计算，直到塑性区贯通整个坡体，即边坡达到了临界的破坏状态，此时的折减系数就作为坡体的稳定系数Fs。

Cf=CFtrial （8）

Φf=tan-1（tanΦFtrial） （9）

式中：Cf为折减后的粘结力；Φf为折减后的摩擦角；Ftrial为折减系数。

3堆积体基本特征

古滑坡堆积体位于贵州某水电站选坝河段，分为四个区域，其中Ⅱ区岩土体结构最为复杂，控制底界以上物质主要以崩坡积块碎石土组成，故以Ⅱ区作为本文研究对象。堆积体所在河段河谷为不对称的“V”型河谷，总体为顺向坡，呈明显圈椅状地貌，顺河长600～700 m，横河宽约730 m，前缘临河，后缘高程为740～790 m。顺河方向上呈上升台阶状，上游侧地势低于下游侧40～70 m。总体上地形坡度约25°。后缘紧靠灰岩陡壁，区域内发育两条浅冲沟。覆盖层物质组成为崩坡积碎块石土为主。下伏基岩主要为J1夹层以下的二叠系下统栖霞组（P1q）中厚层夹薄层含炭质灰岩。

图1Ⅱ区地质分区平面示意图及13条剖面线位置

J1夹层为该区早期层状岩体失稳破坏的控制底界，夹层以上层状座滑体已经完全剥蚀搬运走，目前该区上部覆盖物质来源应为后缘陡壁后期崩塌堆积形成。平面示意图和地质横剖面示意图分别见图1和图2。

图2Ⅱ区典型地质横剖面示意图

4Flac3D数值模拟

41几何模型

模型边坡的主滑方向为S14°E，走向N76°E。模型以主滑方向的反方向为X轴正方向，竖直向上为Y轴正方向，走向NE方向为Z轴正方向。竖直方向坐标采用实际高程坐标，建立直角坐标系，如图3 所示。根据现场提供的地质资料，Ⅱ区堆积体模型选取了13个剖面，分别是Z=0，Z=57，Z=128，Z=215，Z=272，Z=325，Z=380，Z=438，Z=501，Z=561，Z=627，Z=697，Z=786（剖面线位置见图1），左边界略包括了Ⅰ区部分岩体。模型边坡底面高程190 m，坡顶高程850 m，沿走向方向长度为786 m，沿倾向方向长度为1400 m。模型在X轴上的范围是-200～1200 m，在Y轴上的范围是190～850 m，在Z轴上的范围是0～786 m。

图3Ⅱ区堆积体三维模型网格图

图4z=750 m剖面静水压力云图

模型分为5个部分，如图3，对应图上颜色，1，5部分为基岩，2为j1夹层，3为Ⅰ区座滑层状岩体，4为Ⅱ区表面堆积体，其中1，3，4，5部分选用solid45单元类型，夹层选用shell99单元类型。Ⅱ区堆积体模型共有15569个单元，4969个节点。计算模型边界约束形式为：侧边界只对水平方向进行约束，底边界在水平和竖直方向都进行约束，模型的上部边界取为自由面。FLAC程序中，节点速度是主要变量，所以选取模型的边界条件是通过约束模型边界的节点速度实现的，即模型底部边界的水平、竖直方向的速度约束和四周边界水平方向的速度约束。在程序中表达为边界上xvel，yvel设置0. 根据实际地勘报告，计算采用的各层物理力学参数如表1所示。

表1岩、土体物理力学指标

地层岩性密度ρ（g/cm3）天然饱和泊松比抗剪断（岩/岩）（°）c′/MPa饱和（°）饱和c′弹性模量/GPa

基岩26526802445110//200

座滑层状岩体230232022330401//05

堆积体19522002626572024232001

j1泥岩夹层195230025167003//035

地下水面的生成：由于在FLAC3D中直接生成符合勘察资料的地下水面比较困难，所以本文中建立模型时充分利用FLAC3D中的单元界面能自动依附于指定范围内模型表面生成的特性，生成水面。正常河水位为383 m，水库正常蓄水位为400 m。根据勘察报告，本文采用地下水位从坡脚处高程400 m开始，以5°向上倾斜，来近似模拟实际地下水位形态，如图4。沿着水位线，对三维模型采用切割命令，得到地下水面及其下部的几何实体，并剖分网格，导入FLAC3D，并以之为辅助单元，生成与地下水面空间形态一致的单元界面。采用FISH语言[10]遍历界面单元节点，生成水面，同时生成静水压力。

42计算结果分析

421主应力场规律分析

计算表明，坡体的最大主应力及最小主应力基本为层状分布，并在坡脚高程约400～460 m处出现应力集中，坡体大主应力σ1最大值为18139 MPa，小主应力值σ3最大值为10585 MPa。具体应力分布见图5及图6。

图5坡体大主应力σ1分布图

图6坡体小主应力σ3分布图

从剖面Z=650 m的应力分布图（如图7）可以看出，剖面附近的最大主应力（压应力），基本顺着坡面方向，并一直延伸到坡脚。而往边坡内部，最大主应力方向与水平轴的夹角逐步变大，直至铅直；由于岩层分界面的存在，使得其附近区域的最大主应力方向要比其他区域最大主应力方向变化大而且迅速，但并未影响主应力分布的总体走势。这些都表明边坡深部土体主要受铅垂方向的压应力作用，体现为受压屈服。

图7Z=650 m剖面主应力分布图

图8坡体塑性区分布图

422塑性区分布规律分析

从图8来看，坡体的屈服区域较集中，范围较小。坡体后缘（高程约780～850 m）分布有小块拉张塑性区，坡脚处（高程约380～450 m）分布零散的剪切塑性区。堆积体发生张拉剪切破坏的可能性很小，即使发生，也仅是局部区域，不会对堆积体整体稳定性造成重大影响。从剖面塑性区分布图（图9，图10）来看，仅仅在坡体的后缘和坡脚处出现了零星的塑性区，这表明堆积体处于正常的工作状态。

图9Z=210 m剖面塑性区分布图

图10Z=750 m剖面塑性区分布图

坡脚处的剪切塑性区位于表面堆积体内，厚度大约5～10 m。需要强调的是，计算结果显示的是以Mohr-Coulomb屈服准则为依据的塑性区分布情况。该屈服准则则认为材料进入屈服即破坏。

423剪应力（剪应变增量）规律分析

判断堆积体的（潜在）滑动面（带），可根据其剪应力（应变增量）来判断；剪应力较为集中或剪应变增量较大（绝对值）的部位，则为其（潜在）滑动面（带）。变形破坏也都沿此处发生；剪应力分布较为分散（均匀）或者剪应变增量较小或基本上没发生变化的部位，一般不会有潜在滑动面产生，因此，这些部位也不会发生较大的变形和破坏。应力集中带主要出现在3个范围内，坡体上下游侧坡脚处和j1夹层附近，它们是堆积体最有可能发生破坏的部位，见图11。J1夹层附近出现剪应变增量增高带，说明坡体有可能的破坏模式是沿着以j1夹层附近为滑面进行滑动，但由于计算得到的堆积体整体稳定系数Fs为110，所以发生整体性滑动的可能性不大；坡体上下游侧坡脚处的两个剪应变增量增高带，说明坡体可能沿土层内部局部滑弧滑动。Z=730剪应变增量剖面见图12。

图11整体剪应变增量展示图

图12Z=730 m剪应变增量剖面展示图

424位移场规律分析

坡体变形主要限于堆积体中、下部（600 m高程以下）区域。竖直方向位移总体表现为下沉，最大下沉区与水平方向最大位移区一致，局部沉降最大量值为51368 cm，其他下沉部位量值为0～1 cm（如图13）。

图13坡体竖直方向位移云图

出现局部最大沉降的区域有两处，第一处位于Ⅱ区下游侧，高程（Y坐标）范围大约位于420～425 m之间，走向方向（Z坐标）范围大约处于130～150 m之间，表面积小于100 m2；第二处位于Ⅱ区上游侧，高程（Y坐标）范围大约位于410～415 m之间，走向方向范围（Z坐标）范围大约处于760～770 m之间，表面积小于50 m2。水平方向局部最大位移量值为13.483 cm，出现区域与竖直方向局部最大沉降区一致，坡体其他部位的水平方向位移量值大都处于0～2 cm之间（如图14）。坡体最大变形区位置与塑性区位置相吻合。endprint

图14 坡体水平方向位移云图

5结语

1） 根据强度折减法，得到的堆积体整体稳定系数Fs为110，所以发生整体性滑动的可能性不大。文章对滑坡主应力，塑性区，剪应力和位移四个方面的计算结果进行了详细分析，并且四个方面的计算结果都相互吻合相互印证。

2） 可以直观地在模型上看出在蓄水条件下，上下游坡脚处的碎石土堆积层是Ⅱ区的主要不稳区域，具体在两个位置：①下游侧坡脚处，高程（Y坐标）范围大约位于420～425 m之间，走向方向（Z坐标）范围大约处于130～150 m之间，②上游侧坡脚处，高程（Y坐标）范围大约位于410～415 m之间，走向方向范围（Z坐标）范围大约处于760～770 m之间，对后期堆积体的加固治理方案提供了直观详细的位置。

3） 本文在模拟蓄水条件下的地下水位时利用了FLAC3D中的单元界面能自动依附于指定范围内模型表面生成的特性，通过建立辅助单元生成水面并同时生成静水压力，思路清晰明了，适用于比较复杂的坡体模型。

4） 利用FLAC3D软件对堆积体进行稳定性分析具有方便，经济，直观的优势，但是由于建立的模型不可能完全模拟出现实堆积体的所有特性，选取的计算参数也不可能完全与实际相符，所以需要结合其他技术手段综合评价，才能得出更为全面正确的结果。

参考文献：

[1]李素梅. 块体理论在岩石路堑边坡稳定分析中的应用[J]. 云南交通科技， 2003， 19（4）：14-16.

[2]张均峰， 丁桦.边坡稳定性分析的三维极限平衡法及应用[J]. 岩石力学与工程学报， 2005，24（3）：365-370.

[3]杨强， 朱玲， 薛利军. 基于三维多重网格法的极限平衡法在锦屏高边坡稳定性分析中的应用[J]. 岩石力学与工程学报， 2005， 24（增2）：5 313-5 318.

[4]陈昌彦， 王思敬，沈小克. 边坡岩体稳定性的人工神经网络预测模型[J]. 岩土工程学报， 2001， 23（2）：157-161.

[5]贾厚华， 贺怀建. 边坡稳定模糊随机可靠度分析[J]. 岩土力学， 2003， 24（4）：657-660.

[6]陈新民， 罗国煜. 基于经验的边坡稳定性灰色系统分析与评价[J]. 岩土工程学报， 1999， 21（5）：638-641.

[7]祝玉学. 边坡可靠性分析[M]. 北京：冶金工业出版社， 1993.

[8]黄润秋，许强. 显式拉格朗日差分分析在岩石边坡工程中的应用[J]. 岩石力学与工程学报， 1995， 14（4）：346-354.

[9]邹栋， 郑宏. 快速拉格朗日法及其在边坡稳定性分析中的应用[J]. 矿业研究与开发， 2005， 25， （5）：80-83.

[10]陈育民， 徐鼎平. FLAC/FLAC3D基础与工程实例[M]. 北京： 中国水利水电出版社， 2008.

（责任编辑：）endprint