利用Flac3D 进行库区消落带边坡稳定性研究

李望成 冉益铭

(重庆科技学院建筑工程学院，重庆401331)

三峡库区建成以来，每年冬天蓄水发电期和夏季防洪期的水位分别为175 米和145 米，30 米的水位差形成了消落带，消落带边坡每年都会经历一次干湿循环，长期发展下去，会对消落带边坡稳定性造成一定的影响，消落带边坡岩质多为砂岩，砂岩在经过一定次数的干湿循环后，对其基本物理力学性质造成一定的影响，抗劈裂强度、内摩擦角和粘聚力等影响边坡稳定性的参数会随着干湿循环次数的增加而降低[1]，从而影响到边坡安全系数，造成边坡失稳的情况发生。长江三峡大坝建成后，水库蓄水及库岸的二次改造会对长江沿岸岸坡造成深刻的影响，包括崩塌和坍塌、滑动和滑移、蠕移、沉降或沉陷、碎石流与泥石流及洪冲垮塌等[2];随着库区水位升降变化，滑坡体安全系数随着岩体干湿循环的次数呈指数降低，最终造成边坡失稳[3]；前4次干湿循环作用将造成抗剪强度参数劣化幅度占总劣化度的75%左右，在4 次之后劣化幅度将趋于平缓；[4]；随着干湿循环次数增加，坡体下滑力与抗滑力的平衡被打破而造成坡体失稳，而后又恢复平衡状态，坡体间的平衡状态会随着干湿作用的次数增加而更加容易被打破，造成边坡稳定性安全系数降低[5]。

1 Flac3D 的原理与步骤

Flac3D 是美国ITASCA 公司开发的仿真模拟软件，利用源于流体力学的拉格朗日差分法，研究在规定时间内的位移与应力的变化[6]，Flac3D 把模型划分为由多个六面体网格组成，网格数量根据需求来设置，更多的网格会造成更久的计算时间，但可以更精确的得到计算结果；Flac3D 动力学分析基于显式有限差分方案，使用从周围区域的实际密度（而不是用于静态解决方案的虚拟质量） 导出集总网格点来求解完整的运动方程。FLAC3D 采用完全非线性分析方法，加入了“滞后阻尼”的因子对动荷载产生的非线性材料累积变形具有很好的还原能力，动荷载分为Fish 函数和Table 数表两种类型，前者是在面施加荷载，后者是在点上施加荷载。

Flac3D 进行实际工程建模分析步骤如下：（1）根据工程中的地形地貌，进行一定程度的简化和省略，建立出工程中关键地形地貌，根据计算精度不同和需要分组与划分网格；（2）定义模型材料塑性并根据实际地形定义模型的边界条件，施加初始重力，完成初始应力计算；（3）根据实际情况施加动荷载，在循环的过程中改变岩体基本物理参数；（4）设置边坡稳定性安全系数，当破边稳定性安全系数超过设定值时停止循环。

2 数值模拟分析

假设库区边坡水位匀速升降，岩体均质且符合弹性本构关系，基本物理参数在每一次水位升降后出现劣化，不考虑边坡坡面水位升降造成的冲刷与侵蚀以及风化后岩体的损失。本文建立的模型为六面块体和楔形体拼接而成，把模型划分为六面体的网格单元，本模型一共有81000 单元，整个模型采用摩尔库伦弹塑性模型本构关系进行边坡干湿循环位移量的模拟，模型的三个垂直方面底面及顶面固定，坡面施加一个循环的梯度荷载。计算过程中，地面和顶面的横向动荷载随着水位的升降变化，水位会在一年的时间内从年初的最高水位下降到夏季的最低水位在上升到年末的最高水位，每经历一次一年内的水位升降，边坡岩体的基本物理参数就会发生一次劣化，从而影响到后续边坡位移的变化；

3 计算参数

摩尔库伦弹塑性模型材料属性参数包括体积模量、粘聚力、剪胀角、内摩擦角、泊松比、剪切模量、抗拉强度、弹性模量等，本次模拟仅选用其中一部分材料参数进行计算，具体参数如表1：

表1 计算参数

4 计算过程及结果分析

本文主要从以下步骤进行库区边坡模型计算：步骤一：建立边坡模型，按尺寸划分网格，建立弹性本构关系，设置边界条件和初始重力场并进行计算。

图1 初始重力场Z 向位移云图

由图1：边坡受到的重力遵循从下往上越来越大，在边坡下半部分对底面和除了坡面以外的另外三面施加约束后，Z 方向的位移在地面为零，而越靠近顶面位移量的越大，最大位移量在顶面原理坡面一侧；应力方面却是从上到下逐渐增加。步骤二：以最高水位点作为计算的起始点，当边坡水位上涨到最高点时，水作用于边坡坡面一个递增的梯度力，因此在边坡坡面施加一个指向X 轴正方向的梯度力。

图2 起始梯度力Z 向位移云图

由图2 可知：梯度力对边坡沉降有一定影响，边坡在梯度力和重力的影响下产生的沉降相比于在重力的作用下的沉降量更少；而在应力方面，梯度力对边坡坡面产生了负方向的应力，且在边坡底角的应力更加集中了。步骤三：作用在边坡坡面的梯度力会随着水位升降循环增加减小，且在水位升降过程中，边坡也进行了一次干湿循环，而干湿循环会对岩质边坡产生劣化作用，降低岩体的弹性模量和泊松比等基本物理力学参数，所以在每一次水位升降之后对基本力学参数进行调整，尽量还原现场实际情况。如图3 为循环梯度荷载作用下边坡在X 的应力云图，由于梯度荷载循环作用，在X 方向上，边坡顶部的应力向X 的正方向，底部应力向X 的负方向，呈现出自上而下逐渐减小的应力。

图3 Z 向应力云图

图4、图5 为X、Z 两个方向的位移云图，X 方向上的位移整体呈现边坡坡底位移量最大而坡底位移最小，位移由下至上逐渐变大的趋势，在Z 方向的位移为左右分化的趋势，边坡坡面中部位置出现了最大的Z 正方向位移，坡面顶部和底部都仅有少量的位移，而在边坡内部中部位置具有最大的Z 负方向的位移，也就是在内部中部具有最大的沉降量。

图4 X 向位移云图

图5 Z 向位移云图

图6 不同点的应变趋势

5 结论

本文通过FLAC3D6.0 分析了库区边坡消落带水位周期性涨落情形下边坡整体变形，从而得到了以下的结论：通过对库区消落带边坡进行有限元差分数值模拟分析可知，边坡的横向变形由上至下逐渐增加，横向的最大位移量在边坡底面为4.2x10-2，最小位移量在边坡顶部为0.5x10-2；纵向的变形从左至右逐渐减小，边坡坡面中部存在一个向上的变形为5.6x10-3，而在边坡的内部具有最大的向下变形为4.7x10-3。