锦屏一级水电站猴子崖边坡三维非连续变形数值模拟

秦永涛 宋慧 游湘 杨 涛

（1 中国水电顾问集团成都勘测设计研究院，四川 成都 610072；2 西南交通大学土木工程学院，四川 成都 610031）

1 引 言

边坡稳定性是岩土工程学中的一大类问题，国内外众多学者对此都做过研究，应运而生了许多有价值的研究成果，对此很多人做过论述[1]～[6]，为边坡工程的设计、施工起到了很好的指导作用。近几十年来,随着大型水利水电工程、公路工程、铁路遂道工程和采矿工程等的迅速发展, 岩土力学作为岩土工程基础学科得到了较快的发展。特大型的工程实践为岩土力学的发展赋予了巨大的推动力,同时也产生了许多复杂的岩土力学新课题。特别是西部大开发以来，西南地区水利水电工程、公路铁路基础设施建设中遇到的大量复杂地质情况、高地应力、高陡岩质边坡和地下洞室群问题，使得对岩土工程的认识有了更深的要求，出现了更多的工程技术难题，更多的现有理论需要进一步完善。

解决岩土力学问题的方法主要有实验方法、理论分析方法和数值模拟方法三大类。这三类方法相辅相成,互为补充。与其它两种方法相比，数值模拟是解决岩土工程问题较有效的手段,已被学术界和工程界广泛接受,作为一种力学状态的分析工具,它越来越多的应用于岩土体的稳定性分析、岩土工程设计和岩土工程基本问题分析中。但是前面两种方法是数值分析的前提和基础，数值模拟是对实验结果和理论分析成果的佐证。

几十年以来,岩土力学数值分析方法得到了迅速发展, 出现了有限单元法、离散元法、有限差分法、边界元法、无限元法、拉格朗日元法、非连续变形分析方法、无单元法、流形元法及其耦合的数值计算方法和以数值模拟为主的渐进破坏模型等[7]～[14]数值分析技术。其中基于连续介质大变形分析的拉格朗日元法在实际工程中也得到了较好的应用。拉格朗日元法避免了有限元法进行大型矩阵的复杂计算，可以同时考虑岩土体的材料非线性和几何非线性,并能跟踪物体变形的全过程，适于分析岩土力学中的大变形问题。本文运用连续介质大变形分析理论计算了锦屏一级水电站右岸下游河道雾化区猴子崖边坡(简称“猴子坡”)在天然状况和泄洪雾化雨作用下的应力应变变化及位移变化趋势，评价了边坡的稳定性。

2 连续介质大变形理论

三维快速拉格朗日（3D-FLAC）分析是一种基于显式有限差分法的数值分析方法,是非连续变形分析方法中的一种。该法将计算区域划分为若干六面体单元，单元网格可以随着材料的变形而变形，即所谓的拉格朗日算法。这种算法基本原理类同于离散单元法，但它却能像有限元那样适用于多种材料模式与边界条件的非规则区域的连续问题求解。还能针对不同的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为。可以准确地模拟材料的屈服、塑性流动、软化直至大变形,尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点[15]。

2.1 空间导数的有限差分近似

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

其中V为四面体的体积，S为四面体的外表面，nj为外表面的单位法向量分量。

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

式中上标l表示节点l的变量，（l）表示面l的变量。

2.2 运动方程

快速拉格朗日分析以节点为计算对象，在时域内求解。节点运动方程为

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

将（3）式左端用中心差分来近似，则可得

2.3 应变、应力及节点不平衡力

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

有了应变增量，即可由本构方程求出应力增量，进而得到总应力。

2.4 阻尼力

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

阻尼力为

式中 α为阻尼系数，其默认值为0.8，而

3 工程概况

锦屏一级水电站位于四川省凉山彝族自治州盐源县和木里县境内，是雅砻江干流水能资源最富集的中、下游河段五个梯级水电开发的第一级。坝址位于普斯罗沟与手爬沟之间1.5km长的河段上，河流流向N25°E，河道顺直而狭窄。电站以发电为主，兼有防洪、拦沙等作用。坝址控制流域面积10256km2，多年平均流量1200m3/s。水库正常蓄水位1880m，拦河大坝为混凝土双曲拱坝，坝高305m，水库库容77.6亿m3，为年调节水库，电站装机容量3600MW，年发电量166.2亿kW·h。

猴子坡为坝区右岸Ⅳ线下游与Ⅲ线之间一三面临空的凸出山梁。地貌上该斜坡呈SEE向延伸的山脊状，脊状斜坡的SW侧和NE侧均呈陡坎，下部较窄。山脊走向与河流近于垂直，临江坡度50°～80°。坡体由杂谷脑组第二段第6、7层大理岩组成，第6层为薄～中厚层状条带状大理岩，局部夹绿片岩夹层，第7层为厚层条纹状大理岩；岩层产状N40°～55°E，NW∠35°～40°。如图1所示。

图1 猴子坡数值模拟区域

该边坡距大坝约600～780m，根据泄洪雾化模型试验报告，猴子崖边坡处于右岸泄洪雾化区边坡范围内，加上有多种不利结构面组合，需高度重视在泄洪雾化降雨作用下的边坡稳定性分析。

4 三维数值模拟分析

4.1 三维数字模型的建立

根据右岸猴子岩山梁区地形测绘数据和地质资料建立三维几何模型，如图2。模型采用局部坐标系（X、Y、Z），其X轴正方向为NE28°，指向河谷下游，Y轴正方向为NW62°，指向左岸山里，Z轴正方向铅直向上。几何模型建立范围为：上下游X方向，以猴子岩山梁山脊线为基线分别向河流上下游方向各延伸330m和500m，总长830m；垂直河流方向，从河谷向山内侧延伸550m；高程方向，从1480m高程一直到自然坡面，最大高程至2345m。

三维几何模型中真实再现了计算区域复杂的地形地貌、控制性结构面。对猴子坡山体而言，控制性结构面包括：f7、f28、g37-1～g37-5、近SN向结构面。f7、f28的产状及位置均确定，如图2所示。对于层间挤压错动带而言，该区发育g37-1～g37-5等多条层间挤压带，而计算模型中不可能全部真实考虑，只能进行一定程度的概化。根据三维刚体极限平衡分析结果，g37-2作为底滑面的潜在失稳块体稳定性最低、块体体积最大，最可能失稳，是该区的控制性结构面，因此有限元计算中采用g37-2的产状建立层间挤压错动带，而将层间挤压带的厚度设为10m来综合反映其它各条错动带对边坡稳定性的影响。近SN向结构面可能构成猴子坡的下游侧边界，其产状对边坡失稳范围的确定十分重要，计算中采用遍布节理模型来模拟，节理模型的产状采用近SN向结构面的产状。

图2 猴子坡数值计算模型

在几何模型的基础上，建立有限元计算模型。由于猴子坡地形复杂、结构面众多、产状变化大，为较好的模拟这些地质因素，模型全部采用四面体单元进行离散。计算域共剖分了762612个单元，134073个节点。剖分后的有限元网格图如图3所示。对于猴子坡而言，f7、f28、g37～2等结构面是特别关心的，因此进行了重点离散，共离散为62906个单元，如图4所示。

图3 模型网格划分

图4 结构面网格

4.2 计算参数和边界条件

根据计算域的地质情况，按岩体分级分区选择材料，岩体采用遍布节理模型模拟，节理的产状设为近SN向结构面的产状，f7、f28断层、层间挤压错动带g37～2等采用莫尔－库仑材料模拟。材料参数见表1。

表1 计算中采用的岩体力学参数

近SN向结构面的抗剪强度参数按照20%取B1类指标，80%取IV1类岩体指标综合加权确定。

根据模型所处的地形地貌条件、河谷对称性及边坡载荷方向，模型边界顺河向边界(X向)、横河向边界(Y向)和底部边界(Z向)分别取法向支座。

4.3 计算工况

计算分别考虑2种工况：

工况1：天然状态下，猴子坡山梁的应力状态及其稳定性分析；

工况2：泄洪雾化状态下，猴子坡山梁的应力状态及其稳定性分析；

其中，对工况2计算中，考虑最不利的情况，高程1830m以下的岩体均达到饱和状态，岩体容重取饱和容重(γ＝29kN/m3)，在上游侧边界(－X)和山内侧(－Y)1830m高程以下施加随高度变化(梯度为10kN/m)的三角形荷载，以模拟雾化雨形成的水力边界，如图5所示。

图5 雾化雨形成的水力边界示意图

4.4 计算成果分析

(1) 天然状态的计算成果

坡体在自重应力场作用下的应力场分布如图6～图11所示。其中，图6为猴子岩附近区域的最大主应力等值线（拉为正），图7为坡体表面的拉应力区，为分析f28、f7及g37～2等结构面上的应力特征，图8示出了边界面上的最大主应力云图。图9为坡体表面的最小主应力等值线。图10示出了坡面上的塑性屈服区。

由图分析可知，猴子坡潜在失稳块体范围明显。在自重作用下，由于猴子岩山梁为一突出山梁，其上游和低高程（河谷处）由f7和g37～2所围限。在f7和g37～2交界部位，各结构面在坡表面均表现为明显的拉应力区，即使是从最小主应力图来看，这部分坡体表面也呈现出极小的压应力状态，表明猴子岩失稳块体可能沿此界限失稳下滑。从结构面的最大主应力云图来看，f7和f28浅层均出现拉应力区，其下部和g37～2上游侧则表现为压应力，但应力水平较低，而g37～2下游侧压应力水平较高，等值线呈由河谷向山里的展布状态，也揭示了其上块体的失稳趋势。

从塑性区分布来看，猴子坡低高程潜在失稳块体存在失稳趋势，安全裕度不高，但失稳块体还没有完全形成。坡面屈服区仅出现在f7、f28和g37～2的高程局部处，并没有贯穿整个g37～2结构面。结构面上的屈服区分布更能体现这一特征，f7和f28的绝大部分均出现了塑性屈服，表明边坡块体有失稳滑移趋势，后缘边界明显。但g37～2仅坡面附近出现塑性屈服，其分布范围并不大，则又表明块体还没有完全形成，还不可能沿此底滑面出现整体失稳滑移。由于两侧边界的约束作用以及坡脚岩体的抗滑作用，边坡仍能保持整体稳定。但是由于猴子岩山梁底滑面前部（层间挤压错动带）和上游侧边界（f7断层）基本上已经屈服，一旦山梁后缘（近SN向结构面）形成拉裂缝，容易发生顺层的滑移拉裂破坏，安全储备不高。

图6 最大主应力等值线

图7 拉应力区

图8 边界面上的最大主应力

图9 最小主应力等值线

(2) 雾化雨作用下的计算成果

图11～图13所示为雾化雨作用下的坡体应力，图14示出了坡体的位移场变化，图15为坡体结构面上的塑性屈服区。

在雾化雨作用下，坡体应力场进一步恶化，主要表现在拉应力的量值增大，作用范围集中而贯通，块体失稳边界逐步形成。猴子坡山梁上游侧（f7通过区域）拉应力分布明显，整个区域均呈现出受拉状态，沿g37～2向下游，拉应力区集中而呈条带状（如图12），表明g37～2作为上游侧边界、g37～2作为底滑面的失稳范围十分明确，块体已经沿此出现了拉裂破坏。而从河谷应力来看（如图13），河谷出现了极大的压应力集中，势必造成岩体屈服，因此，河谷侧的压应力区也揭示了猴子坡失稳块体的范围所在。从结构面的最大主应力云图来看，除f7和f28浅层均出现拉应力区外，较之天然状态，结构面上的压应力水平大为降低， f7和f28下部和g37～2上游侧压应力水平大为降低，几乎接近0应力状态，考虑到这些结构面位于块体的上游侧，为失稳块体的后缘位置，这些部位出现0应力状态表明块体具有明显的失稳下滑趋势。

在雾化雨作用下，坡体发生了明显变形，变形区域集中在猴子坡前部，揭示了潜在失稳区域。如图14所示，猴子岩前部发生了明显的位移变化，最大位移达到0.055m，向山里逐渐递减，位移等值线方向与近SN向结构面的走向基本一致，表明近SN向的结构面是失稳块体的下游边界面。若以位移大于0.01m的范围作为失稳块体范围，则限定的失稳块体范围如图20所示。可见失稳块体的边界为f7、g37～2及近SN向结构面，而f28由于其出露高程较高，没有形成失稳块体的边界。

图15示出了坡体结构面上的塑性区分布，可见，结构面上的塑性区大为发展，g37～2以上的坡体也沿节理方向发生了塑性流动，由f7、g37～2和近SN向结构面所围限的岩体全部发生了塑性流动，在临空面具备的情况下，块体极有可能出现失稳下滑。

图12 拉应力区

图13 最小主应力等值线

图14 位移大于0.01m的区域

图15 结构面上的屈服区

5 结 论

猴子坡岩体可能的失稳破坏模式为：f7断层、层间挤压错动带、近SN向结构面形成楔形体组合，有限元计算中，考虑了这些结构面的影响，f7断层和层间挤压错动带通过建立单独的岩体材料反映，近SN向结构面为随机分布，采用遍布节理模型加以反映，重点计算分析了天然边坡和雾化雨边坡的稳定性。

由于猴子坡的特殊坡体结构，应力场特征揭示，在天然边坡条件下即具有较明显的失稳块体范围，其范围与地质判断一致。塑性区分布揭示，猴子坡低高程潜在失稳块体存在失稳趋势，安全裕度不高，但失稳块体还没有完全形成。

受雾化雨作用，坡体应力场进一步恶化，应力场特征揭示，失稳块体边界已经逐步形成，通过位移场的判断，可进一步明确失稳块体的范围，该范围内的岩体变形在进一步发展中。结构面和坡体节理方向塑性区大为扩展，在临空面具备的情况下，块体具有明显的失稳下滑趋势。

数值计算表明，猴子坡在雾化雨作用下可能会出现局部失稳，需要对该部分岩体加强支护，做好防排水措施。

[1] 陶振宇. 试论岩石力学的最新进展[J]. 力学进展,1992 ,22 (2) :161～172.

[2] 常士骠.我国岩土工程回顾与展望[J].水文地质工程地质,1993,20(1): 23～26.

[3] 哈秋聆.三峡工程中的若干力学问题[J].力学进展,1994,24(4): 433～439.

[4] 杨志法.关于岩石力学当前发展战略的一些看法[J].岩土工程学报, 1994 ,16 (1) :112～113.

[5] 谢和平,刘夕才,王金安. 关于21 世纪岩石力学发展战略的思考[J].岩土工程学报,1996 ,18 (4) :98～102.

[6] 傅冰骏. 国际岩石力学发展动向[J].岩石力学与工程学报,1997,16 (2) :195～196.

[7] Zienkiewicz O C. The finite element method.London: Mc-Graw-Hill, 1997.

[8] 王泳嘉,邢纪波. 离散单元法同拉格朗日元法及其在岩土力学中的应用[J].岩土力学, 1995, 16 (2):10～15.

[9] 王泳嘉.边界元法在岩石力学中的应用[J].岩石力学与工程学报,1986,5 (2) :205～222.

[10] 应隆安著.无限元方法[M]. 北京:北京大学出版社,1992.

[11] Shi G H. Numerical manifold method[M].Proceedings of the First International Conference on Analysis of Discontinuous Deformation, Chungli,Taiwan,1995:187～222.

[12] Shi G H. Discontinuous deformation analysis: A new numerical model for the static and dynamics of block systems. A dissertation submitted in partial satisfaction of the requirements for the degree of Doctor of Philosophy, Berkeley: University of California, 1988.

[13] Margulies M. Combination of the boundary element and finite element methods. Proceeding in Boundary Element Methods 1981(1): 258～288.

[14] Chang C C. Computational simulation of progressive fracture in fibre composites. Journal of Composite Materials, 1987 (21): 834～885.

[15] FLAC3D Version2.0 Online Manual[Z].Minneapolis: Itasca Consulting Group, Inc, 1997.