基于有限元法的观音岩水电站重力坝12#坝段静力抗滑稳定分析

基于有限元法的观音岩水电站重力坝12#坝段静力抗滑稳定分析摘要重力现作为一种古老的顼型,因其结构简单,设计、施工方便,适宜性好而大量应用在国内外大型水利水电工程中。随着科技的发展,重力坝的设计理论与施工方法得到了快速的发展。尽管如此,重力坝失事时有发生,据统计,全世界失事的大现总数多达15万座,在众多的原因中,坝基失稳是最主要的原因之一。因此,坝基稳定性问题成为重力坝设计与施工关注的重点与难点,其中又以深层抗滑稳定性问题最为突出。重力坝的应力分析目的是为了判断坝在施工期和运用期是否满足强度和稳定方面的要求，同时也为研究与设计、施工有关的其他问题（如确定坝体混凝土的标号分区以及在某些部位配置钢筋等）提供依据。重力坝的应力分析和抗滑稳定分析一直都是判断大坝是否安全的重要指标。目前阶段设计中使用最广泛的是材料力学法和刚体极限法，其计算简单，有完整的规范与之对应，但由于该理论对实际问题做了一些假设，计算结果比较粗略。有限元法可以一定程度上弥补传统理论的不足，该理论可以更全面地反映结构内部各点应力、位移分布规律，找出应力集中部位和抗滑系数较小的危险滑动面。然而该理论往往在应力集中处出现的结果失真，且尚没有相应的规范对其明确的要求。 本文论述了材料力学法、刚体极限平衡法和有限元法基本理论，并分别采用该理论对观音岩水电站重力坝应力、抗滑稳定进行分析。通过对比发现，不同计算理论其计算结果有一定差异性，但其结果均符合规范要求。通过对比分析得出结论应力大小和分布规律与地基刚度有很大关系，距坝基越远，应力与分布规律越相似。通过两种理论综合分析，可以判定该重力坝处于稳定状态。因此，在一些重要工程中，可以分别采用不同计算方法进行综合、全面地分析，以了解该工程真实运行情况。 关键词：浆砌石重力坝；材料力学法；有限元法；刚体极限平衡法；FLAC3D法第1章 绪论1.1工程概况观音岩水电站位于云南省丽江地区华坪县（左岸）与四川省攀枝花市（右岸）的界河塘坝河口附近，为金沙江中游河段规划的八个梯级电站的最末一个梯级，上游与鲁地拉水电站相衔接。该工程以发电为主，兼顾防洪、灌溉、旅游等综合利用的水利水电枢纽工程。电站地理位置适中，对外交通十分便利。电站距攀枝花市直线距离约27km，距华坪县城公路里程约42km，距昆明市公路里程约350km。坝址左岸有简易公路与华坪至攀枝花市的公路相接。成昆铁路支线格里坪站距坝址10km余，小型机动船可从格里坪通至坝址。塘坝河口以上约4km干流河段为云南与四川界河，再往上则全属云南省；塘坝河口以下直至雅砻江口段则全处于攀枝花市境内。观音岩水电站坝址以上金沙江河道长2703km，落差4385m，河道平均坡降1.62，控制流域面积约26万km2，多年平均流量1880m3/s。金沙江流域属典型的季风气候区，每年510月，主要受西南季风控制，同时也受东南季风影响，水汽充沛、降雨集中，为本流域雨季或汛期。年平均降水量呈由南向北递减的总趋势。流域内暴雨较小，其分布中、下游较上游大。混合坝方案坝轴线布置在、勘探线之间，坝顶高程1139m，最大坝高159m。电站正常蓄水位1134m，死水位为1126m，相应死库容为16.89108m3，调节库容为3.83108m3，总库容约20.72亿m3。水库回水至鲁地拉水电站坝址，库长约100km，属峡谷型水库。装机规模为5600MW=3000MW。1.2 重力坝应力分析发展状况重力坝应力分析的主要目的是了解坝体应力分布规律及最大应力，以确定坝体安全系数较低的部位。对结构应力分析一般可通过模型试验、理论计算等方法实现。经过100多年的不断探索、研究，重力坝设计理论日趋完善。随着科学技术提高，国民经济发展，重力坝体型越来越复杂，高度越来越高，需处理问题越来越复杂，因此需要与之配套的设计理论与相关规范。目前常用的理论计算方法包括有材料力学法、弹性理论法及有限元法等。1.2.1 材料力学法材料力学法是设计中一种常用的使用时间最长的计算方法，现行规范均以材料力学法为依据。该理论计算原理简单，可以用于各种荷载和坝体外形，有较长的使用经验和与之配套的应力控制标准，该理论得到了许多实际工程和模型试验的论证，有很强的实用性。材料力学法在对重力坝进行应力分析时，需要做一些基本假设。坝体材料为均质、各向同性的弹性体；视坝体完全固结于坝基上，以悬臂梁形式考虑；材料力学法假定铅直正应力按直线分布，并利用微分体导出其他应力。由于混凝土坝或浆砌石坝体并非是均质、连续的弹性体，在荷载的长期作用下，坝体会发生一定的塑性变形和徐变，使得应力集中部位发生应力重分布，应力近于直线分布。从极限设计观点来看，采用此假设计算出来的结果可以满足设计要求。同时材料力学法计算结果相对粗糙，具有一定局限性，如材料力学法无法用来求解坝址、坝踵、廊道等坝体应力集中问题，不能分析地基应力、变形分布及其弹性模量对坝体应力分布的影响，不能了解坝基各点应力分布等。 材料力学法原理简单，但公式计算较繁琐，特别当实际工程中计算断面数较多时，会大大加大计算工作量，影响工作效率，且容易出错。1.2.2 弹性理论法弹性理论对结构应力分析最早始于19世纪末，由数学家及工程师M.Levy提出。弹性理论最早应用是求解无限楔体在自重和一些边界力作用下的应力分布。以重力坝为例来说明弹性力学应用，在利用无限楔形体求解上游斜率为0 的重力坝剖面时，可考虑坝体的各种集中力，边界力，然后利用叠加法将它们分解或组合成一定标准形式，可通过现有的计算方法，求出应力分布，其结果精度高于材料力学法。但由于往往实际坝体断面为上下游变坡形式或常带有坡角等复杂形式，弹性理论法仍无法完全求出其精确解。通过弹性理论求出重力坝剖面解与材料力学法相比较，可以看出在距坝基1/31/ 4坝高范围内，其分布形式、求解数值结果差别较大，越往上差别越小。该理论可考虑地基条件、荷载条件，得到比材料力学法更精确的解答。弹性理论除了满足材料力学法所满足的平衡条件、边界条件外，还需满足变形相容条件，但往往由于考虑的因素多，实际边界条件的复杂性，造成计算结果过于冗长或根本无法完全考虑实际问题中复杂因素。弹性理论在重力坝应力求解时，没有得到实际应用。1.2.3 数值分析方法有限元法是目前比较流行的一种数值解法。20世纪中叶，经过大量试验研究，学者们提出了有限差分法、拟板法、刚架法等许多数值解法，对坝体应力分布状态有了更全面的掌握，由于计算工程量浩大，而无法普遍使用。随着生产的需要，软件技术发展，有限单元法应运而生。有限元法的出现使弹性理论法进入一个新的局面。有限元法的应用可以更好模拟实际工程问题，通过分析可更好的了解坝体的不同部位应力、位移分布和受力情况。随着软件开发水平不断的发展，许多大型有限元分析软件应运而生，使得有限元法更广泛的应用于实际工程。一般重力坝与地基变位很小，均在线弹性范围内。实际上不论岩石还是混凝土材料，应力在一定范围内本构关系呈线弹性关系，当超过一定范围后，材料将呈非线性变化，出现塑性变形，特别是带有软弱夹层、节理、破碎带的坝基岩体，在荷载作用下出现明显塑性变形；再者材料破坏后，材料不再是连续弹性体。此种情况需采用非线性有限元理论进行考虑。实际计算中，混凝土或岩基的材料参数并不是十分精确；再者非线性有限元需要掌握复杂的本构关系与破坏准则。这是非线性有限元在重力坝工程中未得到广泛应用的主要原因。 采用有限法计算应力时，在坝锺、坝址处往往会出现应力集中现象，单元划分越细，应力集中越剧烈，且应力趋向无穷大，失真现象比较严重。对于中等高度以下并且地址条件比较好的重力坝，其主要是在弹性阶段工作，可按线弹性有限元进行分析。对于高坝或坝址处有软弱夹层等复杂地基，应采用非线性有限元分析。目前规范常以传统材料力学法为依据，对有限元计算结果未做明确规定，世界各国均未制定出相应的应力应变控制标准。因此，有限元法尚处于定性的分析大坝应力、应变阶段，仍不能取代传统的材料应力法。对于重要工程的分析，可通过不同种方法加以比较，进而更全面的了解大坝应力，应变情况，对工程进行客观判断。1.3重力坝抗滑稳定分析发展状况大坝抗滑稳定性是评价大坝安全的重要指标之一。重力坝是悬臂式结构，在各种荷载作用下，靠自重产生的摩擦阻力来维持稳定。目前，比较流行的抗滑稳定分析方法有刚体极限平衡法、有限元分析法和模型试验法等，其中刚体极限平衡法是最基本的、应用时间最长的方法。当遇到重要工程或基础较差需进行深层抗滑稳定分析等复杂地基时，往往需采用另种方法如有限元法或模型试验来进行校核分析，通过不同方法，以对实际问题进行综合分析了解。1.3.1传统刚体极限平衡法传统的刚体极限平衡法是将坝体或坝体与坝基接触面看作一个或若干个滑移面，滑移面上部结构看作整体滑移刚体。该理论使用时间长、原理简单和容易求解计算，得到了许多实际工程反复验证，有完善的控制指标与相应规范，是目前抗滑分析的主要依据，常用的公式有抗剪强度公式和抗剪断强度公式。 抗剪强度公式，又称为“纯摩”公式，该理论认为，大坝抗滑稳定主要靠滑裂面以上结构自重与滑裂面摩擦系数产生的摩擦力来维持稳定。摩擦系数的选取，直接影响到大坝的造价与安全，该参数一般通过实验材料剪断后的实验值，现场实验成果和结合已建类似工程经验等多方面综合考虑确定。由于接触面不可能是完全光滑面，该假定不符合实际，计算出来的抗滑系数有一定的安全储备。20世纪50年代，学者们提出了抗剪断强度公式。该理论认为接触面为胶结面，破坏时只有当胶结面处材料发生断裂、屈服等形式破坏时，才形成滑动通道，这样使得物理概念更为明确，较符合实际情况。1.3.2 数值分析法抗滑稳定分析随着有限元理论和计算机技术的发展，有限元法已广泛应用于坝体抗滑稳定分析。该方法比刚体平衡法更加准确、合理。有限元法可求出坝体及坝基内各点的应力值，可以计算出沿软弱破坏面上的局部抗滑稳定系数，以此了解破坏区的分布范围及其严重程度，更深层次的分析大坝的抗滑稳定。对研究一些软弱夹层带、破碎层带等复杂非线性影响，发挥越来越重要作用。由于现场实验勘测结果的精度有限，非线性问题复杂性，研究中往往带有一定的近似，该计算成果可用于传统方法计算结果的校核和补充。1.3.3 模型实验法模型试验法是用一定比例的实验模型来模拟实际结构的方法，该方法是研究坝体及坝基应力、应变、开裂和失稳等问题的一种重要手段。在模型试验中，模型的几何尺寸，材料的容重、强度、弹模、应变、应力等各种参数与实际物体相应的参数之间应遵循一定的定律，并满足一定的比例关系。在研究复杂的抗滑稳定时，如果想用模型试验来真实地反映地基内的各种情况和性质，往往具有很大的困难性。因为对地质条件的模拟，不仅需要模型材料、岩块的相似，还需能真实反映节理，裂隙组合，断裂层等问题，否则失去试验意义。这种方法比较麻烦，费时，耗资大，常用在重要工程中作为计算结果的校核。1.4研究的目的和意义治理水患、进行水利工程建设是一直以来人们所不断关注的，大约在5000年前人们就开始建造重力坝，它是人类最早使用的一种水坝坝型。近代以来,重力坝工程更是蓬勃发展，在世界各地都能见到它踪影，为人类合理利用自然资源做出了巨大的贡献。 重力坝的应力分析是在坝体断面已初步拟定的情况下进行的。其目的是为了判定坝在运用期和施工期是否满足强度和稳定方面的要求，同时也为研究与设计、施工有关的其他问题（如确定坝体混凝土的标号分区以及在某些部位配置钢筋等）提供依据。重力坝依靠自身重量来维持稳定，所以安全性高就是重力坝设计中最基本、最重要的要求。早期的重力坝断面都比较大，原因是当时对坝体断面的认识是断面越大坝体就越安全。然而，事实上，如果只是一味地加大断面并不能保证大坝的安全。据詹森统计，直到1940年以后水坝失事的情况才很大程度上减少，而这以前全世界已超过15万座水坝失事(其中绝大部分属于小型水坝)。最后，失事数目稳定在水坝总数的1%以内。总结说来，坝基缺陷是水坝失事的原因中最主要的因素之一。1.5本文主要的工作本文以观音岩水电站重力坝为工程实例，运用材料力学法和有限元法分别对重力坝12#坝段进行了应力应变分析，并通过刚体极限平衡法和有限元法分别对坝体进行抗滑稳定分析。分别通过两种不同方法来对重力坝应力及稳定分析，以了解大坝实际运行情况。本文主要工作内容如下： (1)通过查阅大量国内外相关文献，了解重力坝坝的分类、工作特点，并对重力坝应力分析理论、稳定分析理论发展状况进行了总结。 (2) 本文以SOLID45 、SOLID65 6 面体 8 节点空间单元为例，推导了有限元法求解基本过程，为有限元分析提供了理论依据。 (3) 分别采用材料力学法与有限元法对12#坝段进行了应力分析，通过两种方法对比，分析两种计算结果的区别，并对观音岩水电站重力坝12#坝段应力状态进行综合评价。 (4)分别采用刚体极限平衡法、有限元法对观音岩水电站重力坝12#坝段进行抗滑稳定分析，结合两种不同方法分析比较，对重力坝稳定状态进行综合评价。1.6本章小结本章简要介绍了观音岩水电站工程概况，了解了重力坝设计理论的发展状况，阐述了重力坝应力、稳定分析理论的发展，说明了材料力学法、刚体极限平衡法和有限元法理论的优点和不足，结合观音岩水电站重力坝12#坝段这一工程实例，提出了本课题的研究目的和意义。第2章FLAC3D数值分析方法简介及分析过程2.1数值分析方法概述在数值方法中，分为两大类，分别为：连续介质方法，包括有限元、有限差分、边界元、无单元等；非连续介质方法，包括离散元、颗粒元、流形元等。由于计算机计算性能的发展，数值分析方法得到了长足的发展。立足于不同的数学模型，得出了不同的数值分析方法，有限元法（FEM）由线性到非线性的发展，为岩土工程解决了极大的难题。由连续力学到非连续力学的发展出现了离散元（DEM）和不连续变形分析方法（DDA）及流形元法（Manifold Method）的进展，为不连续介质力学的发展做出了很大贡献。边界元法（BEM）、无单元法（EFM）、快速拉格朗日法（FLAC）的出现又为大变形或有限变形分析做出了贡献。有限元法（FEM）发展到今天，已成为求解复杂的岩石力学及岩土工程问题的有力工具，并已为工程科技人员所熟悉。在求解弹塑性及流变、动力、非稳态渗流等时间相关性问题，以及温度场、渗流、应力场的耦合等复杂的非线性问题中其效能已成为岩石力学领域中应用最广泛的数值分析手段。有限元分析中最基本的思想就是单元离散，即将求解域剖分为若干单元，把一个连续的介质换成为一个离散的结构物，然后就各单元进行分析，最后集成求解整体位移（基于最小势能原理的位移法）。为了方便而有效地离散复杂的岩体构造及建于其上的建筑物（如水利水电工程中的大坝结构），用各种实体单元、夹层单元及无穷元，分别模拟岩石及混凝土的连续区、断层及结构面，以及无穷域。就数学概念来说，有限元法是通过变分原理或加权余量法和分区插值的离散化处理把基本的支配方程转化为线性代数方程，把求解域内的连续场函数转化为求解有限个离散点处的场函数值。显然，这种离散化的处理是一种近似，因而只有当单元划分得适当时，才能保证满意的求解精度。有限单元法是建立在固体变分原理基础之上的。用有限元进行分析时，首先将被分析物体离散成为许多个小单元，其次给定边界条件、荷载和材料特性，再求解线性或非线性方程组，得到位移、应力、应变和内力等结果。所得到的结果不是准确解，而是近似解，但足以满足大多数实际工程的需要，因而成为行之有效的工程分析手段。有限元法可以分析形状十分复杂、非均质的各种实际工程结构，可以在计算中模拟各种复杂的材料本构关系、荷载和边界条件；可以进行结构的应力和动力分析，结合前后处理技术可以进行方案的优化和选择，并且可以迅速地用图形直观地表示出来。2.2三维快速拉格朗日分析的数学模型三维快速拉格朗日法是一种基于三维显示有限差分法的数值分析方法，它可以模拟岩土或其他材料的三维力学行为。三维快速拉格朗日分析将计算区域划分为若干四面体单元，每个单元在给定的边界条件下遵循指定的线性或非线性本构关系，如果单元应力使得材料屈服或产生塑性流动，则单元网格可以随着材料的变形而变形，这就是所谓的拉格朗日算法，这种算法非常适合于模拟大变形问题。三维快速拉格朗日分析采用了显示有限差分格式来求解场的控制微分方程，并应用了混合单元离散模型，可以准确的模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。三维快速拉格朗日分析在求解中使用如下3种计算方法：离散模型方法。连续介质被离散为若干六面体单元，作用力均被集中在节点上。有限差分方法。变量关于空间和时间的一阶导数均用有限差分来近似。动态松弛方法。由质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。2.2.1空间导数的有限差分近似快速拉格朗日分析采用混合离散方法，将区域离散为常应变六面体单元的集合体，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体，应力、应变、节点不平衡力等变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均。这种方法既避免了常应变六面体单元常会遇到的位移剪切锁死现象，又使得四面体单元的位移模式可以充分适应一些本构的要求，如不可压缩塑性流动等。如图2-1所示一四面体，节点编号为1到4，第n面表示与节点n相对的面，设其内一点的速率分量为，由高斯公式得：图2-1 四面体单元的面和节点 vvi，jdV=SvinjdS (2-1)其中V为四面体的体积，S为四面体的外表面，nj为外表面的单位法向向量分量。对于常应变单元，vi为线性分布，nj在每个面上为常量。对式(2-1)积分得：Vvi，j=f=14vifnjfSf (2-2)式中，上标f指面f的相关变量值，vi指i速度分量的均值。若速度呈线性变化，则：vif=13l=1，lf4vil (2-3)上标l指节点l的值。将上式代入式(2-1)，有：Vvi，j=13l=14vill=1，lf4njfSf (2-4)在式（2-1）中，若vi=1，应用高斯法则可得：f=14njfSf=0 （2-5）所以，式（2-4）两边同除以V，则有：vi，j=13Vl=14vilnjfSf （2-6）而应变速率张量的分量形式为：ij=-16Vl=14vilnjl+vjlnilSl （2-7）2.2.2节点运动方程一定区域内，静力平衡问题可通过以下的平衡方程求解得到：ij，j+Bi=0 （2-8）式中：为介质密度，Bi=bi-dvidt，bi为介质单位质量的体积力。根据虚功原理，作用于单个四面体上的节点力fi（i=1,4）与四面体应力和等效体力相平衡。引入节点虚速度vi（它在四面体中产生线性速度场v和常应变速率），则节点力Fi和体力B产生的外力功功率等于内部应力ij产生的内力功功率。外力功功率可表示为：E=n=14vinfin+VviBidV （2-9）而内力功功率：I=VijijdV （2-10）由式（2-7），对常应变速率的四面体有：I=-16l=14vilijnjl+vjlijnilSl （2-11）应力张量是对称张量，定义矢量Ti：Til=ijnjlSl （2-12）则：I=-13l=14vilTil （2-13）式（2-8）代入式（2-9），有：E=n=14vinfin+Eb+EI （2-14）Eb和EI分别为体力bi和惯性力所做的外力功功率。若四面体内力功bi为常数，则有：Eb=biVvidV （2-15）EI=-VvidvidtdV （2-16）根据有限差分近似，速度场在四面体内线性变化。为描述它，引进一个参考坐标系（它的坐标原点在四面体的中心上），则有：vi=n=14vinNn （2-17）式中Nn（n=1-4）为一线性函数：Nn=c0n+c1nx1+c2nx2+c3nx3 （2-18）其中，c0n、c1n、c2n、c3n（n=1-4）为下述方程的解：Nnx1j，x2j，x3j=nj （2-19）式中，nj是克罗内克尔增量（Kronecker delta）。通过中心点的定义，所有形如VxjdV的积分均为0，将式（2-18）、式（2-17）代入式（2-14）得：Eb=bin=14vinc0nV （2-20）由克雷姆定律，解式（2-19）得：c0n=14 （2-21）将上式代入式（2-20），有：Eb=n=14vinbiV4 （2-22）同理，将式（2-17）代入式（2-16）得到：EI=-n=14vinVNndvidtdV （2-23）将式（2-22）和式（2-23）代入式（2-14）：E=n=14vinfin+biV4-VNndvidtdV （2-24）对任何虚速度，外虚功率E等于内虚功率I：-fin=Tin3+biV4-VNndvidtdV （2-25）为不变量，则上式可写为：VNndvidtdV=V4dvidtn （2-27）用假想的节点质量mn代替上式中的质量V4：则式（2-25）可写为：-fin=Tin3+biV4-mndvidtn （2-28）对于等效体系，可以建立平衡状态，要求在每个节点上静态等效荷载之和为零。可以写出全部节点上牛顿定律表达式：Fil=Mldvidtl l=1，nn （2-29）式中，nn介质中的所有的节点总数，节点质量定义为：Ml=ml （2-30）不平衡力Fil定义为：Fil=Tin3+biV4l+Pil （2-31）当介质达到平衡时，不平衡力等于0。2.2.3增量形式的本构方程快速拉格朗日分析中，假定时间t内速度为常数，增量形式的本构方程可表示为：ij=Hij*ij，ijt （2-32）式中，ij称为共转应力增量，Hij*为一给定的函数。共转应力速率张量ij等于给定参考系的介质内一点应力的偏导数和以瞬时角速度的转动，数学表达式为：ij=dijdt-wijij+ijwkj （2-33）式中，w为转动速率张量。利用有限差分方程，可以得到转动速率张量的分量形式：wij=-16Vl=14vilnjl-vjlnilSl （2-34）式中符号同前。2.2.4时间导数的有限差分近似由本构方程（式（2-32）和变形速率与节点速率之间的关系（式（2-7），式（2-26）可表示为一般的差分方程：dvildt=1MlFilt，vi1，vi2，vi3，vipl，k l=1，nn （2-35）式中，l是指在计算过程中全局节点l节点速度值的子集（式（2-29）。在时间间隔t中实际节点的速度假定是线性变化的，式（2-35）左边导数用中心有限差分估算。vilt+t2=vilt-t2+1MlFilt，vi1，vi2，vi3，vipl，k（2-36） 类似地，节点的位置也用中心有限差分进行迭代：xilt+t=xilt+tvilt+t2 （2-37）因此，节点位移也有如下关系：uilt+t=uilt+tvilt+t2 （2-38）2.2.5阻尼力为使运动方程获得静态或准静态（非惯性）解，快速拉格朗日分析的静力分析中，在式（2-29）中加入非粘性阻尼力。则式（2-29）变为：Fil+lil=Mldvidtl l=1，n （2-39）式中Fil为阻尼力，lil=-Filsign(vil)，为阻尼系数，其默认值为0.8。sign(y)=+1，if y0-1，if y00，if y=0 （1-40）2.2.6有限差分法岩土工程结构的数值解是建立在满足基本方程（平衡方程、几何方程、本构方程）和边界条件下推导的。由于基本方程和边界条件多以微分方程的形式出现，因此，将基本方程近似用差分方程（代数方程）表示，把求解微分方程的问题改换成求解代数方程的问题，这就是所谓的差分法。差分法由来已久，但差分法需要求解高阶代数方程组，只有计算机的出现，才使该法得以实施和发展。有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以泰勒级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以为显格式、隐格式、显隐交替格式等。目前常用的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用泰勒级数展开方法。其基本的差分表达式主要有4种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分形式的组合，可以组合成不同的差分计算格式。2.2.7 FLAC3D软件FLAC3D（3-D Fast Lagrangian Analysis Code）是由美国ITASCA公司开发的。目前，FLAC有二维和三维计算程序两个版本，二维计算程序V3.0以前的为DOS版本，V2.5版本仅仅能够使用计算机的基本内存（64KB），所以，程序求解的最大节点数仅限于2000个以内。1995年，FLAC2D已升级为V5.0版本，其程序能够使用扩展内存。因此，大大发展了计算规模。FLAC3D是一个三维有限差分程序，目前已发展到V3.0版本。并且其推出的FLAC SLOPE有了Windows界面。FLAC3D的输入和一般的数值分析程序不同，它可以用交互的方式，从键盘输入各种命令，也可以写成命令集文件，类似于批处理，由文件来驱动。因此采用FLAC程序进行计算，必须了解各种命令关键词的功能，然后，按照计算顺序，将命令按先后依次排列，形成可以完成一定计算任务的命令文件。FLAC3D是二维的有限差分程序FLAC2D的扩展，能够对土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析。调整三维网格中的多面体单元来拟合实际的结构。单元材料可采用线性或非线性本构模型，在外力作用下，当材料发生屈服流动后，网格能够发生相应的变形和移动（大变形模式）。FLAC3D采用的显示拉格朗日算法和混合-离散分区技术能够非常准确地模拟材料的塑性破坏和流动。由于无需形成刚度矩阵，因此，基于较小内存空间就能够求解大范围的三维问题。FLAC3D采用ANSI C+语言编写。FLAC3D有以下几个优点：对模拟塑性破坏和塑性流动采用的是“混合离散法”。这种方法比有限元法中通常采用的“离散集成法”更为准确、合理。即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍。采用了一个“显示解”方案。因此，显示解方案对非线性的应力-应变关系的求解所花费的时间，几乎与线性本构关系相同，而隐式求解方案将会花费较长的求解时间求解非线性问题。而且，它没有必要存储刚度矩阵，这就意味着采用中等容量的内存可以求解多单元结构；模拟大变形问题几乎并不比小变形问题多消耗更多的计算时间，因为没有任何刚度矩阵要被修改。当然，它也存在以下几个不足之处：对于线性问题的求解，FLAC3D比有限元程序运行得要慢；因此，当进行大变形非线性问题或模拟实际可能出现不稳定问题时，FLAC3D是最有效的工具。用FLAC3D求解时间取决于最长的自然周期和最短的自然周期之比。但某些问题对模型是无效的。FLAC3D中为岩土工程问题的求解开发了特有的本构模型，总共包括10种材料模型：开挖模型null；3个弹性模型（各向同性、横观各向同性和正交各向同性弹性模型）；6个塑性模型（Drucker-Prager模型、Mohr-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型和修正的cam模型）。FLAC3D网格中的每个区域可以给以不同的材料模型，并且还允许指定材料参数的统计分布和变化梯度。还包括了节理单元，也称为界面单元，能够模拟两种或者多种材料界面不同材料性质的间断特性。节理允许发生滑动或分离，因此可以用来模拟岩体中的断层、节理或摩擦边界。FLAC3D中的网格生成器gen，通过匹配、连接由网格生成器生成局部网格，能够方便地生成生成所需要的三维结构网格。还可以自动生成交叉结构网格（比如说相交的巷道），三维网格由整体坐标系xyz系统所确定，不同于FLAC程序由行列方式确定。这就提供了比较灵活的产生和定义三维空间参数。边界条件和初始条件：定义方式与FLAC相同。在边界区域可以指定速度（位移）边界条件或应力边界条件。也可以给出初始应力条件，包括重力荷载以及地下水位线。所有的条件都允许指定变化梯度。FLAC3D还包括了模拟区域地下水流动、空隙水压力的扩散以及多空隙固体和在空隙内粘性流动流体的相互耦合。流体被认为是服从各向同性的达西定律。流体和空隙固体中的颗粒是可变形的，将稳态流处理为紊态流可以模拟非稳态流。同时能够考虑固定的空隙压力和常流的边界条件，也能模拟源和井。流体模型能够和结构的力学分析独立进行。FLAC3D做计算分析的一般步骤：与大多数程序采用数据输入方式不同，FLAC3D采用的是命令驱动方式。命令字控制着程序的运行。在必要时，尤其是绘图，还可以启动FLAC3D用户交互式图形界面。为了建立FLAC3D计算模型，必须进行以下三个方面的工作：有限差分网格；本构特性和材料特性；边界条件与初始条件。完成上述工作后，可以获得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后进行工程开挖或改变边界条件来进行工程的响应分析，类似于FLAC3D的显示有限差分程序的问题求解。与传统的隐式求解程序不同，FLAC3D采用一种显示的时间步来求解代数方程。进行一系列计算步后达到问题的解。在FLAC3D中，达到问题所需的计算步能够通过程序或用户加以控制，但是用户必须确定计算步是否已经达到问题的最终的解。FLAC3D程序计算方法是美国ITASCA咨询公司Cundall等人提出的显示有限差分法而编制的有限差分软件，具有很强的分析功能，主要特点为：快速朗格朗日求解理论。为物理不稳定系统提供稳定解，时间追赶法辅助用户洞悉节点及网格单元变形历史和应力路径，特别地，能够真实展示与路径相关的材料强度峰后行为。丰富的本构模型库。FLAC3D为连续介质力学分析提供了极为丰富的本构模型库，材料模型的选用取决于计算模式的不同：l 静力分析模式下，FLAC3D提供了12款成熟本构模型，分为弹性、弹性/塑性，如弹性、Mohr-Coulomb理想弹塑性、遍布节理、双屈服、应变软化、修正剑桥和Hoek-Brown模型等；l 特别地，FLAC3D更是提供诸如流体模型、温度模型、蠕变模型和动力模型满足复杂受力条件对介质力学行为影响研究的需要。强大的分析功能。除基本模块所提供的常规分析功能外，还包括若干可选功能模块，包括温度分析、蠕变分析、动力分析和基于C+编译环境的自定义本构模型模块。非连续特征分析能力。针对介质非连续特点（如岩体地质结构面），程序亦提供了界面单元（interface）以模拟沿不连续结构面的滑移/分离响应，如断层、节理和摩擦边界的力学响应。通过三维介质的离散，使所有外力与内力集中于三维网格节点上，进而将连续介质运动定律转化为离散节点上的牛顿定律；时间与空间的导数采用沿着有限空间与时间间隔线性变化的有限差分来近似；将静力问题当做动力问题来求解，运动方程中的惯性项用来作为达到所需静力平衡的一种手段。2.3用FLAC3D分析问题过程FLAC3D是基于命令驱动模式的软件，因为大多数的分析都要用到输入文件。命令语句控制着分析的进程。当然用户交互式的图形控制界面在某些时候还是可以排上用场，比如在控制出图的时候会相对方便一些。这一部分主要讲解用FLAC3D进行简单计算所需要的基本命令的使用方法。以低粘聚力土体中沟渠的开挖问题为例，来讲述解决问题的基本过程。这个过程可以分为几个明显的阶段，每个阶段都分开来论述，每一个部分中都详述了所用到的各种命令，以使得用户更明白每个过程的操作。要建立一个可以用FLAC3D来模拟计算的模型，必须要做以下三步工作：建立模型的有限差分网格；定义本构模型和赋材料参数；定义边界条件、初始条件。由网格来定义所要模拟的几何空间。由本构模型和材料参数来限定模型对于外界扰动做出的变化规律（比如开挖引起的变形反应）。由边界条件和初始条件来定义模型的初始状态（比如说模型在发生变化或扰动前的稳定状态）。做好了以上三步工作，就可以进行模型初始平衡状态的计算了。接着对模型做一些变动（比如开挖或者改变边界条件），然后再对改动后的模型进行计算。FLAC3D为采用显示解法的软件，它的实际求解过程不同于常规的隐式解法。FLAC3D是采用显示时间步推的方法来求解代数方程组的，通过一些时间步的计算，才会得到所要的计算结果。完成计算所需要的时间步可以由软件自动控制，也可以人为地指定计算步数。但最后还是需要用户自己来判断进行了这些时间步的计算后模拟的问题是否已经得到了最终所要的解。第 3 章 重力坝材料力学法应力分析3.1工程基本资料3.1.1 工程地质工程地质情况如图3-1所示，其中，F1表示断层，“砾”代表砾岩，“粉”代表粉砂岩，“砂”代表砂岩，“泥”代表泥岩。图3-1 12#剖面地基岩体力学参数剖面图3.1.2 水文气象特征观音岩水电站坝址以上金沙江河道长2703km，落差4385m，河道平均坡降1.62，控制流域面积约26万km2，多年平均流量1880m3/s。金沙江流域属典型的季风气候区，每年510月，主要受西南季风控制，同时也受东南季风影响，水汽充沛、降雨集中，为本流域雨季或汛期。年平均降水量呈由南向北递减的总趋势。流域内暴雨较小，其分布中、下游较上游大。3.1.3 材料计算参数混凝土容重：24 KN/m3；水容重：10 KN/m3；淤沙高程：1082.000m；淤沙容重：9.5KN/m3；淤沙内摩擦角：24。对于坝基岩体，其相应的计算参数参见表3-1。表3-1剖面材料参数表岩类抗剪断强度变形模量泊松比天然容重fC(MPa)E0(GPa)（g/cm3）泥岩0.810.420.30.32.55砂岩1.131.3015.190.222.60砾岩0.630.220.880.32.20粉砂岩0.981.007.070.242.60断层0.360.0420.30.32.0混凝土（建基面）1.52.025.50.1672.403.1.4 特征水位表3-2 特征水位计算工况上游水位下游水位备注正常蓄水位1134.001017.24一台半机发电尾水位校核洪水位1136.491042.763.1.5 工程规模等级 观音岩水电站总库容约20.72亿m3，装机规模为5600MW=3000MW，参照DJ51802003水电站枢纽工程等级划分及设计安全标准中水利枢纽工程的分等指标（见表3-2）可得，观音岩水电站枢纽工程工程等级为一级，工程规模为大（1）型。表3-3 水电站枢纽工程的分等指标工程等级工程规模水库总库容（亿m3）装机容量（MW）一大（1）型101200二大（2）型 10 1.0 1200 300三中型 1.0 0.10 300 50四小（1）型 0.10 0.01 50 10五小（2）型0.01103.1.6 计算工况及荷载组合根据规范、设计资料以及实际情况，计算工况以及荷载组合参见表3-4表3-4 计算工况及荷载组合计算工况水位荷载自重静水压力扬压力淤沙压力浪压力地震荷载1库空2正常蓄水位3校核洪水位3.2 荷载计算 3.2.1 坝体自重 坝体材料容重：24 KN/m3；坝基材料容重：26 KN/m3；水容重：10 KN/m3；重力加速度：10m3/s。3.2.2 静水压力 静水压力垂直作用于结构表面，可分解包括水平力和垂直自重。水平作用于机构表面的静水压强可表示为：p=H （3-1）式中：p静水压强(KPa ）； 水容重； H计算点处作用水头(m)；3.2.3 淤沙压力淤沙压力是挡水建筑物由于泥沙淤积，而作用在挡水建筑物上的一种土压力，作用在单位长度挡水结构上的水平淤沙压力表达式为：psk=12sbhs2tan245-2 （3-2）式中：psk淤沙压力；sb淤沙的浮重度，大小为9.5 KN/m3；hs坝前泥沙淤积厚度；淤沙的内摩擦角，大小为24；3.2.4 扬压力 扬压力包括渗透压力和浮托力，浮托力为下游水深所引起。渗透力是在上下游水位差作用下，经过时间的推移，在坝体，坝基形成的稳定的渗透场，引起的水压力。影响扬压力分布和大小的因素很多如坝体形状，坝体所采取的防渗、排水设施，坝底与接触面的空隙分布、坝基本身裂隙等有关。工程中一般采用防渗、排水措施减少扬压力。对于本工程，坝基除设置上游主帷幕及主排水孔外，还设置了下游帷幕、副排水孔及抽排系统，基底的扬压力见图3-2。图3-2 坝基扬压力分布图其中a=8.0m，b=11.25m，c=26.25m。3.2.5 浪压力由于设计资料有限，本次毕业设计不考虑浪压力。3.2.6 地震荷载 由于设计资料有限，本次毕业设计不考虑地震荷载。3.3材料力学法简介采用传统的材料力学法，先计算各荷载在计算截面作用的总垂直力Fy，总水平力Fx以及总力矩M，然后再计算在这些荷载的作用下，截面上下游边缘产生的垂