2.7 心墙堆石坝渗流计算理论与方法

2.7.1 渗流分析发展概况

渗流力学是流体力学的一个分支，它是多种科学和工程技术的理论基础。1856年，法国工程师达西 （Henri Darcy）通过试验提出了线性渗透定律，为渗流理论的发展奠定了基础。1889年，H.E.茹可夫斯基 （H.E.Жyкoвcкий）首先推导了渗流的微分方程。此后，许多数学家和地下水动力学科学工作者对渗流数学模型及其解析解法进行了广泛和深入的研究，并取得了一系列研究成果。但解析解毕竟仅适用于均质渗透介质和简单边界条件，在实用上受到很大限制。

随着电子计算机的迅速发展，数值方法（如有限差分法、有限单元法和边界元法等）在渗流分析中应用越来越广泛的。有限差分法是1910年由理查森（L.F.Richardson）首先提出的，经过长期的研究和广泛应用，目前该方法已具有较完善的理论基础和实用经验，有限单元法的基本思想早在1913年由柯朗（R.Courant）提出，1960年克劳夫（R.W.Clough）最先采用“有限单元法”这个名称，以与有限差分法相区别。1965年，津克维茨（O.C.Zionkiewiz）和张（Y.K.Cheung）提出有限单元法适用于所有可按变分原理进行计算的场问题，为该方法在渗流分析中的应用提供了理论基础。随后该方法在渗流分析中逐步推广应用。边界元法建立在经典力学理论基础上，贝蒂（Betti）互换定理及弗雷霍姆（Fredholm）积分方程早在18世纪末及19世纪初就已提出，建立在这些理论基础上的边界元法初见于20世纪60年代后期，当时被称为边界积分方程法（Boundary Integral Equation Method）。直到1978年，边界元法（Boundary Element Method）这个名称才被确立并得到公认。后来许多学者针对所研究问题的不同特点，研究和提出了能集合上述各数值方法优点的杂交元法（Hybrid Element Method），更合理地解决实际的工程渗流问题。

国内的渗流理论及方法的引入发展，均与国家经济建设紧密相连。20世纪50年代毛泽东提出“一定要把淮河治好”，当时南京水利科学研究所（现南京水利科学研究院，以下简称南科院）开始在苏联专家帮助下，引进介绍渗流学科，并全面开展学习、建设研究。南科院以毛昶熙为首的渗流组设计建设了水电模拟试验槽等渗透试验设备并为水利水电工程的渗流控制问题进行了多项试验研究，为水利水电行业做出很大贡献。水利水电行业从“六五”“七五”对100m级土石坝渗流控制问题已开始关注，主要是结合国家科技攻关课题开展土石坝堆石面板坝等关键技术问题的应用研究，为国家水电工程提供了技术支持。在“八五”“九五”“十五”国家科技攻关课题中，更进一步关注了高坝、地下洞室设计建设中遇到的各种岩体渗流问题和高坝岩基及地下洞室的渗流控制问题。从此将纯多孔介质渗流控制的理论观点拓展到裂隙岩体介质的渗流控制研究领域。80年代原北京水利科学院（现中国水利水电科学研究院，以下简称北科院）岩土所以杜延龄、许国安为首的渗流室建设了大型电网络模型，与数值计算并肩前行相互印证，用以研究岩体裂隙各向异性渗流问题，为解决当时一些水电行业岩体工程渗流的问题提供了技术支持。北科院的刘杰在土体渗透变形和反滤层设计方面为工程做了很多有益的贡献。北科院结构所张有天引进创新了岩体渗流数值方法为地下洞室工程渗流控制问题提出了宝贵意见。80年代初河海大学渗流实验室建成了国内最大的三维渗流水电模拟自动测试系统，配置了渗透变形、土工织物等渗透测试设备，为土石坝、尾矿坝、灰坝的渗流控制提供了技术支持。长江水利委员会长江科学院土工研究所和黄河水利委员会水科所（现黄河水利科学研究院）也先后建设了电网络模型和水电模拟试验槽，主要服务于长江三峡工程及长江流域和黄河流域水利水电工程。

国内水电行业开展数值计算起步比国外晚十余年，南科院、黄河水利委员会水科所、北科院等研究机构在20世纪70年代，率先学习引进国外渗流有限元数值计算方法。接下来的10年（20世纪70—80年代）是数值计算方法推广发展的时期，这个时期电模拟技术和有限元数值计算并存，在解决工程渗流问题上两种方法同台共舞。在这期间南科院渗流室同行们，为推广、普及渗流有限单元法做了很多有益教育宣传工作。渗流界（包括科学院所、高等院校、设计工程界）定期联合召开全国水利水电工程渗流学术会议，帮助更多的人去关注学习渗流力学，渗流对水工建筑物安全影响的重要性认识大大提高。90年代，是水电行业开始快速发展的时期，也是渗流数值计算发展完善和进一步推广的重要时期。这时，解决工程渗流问题已开始逐渐依靠初步完善的三维渗流有限元计算方法。

在岩体渗流问题研究方面，国外起步于20世纪50年代，国内起步于70年代，由于各种原因国内进展十分缓慢。90年代开始由于国家水电建设大发展需要解决高坝岩基地下洞室中大量岩体渗流问题，各科学院、高校于80年代末、90年代初连续申请有关岩体渗流问题的国家自然科学基金、省部级开发基金、博士点基金。河海大学在国内率先建设了高压（31.5MPa）多功能渗透仪及渗流应力耦合测试系统、岩体渗流模型试验系统。在此基础上又积极参加了国家科技攻关课题中关于岩体渗流特性、运动基本规律、渗流应力耦合、泄洪雾化雨在岩体岸坡中入渗的非饱和渗流等基础理论问题开展了大量物理模型试验和计算分析研究，得到许多有益成果和启示，发表了一系列论文。河海大学渗流实验室针对国家重大工程，如长江三峡船闸高边坡等工程中的渗流场及渗流控制分析的需要，首次以变分不等式原理为理论基础，研究提出了新的不变网格渗流有限元计算方法，将生产建设的需求和理论紧密结合，以后又紧密结合小湾、溪洛渡、龙滩、锦屏、糯扎渡等重大工程的需求逐步完善开发出大型专业的稳定-非稳定渗流、饱和-非饱和及渗流应力耦合有限元三维计算程序。为了配合设计快速发展的需求，又自主开发了密集排水孔的解析法与渗流有限元相互耦合的快速算法，以满足工程设计中防渗排水多方案优化布置的需要，为复杂岩基上高坝、地下多条洞室的复杂水工枢纽的整体渗流分析及整体防渗排水的优化布置提供了有效的技术支持，取得良好的工程效果。

20世纪50年代以来，越来越多的量测信息可资利用，这成为反分析发展的一个触发点。自1960年Philip提出了一种求解一维非线性渗流反问题的精确解法，国内外众多学者对岩土工程渗流反问题进行了大量的研究。内容涉及各类参数的计算方法、优化方法及适用性研究。我国关于位移反分析的研究开始于20世纪70年代末，基本上与国际同步。冯康院士早在80年代初就大力提倡开展反问题数值解法的研究。此后，反分析方法被拓广应用于岩土工程的各个方面，并发展出各类计算方法。近年来，随着科学技术水平的提高，各种交叉学科、前沿学科不断被引入到岩土力学反分析中，推动了反分析的发展。随着人工智能技术的引入，渗流反分析出现了智能化的趋势。目前，各种渗流反演分析基本方法可概括为以下6种方法：试错法、解析法、脉冲谱法、数值优化算法、随机反演法、人工神经网络法。其中，人工神经网络法自20世纪80年代中后期以来，迅速发展为一个前沿研究领域，被广泛用于各个学科领域。目前，国内外基于人工神经网络的参数反分析研究成果多处可见。与传统的方法相比，若有可能覆盖整个计算域的充足的训练样本对，神经网络反演结果更接近于真实值，而且人工神经网络具有较强的非线性动态处理能力，无需知道状态变量和系统参数之间的关系，可实现相同或不同维数向量之间的高度非线性映射。目前应用较多的是BP网络，在神经网络中需要大量的历史及现存的实测资料对网络进行训练，且需要较长的训练时间，这是该方法在实际应用中的一大障碍。另外神经网络方法还易陷入局部最优，无法达到全局最优值。可采用改进的遗传算法、模拟退火法、牛顿-高斯全局优化算法、交替迭代全局优化算法等和神经网络相结合来寻找全局最优点。

2.7.2 高土石坝坝体坝基系统渗流控制分析计算理论及方法研究

2.7.2.1 渗流有限元基本方程及有限元分析方法

渗流计算是在已知定解条件下求解渗流微分方程，以求得渗流场水头分布、渗流量、渗透梯度等渗流要素，它是工程设计的重要依据。由于无压渗流有渗流自由面（浸润线），且非稳定渗流自由面随库水位升降而变动，加之一般渗流场有不同程度的非均质和各向异性，几何形状和边界条件较复杂，解析求解在数学上存在不少困难，仅能对一些简单流动情况获得解析解。20世纪60年代，电子计算机的普及和数值计算方法的发展，特别是有限单元法的推广应用，促进了渗流数值模型的发展，为渗流计算提供了有效的方法。

有限单元法是把研究区域离散化为有限个单元体的集合来进行研究的。引用变分原理或伽辽金法，针对研究问题建立模型，推导建立求近似解的线性方程组，从而求解出渗流场的分布等渗流要素，供设计使用。

1.渗流有限元基本方程

众所周知，对于饱和稳定-非稳定渗流基本方程为

式中 (x,y,z)——渗透主方向;

H=H(x,y,z,t)——待求水头函数；

K_x、K_y、K_z——主向渗透系数；

S_s——单位储水量或储存率。

当在稳定，渗流状态是均质的条件下，式（2.7-1）退化为

式（2.7-2）为饱和稳定渗流的基本方程。

2.渗流有限元分析基本方法

取8节点等参单元离散渗流场，则单元内的水头分布为

应用Galerkin方法将式（2.7-2）离散，经推导整理得到

显然，当整个区域全部处于承压状态（如混凝土坝坝基渗流），便可直接依式（2.7-4）建立代数方程组，进行求解即得到所要求的渗流水头场。然而，对于具有自由面的渗流问题，其实际渗流区域，往往小于整个渗透介质的区域。由于自由水面正是渗流分析所需要求解的问题，因而，实际渗流域是未知的。这一点便是渗流分析较固体力学问题求解复杂的根本点，从而也决定了具有自由面渗流场不可能一次性直接解出，而必须反复迭代计算获得逼近于真解的数值解。另外，渗流计算问题是典型的边值问题，边界条件对计算结果影响极大。为了得到较准确的渗流场，计算中尽量选用明确可靠的边界条件，如河流边界、分水岭等。这也导致渗流场计算区域比结构计算截取的范围大得多，使计算的前期准备工作量也大大增加。

具有自由面渗流场的数值分析计算，由于在计算中自由面边界是未知的，使得渗流计算变得复杂和困难，通常用迭代逼近的方法来求其近似解。河海大学渗流实验室研究了改进的初流量法和基于不等式求解渗流自由面，并开发了相应的数值程序。

2.7.2.2 裂隙岩体渗流应力耦合分析基本方法

岩体经过长期地质作用一般都发育有裂隙、孔隙和溶隙等不连续结构面和断层等空隙。这些空隙是地下水赋存场所和运移通道，其分布形状、大小、连通性以及空隙的类型都将影响岩体的力学性质和岩体的渗流特性。裂隙是岩体渗流的主要通道，它对岩体的水力行为起控制作用。高坝必然伴随高水头的作用，坝体本身渗透性可人为通过分区设计并加强施工监管进行控制，但基岩断层和裂隙渗透性是天然形成的，往往不以人的意志为转移，在土石坝失事案例中，有相当比例是由于对地质条件认识不够，渗流通道处理不当造成的。因此人们对岩体裂隙的渗流特性给予了很大的关注。

岩体渗流有以下特点：①岩体渗透性大小取决于岩体中结构面的性质及岩块的岩性；②岩体渗流以裂隙导水、微裂隙和岩石孔隙储水为其特色；③岩体裂隙网络渗流具有定向性；④岩体的渗流一般看作非连续介质（对密集裂隙可看作等效连续介质）；⑤岩体的渗流具有高度的非均质性和各向异性；⑥一般裂隙发育岩体中的渗流符合达西线性定律；⑦岩体渗流受应力场影响明显；⑧复杂裂隙系统中的渗流，在裂隙交叉处具有“偏流效应”，即裂隙水流经大小不等的裂隙交叉处时，水流偏向宽大裂隙一侧流动。

1.裂隙岩体渗流分析的等效连续介质模型

(1)达西定律。认为岩体中水流运动为层流，仍服从线性达西定律。

式中 v_i——流速分量;

k_ij——渗透张量；

J——水力坡降；

h——水头。

(2)等效渗透张量k_ij。根据Snow的理论有：假定岩体内的裂隙可以概化为N组，其中第k组裂隙的单位法向量为，平均张开度为b_k。当各组裂隙均为无限延伸时，Snow以立方定律为基础给出等效连续介质的渗透张量为

式中 ρ_k——第k组裂隙的平均体密度；

——该组裂隙的平均等效水力隙宽。

式中 m_(k)——该组裂隙总条数;

——第s条裂隙的等效水力隙宽。

实际裂隙并非无限延伸，为考虑非贯通裂隙，假设裂隙为圆盘形，其平均半径为r_（k）。以式（2.7-6）为基础，渗透张量近似表示为

式中 g——重力加速度；

ν——运动黏性系数；

r——裂隙圆盘半径；

λ——反映裂隙连通性及不同组裂隙相交次数的一个无量纲系数，0≤λ≤1/12，当裂隙无限贯通时，λ=1/12；

ρ——第k组裂隙的体密度，ρ=m_(k)/V。

2.基于等效连续介质的裂隙岩体渗流-应力耦合模型

（1）稳定渗流基本微分方程。将式（2.7-5）代入渗流连续方程v_i,i=0可得：

式中 w——恒定降雨入渗或蒸发量。

(2)根据有效应力原理，裂隙岩体的力学平衡方程为

式中 σ_ij——有效应力;

X_i——体积力;

p=γ_w(h-z)——渗流水压力。

(3)裂隙岩体渗透张量与应力的耦合关系。

裂隙的弹塑性本构关系

裂隙面上应力增量与岩体平均应力增量的关系

式中 ［F_ij］′_k——裂隙的应力集中张量;

——REV上的平均应力增量在裂隙局部坐标系中的投影。

Barton(1985)通过大量的试验，得到机械隙宽b_m、等效水力隙宽b_h与JRC之间的经验关系式

其中，b_h、b_m的单位为μm。

裂隙岩体渗透性与应力的关系是通过裂隙隙宽在应力作用下的改变来实现的。因此下面需讨论一下裂隙等效水力隙宽b_h的计算。当单元的应力{σ_ij}计算出来以后，将其转换到裂隙面局部坐标系上得出，再由式 （2.7-11）计算裂隙面上的应力状态}，然后按式 （2.7-10）计算裂隙机械隙宽的改变量，最后根据机械隙宽与水力隙宽之间的关系式 （2.7-12）求出裂隙等效水力隙宽的大小。最后由式 （2.7-6）或式 （2.7-7）计算岩体的渗透张量。

2.7.2.3 渗流场反演分析基本方法

渗流反演分析法属大坝安全评价范畴，该方法利用高心墙堆石坝渗流场中各测点水位、渗压等的实测值与计算值的最优化拟合准则，开展工程渗流场的反演和反馈分析，以便对工程渗透安全作出正确评价和提出进一步保证工程安全的措施。

土石坝渗流参数反演分析方法，是现代规划理论、数值分析和观测技术的综合运用。它为较为精确地确定土石坝的土体参数提供了有效的手段，因而得到了越来越多的设计和研究人员的重视并开始了这个领域的研究工作。目前发展起来的众多反演方法有其各自的使用范围，而且由于反问题固有的不适定性以及实测资料的多样性，使得不同的反演方法在解决同一问题时应用效果大相径庭。为此，有必要在充分了解各类反演方法的基础上，认真选择适合具体问题的方法，以得到更加可靠的反演结果。

根据渗流观测资料，开展渗流场反演拟合分析，原先的方法之一是把目标函数在设计点处作线性化，但是，求解时往往会出现设计点跳动，收敛较困难。目前改进的主要方法有：①把目标函数序列二次化，引入序列二次规划方法来进行渗流场的反演；②有选择地确定了反演变量类型；③提出更方便的水头函数对反演变量的一阶、二阶导数计算方法，使问题的求解效果好，速度快。

下面简要介绍离散-优化方法的主要求解思路。

1.目标函数

反求渗流参数一般是根据研究的渗流区域内若干已知坐标位置的水头观测值与计算值之间的误差，应用最小二乘法建立目标函数。

式中 K_ij——第i区中的第j个渗流参数；

ω_k——第k个观测点的权函数；

——第k个观测点的水头计算值；

——第k个观测点的水头观测值。

数学上已经证明式（2.7-13）必存在极小值。显然，求解的目标是寻求一组渗流参数，使得式（2.7-13）定义的目标函数值达到最小，这是需要直接解决的问题。应当指出，由式（2.7-13）定义的目标函数中的待求解变量，即地层及断层的渗透参数，还必须满足相应的工程意义。因此，提出以下约束条件。

2.约束条件

反分析渗流参数应在已知渗流参数附近。

这里，α_i为反分析渗流参数与已知的渗流参数的倍比，即反分析渗流参数最多只在已知的参数基础上缩小5被或扩大5倍。

3.反分析变量的正交搜索法

如前所述，由于实际工程问题反分析的极端复杂性，只能采用离散-优化方法来求解式（2.7-13）目标函数在约束式（2.7-14）条件下的极小值问题。尽管如此，原则上可采用诸如单纯形或最速下降等方法求解，但由于偏导数的计算，不仅同样增加正分析计算量，尤其是实际问题的自身复杂性，包括大型数值问题计算自身产生的误差等因素，常常使得偏导数不再单调，产生伪多值情况。为保证反分析的有效性和正确性，提出正交搜索法。实际计算表明，尽管此方法非常原始，但却非常有效。其基本思路是按约束区间作等距划分，从第一个变量开始搜索，根据目标函数值的大小，搜索出第一个变量的相对最优解。在第一个变量最优解搜索完成后，即固定该变量。如此依次对每个变量单独搜索，犹如正交模型实验一样，故称之为正交搜索法。通过与单纯形法比较，不仅最终优化结果相近，而且分析次数减少20%左右。

2.7.2.4 SPGCR-3.FOR计算分析程序介绍

根据上述理论，河海大学渗流实验室编制了大型的三维渗流及渗流应力耦合分析程序SPGCR-3.FOR。该程序应用在大量实际工程中并不断修改完善，用以计算复杂稳定、非稳定三维渗流场以及渗流应力耦合计算，使用效果良好。其功能有：①模拟水利工程的天然渗流场和建坝后的渗流场，进一步对水电站坝体、坝基及地下厂房、洞室群进行三维渗流控制优化分析研究（例如小湾电站坝体、坝基及地下厂房洞室群三维渗流控制优化分析研究）。②在渗流场精细模拟方面，根据地下厂房洞室及主要通道布置，可以模拟地形、地质条件，考虑不同围岩的渗流特性，建立地下厂区的精细模型（网格密度要充分考虑渗流场的水力梯度），还可以进一步用于提出地下厂房防渗排水系统的优化布置方案（例如溪洛渡水电站地下厂区三维精细模型渗流及防排系统分析研究）。③排水孔因尺寸小、排列密集、数量众多导致模拟困难。该程序可以求解密集排水孔的渗控影响问题，它能将排水孔的空间位置、走向及其边界性质的有关集合参量加以描述，从而更适用于许多大型水电工程的复杂渗控系统的优化分析研究。图2.7-1为SPGCR-3.FOR流程图。

2.7.3 渗流计算分析工程实例

工程计算实例对糯扎渡工程区建坝后的渗流场模拟分析，主要包括设计初拟渗控情况、防渗排水系统多组合的渗控优化分析、最大剖面坝段精细模型的计算分析以及溢洪道消力池排水布置渗控计算分析等。以下简单介绍糯扎渡工程区建坝后的整体三维网络模型和部分主要部位的网络模型，以及以大坝精细模型的计算分析为例简略介绍，并给出此例工程计算重要结论的一部分。

2.7.3.1 糯扎渡高心墙坝基本资料

糯扎渡水电站三维渗控分析研究。糯扎渡水电站坝址区地质条件复杂，河床部位基岩面绝大部分为弱风化下部及微风化-新鲜的花岗岩。左岸坝基分布的基岩大部分为花岗岩，高程约815.00m以上为忙怀组沉积岩，左岸节理与断层发育，以F₉、F₁₅断层规模较大。右岸冲沟发育，顺水流方向地形起伏较大，山坡不平整。右岸坝基分布的基岩均为花岗岩，断层发育，岩体蚀变较为强烈，影响较大的有F₁₂、F₁₃、F₁₄、F₅、F₁₁及F₁₆等。在本次渗流控制分析建模工作中，考虑三大地层和12条主要断层。

糯扎渡工程区建坝后的渗控计算分析，经三维有限元离散，共剖分单元总数为32176个，节点总数为30881个。断层、帷幕、坝体、心墙、高塑性黏土、混凝土垫板及防渗排水结构等特殊构造均考虑了过渡、衔接与连续性要求。整体网格模型和各主要部分的网格模型见图2.7-2～图2.7-4。

2.7.3.2 糯扎渡高心墙坝精细模型渗流场

在全面考虑工程细部结构，如坝体心墙、反滤层、过渡层、堆石体、地基各地层、防渗帷幕、固结灌浆、混凝土垫板等的基础上，经过细致的有限元网格划分，最终形成了宽32m、单元10140个、节点12254个的精细网格模型，见图2.7-5。

糯扎渡高心墙坝精细模型坝段各材料渗流梯度见表2.7-1。

糯扎渡高心墙坝精细模型渗流量计算：通过对该坝段（32m）的等效节点渗流量计算分析，得到其坝体坝基渗流量为Q=395m³/d，单宽渗流量为12.35m³/d。和总体模型比较知，最大断面单宽渗流量约为总平均单宽渗流量的3～4倍。计算分析表明，对于糯扎渡高达261m的高坝，在设计渗控布置条件下的渗流量尚较小，满足工程控制要求。

2.7.3.3 糯扎渡工程渗流分析主要成果

应用SPGCR-3.FOR程序针对糯扎渡工程区建坝后设计初拟渗控情况下的渗流分析，得出坝体坝基和边坡的渗流规律，并提出相应的建议。

通过针对糯扎渡工程区防渗排水系统多组合的渗控优化分析，其基本结论如下：

(1)坝基防渗帷幕渗透性。以帷幕渗透系数0.5Lu时的水头势等值线分布最稀疏。单从地下水位控制来看，要求帷幕透水率达到0.5Lu也是未尝不可，为达到这个目标，应该对灌浆材料及灌浆技术做进一步研究。

(2)关于坝基防渗帷幕排数与深度，当时课题组的结论是对于糯扎渡心墙土石坝设计初拟的二排帷幕及其深度是合适的。

(3)两岸防渗帷幕平面延伸长度：①建议帷幕布置按设计初拟方案，即右坝头及右岸帷幕布置为一排，深度为深入微透水岩体；平面长度为设计拟定长度100m。溢洪道左侧帷幕长度83m，斜向山里偏下游转折帷幕120m，防渗帷幕深度为深入微透水岩体。②当两岸按渗控组合表中适当延长后渗流场改变很小，而且渗流量减小均在10%以内。故此，不必延长和加深。

(4)心墙混凝土垫板开裂的敏感性：①当心墙垫板在顺河方向发生开裂最危险。②顺坝轴线局部开裂的具体特性，对工程影响最密切的是高塑性黏土，为此，专门计算分析了此时高塑性黏土的渗流梯度。从关于顺坝轴线方向垫板开裂后其顶面高塑性黏土的渗流梯度分析看，如高塑性黏土可能出现渗透变形，则应当注意采取保护措施。

(5)通过库水位骤降非稳定渗流计算分析，大坝坝基防渗帷幕渗流梯度约为13.0，高塑性黏土渗流梯度约为11.0。对于岸坡区域，渗流梯度均较小，断层的平均渗流梯度也较小约0.3，满足渗控要求。

（6）糯扎渡最大剖面坝段精细模型渗控计算分析表明，在坝体坝基防渗心墙及防渗帷幕作用下，其坝体防渗心墙下游侧相对较低，消刹水头180m左右，坝体中的地下水位总体较低。同时，从渗流场分布可知，防渗系统的效果明显，验证设计布置方案可满足其控制要求。通过对该坝段（32m）的等效节点渗流量计算分析，得到其坝体坝基渗流量为Q=395m³/d，单宽渗流量为12.35m³/d。