有限差分法在高边坡稳定性的应用研究
【摘 要】介绍了有限差分法的基本原理与求解方法，将有限差分算法计算得到的边坡坡体应力分布，与极限平衡条分法结合，得到岩质高边坡的稳定性分析结果，结合广东天然气接收站高边坡计算出坡体内应力场和位移场，结合极限条分法进行稳定性的分析，得到比较符合实际的结果。
【关键词】高边坡；稳定性研究；有限差分法；极限平衡条分法
Application study for the high slope stability in Finite Difference Method
Zhang Ling-shuai,Gao Jian-kang,Su Ying-an
(Shenzhen Gongkan Geotechnical Engineering Co., Ltd Shenzhen Guangdong 518026)
【Abstract】The fundamental principle and method of solution of Finite Difference Method is introduced in this paper. We can get the stability analysis result of high slope if we combine the slope stress on the Finite Difference Method calculation and the limit equilibrium slice method. This paper take the high slope of natural gas station in Guangdong for example, calculating slope stress and displacement and analyzing slope stability with the limit equilibrium slice. The calculation results meet the actual situation more suitably.
【Key words】High slope;Stability study;The Finite Difference Method;The Limit Equilibrium Slice Method
1.前言
高边坡稳定性的分析历来是岩土工程界关心的重点和难点，关于这一问题的研究已有百年的历史，目前已经从均质弹性、弹塑性理论到各向异性的断裂力学、流变力学等多种力学模型，对高边坡的认识深度也从表面宏观组合刚体到了解岩体内部应力场的地质过程机制分析，认为边坡稳定性是一个复杂系统之间的平衡问题，是具有时空效应的系统问题［1］。在此基础上，产生大量的岩质边坡稳定性方法，包括数值计算法、可靠度分析法、系统非线性动力学、突变理论、模糊理论、灰色系统、神经网络、遗传算法等。但是目前还没有一种方法达到完满解决工程实际问题，基于可靠度的分析方法无法对边坡稳定性量化，目前仅有工程地质类比法、极限平衡法和数值计算法大量应用于实际工程［2］［3］。而数值计算法能够较好分析出岩体内部应力机理，是目前能够实际应用的最为先进的方法。而众多的数值计算法中，显式有限差分法能够在较少内存情况下，快速地解决边坡岩体的不连续性和大变形问题，因此受到广泛的关注［4］。
2.有限差分法原理
2.1 基本原理。
有限差分法是国外Cundall等人提出并发展的数值计算方法，该法是把连续的定解区域用有限个网格节点构成的网络来代替，把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似，把原方程和定解条件中的微商用差商来近似，于是原微分方程就近似转化为代数方程组，即有限差分方程组，解此方程组便可以得到原问题在离散点上的近似解。然后再利用插值的方法便可以从离散解得到定解问题在整个区域上的近似解［5］［6］。
对于岩土工程问题，可采用有限差分方法来求解运动方程达到计算岩土体内应力场和应变场。求解时，将计算区域划分为若干相连续单元网格，单元之间采用节点相连，各单元与网格点应满足下列平衡方程：
σij ——应力张量的分量； ρ——物体的密度； g——重力加速度分量；——虚加速度。
当处于静力平衡状态， =0，（1）式左侧等于0，如果（1）式左侧不等于0，系统会出现不平衡力，引起网格节点及单元的系统运动，及≠0，从而得到相应的动力方程：
m +c +ku =0（2）
和u分别为虚速度和虚位移，c为虚阻尼，k为弹簧刚度。通过上述思路，将静力问题转化为拟动力问题，使（2）式所反映的动力响应逐渐随虚拟时间衰减，最终达到系统平衡，当处于平衡状态时，=0，便可以得到相应于（1）式的真实解。
2.2 差分公式。
由于矩形网格的差分格式对于复杂的区域的几何逼近误差大，不能局部任意调整网格，不易处理法向导数边界条件等缺点，因此对于含较多复杂边界的岩土工程，采用三角网的差分格式。该法具有网格灵活，法向导数边界条件容易处理，且能保持质量守恒等优点［7］。
根据Winkins提出的有限差分法，将每个四边形单元分为两个常应变三角形单元，按照高斯散度定理的一般形式可以推导出：
ni ——曲面s的单位法向量；f ——标量、向量或张量；xi ——坐标向量；ds——增量弧长。
在每一个多边形拟合的任意形状面积A上，对于每一无限小的三角形面积 dA上，有
〈f〉 ——第i边上的平均值； △s——三角形第i边的长度。
如果采用每边速度矢量的均值代替 f，每个区域内应变速率 ij用节点速度 表示，根据（4）式有：
（5）和（6）式为节点速度表示的应变速率的差分方程，其中a和b为三角形单元边界两个连续点（见图1）。
2.3 运动方程
如果以节点为计算对象，将力和质量均集中在节点上，三角形节点应力为角点两边作用力的和，表示如下：
每个四边形单元包含两组三角形，每组又有两个三角形，在每组三角形内，对三角形交汇节点力求和，然后对两组三角形的力求平均值，即可得到作用于四边形上的节点力。
网格节点i的力 ，包括外荷载及自重荷载。当物体处于平衡状态，节点i上的合力为0，如果无法到达平衡状态，节点i在时域内将存在加速度，运动方程可表示为如下形式：
式中F i t为t时刻i节点的在运动方向的合力，可由虚功原理导出，mi 为i节点的集中质量，在分析静态问题时，采用虚拟质量以保证数值稳定。将式（8）左端用中心差分来近似，则可得到
对于大变形问题，可以对上式进行积分来确定不断变化的网格点的新坐标：
式（9）和（10）都是基于时间步中心建立的，可以看出，对于中心差分方程，一阶误差将会消失，再对半个时步内的位移和外力进行转换，就可以确定节点i处的速度。
该法在尽可能少的虚拟时间内，划分多个时步，对单元应力和网格节点速率循环往复计算，直至不平衡力消失，位移与应力分别趋向于常数。
2.4 应力、应变求解。
在求得节点速率后，可由速率求某一时步的单元应变增量，如下式：
根据广义的虎克定律，各向同性材料的本构方程为：
式中，λ ，μ为拉梅常数，εii= ε11+ε22+ε33为体积应变，单位张量为 δij =1（ i=j）或δij =0 （i≠j ）。
因此，单元的应力增量可以表达为：
上式中，E为弹性模量，v 为泊松比， I1为应力第一不变量。利用单元应力增量，对各时步的应力增量叠加即可得到单元总应力。
假定某一时刻各个节点的速度为已知，则根据高斯定理可求得单元的应变率，进而根据材料的本构定律可求得单元的新应力。
2.5 计算步骤。
对于二维有限差分法计算边坡稳定性，将模型对象用四边形单元划分成有限差分网格，将应力和密度定义在网格中，速度定义在结点上。在某个节点施加荷载后，该节点的运动方程可以写成时间步长的有限差分形式。在某一小时间步长内，作用于该点的荷载只对相邻的若干节点产生影响。在每一时步如果节点合力不等于零，出现不平衡状态，结点就要运动。于是在合力的作用下，结点产生加速度，在每一个时步中根据单元节点的速度变化及时间步长，可以求出单元之间的相对位移，从而求出单元应变，然后再根据材料本构模型求出单元应力。计算时可以通过计算单元之间的不平衡力，将上一步得出的不平衡力重新分布到各节点，再进行下一步的叠代，直至整个计算体系达到平衡状态。
3. 有限差分法的结果应用
目前的极限平衡条分法计算的边坡安全系数是在高度简化的情况下，在无法得知坡体内应力分布的基础上求得，因此得到的滑动面及安全系数结果比较粗糙。当边坡高度过大，或内部结构面分布过于复杂，计算出的安全系数要么偏大，要么偏小。而有限差分法可以明确坡体内复杂的应力场及位移场，了解边坡逐步破坏机理，跟踪坡体内塑性区发展情况。如果在采用有限差分法求出坡体应力分布的基础上，与极限平衡条分法结合，搜索最危险滑动面，便可以得到较为满意的结果。目前极限平衡条分法的商业软件较多，和专门的有限差分程序结合，便可以解决高边坡稳定判据的问题。
4. 工程应用
4.1 工程概况
广东天然气接收站场地边坡高度53米，长度250米，坡底20m外是两个大型液化天然气储气罐。受附近区域性断裂带的影响，场区地层岩性较复杂。坡体主要岩性为（1）泥盆系石英砂岩，褐黄色，局部为钙质砂岩，裂面铁锰质浸染；（2）碎裂石英砂岩，裂隙面可见绿泥石，岩芯碎块状；（3）构造角砾岩，具破碎角砾结构，砂质或硅质胶结，岩芯柱状。上述岩体均较为破碎，岩层接触面复杂。边坡开挖分为三级放坡，第一级和第二级高度均为15m，坡率1:0.2，第三级为23m高，坡率1:1。由于开挖后较长时间未支护处理，坡面多处崩塌，滑坡。
4.2 计算模型。
对边坡进行二维模拟，分析采用Mohr-Coulomb模型，考虑到边坡开挖后长期裸露，浅表层构造应力已经基本松弛殆尽，因此模型不考虑水平构造应力的作用，只考虑自重应力的作用。本次计算模型岩体物理力学参数如表1所示。
边坡开挖后的典型计算模型如图3所示，根据工程地质条件，利用程序剖分的单元网格，定义不同岩层岩性参数，施加边界条件，生成初始应力场，然后进行逐步分层开挖模拟，得到开挖后的坡体应力场和位移场。
4.3 计算结果。
根据开挖后的主应力矢量图（图4），可以看到，由于泥化夹层的影响，在岩层分界面和软弱夹层附近主应力方向偏转，最大主应力为59.9KPa，最小主应力为-1256KPa，且界面上的应力变化梯度明显高于岩层内部，表明顺坡的岩层倾向对于坡体稳定性有着很大的影响。
分析边坡位移场矢量图（图5），边坡的累计位移位于坡面的强风化岩层或泥化夹层附近，位移变形达到3.8cm，表明开挖后，坡体沿泥化夹层有向坡外滑动趋势。
4.4 边坡稳定性计算结果
将有限差分法得到的应力场导入边坡稳定分析程序，基于坡体实际应力分布采用Morgenstern-Price条分法计算。本次考虑降雨和地震组合工况下的稳定性，水位线设于坡体高度2/3处，地震水平力加速度系数取0.10，计算得到的安全系数为0.679（见图6）。
分析结果表明，坡体不稳定，需要加固防止滑坡失稳，滑动面沿下部软弱夹层剪出，可见滑动面位置较符合实际，如单纯采用极限条分法进行岩质边坡稳定性计算，由于高边坡和结构面组合的坡体应力极为复杂等原因，计算出的边坡安全系数和潜在滑动面与实际情况相差甚远。而基于有线差分法计算结果和极限条分法结合可以考虑实际坡体应力和软弱夹层影响，得到的结果符合实际。根据分析的结果进行边坡加固设计，经过竣工后近两年的跟踪监测，表明采用有限差分法计算的结果应用于坡体稳定性是可行的。
5. 结论
5.1 目前求解岩体边坡的数值方法中，显式的有限差分法把静力问题模拟为动力问题来解，在求解过程中更具简单和快速的特点，尤其适用于高边坡、大变形的非线性材料连续结构体，是比较先进的数值计算方法。
5.2 有限差分法较有限元能更好地考虑岩土的不连续性和大变形性，求解速度快。对具有复杂地貌、地质情况的边坡进行模拟计算，克服了单纯极限平衡分析方法的不足。
5.3 通过广东天然气接收站边坡稳定性分析表明，有限差分法计算结果参与边坡稳定性分析，滑面的形状及位置符合实际情况，结合后续监测表明该方法是合理可行的。
参考文献
［1］ 赵尚毅，郑颖人，时卫民，等. 用有限元强度折减法求边坡稳定安全系数［J］. 岩土工程学报，2002，24（3）：343 ～ 346.
［2］ 连镇营，韩国城，孔宪京. 强度折减有限元法研究开挖边坡的稳定性［J］. 岩土工程学报，2001，23（4）: 407 ～ 411.
［3］ 孙玉科. 中国露天矿边坡稳定性研究［M］. 北京:科学技术出版社，1999.
［4］ 刘娉慧，房后国，黄志全，等. FLAC强度折减法在边坡稳定性分析中的应用［J］. 华北水利水电学院学报，2007，28 （5）：52 ～ 58.
［5］ 马萃林，朱明，王建华，等. 基于有限差分法的边坡稳定性研究［J］，金属矿山，2009，399（9）：173～175.
［6］ 盛谦，任放，等. 岩体工程开挖的显式有限差分法模拟［J］，长江科学院院报，1993，10（2）：34～38.
［7］ 刘波，韩彦辉. FLAC原理、实例与应用指南［M］. 北京：人民交通出版社，2005.
注：本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文