数值模拟作业

中 国 地 质 大 学研究生课程论文课程名称 数值模拟技术 教师姓名 # 研究生姓名 # 研究生学号 研究生专业 地质工程 所在院系 工程学院 类 别: 学术硕士 日 期: 2014年11月 评 语对课程论文的评语:平时成绩：课程论文成绩：总 成 绩：评阅人签名：注：1、无评阅人签名成绩无效；2、必须用钢笔或圆珠笔批阅，用铅笔阅卷无效；3、如有平时成绩，必须在上面评分表中标出，并计算入总成绩。第一部分：数值模拟技术研究文献综述一.前言（理解、体会）近年来，随着中国大规模地进行“西部大开发”和“南水北调”等巨型工程，越来越多的岩土工程难题摆在我们面前，单纯依靠经验显然已不能有效指导工程问题的解决，迫切需要更强有力的分析手段来进行这些问题的研究和分析。自R.W. Clough 上世纪60年代末首次将有限元引入某土石坝的稳定性分析以来，数值模拟技术在岩土工程领域取得了巨大的进步，并成功解决了许多重大工程问题。特别是个人电脑的普及及计算性能的不断提高，使得分析人员在室内进行岩土工程数值模拟成为可能。在这样的背景下，数值模拟特别是三维数值模拟技术逐渐成为当前中国岩土工程研究和设计的主流方法之一，也使得岩土工程数值模拟技术成为当今高校和科研院所岩土工程专业学生学习的一个热点。岩土工程数值模拟正是从岩、土体的受力状态出发，来分析和预测岩、土体破坏情况的一种手段。其基本原理是以典型试样的物理试验（室内试验或现场试验）获得的强度来表征整个地质体的岩、土体强度，以边界条件替代地质体周围所受的约束条件，借由本构关系表达岩、土体在外力作用下的应力-应变特性，最终了解、预测岩、土体变形破坏情况。数值模拟技术具有鲜明的时代特征，以计算机为实现平台，是信息化时代的产物。通过对数值模拟课程的学习，对数值模拟技术有了一定成度的了解，并对FLAC3D软件进行上机学习，初步掌握了一些基本知识。在对数值模拟技术和模拟软件有一定了解的基础上，并结合项目研究的需要，学生在暑假期间对砂土液化变形的数值模拟研究做了较为深入的学习，初步学习了ANASY和FLAC3D等数值模拟软件在重力坝坝基抗滑稳定性研究中的应用。针对这一阶段的学习，在总结和借鉴前人在重力坝坝基抗滑稳定性模拟研究的基础上，做了如下关于重力坝坝基抗滑稳定性非线性有限元研究的读书报告。二、读书报告重力坝坝基抗滑稳定性非线性有限元研究引言：历来人们都在不断地治理水患、建设水利工程，而重力坝l是人类最早使用的一种水坝坝型，其历史可以追溯到公元前2900年的一座高巧米顶长240米的诞生在埃及尼罗河的挡水坝，据历史记载这就是历史上第一座重力坝。从此重力坝一直使用至今，且仍是当今世界水坝中的主要组成部分。概括说来，重力坝具有以下一些优点:(l)相对安全可靠，耐久性好，抵抗渗漏、洪水漫溢、地震和战争破坏能力都比较强;(2)坝体抗冲能力好，坝体中可布置引水、泄水孔口，解决发电、泄洪和施工导流等问题;(3)断面形状简单，结构作用明确，施工方便，安全可靠;(4)对不同的地形和地质条件适应性强，任何形状河谷都能修建重力坝，对地基条件要求相对地说不太高;(5)有丰富的设计建设经验，工作可靠，使用年限长，施工放样、扩建、维护简单。设计、施工技术简单，易于机械化施工，适合于在各种气候条件下修建。综上所述，重力坝以其悠久的历史、稳重的外观、简单的结构、可靠的运行历来博得工程界人士的厚爱。爱者，必使其蓬勃也。近代以来重力坝家族家丁旺盛，在世界各地都能见其踪影，为人类造福绵绵。我国建国初期，全国兴起了修建水库、水电站的热潮5。据统计从1949一1985年间，己建成重力坝达58座，占总数的51%。近二十多年以来，国家加快了水利水电尤其大型特大型水利水电开发的速度，同时加大了开发的力度。实行改革开放以后，已具备更大规模开发和利用水资源的条件，如正在兴建的三峡水利枢纽;正在研究的南水北调工程等都是世界一流的大型水利工程。由于设计理论和施工技术的提高，水利工程建筑速度加快，规模增大，200米以上的高坝不断出现。这些水利工程的兴建为人类带来诸多裨益，全国水库在历年的防洪抗洪中都发挥了不可替代的作用。大坝成为祖国河山坚强的卫士，引导温顺祥和之水，驯服柴赘不逊之水，造福人类，功不可没。而这些工程在造福区域人民的同时也潜在地威胁着他们的生命财产安全。由于大坝是挡水建筑物2，除坝体自重外，还承受很大的水平推力和扬压力，存在向下游滑动的危险性，而重力坝的稳定性全靠坝体自重来维持，当作用在坝体上的全部荷载对坝基任意可能滑动面的滑动力(即对该滑动面的切向分量)大于其阻滑力时，坝基就要发生剪切破坏或滑动，对大坝安全造成很大威胁。在世界坝工史中，不乏由于坝基滑动而导致大坝失事的实例，据西班牙公共工程评论统计的水坝失事原因，地基破坏所占的比例高达40%，1971年英国工程师BiswaS等研究了300多座失事大坝，归纳出因基础问题失事的比率达35%5。坝基抗滑稳定性成为关系到大坝安全的关键问题之一。随着水利事业的不断发展，一方面对水的需求越来越紧迫，另一方面有利的坝址大多数已经开发蓄水。因此越来越多的大坝将修建在复杂地质构造的岩基上，而基岩的性质对大坝深层抗滑稳定性的影响巨大，这就使深层抗滑稳定性问题变得尤为重要，设计时必须对坝基深层抗滑稳定性进行验算，以确保大坝的稳定和安全。深层抗滑稳定性研究对人民生命安全及国民经济都有着极大的影响。坝基滑动破坏的形式有表面滑动和深层滑动两种基本类型2。表面滑动即沿坝基混凝土与岩土体接触面发生的滑动破坏，如格莱诺坝的破坏。深层滑动有两种情况:一种是发生在均质坝基的剪切破坏，其滑动面轨迹近弧形;另一种是沿岩土体中的缓倾角夹层或软弱结构面的滑动，如奥斯汀坝就是由于坝基沿建基面以下60cm的粘土页岩夹层发生滑动和破坏。应用D一P4屈服准则、以塑性屈服区贯通为安全判据、强度折减过程中适当调整弹性模量和泊松比，并考虑不同范围降强、以本文改进的扬压力施加方式考虑水荷载。在此基础上，本章将结合ANsys和ABAQus有限元程序的强大分析功能，对其屈服模式进行扩充，并加以一定的修改，用参数化语言编写适用于重力坝抗滑稳定的强度折减有限元方法计算分析程序，并结合具体工程实例，即金沙江向家坝的两个典型坝段进行非线性有限元针对重力坝深层抗滑稳定性问题的计算分析，对结果进行多方面综合分析，总结规律，提出需要改进之处。1基本原理1.1 广义等 K 法计算原理 Sarma法是岩质边坡稳定分析中常用的一种方法，该方法对土条进行斜分条，并假定在底滑面和倾斜界面上均达到极限平衡，据此推导出计算安全系数的公式。将该方法推广到重力坝深层抗滑稳定分析时，一般假设条块界面为垂直面且对穿越混凝土坝体部分的条块界面的黏聚力c和摩擦系数f值取为0。针对重力坝垂直条分的特殊情况，广义等K法的计算公式为：式中系数按下式计算：1.2 分项系数极限状态设计方法 电力行业混凝土重力坝设计规范规定使用极限状态表达式来进行结构和地基的稳定性验算，要求将建筑物抗滑稳定分析从传统的安全系数表达式改为极限状态表达式，对承载能力极限状态验算的表达式为：在进行重力坝抗滑稳定分析时，可以将分项系数极限状态方法的式改造为使用安全系数K的表达形式：对式中的荷载和抗力各项乘上相应的分项系数后，仍可按传统的极限平衡方法和强度折减意义求解安全系数K，即求解满足K值2 ANSYS 工程建模由于本计算中涉及的软弱结构面较多，表现出很强的非线性特点，因此本文结合ANSYS和ABAQUS各自的优势进行计算分析，即用ANSYS建模，形成几何模型信息、材料信息以及约束和荷载信息，将这些信息转化为ABAQUS命令流文件，进行ABAQUS非线性有限元计算，并进行后处理。计算中要建立的数学物理模型应能模拟大坝结构特点、图中所示坝基的地质构造特征和坝基的主要地质结构面，以满足计算成果精度的需要。在建模时应充分反映结构的几何、材料分区和受力等方面的特征。有限元模拟的第一步就是离散化，这部分以ANSYS为平台来完成。划分网格的疏密对结果有一定的影响，网格划分太疏则造成计算精度很低，如果网格剖分的太密，则会造成计算时间成倍增长，但是对计算结果的精度只有很小的影响，因此，需要对有限单元模型进行收敛性研究以确定合适的网格密度。其中对坝体和坝基进行材料分区，坝基面附近网格较密，以反映材料的突变及大坝在建基面附近的应力集中特点;对于基岩，模拟了不同岩性的岩体，考虑了主要软弱夹层。在本构仿真方面，由于研究重点在于坝基，因此对坝体混凝土采用了弹性模型，对坝基岩体采用弹塑性一断裂模型。在受力仿真方面，按加载路径分别计入初始地应力场、坝体自重、水压力、淤沙压力等作用，以准确反映工程结构受载的物理力学环境和各类荷载的作用效应，初始地应力场仅考虑重力场。图2-1 重力坝网格图图2-2 泄洪道坝段网格示意图3 ABAQUS非线性有限元计算ABAQUS对初始应力场的模拟一般都是采用生死单元法，即在整体建模之后将坝体部分“杀死”，对坝基部分施加自重作用，以模拟初始自重应力场，计算之后，得到初始应力资料，生成输入文件，将初始应力指定为一科，荷载，然后再“激活”坝体部分，施加所有荷载，继续下面的计算。这样坝基部分存在自重引起的初始应力，而位移为零。其中“杀死”单元并不是删除，而是将这部分单元的劲度矩阵乘以一个很小的系数网，并从总质量矩阵消去单元的质量来体现单元的“死”，同时无活性（压力、热应变等）荷载被设置为零。图3-1 坝基未折减时竖向位移等值线图图3-2 坝基未折减时最大主应力等值线图图3-3 坝基未折减时最小主应力等值线图图3-4 天然条件下坝基塑性区图通过分析计算天然状况坝体及坝基的受力状况和塑性区开展情况，发现天然状况下，坝体和坝基变形不大，应力也不大，只在小范围出现拉应力，且数值不大。塑性开展范围不大，大部分位于软弱夹层和坝体与坝基交接面处。这说明在天然状况下，坝基比较安全，满足设计要求。如果按照计算过程中应力和应变不收敛来看，计算出来的安全系数是2.50；如果按照突变的角度来考虑，计算出来的安全系数是2.45。综上所述，以位移和能量突变以及塑性区贯通为失稳判据得到的安全系数较为接近，以安全起见，可认为在工况I情况下坝体及基岩系统的安全系数为2.45。且整个体系破坏主要为压剪破坏，表现为水坝整体与坝基交界处一直到坝址附近质量较差的岩体，最终从下游处滑出，同时，坝址附近的质量较差的岩体也是危险部位。基岩内部从中间向两边有分离断开的趋势。4 结论本文首先详细讨论了目前稳定性计算的方法，指出各种方法的优缺点，然后总结有限元基本理论与分析步骤，分析材料非线性有限元的特点及处理方法，讨论大型商业软件ANSYS和ABAQUS的应用，通过典型算例研究了其中一些关键问题，包括屈服准则类型、在各种程序中的实现、对结果的影响以及选用依据；安全系数的定义、安全判据的形式和对结果的影响分析；单元类型及网格划分等对有限元分析的影响。最后借助大型工程实例金沙江向家坝，利用以上讨论的成果进行较为系统的深层抗滑稳定性计算分析，考察工程的稳定性、失稳破坏的发展过程及轨迹。参考文献：l潘家铮.重力坝设计M，北京:水利电力出版社，1987.2汝乃华.重力坝M，北京:水利电力出版社，1983.3混凝土重力坝设计规范,DL5108一1999中国电力出版社,2000.4陈宪宏.关于岩基混凝土重力坝抗滑稳定的一些问题J.水力发电学报，Vol.24No.1，Feb.，2005.5陈国立.重力坝深层抗滑稳定基本特点的初步分析J.6博弈创作室.ANSYS7.0基础教程与实例详解M,北京:中国水利水电出版社，2003.7Zienkiew O C，HumPheson C and Lewis R W. Associated and non-Assoeiated Viseo-Plastieity and Plasticity in Soil Machanies J. Geoteehnique，1975，25(4):671689.8Sarma，S.K. Stability analysis of embanklllents and slopes J. ASCE Journal of the Geoteehnieal Engineering Division. 1979. 105 (12) : 15111524.9Junuthula Narasimha Reddy. An Introduetion to Nonlinear Finite Element AnalysisM. Oxford University Press，2004.10张弘毅，迟道才.重力坝深层抗滑稳定的可靠性分析J.沈阳农业大学学报，A09期: 5156.第二部分：上机实习报告一、上机实习题目一计算模型分别如图1、2、3所示，边坡倾角分别为30、45、60，岩土体参数为：密度2500 kg/m3, 弹性模量E1108 Pa，泊松比0.3，抗拉强度t0.8106 Pa，内聚力C4.2104 Pa，摩擦角17试用FLAC3D软件建立单位厚度的计算模型，并进行网格剖分，参数赋值，设定合理的边界条件，利用FLAC3D软件分别计算不同坡角情况下边坡的稳定性，并进行结果分析。附换算公式：1 kN/m3= 100 kg/m3剪切弹性模量：=3.85107 pa体积弹性模量：=8.33107 pa图1 倾角为30的边坡(单位：m)图2 倾角为45的边坡(单位：m)图3 倾角为60的边坡(单位：m)1、 代码工况一：角度为30的边坡newgen zone brick size 50 1 10 p0 (0,0,0) p1 (100,0,0) p2 (0,1,0) p3 (0,0,40) p4 (100,1,0) p5 (0,1,40) &p6 (100,0,40) p7 (100,1,40)gen zone brick size 30 1 5 p0 (40,0,40) p1 (100,0,40) p2 (40,1,40) p3 (74.64,0,60) p4 (100,1,40) p5 (74.64,1,60) p6 (100,0,60) p7 (100,1,60)model elasprop bulk 3e9 shear 1e9 density 2500set grav 0, 0, -9.81fix x range x -0.1 0.1fix x range x 99.9 100.1fix yfix z range z -0.1 0.1solveini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0model mohrprop bulk 8.3e7 shear 3.8e7 fric 17prop coh 4.2e4 tens 8e5solve fos associatedsave slope30.sav工况二：角度为45的边坡newgen zone brick size 50 1 10 p0 (0,0,0) p1 (100,0,0) p2 (0,1,0) p3 (0,0,40)gen zone brick size 30 1 5 p0 (40,0,40) p1 (100,0,40) p2 (40,1,40) p3 (60,0,60) p4 (100,1,40) p5 (60,1,60) p6 (100,0,60) p7 (100,1,60)model elasprop bulk 3e9 shear 1e9 density 2500set grav 0, 0, -9.81fix x range x -0.1 0.1fix x range x 99.9 100.1fix yfix z range z -0.1 0.1solveini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0model mohrprop bulk 8.3e7 shear 3.8e7 fric 17prop coh 4.2e4 tens 8e5solve fos associatedsave slope45.sav工况三：角度为60的边坡newgen zone brick size 50 1 10 p0 (0,0,0) p1 (100,0,0) p2 (0,1,0) p3 (0,0,40)gen zone brick size 30 1 5 p0 (40,0,40) p1 (100,0,40) p2 (40,1,40) p3 (51.55,0,60) p4 (100,1,40) p5 (51.55,1,60) p6 (100,0,60) p7 (100,1,60)model elasprop bulk 3e9 shear 1e9 density 2500set grav 0, 0, -9.81fix x range x -0.1 0.1fix x range x 99.9 100.1fix yfix z range z -0.1 0.1solveini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0model mohrprop bulk 8.3e7 shear 3.8e7 fric 17prop coh 4.2e4 tens 8e5ini dens 2500solve fos file slope60.sav associated2、计算结果图2.1建模模型图4 30的斜坡计算模型图5 45的斜坡网格剖分图6 60的斜坡网格剖分2.2各个角度斜坡主应力云图图7 30斜坡最大主应力图8 30斜坡最小主应力图9 45斜坡最大主应力图10 45斜坡最小主应力图11 60斜坡最大主应力图12 60斜坡最小主应力对于FLAC，主应力中负值表示压应力，正值表示拉应力。从图7-图12可以看出，从边坡坡面往坡体内沿深度方向主应力量值逐渐增大，最大及最小主应力量等值线总体平行，在模型底部以高压应力为主，模型上部以低压应力为主，坡面后部上局部地段存在拉应力作用，且随着坡角的增大拉应力的作用范围越大，越靠近坡面，值越大，说明坡角越大边坡越易产生由拉力引起的滑动。2.3各个角度斜坡剪应变增量云图图13 30的斜坡剪应变增量云图图14 45的斜坡剪应变增量云图图15 60的斜坡剪应变增量云图从图13-图15可以看出，剪应变增量在边坡内部较大，呈圆弧状，且剪应变较大的区域呈贯通状态，表明边坡具有滑动的可能性，且可能在坡脚处剪出。随着坡角的增大，剪应变越大，表明坡角越大边坡越不稳定。2.4各个角度斜坡总位移云图图16 30的斜坡总位移云图图17 45的斜坡总位移云图图18 60的斜坡总位移云图由图16-图18可知，滑坡的位移主要发生在坡面并延伸到坡里呈圆弧状，在坡脚处位移量要大于其他位置，表明边坡有滑动的可能性且坡脚处最先滑动。随着坡角增大，最大位移量增大，表明坡角越大边坡越易滑动。FLAC3D采用强度折减法求解边坡的安全系数，即定义安全系数为岩土体的实际抗剪强度与临界破坏时的折减后剪切强度的比值。（FLAC3D中所述安全系数实际为稳定性系数。）基于FLAC3D强度折减法求解的三个不同坡度斜坡的稳定性系数见下表：表1 不同角度斜坡稳定性系数坡角304560稳定性系数1.461.130.97由表1当坡角为30时，斜坡是稳定的，变形主要发生在斜坡内部，临空处未发生大的变形。随着坡角的增大，斜坡坡面处变形增大，并向坡顶扩展，而变形区域却向临空面缩小。由计算出的稳定性系数可以看出，坡角越大斜坡越不稳定。综上所述，采用FLAC3D进行边坡稳定性分析可以取得较好的效果。通过应变云图、位移云图等可以有效地判断边坡潜在滑动面的位置，并能直观地得到边坡中不同位置的变形值，这有助于我们进一步深入分析边坡的稳定问题。且用FLAC3D分析得到的边坡安全系数与用常规的极限平衡法得到的安全系数是较为吻合的，这说明FLAC3D确实是一种分析边坡稳定性的行之有效的手段。二、上机实习题目二边坡开挖，参数与习题一相同图19 边坡开挖示图分六块区域建立模型，如图20所示：图20 建模示意图程序如下：newgen zone brick size 100 1 20 p0 (0,0,0) p1 (100,0,0) p2 (0,1,0) p3 (0,0,40) group 基底gen zone brick size 57 1 10 p0 (43,0,40) p1 (100,0,40) p2 (43,1,40) p3 (63,0,60) p4 (100,1,40) p5 (63,1,60) p6 (100,0,60) p7 (100,1,60) group 基底gen zone brick size 3 1 5 p0 (40,0,40) p1 (43,0,40) p2 (40,1,40) p3 (50,0,50) p4 (43,1,40) p5 (50,1,50) p6 (53,0,50) p7 (53,1,50) group 基底gen zone wedge size 4 1 5 p0 (40,0,40) p1 (50,0,50) p2 (40,1,40) p3 (45.77,0,50) p4 (50,1,50) p5 (45.77,1,50) group 开挖gen zone brick size 3 1 5 p0 (50,0,50) p1 (53,0,50) p2 (50,1,50) p3 (60,0,60) p4 (53,1,50) p5 (60,1,60) p6 (63,0,60) p7 (63,1,60) group 开挖gen zone brick size 5 1 5 p0 (45.77,0,50) p1 (50,0,50) p2 (45.77,1,50) p3 (51.55,0,60) p4 (50,1,50) p5 (51.55,1,60) p6 (60,0,60) p7 (60,1,60) group 开挖fix x y z range z -0.1 0.1fix x range x -0.1 0.1 fix x range x 99.9 100.1fix y model elasticprop dens 2500 bulk 1e10 shear 1e10set grav 0 0 -10solvedelete range group 开挖ini xvel 0 yvel 0 zvel 0 ini xdis 0 ydis 0 zdis 0model mohr prop den 2500 bu shear c 42000 ten 0.8e6 fri 17 dila 20 ran group 基底solve fos file slope fos kaiwa. sav1、构建模型：图21 网格剖分图图22 开挖前的位移云图图23 开挖后的位移云图从位移等值线图中可以很明显地看出，60的边坡经过开挖以后，发生了严重的变形，坡脚部分土体已经剪出，其位移量甚至达到了十几厘米，坡顶土体也产生了较明显的向下位移，整个坡体相对于原来的60边坡来说，由于开挖引起的变形量更大，而且从图中也可以很明显地看到一条贯通的剪切滑动面，如在开挖的过程中不采取任何支护措施，边坡将可能沿此滑动面产生破坏。2、开挖前与开挖后剪应变增量云图和安全系数图图24 开挖前剪应变增量和安全系数图图25 开挖后剪应变增量和安全系数图该边坡开挖计算结果分析见表2：表2 开挖前后稳定性系数状态开挖前开挖后稳定性系数0.991.25从上表可以看出，60的斜坡在开挖后稳定性系数提高。 FLAC3D对于边坡开挖工程的施工具有十分重要的指导意义。从这个开挖实例中可以看出，FLAC3D在边坡工程中具有很好的实用性，它可以分阶段模拟边坡开挖的过程，动态监测边坡开挖过程中坡体产生的变形，从而为我们确定开挖量、开挖角度等及是否应该采取支护措施以及怎样支护提供依据。参考文献1 丁秀丽,付敬,张奇华.三峡水库水位涨落条件下奉节南桥头滑坡稳定性分析J.岩石力学与工程学报，2004，23(17):2913-2919.2 胡新丽.三峡水库水位波动条件下滑坡抗滑工程效果的数值研究J.岩土力学，2006,27(12)3朱大鹏. 三峡库区典型堆积层滑坡复活机理及变形预测研究博士论文 D.武汉：中国地质大学（武汉）.2010.4王洪兵，姚磊华，文海. 基于进行堆载作用下边坡稳定性分析J.勘察科学技术，2012.5-10.5 李晓，张年学，廖秋林，吴剑波. 库水位涨落与降雨联合作用下滑坡地下水动力场分析J.岩石力学与工程学报，2004,23（21）：3714-3720.6刘洋，李世海，刘晓宇. 基于连续介质离散元的双重介质渗流应力耦合模型J.岩石力学与工程学报，2011,30（5）第三部分：数值模拟技术研究应用实例分析：坝基抗滑稳定性的研究一、工程概况小南海水电站中坝址位于重庆市江津区珞璜镇下游约1.4km（图1-1），上距江津市38.5km，下距重庆市40.0km。长江呈NE向流经坝址区，河道长约6-7km，河谷开阔，河床枯水位高程一般为174.5m。至大中坝江心洲长江分流为左、右两支，左支为主河道，河流流向NE26，枯水期河床宽约600m，窄处仅380m；右支为汊河，流向由NE77渐变至NE9与主河道汇合，枯水期河床宽约200m。大中坝江心洲呈梭形，洲顶上游高程较低，多在196m以下；中下游高程在196m以上，最高为203.2m，坝轴线处最高为192.1m。顺水流向长约4.1km，顺坝轴线处宽约1150m。坝址河段属斜交横向谷。图1-1 小南海水电站中坝址地理位置图重庆长江小南海水电站是以发电和航运为主，兼有拦沙减淤、滞洪、灌溉供水等综合利用效益的大型水电水利工程，其主要建筑物包括河床式电站厂房、通航、泄洪、挡水及过鱼设施。主要工程特性指标见表1-1。坝址控制流域面积70.5万km2，多年平均流量8650m3/s，相应多年平均径流量2685亿m3。小南海水电站大中坝坝址（中坝址）水库正常蓄水位为196m，死水位194m，校核水位以下总库容为13.0亿m3；中坝址枢纽布置代表方案的最大坝高为67m，装机容量1680MW。按照防洪标准和水电枢纽工程等级划分及设计安全标准（DL5108-2003）的规定，小南海水电站属1等大（1）型工程，主要建筑物溢流坝、电站厂房、挡水坝、船闸、鱼道为1级建筑物，次要建筑物上下游边坡及导墙等为3级建筑物。二、坝基稳定性问题综述2.1 问题的提出按照岩体工程地质力学的观点，工程岩体稳定问题是重要和关键的工程地质问题之一，而其中最突出的乃是含软弱夹层坝基岩体的抗滑稳定性问题。由于坝基软弱夹层的影响造成工程事故的，国内外不乏其例。因此含软弱夹层坝基岩体的抗滑稳定性分析历来为大家所高度重视。坝基岩体抗滑稳定性分析的力学方法主要有极限平衡法，有限元法和结构地质力学模型法。后两种方法虽具极限平衡法所不及的优点（如能了解全面应力和位移状态，了解失稳的过程和形态等），但他们的精度除受地质工作和力学参数的影响外，结构地质力学模型还受模型材料的模拟正确性和外力作用的相似性等因素影响，所以目前在国内的实际工程中都只作校核用。由于极限平衡法计算简便，成果一般偏于安全，便于设计人员制定设计方案和控制工程量，因此仍得到广发的应用；目前部分学者已应用塑性理论进行极限平衡法计算，国外还将塑性理论推广应用到三维问题。我国、美国、日本等国的设计规范均建议或规定采用极限平衡法。在国内进行重力坝稳定分析时，按极限平衡理论推导公式，不考虑滑动面凝聚力，其计算方法有：被动抗力法、等稳定法和剩余推力法。这些计算方法都是按平面课题考虑的，也就是沿水流方向切取单宽断面或选取一个坝段作为分析对象，不考虑侧向切割面的影响；但是正如潘家铮教授指出的：实际工程的边界条件是很复杂的，深层滑动问题往往是一个空间课题。坝基岩体处于临界失稳状态时，不仅软弱夹层和抗力体将发挥极限抗剪作用，而且侧向切割面也会提供一定的阻滑力。因此，传统的极限平衡理论在重力坝稳定性分析中具有不可忽视的缺陷。2.2 坝基失稳模式研究2.2.1重力坝失稳模式根据过去坝工失事的经验以及室内模型试验的情况来看，大坝失稳形式主要有两种情况：第一种情况是岩基岩体强度远大于坝体混凝土强度，同时岩体坚固完整无显著的软弱结构面，这时大坝的失稳多半是沿坝体与岩基接触处产生，这种破坏形式成为表层滑动破坏（图2-1）。图2-1重力坝表层滑动模式示意图第二种情况是在岩基内部存在着节理、裂隙和软弱夹层，或者存在着其他不利于稳定的结构面，在此情况下岩基容易产生深层滑动破坏（图2-2）。除了上述两种破坏形式之外，有时还会产生所谓混合滑动破坏形式，即大坝失稳时一部分沿着混凝土与岩基接触面滑动，另一部分则沿岩体中某一滑动面产生滑动，因此，混合滑动的破坏形式实际上是介于上述两种破坏形式之间的情况。（a）单斜滑动 （b）双斜滑动 （c）尾岩抗力体挤压破坏 （d）尾岩抗力体隆起破坏图2-2重力坝深层滑动模式示意图根据软弱结构面空间展布性状不同，深层滑动一般可分为4种类型：单斜剪切滑动破坏、双斜剪切滑动破坏、尾岩抗力体挤压破坏和尾岩抗力体隆起破坏。在这4种类型中，当坝趾下游没有临空面和倾向上游的结构面时，坝体的失稳可以归结为后两类，尾岩抗力体挤压破坏和隆起破坏均是因为尾岩的破坏而导致的坝体失稳。当坝基只有倾向下游的结构面时，依据结构面的产状可以将滑移模式分为两种型式。1、缓倾角结构面产状平缓的滑移模式-水平滑移模式。这种滑移模式条件下，结构面的倾角较小，接近水平，此时由于坝体自重产生的下滑力较小，坝趾处尾岩抗力体依据结构面的埋深不同而不同。坝体带动结构面以上岩体整体向下游滑动，沿软弱夹层产生水平滑移，并使尾岩抗力体产生局部压碎和剪断破裂，最终导致坝体失稳。2、结构面倾角较大的滑移模式-倾斜滑移模式。这种滑移模式条件下，由于坝体自重产生的下滑力随结构面倾角的增大逐渐增大，抗力体也逐渐增大，坝体在推力作用下，使下游软弱夹层发生较大的压缩变形，导致坝基岩体沿软弱结构面向下游移动，最终尾岩失效而使大坝失稳破坏。2.2.2 各坝段失稳模式探讨由于大坝影响范围内，地层条件复杂，不同坝段坝基内的剪切带分布状态差异较大，且不同坝段的荷载分布情况也差异较大，因此在进行三维数值计算时，选择分坝段建立模型。根据小南海水电站各坝段坝基下部主要剪切带分布图（图7-47-7，只留下了III类剪切带和部分分布连续的II类剪切带），我们根据剪切带的类型和分布规律，可对各坝段可能的破坏模式进行初步判断，然后用理论计算的方法进行稳定性分析。由图7-4、7-5、7-6、7-7可以看出，坝基剪切带均为缓倾角或近水平，且倾向下游。右岸溢流坝下部的剪切带分布类型为缓倾角、倾向下游剪切带，且剪切带出口在大坝设计的消力池部位，此类剪切带可能发生的破坏类型为单斜滑动；而左岸溢流坝下部的剪切带则是倾向下游的，且离大坝建基面很近，若坝前有反倾向节理裂隙配合，则此类剪切带可能发生的破坏类型为双斜滑动；左岸厂房、右岸厂房坝基岩体内剪切带分布不密集，且倾角较小，但其建基面以上大坝前后剪切带分布较为密集，因此认定该坝段整体稳定性好是合理的，但需要注意在开挖过程中基坑侧壁的稳定性问题，其可能的破坏类型为尾岩抗力体挤压破坏。因此，可以将各坝段可能的破坏型式分为三类：一类为单斜滑动（右岸溢流坝）；另一类为双斜滑动（左岸溢流坝）；第三类为近水平滑移引起的尾岩抗力体挤压破坏（左、右岸电站2.3 坝基稳定性理论计算根据SL319-2005混凝土重力坝设计规范，对于溢流坝坝基深层抗滑稳定计算作用在该坝段上的荷载有混凝土坝体的自重、破坏面上的基岩重、水压力、扬压力，不考虑淤沙压力、库水的动水压力和地震力的作用。2.4 坝基稳定性数值计算2.4.1 数值计算综述传统的坝基稳定性解析法是基于大量的假设和推测的前提下进行的。在取计算参数时，需要将地层参数做加权平均；同时很难模拟地面的起伏和上部荷载的分布情况等；并且由于常规的稳定性解析法不考虑侧向阻滑力。因此，常规的解析法计算结果偏安全，且不能直观的了解坝基岩体内的应力分布状态，需要加以改进。数值计算方法，如有限元法、边界元法、有限差分法等方法在出现之初由于自身算法和计算机水平的限制，只是作为坝基稳定性计算的辅助手段，主要起验证作用。而随着各种数值方法在理论上的日趋成熟，和计算机技术的高速发展，数值计算方法受到人们越来越多的重视。相比传统的解析方法，数值计算方法具有以下的优点：1、能进行复杂的、大量的非线性运算，这一点是解析法望尘莫及的；2、能模拟较为复杂的地层岩性条件，并分别赋予计算参数；3、较为合理的模拟所选取研究区域的边界条件，比如在坝基稳定性计算中就可以将侧向阻滑力作为边界条件施加于所建立的模型上；4、能直观的显示模型内部各单元和节点的应力、应变状态，比如在坝基稳定性计算中能很直观的显示基岩内部应力场在大坝开挖、浇筑直至蓄水整个过程中的变化趋势和大小；5、能考虑坝基内的渗流应力的影响等等。本节将通过通用有限差分软件FLAC3D建立坝基稳定性计算模型。FLAC3D采用“混合离散法”来模拟材料的塑性破坏和塑性流动；采用显式差分法求解微分方程，可以很方便的求得应力增量、不平衡力并跟踪系统的演化过程。值得一提的是，FLAC3D中有一个专门的“接触面单元”来模拟岩体中的节理、裂隙、断层、剪切带等软弱层面。接触面单元可以很好的模拟红层地区大坝稳定性的主要控制因素“层间剪切带”。2.4.2 FLAC3D接触面单元理论简介FLAC3D中接触面单元由一系列三节点的三角形单元构成，接触面单元将三角形面积分配到各个节点中，每个接触面节点都有一个相关的表示面积。每个四边形区域面用两个三角形接触面单元来定义，然后在每个接触面单元顶点上自动生成节点，当另外一个网格面与接触面单元相连时，接触面节点就会产生。图2-3 接触面节点相关面积的分布FLAC3D中接触面是单面的，认识这一点很重要，这点与二维FLAC中所定义的双面接触面不同。可以把接触面看做“收缩带”，可以在指定面上拉伸，从而导致接触面和与之可能相连的其他任何面的相互刺入变得敏感。接触面单元可以通过接触面节点和实体单元表面（称之为目标面）之间建立联系。接触面法向方向所受到的力由目标面方位决定。在每个时间步计算中，首先得到接触面节点和目标面之间的绝对法向刺入量和相对剪切速度，再利用接触面本构模型来计算法向力和切向力的大小。在接触面处于弹性阶段，时刻接触面的法向力和切向力通过下式得到。式中，为时刻的法向力矢量；为时刻的切向力矢量；为接触面节点贯入到目标面的绝对位移；为相对剪切位移增量矢量；为接触面应力初始化造成的附加法向应力；为接触面应力初始化造成的附加切向应力；为接触面单元的法向刚度；为接触面单元的切向刚度；为接触面节点代表面积。图7-9为接触面的本构模型示意图。对于Coulomb互动的接触面单元，存在两种状态：相互接触（intact）和相对滑动（broken）。根据Coulomb抗剪强度准则可以得到接触面发生相对滑动所需要的切向力为：式中：为接触面的凝聚力，为接触面的摩擦角，u为孔压。当接触面上的切向力小于最大切向力（）,接触面处于弹性阶段；当接触面上的切向力等于最大切向力（=）,接触面进入塑性阶段。在滑动过程中，剪切力保持不变=，但剪切位移会导致有效法向应力的增加。式中：为接触面的膨胀角，为修正前的剪力大小。图2-4 接触面单元原理示意图如果接触面上存在拉应力并且超过了接触面上的抗拉强度，接触面就会破坏，切向力和法向力就会为零。默认情况下，抗拉强度为零。节点上计算出的法向力分布在目标面上，节点上的剪切力分布在与节点相连的反方向的面上。把这些力加权平均分布到每个面的节点上。接触面刚度加到了接触面两边节点的计算刚度上，这是为了保持数值计算的稳定性。接触面接触性体现在接触面节点上，并且接触力仅仅在节点上传递。节点应力假定在节点的代表区域上统一分布。接触面参数与每个节点都有联系，并且每个节点也可以有自己不同的参数。接触点上唯一增量的矢量方向可分解为法向和剪切方向，法向和剪切方向的总荷载由下面两式决定：接触面有3中选项可用来模拟实际工程中不同的情况，如下：1、胶合模型。如果接触面是胶合的，不允许张开和滑移，但弹性位移根据给定的刚度仍然会发生。2、库伦剪切模型。库伦剪切强度准则如下式所示：3、抗拉粘结模型。接触面用property关键字设定界面的参数，在所有情况下，必须设定法向刚度kn和剪切刚度ks，对这些参数必须使用一直的单位设置。其他的参数还包括：cohesion：粘结力dialation：膨胀角friction：摩擦角tension：抗拉强度2.4.3 坝基分区数值计算参数数值计算用到的地层、剪切带及结构面的参数均来自重庆长江小南海水电站预可行性研究工程地质勘察报告（2011-03），分别见表2-1，结构面的参数见表2-2。表2-1 数值计算地层参数表岩石名称天然重度（）抗拉强度变模（E0）弹模（Ee）岩体本身抗剪断强度砼/岩体抗剪强度微、新岩体弱风化强风化KN/m3MPa（GPa）（Gpa）FC（MPa）FC（MPa）FC（MPa）FC（MPa）粉砂质粘土岩、粘土岩24.5250.60.751.5222.50.650.30.450.20.300.050.40粘土质粉砂岩、粉砂岩25.52611.22.53.03.03.50.750.50.550.300.350.10.450长石岩屑砂岩25.5456810121.01.00.700.800.450.40.60表2-2 数值计算剪切带及结构面参数表剪切带及结构面类型抗剪断强度抗剪强度特征描述fc（MPa）fc（MPa）类岩块岩屑型0.350.40.10.120.320破碎岩块、岩屑类岩屑夹泥型0.30.350.050.080.280岩屑、部分泥化类1泥夹岩屑0.250.280.030.050.220具泥化特征2泥型0.20.220.010.20泥化层理面0.450.550.150.20.350.40软岩层面取低值缓倾角结构面0.200.250.050.200粘土岩中具泥化0.300.320.080.10.280.300砂岩中破碎带2.4.4 坝基分区稳定性计算二维数值计算2.4.4.1 二维数值模型建立在建模过程中，参考左岸溢流坝横剖面图，按不同岩性分层建立地层单元，在计算时分别赋予不同的计算参数。由于左岸溢流坝建基面下部0-5m的范围内正好有一条剪切带穿过，该剪切带在模型中用接触面单元来模拟。该剪切带倾向大坝下游，在该剖面上的视倾角约为4，属于缓倾角裂隙。在固定这条剪切带的前提下，大坝最可能的破坏滑移方式是大坝下游有缓倾角裂隙与该剪切带配合组成滑移面。以该模式为基础，建立了示意性的模型，通过改变下游缓倾角裂隙的倾角来研究其倾角对坝基稳定性的影响。大坝不同阶段的模型见图2-5、2-6、2-7，以下游反向缓倾角裂隙的倾角为35为例。图2-5 初始地层剖面模型图2-6 建基面开挖后模型图2-7 大坝建立后模型此次数值计算的过程分为4个阶段：1、建立初始地层模型。计算初始应力场和初始位移场；在计算前需先定义模型边界条件和分层赋予模型及剪切带计算参数。2、模拟大坝基坑开挖。在初始应力场的基础上计算开挖后应力场，同时计算由于基坑开挖引起的地层的位移场，在计算此阶段位移场之前需先将初始位移场清零；此阶段的计算结果可反映在剪切带存在的情况下大坝基坑的稳定性。3、模拟大坝建立。在上一阶段应力场和位移场的基础上计算大坝建成后的应力场和位移场；4、模拟水库蓄水。在上一阶段应力场和位移场的基础上，对模型施加静水压力，计算水库蓄水后由于上下游水位差引起的应力场和位移场的变化。2.4.4.2 二维数值计算参数的选取及边界条件的定义为了充分考虑大坝基础开挖及大