河谷边坡形成过程中变形破坏FLAC3D数值模拟.doc

河谷边坡形成过程中变形破坏FLAC3D数值模拟

摘要：本文结合XX滑坡所处地质环境，利用FLAC3D数值分析软件，模拟软硬互层斜坡在河谷边坡形成过程中的应力变化及变形破坏特征。分析得出，初始应力场中，σv垂直应力均为压应力，近水平平行分层分布，由坡顶向深部基本呈线性递增，增加梯度大致与岩体平均重度等同；最大主应力方向近垂直方向，随深度线性递增。随河谷下切，河谷边坡岩体在卸荷作用下引起内部应力重分布，局部应力集中，最大主应力迹线发生偏转，越接近临空面越与之平行，最小主应力则与之正交；坡体由于开挖卸荷，坡表回弹变形明显，根据监测点位移值，坡顶到坡脚变形呈由小变大的趋势。
关键词：FLAC3D；河谷边坡；软硬互层；变形破坏；数值模拟

0 引言
天然斜坡或人工边坡形成过程中，岩体内部原有的应力状态将随着过程的进行而发生变化，引起应力的重分布和应力集中等效应。斜坡岩体为适应这种新的应力状态，将发生不同形式和不同规模的变形与破坏，斜坡变形破坏过程和它所造成的不良地质环境可对人类工程活动带来十分严重的危害，如滑坡、崩塌等的发生[1]。河谷边坡形成过程中，由于河谷下切引起的卸荷回弹，将引起边坡坡体的应力分异及变形破坏，目前河谷下切过程的应力场研究一般采用以实测地应力资料为基础的定性-半定量分析和数值分析两种方法。由于可以处理复杂的边界条件及材料的非均质等问题, 特别是近年来计算机硬件技术的快速发展以及一系列界面友好、计算能力强、操作简单的分析软件, 如3D-σ、FLAC3D等的相继推出, 数值计算方法已经成为河谷( 岩土体) 应力场分析的最主要途径。数值方法在研究特定岩体, 特别是具有复杂边界的计算对象时, 具有明显的优势, 从采用的控制方程、计算过程及结果形式来看, 数值法也是定量的, 但是由于岩体结构的高度复杂性及勘探精度的限制, 参数只能近似取值, 因此从本从质上看, 数值法仍然是定性的[2]。
本文对于XX滑坡所处地质环境的应力场研究以及河谷边坡形成过程中的变形破坏特征的研究，诣在为进一步分析XX滑坡的地应力环境及形成机制以及判断滑坡发展阶段进行基础性研究，同时，河谷边坡最大主应力大小及方向的探明与成坡过程中卸荷回弹水平的初步分析，对于弄清XX滑坡的形成过程紧密相关。
1 工程概况
XX滑坡发生于约30年前，该滑坡位于XX省XX市XX县XX乡，地处XX省东北部，地理坐标为北纬28°19′06.41″，东经108°02′09.20″。据作者2012年4月、11月，2013年4月三次现场地质调查发现，该滑坡目前变形特征明显，主要表现为多处出现拉陷槽，并且在分布上具有规律性，走向均大致为SW-NE走向。该滑坡的变形乃至再次失稳将影响到该地区24户人家，153人的生命财产安全，故对其进行研究具有实际意义。滑坡后缘高程约为962~980m左右，前缘高程870m，垂直高差约为100m，为缓倾外顺层软硬互层岩质边坡滑坡，基岩产状约为300°~310°∠7°~10°。滑坡范围内涉及的岩性，上部为下古生界志留系中下统韩家店群马脚冲组、溶溪组中的黄绿、灰绿色页岩，化石主要为腕足类，数量少，属种单调，下部为下古生界志留系下统石牛栏组，以灰、深灰色中厚层至块状生物骨屑及介屑灰岩为主，时具角砾状及瘤状构造，偶夹泥质灰岩及钙质泥岩，富含珊瑚、腕足类及三叶虫等化石[3]。在滑坡前缘以下，为石牛栏组与松坎组整合接触，呈现厚层灰岩与页岩软硬互层分布，最底部为深切河谷，滑坡后缘与河谷底高差约为280~300m，图1为XX滑坡地质环境示意。

图1XX滑坡地质环境
Figure1 Geological environment of xiangshuping landsilde
2 FLAC3D简介[4~6]
FLAC3D（3-D Fast Lagrangian Analysis Code)是由美国ITASCA公司开发的。目前FLAC有二维和三维计算程序两个版本，FLAC3D是一个三维有限差分程序，是二维的有限差分程序FLAC的扩展，能够进行土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析。调整三维网格中的多面体单元来拟合实际的结构，单元材料可采用线性或非线性本构模型，在外力作用下，当材料发生屈服流动后，网格能够相应发生变形和移动（大变形模式）。FLAC3D采用的显式拉格朗日算法和混合-离散分区技术能够非常准确地模拟材料的塑形破坏和流动。由于无须形成刚度矩阵，因此，基于较小内存空间就能够求解大范围的三维问题。FLAC3D有以下几个优点：（1）对模拟塑性破坏和塑性流动采用的是“混合离散法”。这种方法比有限元中通常采用的“离散集成法”更为准确合理。（2）即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍。（3）采用了一个“显示解”方案。模拟大变形问题几乎并不比小变形问题多消耗更多的计算时间，因为没有任何刚度矩阵要被修改。FLAC3D做计算分析的一般步骤与大多数程序采用数据输入方式不同，FLAC3D采用的是命令驱动方式。命令控制着程序的运行，建立FLAC3D模型必须进行三方面的工作：（1）有限差分网格；（2）本构特性与材料性质；（3）边界条件与初始条件。完成上述工作后，可以获得模型的初始平衡状态，然后进行工程开挖或改变边界条件来进行工程的响应分析。
3 FLAC3D数值计算[7~13]
3.1模型建立
（1） 模型范围：根据现场实际调查和已做出的剖面图，选取穿越滑坡中部及近垂直河谷方向的典型剖面，建立地质模型见图2，在建立模型时，为了尽可能减小底部边界约束对计算结果的影响，底部边界离河谷谷底取约0.7倍的边坡高度200m，在确立横向边界时，为了将整个边坡包括在内并尽量减小边界影响，左侧取河谷谷底处作为边界，右侧取至边坡顶以外约40m处作为边界，取长度方向左端点为坐标原点，水平向右为x轴正方向，竖直向上为y轴正方向。模型在x方向长度为800m，y方向高度为480m，y方向高度由坐标原点向上分别为下部谷底以下深度200m范围，中部河谷边坡高差160m，上部XX滑坡范围高差120m。在CAD中绘制平面模型，导入ANSYS中拉伸Z方向50m厚建成三维模型，并网格化，利用ANSYS-FLAC3D接口程序，按照平面应变建立计算模型，见图3，其中模型长x宽x高=800mx50mx480m（计算范围内相应的坐标为：x方向为0~800m;y方向为0~480m;z方向为0~50m）。
（2）岩性及本构模型：模型最底层为厚层灰岩，向上依次分层为页岩、灰岩、页岩、灰岩，以上全部为页岩（XX滑坡以内），XX滑坡前缘以下为软硬互层岩体结构，见图2；计算模型采用岩土工程中应用最为广泛的Mohr-Coulomb弹塑性模型,该模型包含剪切和拉伸两个准则。
（3）边界条件及初始条件:下部固定约束，左右两侧法向约束,厚度方向前后侧面法向约束，上坡面为自由边界，全场自重应力通过重力加速度g=9.81m/s2施加。
（4）计算模型如图3所示,经过网格剖分,模型共有34112个单元,21663个节点。
（5）如图2所示，首先对模型进行初始应力场计算，计算完成，将位移恢复为零，然后分别按图中Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ的步骤进行开挖模拟计算，每开挖一步计算一步，下一步计算开始前先将位移恢复为零，分析应力场分布及位移变化。

图2 地质模型
Figure2 Geological model
图3 计算模型
Figure3 Calculation model
3.2 岩体物理力学参数
根据现场地质调查，XX滑坡主要由页岩组成，前缘与石牛栏组整合接触带也有少量灰岩，前缘至河谷底为灰岩、页岩软硬互层结构。岩体物理力学参数是综合XX相近滑坡的钻孔、探槽所取得试样的物理力学试验资料以及查阅相关规范，结合工程经验取得，见表1。
表1 岩体物理力学参数
Table1 Physical and mechanical parameters of rock mass
岩体
密度
kg/m3
弹性模量
MPa
泊松比
黏聚力
MPa
内摩擦角
°
抗拉强度
MPa
页岩
2420
1800
0.28
1.67
30
2.4
灰岩
2670
2800
0.24
2.5
40
2.5
3.3 计算分析
（1）初始场分析
边坡工程数值计算中，初始地应力场的模拟是必须首先关注的问题。数值模拟的初始地应力场是否与实际地应力场吻合较好，是决定边坡工程数值模拟能否成功的基本条件通常认为，初始应力场主要由岩体自重和地质构造应力产生，其他影响因素如温度产生的应力在边坡工程中往往忽略不计，而在该边坡的地表浅覆层环境下，地质构造应力很低，仅考虑岩体自重产生的应力场，往往可以满足计算的精度要求[10]。据黄润秋等[9],垂直应力σv随深度线性增加，增加梯度大致与岩体重度等同：σv=γh=ρ平gh≈2500x10x480=12 (Mpa），其中，ρ平为所取页岩与灰岩密度的大致平均值，理论计算值与初始场数值计算结果基本吻合，见图4。由于以自重应力为主，没有考虑水平构造应力，最大主应力方向主要受自重应力影响，为垂直方向，最小主应力与之正交，由于边界效应，在坡顶及两侧主应力迹线略有异常现象，总体上，从边坡表部到坡体内部应力值呈现逐渐增大的趋势，这是因为越往坡体深部，对应的上覆岩体的厚度越大，最大主应力量值最大约为12Mpa，最小主应力值恰好相反，见图5~7。

图4 初始场垂直应力等值云图
Figure4 Vertical stress equivalent cloud of initial field
图5 最大主应力等值云图
Figure5 Maximum principal stress equivalent cloud

图6 最小主应力等值云图
Figure6 Minimum principal stress equivalent cloud
图7 主应力迹线
Figure7 Principal stress trace
（2）河谷边坡形成过程变形破坏特征分析
初始计算完成，现对河谷边坡形成的下切过程进行模拟分析，下切分为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ五级开挖的顺序进行，每开挖一步，分析河谷边坡的应力分布变化与变形破坏情况，分析得出：随着下切，在河谷边坡形成过程中，边坡岩体内的应力状态由于边坡临空面附近的应力重分布及应力集中作用，变得复杂。由于河谷边坡临空面的形成，使得附近原来处于应力平衡状态的边坡体遭受临空面附近局部卸荷的影响，发生边坡体向临空方向的膨胀回弹，主应力迹线也发生明显的变化，具体表现为：
1）主应力迹线在河谷边坡临空面附近发生明显的偏转，越靠近临空面，最大主应力越与之平行，最小主应力则与之正交，见图8；
2）最大主应力表现为由临空面向坡体内逐渐增大，在临空面处达最大值约11.2Mpa，最小主应力则相反，由坡体内向临空面逐渐减小，在临空面处接近为零，见图9；
3）在下切初期（第一步开挖）临空面坡脚出现剪应力集中，随继续下切，剪应力集中区由坡脚呈圆弧形向坡体左下侧转移，最后在边坡坡脚下约40~50m深度处形成一大范围剪应力集中区域，剪应力值最大约6.8Mpa，见图10；
4）临空面附近的卸荷回弹表现为位移值的变化，每一步开挖，坡体向临空面方向均有部分回弹，监测点1到监测点4为边坡坡顶到坡脚依次向下取的点，见图3。由图12~15可以分析得出，在开挖过程中，由坡顶向坡脚，回弹量呈由小变大的规律，开挖到坡脚处时，变形是最大的。同时可以看出，监测点4在每一步开挖中的位移值也几乎是四个监测点中最大的，最后一步开挖完成时位移约为45mm，见图11，说明在边坡卸荷影响范围内，深处的卸荷回弹更明显。

图8 主应力迹线偏转
Figure8 Principal stress trace deflection
图9 开挖完成最大主应力等值云图
Figure9 Maximum principal stress equivalent cloud after excavation

图10开挖完成剪应力等值云图
Figure10 Shear stress equivalent cloud after excavation
图11 开挖完成位移云图
Figure11 Displayment cloud after excavation

图12 监测点1位移图
Figure12 Displayment of monitoring point1
图13 监测点2位移图
Figure13 Displayment of monitoring point2

图14 监测点3位移图
Figure14 Displayment of monitoring point3
图15 监测点4位移图
Figure15 Displayment of monitoring point4
4 结论
本文是为了进行XX滑坡变形破坏过程及机制研究进行的基础性研究，弄清楚滑坡所处软硬互层结构下，在河谷边坡形成历史过程中的应力状态变化具有基础性意义。作者基于FLAC3D数值模拟软件分析得出:
(1)初始场中，垂直应力σv随深度线性增加，增加梯度大致与岩体重度等同；最大主应力方向主要受自重应力影响，为垂直方向，最小主应力与之正交；由于边界效应，在坡顶及两侧主应力迹线略有异常现象；总体上，从边坡表部到坡体内部最大主应力呈现逐渐增大的趋势，最大约为12Mpa，最小主应力值恰好相反。

(2)在下切过程中河谷边坡临空面附近发生了应力分异，主要表现为应力重分布与局部的应力集中，主应力迹线在河谷边坡临空面附近发生明显的偏转，越靠近临空面，最大主应力越与之平行，最小主应力则与之正交；同时，最大主应力由临空面向坡体内逐渐增大，在临空面处达最大值约11.2Mpa，最小主应力则相反，由坡体内向临空面逐渐减小，在临空面处接近为零。
（3）剪应力集中带在河谷形成过程由坡顶向坡脚呈圆弧形转移，最后在坡脚以下深度约40~50m处形成一最高剪应力集中区域，剪应力值最大约6.8Mpa。
（4）临空面附近的卸荷回弹表现为位移值的变化，每一步开挖，坡体向临空面方向均有部分回弹，在开挖过程中，由坡顶向坡脚，回弹量呈现由小变大的规律，坡脚处的变形最大，同时，在边坡卸荷影响范围内，深处的卸荷回弹更明显。
FLAC3D NUMERICAL SIMULATION ON DEFORMATION AND DESTRUCTION
IN THE PROCESS OF INCISED VALLEY
Abstract: This paper combines the geological environment of Xiangshuping landslide in guizhou, using FLAC3D numerical analysis software, simulate the stress changes and deformation failure features of hard and soft interbed rock during the process of incising valley. The analyse shows, in the initial state, sigma v vertical stress is compressive stress, nearly horizontal parallel hierarchical distribution, and almost linear increases from the top to the deep, the increasing gradient is roughly the same with the average unit weight of rock mass; The direction of maximum principal stress is nearly vertical, a linear increasing with depth. Along with the incised valley, the internal stress of the incised valley redistribute and show some local stress concentration under the unloading effect, the trace of maximum principal stress deflect, the closer they get to the surface,the more paralleled to the surface,while the the minimum principal stress vertical to it. Due to excavation,the rebound deformation is obvious,and as the displacement of the monitoring points show,the value increase from the slope top to the toe.
Keywords: FLAC3D;Incised valley;Hard and soft interbed rock;Deformation and destruction;Numerical simulation