基于无单元Galerkin...—非饱和土石坝渗流正演模拟_戴前伟.pdf

第 卷 第 期 年 月水 资 源 与 水 工 程 学 报 ，收稿日期：；修回日期：基金项目：国家自然科学基金项目（）；国家重点研发计划项目（）作者简介：戴前伟（），男，湖南娄底人，博士，教授，博士生导师，主要从事电磁法方法及理论、工程地球物理勘探等方向的研究。通讯作者：雷 轶（），男，湖南邵阳人，博士后，主要从事电磁法和工程地球物理勘探方向研究。：基于无单元 法的饱和非饱和土石坝渗流正演模拟戴前伟，朱泽龙，韩行进，刘 杰，雷 轶（中南大学 地球科学与信息物理学院，湖南 长沙；中南大学 有色金属成矿预测与地质环境监测教育部重点实验室，湖南 长沙；五凌电力有限公司，湖南 长沙；中南大学 土木工程学院，湖南 长沙）摘 要：在土石坝渗流正演模拟中，仅研究饱和区域的渗流场并不能全面真实地反映出土石坝地下渗流状态，非饱和渗流区域的研究同样至关重要，利用无单元 法对饱和非饱和渗流域进行求解。首先从达西定律出发，推导了渗流方程及边界条件，其中详细推导了无单元 法 通过滑动最小二乘法构造形函数，同时利用罚函数的方法计算边界条件。然后，通过与 软件计算出的仅饱和渗流场以及饱和非饱和渗流场进行比较，证明了该方法的准确性及有效性。最后，通过不同的均质和非均质模型，研究了坝体中零压力线、水头值、孔隙压力、含水率的分布。无单元 法只需要通过节点来实现对全域渗流场的精确逼近，解决了对网格单元的依赖问题，与其他数值方法相比，其具有前期处理数据简单和精度高的优点，更适合稳定饱和非饱和渗流场这种复杂情况的正演模拟。关键词：饱和非饱和渗流；无单元 法；滑动最小二乘法；罚函数；土石坝中图分类号：文献标识码：文章编号：（），（，；，；，；，）：，：；研究背景饱和非饱和渗流分析对于岩土工程、地下水水文、水文地球物理以及高土石坝的稳定性等许多应用都具有重要意义。在渗流分析领域中，精确确定零压力面和溢出点的位置仍然是一项棘手的任务，事实上，在渗流问题的迭代求解过程中，零压力面的位置和溢出点的位置是不断变化的。目前，在渗流问题的研究中，很多学者仅将研究区域限定在饱和区，而非饱和区内的渗流问题对工程也有重要影响。非饱和区的渗透系数会随着孔隙压力的变化而变化，从而影响零压力面与溢出点的位置。因此，在求解土石坝渗流问题时仅考虑饱和区内的渗流并不全面，可能造成计算结果与实际情况偏差较大的问题。不同的数值模拟方法已经被国内外水文地质学者们应用于渗流研究中。在饱和非饱和渗流分析中，数值模拟方法主要分为有限元法（，）和有限差分法（，）两大类。在有限元计算中，网格离散化是必不可少的，但该方法存在一些缺点：（）有限元法不太适合于涉及几何变化、严重畸变、内部边界、裂纹扩展、地形起伏、近场源、形状设计优化等问题的应用。（）虽然区域划分的思想非常巧妙，但是分析区域的网格划分可能非常费力和耗时。（）通常很难提高局部精度。（）在模型中加入新节点时，需要重新网格化。（）对于不同类型的单元网格，由于有限元法中形函数的约束条件复杂，使形函数的确定较为困难。无单元 法（，）具有不受网格限制的优势，被广泛应用于各种工程和地球科学模拟中。最初由 等提出，用于解决弹性问题，而后在计算电磁学、计算力学、地球物理学、渗流场 等多个领域被越来越多的学者关注。在 中，节点和背景单元相互独立，消除了有限元法中严格要求的网格离散化所带来的局限性。不依赖网格使该方法变得更加灵活，其形函数的构造采用局部领域中可自动生成的滑动最小二乘法（，），形函数对于任意不规则的节点分布具有灵活性。等运用变分多尺度插值方法改进了 中无插值滑动最小二乘法构造具有插值性质的形函数，使渗流场数值模拟能够达到更高的精度和稳定性。曹阳等在运用 中对目前已提出的数十种关于如何方便准确地处理边界条件的方法进行了归纳总结，比较了这些方法的优缺点，最后提出了展望。王成华等总结分析了大量的全域渗流数值模型，介绍了饱和非饱和渗流数值分析的研究现状，讨论了现有分析方法存在的问题。通过提出合理的边界条件提高了饱和非饱和渗流数值分析的准确性。本文运用 来解决饱和非饱和渗流场问题。减少了对网格单元的依赖，提高了求解精度，其求解过程通过滑动最小二乘法构造水头场函数的形函数，再利用罚函数的方法求解边界条件，最后通过求解控制方程，得到场节点的水头值分布情况，然后利用边界流量控制和水头值更新溢出点。在迭代过程中修正边界溢出点的位置，规定各节点与高程的容差，直至达到收敛标准，得到最终稳态的零孔隙水压力面（即零压力面）和溢出点。饱和非饱和渗流场的数学描述 达西定律达西定律可以用来表示饱和土中水的渗透状态，达西认为土体中水的流速受其水力梯度的影响，二者之间成正比关系：（）式中：为渗流速度，；为渗透系数，；为渗流水头值，；为垂直方向坐标，。由公式（）可得出饱和域内的渗流速度为：|（）式中：为水平方向渗流速度，；为垂直方向渗流速度，；为水平方向渗透系数，；为垂直方向渗透系数，。饱和域中的流体连续性方程为：（）水 资 源 与 水 工 程 学 报 年在这种渗透条件下，将公式（）代入公式（）可以得到二维稳定渗流控制方程：()()（）达西定律最初仅针对饱和渗流，但经过大量的研究发现，其同样适用于饱和非饱和渗流，唯一的区别在于非饱和区渗透系数与孔隙压力相关。边界条件为了方便描述饱和非饱和渗流问题，以图 为例，其中 为零压力面，零压力面以上的部分 为非饱和域，以下的部分 为饱和域。图 土石坝中的饱和非饱和渗流示意图如图 所示，饱和非饱和定解问题一般分为以下 类边界条件：（）图 中、段，为已知水头边界，需满足：（）（）式中：、分别为已知上、下游水位高度，。（）图 中 段，为不透水边界，需满足：（，）（）式中：为不透水边界 的外法线方向；为 的水流单宽流量，；边界 为 段；为高程，。（）图 中 段，为溢出面边界，此边界为疑似边界，在求解溢出点以下位置时，其渗流量，水头值 。反之，溢出点以上位置的渗流量，。（）图 中 段，为零压力面，需满足 。无单元 法的数学基础 建立在 的数学基础之上，的基本原理为利用加权最小二乘法对构造的近似场函数进行加权离散处理。因此无单元 法具有局部函数紧支性的特点。基本原理如图 所示，将二维渗流场模型分成 个节点（实心点）。假定二维渗流场场函数为（），在二维空间内，坐标函数为场节点，给定场函数（）在个实心点的值为：（）（，）（）式中：为每个场节点的场函数。图 滑动最小二乘法近似示意图通过可以将二维渗流场模型分成个节点场值，于是可以得到一个近似场函数（）：（）（）（）（）（）（）式中：（）为维基函数矢量；（）为维系数矢量。对于二维情况，（）有以下 种选取方案：（）（，）（线性基，）（）（，）（二次基，）（）（，）（三次基，）|（）式中：，为高斯积分点位置坐标，。本方法中的（）是一个需要求解的变量，为空间场节点 的函数，即：（）（），（），（）（）（）的求解需要构建一个加权的离散范数：（）（）（）（）式中：（）为点的权函数；为点一个有限范围内的场节点；为权函数影响的节点个数。将 取极小值求解（）得：（）公式（）通过公式（）展开可得：（）（）（）（）其中：（）（）（），（）（）（）（），（）（），（）（），|（）第 期 戴前伟，等：基于无单元 法的饱和非饱和土石坝渗流正演模拟由公式（）可导出：（）（）（）（）式中：为二维渗流场的各节点场函数的矢量矩阵。将公式（）代入近似场函数式（）可得：（）（）（）（）（）（）式中：（）为节点 的形函数。将公式（）代入由式（）可得（）的表达式为：（）（）（）（）（）上式可简化为：（）（）（）（）（）式中：、均为矢量形式。上述这种构建形函数的方法称为滑动最小二乘法（）。权函数的选择由于选取的权函数不同，的精度及其计算复杂程度会发生相应变化，所以如何选取权函数是运用 的重点。权函数选取方式应该遵循以下规律：（）权函数的值必须为非负数。（）权函数的最大值应该取在每一个场节点上，并且离场节点越远其值越小，而在场节点的影响半径外其值为。（）权函数应是一个连续可导的函数。因此，可将权函数表示为两点之间的距离：（）（）（）式中：（）为 与 的距离，。权函数的选择是随机的，但需满足以上条件。本文将选取指数函数作为权函数：（）（）（）（）（）（）|（）式中：为节点 的影响域半径，；，其中 为相连两节点的距离，值范围为，本文取；为非负的整数，本文取 。在均匀分布的节点上 取值为：（）式中：根据节点分布的密度取值；由（）基函数确定；为常数，一般取 ，本文取 。无单元 法在渗流场中的数值模拟 基本方程由第 部分理论可知，水头场函数的数学表达式为：（）式中：为节点水头值，；为每个高斯积分点需要求解的节点个数；为节点形函数。整个系统的总能量可以表示为：（）根据变分原理可知：（）由公式（）可得整体平衡方程为：（）式中：为渗透系数矢量形式；为每个节点的水头值矢量形式；为右端项。渗流平衡方程渗流场的基本微分方程的能量泛函可表示为：()()（）通过变分原理对公式（）求极值可得：（）（）（）因此有：（）（）（）（）（，）（）其中：（）（），（），（）（）|（）水头边界条件采用罚函数方法处理已知的水头边界，设水头边界上含有一个强透水层，通过强透水层可以表示出水头边界上的能量泛函：()（）（）（）水 资 源 与 水 工 程 学 报 年式中：为强透水层的渗透系数，；为强透水层厚度，；为强透水层的测压水头，；为水头边界的测压水头，。根据变分原理，对公式（）求极值，即：（）其中：（）式中：为常数，取值范围为 ；为边界上节点最小间距的负值，。即得：（）（）|（）由公式（）和（）可得：（）（）（）（）（）即：|（）|利用高斯积分法可求得 和 值。数值实现方法本文采用四边形网格，选取四边形网格特有的高斯积分点及特定的加权系数。在每个矩形网格的 个场节点内，需要布置多个高斯积分点，这里确定为 的平方。当高斯积分点过少时将会影响计算精度；当高斯积分点过多时，将会增大计算时间，影响计算效率。四边形网格内包含场节点的背景网格需要进行高斯积分计算，而没有场节点的高斯积分无需计算。无单元 法求解如图 所示。图 无单元 法求解示意图由高斯积分计算方法可知任一函数（，）在二维域内积分为：（，）（，）（）式中：为全域面积，；为全域内高斯积分点的总个数；为每个高斯点的积分权。同理，由高斯积分计算方法可知任一函数（，）在边界上的积分为：（，）（，）（）式中：为整个边界，；为全边界内高斯积分点的总个数；为每个高斯点的积分权。饱和非饱和渗流问题的数值程序饱和非饱和渗流问题的数值程序模拟运算步骤如下：（）生成高斯积分点，覆盖整个求解域 。（）在全域范围内生成场节点。（）循环高斯积分点：确定当前积分点的材料区。如果当前积分点不属于任何材料区域，则退出当前点并移动到下一个积分点；确定当前积分点的影响域；计算形函数及其导数；计算当前积分点的渗透矩阵；通过公式（）将当前积分点的渗透矩阵组合成总渗透矩阵；将当前积分点的右端项矩阵组合成总右端项。（）利用公式（）实现边界条件，形成全局方程组。（）求解全局方程，得到节点水头值。（）判定结果是否收敛。如果收敛，则进行下一步；如果不收敛，则根据零压力面及溢出点位置更新渗透系数重新求解零压力面。（）得到零压力面及其他相关结果。方法更新非饱和区渗透系数 于 年提出了相对简单的土壤含水率压头曲线方程，并通过方程的特殊形式推导出了非饱和区渗透率的封闭解析表达式。根据 开发的技术，可以使用需要曲线拟合参数的封闭形式方程来生成非饱和域体积含水率函数：（）（）式中：、均为控制含水率函数的参数，第 期 戴前伟，等：基于无单元 法的饱和非饱和土石坝渗流正演模拟，这些参数可以通过土壤水分保持模型与实验数据拟合得到；为饱和含水率；为残余含水率；为孔隙水压力，。从饱和渗透率和体积含水率函数来估计非饱和区渗透率函数。非饱和区渗透率的封闭形式方程为：（）（）（）（）（）式中：（）为非饱和区渗透率，；为饱和渗透率，。算例分析该部分利用 来求解饱和非饱和全域渗流问题经典模型的零压力面位置。研究全域渗流场需要对非饱和区的渗透率进行不断迭代更新，通过 方法，利用每一次迭代计算的背景网格单元平均孔隙水压力、水头和输入饱和渗透率，对非饱和区的渗透率进行更新迭代。并通过与其他的土石坝渗流模拟方法进行对比分析论证本方法的准确性及实用性。经典矩形均质模型模型 为矩形均质坝模型，设置坝高为 ，坝宽为 ，背景网格数为 ，节点均匀布置。节点间距为 ，上、下游水位分别为 和，背景饱和渗透率为 。背景饱和含水率 ，残余含水率 ，为控制含水率的参数（公式（）。本文方法与 软件对模型 渗流零压力线的计算结果对比见图。由图 可知，本文采用的 与 软件计算得到的饱和非饱和渗流零压力线基本一致，且在相同的网格情况下，由于 具有解高次连续的优点，其迭代次数比 的迭代次数少，提高了计算效率。根据图 中溢出点的差异可看出，当仅考虑饱和域时，上部分非饱和域不参与计算，零压力面及溢出点明显上移，计算结果验证了非饱和域对零压力面的影响及本方法的正确性。图（）、（）为 计算出的孔隙水压力与水头等值线图。由图（）中可以清晰地看到零压力面的具体位置。图（）、（）为渗透率、含水率分布图，本文采用 公式，通过孔隙水压力与渗透率、含水率之间的关系，计算得到渗透率分布和含水率分布。根据公式（）、（）可知，非饱和域内的渗透率和含水率随着孔隙水压力的改变而改变。图 不同方法对模型 渗流零压力线的计算结果图 及 公式计算的模型 水文状态变量分布 多个异常体矩形坝模型对于土石坝渗流反演，其实质是利用水头、流量等观测数据求解渗流介质的初始条件、渗透系数和边界条件等的反演问题。土石坝中介质的非均匀性会影响零压力面的位置和形状，而在目前的无压渗流研究中，很少有针对穿过零压力面的异常体的分析。实际情况表明，靠近零压力面和穿过零压力面的异常体对零压力面的影响最大。因此，针对坝体内含异常体和异常体穿过零压力面的渗流研究是非常必要的。水 资 源 与 水 工 程 学 报 年模型 为含异常体（不同渗透系数）的矩形坝模型，如图（）所示。模型 中背景饱和渗透率 ，；红色异常体为高渗透率异常体，；黑色异常体为低渗透率异常体，。其余参数取值与模型 相同。模型 渗流零压力线、孔隙水压力及水头计算结果如图（）（）所示。图（）显示，与 软件计算结果基本一致，含异常体的非饱和域对零压力面位置的影响与仅饱和域相比有明显差异。从图（）异常体位置和图（）零压力面位置可知，零压力面是从上部两块异常体中间穿过，零压力面受到其附近异常体的影响明显。从图（）和图（）可以看出，高渗透率异常吸引渗流，水头等值线稀疏；低渗透率异常排斥渗流，水头等值线密集。且异常体渗透率与背景渗透率差异越大，对零压力面的影响越大。图 为 公式计算的模型 水文状态变量分布图。从图（）可以看出，非饱和区内的渗透率随孔隙水压力的改变而改变。图 模型 示意图及计算结果图 公式计算的模型 水文状态变量分布 含排水沟的等腰梯形均质模型为了更好地与工程实际相结合，模型 设计为含排水沟的均质等腰梯形坝，背景网格数为 ，节点均匀布置。纵坐标节点间距为 ，横坐标节点间距为 。上游水位为，排水沟长度为，饱和渗透系数 ，饱和含水率 。其他参数取值与模型 相同。本文方法与 软件对模型 渗流零压力线的计算结果对比见图（圆点标记部分为排水沟区域）。模型 水文状态变量分布见图。图 不同方法对模型 渗流零压力线的计算结果第 期 戴前伟，等：基于无单元 法的饱和非饱和土石坝渗流正演模拟由图 可以看出，与 软件计算方法得到的饱和非饱和渗流零压力线基本一致。分析图 可知：（）孔隙水压力和水头计算结果不受坝体形状影响，计算过程中节点分布均匀，不需要采用不规则化网格（图（）、（）；（）背景网格渗透率及含水率的更新与网格内是否包含水头节点有关，当背景网格无水头节点分布时，渗透系数不更新（图（）、（）。图 及 公式计算的模型 水文状态变量分布图 结 论本文采用无单元 法（）求解土石坝饱和非饱和渗流问题，通过与 软件的计算结果对比验证了所提方法的有效性。针对不同类型的土石坝模型进行求解，得到如下结论：（）只需要通过节点来实现对全域渗流场的精确逼近，解决了对网格单元依赖的问题，与其他数值方法相比，其前期数据处理步骤简单且计算精度高，更适合饱和非饱和渗流场的数值模拟。（）通过与 软件计算结果比较，验证了本文方法的正确性。坝体内含异常体模型的计算结果表明，靠近零压力面和穿过零压力面的异常体对零压力面的影响最大。梯形坝模型计算结果表明，本方法不受坝体形状的影响，计算过程中节点分布均匀，不需要采用不规则化网格。（）通过对饱和非饱和渗流的求解，可以获取全域的水文变量信息（如含水率、饱和度等），这将为水文地球物理研究提供理论基础。参考文献：，：，：，（）：，：，（）：，：，水 资 源 与 水 工 程 学 报 年 ，戴前伟，孔重阳，雷 轶，等 基于光滑有限单元法的土石坝无压渗流场数值模拟 水资源与水工程学报，（）：严绍洋，李亮辉，高燕希，等 公路隧道开挖渗流场的有限差分法分析 中外公路，（）：豆海涛 基于渗漏水的隧道渗流场有限差分法分析 铁道勘察，（）：，：，（）：，（）：，（）：，（）：，（）：，（）：，（）：麻昌英，柳建新，郭荣文，等 耦合有限单元法扩边的直流电阻率勘探无单元 法正演 地球物理学报，（）：李俊杰，严家斌，皇祥宇 无单元 法大地电磁三维正演模拟 地质与勘探，（）：刘昌军，徐甲存，唐 波，等 尾矿坝二维渗流场的无单元模拟方法研究 西部探矿工程，（）：，：曹 阳，陈莹婷，姚林泉 无单元 方法施加本质边界条件研究进展 力学季刊，（）：王成华，张燕青 饱和 非饱和土渗流数值分析方法综述 建筑技术，（）：，（）：邓高阳 地下洞室三维渗流场溢出边界模拟与渗流反演分析 武汉：武汉大学，第 期 戴前伟，等：基于无单元 法的饱和非饱和土石坝渗流正演模拟