代古寺面板堆石坝渗流场有限体积法数值分析

张明亮，宿晓辉，吴正桥，张建涛，李浩瑾

(1.大连理工大学 建设工程学部水利工程学院，辽宁 大连 116033； 2.中水北方勘测设计研究有限责任公司，天津 300222)

渗流分析和渗流控制是面板堆石坝工程中极为重要的研究内容，国内外有许多由于渗流控制不当而引起破坏或溃决的工程，1993 年青海沟后面板堆石坝因渗透破坏导致大坝溃决[1]，墨西哥的Aguamilpa[2]大坝由于面板上的水平裂缝，渗漏量达到了257.7 L/s。由此可见，面板堆石坝中由面板和防渗帷幕组成的防渗系统对坝体和坝基的安全至关重要。目前大多数渗流数值分析采用有限单元法，张凤财等[3-4]采用饱和稳定渗流有限元法研究混凝土面板坝在面板产生裂缝后的渗流场分布及渗流量的变化；潘正阳等[5]通过模拟不同深度、类型、位置、大小及分布的面板裂缝，获得了面板裂缝对大坝渗流量及浸润线等渗流要素的影响规律；何玲丽等[6]以 Richards 方程、运动波方程和有限元法为基础，提出了一种降雨时考虑径流流量补给的滑坡渗流数值模拟方法，正确反映降雨入渗对基质吸力大幅降低的作用。然而，有限元方法对计算机存储需求高，计算工作量大，应用到大型工程渗流计算时耗时严重。有限体积法与有限元法相比，计算工作量小，并且在控制体积内局部绝对守恒的特性对流动问题处理具有天然的优势。陈国芳等[7]提出了一种修正的中心型有限体积法，避免了标准有限体积方法的 “数值障碍”现象； Manzini 等[8]在非结构网格上用二阶精确单元中心有限体积法数值求解二维渗流方程，提高了数值精度；郭影等[9]提出了一种渗流吸水诱发岩体强度弱化的有限体积数值计算方法，实现了渗流作用下孔隙内自由水向基质吸附水转化过程的模拟，并实现了岩体模量与强度随渗流时间逐渐弱化过程的模拟。

本文通过自主开发的三维有限体积法渗流计算软件DUT-Seepage[10]，采用无系数矩阵数值迭代求解方法，对代古寺混凝土面板堆石坝三维复杂渗流场进行渗流特性分析，通过避免数值计算过程中的矩阵操作，减少计算工作量，并采用GPU-OpenACC 并行技术加快求解速度[10]，实现大型工程渗流模拟的快速求解。本文重点研究面板和防渗帷幕联合作用下的渗流控制效果，并开展基岩渗透敏感性分析，研究基岩渗透系数对渗流量的影响，最后探讨坝基挤压破碎带对面板堆石坝渗流场特性的影响，为其他存在地质缺陷的工程提供参考。

1 工程概况

代古寺水库河道走向呈大“S”形，河谷狭窄，呈深切“V”形，为横向或斜向河谷。两岸坝肩边坡岩体均为志留系板岩，均倾向坡内，自然边坡稳定，整体稳定性较好。坝址区滑坡、泥石流不发育。坝址区未发现区域性断层通过，构造行迹以裂隙及小规模断层为主，多倾向坡里，对坝体抗滑稳定影响不大。为了加强抗渗性，沿左右基岩向下进行帷幕灌浆，减少渗漏量。工程拟定正常蓄水位1 804 m，最大坝高151 m。

本工程地质条件复杂，坝址区岩体以物理风化为主，JI38 挤压破碎带自坝（0+115）断面往右岸延伸。坝址区风化程度不均，左、右岸强风化带厚9.4～26.3 m，弱风化带厚29.5～66.5 m。右岸岩体强卸荷带水平深度15.0～18.0 m，弱卸荷带水平深度大于35.0 m，左岸岩体强卸荷带水平深度5.0～8.0 m，弱卸荷带水平深度15.0～20.0 m。大坝最大断面如图1 所示。根据坝体填料分区原则，该面板堆石坝分区自上游至下游依次为：上游盖重区（1B）、上游铺盖区（1A）、混凝土面板（F）、垫层料区（2A）、特殊垫层区（2B）、过渡料区（3A）、主堆石区（3B）、次堆石区（3C）和下游坝基排水区（3F）。

图1 代古寺面板堆石坝最大断面（单位：m）Fig.1 The maximum cross section of Daigusi concrete face rockfill dam (unit: m)

2 渗流计算方法和原理

三维稳态渗流控制微分方程[10]为：

式中：H为总水头（m）；直角坐标轴x,y,z为渗透主方向；Kxx,Kyy,Kzz为沿渗透主方向的渗透系数。在本次工程渗流计算分析中，假定坝体各部分材料均为各向同性材料。

本文提出采用基于非结构化网格有限体积法对渗流控制方程进行离散求解，采用预处理共轭梯度法（PCG）+ GPU 并行技术对求解过程进行加速。采用Cell-Vertex（CV）方法形成控制体，并在控制体上对控制方程进行积分求解。控制体如图2 所示，控制体由围绕节点P的边中心点与相应的形心点连接而成，此控制体形成方法会在一定程度上降低网格奇异性的影响，保证计算结果的准确性。

图2 围绕节点P 的控制体示意Fig.2 Diagram of control volume around node P

采用有限体积法对渗流控制方程进行离散求解：

对方程（2）应用高斯散度定理：

式中：积分区域cv为控制节点i的控制体；积分区域S为控制体的闭合曲面；K=Kxx=Kyy=Kzz为渗透系数；H为总水头；V为控制体体积；ncell为与控制节点P相关联的单元数目；ΔSci为控制体在四面体单元i中的对应边界面。对于给定的四面体单元，可以采用闭合曲面积分公式dS=0计算nΔSci。

根据对封闭控制体所有边界面法向量积分为0，可以得到：

将式（4）代入式（3）可以得到：

式中：nΔSpi为含有控制节点P的四面体单元中P点所对应的面和法向量的乘积；(∇H)i为含有控制节点P的四面体单元的中心梯度；Ki为含有节点P的四面体单元的中心渗透系数，其四面体计算单元如图2 所示。

四面体单元的中心梯度采用格林公式计算：

四面体单元的中心渗透系数通过体积加权法计算：

式中：Hi为四面体单元顶点i处的变量值；Si为四面体单元顶点i所对应的面与法向量之积；V为四面体单元体积；Km为四面体单元顶点渗透系数；Vm为单元顶点m控制体的体积。

最终得到渗流控制方程的离散形式：

采用预处理共轭梯度法对式（8）进行迭代求解，其求解步骤为：

（1） 给定求解域初始值H0；

（2） 通过初始值H0，求得初始残差向量r0，并利用预处理矩阵M−1得到p0；

式中：M=diag(a11,a22,···,ann)

（3） 循环迭代计算，当计算收敛误差小于给定允许值时，循环迭代停止。

3 三维渗流场分析

3.1 渗流计算模型边界及参数

渗流计算模型边界类型主要包括三类：（1）水头边界，包括坝址区上游水位线以下的上游坝坡、水库岸坡和库底及坝址区下游水位线以下的下游坝坡、岸坡和河道；（2）逸出边界，包括坝体下游坡面在下游水位线以上的区域、坝址区下游水位线以上的左、右岸边坡面、模型下游侧截取边界面在下游水位线以上的区域；（3）不透水边界，包括模型上游侧截取边界面、模型左右岸两侧截取边界面和模型底部截取边界面。该面板堆石坝的三维渗流计算模型如图3 所示。

图3 三维渗流计算模型Fig.3 Three-dimension seepage calculation model

根据渗透及反滤试验，并结合地质、水文及钻孔压水实验资料，对模型中各介质的渗透系数取值如表1 所示。正常蓄水工况下，坝体上游水位1 804.00 m，下游水位1 675.00 m。

表1 面板堆石坝各材料分区渗透系数Tab.1 Permeability coefficient of each material partition of concrete face rockfill dam单位： cm/s

3.2 三维渗流场计算结果分析

图4 为渗流计算域的浸润面，图5 为沿坝轴线方向孔隙压力分布，图6 为坝防渗结构体的总水头分布。由计算结果可以看出：

图4 计算域渗流浸润面Fig.4 Seepage infiltration surface of the calculation domain

图5 沿坝轴线方向不同典型剖面压力水头分布Fig.5 Pressure head distribution of different typical cross sections along dam axis

图6 面板堆石坝不同部位的总水头分布Fig.6 Total water head distribution diagram of different partitions of concrete face rockfill dam

（1）面板坝渗流场水头分布规律合理，总水头等值线在面板、趾墙、帷幕等处较为密集，上游水头由这些部位承担，渗流控制系统具备很好的防渗效果。

（2）在混凝土面板堆石坝内，面板上下游侧水位下降最大，上游水流经过混凝土面板后，总水头由上游正常蓄水位1 804.00 m 降到1 680.50 m 左右，经垫层区、砂砾堆石区至排水区时，总水头降至 1 675.15 m，经过混凝土面板后，坝体最大断面逸出点水力坡降为0.008，小于允许水力坡降0.1，满足正常运行要求，其余各填筑分区水力坡降均小于允许水力坡降。

（3）通过计算坝轴线断面的流量得到最终坝基渗流及绕坝渗流量，其中通过坝基和帷幕的渗流量为586 m3/d，两岸绕坝渗流量为82 m3/d。计算结果表明左、右岸渗漏量很小，绕坝渗流并不明显。

（4）通过GPU-OpenACC 并行技术[10]，对本案例千万级计算单元网格的渗流计算仅需3～5 h，能够满足工程渗流时间要求。

4 基岩渗透系数对渗流量的影响

由于面板堆石坝基岩实际渗透系数难以确定，为了分析基岩渗透系数的变化对面板堆石坝渗流量的影响，本文另外选取了2 组不同的基岩渗透系数，计算得到不同基岩渗透系数对应的渗流量，计算结果见表2。

表2 正常蓄水位下不同基岩渗透系数坝体渗流量Tab.2 Seepage flow with different bedrock permeability coefficients under normal storage level

从表2 可知，随着基岩渗透系数的增大，左、右岸及坝基渗透量也呈现逐渐增大的趋势，尤其是当微风化基岩渗透系数增大时，坝基渗透量出现急剧增大现象。

通过与相近坝高的国内外面板堆石坝工程进行渗流量[11-17]比对发现，该面板堆石坝的渗流量处于中等偏低范围，总体偏于安全，对比结果如图7 所示。

图7 相近工程渗流量对比Fig.7 Comparison of engineering seepage flow

5 JI38 挤压破碎带影响分析

本工程基岩中存在JI38 挤压破碎带，宽度自左岸向右岸逐渐加宽，坝（0+285）断面的破碎带宽度较大，约11 m。由于JI38 挤压破碎带的渗透系数未知，故选取3 组不同的渗透系数分析JI38 挤压破碎带对浸润线位置的影响，并将计算结果与GEO-Studio 计算结果进行对比。

建立面板堆石坝（0+165）二维渗流分析模型，针对坝（0+165）断面处存在的JI38 挤压破碎地质缺陷，在正常蓄水位1 804.00 m、尾水位1 675.00 m 工况下，开展JI38 破碎带渗透敏感性分析。JI38 挤压破碎带的渗透系数分别取1.00×10−1、1.00×10−4和1.00×10−5cm/s，坝体其余部分渗透系数参照表1。

分别采用DUT-Seepage 与GEO-Studio 程序，分析不同工况下JI38 破碎带渗透系数对压力水头和总水头分布的影响，计算结果见图8～10。对比可发现DUT-Seepage 与GEO-Studio 的计算结果较一致；不同渗透系数下的大坝浸润线位置几乎不变，但当渗透系数为1.00×10−1cm/s 时，坝基的压力水头及总水头分布产生较大变化，并且JI38 破碎带周边的水力坡降从0.40 增至0.75，增大了坝基渗透破坏风险，建议对 JI38 破碎带进行加固处理，降低渗透系数，避免坝基发生渗透破坏。

图8 不同工况下JI38 破碎带渗透系数对孔隙水压力分布的影响Fig.8 The effect of permeability coefficient of JI38 crush zone on pore water pressure distribution under different cases

图10 不同工况下JI38 破碎带渗透系数对水力梯度的影响Fig.10 The effect of permeability coefficient of JI38 crush zone on hydraulic gradient distribution under different cases

6 结 语

本文利用三维有限体积法渗流计算分析软件DUT-Seepage，采用无矩阵数值算法及GPU-OpenACC 并行技术对代古寺混凝土面板堆石坝进行了三维渗流数值模拟分析。三维渗流结果显示代古寺面板堆石坝工程的防渗系统，包括面板与防渗帷幕有效阻断了上游水流的渗透，具有较好的防渗效果，堆石坝主体无渗透破坏的风险；当基岩渗透性较大时，大坝整体渗漏量约为2 100 m3/d，处于同类工程的偏低水平；若挤压破碎带渗透系数过大时，存在地基渗透破坏风险，建议对坝基破碎带进行加固处理。