ANSYS有限元数值模拟在某开挖高边坡稳定性分析中的应用

许容 王辉

摘要：ANSYS有限元数值可以很好地模拟岩土体的力学性能，甚至通过选取适当的输入参数和计算模型，也能模拟出节理裂隙、软弱夹层、活动性断裂等地质情形。文章结合某开挖高边坡工程实例，利用ANSYS有限元软件，通过应力、变形以及能量计算分析边坡的稳定性。

关键词：ANSYS有限元;数值模拟;稳定性

0 引言

边坡稳定分析方法很多，常用的方法主要是通过应力、变形以及能量计算分析边坡稳定性。其采用矩阵分析模式，更加方便编程分析。应当注意的是岩土工程在做有限元分析时，不应当过度注重其具体量值的大小，应更加注重应力应变的分布情况与相对的变化情况。

1 边坡工程地质概况

该边坡位于云南省境内，构造主要受北西向和近南北向构造控制，区域地质构造背景十分复杂。该边坡区2#倾倒蠕变岩体在边坡冲沟下游侧1 570～1 840 m高程处发育。变形体表层覆盖有块石、碎石混粉土的坡积层;倾倒岩体主要为板岩（J2h）;下伏基岩为变质火山角砾岩夹片理化变质凝灰岩（T3xd）。

2 模型建立及数值模拟

2.1 有限元模型建立

选取开挖边坡横剖面，该剖面位于正面开挖边坡中段，剖面线总长535 m，前缘高程1 540 m，后缘高程1 960 m。右岸边坡开挖至缆机平台，高程为1 678 m（见下页图1）。

本文通过将剖面的二维极限平衡的渗流场通过插值的方式赋予网格中的各个节点，具体步骤如下：

首先将对各个工况下的渗流场进行模拟，得出某个时刻的孔隙水应力图（见下页图2）。提取图中各个坐标和孔隙水压力值，通过插值的方法将坐标输入ANSYS FLUENT模块中，使用TB、PM命令制定的渗透系数和BIOT系数定义孔隙介质属性。然后，通过fluent输出.cdb格式文件，最后再导入ANSYS中。

ANSYS有由頂向下和由底向上两种建模方式。本论文主要使用由底向上方式建立边坡的主体模型：首先提取坝址区CAD地质图中的等高线数据，存为XYZ坐标文件，然后导入到Surfer，选取建模区域内的坐标，通过克里格插值法将这些不规则高程点离散为均匀的网格点，导出这些网格数据为文本，最后在Excel中处理为ANSYS可以识别的APDL格式文件。再根据此初步模型，将上述APDL文本输入ANSYS，生成地表关键点，再通过“由底向上”的方式建立包括下切地层的完整边坡模型，然后依次分割出各个下切地层、分级开挖体和不同卸荷层等模型分区，最后对几何模型划分网格。采用智能剖分（SmartSize）和限定单元尺寸的方法优化网格，多次尝试后划分出合适的网格模型。

一般情况下，ANSYS模型越大、越精细，其求解精度就会越高，但是求解的时间也越多。根据本论文的需要和求解时间的限制，本论文建立模型的大小为沿河流方向（X轴）为364 m，垂直于河流方向（Y轴）为363 m，由于边坡在1 678 m高程，山顶高程为1 906 m（见图3），考虑到计算的合理性，建模的最低高程应在边坡高程以下，以防止在模型底部产生的应力集中影响计算结果。模型的顶部应至坡体顶部，能更合理地计算岩土体自重产生的应力应变。综上所述，建模高程（Z轴）取1 500～1 886 m中间的386 m高程。考虑模型开挖的过程，在模型开挖面以及开挖高程的同时建模。

模型使用SOLID95单元四面体剖分岩土体部分，SOLID95单元具有塑性、蠕变、单元生死、初始应力输入等适合岩土计算的特性。设置法向接触刚度为1.0，防止过大的接触面嵌入并且保证数值计算的稳定性;使用拉格朗日乘子法计算，设置侵入容差为0.05，以定义允许的最大侵入值;设置力学表面作用模式为法向单相接触，使得模型在计算时，如果接触处于分离状态，法向压力为零;单元网格划分时，依然综合考虑了计算精度与计算效率的问题，在坡体表面、开挖区域以及断层区域使用较小网格剖分，在模型的底部，由于形状规则，使用较大的网格剖分，模型单元数为10 103，节点数为37 460（见图4）。

模型分析对于力的容差为0.005，对于位移的容差为0.05，设置使用L2范数检查SRSS值，位移允许的最小值为1.0。

2.2 数值计算参数

综合前期勘察报告和作者现场跟踪调查成果，确定了不同岩层的数值计算参数取值，具体见表1。

3 有限元模拟计算分析

3.1 天然边坡应力应变分析

边坡的应力重分布即二次应力是天然状态边坡应力场特征的研究关键，对于高地应力下的变形破坏具有很高的研究意义。原因在于天然应力场对边坡开挖后应力重分布的趋势和结果起决定作用;同时在天然应力场的影响下形成的边坡内部应力增高带和应力集中也是影响边坡稳定性的因素。

3.1.1 应力场分析

天然状态下边坡大主应力场主要来源于岩体自重力，随着岩体深度的变化，分布较均匀，主要为压应力，局部位置为拉应力。边坡的坡面为临空面，初始的最大主应力为沿坡面切线方向受重力影响指向坡面以下，表现为拉应力集中在边坡顶部附近下游部位，此处为J2h千枚状泥质板岩，且后缘部位为2#倾倒变形岩体，最大主应力与最小主应力都表现为较小应力，与现实相符合。2#倾倒变形岩体附近最大主应力大小为0.82～1.33 MPa，坡底部表现出由重力引起的较大的压应力，大小为1.25～3.23 MPa;坡面最小主应力为拉应力，大小为0～0.27 MPa，在边坡底部有较大的压应力。

3.1.2 位移场分析

天然状态下，边坡模型在施加自重应力后会产生初始位移（见图5～7）。后面的工况中，ANSYS计算出的位移包含了初始位移，因此在绘制位移云图时从现位移中减去该初始位移，即可得到因荷载改变产生的增量位移，即真实位移。

3.2 开挖边坡应力应变分析

3.2.1 应力场分析

天然边坡坡面主应力较小，在开挖后J2h板岩倾倒变形体开挖面部位在水平方向上位移较大，开挖后的临空面应力释放，最大主应力由压应力变为拉应力，大小为0～0.29 MPa，最小主应力为0～0.21 MPa，边坡下游圆弧开挖面应力集中现象不明显，上游折线开挖面在坡脚与开挖折线处产生较大的应力集中，最大主应力大小为3.30 MPa，此开挖脚处产生较大的塑性变形（见图8和图9）。坡体内部最大主应力在F12断层附近出现应力松弛效应，在F12断层附近，断层上盘的主应力大于断层下盘的主应力。边坡水平向位移最大位置在开挖面处J2h板岩倾倒变形体F12断层上盘附近，最大主应力在断层上盘的应力值大于下盘，连续的应力等值线在此处错开。在J2h千枚状泥质板岩与T3xd火山角砾岩、碎屑岩夹凝灰岩岩层分界面处，T3xd岩层最大主应力应力值大于J2h岩层应力。

3.2.2 位移场分析

水平向的位移仍然以上部J2h板巖区域倾倒变形体最大。在坡体内部断层与岩层界面之间的岩体主应力值相较两侧岩体小，断层具有应力松弛效应，倾倒变形体与下伏基岩界面上部变形体应力值小于下部基岩。2#倾倒变形体的后部J2h千枚状泥质板岩在表层有较大的指向坡外的位移，基本与现场监测情况一致。开挖期间，边坡向临空方向卸荷回弹变形以及沿各缓倾断层滑移变形继续增大（见图10）。

4 结语

本文对云南省某开挖边坡进行了三维模拟，计算边坡天然应力，对边坡进行开挖模拟，计算开挖对边坡岩体的应力应变影响，并对比开挖后的应力应变、位移，得出开挖边坡的稳定性分析结论：

（1）针对边坡开挖的卸荷变形问题，介绍了使用有限单元法模拟边坡开挖的方法。

（2）自然边坡在人工开挖后，边坡开挖面在卸荷作用下，开挖面附近坡体原有应力场将发生变化，应力场的变化导致了在开挖面附近岩土体产生变形。变形在开挖面上最大，向坡体内部逐渐减小。在开挖过程中，开挖面附近岩体卸荷，由于坡面没有约束，应力释放明显。随着向坡体内部的深入，随着岩土体间的内摩擦力与粘聚力的作用越来越强，应力释放就越来越小，但是这种变化是有限制的，开挖边坡存在一个卸荷影响范围，当超出此范围以后开挖卸荷影响就会减小，此处边坡岩土体的应力就基本仍为原边坡的应力状态。

（3）边坡的人工开挖对开挖范围附近的应力分布会有较大影响，其程度会因开挖面的深度、开挖面的形状与开挖方式而不同。开挖后开挖面会发生卸荷，随着初始应力由坡顶向坡底逐渐增加，开挖的卸荷量也逐渐增加。在开挖的台阶处发生了较大的应力集中，会导致在开挖的过程中就有岩体的卸荷变形现象。顺坡向开挖坡脚出现向外的位移，同时竖直方向在开挖面上出现指向坡外的位移，上游侧折线开挖面在坡脚处常常会发生应力集中，具有较大的张拉应力，这些都是由于岩体倾倒蠕变所导致的，如果不及时加固会发生岩体卸荷变形。所以开挖时应避免较大折角的产生。

（4）现场监测分析表明该开挖边坡稳定，理论分析与实际结论一致，可为今后类似工程模型的建立提供依据，具有一定的理论意义及工程运用价值。

参考文献：

[1]王敏强，许 原. 有限元分析中开挖释放荷载计算的讨论[J]. 武汉大学学报（工学版），2001，34（1）：56-59.

[2]章 青.有限元分析中开挖释放荷载的正确计算[J].河海大学学报，1999，27（3）：112-115.

[3]李 围. ANSYS土木工程应用实例[M]. 北京：中国水利水电出版社，2006.

[4]王 勇，殷宗泽.有限元计算深开挖中挖方等效荷载的分析[J].河海大学学报，1998，26（5）：71-74.

收稿日期：2020-06-02