应用FLAC3D模拟分析输变电工程中的边坡稳定

郑晏超 索亮 万正东

摘要：近年来，我国电网工程发展迅速，输变电工程设计中应用FLAC3D（三维快速拉格朗日分析程序）模拟分析杆塔定位和变电站选址过程中遇到的边坡稳定问题，可帮助设计人员对边坡的稳定性作出快速有效的判断，以便及时避让或提前做好防治预案。

关键词：FLAC3D；输电线路；边坡稳定

1 引言

连续介质快速拉格朗日差分法(Fast Lagrangian Analysis of Continua，简写FLAC)是近年来发展起来的一种新型数值分析方法，它与通常的有限元法和边界元法的不同之处在于：前者是显式的方法，而后者则是隐式的。显式差分法求解时未知数集中在方程式的一边，无须形成刚度矩阵，不用求解大型联立方程，因而占用内存少，便于微机求解较大的工程问题。另外，由于该方法采用的是随流观察法，即研究每个流体质点随时间变化的情况，所以适合于解决非线性大变形问题，这一点恰好可以满足岩土工程中力学分析的需要。

2 FLAC3D的基本原理简介

FLAC3D是美国Itasca公司开发的三维快速拉格朗日分析程序，该程序使用了如下三种计算方法[1]：(1)离散模型方法。 (2)有限差分方法。 (3)动态松弛方法。

(1) 空间导数的有限差分近似

在快速拉格朗日分析中采用混合离散方法[2]，即将区域离散为常应变六面体单元的集合体，而在计算过程中，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体，应力、应变、节点不平衡力等变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均，六面体内四面体的应力应变第一不变量由该六面体内所有四面体的体积加权平均得到，然后重新得到四面体的应力、应变。如图4.2所示一四面体，节点编号为1到4，第n面表示与节点n相对的面，设其内任一点的速率分量为vi，则可由高斯公式得：

（1）

图1四面体

式中：V为四面体的体积；S为四面体的外表面；nj为外表面的单位法向向量分量。

对于常应变单元，vi为线性分布，nj在每个面上为常量，由式(1)可得：

（2）

式中：上标l表示节点l的变量；(l)表示面l的变量。

(2) 运动方程

快速拉格朗日分析以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可表示为如下形式：

（3）

式中：Fil(t)为在t时刻l节点的在i方向的不平衡力分量，可由虚功原理导出。ml为l节点的集中质量，在分析静态问题时，采用虚拟质量以保证数值稳定，而在分析动态问题时则采用实际的集中质量。

将上式左端用中心差分来近拟，得到：

（4）

(3) 应变、应力及节点不平衡力

快速拉格朗日分析由速率来求某一时步的单元应变增量

（5）

式中导数可由式(2)近似。

有了应变增量，即可由本构方程求出应力增量，各时步的应力增量叠加即可得到总应力，在大变形情况下，还需根据本时步单元的转角对本时步前的总应力进行旋转修正。然后即可由虚功原理求出下一时步的节点不平衡力，进入下一时步的计算，其具体公式这里不再赘述。

(4) 阻尼力

对于静态问题，在式(3)的不平衡力中加入了非粘性阻尼，以使系统的振动逐渐衰减直至达到平衡状态(即不平衡力接近零)。此时式(3)变为：

（6）

阻尼力为：

(7)

式中：为阻尼系数，其默认值为0.8.

（8）

(5) 计算循环

以上介绍了三维快速拉格朗日分析（FLAC3D）的基本原理[3]，可以看出其计算循环如图2所示。

图2计算循环

3实例模型的创建

以某输电线路工程拟选杆塔塔位附近边坡体为例创建FLAC模型如图3所示：

图3边坡模型的建立

由于土坡垂直高度不大（8m），故本模型尺寸采用1：1比例。材料本构模型为摩尔-库仑模型，选用大变形理论计算。

边界条件：除了上表面及斜坡面不加约束外，其余表面均施加固定约束。

荷载条件：本例仅考虑土坡体在天然与浸水饱和两种状态下，在重力场中的变形模拟，暂不考虑施加杆塔荷载的情况。

计算参数：弹性模量E为1.8e6Pa，泊松比λ取0.3，抗拉强度为1.78e3Pa，本例中边坡体全部处于黄土层，土的重度以及抗剪强度指标选用表1中的数值。

表1 边坡土抗剪强度指标

4 模拟结果分析

(1) 塑性区分析

塑性区图能反映出边坡最易破坏的的位置，本例中边坡在天然状态和浸水饱和状态两种工况下的塑性区图见图4和图5。

图4 天然状态下塑性区图

图5 浸水饱和状态下塑性区图

由图4可知，在天然状态下，滑坡体的塑性区域仅出现于坡体前缘及坡脚处，后缘也有一小部分单元出现塑性破坏。其中，滑体后缘及前缘上部塑性区以拉伸破坏为主，坡脚处主要为剪切破坏。滑坡体的后緣和坡脚的塑性破坏区远没有贯通，因此，在天然状态下，边坡处于稳定状态。随着地表水及雨水的浸入，边坡体含水量增加，其塑性区也随之扩大。从图5可以看出，土坡处于浸水饱和状态时，塑性区面积几乎充满了整个土坡模拟体，且滑体大部以张拉破坏为主，张拉区域发生在滑坡体表层易使滑坡体顶部形成裂隙，导致坡体后缘水的涌入，最终导致滑坡失稳。此时，由于地下水和表面降水入渗以及人为因素的影响，边坡体塑性破坏区完全贯通，坡体将发生滑动破坏。

(2) 位移分析

由位移等值线云图可以看出边坡产生滑动后的位移矢量及坡体变形趋势。由图6和图7可分析边坡的变形特征如下：

随着土坡由天然状态转变为浸水饱和状态，滑坡坡角的粘聚力变低，各方面的力学性变差，滑坡坡面在z方向（垂直方向）的位移逐渐增大。图6与图7中深色

图6 天然状态下z方向位移云图

图7 浸水饱和状态下z方向位移云图

部分为垂直位移最大部分，位移由5.76cm（天然状态）增至75.6cm（浸水饱和状态,坡体已破坏）。此外，由图7能清晰地看出边坡最可能滑动面的位置。

(3) 剪应变增量分析

岩土体失稳，大都是沿剪应变最大的部位发生，大量工程实例分析结果也证明了这一点。判断边坡体的(潜在)滑动面(带)，可根据其剪应变增量的大小来判断：剪应变增量较大(绝对值)的部位，则为其(潜在)滑动面(带)，变形破坏也多沿此处发生。我们可以通过FLAC3D计算出的剪应变增量来找出滑坡体内的最易剪切破坏的区域，亦即边坡实际滑动面的位置所在。

土坡在天然与浸水饱和两种状态下的剪应变增量云图及剪应变增量集中带分布见图8和图9。由于浸水饱和状态下滑坡体最终破坏，故此处选取其破坏前的某一临近状态时刻的剪应变增量云图。

图8天然状态下剪应变增量等值线图

图9浸水饱和状态下剪应变增量等值线图

由图8可知，天然状态时，坡体最大剪应变增量仅为2.34e-02，剪应变集中带没有贯通，边坡体处于稳定状态；当坡体含水量增大时，剪应变增量的量值、范围也将随之增大，由图9可以看出，当土坡体处于浸水饱和状态时，图中所示时刻的剪应变增量最大值已达3.68e-01，剪应变增量带已经发展到即将贯穿整个边坡体，此时的边坡或已发生失稳破坏。

剪应变增量集中带的发展趋势反映了滑坡体在水的作用下土坡破坏的渐进过程，剪应变增量集中带的分布位置则显示了滑坡可能存在的潜在滑动面，由图9我们可以看出本例中边坡体潜在滑动面的位置。

除从上述几个方面对边坡稳定进行模拟分析外，FLAC3D还提供如滑移速率、应力变化等角度的模拟分析， 不管从哪个角度对边坡体进行分析，最终的结论都是一致的。本例中，我们便可得知，该边坡在天然状态下基本是稳定的，但当遭遇暴雨、地下水位升高、地表水侵入以及人类活动等因素影响时，边坡体就会在浸水情况下处于失稳状态，对于输变电工程而言，此类边坡是不易立塔或作为变电站选址的。

* 结语

本文通过介绍应用FLAC3D（三维快速拉格朗日分析程序）模拟分析杆塔定位和变电站选址过程中遇到的边坡稳定问题，从而可帮助设计人员对边坡的稳定性作出快速有效的判断，以便及时避让或提前做好防治预案。

参考文献(References)：

[1] Itasca Consulting Group, Inc.USA.FLAC-3D[M]。Fast Lagrangian Analysis of Continua in 3 Dimensions, Version 2.0, User's Manual。

[2] Marti J P Cundall. Mixed Discretization Procedure for Accurate Modelling of Plastic Collapse[J]。Int.J.Num.& Analy.Ethods in Geomech.,1982,(6).

[3] 寇晓东，周维垣，杨若琼.三维快速拉格朗日法进行小湾拱坝稳定分析[J].水利学报，2000，9：8-14。

注：文章内所有公式及图表请以PDF形式查看。