2.1 无限深透水地基上带有防渗斜墙和水平铺盖的土石坝有限元渗流计算
在渗透系数各项异性的深厚覆盖层或非均质无限深透水地基上，存在着很多限制条件和技术难度。而且在实际工程中，水平防渗体无论是利用黏土还是土工膜，都会因为人为因素，或不可抗力等原因而存在渗漏的问题。因此，本节对非均质无限深地基上带有不同长度的微透水水平铺盖和斜墙防渗体的土石坝，进行了有限元渗流分析计算，并得出以下结论：当铺盖长度是坝前水深的8倍时，土石坝坝后渗透坡降小于0.1，满足渗透稳定的要求；当铺盖长度是坝前水深的26～30倍时，土石坝坝体、坝基渗流量基本趋于稳定，此时若再通过增加铺盖长度的方法减少渗流量，其效果已经不明显；当水平铺盖长度是坝前水深30倍时，比规范规定5～8倍水头的设计节水35%～42%。
静力有限元法一般可以用来处理土石坝非均值无限深透水地基上关于渗流的问题。这种方法是通过综合考虑上游水头、外部静荷载等因素的影响，从而分析坝体有关渗透量和渗透特性方面的问题。
2.1.1 有限元法概述及在渗流问题中的应用
有限元法是伴随着计算机技术的发展而形成的一种数值分析方法，它是应用数学、力学以及计算机技术相结合的产物。有限元法的基本思想是：将连续的求解域离散为一组单元的组合体，用在每个单元内假设的近似函数来分片地表示求解域上待求的未知场函数，近似函数通常由未知场函数及其导数在单元各节点的数值插值函数来表达。从而使一个连续的无限自由度问题变成离散的有限自由度问题。
作为一种离散化的数值解法，有限单元法首先在结构分析，然后又在其他领域中得到广泛应用。离散化的思想可以追溯到20世纪40年代，1941年A.Hrennikoff首次提出用构架方法求解弹性力学问题，当时称为离散元素法，仅限于用杆系结构来构造离散模型。如果原结构是杆系，这种方法是精确方法，发展到现在就是大家所熟知的结构分析的矩阵方法，究其实质这还不能够说是有限单元法的思想。1943年Rcourant在求解扭转问题时为了表征翘曲函数而将截面分成若干三角形区域，在各三角形区域设定一个线性的翘曲函数，这是对里兹法的推广，实质上就是有限单元法的基本思想，这一思想真正用于工程中是在电子计算机出现后。20世纪50年代，由于航空业的需求，波音公司的专家第一次用三个结点的三角形单元，把矩阵位移法用在一个平面的问题上。与此同时，联邦德国斯图加特大学的J.H.Argyris教授发表了一些能量原理与矩阵分析的论文，这种方法对理论基础研究作出了杰出的贡献。1960年美国的R.W.Clough教授首次在一篇题为“平面应力分析的有限单元法”的论文中提到有限单元法（The Finite Element Method）一词，后来这一名称得到广泛承认。20世纪60年代，随着有限元方法的迅速发展，除了力学学者，许多数学家也参与了这项工作，奠定了有限元方法的理论基础，找出有限元法和变分法之间的关系，开发出了许多的单元模式，使有限单元法的应用范围得到扩大。20世纪70年代之后，有限元法进一步迅速发展，其应用范围扩展到工程的各个领域，成为连续介质问题的数值解法中最常用的分支，由变分法有限元扩展到加权残值法和能量平衡法有限元，从弹性力学的平面问题扩展到空间问题、板壳问题，从静态平衡问题扩展到稳定问题，动力问题及波动问题，从线性扩展到非线性问题，研究分析的对象由原来的弹性材料扩展至塑性、黏弹性、黏塑性及复合材料，研究目标由单纯的结构分析扩展到结构优化甚至设计自动化，研究领域由固体力学领域扩展到流体力学、传热学、电磁学等领域。伴随着计算机的发展，有限元法中的繁琐的计算过程逐渐由人工计算转为计算机代替计算，这大大节省了人力和时间，科技人员又将计算机技术、数学理论与有限元理论相结合，开发出许多有限元通用软件，常见的有限元软件有ANSYS、ADINA、ABAQUS及NASTRAN等。
关于有限元法的基本原理，简单而言，有限元法研究就是将连续体离散化为有限个单元的集合体来进行的。通过变分原理建立模型，推导出近似解的一组方程，最终转化为求解多阶系数矩阵的线性方程组。这类计算主要利用计算机进行。
对于渗流计算，由于无压渗流具有渗流自由面（浸润线），且非稳定渗流自由面随着上游库水位升降而变动，而且一般的渗流场都有不同程度的非均质和各向异性，几何形状和边界条件较复杂，求解解析解在数学上存在着不少困难，仅仅能对一些简单的流动情况获得解析解。20世纪60年代后，电子计算机的普及和数值计算方法的发展，特别是有限元法提出后，推进了渗流数学模型的发展。有限元法能系统地编成计算程序，很方便地处理复杂的边界条件，非均质土层，以及第二类流量边界不需专门处理而能自动满足，故目前在渗流计算中得到广泛的应用。
用有限元法求解渗流问题，一般有以下几个步骤：
（1）连续区域离散化。将需要研究的渗流区域分解成若干相互连接的单元，单元的划分以计算机性能和精度要求为基础。
（2）插值函数的选择。单元内的水头分布采用一些简单而适用的函数表示。
（3）推导特性矩阵方程。
（4）组合各个单元的渗透矩阵集合形成整个渗流区域的总体渗透矩阵［K］，其代数方程的形式为
式中 {H}——渗透区各节点水头值的列向量；即
式中 n——内节点和第二类边界条件上的节点数（即未知节点数）；
{F}——自由项（已知项）组成的列向量。即
（5）代数方程组的求解。利用消元法和迭代法求解各个节点的未知水头值。
（6）由节点水头值计算水力比降、渗流量以及其他需要计算的物理量。
2.1.2 数值模型建立
某水库为灌注式平原水库，坝基为渗透系数不同的深厚粉细砂覆盖层，渗流设计按照无限深透水地基考虑。坝前最大水深为8m，坝基70m之上和100m之下为较强透水层，渗透系数k1为6.26×10-3cm/s；70～100m为较弱透水层，渗透系数k2为2.50×10-4cm/s，允许渗透坡降为0.1。水库的坝体为土石坝，坝高为12m，坝顶宽度为6m，坝底最大宽度54m，上下游边坡比均为1∶2，采用复合土工膜防渗结构作为水平防渗体和斜墙防渗体，坝后地下水位与地面同高。
以水平防渗体长度为32倍水头为例，建立非均质无限深透水地基微透水水平铺盖的有限元计算模型。坝基为长方体（800m×100m×224m），上游坝脚至上游边界600m，下游坝脚至下游边界146m。坝基70m之上为第Ⅰ层，100～224m之下为第Ⅲ层，渗透系数k1为6.26×10-3cm/s；70～100m为第Ⅱ层，渗透系数k2为2.50×10-4cm/s；考虑到施工时人为破坏的影响，水平防渗体和斜墙防渗体的渗透系数k3取1×10-7cm/s，具体渗透系数和密度见表2.1。
表2.1 坝体和坝基的参数
由于初次对该问题进行研究，且该模型参数在X、Y方向上是均匀的，因此只对模型进行了二维的分析。对全局每5m进行一个单元的划分。坝体和坝基共有7330个节点，7239个四面单元体。模型的建立及单元剖分如图2.1所示。
图2.1 土石坝的网格剖分图
该模型分析中一共分析了16种工况，分别是水平铺盖的长度是坝前水头的2～32倍（代表铺盖长度分别为16m、32m、48m、64m、80m、96m、112m、128m、144m、160m、176m、192m、208m、224m、240m、256m）的工况。
2.1.3 计算结果与分析
此计算设定在稳定渗流的条件下，不同铺盖长度计算终止时的情况。以水平铺盖长度为32倍坝前水头为例的孔隙水压力等值线云图、坝体浸润线图、总水头等值图、流速分布云图分别如图2.2～图2.5所示。
图2.2 水平铺盖长度为32倍上游水头计算终止时的孔隙水压力等值云线图（单位：kPa）
由图2.2和图2.3可以看出，在坝体的水平面以上的部分均存在大小不同的负孔压区。意味着该区域是非饱和的，而孔压为零的位置即为浸润线（面）。
由图2.6可以看出，当假设水平防渗体和斜墙防渗体的渗透系数在人为破坏的情况下，保守取1×10-7cm/s时，放大可以明显看到水流渗透防渗体的情况。证明了此类假设的可行性。
图2.3 坝体浸润线图
（a）水平铺盖长度为2倍上游水头计算终止时的浸润线图；（b）水平铺盖长度为32倍上游水头计算终止时的浸润线图
图2.4 水平铺盖长度为32倍上游水头的总水头等值图（单位：kPa）
图2.5 水平铺盖长度为32倍上游水头的流速云线图（单位：m/s）
图2.6 水平铺盖长度为32倍上游水头的防渗体放大图
通过不同的防渗体长度对饱和渗流模型浸润线位置的分布比较可以看出，在饱和-非饱和渗流中，水是可以通过浸润线流入非饱和区域的。同时可以发现，在大坝上游面，水流通过不饱和渗流线会继续流动。这些结果表明，饱和线不同于传统的渗流分析假设。
分别选取不同的铺盖长度进行模型计算，可以通过流速线得到计算终止时的流速，从而得到在此铺盖长度下模型的流量。选择监控点，可以得到坝后渗流坡降，如图2.7所示。
将所有模型该位置的流量提取并换算汇总得出的数据见表2.2。
图2.7 水平铺盖长度为32倍上游水头计算终止时的流速及各项数据
表2.2 铺盖长度与坝体坝基单宽流量、坝后渗透坡降关系
根据表2.2绘制出水平铺盖长度与单宽渗流量关系曲线图，如图2.8所示。
图2.8 铺盖长度与单宽渗流量的关系曲线
由图2.8可知，当土工膜水平铺盖长度是坝前水深的26～30倍的时候，流量减小的趋势明显变缓，此时若再通过增加水平铺盖的长度来减小渗流量，其效果已经不明显。当水平铺盖长度是坝前水深30倍时，比规范规定5～8倍水头的设计节水35%～42%。