【《两河口心墙堆石坝枢纽布置、大坝剖面设计和渗流坝坡稳定性分析》16000字】

绪论1.1论文的研究背景由于心墙土石坝具有对坝址所在的地质条件适应性强、施工工艺易于掌握且造价便宜、坝体材料可就地取材来源丰富以及良好的抗震性能等优点而被广泛利用，其建设历史悠久，有着丰富的实践经验。随着堆石体材料性质的深入研究和施工工艺的不断进步，堆石坝工程建设技术逐渐提升，心墙堆石坝成为当下高坝建设的最常用的坝型之一[1]。堆石坝坝体土石料应力应变关系的非线性会随着大坝增高而更为复杂，致使分析难度也随着增加。因此，对高土石坝的应力变形特性进行分析研究十分必要,可以利用有限元分析方法对坝体在施工、蓄水及运行过程中应力应变关系进行三维非线性分析[2]。近年来我国大力推进西部大开发和西电东送的发展战略，我国西南部的水利水电工程建设得到了迅速的发展。在这样的战略背景之下，一批高堆石坝工程项目陆续上马，开始规划和建设，如坝高为314m的双江口心墙堆石坝，356m的其高心墙堆石坝、261.5m的糯扎渡心墙堆石坝以及高达295m的雅砻江两河口心墙堆石坝等工程[3]。超高心墙堆石坝需要防止心墙体发生较大变形，因为心墙较大的变形可能会导致防渗结构产生裂缝，裂缝的产生就会导致防渗结构防渗能力大大降低，甚至坝体的稳定也会受到影响[4]。因此有必要利用三维非线性的有限元计算方法对坝体在施工、蓄水过程中的应力变形特性进行分析分析坝体和心墙的应力变形特性，从而保证高心墙坝的安全稳定运行[5]。本论文以两河口心墙堆石坝为背景，通过对有限元和心墙堆石坝的设计的基本理论和应力变形计算分析研究的学习，采用有限元方法，按照设计拟定的实际施工程序和蓄水情况施加荷载，使用软件ANSYS进行建模以及使用ABAQUS进行分析计算，对两河口心墙堆石坝的应力变形进行分析，论证其运行稳定性，为设计提供参考。1.2国内外研究进展近年来，通过不断成熟的堆石坝工程设计技术与技术创新、施工机械的不断革新以及先进的施工管理技术，人们得以建成坝高达300m级的心墙堆石坝。如今坝工技术仍在飞速持续的发展进步，将会有更多的，地质条件更为复杂的并且工程规模更大的高心墙堆石坝工程开始规划与建设。国内外的研究人员都在心墙堆石坝应力变形分析方面开展了大量的研究，进行了大量的工作。研究人员通过大量的试验和对试验结果的理论分析，在对堆石料本构关系的研究方面取得了长足的进步，各种不同的堆石料本构模型被研究人员们陆续提出。而大量的研究与实践表明，堆石料的变形特性可以被非线性弹性模型和弹塑性模型较好地反映。自此有限元分析的方法在土石坝的应力应变分析中开始被大量采用[6]。我国自1985年以来，经过多年的不断耕耘，已完成了在心墙堆石坝应力应变分析方面从探索道路到慢慢走向成熟的关键过程，现在我国在堆石坝应力应变分析的研究方面，已经走在了世界前列。目前，在我国，将土石坝的应力变形看作是是典型的位移场和渗流场耦合问题，其设计方法仍处在半理论半经验阶段[7]。以有限元为代表的数值分析方法被广泛地应用在对土石坝地应力应变分析中，其中堆石体的本构模型在有限元分析中十分关键。现在国内已经有许多堆石体本构模型，这些本构模型也被经常采用，主要有邓肯-张非线性弹性模型、四川大学KG模型、清华非线性解耦KG模型以及沈珠江双屈服面弹塑性模型等本构模型。在高塑性心墙及软土地基的数值计算中也有使用修正剑桥模型的计算程序[8]。这些模型各具优缺点，但总体来说已能较好地模拟堆石坝的真实工作情况。当下对高土石坝工程中的本构模型的研究仍然很活跃。例如，有关在高围压以及复杂应力路径下堆石料的变形特性，对考虑颗粒破碎及复杂应力条件下的堆石体应力应变模型改进，堆石体长期破损特性试验及计算模型等[9]。在土石坝的变形控制方面，我国近年来在高土石坝工程中的实践表明，在采用相同的变形控制标准的条件下，土石坝变形控制难度会随着坝高的提高而大大增加，有时甚至需要花费巨大的代价。并且当下高土石坝变形控制的标准是笼统的，并没有与具体的大坝坝体破坏形式相关联，因此有必要针对坝体裂缝等可能的表现行为探讨大坝坝体的变形控制标准。1.3论文的主要内容及研究方法1.3.1主要内容本论文以两河口心墙坝为背景，采用有限元分析方法，先整体建模，根据施工填筑过程和坝体的剖面图对单元进行划分，并定义单元的材料参数。然后按照设计拟定的实际施工程序和蓄水情况施加荷载，利用邓肯-张E-B模型对两河口心墙堆石坝在竣工期、蓄水期以及运行期进行计算，分析大坝的应力变形特性，论证其运行稳定性，为设计提供参考。其中主要内容主要包括常规设计和专题计算两部分：（1）常规设计部分：根据所给资料，规划大坝枢纽布置，设计大坝体型。常规渗流计算是利用水力学法计算渗透流量和确定浸润线位置，为下游坝坡稳定性计算提供计算依据，然后利用瑞典圆弧法分析坝体的抗滑稳定性。这里采用上游正常蓄水位与下游相应的最低水位作为计算工况。（2）专题计算部分：根据两河口心墙堆石坝已有的工程资料，选取坝典型剖面，在CAD中对典型坝剖面进行一定的简化并划分网格，然后应用有限元软件ANSYS将坝体的三维模型建立出来，从而得到坝体的三维有限元模型，先计算稳定渗流，再考虑静力、流变两种情况，按照实际工程建设施工中坝体的填筑、蓄水过程和运行情况分别对使用ANSYS软件建成的堆石坝三维有限元模型施加荷载，然后进行计算分析，得到应力变形计算成果。然后再单独对心墙，根据计算结果确定心墙在静力、流变种情况下施工、蓄水和运行期的应力变形。最后对计算结果进行检验，看计算结果是否满足设计所允许的范围，并与国内外实测的数据相比较，分析其变形规律是否一致，检验成果是否满足大坝安全稳定运行的要求，从而保证大坝稳定的运行。1.3.2研究方法本论文通过对文献的查阅研究，依据相关规范中提供的计算方法，对堆石坝采用有限元方法进行计算，对计算结果进行分析和研究。具体方法为：利用CAD和ANSYS软件对两河口心墙堆石坝进行三维建模，利用ABAQUS软件并根据邓肯-张E-B模型对两河口心墙堆石坝做有限元分析计算，利用中点增量法模拟施工步骤的实现过程，一直到坝体竣工并且蓄水过程完毕。根据计算所得的应力应变论证坝体运行的稳定性。1.3.3技术路线根据两河口心墙堆石坝工程提供的资料，用ANSYS和ABAQUS软件为建模计算分析的工具，采用有限元方法，按照设计拟定的实际施工程序和蓄水情况施加荷载，通过对两河口心墙堆石坝的应力变形的分析，论证其运行稳定性，为设计提供参考。1.3.4可行性分析有限单元法建模可以将复杂的几何体通过相应的有限单元进行简化，有限单元之间还相互关联，因此可以建模分析计算出结果，使复杂的工程问题简单化[10]。邓肯E-B模型具有许多优点，首先该模型能够对堆石的变形性状进行较好地模拟，而且该模型参数的测定较为简单并且当下已经有较为成熟的经验[11]。邓肯E-B模型不仅公式十分简单，还是堆石坝中使用频率最高的模型之一。通过三轴试验研究结果可以看出，邓肯模型能够较好反映土体应力应变的非线性特性[12]。2两河口心墙堆石坝工程资料2.1大坝基本情况两河口水电站工程地理位置处在四川省甘孜州雅江县境内的雅砻江干流上，两河口水电站工程坝址处的控制流域面积为6.57万km2，坝址处的河流多年平均流量为666m3/s。此工程为一等大（1）型工程，该工程中永久性主要建筑物如挡水、泄洪、引水及发电等建筑物都是为1级建筑物。两河口水电站工程的主要开发任务是发电，同时还兼具防洪的任务。工程采用坝式开发，水库正常蓄水位高程2865m，水库总库容为107.67亿m3，调节库容65.6亿m3，该工程具有多年调节能力。两河口大坝为砾石土心墙堆石坝，坝顶高程2878.00m，坝顶宽度为16.00m，大坝心墙在河床部位的底部高程为2582.00m，坝体可根据填筑材料的不同划分为四大区，即防渗体、反滤层、过渡层和堆石区四大区。其中防渗体采用砾石土直心墙型式，坝壳采用堆石填筑，反滤层以及过渡层设在心墙与上、下游堆石之间。2.2施工填筑过程根据施工组织设计，两河口大坝坝体填筑工期约为69个月，计划于2016年7月开始浇筑，于2022年3月完成大坝填筑施工。将大坝主体施工分为六期。大坝上下游堆石、上下游过渡层、上下游反滤层填筑过程如表2.1~2.4所示：表2.1大坝反滤层及心墙填筑过程分期期初高程（m）期末高程（m）填筑时间时长（天）一期258026052016.10.01~2017.03.31182二期260526202017.04.01~2017.07.31122三期262026692017.08.01~2018.03.31243四期26692718.52018.04.01~2019.02.28334五期2718.527842019.03.01~2020.04.30427六期278428762020.05.01~2022.03.31700合计2581.528762016.10.01~2022.03.312008表2.2大坝上下游过度料填筑过程分期期初高程（m）期末高程（m）填筑时间时长（天）一期258026052016.10.01~2017.03.31182二期260526232017.04.01~2017.07.31122三期262326692017.08.01~2018.03.31243四期26692718.52018.04.01~2019.02.28334五期2718.527842019.03.01~2020.04.30427六期278428762020.05.01~2022.03.31700合计2581.528762016.10.01~2022.03.312008表2.3大坝上游堆石料填筑过程分期期初高程（m）期末高程（m）填筑时间时长（天）一期259626002017.01.01~2017.03.3190二期260026392017.04.01~2017.07.31122三期263926692017.08.01~2018.03.31243四期266927382018.04.01~2019.02.28334五期273827862019.03.01~2020.04.30427六期278628672020.05.01~2021.12.31640合计258028672016.10.01~2021.12.311826注：大坝上游堆石已于2017年1月填筑至高程2596m。表2.4大坝下游堆石料填筑过程分期期初高程（m）期末高程（m）填筑时间时长（天）一期26062629.52016.07.01~2017.03.31274二期2629.526402017.04.01~2017.07.31122三期264026732017.08.01~2018.03.31243四期267327272018.04.01~2019.02.28334五期272727862019.03.01~2020.04.30427六期278628672020.05.01~2021.12.31640合计260628672016.10.01~2021.12.312010注：大坝下游堆石已于2016年7月填筑至高程2606m。2.3蓄水方案两河口工程的蓄水阶段划分为初期蓄水期和初期运行期，结合导流规划，两河口蓄水时段划分及蓄水规划如下：表2.5蓄水方案蓄水阶段起始时间末时间初高程（m）末高程（m）蓄水量（亿m3）初期蓄水期2020.11.012020.11.112605.52678.53.0862020.11.122021.05.312678.52686.52021.06.012021.08.072686.52785.0（死水位）31.9初期运行期2022.07.092022.07.302785.02800.09.12022.08.092022.09.082800.02815.010.42022.09.182022.10.112815.02830.011.82022.10.212022.10.312830.02835.54.772022.11.102022.12.102835.52850.02022.12.202023.01.202850.02865.0注：考虑到两河口为高坝大库，建议2023年继续控制蓄水速率，其蓄水过程与2022年一致。后续运行可按水库调度规程运行。3常规设计3.1枢纽总体布置（1）挡水建筑物（堆石坝）：大坝建于主河床，布置呈直线。（2）泄水建筑物（溢洪道）：采用洞式溢洪道，建于左岸山体内。（3）水电站建筑物：电站厂房和开关站等布置在右岸，在开挖的基岩上布置厂房，将开关站布置在厂房旁边。（4）施工导流洞及水库放空洞：洞式溢洪道和导流洞均位于在左岸，布置在山体之中。枢纽总体布置见图3.1。图3.1枢纽布置图3.2大坝剖面设计堆石坝的剖面设计主要是设计出大坝的典型剖面图，在设计过程中需要确定坝顶高程、坝顶宽度、上下游坝坡及心墙的断面尺寸。大坝的剖面设计图见图3.2。图3.2大坝剖面设计图3.2.1坝顶高程大坝坝顶高程的确定取决于水库静水位与超高，两者加起来即为坝顶高程。其中坝顶超高的确定主要依据最大波浪在坝坡上的爬高、安全加高以及最大风壅水面高度。本设计中坝顶高程取为2878.00m。3.2.2坝顶宽度本设计中，综合考虑运行、施工、构造、交通和人防等方面的要求后，确定堆石坝坝顶宽度为16m。3.2.3上下游坝坡坡度对心墙坝，上下游坝坡一般在1:2~1:4之间选取，而且一般上游坝坡的坡度比下游坝坡的稍缓。沿坝坡每一定高度改变一次坝坡，并设置一级马道。本设计中，最大坝高约为298m。本设计中，大坝坝高较高，需要布置多级马道，故大坝下游坝坡采用了五级变坡，第一级到第四级马道的高程依次为2640.0m、2700.0m、2760.0m、2820.0m。每级马道宽取为5.0m。大坝下游坝坡坡度从坝顶至坝趾依次取为1:1.6；1:1.8；1:1.6；1:1.8；1:2.0。上游采用两级变坡，马道的高程为2790.0m，坡度为1:1.8；1:2.0。3.2.4心墙断面尺寸（1）心墙位置：心墙位于大坝断面的中心线。（2）心墙顶宽及坡度：确定心墙顶部的水平宽度需要综合考虑多种因素，首先应当考虑满足控制渗流量和渗透比降的要求，且由于施工填筑过程中会使用机械碾压，因此心墙顶宽也应该满足施工机械碾压的要求，依据规范心墙顶宽一般不能小于3m。本设计中取心墙的顶宽为6m。通常，心墙体上、下游坡度取为1:0.15~1:0.3。在本设计中取上下游坡度为1:0.2。（3）心墙顶部高程：为了防止漫顶，在设计时通常要使心墙顶部在水库静水位以上存在超高。在正常运用情况下，心墙顶部超高取0.3~0.6m，本设计中取超高为0.5m，则心墙顶部高程取2865.5m。在非常运用情况下，心墙的顶部高程则不能低于此情况下的静水位。综合考虑这两种情况，最终心墙顶部高程取为2876m。（4）心墙底宽：心墙的底部宽度首先取决于防渗。为了满足防渗要求，心墙底宽不宜小于水头的四分之一，本工程中上下游最大作用水头为249.11m，故心墙厚应当满足T≥62.28m。本设计中心墙底部宽度为140m，能够满足上述要求。（5）心墙保护层：心墙顶部应设保护层，防止冰冻和干裂。保护层可采用砂、砂砾或碎石，其厚度不小于该地区的冻深或干燥深度，此处取2.00m。（6）反滤层：本设计设置两层反滤层，坡度取1:0.2，上游宽度为4.0m，下游取6.0m。3.3渗流计算3.3.1渗流计算的基本假定渗流计算中，基本假定如下：（1）本设计中，心墙粘土料的渗透系数为K＝5×10-6cm/s，坝壳堆石体的渗透系数为K＝3×10-1cm/s，可以发现，心墙粘土料与堆石体的渗透系数相差了105倍，因为渗透系数相差很大，粘土心墙透水量相对堆石体非常小，故把心墙看作是相对不透水层，因此在计算中，可以对上游楔行降落水头的作用不予考虑。在大坝下游设有棱体排水，可近似地假定下游水位与堆石棱体内坡的交点就是浸润线的逸出点。（2）渗流在土体中的流动速度比较小而且渗流处于层流状态，且满足达西定律。（3）渗流为连续的，且当渗流发生时土体孔隙尺寸以及饱和度不发生变化[13]。3.3.2计算条件在进行渗流计算时，对计算条件的选取需要考虑如下几种情况，综合对比考虑后选择这几种情况中最不利的情况作为计算时的控制条件：（1）上游正常高水位，下游相应的最低水位；（2）上游设计或校核洪水水位，分别相应的下游水位；（3）对山游坝坡稳定最不利的库水降落后的水位。本设计取定计算条件为上游设计洪水位+相对应的下游水位根据工程的基本资料可以得到渗流计算的计算条件为：上游设计洪水位为2867.11m，下游水位取对应的水位，为2618.00m。3.3.3计算方法根据工程资料给出的大坝图纸，将如美心墙堆石坝最大断面进行简化，近似将坝基看作不透水地基，且坝下游无排水设备，故可用水力学方法进行渗流计算。其中有几种不同的方法，本次计算选用分段法进行计算[14]。分段法计算简图如图3.3图3.3两河口坝渗流计算简图分段法，顾名思义是将坝内渗流分为两段，然后根据达西定律和杜平假设，可以列出各段的运动方程，又因为水流具有连续性，故可以据此对渗透流速、渗透流量和浸润线等进行计算求解。为计算方便，可以对心墙断面进行简化，即将心墙简化为厚度为的bc等厚矩形断面。由于心墙的渗透系数K远小于坝壳的渗透系数K，故当渗流通过心墙时会产生巨大的水头损失，浸润线在心墙中发生很大跌落。因此，可以近似假定心墙上游浸润线与水库水位处于相同位置。当下游有堆石排水时，可以近似假定将下游水面与堆石内坡的交点处作为浸润线逸出点。3.3.4计算过程及计算结果分析假设坝址处地基为不透水地基。渗流的计算公式如下：q=qq=kh

式中：k——坝壳料渗透系数，取3×10-1cm/s；k1——防渗心墙渗透系数，取k=5×10-6cm/s；H1——上游设计水深，H1=287.11m；H2——下游相应水深，H2=38.0m；bc——防渗心墙平均厚度；bc=（b1+b0）/2=73m；L——防渗心墙右侧到溢出点的水平距离。L=D-bc/2，D为心墙重心到溢出点的水平距离:D=459.5m。故L=423m；h——心墙下游浸润线高度,m；T——透水地基厚度，取T=0m；联立计算可得h=38.103m、q=2.736×10-5m3/(s·m)大坝沿坝轴线方向的长度可以根据地形地质图测量得到，为L=668m，沿整个坝段的总渗流量Q=Lq,式中m是折减系数，因为坝宽、坝厚、渗流量沿坝轴线并非是均匀的，所以增加了该折减系数，取m=0.8，得Q=1.462×10-2m3/s。则浸润线方程为：y=h3.4坝坡稳定性计算3.4.1滑动面的选择堆石坝下游坝坡的堆石体内摩擦角=49.2○，堆石粘聚力忽略不计取为0。坝坡抗滑稳定计算采用简化的Bishop法，首先需要做的就是确定最危险滑动圆弧位置。本设计只取一个滑动面做稳定计算（滑动面的起点在坝顶、与心墙相交、与坝基接近或切入坝基、终点在坝坡脚附近），取圆弧半径为R=800m，然后对滑动面以上土体进行分条并编号，大坝条分如图3.4所示。图3.4大坝条分示意图3.4.2计算条件本设计选取为上游设计洪水位+下游相应水位作为坝坡稳定计算的计算工况。即取上游水位高程2867.11m，下游相应水位为2618m。3.4.3计算过程及结果分析简化的Bishop法计算公式如下：K=W在本设计中不考虑坝区地震，故可以将Q、V、Mc忽略不考虑，用容重代替法考虑孔隙水压力的影响，则将计算公式（3.3）简化为：K=W'tan计算稳定安全系数K时，第一步需要预先估计一个K值代入公式（3.4）的等号右部的方程计算出K值，将新计算出的K值继续代入（3.4）重新计算出新的K值，如此重复计算下去，直到代入的K值与新计算得的K值相等或者两者相差不大为止，此时所得的K值即代表计算滑动圆弧的稳定安全系数值。滑弧面的稳定计算结果见表3.1。表3.1坝坡抗滑稳定性系数计算表土条编号S(m2)W(KN)α土条长度l(m)W'*tanφ'+c'b

(KN)W*sinα(KN)1.00849.9018077.3749.0075.6020942.6414830.732.001760.2037439.4543.0068.8043373.6127986.353.002559.4054438.4439.0064.0063066.9337730.664.003011.3064050.3534.0060.4074202.3339653.265.003249.5069116.8730.0057.7080071.8938399.356.003037.8064614.0126.0055.8074855.3331571.847.002765.4058820.0622.0053.9068143.0424625.648.002220.4047227.9118.0052.6054713.5316346.4290.201211.3918.001403.39419.289.001042.7022178.2314.0051.6025693.486020.08704.309458.7514.0010957.962567.4910.00110.602352.4611.0050.902725.33504.161195.2016051.5411.0018595.703440.0011.00900.9012099.096.0076.5014016.791422.10总和245517.38预估K=1，代入公式（3.4）进行第一次计算，迭代过程见表3.2。表3.2坝坡抗滑稳定系数迭代表土条编号ABABABABAB10.6213036.240.7916493.800.8417686.880.8618027.740.8718120.2120.6327153.610.7733609.980.8235773.290.8436385.450.8436551.0830.6339913.550.7748678.420.8251557.940.8352367.660.8352586.3540.6547988.310.7757426.510.8160450.290.8361293.920.8361521.2650.6653039.890.7862478.150.8265439.770.8366260.790.8366481.6560.6851100.580.7959207.480.8261696.850.8362382.450.8462566.5570.7148253.040.8154938.720.8456945.400.8457494.330.8557641.4680.7440467.580.8345219.430.8546611.460.8646989.560.8647090.710.741037.980.831159.870.851195.570.861205.270.861207.8790.7819998.660.8521897.790.8722439.760.8822585.870.8822624.880.788529.190.859339.150.879570.290.889632.600.889649.24100.812215.880.882386.590.892434.290.902447.070.902450.480.8115119.630.8816284.410.8916609.850.9016697.070.9016720.32110.8812396.680.9212949.240.9313097.720.9413137.090.9413147.56总和380250442069461509466906468359K1.551.801.881.901.91其中，A1/cossintan`/K，B=A×(W'tanφ'+c'b)。综上，最终计算得到的抗滑稳定性系数为K=1.91。大于规范要求的1.5，故可以初步判定，正常水工况下，大坝坝坡的稳定性可以满足要求。

4渗流计算4.1有限元模型建立本次设计中，对两河口心墙堆石坝整个坝体进行了完全的模拟，根据坝体最大断面，在CAD中对大坝剖面进行适当的简化，在ANSYS中拉伸建出整个坝体的模型并再对坝体进行适当的简化。本设计使用ANSYS软件对堆石坝进行了三维建模，得到的三维有限元模型一共有148266个单元，146961个节点，该有限元模型单元的选择为，主要采用的单元为8结点6面体，部分采用了棱柱体单元，采用棱柱体单元主要是为了适应边界的过渡。在堆石坝的应力变形仿真计算中，本设计大坝堆石体材料的本构关系采用的模型是邓肯E-B模型。通过ANSYS软件建立的两河口心墙堆石坝的三维有限元模型如图4.1所示，心墙体三维有限元模型如图4.2所示，三维坝体剖面的填筑材料分区如图4.3所示。图4.1两河口心墙堆石坝三维有限元模型图4.2心墙体三维有限元模型图4.3三维坝体剖面材料分区4.2渗流计算原理假设在非饱和土中，达西定律仍然是存在的，并且当时间变化时，孔隙气压力不会因此而产生变化，并且对于土结构平衡条件和不同流体流动之间的相互作用可以不用考虑，综上，就可以得到当地下水非稳定渗流流经饱和区与非饱和区时的控制方程[14]：∂∂xkx稳态渗流即水头H不会随着时间变化而变化时，则根据（4.1）能够得到地下水稳态渗流的控制方程：∂∂xk式中，H为水头，y为位置水头；kx、ky为x、y方向的渗透系数；Q为微元体边界流量；mw为体积含水量变化系数；4.3计算参数本设计中渗流计算的计算参数如表5.1所示表4.1材料渗透系数材料名称渗透系数K（cm/s）心墙料5×10-6反滤层15×10-3反滤层25×10-2过渡层3×10-1堆石区3×10-14.4计算工况水体在土体中流动时会产生水头损失，主要是因为水体在土体中流动时，土颗粒会对水体有阻碍作用，对水体产生了阻力，而由牛顿第三定律即作用力与反作用力的原理可以知道，相应于水体受到的来自土颗粒的阻力，水体流经的地方也会对土体颗粒施加一种渗透体积力，这种渗透体积力会对坝体的稳定和应力变形起很大的作用，为了真实地模拟坝体的应力变形状态，考虑水库蓄水全过程，为后期进行坝体流固耦合分析计算提供基础数据，需要选取多个蓄水高程，分别对这些蓄水高程进行稳定渗流分析。在上游水位中选取13个水位高程作为本设计有限元稳定渗流计算工况，具体见表4.2。下游水位取2595m。表4.2上下游水位对应表序号上游水位下游水位12605259522678259532685259542693259552742259562785（死水位）2595727902595828082595928202595102832259511284425951228532595132865（正常蓄水位）2595针对蓄水过程，列出2678m、2785m、2865m三个水位的计算结果。4.5计算结果经三维有限元法计算，选取有限元网格z方向处于最中间的一层网格作为最大断面进行分析整理，各水位下心墙削减水头百分率如表4.3所示，绘制剖面等压力水头线如图所示4.1~4.3。表4.3各库水位高程心墙削减水头值库水位高程心墙内浸润面高程削减水头百分率上游下游差值26782674.32595.878.596.0%27852782.42598.6183.896.1%28652860.62600.5260.196.3%注：表中削减水头百分率=（H上心墙−H由坝体剖面压力水头线分布图可见，浸润面在心墙上下游发生了突降。由表4.3可见，心墙削减水头随着库水位的增高而变大，库水位高程2678m、2785m、2865m下，心墙削减水头分别为78.5m、183.8m、260.1m，分别占上下游水头差的96.0%、96.1%、96.3%。由此可见，心墙的防渗效果十分显著。图4.4上游水位2678m高程渗流场压力水头图4.5上游水位2785m高程渗流场压力水头图4.6上游水位2865m高程渗流场压力水头

5三维有限元应力分析5.1邓肯E-B本构模型及参数切线弹性模量Et和卸荷下回弹性模量Eur计算方法为[15]：EtEur式中，Rf为材料参数；S表示应力水平；σ3为侧限压力；Pa为单位大气压力；Kur为卸荷模量；nur为卸荷模量指数。切线体积模量Bt和粗粒料的内摩擦角φ的计算方法为：Btφ=式中：Kb为无因次的体积模量数；m为模量指数；φ0为σ3等于单位大气压力时的内摩擦角；∆φ为φ值σ3随本次不考虑流变和湿化的静力计算所采用的邓肯E-B模型材料参数具体数值如表5.1所示。表5.1两河口心墙坝E-B材料参数材料分区KnKbmRfC（kPa）φ（度）φ0（度）△φ（度）KurNur堆石Ⅰ区9200.292900.20.7518338.549.25.018400.29堆石Ⅱ区8530.252800.160.7526035.349.28.117060.25堆石Ⅲ区7830.272870.190.7519737.446.85.315660.274心墙4500.413500.390.8812223.0\\9000.41反滤层18300.282950.230.7720035.747.57.016600.28反滤层29150.253160.210.7825738.248.05.318300.25过渡层9040.282480.20.7424637.350.88.418080.28接触粘土1510.2475.90.260.864219.7\\3020.245.2流变本构模型及参数通过对流变试验的研究，选择用以指数型衰减的Merchant模型来模拟常应力下的ε~t衰减曲线。Merchant模型的蠕变曲线表达式[16]如下：εt=εi式中：εi=σ对式（5.5）求导，可以得到：ε=αεfe-αt（由上式可见αεf为t=0时的初始变形率，而α则为初始相对变形率。故可以ε=αεf在Prandtl-Reuss的假设下，可以将应变率的张量表示如下式：ε=13式中s为偏应力，σs沈珠江对体积和剪切变形假定的变形速率如下：εV=αεγ=αγf式中εVf和γf表示εVf=bσγf=dSl1−式中εVt和γtεvt=εγt=γΔt广义剪应变和广义体积应变可表达为：γ=29ε在轴对称的情况下，即ε2=ε3γ=ε1−在使用有限元方法做应力应变的分析中，常把流变作为初应变来进行考虑[17]。为了能够让因流变变形而引起的应变增量在对坝体应力应变的分析过程中被考虑到，需要把堆石体流变变形的体积流变εVf和剪切流变γf的计算公式进行修正，εvf=bσ3式中，b、c、d、m1、m2和m3为模型参数，S本次心墙堆石坝进行流变计算所采用的材料参数如表5.2。表5.2两河口心墙坝材料流变模型试验参数材料分区αbcdmmm堆石Ⅰ区0.000720.019430.005560.029580.2010.2810.552堆石Ⅱ区0.000720.019430.005560.029580.2010.2810.552堆石Ⅲ区0.000720.019430.005560.029580.2010.2810.552心墙0.000720.019430.005560.029580.2010.2810.552反滤层10.000720.019430.005560.029580.2010.2810.552反滤层20.000720.019430.005560.029580.2010.2810.552过渡层0.000720.019430.005560.029580.2010.2810.552接触粘土0.000720.019430.005560.029580.2010.2810.5525.3大坝填筑蓄水过程模拟本设计大坝主体施工共分6期，坝体在填筑到6期中时开始蓄水，然后边填筑边蓄水，开始蓄水后第三年为运行期，具体有限元荷载步见表5.3（大坝填筑过程共占58步）。表5.3大坝有限元模拟荷载步工程项目高程（m）有限元荷载步坝体填筑28431~48坝体填筑和蓄水287849~58蓄水287859~67坝体运行287868~77注：静力计算只考虑坝体填筑和蓄水，不用考虑坝体运行，而流变计算需要考虑坝体运行期。5.4静力计算两河口心墙堆石坝坝体的静力计算结果如表5.4。表5.4两河口心墙堆石坝应力变形最大值坝体计算工期竣工期蓄水期坝体变形（cm）向上游水平位移-40.6-31.8向下游水平位移101.4130.2竖向位移-353.6-369.1坝体应力（MPa）小主应力5.045.18大主应力2.362.48注：竣工期是指坝体填筑完毕的时刻，在有限元模拟荷载步中为58步；蓄水期是指大坝蓄水至正常高水位的时刻，在有限元模拟荷载步中为67步。规定向下游方向和竖直向上为正方向。5.4.1坝体位移两河口心墙堆石坝的坝体最大断面位移图如图5.1所示（取有限元网格z方向处于最中间的一层网格为最大断面，即z=135m处的断面）。由计算结果知，心墙体的土料相对来说比较软，所以在心墙中部部位会发生最大沉降，由计算结果图可以看到，最大沉降在竣工期和蓄水期时产生的位置是基本一致的，而相较于大坝竣工期，在蓄水期，坝体最大沉降值会增大一些。在蓄水期，因为水压力的作用，坝体向上游位移减小，向下游位移增大，坝体整体向下游变形较为明显。大坝竣工期和蓄水期的铅直位移规律也相似，最大铅直向下位移都出现在心墙中部位置。大坝上游坝体表面局部位置出现了铅直向上的位移，该现象的产生是由于竣工期水位较高，浮托力使坝体产生上抬，但是这种现象发生的范围很小，对坝体的影响几乎没有。在蓄水期时，出现这种现象的范围随着浮托力的增大而变大。（a）坝体竣工期顺河向水平位移图（b）坝体竣工期垂直沉降图（c）坝体蓄水期顺河向水平位移图（d）坝体蓄水期垂直沉降图图5.1坝体位移图（单位：cm）5.4.2坝体应力两河口心墙堆石坝的坝体最大断面应力图如图5.2所示。由图可以发现，心墙的变形比两侧堆石变形是要大一些，这是因为心墙的土料要软一些，这也会导致心墙与堆石之间产生了明显的变形不协调现象，引起心墙和堆石体之间应力的重分布，从而出现了明显的拱效应[16]。由图可知，竣工期、蓄水期的坝体应力规律相似。这是由于在竣工期和蓄水期，两河口心墙堆石坝被施加的荷载主要是自重荷载、渗透体积力和浮托力，其中大坝的自重荷载是大坝的主要荷载，且采用边填筑边蓄水的方案，竣工期的水位较高，相比正常蓄水位，蓄水期水位上升较小，因此荷载在竣工期和蓄水期的差别不大，故由于这两个时期的荷载的变化很小，在这两个时期坝体的大主应力虽然随着水位的升高有所增大，但是变化比较小。由于库水推力的作用，使大坝坝体向下游变形，同时在心墙底部也存在一定弯曲，故在心墙底部下游侧的应力会增加，这个部位就是会出现应力的最大值的位置。（a）竣工期大主应力图（b）竣工期小主应力图（c）蓄水期大主应力图（d）蓄水期小主应力图图5.2坝体应力图（单位：MPa）5.4.3心墙位移坝轴线心墙剖面位移图如图5.3。由图可知，由于心墙土料较软，坝体的最大沉降发生在心墙内，在竣工期，心墙的最大沉降发生部位为心墙内高程约2740~2750m处，心墙向两岸坝肩沿坝轴线方向发生变形，变形基本呈对称分布，向两岸最大位移为56.0cm。在蓄水期，心墙蓄水期和竣工期的位移规律相似，蓄水期心墙向两岸最大位移较竣工期有所增加，心墙向两岸最大位移大致呈对称分布。由于堆石体和心墙之间存在黏滞摩擦作用，使蓄水期心墙的沉降值减小。（a）心墙竣工期横河向位移图（b）心墙竣工期铅直位移图（c）心墙蓄水期横河向位移图（d）心墙蓄水期铅直位移图图5.3心墙剖面位移图（单位：cm）5.4.4心墙应力坝轴线心墙剖面应力图如图5.4。由图可知，在竣工期，心墙的大小主应力的最大值出现的位置，都在心墙的底部，等值线基本呈水平状；黏土与岸坡混凝土之间存在摩擦作用，使得在靠近岸坡处出现局部弯折；同时岸坡具有顶托作用，使得靠近岸坡部位的小主应力比同高程心墙内部的应力小。心墙在蓄水期和竣工期这两个时期的应力分布规律没有太大的差别，是基本一致的。心墙在蓄水期时，大主应力最值为5.0MPa；小主应力有所减小，最值为2.56MPa。心墙在竣工期、蓄水期均没有出现拉应力，故在心墙内部产生裂缝的可能性较小。（a）心墙竣工期大主应力图（b）心墙竣工期小主应力图（c）心墙蓄水期大主应力图（d）心墙蓄水期小主应力图图5.4坝轴线心墙剖面应力图（单位：MPa）5.5流变计算考虑流变效应对大坝应力变形的影响，两河口坝的流变计算结果见表5.5。表5.5两河口心墙堆石坝应力变形最大值坝体计算工期竣工期蓄水期运行期坝体变形（cm）向上游水平位移-58.1-53.0-55.8向下游水平位移117.9159.8173.7竖向位移-491.5-519.7-558.6坝体应力（MPa）小主应力4.414.524.54大主应力2.12.052.00注：竣工期是指坝体填筑完毕的时刻，在有限元模拟荷载步中为58步；蓄水期是指大坝蓄水至正常高水位的时刻，在有限元模拟荷载步中为67步；运行期是指大坝蓄水至正常高水位一年后时刻，在有限元荷载步中为77步。5.5.1坝体位移考虑流变效应的两河口心墙堆石坝的坝体最大断面位移图如图5.5所示。由计算结果知，在考虑了流变效应后，在水平位移，大坝竣工期、蓄水期和运行期规律相似，向上游最大位移出现在上游堆石区1/3高程处，向下游最大位移出现在心墙下游面与反滤层交接的位置，约坝高1/3高程处在。蓄水期大坝向上游的位移减小，向下游的位移增大，坝体整体向下游变形较为明显。运行期大坝向上游位移与蓄水期相近，向下游位移则进一步增大。在铅直位移，大坝竣工期、蓄水期和运行期铅直位移的规律也相似，最大铅直向下位移都出现在大坝心墙中部位置。大坝的上游坝体表面局部位置出现了铅直向上的位移，该现象产生的原因与静力计算中的情况相同。但是这种现象发生的的范围很小，对坝体的影响几乎没有。在蓄水期时，出现这种现象的范围随着浮托力的增大而变大。由于考虑了流变效应，与静力计算结果相比，竣工期和蓄水期的最大沉降值的变化减小；而大坝坝体在竣工期和蓄水期这两个时期时，向上游和下游的位移都增加了。（a）坝体竣工期顺河向水平位移图（b）坝体竣工期垂直沉降图（c）坝体蓄水期顺河向水平位移图（d）坝体蓄水期垂直沉降图（e）大坝运行期顺河向水平位移图（f）大坝运行期垂直沉降图图5.5考虑流变效应坝体位移图（单位：cm）5.5.2坝体应力考虑流变效应的两河口心墙堆石坝的坝体最大断面应力图如图5.6所示竣工期、蓄水期坝体应力规律计算结果与静力计算的应力结果相同。从竣工期到蓄水期大主应力随着水位的升高有所增大，而竣工期和蓄水期荷载的变化较小，致使大主应力在竣工期、蓄水期变化较小。且运行期的水位不变，故大坝运行期的大主应力变化很小。因为流变效应，坝体的小主应力从竣工期到运行期均有所减小。（a）竣工期大主应力图（b）竣工期小主应力图（c）蓄水期大主应力图（d）蓄水期小主应力图（e）运行期大主应力图（f）运行期小主应力图图5.6考虑流变效应的坝体应力图（单位：MPa）5.4.3心墙位移考虑流变效应的坝轴线心墙剖面位移图如图5.7。由图可知，竣工期，蓄水期和运行期的心墙位移规律相似，心墙向两岸坝肩沿坝轴线方向发生变形，变形大致呈对称分布，心墙向两岸的最大位移为146.2cm，蓄水期较竣工期增大但差别不大，由于流变作用，运行期较蓄水期增大。由于堆石体和心墙之间存在黏滞摩擦作用，心墙沉降值在蓄水期较竣工期有所减小，而由于流变效应，使得运行期心墙最大沉降值增大。（a）心墙竣工期横河向位移图（b）心墙竣工期铅直位移图（c）心墙蓄水期横河向位移图（d）心墙蓄水期铅直位移图（e）运行期横河向位移图（f）运行期铅直位移图图5.7考虑流变效应的心墙剖面位移图（单位：cm）5.4.4心墙应力考虑流变效应的坝轴线心墙剖面应力图如图5.8。由图可知，考虑流变效应后，心墙大小主应力的最大值出现的位置仍然是心墙的底部，等值线也还是基本呈水平状，且仍然由于黏土与岸坡混凝土之间的摩擦作用，会在靠近岸坡处出现局部弯折；岸坡的顶托作用，还是会使靠近岸坡部位的小主应力比同高程心墙内部的应力小。心墙在竣工、蓄水和运行这三个时期的应力分布规律相差不大，基本是一致的。在蓄水期，心墙的大主应力最值为4.51MPa；小主应力有所减小，最值为2.39MPa。运行期心墙大主应力最值为4.52MPa；小主应力有所减小，最值为2.34MPa。竣工期、蓄水期和运行期均没有出现拉应力，故在心墙内部产生裂缝的可能性较小。（a）心墙竣工期大主应力图（b）心墙竣工期小主应力图（c）心墙蓄水期大主应力图（d）心墙蓄水期小主应力图（e）心墙运行期大主应力图（f）心墙运行期小主应力图图5.8考虑流变效应的心墙剖面应力图（单位：MPa）5.6小结综上在静力和流变计算条件下不同时期的大坝应力变形计算结果汇总如表5.6~5.7所示。表5.6竣工期、蓄水期和运行期应力变形最大值坝体计算工况竣工期蓄水期运行期静力流变静力流变流变坝体变形（cm）向上游水平位移-40.6-58.1-31.8-53.0-55.8向下游水平位移101.4117.9130.2159.8173.7竖向位移-353.6-491.5-369.1-519.7-558.6坝体应力（MPa）小主应力5.044.415.184.524.54大主应力2.362.12.482.052表5.7蓄水完成时如美坝体最大沉降与已建工程实测资料结果对比观测数据小浪底瀑布沟糯扎渡两河口最大坝高（m）154186261.5295沉降最大值（m）3.71.963.555.586占最大坝高百分比（%）2.401.061.361.89由以上表格可知：（1）流变效应会让大坝大主应力和小主应力都有所减小；（2）流变效应都会让大坝向上游和下游的位移以及铅直位移都有所增加；（3）大坝坝体的大主应力因为竣工期和蓄水期时的荷载变化较小，产生的变化也并不大。而又因为运行期的水位不变，致使大坝大主应力在运行期的变化也还是很小。（4）与已建工程实测资料对比可知，两河口心墙堆石坝运行期最大沉降占最大坝高百分比（1.89%）介于小浪底心墙坝（2.40%）与糯扎渡心墙坝（1.36%）之间，计算结果是合理的。6总结与展望6.1总结本设计首先在常规设计中利用分段法对两河口心墙堆石坝进行渗流计算，确定大坝心墙下游堆石体内浸润线所在的位置，并利用瑞典圆弧法对两河口心墙堆石坝的坝坡抗滑稳定性进行了计算，然后在专题计算中，利用ANSYS和ABAQUS软件，对两河口心墙堆石坝进行了三维有限元计算分析，对不同工况下大坝坝体和心墙的应力变形规律进行了计算分析。最终主要结论如下：（1）本文常规设计中对两河口心墙堆石坝的坝坡抗滑稳定的计算采用了简化的Bishop法，最终计算得到的结果，是符合一般规律的，坝坡抗滑稳定性满足要求。（2）本文专题计算中通过建立的两河口心墙堆石坝的三维有限元模型，模型中考虑了坝体各种材料分区以及心墙靠近两岸坝基的渐变形状等细节；在计算过程中真实地模拟了坝体的填筑过程和蓄水过程。（3）两河口心墙堆石坝在竣工期和蓄水期和运行期等不同阶段，随着时间推移，坝体最大沉降量不断增加；大坝向上游的位移随着时间推移而增大，但趋势减缓；随着水位上升，在水压力带来的作用之下，相比竣工期，坝体在蓄水期向上游位移逐渐开始减小。大坝向下游位移明显增大，其整体向下游变形明显。向上游最大位移出现在上游堆石区1/3高程处，向下游最大位移出现