基于光滑有限单元法的土石坝无压渗流场数值模拟

2021-12-22戴前伟孔重阳韩行进

水资源与水工程学报 2021年5期

戴前伟，孔重阳，雷轶，张彬，韩行进

(1.中南大学地球科学与信息物理学院，湖南长沙 410012； 2.中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙 410012； 3.五凌电力有限公司，湖南长沙 410004)

1 研究背景

水利工程中水坝的修建能调节控制下泄水量，从而降低洪涝等灾害带来的影响，并且提高水资源的利用率，满足人们生产和生活的需要。无压渗流问题在岩土工程、水利工程、地下水等领域广泛存在，渗流或者孔隙压力过大会导致水坝失稳[1]。土石坝作为最常见的一种水坝类型，其坝体具有一定的透水性，在水的渗流作用下容易发生渗漏而遭到破坏[2-4]，因此需要对土石坝的渗流安全问题进行更深入的研究。

目前，坝体的工程设计、施工、检测等常常要求探测出自由面的位置，从而对坝体的渗漏情况进行初步判断并提前做好相应的防渗措施，然而自由面的位置和形状是事先未知的，并且在求解过程中会不断发生变化[5]，如何能高效准确地确定自由面以及渗流溢出点的位置正是渗流分析中需要解决的关键问题[6]。渗流分析中最常用的方法是数值模拟，不同的数值模拟方法已经被国内外的水文地质学者们应用于相应的研究之中[7]。Neumann[8]率先利用固定网格的有限元法对渗流自由面问题进行了求解，其主要特点是在求解渗流自由面的迭代过程中不对原始单元网格进行修改和调整。付延玲等[9]将整个求解域的网格单元划分为3种类型，在处理过渡单元的渗透矩阵时，引入了双参数罚函数，利用流量等效法对渗透系数进行调整，提高了计算过程的效率以及结果的精度。唐红影[10]对弃单元法进行改进，使计算域逐渐靠近实际渗流区域，并对复合单元进行二次细分，降低了网格单元对计算精度的影响。Darbandi等[11]提出的有限体积法确保了单元的质量守恒，但要求网格与自由面相匹配。Chugh等[12]的边界元法只对域边界进行离散，通过移动边界节点可以简单地改变域的几何结构，虽然可以降低问题的维数，简化求解，但对于非均匀区域和非线性问题，边界元法很难求解。Zheng等[13]将数值流形方法与无网格Galerkin 方法相结合，对非均质坝体的自由面进行了数值模拟，通过增加积分点支撑域的场节点，计算精度得到了很大提高。Zhang 等[14]用移动Kriging 无网格法和蒙特卡罗积分法分析了无约束渗流问题，但在迭代过程中会移动场节点，这给有限元分析中自适应网格法的使用带来了一定的困难。Dai等[15]提出了一种新的综合方法，将自适应网格策略与Galerkin有限元法相结合，对自由曲面的位置进行搜索，取得了很好的效果。

本文将光滑有限元法应用于无压渗流问题的分析中，建立二维无压渗流问题的光滑有限元模型，在饱和渗流区域计算求解自由面的位置。通过对矩形均质坝和直角梯形均质坝经典算例的计算，并与其他数值方法的结果进行对比来验证本文方法的有效性。同时模拟了土石坝中不同的渗透异常体对渗流速度、流体压力、水头分布等渗流参数的具体影响，为之后的渗流反演提供相应的依据。

2 无压渗流问题的数学描述

为了方便描述无压渗流问题，以图1为例示意土石坝中的无压渗流。图1中的问题域被自由面AE分为两部分，自由面以上的部分ΩU为不饱和域，以下的部分ΩS为饱和域。假设问题域中渗流仅仅流经饱和区域。无压渗流分析中最重要的问题即为自由面的求解，只要寻找到自由面的位置，饱和域与不饱和域便可区分开来。

图1 土石坝坝体无压渗流示意图

在坝体饱和域中，任意一点的水头值的数学表达式为：

(1)

式中：φ为渗流测压水头，m；p为流体压力，Pa；ρ为流体密度，kg/m3；g为重力加速度，m/s2；y为竖直方向坐标，m。

饱和域中流体连续性方程为：

(2)

根据达西定律，饱和域中的渗流速度v为：

(3)

式中：vx、vy分别为水平向、竖直向的渗流速度，m/s；kx、ky分别为水平向、竖直向的渗透系数，m/s。

将公式(3)代入公式 (2)中，可得到二维无压稳定渗流的微分方程[16]：

(4)

在图1中，相应边界条件如下：

(1) 上下游库水AB、尾水CD水头边界条件：

φAB=H1

(5)

φCD=H2

(6)

(2) 底部边界BC为不透水边界，其外法向方向流量qn为0，即满足流量边界条件：

qn=v·n=0

(7)

式中：n为单位外法向向量。

(3) 自由面AE需同时满足两个边界条件：

qn=0,φ=y

(8)

(4) 渗流面DE边界条件：

qn≤0,φ=y

(9)

即在渗流面上有流量流出。

3 无压渗流场的光滑有限元模型

本部分简要介绍应用在无压渗流问题中的光滑有限元法。利用固定网格将整个问题域进行离散化，每个节点的水头值可以通过插值得到：

uh(x)=N·H

(10)

式中：N为形函数矩阵；H为网格节点的水头值。

求解无压渗流问题控制方程的过程中，需要得到水头梯度，而该值的获得需要计算出形函数梯度。在经典有限元法中，形函数梯度矩阵由形函数求导得到，而在光滑有限元法中，利用梯度光滑技术得到形函数梯度矩阵：

(11)

本文中采用的光滑函数为Heaviside-type型函数：

(12)

利用高斯散度定理(或者是分部积分法)，将面积分转化为线积分,并将公式(12)代入公式(11)得：

(13)

由公式(13)可知，光滑梯度技术将原先二维的面积分转换成沿着边界进行的线积分。将公式(13)代入节点水头公式(10)中得到：

(14)

(15)

其中：

(16)

将控制方程和边界条件改写为积分弱形式，并将水头插值代入得到无压渗流问题的Garlerkin弱形式：

(17)

(18)

弱形式中，R与自然边界条件相关，由前文边界条件可知边界沿外法向方向流量为0，即R=0[17]。

利用固定网格法求解无压渗流问题，可以降低对网格的依赖性，并且不需要在迭代过程中进行网格的修改。由于网格是固定的，每次迭代过程中得到的自由面会将整个求解域分为上下两部分，产生3种类型的光滑单元，如图2所示，自由面以上的单元称为外部光滑单元，以下的称为内部光滑单元，自由面经过的单元称为相交光滑单元。由于采用的是四边形单元，自由面经过相交光滑单元时会将该单元切分成不同的形状，可能出现的切分形状如图3所示。

图2 3种类型的光滑单元示意图

图3 边界相交单元被自由面可能切分的形状

(19)

式中：Ω1为内部光滑单元区域；Ω2为相交光滑单元区域。

在每个单元上运用梯度光滑技术得到单元的水头值梯度，公式(19)可改写为单元积分和的形式：

(20)

在计算过程中，认为每个网格中的渗透系数为常数，而且光滑形函数梯度矩阵也为常数矩阵，则总体刚度矩阵为：

(21)

4 迭代求解算法

利用光滑有限元法求解无压渗流问题，根本目的是寻找自由面的位置，根据自由面满足的边界条件，不断在迭代过程中更新自由面位置和形状，直至达到一定的收敛条件。具体的迭代过程如图4所示。

图4 自由面迭代求解过程示意图

(1) 选择1个线段作为初始自由面，如图4中AD。自由面的溢出点需要在CE上，而初始点位于BF上。

(2) 在CE上选择一点D作为自由面溢出点的初始解。

(3) 连接点A与D，形成的线段AD作为自由面的初始解。

(4) 将初始自由面AD切分成一系列控制点H1，H2，…，Hn。选择如图4所示的BC线段做为基线，并将基线BC划分为一系列基点G1，G2，…，Gn。

(5) 分别连接H1点与G1点，H2点与G2点，…，Hn点与Gn点，得到一系列射线H1G1、H2G2、…、HnGn。在迭代过程中，控制点H1，H2，…，Hn组成了自由面，并且这些控制点会沿着相应的射线移动，形成新的自由面。

上述准备工作完成后，利用假设的初始自由面AD解上述方程，可以得到求解域各个节点的水头值。在这一步骤中，控制点的位置必须沿着射线进行修改以满足自由面上的边界条件。控制点的垂直坐标表示为：

yinew=yiold+λσi

(22)

σi=hi-yiold

(23)

式中：yiold和yinew分别为第i个控制点在前一次和下一次迭代的高程值；hi为在控制点处计算得到的测压水头值；λ为给定的步长，本文取λ=0.8。若步长过小，则迭代次数增多，计算消耗较大；若步长过大，则难以到达收敛状态。在能达到收敛状态的情况下，选择适合的步长来减少迭代次数，降低计算消耗。每次进行迭代前进行如下判断：

(24)

式中：n为控制点个数；ε为非常小的参数，本文中取值为ε=10-3，该值越小，说明计算结果中控制点水头值与高程值越接近，但迭代次数会增加、计算消耗增大，也有可能出现迭代难以收敛的问题；该值越高，则计算结果的误差越大，难以满足计算精度的要求。如果判断式(24)得以满足，则迭代过程终止，否则根据公式(22)、(23)对控制点位置进行更新。

5 算例分析

算例分析即利用本文提出的算法来求解无压渗流问题经典模型和异常体模型的自由面位置，通过与其他方法[11,15,18-20]在标准模型应用得到的结果进行对比来验证本文算法的有效性，在异常体模型中研究异常体的存在对渗流速度、水头分布、流体压力等渗流场参数的影响。

5.1 矩形均质坝模型(示例1)

示例1采用的经典模型为高11 m，宽5 m的矩形均质坝模型。上游库水的水头值H1=10 m，下游水头值H2=2 m。根据文献[17]、[21]中选取的渗透系数值，本部分渗透系数设置为1 m/s，底部边界设置为不透水边界。所采用的网格如图5(a)所示，初始自由面和计算得到的自由面位置如图5(b)所示。假设该模型的渗流区域为均质各向同性介质，初始自由面可以由Dupuit公式得到：

图5 示例1中网格尺寸及初始自由面与最终自由面结果对比

(25)

式中：H1、H2分别为上游库水、尾水的水头值，m；L为底部边界的长度，m。

示例1中本文方法对自由面计算结果与其他方法计算结果的对比如图6(a)所示，本文方法计算得到的水头等值图、渗流速度图及流体压力等值图如图6(b)～6(d)所示。由图6(a)可看出，本文方法计算得到自由面位置与其他方法得到的自由面位置具有很好的一致性，因此本文方法可以很好地预测自由面以及溢出点的位置。图6(d)中的流体压力可通过公式(1)计算得到，其零压力等值线与自由面位置相吻合。

图6 示例1中自由面本文方法与其他方法计算结果对比及渗流场各参数计算结果

5.2 具有单异常体矩形坝模型(示例2和示例3)

在土石坝中，由于介质渗透系数的不同，自由面的位置会受到一定的影响[21]。实际中常常需要根据自由面的位置和形状推测渗透异常区域的位置，因此研究异常体对自由面所产生的影响是十分有必要的。示例2和示例3的模型尺寸、边界条件和网格与示例1相同，但在模型中加入了不同渗透系数的异常体，将具有不同异常体位置(近自由面和远自由面)的两个模型分别设为示例2和示例3。分别对示例2和示例3计算表1中列出的4种不同的渗透系数情况下的自由面位置及流场分布，示例2和示例3背景区域的渗透系数设置为1 m/s。

表1 示例2和示例3模型中异常体不同渗透率的溢出点高程

与示例1相同。示例2近自由面渗透异常体位置及异常体不同渗透情况下的自由面计算结果见图7，模型流场分布见图8；示例3远自由面渗透异常体位置及异常体不同渗透情况下的自由面计算结果见图9，模型流场分布见图10。

图7 示例2近自由面渗透异常体位置及异常体不同渗透情况下的自由面计算结果

图8 示例2异常体不同渗透情况下的模型渗流场

图9 示例3远自由面渗透异常体位置及异常体不同渗透情况下的自由面计算结果

图10 示例3异常体不同渗透情况下的模型渗流场

由图7(b)可以看出，当异常体位置靠近自由面时，自由面受异常体影响较大。当异常体渗透系数高于背景区域渗透系数时，自由面向异常体方向弯曲；而异常体渗透系数低于背景区域时，自由面向异常体相反方向弯曲。且异常体渗透系数与背景区域渗透系数差值越大，则自由面受影响程度越大。由图9(b)可以看出，当异常体位置距离自由面较远时，自由面受异常体影响较小。由表1可知，溢出点位置受异常体的影响，异常体离自由面越近，对其位置的影响越大，异常体位置较远时，溢出点位置的变化不大。

图8和10显示了异常体不同渗透系数影响下的渗流速度变化。当异常体的渗透系数由低于背景区域逐渐增大至高于背景区域时，渗流速度逐渐增大，流向由被低渗透异常体排斥变为向高渗透异常体汇聚。

通过示例2和示例3计算结果的比较，发现异常体位置对自由面和渗流速度均有影响。异常体离自由面越近，对自由面影响越大。然而，由于自由面附近的渗流速度较小，渗流速度受异常位置的影响相对较小。异常体离自由面越远，对自由面影响越小，异常体位置对渗流速度的影响越大。

5.3 具有多异常体矩形坝模型(示例4)

示例4旨在探讨高渗透和低渗透多种不同渗透异常体对水头分布、渗流速度和流体压力的影响。各异常体的具体位置如图11(a)所示。背景区渗透系数仍设为1 m/s，低渗透系数和高渗透系数分别设为0.1和10 m/s，高渗和低渗异常体的渗透系数均与背景区域相差10倍。模型的边界条件和网格与示例1一致。

具有多种不同渗透异常体的自由面计算结果如图11(b)所示，水头等值图、渗流速度图及流体压力等值图如图12所示。由图12(a)、12(b)可看出，高渗透区域渗流速度增大，水头等值线变稀疏，低渗透区域排斥周围渗流，通过低渗透区域的渗流速度减小，低渗透区域水头等值线变密集。两个低渗透区域对水头和渗流流速的影响程度也不相同，深部的低渗透区域对周围水头和渗流流速的影响更加明显。在图12(c)中，流体压力随异常渗透率的不同而发生变化，流体压力等值线在高渗透区域变得稀疏，而在低渗透区域变得密集。