基于有限元应力的三维边坡稳定分析系统

2022-09-28韦立德刘文连陆得志李江龙

计算机仿真 2022年8期

韦立德，刘文连，陆得志，李江龙

(1. 中国科学院武汉岩土力学研究所，湖北武汉 430071；2. 中国有色金属工业昆明勘察设计研究院有限公司，云南昆明 650051)

1 引言

自然界发生的滑坡绝大多数呈三维形态，越来越多的工程实际提出了建立三维边坡稳定分析的要求。三维边坡稳定分析可以更加真实地反映边坡的实际形态，特别是当滑裂面已经确定时，使用三维分析可以恰当地考虑滑体内由于滑裂面的空间变异特征对边坡稳定安全系数的影响。边坡稳定性分析方法有多种，但应用最为广泛和有效的主要有两种，即基于刚体极限平衡理论的传传统计算方法具有简单、方便、快速的特点，目前二维刚体极限平衡法已经成熟，在边坡方面应用较为广泛，且三维刚体极限平衡法在工程上也得到了应用[1-4]。然而，传统计算方法在条间力、滑动面上假设过多，且不考虑边坡的应力-应变关系，所得到的底滑面上的力不尽符合实际，也不能求得边坡本身的变形对边坡变形及稳定性的影响。有限元法不仅满足力的平衡条件，而且还考虑了应力、变形关系，能够得到边坡在荷载作用下的应力、变形分布，所得到坡体的应力、应变与实际情况比较接近，比刚体极限平衡法更为精确合理[4，5]。有限单元法在分析边坡计算安全系数时的困难之一是如何将有限元计算的结果转化为工程上适用的安全系数。目前有两类基本的有限元方法，即建立在强度折减基础上的有限元法和建立在滑裂面应力分析基础上的有限元法。强度折减法因通过计算能同时得到边坡体安全系数及临界滑面位置，因此在工程实践中得到了较广泛的应用。然而该方法需要反复试算，在接近极限状态时，收敛很慢，尤其对于复杂的三维动力问题，其计算时间让人无法忍受。除此之外，由于影响强度折减法的因素较多，到目前为止，还没有一种被工程界和学术界都广为接受的通用标准，一定程度上限制了该方法在实际工程中的广泛应用。建立在滑裂面应力分析基础上的有限单元法，是在边坡中定义一个潜在的滑裂面，根据有限元得出的应力分布，计算滑裂面上各点应力水平，再定义安全系数；基于有限元计算的三维边坡滑裂面搜索方面有一定的研究[5-10]，但主要采用遗传算法等方法搜索且搜索变量个数比较少，当前边坡工程商业软件往往采用穷举法搜索因为普遍认为该法更可靠[11]，基于有限元计算的三维边坡滑裂面搜索方法在易行性、应用于复杂工程方面还有待于进一步地研究和完善。

本文作者曾致力于基于有限元计算的边坡三维滑裂面搜索方法及工程应用探索[12]，本文是该工作的一部分。本文建议了一个更加易行的三维边坡分析系统，总结基于有限元计算的边坡三维滑裂面搜索方法，然后用穷举算法搜索出最危险滑裂面，得出相应的安全系数作为边坡的安全度评价指标，建议滑动体的Flac3D建模与显示技术，并研制了基于该法的程序，用于边坡稳定分析。本文研究成果对三维边坡稳定性分析与治理具有借鉴意义。

2 安全系数的定义

在此边坡稳定性安全系数采用基于剪应力计算方法，即坡体沿整个滑动面的安全系数定义为

(1)

式中：σn、τ为滑动面上的法向应力和切向主应力；c为材料的黏聚力；φ为材料的内摩察角；A为滑动面的面积。

3 基于有限元软件计算应力的三维边坡极限平衡分析方法

利用商业有限元软件的强大建模能力建立边坡几何模型和有限元网格模型，利用商业有限元软件可以计算有限元网格模型节点多的优势进行三维边坡有限元计算，利用有限元计算中的有限元网格和有限元单元应力信息，滑动面采用很多三角形模拟，计算边坡稳定性安全系数、下滑力和抗滑力，为边坡稳定评估和加固设计提供依据。

3.1 三维边坡极限平衡分析系统中边坡模型构建

由边坡有限元网构成边坡几何模型，由有限元软件计算单元高斯点应力分布构成边坡应力场分布。为了更好模拟自然边坡地质构成的复杂性，边坡有限元网模型由四面体单元构成，单元高斯点应力为四面体单元的高斯点应力平均值。

3.2 滑动面参数设置和由三角形组成的滑动面构建

对三维问题，一般假设滑动面为椭球面，本文采用此假设，由三个球心坐标和三个半径6个变量(xi，yi，zi，ai，bi，ci)及滑动角θi共7个变量确定潜在椭球形滑动面，其方程如下面所示

(2)

在确定7个变量变化范围后基于穷举算法通过循环语句给出一组组7个具有定值的变量组成的数组，每个数组代表一个潜在椭球形滑动面的确定参数；确定哪个四面体与椭球面相交，基本搜索方案是遍历有限单元列表，对选定单元的所有边(线段)的两个端点坐标(x1,y1,z1)和(x2,y2,z2)带入椭球面方程，以下列不等式成立表示椭球面与线段相交，满足选定单元三条线段与椭球面相交表示椭球面与该四面体单元相交，根据交点建立三角形模拟在该选定单元范围内的潜在滑动面。

(3)

遍历所有有限单元列表，所有建立的三角形滑动面的集合就是采用的以三角形模拟的对应(xi，yi，zi，ai，bi，ci)和θi确定的滑动面，如图1所示。

图1 三维边坡滑动面和滑动角

3.3 形心所在单元搜索

为求解潜在滑动面应力，首先需要确定潜在滑动面每个三角形形心在哪个单元内。基本搜索方案是遍历有限单元列表，以确定形心是在当前单元的“内部”或“外部”。如果形心是在“内部”，单元被找到，搜索完成。如果形心是在“外部”，继续搜索列表中的下一个元素，直到找到包含该形心的单元。在此采用矢量点积法判断三角形形心是否在某个有限元单元内。对单元的平面定义指向外部的法线向量，三角形形心到该单元平面上任意一节点定义第二个向量，这两个向量的点积小于或等于0说明三角形形心在单元内部，据此找到三角形形心所在单元。

3.4 滑面三角形内的应力分析

有限元分析结果给的应力是单元高斯点应力，取高斯点应力平均值为单元平均应力。三角形滑面应力近似取为其形心所在有限元单元平均应力。

图2 应力分量图

S1=cosθ，S2=sinθ，S3=-cosθ·n1/n3-sinθ·n2/n3

(4.1)

特殊情况在n3=0时，该小滑块处于竖直方向，那么该滑块的滑动方向一定在沿z轴负向，那么该滑块的滑动方向为

S1=0，S1=0，S1=-1，s1=S1/S，s2=S2/S，

(4.2)

沿着滑动方向的切向应力为ts=txs1+tys2+tzs3。

3.5 安全系数计算

已知组成滑动面的三角形形心处的正应力tn和切应力ts之后，由摩尔-库仑准则可知，该滑块的抗滑力为fn=(c+tn·tanφ)A。而该滑块的下滑力为fs=tsA。通过叠加计算，就可以得到滑面的抗滑力和下滑力，那么，在整个滑面上的安全系数可以表示为

(5)

其中ci和φi分别为滑块三角形i形心处的凝聚力和内摩擦角，Ai为滑块三角形i的面积。

3.6 优化计算

采用“穷举法”进行优化计算，通过使用嵌套 7 个循环并依次改变循环变量xi、yi、zi、ai、bi、ci和θi的值，来搜寻边坡稳定性安全系数Fn最小值，搜寻到最小值Fn及对应变量xi、yi、zi、ai、bi、ci和θi的值，即为要求解边坡稳定性安全系数Fn和对应的滑动面。

3.7 滑动体图形显示

在获取含滑面穿过的四面体单元和滑动面内部包含的单元组成的单元集后，即可用该单元集中的每个单元的有限单元信息和包含的节点的坐标编写出Flac3D的对应该单元的块体建模语句，对该单元集中的全部单元实施从而完成该单元集Flac3D建模，进而输入Flac3D建模和利用Flac3D显示出近似最危险滑动体。

4 程序验证

本文的方法编入了软件FEMSTAB3D计算。ZHANG Xing曾发表文章，对一系列具有简单体形的椭球滑面的三维安全系数提供了解答，本文选用其中一个算例(见图3)对程序进行验证。

对此问题采用本文方法FEMSTAB3D计算分析，有限元网格见图4，只求解(xi，yi，zi，ai，bi)共5个变量解，采用步差为0.25米的穷举算法求得边坡稳定性安全系数计算结果为2.074，最危险滑动面球心坐标点对应图3坐标系的最危险椭球面滑面的球心坐标点是(35.10，28.15)，椭球面在各轴上的截距为：b=69.9，a=c=25.90，滑动体大致如图5所示。由于如果滑动体只包含滑面范围以内岩土体的显示比较麻烦，为了图显示方便这里的滑动体把滑面穿过的4面体单元全部包括到滑动体了，因此滑面显示不是曲面，通过这样大致滑动体包含滑面穿过的四面体单元和滑动面内部包含的单元。各种方法计算结果见表1，本文方法和其它方法安全系数计算结果一致，在工程上是可以接受的。

表1 各种方法计算得到的安全系数对比

图3 椭球面边坡模型

图4 有限元网格模型

图5 潜在滑动体

5 工程应用

攀钢西渣场位于金沙江东岸，在弄沟口北侧堆放。渣场边坡稳定问题较为突出，渣场边坡的稳定性对攀枝花市的安全意义重大，对渣场边坡的稳定性分析就非常必要了[12]。渣场高边坡有限元模型见图6，该模型模拟了 10种主要的组成材料，其主要参数见表2。本文关注渣场边坡右半部分(接近下游方向部分)边坡的稳定性。对此问题采用本文方法FEMSTAB3D计算分析，采用步差为0.25米的穷举算法搜索到危险滑面为对应球心(630.5，-2.0，115.5)和半径(63.5，108.0，108.0)的潜在椭圆球滑动体大致如图7所示，边坡稳定性安全系数计算结果为1.400，边坡稳定。对穿过该滑动体的一个剖面采用Slide软件极限平衡法分析，求得稳定性安全系数为1.321(简化毕肖普法)。采用开发软件进行边坡稳定性安全系数计算结果与采用二维极限平衡法计算结果相差大约6%，分析结果与一般规律一致。