洮河海甸峡堆石面板坝在稳定状态和非稳定状态下渗流分

洮河海甸峡堆石面板坝在稳定状态和非稳定状态下渗流分 1绪 论 1.1本文研究得意义 大坝作为一种挡水构筑物一般都会涉及到渗流稳定性问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ,并且它们的破坏主要都是从渗流开始。如若不及时引起重视进行处理,必将危及大坝的整体稳定性。坝体发生渗流破坏的原因,通常都是由于水力梯度引起的。因此，要研究确定水工构筑物的水力梯度的分布,同时必须通过其进行渗流分析,确定浸润自由面，再对坝体进行渗流稳定分析，以此对大坝的安全运行进行作出合理评估。 水利枢纽工程在蓄水期，由于库区内水位的上升，在水压力的作用下会在大坝的坝体内部和大坝基底产生渗流的现象，这对水工建筑物的正常运行会产生比较不利的影响，可是又是不可避免的。根据统计，我国已建成上万座各种类型水库，经过长年累月的运行，很多坝体在不同程度上都存在严重的病害或者隐患，这种病险危害的存在会从根本上影响水工建筑物的渗流稳定性及其变形稳定。近年来，在水利枢纽工程的正常运行过程中发生了很多具有灾难性的的破坏，经过广泛分析研究后，分析得出在可能对坝体造成破坏的所有模式中，渗流都在其中扮演着重要的作用。因此，对大坝渗流的研究是大坝安全检测的一项重要内容，深入的分析和研究坝体渗流状况，这对于我们水工建筑物的安全运行以及应对具有破坏性事故的威胁和保护人民生命财产具有重大而深远的意义。 目前较为理想的一种计算渗流场的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 是有限单元法。在渗流场的基础上对其进行整体稳定性分析,才能得到较为可靠的结果。因此,用有限单元法计算其渗流场，对挡水构筑物安全、维护、治理等都有巨大的社会效益和经济效益。 有限单元法的基本思路和原理早就有人提出，它是一套求解微分方程的数值计算法，与传统的解题法相比，有物理意义直观明确，解题效能高方法本身的适应性强，解题步骤易规范化，电子计算机配合使用后，有限单元法得到了迅速的应用和发展。 1.2国内外渗流计算研究的发展现状 法国工程师达西(H.Dracy)，在1986年通过试验研究提出了线性渗透定律，为渗流理论的发展莫定了基础。达西的试验是通过液体在均质沙土中作匀速流动的情况下进行的。这个研究结果已被后来的学者逐渐推广至整个渗流计算领域中 去，达西定律成为渗流分析计算中最基本、最重要的公式,它解决了土体的基本渗透计算理论。达西定律研究得出的公式是对均匀渗流平面上计算平均流速以及渗流区域内任意一点的渗流流速。1857年，法国学者杜比推导出杜比公式来研究非均匀渐变渗流分析的规律。1889年，H.E茹可夫斯基研究发现土体中渗流规律分析计算过程符合拉普拉斯方程，首次推导出了渗流的微分方程。但是，相对于许多数学模型和解析解的研究只能适用于均匀的渗透介质和边界比较简单的条件基础上，因此在实际工程的运用中就要受到一定限制。1914年，P. Forehheiness提出利用绘制流网的方法进行图解渗流试算。流网针对分析整个渗流场这个二维平面问题给出了比较可靠全面的图案。绘制出了流网，整个渗流场问题就会得到解决。直到1922年，为解决边界条件和渗流条件比较复杂的渗流问题，H.H?巴甫洛夫斯基研究出了一个有效的求解渗流场的工具—电拟法。所谓电拟法，就是采用电导液模型进行渗流模拟试验，但仍有一些弊端就是无法对非均质各向异性渗透介质的渗流进行模拟实验，当遇到一些比较复杂的地质和复杂的边界条件时就是去了适用性。后来为了对渗流有进一步的研究和提高，针对电拟法改进而形成了电网络法。 近年来，在计算机电子技术日新月异发展的基础上，很多渗流计算方法也随之出现。其中典型的数值模拟方法有边界元法、有限差分法、有限单元法三种，而且以围绕这三种方法为主的研究大坝坝体渗流分析计算在国内得到了迅速发展。 边界元法是在贝蒂互换定理经典力学理论及弗雷德霍姆积分方程基础上建立起来的，该方法提出于20世纪九十年代末、21世纪初。边界元法的基础理论形成的相对较晚大约在六十年代末。边界元法在最初被称为边界积分仪法。边界元法这个名词直到1978年才得到渗流研究领域的认可 有限差分法是1910年理查德森提出来的，该方法由于在渗流分析研究中的长期应用已具有一定的研究基础。目前，有限差分法的理论和和研究 经验 班主任工作经验交流宣传工作经验交流材料优秀班主任经验交流小学课改经验典型材料房地产总经理管理经验 已在渗流研究领域有了广泛的推广。1960年，克劳夫为了区别有限差分法，而最早提出并采用“有限单元法”这个名词。而在1965年，由于研究后发现有限单元法具有可以适用很多以变分进行的的渗流场计算的研究分析问题使该种方法在渗流领域的发展有向前迈进了一步，而且其成果也十分的显著有效。 随着水能源的大力开发，为了使水利工程更好的服务于人类，在相关坝体渗流的研究和应用在国内外对都引起了高度重视，特别是在土石坝的渗流研究领域，各国学者和专家一直在致力于对坝体渗流机理的探索。 达西定律作为渗流理论最基本的理论基础以及渗透系数的提出就为在水头作用下会产生渗流水量的函数关系提供了强大的应用基础。在对水工模型进行数值模拟的基础上利用有限单元法对坝体非饱和区域的渗流状况进行了实验。研究表明数值模拟能够对一些复杂的边界条件计算并具有较高的精度。随着渗流理论的逐步发展自适应理论在渗流分析计算中也得到了较为广泛的推广，并将在二维 这种自适应理论在无压渗流分非线性的渗流分析中吸收了线性误差估计的方法。 析方面也进行了验证分析。通过划分单元网格的方法，在渗流分析过程中让单元自己根据渗流场的变化去适应，这样就会大大提高其计算进度和效率。 随着现代计算机技术的高度发展和应用，有限元法也得到了迅速发展。近几年来，国内外众多水利研究专家和水利水电工程技术人员在从事水利枢纽工程研究和建设过程中，从渗流的基本理论出发，充分利用现代计算机技术的基础上结合实践都深入的就渗流方面做了探究。在解决渗流场的研究计算方面对有限元法的应用有了相对较大的进展，有限单元可以对一些简单的坝体渗流状况作出合理的分析。但是水工建筑物作为一种实际的工程结构以及所处环境条件十分复杂，要寻找一种简单、可靠、适用以及能对多种复杂边界条件的渗流问题以上研究方法还有待实践的验证。本文在研究堆石面板坝的坝体渗流曲线方面做了简单的尝试。 1.3有限元法在渗流分析中的发展现状和存在问题 早期从事渗流分析方面的研究学者，采用有限元方法分析具有自由面的无压渗流问题。首先对自由面的位置进行假定，然后进行试算过程，通过不断的的迭代计算一步步的逼近结果值，最后获得较为近似的结果。可是往往由于渗流流域的边界具有的复杂性，渗流区域的范围是未知的，这类渗流问题需要通过反复不断的分析，再反复不断的修改自由面的位置，直至求解、迭代逐步确定。其中变网格法和固定网格法是修改渗流自由面最常用到的分析方法。 在渗流计算过程中由于渗流自由面的位置是未知的，由于渗流边界是未定的这就使得求解边值问题存在一定的难度。这就使得渗流自由面和渗流溢出面更加 困扰渗流场的复杂求解。在求解不定边界问题的渗流区域中未定边界部分通过迭代逼近的方法得到合理的结果。 (1) 变网格法也是一种利用单元划分方法离散研究区域，在求解渗流问题 的过程逐渐迭代直至收敛的过程，在求解坝体渗流自由面与坝体渗流溢出点的方面其优点是迭代过程会随着求解发展自行稳定形成。 为了突破这个局限提出了新的研究方法。在变网格发的基础上，渗流科学研究者提出了固定网格法，其核心思想就是在渗流分析过程中划分的不变网格。 (2)固定网格在渗流分析的研究中是运用的较为成熟的理论方法之一。其中纽曼提出的Galerkin方法对固定网格在渗流分析中的作用意义深远。 固定网格法是将整个渗流区域纳入计算求解的范围。由于渗流的作用会在自由面附近形成饱和区域和非饱和区域，固定网格法会将这两个区域视为一个整体划分计算的方法。 1.4目前国内外渗流计算方法的发展现状 1.4.1常用的渗流计算方法 常用的渗流计算方法分为两大类;如图所示。 图1-1渗流计算方法 Figure 1 The methods of seep computing 一理论分析法。 (1)解析法中水力学法和流体力学法是理论分析法的两大类。 水力学法，它是一种比较简单的渗流分析方法，主要针对一些渗流条件比较简单的分析计算，由于计算分析过程简练，因此在一些坝体材料和结构较为单一的渗流问题普遍运用水力学法。众所周知，坝体的渗流计算需要求解许多的渗流要素，而水力学法不能满足这一要求只能求解坝体某一界面渗流要素的平均值，但当细化到每一个具体的点时就会出现制约，再加上水力学法求解的前提是对在渗流场假定的基础上的，可能会与现实的情况有一定的差距，而正是各种误差的累积导致我们得出的结果精度不高。为了解决水力学法中的不足，通常会用流体力学的方法加上一些实验取得的结果对采用水力学法所得值进行修正，这种方法在一定程度上会提高计算精度。这种方法的原理是基于流体力学在应对边界条件复杂基础上对反渗流分析研究的能够较大提高计算结果精度建立起来的。在渗流问题的计算中对任意一点的各个渗流要素和任一截面处的渗流要素如水头、压力、流速等等。由于要利用实验所得结果和流体力学法同时对水力学所得结果修正，这就会使问题求解变的比较复杂，而且其对渗流边界的研究还不能满足坝体渗流区域复杂的情况。 (2)数值方法。边界元法、有限单元法和差分法是目前最常用的三种数值方法。 边界元法:首先将渗流区域复杂边界通过单元划分的方法离散为能适应复杂边界的子域，然后通过区域积分形成积分方程，最后再化为线性求解方程进行计算。边界元法相对于有限元法而言具有以下优点:便于对半无限或无限渗透介质的处理、渗流奇异和自由面等问题。一般说其计算精度也高于有限单元法。边界元法的缺点:当面对三维非均质各向异性的渗透介质问题时，在应用方面尚存在相当大的困难。 有限单元法:它归纳总结了对古典近似计算方法，既吸收了变分计算中选择试探函数的内核，又继承了有限差分法中离散处理的精髓，并对分析区域进行积分，这一过程是对每个子单元对整个区域影响的直接表达。相对于边界元法该法在面对复杂的渗流边界时更好的模拟渗流曲线以及对非均质的坝体渗流方面都体现出了一定的优越性。 差分法:在渗流研究中差分法的基本思想是进行区域离散，离散的各个子域 会形成渗流的微分方程。再通过控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点水头值为未知量的方程组。将微分问题转化为函数运算问题。但是常常由于差分网格的规则受到一定的局限;其次，在对各向异性的曲线边界进行模拟分析时会有一定的困难。 (3)图解法:在一定理论分析和实验结果的基础上对所研究渗流区域进行流网绘制，利用中线法可求得断面处的渗流流速，也可以根据绘制的单元网格绘制坝体内的各个节点的渗流等水头线。 二试验分析方法。 由于大多数渗流基本理论是在试验基础上总结出的结果，所以采用试验分析的方法也是渗流计算研究得一种途径。缝隙水槽法、砂槽模型法、水力积分仪法、粘滞流模型法、电阻网模型法和电拟试验法。现只对以下三种实验分析方法 进行简单介绍。 (1)缝隙水槽模型试验:是指利用两块玻璃板或者有机玻璃板做成的水槽模型、然后将其平行地安装在试验水槽内，并在两块模型的板之间保持一定的缝隙宽度。试验的观测方法:一般是在玻璃板上顶先画上方格，然后将不同颜色的液体注入上游的不同高程。这样能使模型中的渗流线能明显地显示出来，便于观侧。通常在槽中注入的是纯甘油。 (2)砂槽试验:砂槽试验利用砂或者砂性土做模型，将模型的上下游水位保持与实际建筑物上下游水位相对应的时，在上下游水位水头差的作用下产生与实际相当的渗流，并且形成自由渗流的水面曲线。通过量测设备量测出水面线高度。 (3)电拟法:是在渗流和电场的控制方程符 合同 劳动合同范本免费下载装修合同范本免费下载租赁合同免费下载房屋买卖合同下载劳务合同范本下载 一种形式的基础上而进行求解的。电网络和导电液模型是目前电拟法所具有的两种模型。由于在模型中其材料是各向同性的均匀介质，只能对求解一些急变的渗流区域。 目前，在研究渗流问题分析应用中对坝体材料为非均质各向异性渗透介质的研究较为苦难，尤其在地质条件比较复杂的坝址条件下电网络法也不尽适用，电网络法在一定程度上优化了其它渗流研究的缺陷，其所得结果由于没有误差累积，又有一定的稳定性，一直以来也视作是一种应用性较强的求解渗流场复杂问题的计算工具。电网络法的特点:在求解渗流问题的过程中不产生累积误差以及具有容量、稳定性基本不受限制在渗流分析中沿用已久，在求解大型复杂渗流场 问题时任是有效的工具。 1.5主要研究内容与研究方法 本文在大坝渗流曲线的研究中有效结合有限元分析软件ANSYS中的生死单元技术，通过渗流场的实际条件决定单元生或死。这种方法对于跨自由面单元固定网格法有一定的弊端，这必将会影响其计算结果的准确性。在有限单元法的求解中，通过网格划分的各个单元会形成单元渗透矩阵，而所有子单元最终组合成总渗透矩阵，处于饱和区与非饱和区域的自由面单元通过调整可视为一种复合材料单元。由于各个单元彼此相互作用而节点与节点也相互联系，在形成的结点处渗流要素以多项插值连续函数表达出来。如果由于单元渗透矩阵的不稳定性会导致总渗透矩阵计算结果值精度降低，何况在应用更加复杂的三维渗流计算方面。 总之，从以上各个研究分析方法可以看出，对渗流自由面的确定需要进一步研究，提出更新更全面的更有效的分析方法。 为了对大坝的安全运行、社会效益及经济效益作出积极的 评价 LEC评价法下载LEC评价法下载评价量规免费下载学院评价表文档下载学院评价表文档下载 ，就要求我们对建成坝体的渗流状况进行掌握，制定出有效的工程措施控制坝体渗流。而这些都是建立在对坝体渗流量、流速大小、水头等充分计算研究的基础上。本文通过众多资料的学习，初步探讨性地提出利用ANSYS有限元法对大坝进行渗流分析。以有限单元法的基本原理，利用大型有限单元计算软件ANSYS为平台，对坝体的渗流场进行仿真模拟计算。将坝体作为一个连续体对渗流区域部分进行离散，通过有限单元法把每一个单元的微分方程按边界条件要求在变分原理的基础上，形成线性方程组把泛函求极值的问题通过辅助设备精确计算求解。 2渗流分析的基本理论和计算原理 2.1渗流的基本理论 2.1.1渗流的基本概念 水或其它具有水流运动特点的流体在多孔介质中的流动，统称渗流。 根据渗流所处渗流场在渗流自由面的位置不同将渗流分为饱和渗流与非饱和渗流两种。稳定渗流和非稳定渗流的区别在于渗流场的渗流要素的大小和方向是否会随着时间的推移发生变化。其基本渗流要素有:渗流的水头和流速及水力梯度等。 由于要受到地下水的补给、径流和排泄等天然和人为因素的影响，引起的渗流要素中的中的任一参数随着时间而引起变化的，称此为非稳定渗流。处于河床中的坝体由于会受到库区水位的影响会产生变化，这些渗流基本上都是非稳定渗流，若果库区水位变化不大，引起的坝体渗流较为稳定时，可将非稳定渗流视为稳定渗流来考虑，对渗流计算量会大幅减小。 2.1.2多孔介质的渗透系数 多孔介质是指一个划分成很多微小体积的介质，这种介质中都存在微小的固体体积和液体体积。液体和固体构成体积的空隙和骨架，在有效空隙中，流体可以从一点连续运动到另外任意一点。 渗透系数是体现坝体渗流重要特征的一个参数，也称为水力传导系数。它是对多孔介质渗透强弱程度的一种反映。渗透系数大小与多孔介质的类型、孔隙度、空隙比以及其中渗流液体的密度和介质的粘滞系数都有一定的关系。 ,kgk,其表达式为: (2—1) , ρ——液体密度 γ——液体粘度; 2——9.8m/s; g ——渗透系数。(是指介质材料本身具有的属性与介质的大小、形状、排k 列相关)。 渗透系数决定于多孔介质的结构和流体的性质两种因素。不同的土类都有不同的渗透系数。 表2-1The permeability coefficient of several porous media 各种多孔介质的渗透系数 土类 渗透系数 -1-2 砂质砾 3×10,5×10 -1-2砾质砂 1×10,2×10 -2-2中砂 4×10,1×10 -2-2粉砂 2×10,5×10 -3-3砂壤土 5×10,1×10 -4-6粘壤土 5×10,1×10 -6 ?10粉质粘土 2.2渗流基本微分方程(控制方程) 2.2.1达西定律 考虑不可压缩流体通过均质各向异性的孔隙介质的稳定渗流，根据达西定 律，x、y、z方向的渗流速度可表示为: ,,V,,K (2—2) xx,x ,,V,,K yy,y ,,V,,K zz,z p,,y， (2—3) ,式中:—x方向上的流速; Vx V—y方向上的流速; y —z方向上的流速; Vx K、K、K 分别为沿x、y、z三个方向的渗透系数; xyz ,—水头函数; —流体压力; p —流体的容重; , 2.2.2渗流连续性方程 根据质量守恒原理见图2-1，推得渗流连续方程为: z PVxdy dx ,()PVx,PVdxx,xdz yx 图2-1流量守恒原理 Figure 2-1 The principle of flow conservation ,V,V,Vyxz，，,0 (2-4) ,x,y,z 将式(2—1)流速分量Vx、Vy、Vz代入式(3—3)得下式: ,,,,,,,,,(K)，(K)，(K),0 (2-5) xyz,x,x,y,y,z,z 当时，即孔隙介质为均质各向同性时，式(2—5)变为: K,K,K,Kxyz 222,,,,,, (2-6) ，，,0222,x,y,z 上式即为稳定渗流的拉普拉斯方程，即稳定渗流的基本微分方程。 在具体的坝体稳定渗流分析计算中会涉及边界条件的问题。这两类边界一种边界已知水头值，如处于上下游正常水位以下的坝体外部、坝基不透水层;另一类边界是渗流量大小能通过计算得出的。 ,,y已知水头边界: ,,已知流量边界: K,,q,n 结合土坝稳定渗流情况说明边界条件如图(2-1)所示，在上游水位以下坝面边界为已知水头值的边界，同样下游坝面下游水位以下也是水头已知的边界。 ,,而经过处理的坝基由于渗流量为0 作为一种特例出现，即，浸润面(自K,0,n 由面)也是边界，它应同时满足下式: ,,,,y ， (2-7) K,0,n 渗出点以下，下游水位以上的坝坡称为渗出段边界，它应满足下列条件:,,y。 ,,y,,,,, K,0,,n, 渗出点c 浸润面,自由面,4φ=y(p=0)1BjyH5iφ=Hm,,h2B,0c,h3x 坝基不透水层 图2-2 坝体渗流的边界条件示意图 Figure2-2 The diagram boundary condition of seepage dam 2.2.3二维渗流的基本微分方程 22,,,, (2-8) ，,022,x,y 将上式(7—6)当作欧拉方程，运用数学手段可得出相应的泛函积分表达式: ,,1,,22,，I[(x,y)][K()K()]dxdy (2-9) ,xy,,,,2xy, 在全部边界上。 ,,,B 有限单元法是吸取里慈法和差分法的思路和格式而发展起来的一种有效的数值解法。 在泛函积分表达式(2—9)中的积分区上，用三角形单元或其它形状的单元将Ω域划分成有限个子域(单元)?Ω。有限元法区域的离散要比差分网格的划分灵活得多，适应复杂形状的边界能力更强。 3有限单元法在渗流中的应用 3.1有限单元法概述 有限单元法的基本思路和原理早就有人提出，它是一套求解微分方程的数值计算方法，与传统的解题法相比，有物理意义直观明确，解题效能高、方法本身的适应性强，解题步骤易规范化，近年来，有限单元法随着电子计算机的广泛应用和数值分析方法发展得到了迅速的应用和发展。有限单元法之所以受到人们的青睐是因为在渗流区域研究中能对稳定渗流和非稳定渗流问题得以有效的解决。到目前为止有限单元法已可解决各种复杂边界的平面渗流问题及一些空间渗流 问题。随着计算机技术在全球的高速发展，其相关应用也得以长足发展，有限单元法就是在这种优势下快速发展起来的。 有限单元法以分块插值和剖分离散为指导思想，将渗流区域离散划分，对求解单元通过单元共用结点组合起来，单元参数之间形成一定的函数关系，在进行分段插值、分区求解，形成整个渗流场的场函数，再将未知量通过近似函数插值计算。为使解的精度不断提高，随之应增加单元自由度和单元数目这有利于结果值的收敛也有助于插值函数精度的提高。 在数值计算的方法中有限单元法与有限差分法相比在求解实施的步骤有些类似，但其着手研究的思路去不同，有限单元法是在变分原理的基础上求解泛函的积分表达式，通过微分方程将渗流边界条件转化成对泛函求极值的问题，而差分法的原理首先对微分方程的解析，以迭代逼近的方法获得离散方程的导数。 3.2有限元法计算原理及步骤 3.2.1定界区域的离散化 定解区域在渗流计算中就是指所研究的渗流场及边界，根据渗流场的形状与计算要求的精度来选定单元的形状。 3.2.2单元分析 单元分析有两方面的工作;第一是选择插值多项式代替原待定的水头函数φ(x,y)进行分析。第二是建立单元的渗透矩阵，这是有限元法的核心。 在平面渗流计算时，常用单元型式有三角形单元和四边形单元，本研究采用三角形单元。在三角形单元中多采用如下的线性多项式来近似地表达待求的水头函数φ: (3-1) ,(x,y),,，,x，,y123 式中:为待求的参数。 ,、,、,123 3-1三角形单元节点及坐标 Figure 3-1 A triangular element node and coordinate 根据上式(3-1)及图3-1所示，任意一个单元的节点水头值可表示为,,,x,y,,，，i12i3i ,,,,x,y,，，j12j3j, (3-2) ,,,,x,y,，，m12m3m, 代回式得单元水头函数的表达式: 将,、,、,,(x,y),,，,x，,y123123 1,,(x,y),(xy,xy)，,,jmmjixy，xy，xy,xy,xy,xyjmijmijiimmj ,,,, (3-3) (xy,xy)，(xy,xy)，(y,y)，(y,y),,miimjijjimjmimij ,,,,，(y,y,)x，(x,x,)，(x,x),，(x,x),ijmmjiimjjim 因此，可以将单元水头函数表达式写为: 1,，，, (3-4) ,(x,y),a，bx，cy,，(a，bx，cy),，(a，bx，cy),i1iijjjji1im2A A—为单元三角形面积。为使A为正值，单元节点次序必须顺着逆时针方向 变换。 简化水头函数的表达式，将以 a，bx，cy,iiiN,i,2A,a，bx，cy,jjjN, (3-5) ,j2A,a，bx，cymmm,N,m,2A, N、N、N表示,(x,y) 用得如下矩阵: ijm ,,,i,,e (3-6) ,,,,,,,,,(x,y),N,，N,，N,,NNN,,N,,,iijjmmijmj ,,,m,, 单元边界上水头线性变化，两个相邻单元公共节点上具有相同的水头值，因 此它们的公共边界上，应具有相同的水头分布值，即所选择的水头函数保证了相 邻单元之间水头分布的连续性。 3.2.3建立单元的特性方程 设整个定解域Ω内的泛函是，则任一单元可作为求解区域Ω的一I,,,(X,Y) ?ΔΩ 个子域ΔΩ。 Ω, 则在这个子域上的泛函表达式为 ，，,,1,,e22 (3-7) ,，I,,(x,y)K()K()dxdy,xy,,,,,,2xy，，,, 对单元节点值求导 ,,，，,,,,e,,,,,,,,Iy,x,, ,， (3-8) KKdxdyxy,,,,,,,,,,xy,,iii,,,,，， ,,对单元, 求导得 ,(x,y),N,，N,，N,iijjmm ,N,,,N,Njim,，， (3-9) ,,,ijm,x,x,x,x ,,,,,,,,,,,,,,,,,y,N,N,x,,,,ii,,,,,,x,,yii ,,,,,,,,,,,,,,,,,N,N,y,xjj,,,, (3-10) ,,,,,,y,,yjj ,,,,,,,,,,,,,,,,,y,N,N,x,,,,mm,,,,,y,,,ymm e,I将(3-9)(3-10)以上各式代入得: ,,i e,I1ee (3-11) ,,,,,,,Kbbb,,,b，Kccc,,,cxijmiyijmi,4A,i 同理计算简化可得: e,I1ee,,,Kbbbb，Kcccc,,,,,,,,,,xijmjyijmj,,4Aj (3-12) e,I1ee,,,,,,,,,,,Kbbb,b，Kccc,cxijmmyijmm,4A,m 将以上各式写为矩阵形式: e,,,I,,ee (3-13) ,,,,,K,,,,,,,,, 其中: ,,,i,,,,bbbbbbcccccciiijimiiijim,,eKK,,,,yex,,,,,,,jK,bbbbbb，cccccc,, ,,,,jijjjmjijjjm4A4A,,,,,,,mbbbbbbcccccc,,mimjmmmimjmm,,,, 上式即为单元特性方程。 e系数矩阵[K]只与单元的节点坐标及渗透系数有关，故又称单元渗透矩阵。 3.2.4建立定解区的总体方程 I,泛函,0的极值条件为是对整个的求解区域而言的。 ,,I,(x,y),,l e,,,I,,e式仅是某个单元域上的表达式，对某个节点l来说可能为多,,,,,K,,,,,,,j,, 个单元所共有， 值也必然分配在多个单元泛函之中，所以应当将有关的单元,l e,I特性方程中的叠加以后，才能令其为零，即: ,,l e,,II (l=1、2、3……n) (3-14) ,,0,,,,,ell 上式即为单元特性方程的结合式，式中n为节点总数，?表示对l节点有关的全部单元求和。对于求解区域内各个单元的未知节点的特性方程过程基本类似，但是由于在渗流场中不能对任一单元的节点无满足上式的条件，同样就,K 不存在变分。在整个区域中未知节点的总数n同形成的线性方程的数目是一致的,通过所有的未知结点的线性方程以及已知节点值建立起有自由项存在的数目为n的总方程组，一般以矩阵的形式表达，其矩阵又称总渗透矩阵: (3-15) ,,,,,,K,,F,0 ,F,为形成总体矩阵时，由已知结点构成常数项的列向量。 , ,,为未知结点的列向量。 [K]为n×n的方阵，称为总渗透矩阵。 3.2.5总体方程的数值求解 方程组的系数矩阵一般是正定的，这就保证了解答的唯一性。总体方程解出后，就可得到各未知节点的水头值，求得坝体渗流浸润线，绘出流网，通过流网求得渗流的其它水力要素。 4 ANSYS软件在大坝渗流分析应用 4.1 ANSYS分析软件简介 有限元法是工程领域应用最广泛的数值模拟方法，而ANSYS是目前应用十分广泛的通用且成熟的有限元工程分析软件之一。 ANSYS软件是一种应用范围十分广泛的大型通用软件，它在力学、流体力学、热分析等等方面的仿真模拟技术都较为成熟。ANSYSY具有的前处理模块和后处理模块以及计算分析能力功能十分强大，能够同时模拟结构、热、流体、电、磁、声学、压电以及多种物理场间的耦合效应，目前在全世界拥有广泛的用户，大量应用于土木工程、水利水电工程、汽车工程、机械、采矿、核工业、船舶等研究领域。近年来，国内外在有限元软件技术的研究开发、分析、交流、应用领域方面都大多以ANSYS为技术平台。 4.1.1ANSYS软件分析类型 作为一种大型的通用有限元软件其分析类型有:声场、压电、结构静力和动力、热传导、结构线性和非线性、电磁场以及结构屈服等分析类型。本文在研究ANSYS热传导分析的基础上模拟坝体的稳定和非稳定渗流。 4.1.2ANSYS软件的操作方式 ANSYS有两种操作方式:GUI(用户图形界面)方式、命令流方式。GUI (Graphical User Interface)即图形拾取方式，通过点击菜单项，在弹出的对话框中输入参数并进行相应的设置，从而进行问题的分析与求解。命令流方式是指在ANSYS的命令流输入窗口输入求解所需的命令，通过执行这些命令来实现问题的解答。对于操作不太娴熟的初学者来说用户界面操作方式比较适用，通过直观简单的操作步骤重复适用界面按钮执行命令，不过这种方式对于计算量大，过程复杂的求解问题是会显得繁琐和效率低。APDL是命令流方式特有的一种程序化语言操作方式，类似于一些编程语言。作为一种解释性语言能快捷有效的实现命令的重复和纠错、交互输入，通用性相对较强。依据参数可快速建立有限元模型定义所需的宏、矩阵、标量分支(if-then-do)、重复操作以及优化设计模块和网格自适应划分方面提供了一系列的便捷方式。 4.2 ANSYS模块化结构 ANSYS在仿真模拟计算中所体现出超越其它软件的优越性与其所固有的模块化结构有着必然的联系。模块化结构如图4-1所示。 预处理模块(PREP7) 求解模块(SLOUTION) 通用后处理器 (POST1) 后处理模块 时间历程后处理器 ANSYS的模块结构 优化设计模块(OPT) (POST26) 估计分析模块 (RUNSTAL) 其他模块(OTHER) 图4-1 ANSYS模块结构示意图 Figure 4-1 The structural representation of ANSYS modules 其中以下3个模块是在有限元分析过程中会经常用到的基本模块。 4.2.1预处理模块 预处理模块可以简单快捷地构建有限元模型，通过对参数的定义、几何模型 的创建及网格划分实现预处理模块的基本功能，实现有限元模型的求解模块的准确定义。 4.2.2求解模块 求解模块是对所要进行的材料基本参数及分析类型的一个定义过程，包括:定义分析类型、定义实常数、定义分析特性、定义分析选项、载荷等，接着对模型结构分析求解。 4.2.3后处理模块 POST1和POST26分别是后处理模块所具有的两大组成部分。 POST1通常称为通用后处理器，起作用是可以查看某个特定时间步内的结果，其输出的结果形式可以是数据列表形式，也可以是等值线图。时间历程后处理器(POST26)主要功能是查看建立的模型特征点，获得结果数据对时间(或频率)的列表或关系曲线。 4.3 ANSYS热-结构耦合分析简介 ANSYS热分析适用于分析材料内部或外部的温度传导状况，并导出热流场分布的一些重要参数，比如温度变化与速度的方向(thermal gradient)温度梯度、材料截面单位时间内通过的热量(Thermal Flux)热通量以及热传导过程中热量的损失和获取。能量守恒方程是ANSYS热结构耦合分析的基本原理。用有限单元法进行的网格划分，热量通过节点传导，再计算出各节点的温度，并导出其它热分析物理参数。 利用ANSYS进行大坝渗流分析主要是通过ANSYS分析的热模块去分析坝体的渗流状况。常用到的大坝渗流分析的二维或三维的单元类型有。 表4-1 二维及三维单元类型 Table 4-1 The element types of Two dimensional and three-dimensional 4.4 ANSYS热分析中热流场与渗流场的相似性 渗流场 热流场 水头H 温度T 渗透系数k，k，k，k，k热传导系数kxyz xyz 渗流速度v 热流速度v ,T,,i,,KV,,K达西定律 热流定律 ,l,l基本微分方程: 基本微分方程: ,,Z,,Z,,Z,,,,,,,,,(K)，(K)，(K),0(K)，(K)，(K),0xyzxyz,x,x,y,y,z,z,x,x,y,y,z,z 边界条件: 边界条件: ,T,,不透水层面 绝热面 ,0,0,n,n ,,f(x,y,z)T,f(x,y,z)水头函数 温度函数 等水头面φ=常数 等温面T=常数 渗流途径长度L 热流途径长度S 渗流通过的过水断面A 热流通过的过水断面A TT,,,,1212EJ水力坡降 水力坡降 ,,LS TT,,,,1212QKAQKA渗流量 热流量 ,,LS 5 ANSYS生死单元技术 ANSYS软件所具有的生死单元技术可以让已建模型中的单元根据求解要求在求解过程中的作用“消失”或“存在”，也即单元的“生”或“死”。 生死单元技术利用模型中材料属性的改变删除或加入一定的材料，促使相应类型的单元满足或不满足某一特定状况下的计算要求而选择杀死或激活所需单元，再通过单元生死选项的设置使得对单元在仿真模型中的作用做出合理判断。 5.1生死单元实现原理 在ANSYS生死单元技术实现的原理过程中，根据单元的一些在渗流场中所处位置不同对单元生死状态做出判别。其依据是会对分析区域的结果产生较为显著影响。在渗流场中，如单元内水头，载荷等。单元从有限元模型中消失，是将其特性乘以一个很小缺省因子“ESTIF”，从而忽略其在进行模型分析计算的过程中所发挥的作用。这样被乘以缺省因子的单元在载荷求解过程中发挥的作用可以忽略不计。 若果在计算中要重新启用单元对整个模型的影响作用必须重新激活他们在模型中的作用，即就是指单元的“生”。由于在后期的求解和后处理的过程中不能生成单元，所以在最初建立模型的预处理模块中生成所有单元。 5.2生死单元特性使用 ANSYS有限元分析软件特有的生死单元技术，在本文中求解渗流场的稳定渗流分析和非稳定渗流分析中做了初步尝试性的应用，相比较与其它一些数值模拟生死单元特性的使用显现出了特有的优势。 这一过程主要包括以下三个步骤: (1)建模 要在在预处理模块中，将所有单元都生成，同样包括那些可能会在以后载荷步中才会激活的单元。 (2)在所建模型上加荷载。这一过程须在求解模块中完成: ?首次对载荷步的定义: 指定分析类型 GUI:Main Menu>Solution>-Analysis Type-New Analysis-Thermal 在结构分析中打开大变形效果 GUI:Main Menu>Solution>Analysis Options 设定牛顿,拉夫森选项 Command:NROPT 后继载荷步中将要加入的单元杀死 GUI:Main Menu>Solution>-Load Step Opts-Other>Kill Elements 指定缩减因子数值(一般为10E-6) GUI:Main Menu>Solution>Other>StiffnessMult ?建立后继载荷步 为了充分发挥单元的自由度，操作者可以按自己的意愿在后继载荷步中选择 所需激活的单元，杀死载荷与约束作用为零的的单元。 (3)查看计算结果 包括单元显示、图形输出、列表输出等操作: GUI: Main Menu>General Postproc>Plot Results>Nodal Solu 5.3 ANSYS结果控制生死单元的使用 根据ANSYS计算出的温度确定单元的“生”“死”。使用命令识别并选择关键 单元，依靠ANSYS计算结果(如温度)决定杀死或需要重新激活的单元。 关键单元识别 GUI:Main Menu>General Postproc>Element Table>Define Table 关键单元选择 GUI:Utility Menu>Select>Entities 使单元“死”的效果 GUI: GUI:Main Menu> Analysis Options- Postproc >Kill Elements 使单元“生”的效果 GUI:Main Menu>Solution>-Load Step Opts-Other>Activate Elem 5.4生死单元技术确定浸润线 由于大坝坝体渗流都会存在自由水面，在浸润面以下的土体为饱和土，其渗透系数一般为常数，而处于浸润面以上的土体都处在非饱和状态。在渗流曲线附近各单元结点渗流流速为0 ，处于非饱和状态下的坝体材料其渗透系数与土体含水量形成具有一定关系的函数表达式，生死单元技术以在假定浸润线以上土体渗透系数为0 的前提下，对其不纳入分析计算，也即对坝体浸润曲线上部的单元杀死，对于要重新考虑计算的下部单元进行激活，其原理是对单元渗透矩阵的调整，这个过程是不断试算分析，不断修正提高计算精度的过程。通过生死单元技术求解坝体浸润线的方法有一定的可行性。 浸润线计算求解步骤: 第一步，对整个渗流区域渗透系数进行假定，求出在正常水位情况下，坝体内各个点的总水头。 ,n 第二步，对坝体内划分的单元结点在笛卡尔坐标系中的大小与比较，若,n果坐标值小于说明单元结点位于浸润自由面以上，通过生死单元技术给单元,n 乘以缺省因子，反之对单元进行激活处理。对模型中所有单元分析比较。 第三步，为了提高计算精度，在第二步结果的基础上重新确定渗流区域各点水头。重复以上分析过程，直到各个结点第n+1次所得水头与第n次结果差,n，1 值满足一定的精度要求即终止计算。这个反复迭代的计算过程必须通过计算机技术来实现。 生死单元技术确定浸润线以网格划分的单元为对象，通过ANSYS的网格自适应功能对坝体水头变化较大的区域以及在浸润线附近较大区域的进行局部网格加密的方法，来充分提高计算精度，获得较为理想的光滑的浸润线。这种不断迭代的方法需要反复进行大量重复步骤操作，需要较长的时间和代价。但是，采用ANSYS提供的命令流(APDL)语言将模型参数化，会使问题得到进一步的 简化。 6基于ANSYS软件的坝体稳定渗流分析过程 坝体的稳定渗流状态分析属于ANSYS中的静态(static)分析问题。在进行坝体非稳定渗流分析以前，通常先进行坝体的稳定渗流分析，确定坝体内部初始水头分布。ANSYS分析坝体稳定渗流的一般过程分为建模、施加载荷、后处理三大模块。 首先在建立坝体仿真模拟结构的同时选择单元类型，设定分析单元大小选项。定义坝体材料参数，在稳定渗流分析中则可设置材料属性的渗透系数，最后创建几何模型进行网格划分。其次，在已建立的坝体模型上施加边界条件:上下游水位恒定的水头、以自有约束施加的边界流量、在渗流时段内库区外界条件的体载等作为载荷直接作用与实体模型上。在后处理中AYSYS计算结果会直接将渗流分析所得数据写入*.rth文件中，其所得结果可有三种查看方式:彩色云图显示(等水头图)、坝体内流速矢量图、结果列表显示;输出结果包括:坝体内各个节点水头、渗流密度、渗流梯度等。 7利用简单算例对ANSYS软件求解渗流问题进行求证 7.1利用有限单元法求解渗流场 如图所示的一个矩形土坝算例来示范单元渗透矩阵及总渗透矩阵形成，求解在上下游固定水位条件下的稳定渗流场，通过划分少数单元来具体说明计算过程。 z 13(0,6)试验点 5(1,5.5) 试验点8 6(2,4.5)9 10 117 7(2,3.5)2(4,1)6 1213 12(0,3)8(4,2.5)3(2,3) 4(3,2)14 54 3 9(4,1)1 7-1矩形土坝网格划分 2 x Figure 7-1 The mesh of Rectangular earth dam11(0,0)10(4,0)1(2,0) 将渗流场以三角形单元进行划分，其中1至8号节点为未知水头结点，9至14号为已知水头结点。为使形成总矩阵的带宽最小将自由面上的结点和单元编成连号，三角形单元和结点数据如下表。 表7-1 三角形单元与结点数据 Table 7-1 The Triangular element and Node data 三个结点 坐标 单元号 结点号 i j m x y 1 1 12 11 1 2 0 2 1 10 9 2 1 4 3 1 9 4 3 2 3 4 1 4 3 4 3 2 5 1 3 12 5 1 5.5 6 2 12 3 6 2 4.5 7 2 13 12 7 3 3.5 8 2 5 13 8 4 2.5 9 2 6 5 9 4 1 10 2 3 6 10 4 0 11 3 7 6 11 0 0 12 3 4 7 12 0 3 13 4 8 7 13 0 6 14 4 9 8 通过上面的数据，即可计算各个三角形单元的矩阵。再将各单元矩阵中的重复单元的渗透矩阵系数相加，得到总渗透矩阵的各个元素，最终形成总渗透矩阵〔K〕。其次计算方程组中的常数项。最后求得未知结点水头的矩阵形式: ,240,2,80000-59,,,,,,1,,,,,,,52,90,10,1200-1262,,,,,,,,,,,,,51,150,4,120-543,,,,,,,5000,4,12-11,,,,,,4 ，,0,,,,,,,对25,300-725,,,,,, ,,,,,,,称22,3006,,,,,,,22,307,,,,,, ,,,,,,,17-28,,,,,, 求解以上方程组就可得出1至8结点的水头值。按照常用的直接消去法和迭代计算得出的计算结果与甘油模型试验结果的自由面位置相比较，如下表。 表7-2 计算值与试验值结果 Table 7-2 The results of calculation and experimental 结点水头 ,,,,, ,,,35678124 3.69 5.47 4.58 3.37 5.67 5.24 4.47 3.40 直接消去法 计算结果 3.94 5.46 4.54 3.31 5.65 5.22 4.44 3.33 迭代法计算 结果 5.63 5.10 4.33 3.25 甘油模型试 验结果 7.2利用ANSYS软件对渗流场的仿真模拟计算 在ANSYS软件中利用渗流场与热流场具有相似特性的关系，仿真模拟矩形土坝渗流场，得出结果与传统的有限单元法计算结果以及试验研究结果进行对照。 7-2矩形土坝有限元网格划分 7-3矩形土坝坝体内水头等势图 The finite element grid division of Rectangular earth dam The potentionmetric diagram of water head H in the dam 在利用有限单元法以矩形土坝为简单算例求解其渗流场，通过迭代法和直接 比较，其获得的未知结点水头值法得出的未知结点水头值与试验所得结果相分析 基本一致，误差较小，而且若果增加划分的单元并合理安排，会得到更高精度的计算结果。再将以上结果与通过有限元ANSYS软件模拟计算求解值对应分析，所得水头值也是一致的，而且模拟计算更容易进行单元网格的划分，在边界条件比较复杂时可以自动加密，使得计算精度单位提高。 8基于ANSYS的稳定渗流分析求解实例 本文利用甘肃省临洮县海甸峡水电站工程为例。水电站工程由引水枢纽、输水系统及发电站厂房三部分组成。枢纽为闸坝结合型式，右岸挡水坝段采用混凝土面板堆石坝，左岸挡水坝段采用混凝土重力坝，坝顶高程2005.00m，最大坝高49m;水库正常蓄水位2002.00m，校核洪水位2003.70m，水库总库容2200 33万m,调节库容200万m。实例计算主要为右岸挡水坝段(混凝土面板堆石坝)，坝段长160m,坝顶标高2004.50m，防浪墙顶高程防浪墙顶高程2005.20m，趾板建基高程1956.00m，设计坝高45m，坝顶宽5.0m。坝体上游边坡比为1:1.4，下游边坡比为1:1.3—1:1.4，下游坝坡坡比分为两级，高程1983.10m以上为1:1.3，该高程以下为1:1.4。下游坝坡以0.3m厚的干砌块石护坡，在高程1983.10m处设置马道，马道宽度2.0m。上游坡面C25钢筋混凝土面板，厚0.45m;下铺垫层料和过渡料，垫层区采用导流洞开挖石渣和天然河砂制备，根据施工要求其水 平宽度为3.0m，过渡区填料以洞挖石料为主，水平宽度3.0m。混凝土面板上游设土质斜墙铺盖及碎石土盖重体，其中:土质斜墙铺盖厚度为1.50m，土质斜墙铺盖及盖重体顶部标高1975.0m。坝顶填2.0m厚砂碎石，表面铺30cm厚混凝土路面。坝顶上游侧设“L”型钢筋混凝土防浪墙，墙高2.70m，墙顶高程2005.20m，墙体厚度0.35m。其面板堆石坝段桩号0+151坝体断面形状见图8-3。 表孔溢流坝段采用上游面为铅直的基本三角形断面，堰顶高程为1994.00m，建基面高程1962.10m，最大底宽38m。溢流坝为2孔，孔口净宽10m。堰顶上游堰头为双圆弧曲线，下游堰面采用WES幂堰面曲线。坝体采用现浇C15混凝土，堰面设0.2m厚的C50硅粉混凝土。闸墩采用现浇C20钢筋混凝土，中墩厚3m，边墩厚2m。工作闸门及检修闸门均采用钢制平板闸门。消能形式为底流消能，消力池长46m，宽23m，底板高程1970.20m，采用现浇C20钢筋混凝土，厚2m，面层设0.2m厚的C50硅粉混凝土防冲刷。 坝体主要材料分区为:垫层、过渡层、主堆区、次堆区。上游水位为2002.00m，下游水位为1972.00m;混凝土面板，渗透系数为1.0E-7;垫层，渗透系数为2.0E-3;过渡层，渗透系数为2.0E-3;主堆区，渗透系数为5.0E-3;次堆区，渗透系数为1.0E-5;基岩弱透水层5.0E-4。 工程的溢流坝、泄洪冲沙闸和挡水坝等主要建筑物的平面布置如图所示8-1，上游立视图如图所示8-2。 124708dlQkstk420480 52889k1ssh-ss1 1000500053481左岸截渗墙段左岸副坝段泄冲闸段溢流坝段1ssh-ss25002500250030002500250030002500100001000080008000800020000240001450018500137443004008004000800250010003000116002500940094003750R150003000790012000143.2555?F41cgk15000dlQkst53.2555?4alQfzrtk1ssh-ss41 5000k1ssh-ss1k1ssh-ss1alQspgr41cgk1dlQ450001cgk1k1ssh-ss1R7500 kst 8-1海甸峡水电站平面布置图 Figure 8-1 The Floorplan of Hai Dian Xia Hydropower 8-2海甸峡水电站上游立视图 Figure 8-2 The upstream elevation of Hai Dian Xia Hydropower 图8-3洮河海甸峡大坝0+151坝体断面示意图 8.1图形界面操作方式(GUI) (1)在工作名框中输入:Haidianxia; (2)在下边ANSYS标题框中标题输入:Analyzing the seep-field of the Taohe Haidianxia Dam Using ANSYS; (3)在分析计算中使用国际统一单位:/units.SI; (4)在分析类型中，选择热分析类型(Thermal); (5)选择的单元类型如下图所示; 在选择单元类型时选择热分析实体单元plane35(三节点三角形单元)，理由是在对渗流场复杂边界较强的适用性，而且也好控制。 (6)MainMenu:Preprocessor>Material Prop> Material Models>Define Material Model Behavior>Materrial >New Model >Define Material ID ;新建6种材料编号。在Material Models Available;Thermal>Conductivity>Isotropic>Conductivity for Material。在热流传导框中输入面板材料的渗透系数，KXX框中输入1.0E-7的渗透系数，选择应用;增加第二种材料编号2，KXX框中输入垫层的渗透系数2.0E-4，选择应用;增加第三种材料编号3，KXX框中输入过渡层的渗透系数2.0E-3，选择应用;增加第四种材料编号4，KXX框中输入主堆区渗透系数5.0E-3，选择应用;增加第五种材料编号5，KXX框中输入次堆区渗透系数1.0E-4，选择应用;增加第六种材料编号6，KXX框中输入5，基岩弱透水层的渗透系数5.0E-4， 选择结束; (7)Main Menu:Preprocessor>Modeling>Create>Keypoints>In Active CS里面定义关键点; (8)Main Menu:Preprocessor>Modeling>Create>Lines> Lines>Straight Line，由关键点连成线; (9)Main Menu:Preprocessor>Modeling>Create>Areas>Arbitrary>By Lines, 由线创建面; (10)Main Menu:Preprocessor>Numbering Ctrls>Compress Numbers,选择All,单击OK,压缩模型元素的编号; (11)Main Menu:Preprocessor>Modeling>Operate>Booleane>Glue>Area,选择PICK ALL，将创建的的相邻两个面的边界粘合; (12)在主菜单Preprocessor中的单元面积大小进行控制，先选择堆石面板坝混凝土面板,在尺寸框中输入单元边长2,选择APPLY,选择心墙上垫层,在SIZE框中输入单元边长3,选择APPLY,选择过渡层,在SIZE框中输入单元变长3,选择APPLY,选择主堆区,在SIZE中输入单元变长4，选择次堆区,在SIZE框中输入3, 选择基岩若透水层,在输入框中定义的单元边长4选择OK; (13)在主菜单界面的Preprocessor中定义材料的属性，对已经编号的坝体分区材料进行划分，选择面板,在MAT框中输入1,选择APPLY,选择垫层,在MAT框中输入2,选择APPLY,选择过渡层,在MAT框中输入3,选择APPLY,选择主堆区,在MAT中输入4，选择次堆区,在MAT框中输入5, 选择基岩若透水层,在MAT框中输入6,选择OK; (14)Main Menu:Preprocessor>Meshing>Mesh>Areas>Free,选择PICK ALL; (15)在求解模块中对坝体单元加入载荷计算，在上游坝面选择正常水位以下的节点单元进行加载，上游水位2002m;选择下游正常淹没水位节点单元输入水头1972,进行求解; (16)在求解过程中会出现Current LS的问题，直接结束便可; (17)在通用后处理器中对所得结果处理输出，一般为结点温度(水头)以彩色云图显示—等势图; (18)Main Menu>General Postproc>Plot Results>Pre-defined or Userdefined， 选择OK;矢量图显示—流速矢量图; 8.2命令流APDL方式 利用/CLEAR命令对界面进行清除 /定义的文件名为:Taohe ANSYS操作界面定义渗流分析的标题:Analyzing the Seep-field of the TaoheHai dian xia Dam Using ANSYS uplevel=2002.00! downlevel=1972.00! kx1=1.0e-7 ky1=1.0e-7 kx2=2.0e-4 ky2=2.0e-4 kx3=2.0e-3 ky3=2.0e-3 kx4=5.0e-3 ky4=5.0e-3 kx5=1.0e-4 ky5=1.0e-4 kx6=5.0e-4 ky6=5.0e-4 创建坝体有限元模型 进入预处理模块 在这一步中选择适宜单元类型划分的模型网格:plane35(Triangl 6node 35) 能够适应复杂的渗流边界条件而且比较容易控制。 mp,kxx,1,kx1 mp,kyy,1,ky1 mp,kxx,2,kx2 mp,kyy,2,ky2 mp,kxx,3,kx3 mp,kyy,3,ky3 mp,kxx,4,kx4 mp,kyy,4,ky4 mp,kxx,5,kx5 mp,kyy,5,ky5 mp,kxx,6,kx6 mp,kyy,6,ky6 坝体模型关键点的创建以及由坝体模型关键点创建线、面; a,p1,p2,p3,,,,,～面板 a,,,,,～垫层 a,,,,,～过渡层 a,,,,,～主堆区 a,,,,,～次堆区 a,,,,,～基岩弱透水层 aglue,all /pnum,area,1 网格控制与划分 aesize,5,3 amesh,5 mat,6 aesize,5,4 amesh,6 NROPT,Full OUTPR,BASIC,All 加载荷并且求解 nsel,s,ext nsel,r,loc,x,0,1972 nsel,r,loc,y,(271.50+1e-5),downlevel d,all,temp,downlevel solve 图形查看计算结果 /post1 plnsol,temp FINISH 8.3计算结果分析 坝体有限元计算模型材料分区如图，坝体有限元计算模型网格划分如图;计 算结果参见坝体内水头等势图和坝体内流速矢量图 图8-4 坝体材料分区图 Figure 8-4 The materials division of dam 图8-5坝体有限元网格划分图 Figure 8-5 The finite element grid division diagram of dam 图8-6 大坝渗流浸润曲线图 图8-7 坝体内流速大小与方向示图 本文大坝坝体有限元模型共划分单元总数是371，节点总数是809 。本文通过ANSYSY软件自带的生死单元技术杀死压力值小于零的非饱和区域中的单元，将饱和区域大于零的单元进行激活。这样就会在非饱和区域与饱和区域的界面上存在一个临界区域，在这个区域面内各个单元节点的水头线是等势的，这就是所谓的渗流曲线。在上图8-6中，堆石面板坝体内的水头在通过面板时下降过程不太明显，水头有可能会有个增幅，在次堆区内其降幅明显 说明坝体构筑材料渗透系数比较小的坝体区域会对大坝的防渗有着非常显著的作用。通过图8-7可以看出，透过上游坝面，即混凝土面板和主堆区的流速相对于下游坝体时流速都较小，而在下游坝面与下游水位交界点处有最大流速出现，这对水工建筑物的运行较为有利的，可以将积存在坝体内的水位得到有效的降低，降低扬压力，确保大坝的正常运行。 9利用ANSYS软件的坝体的非稳定渗流计算 在基于上面稳定渗流分析成果的基础上，研究在上游不同水头的作用下载荷变化对渗流场的影响，其主要区别是载荷(库水位)随时间的变化引起的。 图9-1载荷—时间关系曲线 在上述关系曲线中随着时间的变化形成一个载荷步，在上图9-1中任何一个突变处就是载荷发生变化的地方。在模型分析计算的每一个载荷步中，每一个载荷值都会有一个对应的的时间值，同时必须指定的施加载荷步的方式为渐变或阶越。 9.1非稳定渗流分析的步骤 9.1.1建模 1确定工作名(JOBMANE)、标题(TITLE)、UNITS后，进入预处理模块结构(PREP7); 2在稳定渗流计算中使用的是三节点三角形单元(plane35)可以在非稳定渗流中继续定义。 3定义单元实常数;(除定义和稳定渗流分析相同是常数外，还需定义坝体材料的弹性模量、泊松比等) 4定义坝体主要材料渗流性能:非稳定渗流分析要定义给水度及渗透系数; 5建立坝体主断面的有限元模型; 6对生成的几何模型进行单元网格划分。 9.1.2加载求解 第一步:分析类型的定义(如图9-2)，指定分析类型为瞬态分析(即指渗流场中的非稳定渗流分析)，并且可以对这个分析过程进行重新的启动。 ?首先在主菜单的求解模块中定义分析类型。 ?如果在增加其它载荷的情况下接着上次的分析继续进行。 GUI:Main Menu>Solution>Analysis Type>Restart 图9-2定义非稳定分析类型(GUI方式) 第二步:对于非稳定渗流分析获得初始条件 在进行了稳定渗流分析后，将其稳定渗流结果作为非稳定渗流分析的初始条件进行的渗流分析。坝体的初始水头分布是通过假定先初设以某一值，在计算时可以通过命令对初始条件选择。 通过以下命令来施加初始条件 Main Menu>Solution>Define Loads>Apply>Initial Condition>Define命令来施 非稳定分析有以下方法指定具体的初始条件: ?以水头、渗流量等设定合适的稳定载荷 ?在载荷步的设定时，首先考虑子步的载荷的为一极小值。在主菜单Preprocessor子目录的载荷Loads中设定Time/Frequence和Time and Substps。 ?在进行程序菜单命令操作时选择关闭时间积分的设置。 ?通过时间—载荷命令(KBC)指定载荷的变化类型:渐变或阶越。 要是载荷值随时间呈线性变化过程，就将其设定为渐变选项;要是载荷在这个载荷步是恒定的，就将其设定为阶越选项;通过以下命令实现: Command:KBC ?进行求解 Command:SOLVE ?写入载荷步文件: Command:LSWRITE 这个过程中需要注意的是:如果在稳定渗流分析中设定的水头有与后边非稳定渗流分析中相同的节点水头的将其保留，否则删去所有设定的水头。 第三步:设定时间—载荷步选项 在进行非稳定渗流分析需要指定时间—载荷步选项主要有以下三项:普通选项、非线性选项以及输出控制选项。 ?普通选项;时间:本选项设定每一载荷步结束时的时间。 Command:TIME 在非稳定渗流分析过程中，随着时间的变化上游水位在下降过程中的载荷会形成许多子步，而且载荷子步越多，步长越小计算输出的结果值精度就越高，反之步长越大计算进度就越低。 阶跃载荷选项(KBC)，根据载荷变化情况分为阶变和渐变。当载荷值随时间线性变化时，将其设定为渐变选项，要是载荷在这个载荷步是恒定的，需要设 为阶越选项。 可以用下面的命令或GUI操作方式来实现。 Command:KBC ?非线性选项:对于大多数非稳定分析来说，系统默认的迭代次数已经能够满足计算要求。要是分析的问题比较难于收敛时，也可以通过以下命令来指定增加迭代次数。 Command:NEQIT 时间积分效果—要进行坝体的非稳定渗流分析，将非线性选项设置为ON。 ?输出控制选项:在定义完成所有求解分析选项后，进行求解，输出控制结果。 9.1.3后处理(POST1、 POST26) 对于非稳定渗流分析和瞬态热分析一样，ANSYS通同样提供通用后处理器(POST1)和时间—历程后处理器(POST26)两种后处理: (1)通用后处理器(POST1)，在对坝体仿真模型的某一时间段内的渗流成 果加以不同载荷并对后续结果进行处理。 Command:POST1 (2)时间-历程后处理器(POST26)，在非稳定渗流中对水位和时间变化对坝体内渗流结果进行优化。 第一步采用通用后处理器进行后处理 ?进入通用后处理器后，就可以将某一时间的结果读出。并且AYSYS会自动通过线性插值，可以设定不在任何一个子步上的时间点。 Command:SET ?还可以输出某一载荷步的结果，同上面实例中的输出形式一样，对非稳定渗流分析结果可以有三种方法进行输出。 第二步采用时间—历程后处理器进行后处理 ?首先定义变量:上游水位随时间变化。 Command::NSOL / ESOL / RFORCE ?输出列表: Command:PRVAR 或绘制出这些变量随时间变化的曲线图: Command:PLVAR 9.2上游水位变化过程的表型数组定义过程 当上游水位变化时在ANSYS软件对非稳定渗流分析必须定义变量，而变量的定义一般数组的形式出现，如以下几种数组参数: 9-3数组的定义 Figure 9-3 The definition of array parameter 9-4 数组的赋值 Figure 9-4 The valuation of array parameter (2)表数组特点 表数组类型的参数是数值型数组参数和数值型数组类似。 表述组的下标值可以是实数或者为递增的整数。在表数组的赋值过程中根据数组的“列”和“行”的下标值进行赋值，如果表述组不需要赋值的则程序会自 动赋一极小值。 表数组的最大特点是在计算过程中允许在已定义的表数组元素之间通过线性插值插入任意的值。因此在进行线性插值计算的功能是按照行和列的下标进行的插值，这就允许在非稳定渗流中用于定义随时间变化的载荷或者复杂多变的边界条件。 本文就是充分利用这一特点的基础上，在表数组中定义各个时间段得上游水位下降过程，将其定义为UPLEVEL，来实现非稳定水头变化的载荷施加，来实现非稳定渗流分析的操作。 10非稳定渗流分析计算求解实例 以下利用海甸峡水电站工程实例来说明ANSYS求解非稳定渗流问题过程。 例如:海甸峡水电站的库水位在蓄水期一般都是处于正常水位2002.0m高程，需要水库下泄流量时库区水位会以每天大约2米的速度下降，在上游水位下降的过程中保持下游水位变幅保持不变为1972.0m，在研究过稳定渗流的基础上对坝体非稳定渗流状况分析。水位下降过程如图10-1。 上游水位下降过程 2004 2002库水位2000线性 (库水位)1998 1996 1994 1992 1990 1988上游水位()m1986 1984 1982 1980 024487296120144168192216240264288 时间(h) 10-1水位下降与时间关系曲线 Figure10-1 The relation curve of water level declining and time 10.1图形菜单GUI方式与APDL方式交互求解非稳定渗流过程 (1)在实用程序菜单文件夹的工作名中输入:Haidianxia; (2)实用程序菜单文件夹的标题栏中输入:Analyzing the seep-field of the Taohe Haidianxia Dam Using ANSYS; (3)计算过程一致采用国际通用单位，输入命令:/units,SI; (4)上游水位UPLEVEL表数组的定义:Utility Menu> Parameters> Array Parameters> Define/Edit>Add Array new Parameters 10-2 UPLEVEL数组的定义 Figure 10-2 The definition of table parameter UPLEVEL (5)上游水位UPLEVEL表数组的赋值:Utility Menu> Parameters> Array Parameters> Define/Edit>Table Array:UPLEVEL=f (TIME，LEVEL) 10-3 UPLEVEL数组的赋值 Figure 10-3 The valuation of table parameter UPLEVEL 将上游水位数据从文件uplevel.txt中读入到数组(12,1) (6)Main Menu:Preprocessor，选择Thermal; (7)在主菜单Preprocessor中单元类型中添加或编辑渗流分析所需的单元类型，与稳定渗流相同都选择Plane35。 (8) Model Behavior>Materrial >New Model >Define Material ID ;新建6种材料编号。在Material Models Available;Thermal>Conductivity>Isotropic>Conductivity for Material。 对坝体各个分区材料的渗透系数的定义，在上图所示的热流传导框中输入面板材料的渗透系数，KXX框中输入1.0E-7的渗透系数，选择应用;增加第二种材料编号2，KXX框中输入垫层的渗透系数2.0E-4，选择应用;增加第三种材料编号3，KXX框中输入过渡层的渗透系数2.0E-3，选择应用;增加第四种材料编号4，KXX框中输入主堆区渗透系数5.0E-3，选择应用;增加第五种材料编号5，KXX框中输入次堆区渗透系数1.0E-4，选择应用;增加第六种材料编号6，KXX框中输入5，基岩弱透水层的渗透系数5.0E-4，选择结束; (9)Main Menu:Preprocessor>Modeling>Create>Keypoints>In Active CS里面定义关键点; (10)Main Menu:Preprocessor>Modeling>Create>Lines> Lines>Straight Line，由关键点连成线; (11)Main Menu:Preprocessor>Modeling>Create>Areas>Arbitrary>By Lines, 由线创建面; (12)Main Menu:Preprocessor>Numbering Ctrls>Compress Numbers,选择All,单击OK,压缩模型元素的编号; (13)进行布尔运算，对不同材料的面进行结合。 ，选择所有创建的面域，将相邻两个面的边界粘合。 (14)在菜单操作中选择不同坝体材料输入相对应的单元变长的大小，首先选择面板,在SIZE框中输入单元边长2,选择APPLY,选择心墙上垫层,在SIZE框中输入单元边长3,选择APPLY,选择过渡层,在SIZE框中输入单元变长3,选择 ，选择次堆区,在SIZE框中输入APPLY,选择主堆区,在SIZE中输入单元变长4 3, 选择基岩若透水层,在SIZE中输入单元变长4，最后结果一上定义过程; (15)对以划分的坝体各个分区进行不同的编号进行选择划分网格，选择面板,在MAT框中输入1,选择APPLY,选择垫层,在MAT框中输入2,选择APPLY,选择过渡层,在MAT框中输入3,选择APPLY,选择主堆区,在MAT中输入4，选择次堆区,在MAT框中输入5, 选择基岩若透水层,在MAT框中输入6,选择OK; (16)Main Menu:Preprocessor>Meshing>Mesh>Areas>Free,选择PICK ALL; 10.2非稳定渗流分析 APDL命令流: (1)antype,transient (2)trnopt,full(非稳定完全分析) (3)打开大几何变形开关:!NLGEOM,1 (4)做稳定分析，确定非稳定分析的初始条件:timint,off/ time,1e-6/ deltim,1e-6(设定只有一个子步的时间很小的载荷步) (5)设定上游水位节点载荷: Command:NSEL,Item-S,Comp-EXT Command:TEMP -r,Loct,X,Y,-2002,0 Command: y,(271.50+1e-5),UPLEVEL! Command; D,ALL,TEMP,UPLEVEL! (6)设定下游水位节点载荷: Command:NSEL,Item-S,comp-EXT Command:R,Loc,X,0,1972,0 Command:NSEL,R,Loc,Z,(271.50+1e-5),DownLevel Command:D,All,TEMP,DownLevel solve(求解出渗流场的分布) (7)打开时间积分:timint,on (8)关闭自动时间步长:autots,0 (9)记录每一子步数据文件: Command: kbc，0 Command: outres,all,all (10)读取上游水位数组值及上下游水位节点载荷 将上游水位定义为UPLEVEL表数组，定义每一载荷步的时间值TIME Command: *DOF,Lab-1,Lab-1,Lab-12,Lab-1 定义随时间—水位的变化 Command: nsel,r,loc,x,-2002,0 将*rth文件中的上游水位表数组读入计算模型 Command: nsel,r,loc,y,(271.50+1e-5),UPLEVEL(1,2) 上游水位下降过程中下游水位高程 Command: nsel,r,loc,x,0,1972 将*rth文件中的下游水位表数组读入计算模型 Command: nsel,r,loc,y,271.50+1e-5),downlevel (11)删除节点温度(稳定分析) delet,all,temp (12)求解 solve (13)结束求解并退出 Command: *ENDDO Command: finish (14)进入后处理:与预处理相对应，对坝体非稳定渗流的后序工作 定义结点变量:Command:ESOL 绘制曲线或输出列表:Command:PLVAR/PRVAR 10.3堆石面板坝非稳定渗流计算结果分析 10-4 坝底11个点的水头响应随时间的变化曲线 Figure 10-4 The 11 points of fluctuation curve under waterhead and time 10-5 水头沿路径的变化曲线 Figure 10-5 The fluctuation curve of water head and time 在对大坝进行非稳定渗流分析时，随着上游库水位以每天2m的速度下降的 过程中坝体内的渗流水头响应随时间的变化称曲线形式。通过上图中坝底11点 渗流水头曲线的变化变化可以看出，处于坝体上游迎水面坝壳上的2个点的水头 曲线的斜率变化比较明显，与上游水位的变化基本相近，介于坝体上游坝壳与下 游坝壳之间的4个点的水头的曲线变化斜率比较小，逐渐过渡到下游坝壳5个点 时水头曲线的斜率几乎不发生变化保持为0，总体来看以上11个点的曲线斜率 在库水位逐渐下降的过程中从上游到下游方向渐渐变小。 结 论 本文主要通过有限单元法渗流分析的基本原理，根据渗流场与其温度场的相似性，利用大型有限元分析软件ANSYS，并以其热结构耦合分析模块对洮河海甸峡堆石面板坝在稳定状态和非稳定状态下渗流分析的研究。研究结果表明利用ANSYS软件中的热分析模块对大坝进行渗流分析的研究方法和思路具有一定可行性、可靠性。同时，在软件的应用过程中发现其命令流方式(APDL)形成的参数化语言要比通过窗口菜单模式(GUI)有易保存、易修改的优点。 对堆石面板坝的渗流研究以渗流的基本理论出发，介绍多孔介质渗流概本文 念以及不同坝体材料渗透系数。同时对渗流分析计算的基础理论、基本微分方程、连续性方程(控制方程)以及渗流有限单元法的基本原理进行了分析、推导，得出可以进行稳定渗流场理论计算依据。这为利用大型通用有限元分析软件ANSYS的仿真模拟计算、分析坝体渗流问题打下了坚实的理论基础。 在本文中还介绍了有限单元法的基本思路和原理，并对有限单元法如何求解微分方程重点进行了推导演算，论述了有限单元法求解二维平面渗流问题的基本过程。通过一个简单算例—矩形土坝的渗流分析对ANSYS热分析计算坝体渗流问题可行性进行了验证。在对洮河海甸峡面板堆石坝实例进行的稳定渗流分析计算时，利用ANSYS特有的生死单元技术来确定坝体内渗流水头位置及坝体内各点的渗流流速，其结果直观、精确并以彩色云图或列表形式给出。通过对海甸峡面板堆石坝非稳定渗流的定性分析，当上游库水位以2m/d的速度下降时，研究了大坝底部11个位置水头曲线斜率的变化，处于坝体上游迎水面坝壳上的2个点的水头曲线的斜率变化比较明显，与上游水位的变化基本相近，介于坝体上游坝壳与下游坝壳之间的4个点的水头的曲线变化斜率比较小，逐渐过渡到下游坝壳5个点时水头曲线的斜率几乎不发生变化保持为0。以上实例求解结果表明利用有限元分析软件ANSYS的人结构耦合模块对面板堆石坝的渗流分析具有的可行性，这为我们以后在面板堆石坝的坝体渗流分析和研究提供了新的研究思路和方法，具有一定的参考价值。