饱和无压渗流模型的虚拟单元法求解

Virtual Element Method for Saturated Unconfined Seepage Flow Analysis

摘要

局部修改网格是固定网格法求解饱和无压渗流模型的重要途径之一。以有限单元法为计算工具时，单元形状受到限制导致网格修改过程复杂，阻碍了局部修改网格法的发展。为解决此问题，引入适应于任意形状多边形单元的虚拟单元法，采用局部修改网格法求解二维饱和无压渗流模型。采用1阶虚拟单元表达单元的水头变化并建立单元基本方程，推导单元渗透矩阵，并揭示1阶虚拟单元与有限单元数学上的联系。利用虚拟单元法对单元形状的包容性设计网格修改规则，建立整体渗透矩阵和整体流量向量的更新算法以减少计算成本。采用均质矩形土坝、均质梯形土坝和非均质矩形土坝等3个数值算例检验本文算法的计算精度和效率。结果显示：本文算法具有可靠的计算精度，获取的均质矩形土坝自由面与解析解相对误差为0.45%，均质梯形土坝和非均质矩形土坝自由面与改进流形单元法结果基本吻合；迭代过程中自由面节点将加密，在提升计算精度的同时，也使得自由面形态更为平滑，但导致计算效率与传统局部修改固定网格方法相比略微下降。

Abstract

Objective Local mesh modification is one of the important approaches for solving the saturated unconfined seepage model when the fixed mesh method is used. When the finite element method is employed as the calculation tool, the complexity of the mesh modification process, resulting from the restricted shapes of the elements, hinders the development of local mesh modification methods. Therefore, the virtual element method, which is adaptable to arbitrarily shaped polygonal elements, is introduced to solve the two-dimensional saturated unconfined seepage model using the local mesh modification method. Methods First, the basic mathematical equation of the saturated unconfined seepage model was described. Then, the change in water head within an arbitrary polygonal element was expressed using the first-order virtual element, leading to the establishment of the first-order virtual element formulation for seepage analysis. The trial function space V₁(Ω^e) governed the water head in each element, with the water heads at the vertices defined as the degrees of freedom of the element. The element conductivity matrix and flow vector were derived using the Galerkin method, and the procedure for computing the element conductivity matrix was deduced based on the principles of the virtual element method. The mathematical relationship between the first-order virtual element and classical finite element methods was elucidated. Next, the versatility of element shapes within the virtual element method was utilized to redesign the mesh modification rules of the local mesh modification method, developing a comprehensive solution flow for the saturated unconfined seepage model. The mesh was modified by removing elements located above the free surface while retaining those below it. In cases where an element was intersected by the free surface, it was first removed, and the portion below the free surface was treated as a new element. Update algorithms for the global conductivity matrix and the global flow vector were established to reduce computational costs. Finally, three numerical examples, homogeneous rectangular earth dams, homogeneous trapezoidal earth dams, and heterogeneous rectangular earth dams, were employed to evaluate the calculation accuracy and efficiency of the proposed method. The performance of the proposed method was also compared to other methods based on the finite element method and the numerical manifold method. Results and Discussions The numerical accuracy of the proposed method was evaluated using the analytical solution of the rectangular homogeneous earth dam model as a benchmark. The maximum error identified along the free surface, which occurred at the exudation point on the right side of the earth dam, was a relative error of 0.45%. This result indicated that the proposed method exhibited good numerical accuracy. The method was applied to solve the trapezoidal homogeneous earth dam model and the inhomogeneous rectangular earth dam model. Due to the absence of analytical solutions for these two models, the results were compared to those obtained using Seep software, the local mesh modification method based on the finite element method, and the improved numerical manifold method, which is widely recognized for its high computational accuracy. For the trapezoidal homogeneous earth dam model, Seep software produced the lowest free surface position, while the local mesh modification method based on the finite element method produced the highest. The proposed method yielded a free surface position comparable to that obtained using the improved numerical manifold method, although the exudation point on the right side of the dam was slightly higher. For the inhomogeneous rectangular earth dam model, all methods except the local mesh modification method based on the finite element method exhibited a rapid decline in free surface levels as they approached the sudden change line of the permeability coefficient, followed by slight fluctuations after crossing this line. When ranking the free surfaces of the inhomogeneous rectangular earth dam model obtained using the four methods, the results closely matched those of the trapezoidal homogeneous earth dam model. Under conditions of comparable numerical accuracy, the proposed method employed the virtual element method, which is highly compatible with the finite element method, avoiding the introduction of complex mathematical concepts such as numerical manifolds. Under the same computer configuration, the running times of the proposed method for solving the trapezoidal homogeneous earth dam and the inhomogeneous rectangular earth dam were 14.4 and 8.1 seconds, respectively. In comparison, the running times of the local mesh modification method based on the finite element method were 13.6 and 7.2 seconds, respectively. Therefore, the computational efficiency of the proposed method was slightly lower than that of the local mesh modification method based on the finite element approach. This difference was attributed to the mesh modification strategy of the proposed method, which directly cut the fixed mesh using the free surface, reducing computational time. However, the encryption and intervention of free surface nodes during the iterative process were relatively time-consuming. The encryption of free surface nodes improved calculation accuracy and produced a smoother free surface profile. Therefore, the slight reduction in computational efficiency of the proposed method was considered acceptable. Conclusions The results demonstrate the validity and accuracy of the proposed method for solving saturated unconfined seepage. The method overcomes the limitations associated with the restricted element shapes in the local mesh modification approach within the finite element method, establishing the local mesh modification method as a reliable tool for analyzing saturated unconfined seepage. In addition, this advancement expands the application domains of the virtual element method.

Graphical abstract

关键词

虚拟单元法 / 饱和无压渗流 / 自由面 / 局部修改网格 / 计算精度

Key words

virtual element method / saturated unconfined seepage / free surface / local mesh modification / computational accuracy

本刊网刊

无压渗流问题是水利工程领域的经典问题，可采用饱和无压渗流模型和饱和‒非饱和渗流模型描述后进行数值求解。尽管饱和无压渗流模型缺乏描述非饱和区中水分迁移过程的能力，但其所需参数为易得到的渗透系数，可解决工程实践中确定土‒水特征曲线等饱和‒非饱和渗流模型参数的困难^[1]。因此，饱和无压渗流模型求解的非线性程度更高^[2]，其求解方法一直备受关注^[3‒4]。

有限元法（FEM）、流形元法^[5‒7]和无网格法^[8]等均被用于求解饱和无压渗流模型，其中，FEM应用最多。由于自由面的位置预先未知，饱和无压渗流模型一般须迭代求解。根据迭代过程中网格处理的差异，FEM求解存在调整网格法和固定网格法两类方法。调整网格法视自由面为可变边界，迭代过程中修改自由面位置使网格发生变形，直至自由面位置稳定为止。调整网格法存在计算量大、网格可能畸变和非均质土体难收敛等局限，虽然引入等几何分析等^[9‒11]技术后，性能有所改善，但与固定网格法相比，目前所受关注较少。

固定网格法的核心在于计算过程中网格基本保持不变，求解方法可分为变分不等式法和直觉化方法^[12]。变分不等式法具有严密的数学基础，无须关注自由面的位置，数值精度和求解稳定性表现优异^[13‒15]，但由于理论过于复杂，尚未被工程师广泛接受。直觉化方法可基于调整流量^[16‒17]、调整单元渗透矩阵^[18]和局部修改网格^[19‒20]3种路径实现，均需考虑自由面的位置：调整流量需要根据高斯点和自由面位置的关系确定该点对单元流量的贡献；调整单元渗透矩阵需要对自由面经过的单元用Heaviside函数进行人为微调；局部修改网格则直接修改局部的单元和节点信息，使其准确地反映自由面位置，然后将自由面以下作为计算区域。实际操作过程中，由于FEM对单元形状的限制，局部修改网格往往须精心设计修改规则，网格修改过程复杂。受制于此，尽管局部修改网格无需引入额外的物理或数学概念，且最为简单直观，但近年来未能进一步发展。

Beirão等^[21‒22]于2013年提出虚拟单元法（VEM），一种新型数值方法，可适用于任意形状的多边形单元^[23‒25]，鉴于其在处理悬挂节点、接触和晶体变形等问题上具有优势，受到国内外计算力学学者的广泛关注^[26‒28]。VEM与FEM高度兼容^[29]的同时，极大放松了对单元形状的限制，因此有破解局部修改网格法现有局限的潜力。本文尝试采用1阶VEM求解二维饱和无压渗流模型：首先，按VEM理论基本框架推演饱和无压渗流问题的单元基本方程及相关矩阵和向量的计算方法，建立VEM求解格式；然后，采用局部修改网格法迭代求解模型，利用VEM对单元形状的包容性设计网格修改规则，并建立整体渗透矩阵和流量更新算法以减少计算成本；最后，采用数值算例对本文算法的计算精度和效率进行检验。

1 饱和无压渗流模型

土坝饱和无压渗流示意图如图1所示。图1中，Ω_w为饱和渗流区域，AB和CD分别为上游面和下游面，BC为不透水地基，AE为自由面，DE为出渗面，H_up为上游水位高度，H_down为下游水位高度。

流动区域中任意一点P的总水头H等于该点的位置水头与压力水头之和，具体的表达式为：

H=y+p/{\gamma }_{\mathrm{w}}（1）

式中：y为该点的位置水头，等于该点的纵坐标；p为孔隙水压力；γ_w为水的重度。

假定P点的流速 v =(v_x,v_y )^T满足达西定律，其中，v_x 和v_y 为渗流速度 v 在x方向和y方向上的分量，则：

\mathit{v}=-\mathit{D}\nabla H（2）

式中， D 为P点处的2阶渗透张量，\nabla为哈密顿算子。

当水与土体都不可压缩时，流体在多孔介质中的连续性方程为：

\nabla \mathit{v}=\mathrm{0}（3）

将式（2）代入式（3），即可得到渗流微分方程：

\nabla \mathit{D}\nabla H=\mathrm{0}（4）

第一类边界条件直接给出水头在边界上的值，称为水头边界条件。在上游面AB上应满足：

H=H_{\mathrm{u}\mathrm{p}}（5）

在下游面CD上应满足：

H=H_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}（6）

在自由面AE上应满足：

H=y（7）

在出渗面ED上应满足：

H=y（8）

第二类边界条件给出水头在边界上法向导数的值，称为流量边界条件。在不透水地基BC上应满足：

q_{\mathit{n}}=-{\mathit{n}}^{\mathrm{T}}\mathit{v}=\mathrm{0}（9）

式中：q_{\mathit{n}}为边界法向流量； n 为边界BC的单位外法线方向向量， n =(n_x,n_y )^T，其中，n_x 为单位外法线向量在x方向分量，n_y 为单位外法线向量在y方向分量。

在自由面AE上应满足：

q_{\mathit{n}}=\mathrm{0}（10）

由于自由面AE的位置事先未知，因此，确定其位置就成了饱和无压渗流模型求解的主要内容。

2 1阶虚拟单元法求解渗流问题的计算格式

2.1 基于虚拟单元法的水头插值与方程

数学原理上，VEM具有与FEM相同的变分背景^[22]，可视为FEM向一般的多边形/多面体单元的拓展。二维分析时FEM单元形状限定为三角形、矩形和四边形，而VEM则不受此限制。平面域Ω的多边形网格离散如图2所示，平面域Ω可采用一般的多边形网格离散后用VEM求解，形成单元域Ω^e，其边数l\ge \mathrm{3}，经典的FEM网格可视为一般多边形网格的特例。

VEM的特点在于物理量测试函数仅有形式化定义，而无显式表达，因此被称为“virtual”。求解区域离散后，VEM定义单元域Ω^e上物理量流速v的测试函数空间V_k (Ω^e)为：

V_k({\Omega }^{\mathrm{e}})=\left\{\begin{array}{l}V\in H^{\mathrm{1}}({\Omega }^e)\bigcap C^{\mathrm{0}}(\partial {\Omega }^e),\\ \mathrm{\Delta }V\in P_{k-\mathrm{2}}({\Omega }^e),\\ V_{\Gamma }\in P_k(\mathrm{\Gamma })\forall \mathrm{\Gamma }\in \partial {\Omega }^e\end{array}\right.（11）

式中：V为测试函数；V_{\mathrm{\Gamma }}为测试函数在边界上的值；∂Ω^e代表单元边界；H¹为1阶Sobolev空间；C⁰为连续函数空间；\mathrm{\Delta }为拉普拉斯算子；Γ构成\partial {\Omega }^{\mathrm{e}}的一条边；P_k 为由不超过k阶的多项式构成的函数空间，k为不小于1的整数。

使用V_k (Ω^e)的单元可称为k阶VEM单元。根据定义，必然存在P_k (Ω^e)⊂V_k (Ω^e)。

k阶VEM单元的自由度包括：1）单元每个顶点处的v值；2）在每条边上执行k+1点Gauss‒Lobatto积分所需的k‒1个内积分点处的v值；3）在单元域上v的0至k‒2阶矩。k阶VEM单元可精确表达单元域内k阶多项式分布的物理量，因此，增大k取值可得到单元内物理量的更高阶分布情况，但其代价是单元自由度数量的增加。

本文拟采用1阶VEM求解二维饱和无压渗流模型，单元域Ω^e内水头H∈V₁(Ω^e)，单元每条边上水头均为线性分布。此时单元的自由度仅包含各顶点处水头值，自由度向量 H^e表达为：

{\mathit{H}}^{\mathrm{e}}={(\begin{array}{ccc}{H}_{\mathrm{1}}& H_{\mathrm{2}}& \cdots \end{array}\begin{array}{ccc}{H}_i& \cdots & H_n\end{array})}^{\mathrm{T}}（12）

式中，H_i 为单元第i个顶点处水头值，n为单元顶点数。

单元域Ω^e内任意一点的水头值可表达为如下插值形式：

H={\mathit{\phi }}^{\mathrm{e}}{\mathit{H}}^{\mathrm{e}}（13）

式中， φ^e为规范化的单元基函数向量，可表示为：

{\mathit{\phi }}^{\mathrm{e}}=(\begin{array}{cccc}{\phi }_{\mathrm{1}}& {\phi }_{\mathrm{2}}& \cdots {\phi }_i\cdots & {\phi }_n\end{array})（14）

式中，φ_i 具备Dirac‒Delta性质：

{\phi }_i({\mathit{x}}_j)={\delta }_{ij}（15）

式中：{\mathit{x}}_j为单元内的坐标；δ_ij 为克罗内克函数，当i=j时，{\phi }_i({\mathit{x}}_j)={\delta }_{ij}=1，当i\nej时，{\phi }_i({\mathit{x}}_j)={\delta }_{ij}=0，其中，j表示非i节点。

根据H∈V₁(Ω^e)必然有φ_i ∈V₁(Ω^e)，可知φ_i 在每条边上线性分布。

使用伽辽金格式加权余量法，取式（13）中的单元基函数为权函数，可建立VEM单元的方程：

{\mathit{K}}^{\mathrm{e}}{\mathit{H}}^{\mathrm{e}}={\mathit{F}}^{\mathrm{e}}（16）

式中： F^e为单元流量向量； K^e为单元渗透矩阵，可表示为：

{\mathit{K}}^{\mathrm{e}}=\left[\begin{array}{cccc}{k}_{\mathrm{11}}& k_{\mathrm{12}}& \cdots & k_{\mathrm{1}n}\\ k_{\mathrm{21}}& k_{\mathrm{22}}& \cdots & k_{\mathrm{2}n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{n\mathrm{1}}& k_{n\mathrm{2}}& \cdots & k_{nn}\end{array}\right]（17）

单元渗透矩阵中元素k_ij 的计算式为：

k_{ij}={\int }_{{\Omega }^{\mathrm{e}}}(\nabla {\phi }_{i}D\nabla {\phi }_j)\mathrm{d}\mathit{x}（18）

式中， x 为单元内任意点坐标。

F^e的计算式为：

{\mathit{F}}^{\mathrm{e}}={\int }_{{\mathrm{\Gamma }}^{\mathrm{e}}}(({\mathit{\phi }}^{\mathrm{e}}{)}^{\mathrm{T}}{q}_{\mathit{n}})\mathrm{d}s（19）

式中，Γ^e为单元边界，s为弧长参数。

由式（9）和（10）可知，单元流量向量 F^e为0向量。由于φ_i 在每条边上线性分布且有Dirac‒Delta性质，水头边界条件的施加也不存在任何困难。因此，采用1阶VEM求解二维饱和无压渗流模型时，由于基函数φ_i 缺少显式表达，单元渗透矩阵元素k_ij 的积分表达无法直接求取，因此需要构造相应的数值计算方法。

2.2 基函数在多项式函数空间的投影

为计算类似于式（18）的内积积分，VEM基于内积规则构造V_k (Ω^e)中任意元素v向多项式函数空间P_k (Ω^e)的投影Πv∈P_k (Ω^e)。

对于任意元素v∈V₁(Ω^e)，本研究引入正交条件：

{\int }_{{\Omega }^{\mathrm{e}}}(\nabla v\nabla p)\mathrm{d}\mathit{x}={\int }_{{\Omega }^{\mathrm{e}}}(\nabla \mathrm{\Pi }v\nabla p)\mathrm{d}\mathit{x},\forall p\in P_{\mathrm{1}}({\Omega }^{\mathrm{e}})（20）

构造投影Πv为V₁(Ω^e)→P₁(Ω^e)。

由于Πv∈P₁(Ω^e)，可假定其表达式为：

\mathrm{\Pi }v=s_{\mathrm{1}}{p}_{\mathrm{1}}+s_{\mathrm{2}}{p}_{\mathrm{2}}+s_{\mathrm{3}}{p}_{\mathrm{3}}（21）

式中：s₁、s₂和s₃为相应系数；p₁、p₂和p₃为P₁(Ω^e)的3个基函数，文献[22]中，在p₁、p₂和p₃的表达式中引入了单元形心位置等参数，但其本质上对应一次多项式函数空间的3个经典基函数1、x和y，本研究定义为：

p_{\mathrm{1}}=\mathrm{1},p_{\mathrm{2}}=x,p_{\mathrm{3}}=y（22）

利用式（20）中p的任意性建立求解s₁、s₂和s₃的方程时，基函数p₁梯度为0向量，导致方程欠定。为解决此问题，VEM约定任意元素v∈V₁(Ω^e)向0阶多项式空间P₀(Ω^e)的投影为：

P_{\mathrm{0}}(v)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}v({\mathit{x}}_i)（23）

式中，{\mathit{x}}_i为单元第i个顶点的坐标向量。

联立式（20）、（21）和（23），可建立方程：

\begin{array}{c}\left[\begin{array}{ccc}{P}_{\mathrm{0}}(p_{\mathrm{1}})& P_{\mathrm{0}}(p_{\mathrm{2}})& P_{\mathrm{0}}(p_{\mathrm{3}})\\ \mathrm{0}& {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{2}}\nabla p_{\mathrm{2}})\mathrm{d}\mathit{x}& {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{2}}\nabla p_{\mathrm{3}})\mathrm{d}\mathit{x}\\ \mathrm{0}& {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{3}}\nabla p_{\mathrm{2}})\mathrm{d}\mathit{x}& {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{3}}\nabla p_{\mathrm{3}})\mathrm{d}\mathit{x}\end{array}\right]\cdot \\ \left(\begin{array}{c}{s}_{\mathrm{1}}\\ s_{\mathrm{2}}\\ s_{\mathrm{3}}\end{array}\right)=\left(\begin{array}{c}{P}_{\mathrm{0}}(v)\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{2}}\nabla v)\mathrm{d}\mathit{x}\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla p_{\mathrm{3}}\nabla v)\mathrm{d}\mathit{x}\end{array}\right)\end{array}（24）

方程中涉及\nablav的相关积分，可采用高斯公式转化为边界积分，然后利用v在每条边上为线性分布的特点正确计算。由式（24）解出s₁、s₂和s₃后，整理可得：

\begin{array}{l}\mathrm{\Pi }v=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{v}_i+\left(x-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{x}_i\right)\frac{\mathrm{1}}{A}{\int }_{\mathrm{\Gamma }}v{n}_x\mathrm{d}s+\\ \left(y-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{y}_i\right)\frac{\mathrm{1}}{A}{\int }_{\mathrm{\Gamma }}v{n}_y\mathrm{d}s\end{array}（25）

式中，v_i 为第i个顶点处v的值，x_i 和y_i 分别为第i个顶点在x和y坐标的值，A为单元面积。

φ_i 为V₁(Ω^e)的元素，且具备Dirac‒Delta性质，针对φ_i 执行式（25）可得：

\begin{array}{l}\mathrm{\Pi }{\phi }_i=\frac{\mathrm{1}}{n}+\frac{l_i{n}_{ix}+l_{i+\mathrm{1}}{n}_{i+\mathrm{1},x}}{\mathrm{2}A}\left(x-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{x}_i\right)+\\ \frac{l_i{n}_{iy}+l_{i+\mathrm{1}}{n}_{i+\mathrm{1},y}}{\mathrm{2}A}\left(y-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{y}_i\right)\end{array}（26）

式中，l_i 、l_i+1为与节点i关联的两条邻边的长度，n_{ix}、n_{iy}为与节点i关联的第i条边的单位外法线向量的x、y分量，n_{i+\mathrm{1},x}、n_{i+\mathrm{1},y}为与节点i关联的第i+1条边的单位外法线向量的x、y分量。

计算Πφ_i 所使用的几何信息如图3所示。

2.3 单元渗透矩阵元素的计算

建立基函数φ_i 在P₁(Ω^e)的投影后，可将基函数φ_i 分解为：

{\phi }_i=\mathrm{\Pi }{\phi }_i+({\phi }_i-\mathrm{\Pi }{\phi }_i)（27）

由于φ_i ∈V₁(Ω^e)，Πφ_i ∈P₁(Ω^e)，且P₁(Ω^e)⊂V₁(Ω^e)，则有投影残差(φ_i-Πφ_i )∈V₁(Ω^e)。式（18）可计算为：

\begin{array}{l}{k}_{ij}={\int }_{{\Omega }^{\mathrm{e}}}(\nabla (\mathrm{\Pi }{\phi }_i))D(\nabla (\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}+\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla (\mathrm{\Pi }{\phi }_i))D(\nabla ({\phi }_j-\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}+\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla ({\phi }_i-\mathrm{\Pi }{\phi }_i))D(\nabla (\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}+\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla ({\phi }_i-\mathrm{\Pi }{\phi }_i))D(\nabla ({\phi }_j-\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}\end{array}（28）

根据构造投影时所采用的正交条件式（19），可知式（27）右侧的第2、3项为0，式（28）进一步简化为：

\begin{array}{l}{k}_{ij}={\int }_{{\Omega }^{\mathrm{e}}}(\nabla (\mathrm{\Pi }{\phi }_i))D(\nabla (\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}+\\ {\int }_{{\Omega }^{\mathrm{e}}}(\nabla ({\phi }_i-\mathrm{\Pi }{\phi }_i))D(\nabla ({\phi }_j-\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}\end{array}（29）

根据式（26）、（29）右侧的第1项是可准确计算的，VEM称其为k_ij 的一致项；式（29）右侧的第2项在VEM称为k_ij 的稳定项，只能采用近似方式进行处理。

应用VEM求解各种不同的问题时，合适的稳定项是国内外学者研究的热点。文献[22]针对泊松问题的求解，使用投影残差在单元顶点处的值构建稳定项：

\begin{array}{l}{\int }_{{\Omega }^{\mathrm{e}}}(\nabla ({\phi }_i-\mathrm{\Pi }{\phi }_i))(\nabla ({\phi }_j-\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}\approx \\ {\displaystyle \sum _{m=\mathrm{1}}^n}(\lambda ({\phi }_i(x_m)-\mathrm{\Pi }{\phi }_i(x_m))({\phi }_j(x_m)-\mathrm{\Pi }{\phi }_j(x_m)))\end{array}（30）

式中：λ为用于保证稳定项与一致项数量级可比的近似乘子，文献[22]中λ取1；m为求和指标，遍历单元的所有顶点，m=1,2,…,n。

鉴于k_ij 由一致项和稳定项组成，为考虑渗透张量 D 的影响，将单元渗透矩阵 K^e分为与其同阶的两部分：

{\mathit{K}}^{\mathrm{e}}={\mathit{K}}^{\mathrm{e}\mathrm{c}}+\lambda {\mathit{K}}^{\mathrm{e}\mathrm{s}}（31）

式中：{\mathit{K}}^{\mathrm{e}\mathrm{s}}为组装稳定项对 K^e的贡献； K^ec为组装一致项对 K^e的贡献，其元素k_{ij}^{\mathrm{c}}为：

k_{ij}^{\mathrm{c}}={\int }_{{\Omega }^{\mathrm{e}}}(\nabla (\mathrm{\Pi }{\phi }_i))D(\nabla (\mathrm{\Pi }{\phi }_j))\mathrm{d}\mathit{x}（32）

K^es的元素k_{ij}^{\mathrm{s}}为：

k_{ij}^{\mathrm{s}}={\displaystyle \sum _{m=\mathrm{1}}^n}[({\phi }_i(x_m)-\mathrm{\Pi }{\phi }_i(x_m))({\phi }_j(x_m)-\mathrm{\Pi }{\phi }_j(x_m))]（33）

式（31）中λ的计算思路为：既然稳定项在VEM中用投影残差在单元顶点处的值估算，那么一致项也可尝试用投影在单元顶点处的值估算。按此种方式估算一致项后进行组装，将得到参考矩阵 K^refer，其元素为：

k_{ij}^{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}}={\displaystyle \sum _{m=\mathrm{1}}^n}[\mathrm{\Pi }{\phi }_i(x_m)\mathrm{\Pi }{\phi }_j(x_m)]（34）

K^refer与 K^es数量级应大致接近。要保障 K^es的数量级与 K^ec具有可比性，乘子λ可确定为：

\lambda =\frac{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}{\mathit{K}}^{\mathrm{e}\mathrm{c}}}{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}{\mathit{K}}^{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}}}（35）

式中，trace为对方阵求迹。

按照式（32）和（20），可计算任意多边形单元的单元渗透矩阵和单元流量向量，进而组装得到整体渗透矩阵和整体流量向量。三角形单元、矩形单元和四边形单元是经典FEM二维分析的常用单元，在FEM中其数学处理方法存在差异；但在VEM中，上述单元只是特定边数和形状的多边形单元，单元渗透矩阵和单元流量向量均按照式（31）和（19）执行。

2.4 1阶VEM与FEM单元的数学联系

1阶VEM采用单元顶点处物理量的值为自由度，当单元拥有不超过4条边时，其可与FEM中具备同样自由度的单元建立联系。具体而言，当单元为三角形时，1阶VEM可与FEM的3节点三角形单元关联；当单元为矩形时，1阶VEM与FEM的4节点矩形单元存在数学关联。

单元为三角形时，1阶VEM基函数φ_I 的投影Πφ_I 可具备Dirac‒Delta性质，其证明如图4所示。

根据式（25），可整理得基函数φ_I 的投影Πφ_I 为：

\mathrm{\Pi }{\phi }_I=\frac{\mathrm{1}}{\mathrm{3}}+\frac{l_I{\mathit{n}}_I^{\mathrm{T}}}{\mathrm{2}A}(\mathit{x}-{\mathit{x}}_{\mathrm{o}})+\frac{l_J{\mathit{n}}_J^{\mathrm{T}}}{\mathrm{2}A}(\mathit{x}-{\mathit{x}}_{\mathrm{o}})（36）

式中， x_o为单元中心。

根据矢量点乘的几何意义，当 x 位于顶点I处时，存在Πφ_I =1/3+1/3+1/3=1；当 x 位于顶点J处时，存在Πφ_I =1/3-2/3+1/3=0；当 x 位于顶点K处时，存在Πφ_I =1/3+1/3-2/3=0；因此，Πφ_I 具备Dirac‒Delta性质。类似地，可进行Πφ_J 和Πφ_K 的性质证明。

单元为三角形时，式（14）中的基函数及其投影均具有Dirac‒Delta性质，意味着稳定项为0，仅由一致项构成的 K^e与FEM 3节点三角形单元的结果将完全一致。以上特征存在的根本原因，是FEM 3节点三角形单元的域内物理量采用一次完备多项式近似，其测试函数空间P₁(Ω^e)为V₁(Ω^e)的子空间，同时也正好是V₁(Ω^e)的投影空间。

单元为矩形时，FEM 4节点矩形双线性单元对域内物理量用双线性近似，由式（10）可知双线性函数满足V₁(Ω^e)的全部规定。此测试函数空间为V₁(Ω^e)的子空间，但非V₁(Ω^e)的投影空间，因此，对矩形单元必须进行稳定项计算。当渗透系数非各向异性时，由式（30）计算得到的单元渗透矩阵与FEM 4节点矩形双线性单元具有良好的可比性。

值得注意的是，单元为任意四边形时，1阶VEM与FEM四边形4节点等参单元的自由度定义同样一致，但后者对域内物理量采用等参近似，测试函数不满足V₁(Ω^e)规定，因此两者之间不存在直接的数学联系。

3 求解流程

局部修改网格法求解饱和无压渗流模型时，求解流程一般为：1）对求解域进行固定网格剖分；2）施加边界条件，数值求解渗流场水头值，此时自由面仅施加流量边界条件；3）判断渗流自由面上水头值与其纵坐标的差是否小于收敛阈值，若已满足则迭代结束，否则根据求得的自由面水头值确定新自由面；4）根据新自由面局部修改网格；5）重复步骤2）～5），迭代计算。在执行步骤4）时，由于VEM适应于新自由面切割固定网格可能形成的多边形单元，因此局部修改网格可更直接。

3.1 局部修改网格的单元和节点更新

已有局部修改网格法与本文方法的区别如图5所示。图5中，数字为原节点编号，字母为与自由面相交的新的节点，黑色实线为固定网格，红色实线为本次迭代求出的自由面。以图5的固定网格和自由面为例，阐述本文局部修改网格与现有FEM方法的根本区别。

现有FEM方法在局部修改网格时，一般先根据单元边数是否符合FEM要求进行初步修改，再根据单元形状是否病态继续调整，以吴梦喜等^[19]的方法为例，区域16‒b‒c‒d‒14内最终形成3个单元。本文利用VEM对单元形状的包容性，直接计算自由面与固定网格的交点形成节点，区域16‒b‒c‒d‒14内最终形成4个单元，其中单元15‒8‒c′‒d‒14为FEM无法处理的五边形单元。

本文局部修改网格的实现过程：用自由面切割固定网格的第1步判断单元的类型，通过计算单元节点位置与自由面的距离关系实现。本文方法自由面切割单元的处理如图6所示，根据切割时需执行的操作，分为3种：1）删除单元，全部节点在自由面之上，设定操作类型号为0；2）原生单元，全部节点均在自由面之下，设定操作类型号为1；3）贯穿单元，部分节点在自由面之下，部分节点在自由面之上，设定操作类型号为2。

对贯穿单元进一步处理，计算自由面与贯穿单元的交点，这些交点与贯穿单元位于自由面以下的节点构成新生单元。以图6（b）为例，贯穿单元等同于删除单元加上新生单元，新生单元继承贯穿单元的操作类型号2，同时将原单元操作类型号修改为0，代表删除。以图6（a）中虚线框内的处理为例，6‒16‒15‒8为操作类型号1的原生单元，b‒6‒8‒c、c‒8‒c′和c′‒8‒15‒14‒d为操作类型号2的新生单元，5‒6‒8‒7、7‒8‒10‒9和8‒15‒14‒10为操作类型号0的删除单元。

根据切割时节点的操作方式，节点也分为删除节点、原生节点和新生节点3大类。以图6为例，3、5、7、9、10、11和12为操作类型号0的删除节点，2、4、6、8、13、14、15和16为操作类型号1的原生节点，a、b、c、c′、d和e为操作类型号2的新生节点。

以固定网格信息为基础，出现1个新生单元和节点则单元和节点的总数增加1，并记录新生单元的节点组成和新生节点的坐标等信息，最终形成新渗流区域网格的单元节点信息。虽然新的渗流区域网格上仅由原生单元和新生单元组成，但删除单元和删除节点的相关信息仍保存于固定网格信息库中。

值得注意的是，采用本文局部切割网格方法时，可能出现以图6中节点c′为代表的自由面新节点。因此，随着迭代求解的进行，本文自由面上节点数量可能逐步增加，从而对自由面进行节点加密处理。

3.2 整体渗透矩阵和整体流量向量的更新

为减少计算成本，本研究将固定网格的单元、节点信息及整体渗透矩阵 K 和整体流量向量 F 等全部存储。每次迭代计算确定新自由面之后，均用新自由面对固定网格进行切割执行网格局部修改，以便在计算新渗流区域网格的整体渗透矩阵和整体流量向量时，继承固定网格的诸多信息。

在新的渗流区域网格切割形成之后，考虑网格库中节点数量多于固定网格，根据网格库中节点总数将 K 和 F 分别扩充为 K ′和 F ′。接着，按照操作类型逐个单元计算分析，对 K ′和 F ′进行更新：操作类型号1的原生单元无须对 K ′和 F ′进行任何修改；操作类型号0的删除单元，计算其单元渗透矩阵和单元流量向量并从 K ′和 F ′中剔除；操作类型号2的新生单元，计算其单元渗透矩阵和单元流量向量并添加至 K ′和 F ′中。所有单元执行完毕则新渗流区域网格的整体渗透矩阵和整体流量向量计算完成。更新完毕之后，检查 K ′和 F ′中删除节点所对应的行列信息，正确情况下其应均为0，将删除节点对应的行列从 K ′和 F ′中删除。

综合上述内容，VEM迭代求解饱和无压渗流模型的流程如图7所示。收敛标准为自由面上节点处的水头值与纵坐标之差小于阈值。阈值取值根据问题尺寸规模设定，本研究中取0.001倍的问题分析区域纵向尺寸。

4 数值算例验证

编制本文算法Matlab程序，求解矩形均质土坝、梯形均质土坝和非均质矩形土坝等4个数值算例，并与其他算法结果比较，以验证本文算法的计算精度和计算效率。由于本文研究未涉及应力场的计算，因此，模型初始自由面设定为与上游水头齐平的水平线，以其为固定网格的边界可保证真正自由面处于固定网格范围内。为保证程序的一致性，固定网格的渗流计算和局部修改网格后的渗流计算均使用VEM进行。鉴于1阶VEM三角形单元（矩形单元）与FEM 3节点三角形单元（4节点矩形单元）的可比性，固定网格剖分仅使用三角形单元和矩形单元，以便与吴梦喜等^[19]基于FEM的局部修改网格法进行结果对比。需要说明的是，计算程序得到的自由面为多个节点相连的折线，后续图形展示中未做平滑处理。

4.1 算例1

矩形均质土坝算例为唯一有解析解可参考的算例^[30]，用其验证本文算法的计算精度。图8为矩形均质土坝模型及计算结果，上游水位H₁=1.0 m，下游水位H₂=0.5 m，坝体底部宽L=0.5 m，渗透系数k=1 m/d。对整个土坝区域用尺寸0.05 m×0.05 m的4节点矩形单元进行划分，初始自由面设定为坝顶水平线，采用本文方法进行求解，迭代收敛阈值为0.001 m。

由图8可知，7次迭代后的结果与解析解基本一致。迭代中，0.1、0.2、0.3、0.4和0.5 m这5个水平刻度处渗流自由面位置坐标如表1所示。

本文算法计算结果与解析解的绝对误差随迭代的变化趋势如图9所示。以解析解为参考，经历7次迭代后水平刻度0.5 m的出渗面处误差最大，绝对误差为0.003 m，相对误差为0.45%。总体上误差量级很小，说明本文算法具有良好的计算精度。

本文算法迭代求解过程中自由面节点会自然加密。此算例求解时初始自由面仅有节点11个，7次迭代结束时自由面节点数量增加为22个，最终获取的自由面折线与解析解曲线贴合程度良好。但是自由面节点数量的增加会在一定程度上影响计算效率。求解在便携式笔记本上执行（型号Legion Y70002021，CPU为i7‒11800H，内存16 G），第1次迭代计算耗时0.54 s，第7次迭代计算耗时0.79 s。

4.2 算例2

梯形均质土坝模型及计算结果如图10所示，上游水位H₁=5 m，下游水位H₂=1 m，坝体底宽L=7 m，渗透系数k为1 m/d。对整个土坝区域用尺寸为0.5 m×0.5 m的4节点矩形单元和3节点三角形单元（尺寸为矩形单元一半）进行划分，初始自由面设定为坝顶水平线，采用本文算法进行求解，迭代收敛阈值为0.005 m。

由图10可知，17次迭代后自由面趋于稳定。迭代中，1、2和3 m这3个水平刻度处及出渗点的渗流自由面位置坐标如表2所示。表2中，“\”表示前5次迭代过程中出渗点尚未达3 m水平刻度位置。必须说明的是，迭代求解过程中若不加干涉，自由面节点自然加密可能导致节点过于接近出现单元病态，因此，算法中当自由面上两个节点距离小于1/5的固定网格单元尺寸时，将两个节点合并处理。本例迭代结束时自由面节点数量为11个。

因无解析解可参照，用Seep软件、Zheng^[5]和吴梦喜^[19]等方法对比分析本文算法的计算精度。3种方法的模型和求解工具各异：Seep软件用饱和‒非饱和渗流分析模型，工具为FEM；Zheng等^[5]用饱和无压渗流模型，工具为改进的流形元；吴梦喜等^[19]用饱和无压渗流模型，工具为FEM。使用单元尺寸大致相同的网格进行求解，不同算法的自由面结果比较如图11所示。Seep计算时，材料模型为“饱和/不饱和”，使用数据点拟合创建体积含水量函数，饱和土水含量设置为0.5，水力传导率函数中设置渗透系数为1 m/d，其余参数皆采用默认设置。图11中，Seep的自由面位置最低，而吴梦喜等^[19]的结果最高，本文结果与Zheng等^[5]大致吻合，但出渗点位置略高。鉴于Zheng等^[5]算法的精度已得到充分检验^[7]，因此本文算法也具备足够计算精度。在计算精度相近的情况下，本文所采用的VEM与FEM高度兼容，同时避免了流形元方法中流形覆盖Zheng等^[5]复杂数学概念的引入。

本文算法与吴梦喜等^[19]方法同属局部修改网格法，比较两者的计算效率有助于理解本文算法的优势与不足。使用与第4.1节同样的电脑配置，吴梦喜等^[19]方法分析模型耗时13.6 s，本文算法耗时14.4 s，可见本文算法计算效率略有下降。因为本文算法用自由面切割固定网格时更为直接，相比吴梦喜等^[19]方法节省时间，但迭代过程中自由面节点加密和干预则相对耗费时间。考虑自由面节点加密有利于提高计算精度，且使自由面形态更为平滑，本文算法计算效率相比吴梦喜等^[19]传统局部修改网格法的略微下降可以接受。

4.3 算例3

非均质矩形土坝模型及计算结果如图12所示。使用图12的非均质矩形土坝算例对本文算法进行考核。该土坝上游水位H₁=60 m，下游水位H₂=20 m，坝体底部宽L=100 m，渗透系数在AB线处发生改变，左侧的渗透系数k₁为1 m/d，右侧的渗透系数k₂为10 m/d。对土坝区域用尺寸为10 m×10 m的4节点矩形单元进行划分，初始自由面设定为坝顶水平线，采用本文算法进行求解，迭代收敛阈值为0.06 m。

由图12可知，7次迭代后自由面趋于稳定。迭代中20、50、80和100 m等4个水平刻度处的渗流自由面位置坐标如表3所示；渗透系数变化处自由面呈迅速下降形态，在AB线右侧10 m范围内自由面呈现出轻微的下凸特征，与Zheng等^[5]的结论大致吻合。本算例求解时，对自由面节点自然加密过程同样施加干预，迭代结束时自由面节点数量为21个。

同样采用Seep软件、Zheng^[5]和吴梦喜^[19]等方法分析本算例，使用单元尺寸大致相同的网格进行求解，不同算法的自由面结果比较如图13所示。Seep计算时，体积含水量函数创建方式与算例2相同，水力传导率函数中分别将左右两区域的渗透系数设置为1 m/d和10 m/d。图13中，右侧出渗点附近自由面无明显差别，左侧区域和渗透系数变化处的自由面则存在差异。首先，按自由面位置的高低排序，吴梦喜等^[19]最高，本文算法、Zheng等^[5]和Seep软件结果逐次下降。其次，除吴梦喜等^[19]之外，自由面在接近渗透系数突变线时均迅速跌落，越过突变线后则相对波动较小。本文算法结果与Zheng等^[5]的结果吻合程度较好，进一步证实了本文算法的计算精度。

计算效率方面，本文算法与吴梦喜等^[19]方法相比计算耗时仍略有增长。使用与第4.1节同样的电脑配置，吴梦喜等^[19]方法分析模型耗时7.3 s，本文算法耗时8.1 s。

5 结 论

本文引入虚拟单元法求解饱和无压渗流模型，构建了1阶VEM求解此问题的计算格式，提出相应的网格切割方法和整体渗透矩阵、整体流量向量更新算法，编制求解程序，并通过算例验证了本文提出方法的有效性。相关结论如下：

1）采取VEM可有效解决局部修改固定网格方法求解饱和无压渗流模型时因单元形状受FEM限制需精细处理的局限。

2）本文提出的网格切割方法和整体渗透矩阵、整体流量向量更新算法，可有效继承固定网格的计算信息，从而节省计算成本。

3）数值算例结果表明本文所提方法可有效地求解饱和无压渗流模型，计算结果与解析解或改进流形单元法的结果具有很好的可比性，具备足够的计算精度。

4）迭代过程中自由面节点的加密，有利于提高计算精度，且使自由面形态更为平滑，但导致与传统局部修改固定网格方法相比计算效率略微下降。

参考文献

[1] Chen Yifeng, Zhou Chuangbing, Hu Ran,et al.Key issues on seepage flow analysis in large scale hydropower projects[J].Chinese Journal of Geotechnical Engineering,2010,32(9):1448‒1454.

[2] 陈益峰,周创兵,胡冉,等.大型水电工程渗流分析的若干关键问题研究[J].岩土工程学报,2010,32(9):1448‒1454.

[3] 郑宏.数值流形法[M].北京:科学出版社,2022.

[4] Zhao Lanhao, Zhang Hairong, Mao Jia,et al.An ICLS-based method for solving two-phase seepage free surface considering compressible gas in porous media[J].Computers and Ge-otechnics,2022,141:104528. doi:10.1016/j.compgeo.2021.104528

[5] Wei Wei, Jiang Qinghui, Ye Zuyang,et al.Equivalent fracture network model for steady seepage problems with free surfaces[J].Journal of Hydrology,2021,603:127156. doi:10.1016/j.jhydrol.2021.127156

[6] Zheng Hong, Liu Feng, Li Chunguang.Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method[J].Applied Mathematical Modelling,2015,39(2):794‒808. doi:10.1016/j.apm.2014.07.007

[7] Jia Zhen, Zheng Hong.A new procedure for locating free surfaces of complex unconfined seepage problems using fixed meshes[J].Computers and Geotechnics,2024,166:106032. doi:10.1016/j.compgeo.2023.106032

[8] Li Xilong, Zhang Hong.Analyzing unconfined seepage fl-ow with corner singularity using an enhanced second-ord-er numerical manifold method[J].Computers and Geotechnics,2024,167:106101. doi:10.1016/j.compgeo.2024.106101

[9] Yuxin Jie, Liu Lizhen, Xu Wenjie,et al.Application of NEM in seepage analysis with a free surface[J].Mathematics and Computers in Simulation,2013,89:23‒37. doi:10.1016/j.matcom.2013.03.006

[10] Kazemzadeh‒Parsi M J.Isogeometric analysis in solution of unconfined seepage problems[J].Computers & Mathe-matics with Applications,2019,78(1):66‒80. doi:10.1016/j.camwa.2019.02.011

[11] Wang Zhaoqing, Li Shucai, Li Shuchen.Multi-node finite element approaches for unconfined seepage flow problem[J].Rock and Soil Mechanics,2008,29(10):2647‒2650. doi:10.3969/j.issn.1000-7598.2008.10.010

[12] 王兆清,李术才,李树忱.无压渗流问题分析的多节点有限元方法[J].岩土力学,2008,29(10):2647‒2650. doi:10.3969/j.issn.1000-7598.2008.10.010

[13] Dai Qianwei, Lei Yi, Zhang Bin,et al.A practical adaptive moving-mesh algorithm for solving unconfined seepage problem with Galerkin finite element method[J].Scientific Reports,2019,9:6988. doi:10.1038/s41598-019-43391-4

[14] Zheng Hong, Dai Huichao, Liu Defu.Improved Bathe's algorithm for seepage problems with free surfaces[J].Rock and Soil Mechanics,2005,26(4):505‒512. doi:10.3969/j.issn.1000-7598.2005.04.001

[15] 郑宏,戴会超,刘德富.改进的有自由面渗流问题的Bathe算法[J].岩土力学,2005,26(4):505‒512. doi:10.3969/j.issn.1000-7598.2005.04.001

[16] Chen Yifeng, Hu Ran, Zhou Chuangbing,et al.A new parabolic variational inequality formulation of Signorini's condition for non-steady seepage problems with complex see-page control systems[J].International Journal for Numerical and Analytical Methods in Geomechanics,2011,35(9):1034‒1058. doi:10.1002/nag.944

[17] Luo Guanyong, Pan Hong.Using Bathe algorithm and Sig-norini condition to solve unconfined saturated‒unsaturated seepage problems[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(11):2275‒2282. doi:10.3969/j.issn.1000-6915.2013.11.013

[18] 骆冠勇,潘泓.结合Bathe算法及Signorini条件求解饱和‒非饱和无压渗流问题[J].岩石力学与工程学报,2013,32(11):2275‒2282. doi:10.3969/j.issn.1000-6915.2013.11.013

[19] Alnashri Y, Droniou J.A gradient discretization method to analyze numerical schemes for nonlinear variational inequ-alities,application to the seepage problem[J].SIAM Journal on Numerical Analysis,2018,56(4):2375‒2405. doi:10.1137/16m1105517

[20] Pan Shulai, Wang Quanfeng, Yu Jin.Improvement of analysis of free surface seepage problem by using initial flow method[J].Chinese Journal of Geotechnical Engineering,2012,34(2):202‒209.

[21] 潘树来,王全凤,俞缙.利用初流量法分析有自由面渗流问题之改进[J].岩土工程学报,2012,34(2):202‒209.

[22] Zhou Bin, Yan Jun, Liu Sihong,et al.Theoretical interpretation of nodal virtual flux method and its optimized algori-thm[J].Rock and Soil Mechanics,2018,39(1):349‒355.

[23] 周斌,严俊,刘斯宏,等.结点虚流量法理论基础阐释及改进算法[J].岩土力学,2018,39(1):349‒355.

[24] Fu Yanling, Zhou Zhifang, Wu Yongxia.Improved adjustm-ent method of compound element conductivity matrix for calculating 3D seepage field with free surface[J].Chinese Journal of Geotechnical Engineering,2009,31(9):1434‒1439.

[25] 付延玲,周志芳,武永霞.改进复合单元渗透矩阵调整法求解自由面三维渗流场[J].岩土工程学报,2009,31(9):1434‒1439.

[26] Wu Mengxi, Zhang Xueqin.Imaginary element method for numerical analysis of seepage with free surface[J].Journal of Hydraulic Engineering,1994(8):67‒71.

[27] 吴梦喜,张学勤.有自由面渗流分析的虚单元法[J].水利学报,1994(8):67‒71.

[28] Liang Yeguo, Xiong Wenlin, Zhou Chuangbing.Subeleme-nt method for seepage analysis with free surface[J].Journ-al of Hydraulic Engineering,1997,28(8):34‒38.

[29] 梁业国,熊文林,周创兵.有自由面渗流分析的子单元法[J].水利学报,1997,28(8):34‒38.

[30] Beirão da Veiga L, Brezzi F, Cangiani A,et al.Basic principles of virtual element methods[J].Mathematical Models and Methods in Applied Sciences,2013,23(1):199‒214. doi:10.1142/s0218202512500492

[31] Beirão da Veiga L, Brezzi F, Marini L D,et al.The hitchhiker's guide to the virtual element method[J].Mathematical Models and Methods in Applied Sciences,2014,24(8):1541‒1573. doi:10.1142/s021820251440003x

[32] Benvenuti E, Chiozzi A, Manzini G,et al.Extended virtual element method for two-dimensional linear elastic fracture[J].Computer Methods in Applied Mechanics and Engine-ering,2022,390:114352. doi:10.1016/j.cma.2021.114352

[33] Marfia S, Monaldo E, Sacco E.Cohesive fracture evolution within virtual element method[J].Engineering Fracture Mechanics,2022,269:108464. doi:10.1016/j.engfracmech.2022.108464

[34] Lin Shan, Yang Yongtao, Sun Guanhua,et al.Elastoplastic mechanical analysis based on the virtual element method[J].Chinese Journal of Solid Mechanics,2020,41(1):30‒40. doi:10.19636/j.cnki.cjsm42-1250/o3.2019.042

[35] 林姗,杨永涛,孙冠华,等.弹塑性力学问题的虚单元法[J].固体力学学报,2020,41(1):30‒40. doi:10.19636/j.cnki.cjsm42-1250/o3.2019.042

[36] Lin Shan, Guo Yukui, Sun Guanhua,et al.Virtual element st-rength reduction method for slope stability analysis[J].Ch-inese Journal of Rock Mechanics and Engineering,2019,38(S2):3429‒3438.

[37] 林姗,郭昱葵,孙冠华,等.边坡稳定性分析的虚单元强度折减法[J].岩石力学与工程学报,2019,38(S2):3429‒3438.

[38] Jiang Wei, Xu Jiancheng, Wang Lehua,et al.A novel formulation of discontinuous deformation analysis enlightened by virtual element method[J].Chinese Journal of Rock Mechanics and Engineering,2022,41(1):106‒119.

[39] 江巍,徐建城,王乐华,等.基于虚单元法的非连续变形分析方法新格式[J].岩石力学与工程学报,2022,41(1):106‒119.

[40] Jiang Wei, Yin Hao, Wu Jian,et al.Virtual element method for solving 2d geometric nonlinear problems upon S‒R decomposition theorem[J].Engineering Mechanics,2024,41(8):23‒35.

[41] 江巍,尹豪,吴剑,等.基于S‒R和分解定理的二维几何非线性问题的虚单元法求解[J].工程力学,2024,41(8):23‒35.

[42] Liu Chuanqi, Xu Guangtao, Wei Yujie.Virtual element met-hod:Theory and applications[J].Advances in Mechanics,2022,52(4):874‒913. doi:10.6052/1000-0992-22-037

[43] 刘传奇,许广涛,魏宇杰.虚单元计算方法的最新理论与应用进展[J].力学进展,2022,52(4):874‒913. doi:10.6052/1000-0992-22-037

[44] Liu Fangxue, Lei Guohui, Wang Weiyu,et al.Charts for free surfaces in steady-state seepage flow through homogen-eous isotropic rectangular dams[J].Journal of Hydrology,2022,612:128082. doi:10.1016/j.jhydrol.2022.128082

基金资助

国家自然科学基金项目(U22A20600)

国家自然科学基金项目(52079070)

国家自然科学基金项目(52079002)

防灾减灾湖北省重点实验室开放基金项目(2022KJZ07)