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

科学技术创新2020.35
利用Flac3D 进行库区消落带边坡稳定性研究
李望成冉益铭
(重庆科技学院建筑工程学院，
重庆401331)三峡库区建成以来，每年冬天蓄水发电期和夏季防洪期的
水位分别为175米和145米，30米的水位差形成了消落带，
消落带边坡每年都会经历一次干湿循环，
长期发展下去，会对消落带边坡稳定性造成一定的影响，
消落带边坡岩质多为砂岩，砂岩在经过一定次数的干湿循环后，
对其基本物理力学性质造成一定的影响，抗劈裂强度、内摩擦角和粘聚力等影响边坡稳定性的参数会随着干湿循环次数的增加而降低[1]，从而影响到边坡
安全系数，造成边坡失稳的情况发生。

长江三峡大坝建成后，
水库蓄水及库岸的二次改造会对长江沿岸岸坡造成深刻的影响，
包括崩塌和坍塌、滑动和滑移、
蠕移、沉降或沉陷、碎石流与泥石流及洪冲垮塌等[2];随着库区水位升降变化，
滑坡体安全系数随着岩体干湿循环的次数呈指数降低，最终造成边坡失稳[3]；前4
次干湿循环作用将造成抗剪强度参数劣化幅度占总劣化度的
75%左右，在4次之后劣化幅度将趋于平缓；[4]
；随着干湿循环次
数增加，坡体下滑力与抗滑力的平衡被打破而造成坡体失稳，而后又恢复平衡状态，坡体间的平衡状态会随着干湿作用的次数增加而更加容易被打破，造成边坡稳定性安全系数降低[5]。

1Flac3D 的原理与步骤
Flac3D 是美国ITASCA 公司开发的仿真模拟软件，利用源
于流体力学的拉格朗日差分法，
研究在规定时间内的位移与应力的变化[6]，Flac3D 把模型划分为由多个六面体网格组成，网格数量根据需求来设置，更多的网格会造成更久的计算时间，
但可以更精确的得到计算结果；Flac3D 动力学分析基于显式有限差分方案，使用从周围区域的实际密度（而不是用于静态解决方案的虚拟质量）导出集总网格点来求解完整的运动方程。

FLAC3D 采用完全非线性分析方法，加入了“滞后阻尼”的因子
对动荷载产生的非线性材料累积变形具有很好的还原能力，
动荷载分为Fish 函数和Table 数表两种类型，前者是在面施加荷载，后者是在点上施加荷载。

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

2数值模拟分析
假设库区边坡水位匀速升降，岩体均质且符合弹性本构关
系，基本物理参数在每一次水位升降后出现劣化，
不考虑边坡坡面水位升降造成的冲刷与侵蚀以及风化后岩体的损失。

本文
建立的模型为六面块体和楔形体拼接而成，把模型划分为六面
体的网格单元，本模型一共有81000单元，整个模型采用摩尔库
伦弹塑性模型本构关系进行边坡干湿循环位移量的模拟，模型
的三个垂直方面底面及顶面固定，
坡面施加一个循环的梯度荷载。

计算过程中，地面和顶面的横向动荷载随着水位的升降变化，水位会在一年的时间内从年初的最高水位下降到夏季的最低水位在上升到年末的最高水位，每经历一次一年内的水位升
降，边坡岩体的基本物理参数就会发生一次劣化，
从而影响到后续边坡位移的变化；
3计算参数
摩尔库伦弹塑性模型材料属性参数包括体积模量、
粘聚力、剪胀角、内摩擦角、泊松比、剪切模量、
抗拉强度、弹性模量等，本次模拟仅选用其中一部分材料参数进行计算，
具体参数如表1：表1计算参数
4计算过程及结果分析
本文主要从以下步骤进行库区边坡模型计算：
步骤一：建立边坡模型，按尺寸划分网格，
建立弹性本构关系，设置边界条件和初始重力场并进行计算。

图1初始重力场Z 向位移云图
由图1：边坡受到的重力遵循从下往上越来越大，
在边坡下半部分对底面和除了坡面以外的另外三面施加约束后，Z 方向的位移在地面为零，而越靠近顶面位移量的越大，最大位移量在
摘要：运用Flac3D 对库区消落带边坡进行稳定性分析，根据消落带边坡水位升降变化情况，结合室内试验得到不同干湿循
环次数后岩石的物理力学参数，模拟消落带边坡随着水位升降情况下的位移状态。

结果表明：
通过对库区消落带边坡进行有限元差分数值模拟分析可知，边坡的横向变形由上至下逐渐增加，横向的最大位移量在边坡底面为4.2x10-2，最小位移量在边坡顶部为
0.5x10-2；纵向的变形从左至右逐渐减小，边坡坡面中部存在一个向上的变形为5.6x10-3，而在边坡的内部具有最大的向下变形为4.7x10-3。

关键词：边坡；循环；稳定性中图分类号:U416.1+6文献标识码:A 文章编号:2096-4390(2020)35-0122-02密度/(3kg m ) 2500 粘聚力/(kPa ) 19 内摩擦角/(°) 30.53 体积模量/(MPa ) 100 弹性模量/(MPa ) 17.5 剪切模量/(MPa ) 77 抗拉强度/(MPa )
0.5
122--
2020.35科学技术创新
顶面原理坡面一侧；应力方面却是从上到下逐渐增加。

步骤二：以最高水位点作为计算的起始点，当边坡水位上涨到最高点
时，水作用于边坡坡面一个递增的梯度力，
因此在边坡坡面施加一个指向X 轴正方向的梯度力。

图2起始梯度力Z 向位移云图
由图2可知：梯度力对边坡沉降有一定影响，
边坡在梯度力和重力的影响下产生的沉降相比于在重力的作用下的沉降量更少；而在应力方面，梯度力对边坡坡面产生了负方向的应力，且在边坡底角的应力更加集中了。

步骤三：作用在边坡坡面的梯度力会随着水位升降循环增加减小，且在水位升降过程中，边坡也进行了一次干湿循环，而干湿循环会对岩质边坡产生劣化作用，降低岩体的弹性模量和泊松比等基本物理力学参数，
所以在每一次水位升降之后对基本力学参数进行调整，
尽量还原现场实际情况。

如图3为循环梯度荷载作用下边坡在X 的应力云图，由于梯度荷载循环作用，在X 方向上，边坡顶部的应力
向X 的正方向，底部应力向X 的负方向，呈现出自上而下逐渐
减小的应力。

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

图4X 向位移云图
图6为边坡坡面上（0，0，0）、（7.5，0，15）、（15，0，30）三个点在X 和Z 方向上的应变情况，从图9可以看出（15，0，30）点Z 方向的位移与位移云图结果相同，从三点X 方向变形曲线可以看出在计算到10000多步之后边坡就已经趋近平衡了；而三个
点的X 方向变形也可以从图中看出从上到下是逐渐增加的，
最大的变形接近，顶点和中间点的变形分别达到了、；顶点位置的Z 方向的变形由于没有Z 方向的荷载变形为0。

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

参考文献
[1]傅晏,刘新荣,张永兴,胡元鑫,谢应坤.水岩相互作用对砂岩单轴强度的影响研究[J].水文地质工程地质,2009,36(06):54-58.[2]Zeng Sheng,Li Zhencun,Chen Hanjie,Zhao Jian.Strength
attenuation law of Red Sandstone under dry-wet cycle and its engineering application J.Proceedings of the Changsha University of Science and Technology,2011,804:18-23.
[3]邓华锋,李建林,朱敏,王孔伟,王乐华,邓成进.饱水-风干循环作用下砂岩强度劣化规律试验研究[J].岩土力学,2012,33(11):3306-3312.
[4]Deng Huafeng,Li Jianlin,Zhu Min,Wang Kongwei,Wang Lehua,Deng Chengjin.Experimental study on strength degradation of SANDSTONE UNDER WATER -AIR -DRYING CYCLES J.Geomechanics,2012,3311:3306-3312.
[5]张振华,王野.水库运行期岸坡消落带红砂岩抗剪与抗压强度劣化机制[J].岩土工程学报,2019,41(07):1217-1226.
[6]Zhang Zhenhua,Wang Ye.Mechanism of Shear and compressive strength deterioration of red sandstone in the
fluctuating zone of bank slope during reservoir operation J.
Geotechnical engineering,
2019,4107:1217-1226.
图5Z 向位移云图
24.210 20.510 2
2.710
123--。