基于弱有限元法的渗流分析

何朝葵; 速宝玉; 盛金昌; 罗玉龙

doi:10.11779/CJGE20220635

基于弱有限元法的渗流分析 English Version

1.
河海大学理学院, 江苏南京 211100
2.
河海大学水利水电学院, 江苏南京 210024

详细信息

作者简介:
何朝葵(1974—)，男，硕士，讲师，主要从事有限元方法方面的教学和科研工作。E-mail: hzk@hhu.edu.cn

中图分类号: TU139.14
计量
- 文章访问数: 301
- HTML全文浏览量: 61
- PDF下载量: 86
出版历程
- 收稿日期: 2022-05-16
- 网络出版日期: 2023-02-23
- 刊出日期: 2023-06-30

Seepage analysis based on weak finite element method

1.
College of Science, Hohai University, Nanjing 211100, China
2.
College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210024, China

摘要

摘要: 在坝、闸和堤防等水工建筑物的设计和运行管理中，渗流计算占有非常重要的地位。目前渗流分析的方法以一般有限元方法为主。在渗流分析时往往希望通过加密网格来追求更高精度的数值结果，但是由于一般有限元法属于协调有限元方法，它不适合处理混合剖分网格，加密网格容易导致计算格式不稳定、计算结果不收敛。而弱有限元方法是一种非协调有限元方法，用弱函数的弱梯度算子替代一般有限元方法变分式中的经典梯度算子，并在变分式中增加稳定子项，保证了计算格式的绝对稳定性，并且它适用于混合网格。对于同一剖分网格，弱有限元法的总体代数方程组的自由度远远比一般有限元方法的要大很多，为了降低弱有限元法的自由度，探讨根据弱有限元方法适用于混合剖分网格的特点，采用有针对性的混合剖分网格，运用弱有限元法求解渗流自由面和进行有防渗帷幕的闸基渗流分析，数值结果表明弱有限元方法能灵活处理混合剖分网格，并且计算结果具有高精度。
- 弱有限元方法 /
- 混合剖分网格 /
- 渗流 /
- 自由面 /
- 闸基
Abstract: In the design and operation management of hydraulic structures such as dams, gates and embankments, the seepage calculation plays a very important role. The current methods for seepage analysis are mainly based on the general finite element method. It is often hoped, by refining the mesh, to pursue the calculated results with higher-precision in the seepage analysis. However, the general finite element method, which belongs to the coordinated finite element method, cannot handle hybrid grids, and the refined mesh may lead to unstable calculation formats and the divergence of the calculated results. The weak finite element method is a non-coordinated one, which replaces the classical gradient operator in the variational formula for the general finite element method with the weak gradient operator of the weak function, and adds a stable sub-term in the variational formula, to obtain an absolute stability of the calculation format, and it can handle hybrid grids. For the same mesh grid, the degree of freedom of the overall algebraic equations for the weak finite element method is much larger than that of the general finite element method. In order to reduce the degree of freedom of the weak finite element method, the targeted mixed meshing grid is used to establish the weak finite element method which is used to solve the free surface of seepage and to analyze the seepage field of the gate foundation with anti-seepage curtain. The numerical simulation shows that the weak finite element method can be used to process flexibly hybrid grids, and it is of high accuracy.
- weak finite element method /
- hybrid grid /
- seepage /
- free surface /
- brake base

HTML全文

0. 引言

渗流计算是评价水工建筑物渗流安全及进行渗控方案设计的依据。渗流计算的实践表明有限元方法是一种非常成熟、非常有效的数值计算方法。有限元方法分为协调有限元方法和非协调有限元方法，要求近似函数在剖分单元边界上连续的称为协调有限元方法，允许近似多项式函数在单元边界上不连续的称为非协调有限元方法，又称为间断有限元方法。通常讲的有限元方法是指协调有限元方法，称为一般有限元方法。一般有限元方法连续性的要求，使得剖分网格只能使用单一的单元类型，同一次数的近似多项式函数，因而一般有限元方法不适用于混合剖分网格。在渗流计算中，往往希望通过加密网格的方法来提高计算精度，而这有可能导致计算格式不稳定、计算结果不收敛。国内外众多学者一直在探索适用于混合剖分网格的间断有限元方法，主要成果有：内罚间断有限元法（IPDG）^[1]，局部间断有限元法（LDG）^[2]和可杂交的间断有限元法（HDG）^[3]。与一般有限元法比较，非协调有限元方法有两方面的优势：①近似函数灵活，可以在不同剖分单元采用不同次数的近似多项式函数；②网格生成灵活，适用混合剖分网格。但是，间断有限元方法的变分形式比协调有限元方法的变分形式要复杂得多，计算格式的稳定性往往需要靠调整参数来保证、而辅助变量的引入和变分的复杂化则使得单元矩阵非局部化，这些困难使得间断有限元方法难以在渗流分析中充分发挥其优势^[4]。弱有限元法（WG）是Wang等^[5]提出的，它属于非协调有限元方法，通过在剖分单元上定义弱函数和弱梯度，从而用弱梯度算子替代经典的梯度算子，并在计算格式中引入了一个不含有参数的稳定子或光滑子，保证了格式的绝对稳定性^[5-7]。因此，弱有限元方法不但具有间断有限元方法的优点，还克服了间断有限元方法的缺点，目前，弱有限元方法可以用来求解Stokes方程、Helmholtz方程、Maxwell方程、线弹性方程等^[8-12]。

本文主要研究弱有限元法针对混合剖分网格在渗流分析上的具体应用。安排如下：在第一节中将介绍弱有限元方法的原理，并通过一个简例介绍弱有限元方法处理混合网格的流程；在第二节中运用弱有限元方法，采用混合剖分网格求解均值土坝渗流自由面和闸基渗流场。

1. 弱有限元方法

1.1 变分形式

稳定渗流的数学模型为椭圆型方程边值问题：

$\left\{ {\begin{array}{*{20}{c}} { - \nabla (p\nabla u) = f}&{{\text{in }}\Omega } \\ {u = c}&{{\text{on }}{\Gamma _1}} \\ {\frac{{\partial u}}{{\partial {\bf{n}}}} = g}&{{\text{on }}{\Gamma _2}} \end{array}} \right. \text{，}$

(1)

式中： ${\Gamma _1} \cup {\Gamma _1}_2 = \partial \Omega$ ，n为 ${\Gamma _2}$ 上的外法线单位向量；u为水头函数，p为渗透系数矩阵。此问题的一般有限元法的变分形式为：求 $u \in {H^1}(\Omega )$ ，使得在 ${\Gamma _1}$ 上 $u = c$ ，且满足

$a(u, v) = (f, v) + {\left\langle {g, v} \right\rangle _{{\Gamma _2}}}{\text{ }}\forall v \in H_0^1(\Omega ) 。$

式中， $a(u, v) = (p\nabla u, \nabla v)$ 。在弱有限元法中，允许近似函数在单元边界上不连续，为此引入弱函数和弱函数的梯度算子。假设 ${\mathcal{T}_{\text{h}}}$ 是区域 $\Omega$ 的一个剖分，剖分单元可以是三角形、四边形或其它多边形。单元T上的弱函数v定义为

$v = \left\{ {{v_0}, {v_{\text{b}}}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{v_0}}&{{\text{in }}T} \\ {{v_{\text{b}}}}&{{\text{on }}\partial T} \end{array}} \right. 。$

即，弱函数 $v$ 在单元 $T$ 上包含内部函数 ${v_0}$ 和边界函数 ${v_{\text{b}}}$ 两部分。用 ${P_m}(T)$ 表示定义在T上且次数不超过m的多项式函数，弱函数v的弱梯度算子 ${\nabla _w}v$ 定义为对 $\forall {\bf{ \pmb{\mathit{ τ}} }} \in {[{P_m}(T)]^2} = \left( {{P_m}(T), {P_m}(T)} \right)$ ，有

$\left( {{\nabla _w}v, {\bf{ \pmb{\mathit{ τ}} }}} \right) = - {\left( {{v_0}, \nabla \cdot {\bf{ \pmb{\mathit{ τ}} }}} \right)_T} + {\left\langle {{v_b}, {\bf{ \pmb{\mathit{ τ}} }} \cdot {\bf{n}}} \right\rangle _{\partial T}} 。$

(2)

单元 $T$ 上的弱有限元空间定义为

${W}_{j, l}(T)=\{v=\left\{{v}_{0}, {v}_{\text{b}}\right\}|{v}_{0}(T)\in {P}_{j}(T), {v}_{\text{b}}(e)\in {P}_{l}(e)\text{，} \left. {e \in \partial T} \right\} \text{，}$

弱梯度算子是由 ${W_{j, l}}(T)$ 到 ${[{P_m}(T)]^2}$ 上的一个映射。剖分 ${\mathcal{T}_{\text{h}}}$ 上的弱有限元法空间定义为

${V_{\text{h}}} = \left\{ {\left. v \right|v(T) \in {W_{j, l}}(T), \forall T \in {\mathcal{T}_h}} \right\} 。$

令

$V_{\text{h}}^{\text{0}} = \left\{ {v = \{ {v_0}, {v_{\text{b}}}\} \in {V_{\text{h}}}, {{\left. {{v_{\text{b}}}} \right|}_{{\Gamma _1}}} = 0, } \right\} \text{，}$

问题(1)的弱有限元法的变分形式为，求 ${u_{\text{h}}} \in {V_{\text{h}}}$ ，使得在 ${\Gamma _1}$ 上 ${u_{\text{h}}} = c$ ，且满足

${a_{\text{h}}}({u_{\text{h}}}, v) + s({u_{\text{h}}}, v) = (f, v) + {\left\langle {g, v} \right\rangle _{{\Gamma _2}}}, {\text{ }}\forall v \in V_{\text{h}}^{\text{0}} 。$

(3)

其中 ${a_{\text{h}}}({u_{\text{h}}}, v) = (p{\nabla _{\text{w}}}{u_{\text{h}}}, {\nabla _{\text{w}}}v)$ ， $s({u_{\text{h}}}, v)$ 是为了保证格式稳定的一个补充项，称为稳定子或光滑子，且

$s({u_{\text{h}}}, v) = \sum\limits_{T \in {\mathcal{T}_{\text{h}}}} {{{\left( {{J_\partial }({u_{\text{h}}}), {J_\partial }(v)} \right)}_{\partial T}}} 。$

(4)

式中， ${J_\partial }(v) = h_{\text{T}}^{ - 1/2}(Q({v_0}) - {v_{\text{b}}})$ ， ${h_{\text{T}}}$ 表示单元T的直径，即T上任意两点间距离的最大值（编程计算时一般取为多边形单元任意两节点间距离的最大值）， $Q({v_0})$ 表示内部函数 ${v_0}$ 在边界上的投影。 $\Theta (\nabla {u_0})$ 表示 $\nabla {u_0}$ 在单元边界 $\partial T$ 上的投影。根据弱梯度的定义和分部积分公式计算得

${({\nabla _{\text{w}}}{u_{\text{h}}}, {\nabla _{\text{w}}}v)_{\text{T}}} = {(\nabla {u_0}, \nabla {v_0})_{\text{T}}} - {\left\langle {\Theta (\nabla {u_0}) \cdot {\bf{n}}, Q({v_0}) - {v_{\text{b}}}} \right\rangle _{\partial T}} - {\left\langle {Q({u_0}) - {u_{\text{b}}}, {\nabla _{\text{w}}}v \cdot {\bf{n}}} \right\rangle _{\partial T}} \text{，}$

注意到

$\Theta (\nabla {u}_{0})\cdot n\approx (Q({u}_{0})-{u}_{\text{b}})/{h}_{\text{T}}, \text{ }{\nabla }_{\text{w}}v\cdot n\approx (Q({v}_{0})-{v}_{\text{b}})/{h}_{\text{T}}，$

故

${(\nabla {u}_{0}, \nabla {v}_{0})}_{\text{T}}\approx {({\nabla }_{\text{w}}{u}_{\text{h}}, {\nabla }_{\text{w}}v)}_{\text{T}}+2{h}_{\text{T}}^{-1}{\langle Q({u}_{0})-{u}_{\text{b}}, Q({v}_{0})-{v}_{\text{b}}\rangle }_{\partial T}。$

因此在变分式中用 ${a_{\text{h}}}({u_h}, v)$ 近似 $a(u, v)$ ，把 $s({u_{\text{h}}}, v)$ 作为补偿项加入变分式中是合理的^[7]。

1.2 计算格式

弱有限元空间 ${W_{j, l}}(T)$ 表示在单元T内部近似函数采用 $j$ 次多项式，在单元边界 $\partial T$ 上近似函数采用 $l$ 次多项式。 ${[{P_m}(T)]^2}$ 表示弱函数的弱梯度是分量为m次多项式的向量。选取不同的 $j, l$ 和m，会得到单元局部矩阵不同的计算格式， ${W_{j, l}}(T)$ 和 ${[{P_m}(T)]^2}$ 简记为 ${W_{j, l}}(T)/{[{P_m}(T)]^2}$ ，下面假设单元T是一个三角形单元（图 1）。

图 1 三角形单元图

Figure 1. Triangular element

下载: 全尺寸图片幻灯片

若选取弱有限元空间 ${W_{1, 0}}(T)/{[{P_0}(T)]^2}$ ，取基函数为

${\varphi _1} = \left\{ {\begin{array}{*{20}{c}} {1, }&{{\text{in }}T, } \\ {0, }&{{\text{on }}\partial T.} \end{array}} \right. \text{，} {\varphi _2} = \left\{ {\begin{array}{*{20}{c}} {x, }&{{\text{in }}T, } \\ {0, }&{{\text{on }}\partial T.} \end{array}} \right. \text{，} {\varphi _3} = \left\{ {\begin{array}{*{20}{c}} {y, }&{{\text{in }}T, } \\ {0, }&{{\text{on }}\partial T.} \end{array}} \right. \text{，}$

${\phi }_{3+i}=\{\begin{array}{cc}0, & 其它, \\ 1, & \text{on }{e}_{i}.\end{array}\text{ (}i=1, 2, 3)\text{ }。$

则，

${W_{j, l}}(T) = {\text{span}}\left\{ {{\varphi _1}, {\varphi _2}, \cdots , {\varphi _6}} \right\} \text{，}$

${u_{\text{h}}}(T) = \left\{ {{u_{{\text{h0}}}}(T), {u_{{\text{hb}}}}(T)} \right\} = \sum\limits_{i = 1}^6 {{a_{Ti}}{\varphi _i}} 。$

变分式（3）右端项，

${(f, v)_{\text{T}}} = {(f, {v_0})_{\text{T}}} \text{，}$

(5)

${\left\langle {g, v} \right\rangle _{\partial T \cap {\Gamma _2}}} = \sum\limits_{e \in \partial T \cap {\Gamma _2}}^{} {{{\left\langle {g, {v_{\text{b}}}} \right\rangle }_e}} \text{，}$

(6)

依次取 $v = {\varphi _i}, i = 1, 2, \cdots , 6$ ，由式（5）和（6）可得单元常数项矩阵.单元局部系数矩阵A_T包含两部分，

$\boldsymbol{A}_T=\boldsymbol{M}_T+\boldsymbol{S}_T。$

(7)

其中，

$\boldsymbol{M}_T=\left(\nabla_{\mathrm{w}} \varphi_i, \nabla_{\mathrm{w}} \varphi_j\right)_{6 \times 6}，$

${\mathit{\boldsymbol{S}}_T} = {\left( {\sum\limits_{i = 1}^3 {{{\left\langle {{J_\partial }\left( {{\varphi _i}} \right),{J_\partial }\left( {{\varphi _j}} \right)} \right\rangle }_{{e_i}}}} } \right)_{6 \times 6}}{\rm{ }}。$

根据弱梯度的定义（2），由 $\boldsymbol{\tau} \in\left[P_0(T)\right]^2$ ，有 ${\nabla _{\text{w}}}{\varphi _i} =$ 0， $\nabla_{\mathrm{w}} \varphi_{3+i}=\frac{\left|e_i\right|}{|T|} \boldsymbol{n}_i, i=1,2,3$ ，其中， $|{e_i}|$ 为单元边界 ${e_i}$ 的长度， $|T|$ 为单元T的面积，n_i为单元边界 ${e_i}$ 的外法线单位向量。令

$\boldsymbol{M}_T=\left(\begin{array}{cc} O & O \\ O & \boldsymbol{M}_{b b} \end{array}\right)$

(8)

这里 $M_{b b}=\left(\frac{\left|e_i \| e_j\right|}{|T|} \boldsymbol{n}_i \boldsymbol{n}_j\right)_{3 \times 3}$ 。S_T通过数值积分计算，用 $({x_i}, {y_i}), i = 1, 2, 3$ 表示单元边界 ${e_i}$ 的中点坐标，则

${\mathit{\boldsymbol{S}}_T} = \sum\limits_{i = 1}^3 {\left| {{e_i}} \right|\mathit{\boldsymbol{N}}_i^{\rm{T}}{\mathit{\boldsymbol{N}}_i}}$

(9)

其中

$\boldsymbol{N}_1=h_T^{-1 / 2}\left(1, x_1, y_1,-1,0,0\right) \text {, }$

$\boldsymbol{N}_2=h_T^{-1 / 2}\left(1, x_2, y_2, 0,-1,0\right) \text {, }$

$\boldsymbol{N}_3=h_T^{-1 / 2}\left(1, x_3, y_3, 0,0,-1\right) 。$

式（5）~（9）称为弱有限元法的计算格式。对于四边形单元的计算格式参考文献[6]，其他多边形的计算格式按照这个过程也易于推得。

由于相邻单元内部函数之间没有联系，所以可以选取不同次数的近似多项式函数，因而弱有限元方法适用于混合剖分网格，对于边中点网格，如剖分图 2中第4号单元是个四边形，一条边中间有个节点，那么这条边上的近似函数可以理解为分段近似，这在局部矩阵计算格式上相当于两条边，这个点把该边一分为二，这个单元就可以被认定为五边形。类推，如果一条边中间有多个节点，就可把该条边分成多条边，该单元就被认定为相应类型的多边形单元。因此，对于弱有限元方法，网格剖分比较自由。

图 2 剖分网格图

Figure 2. Diagram of partial grids

下载: 全尺寸图片幻灯片

1.3 混合剖分网格应用简例

下面通过一个简单的例子来了解弱有限元法处理混合网格的计算过程和计算效果。

实例1：考虑

$\left\{\begin{array}{rc} -\nabla \cdot(\nabla u)=0, & (x, y) \in[0,2] \times[0,2] \\ u=1, & x=0 \\ u=0, & x=2 \\ \frac{\partial u}{\partial \boldsymbol{n}}=0, & y=0 \text { 或 2。 } \end{array}\right.$

(10)

对正方形区域 $\Omega = [0, 2]$ × $[0, 2]$ 做剖分，如图 2所示。区域被分成4个单元，其中第1、2号单元为三角形，第3号单元为四边形，第4号单元有一条边中间有一个节点，因而该单元被认定为五边形，节点共8个，单元边界共11条。单元信息见表 1，各节点坐标信息见表 2。

表 1 剖分单元信息

Table 1. Information of elements

单元编号	单元类型	节点编号
1	3	3	2	4
2	3	3	1	2
3	4	3	4	5	6
4	5	3	6	7	8	1

下载: 导出CSV

| 显示表格

表 2 节点坐标信息

Table 2. Information of node coordinates

节点	1	2	3	4	5	6	7	8
x	1	2	1	2	2	1	0	0
y	0	0	1	1	2	2	2	0

下载: 导出CSV

| 显示表格

每个单元的弱有限元空间取 ${W_{1, 0}}(T)/{[{P_0}(T)]^2}$ ，根据式（7）~（9）计算各单元的局部系数矩阵A_T，用 $\mathit{\boldsymbol{A}}_T^i$ 表示第 $i$ 号单元的局部系数矩阵，则 $\mathit{\boldsymbol{A}}_T^1，\mathit{\boldsymbol{A}}_T^2$ 为6阶方阵， $\mathit{\boldsymbol{A}}_T^3，\mathit{\boldsymbol{A}}_T^4$ 分别为7阶、8阶方阵。为了把单元局部矩阵组合成整体矩阵，需要把所有自由度排序，每个单元的内部近似函数有3个自由度，每条单元边界对应一个自由度，为此先给所有单元边界编号如下（图 3）。

图 3 单元边界编号

Figure 3. Number of element boundary

下载: 全尺寸图片幻灯片

如果把所有单元内部近似函数的自由度按单元序号排前面，即，第 $i$ 号单元内部近似函数的3个自由度序号依次为 $3(i - 1) + 1$ ， $3(i - 1) + 2$ ， $3(i - 1) + 3$ ，把边界自由度按序号排在所有单元内部近似函数自由度序号的后面，比如，第1号边界的自由度排在第3×4+1=13号，第（8）号边界的自由度排在第3×4+8=20号。对于这个剖分网格，一共有23个自由度。依据自由度序号把单元局部矩阵组合成整体系数矩阵A。例如，第2号单元的局部矩阵的元素 $\mathit{\boldsymbol{A }}_T^2(3, 6)$ ，它是第2号单元上的基函数 ${\varphi _3}$ 和 ${\varphi _6}$ 作用的结果，它们对应的自由度序号分别为6和13，所以应该累加到总体系数矩阵的A（6，13），把所有局部单元矩阵组合成整体系数矩阵A，A是一个23阶的方阵，此问题的总体代数方程组 $\mathit{\boldsymbol{A X}}{\rm{ = }}\mathit{\boldsymbol{b}}$ ，由于此例 $f = 0$ ，且在第二类边界外法线方向上的方向导数为零，所以此时常数项矩阵b是一个23×1阶的零矩阵，未知数X是由所有自由度构成的一个23×1阶矩阵。问题的第二类边界条件属于自然边界条件，已在变分式中体现了，接下来需要处理第一类边界条件。处理第一类边界条件有不同的方法，这里采用如下处理方法：由于第2号、第6号和第10号单元边界属于第一类边界，它们对应的自由度序号分别为14，18，22，第2号和第6号边界上函数值为0，第10号边界上函数值为1，故对总体系数矩阵和常数项矩阵做如下修正：

$\begin{aligned} & \boldsymbol{A}(14, :)=0, \boldsymbol{A}(14, 14)=1, \boldsymbol{b}(14)=0, \\ & \boldsymbol{A}(18, :)=0, \boldsymbol{A}(18, 18)=1, \boldsymbol{b}(18)=0, \\ & \boldsymbol{A}(22, :)=0, \boldsymbol{A}(22, 22)=1, \boldsymbol{b}(22)=1 \text { 。 } \end{aligned}$

其中， $\mathit{\boldsymbol{A}}(i, :)$ 表示矩阵A的第i行的所有元素。处理边界条件后就得到该问题的弱有限元法的最终的总体线性方程组 $\mathit{\boldsymbol{A X}}{\rm{ = }}\mathit{\boldsymbol{b}}$ 。解总体线性方程组，得每个单元的内部近似函数，结果都为

${u_{h0}}(T) = 1.002 - 0.5001x 。$

而此问题的解析解为 $u = 1 - 0.5x$ 。这说明弱有限元法处理混合剖分网格是可行，且有效的。

2. 渗流分析

2.1 自由面问题

在岩土边坡、土坝、地下洞室及地下水运动等渗流分析中，均存在有渗流自由面问题。渗流自由面的确定是一个重点，也是个难点。求解渗流自由面的方法有虚单元法^[13]、截止负压法^[14]、高斯点法^[15]、变渗透系数法^[16]等，这些方法都是以一般有限元方法为基础，因而精度都受制于一般有限元方法的不足。现在用弱有限元法来求解，看看计算结果的精度如何。

计算模型1：假设有一10 m×10 m的均质土坝，上游水位10 m，下游水位2 m，底部为不透水边界。现把分析区域进行剖分，剖分网格为网格Ⅰ（）：区域被均匀划分成200个直角边为1 m的等腰直角三角形单元，一共有320条单元边界。坝体上游和下游水位以下部分为第一类边界，初始计算时假定溢出段为： $x=10，\text{ }2\le y\le 5$ 。弱有限元空间取 ${W_{1, 0}}(T)$ / ${\left[ {{P_0}(T)} \right]^2}$ ，由于溢出段也是第一类边界，溢出段上水头 $u = y$ ，对于落在溢出段上的单元边界 ${e_i}$ ，需做特殊处理，该单元边界上对应的基函数应取为

图 4 坝体剖分网格Ⅰ

Figure 4. Diagram of dam body meshing Ⅰ

下载: 全尺寸图片幻灯片

$\varphi_{3+i}=\left\{0, \varphi_{b i}\right\}=\left\{\begin{array}{cc} y & (x, y) \in e_i \\ 0 & \text { 其他 } \end{array}\right. \text { 。 }$

通过计算得

$\nabla_{\mathrm{w}} \varphi_{3+i}=\frac{\left|e_i\right|}{|T|} y_i \boldsymbol{n}_i, \quad J_{\partial}\left(\varphi_{3+i}\right)=-h_T^{-1 / 2} y_i \text { 。 }$

所以单元局部矩阵中的M_bb和S_T部分的计算公式需做相应调整。比如，假设三角形单元（）的边界 ${e_2}$ 属于溢出段上， ${e_1}$ 和 ${e_3}$ 属于内部边界，则该单元局部矩阵

$\boldsymbol{A}_T=\boldsymbol{M}_T+\boldsymbol{S}_T,$

其中

$\boldsymbol{M}_T=\left(\begin{array}{cc} O & O \\ O & \boldsymbol{M}_{b b} \end{array}\right),$

$\begin{aligned} & \boldsymbol{M}_{b b}=\frac{1}{|T|}\left(\begin{array}{ccc} \left|e_1\right|^2 & \left|e_1 \| e_2\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_1 & \left|e_1 \| e_3\right| \boldsymbol{n}_3 \boldsymbol{n}_1 \\ \left|e_1 \| e_2\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_1 & \left|e_2\right|^2 y_2^2 & \left|e_2 \| e_3\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_3 \\ \left|e_1 \| e_3\right| \boldsymbol{n}_3 \boldsymbol{n}_1 & \left|e_2 \| e_3\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_3 & \left|e_3\right|^2 \end{array}\right) \\ & \boldsymbol{S}_T=\sum\limits_{i=1}^3\left|e_i\right| \boldsymbol{N}_i^T \boldsymbol{N}_i, \end{aligned}$

其中，

$\boldsymbol{N}_1=h_T^{-1 / 2}\left(1, x_1, y_1,-1,0,0\right) \text {, }$

$\boldsymbol{N}_2= h_T^{ - 1/2}\left(1, {x_2}, {y_2},0,- {y_2},0\right) \text{，}$

$\boldsymbol{N}_3 = h_T^{ - 1/2}\left(1, {x_3}, {y_3}, 0, 0, - 1\right) \text{。}$

处理第一类边界条件时，该边界 ${e_2}$ 上自由度的值为1。根据变单元渗透系数法^[16]的原理，先假设整个坝体渗透系数为1，用弱有限元法初始计算后，根据计算结果把自由面以上单元的渗透系数赋值为1/1000，再进行计算，循环计算3次后得数值结果（表 3），表中数据精确到0.01。

表 3 计算结果对比Ⅰ

Table 3. Comparison of calculated results Ⅰ

横坐标	解析解	FEM	误差1	WG	误差2
0	10.00	10.00	0.00	10.00	0.00
1	9.59	9.73	0.14	9.62	0.03
2	9.17	9.39	0.22	9.26	0.09
3	8.72	8.99	0.27	8.79	0.07
4	8.25	8.53	0.28	8.30	0.05
5	7.75	8.03	0.28	7.78	0.03
6	7.21	7.46	0.25	7.23	0.02
7	6.63	6.83	0.20	6.63	0.00
8	6.00	6.11	0.11	6.04	0.04
9	5.29	5.21	-0.08	5.34	0.05
10	4.47	4.50	0.03	4.50	0.03

下载: 导出CSV

| 显示表格

确定溢出点的方法：根据计算结果得到各节点的水头值后，搜索在坝体下游坝面（ $x = 10, y \geqslant 2$ ）上的单元边界 ${e_i}$ ，计算 ${e_i}$ 的两个节点处的纵坐标值与水头值的差值，如果两个差值异号，则表明出溢点在该单元边界内，再利用计算结果求该单元内部近似函数在边界 $x = 10$ 上的投影，解方程 $y = {Q_{\text{b}}}({u_{\text{h}}}(T)) =$ ${a_{T1}} + {a_{T2}} \times 10 + {a_{T3}}y$ 即可求得溢出点的位置。自由面图见。该问题的自由面有解析解： ${y^2} = 100 - 8x$ ，参考文献[17]中有针对剖分网格Ⅰ，利用虚单元法求解的自由面数据结果，虚单元法是以有限元方法为基础，所以在表 3中对应FEM列。

图 5 自由面

Figure 5. Free surface

下载: 全尺寸图片幻灯片

结果对比显示，在相同剖分网格下，弱有限元法的最大绝对误差仅为0.09，最大相对误差为0.0098，而一般有限元法的最大绝对误差为0.28，最大相对误差为0.0361。这表明弱有限元法的计算结果比一般有限元法的精度明显要高。弱有限元法不足的地方在于：按照1.2节算例中第一类边界条件的处理方法，一般有限元法的自由度个数为121个，而弱有限元法的自由度个数为3×200+320=920个。为了发挥弱有限元法适用于混合剖分网格的优势，减小弱有限元法的自由度个数，根据自由面的预估位置，对区域进行有针对性的剖分（见），在自由面附近区域采用边长为1 m的正方形网格，其它部位采用粗网格单元，一共剖分成57个单元，其中52个四边形，2个五边形和3个六边形，单元边界有133条，弱有限元空间仍然选取 ${W_{1, 0}}(T)/{[{P_0}(T)]^2}$ 。自由度个数减少为3×57+133 =304个。仍然循环计算3次后得计算结果（见表 4）

图 6 坝体剖分网格Ⅱ

Figure 6. Diagram of dam body meshing Ⅱ

下载: 全尺寸图片幻灯片

表 4 计算结果对比Ⅱ

Table 4. Comparison of calculated results Ⅱ

横坐标	解析解	WG	误差
0	10.00	10.00	0.00
1	9.59	9.59	0.00
2	9.17	9.18	0.01
3	8.72	8.77	0.05
4	8.25	8.23	-0.02
5	7.75	7.74	-0.01
6	7.21	7.16	-0.05
7	6.63	6.61	-0.02
8	6.00	6.00	0.00
9	5.29	5.34	0.05
10	4.47	4.50	0.03

下载: 导出CSV

| 显示表格

总体方程组的自由度大大地减少了，而数值解的精度却没有降低，最大绝对误差仅为0.05，最大相对误差为0.0095。这说明，弱有限元法能处理混合网格这一优势是非常显著的。

2.2 闸基渗流场

弱有限元法适用于混合剖分网格，那么在进行复杂的渗流场分析时，可以根据需要采用有针对性的剖分网格，这样可以用较少的自由度得到精度较高的数值计算结果。下面用弱有限元法来分析有防渗帷幕的闸基渗流场。

计算模型2：模型长和高分别为60，20 m，闸基底板长20 m，上游20 m，下游20 m，上游水位深4 m，下游水位默认为0 m，闸基底部默认为不透水边界，在闸基底部靠上游端下方有1 m厚防渗墙（见图 7）。

图 7 闸基模型

Figure 7. Model for brake base

下载: 全尺寸图片幻灯片

闸基渗流分析重点关心防渗墙下端处和闸基底板下游出口处的水力梯度，因此采用有针对性的混合剖分网格（图 8）。

图 8 闸基剖分网格图

Figure 8. Meshes of gate base

下载: 全尺寸图片幻灯片

防渗墙周围和闸基底板下游出口处两个部位采用细网格的四边形单元，其它部位采用大网格四边形单元，大网格与细网格相邻的大网格被认定为八边形单元，因此，单元类型共两种：四边形和八边形，剖分单元一共有192个单元，其中四边形180个，八边形12个，节点共241个，有单元边界432条。弱有限元空间取 ${W_{10}}(T)/{[{P_0}(T)]^2}$ ，闸基渗流不含有需要特殊处理的边界，所以计算格式按1.2节中介绍的公式计算即可。自由度总数为3×192+432=1008个。

防渗墙和地基黏土的渗透系数比假定为1∶1000，对防渗墙深5，10，15 m 3种工况采用弱有限元法进行了渗流计算，计算结果的水头等值线分布如图 9~11（单位：m），防渗墙底端和闸基出口处的水力梯度值见表 5。

图 9 防渗墙深5 m时水头等值线分布图

Figure 9. Distribution of water head contour line at depth of anti-seepage wall of 5 m

下载: 全尺寸图片幻灯片

图 10 防渗墙深10 m时水头等值线分布图

Figure 10. Distribution of water head contour line at depth of anti-seepage wall of 10 m

下载: 全尺寸图片幻灯片

图 11 防渗墙深15 m时水头等值线分布图

Figure 11. Distribution of water head contour line at depth of anti-seepage wall of 15 m

下载: 全尺寸图片幻灯片

表 5 水力梯度值

Table 5. Values of hydraulic gradient

防渗墙深/m	防渗墙底端处水力梯度	闸基出口处水力梯度
5	0.323	0.332
10	0.333	0.262
15	0.358	0.190

下载: 导出CSV

| 显示表格

计算结果表明随着防渗墙的加深下游水头值逐渐降低，防渗墙的防渗效果得到充分体现，也说明了采用有针对性的混合剖分网格，运用弱有限元法来处理类似闸基这种比较复杂的渗流场问题是非常有效的。

3. 结论

计算渗流自由面模型的数值结果表明弱有限元法的计算结果的精度比一般有限元法的要高，并且采用较稀疏的混合剖分网格也能得到精度很好的数值结果。利用弱有限元法适用混合剖分网格的特点，对较复杂的闸基渗流场进行分析，根据数值结果绘制的等水头线表明结果符合闸基渗流的特征。高精度和能适用于混合剖分网格这两个优势，使得弱有限元法在渗流分析中有很好的应用前景。

图 1 三角形单元图

Figure 1. Triangular element

下载: 全尺寸图片幻灯片

图 2 剖分网格图

Figure 2. Diagram of partial grids

下载: 全尺寸图片幻灯片

图 3 单元边界编号

Figure 3. Number of element boundary

下载: 全尺寸图片幻灯片

图 4 坝体剖分网格Ⅰ

Figure 4. Diagram of dam body meshing Ⅰ

下载: 全尺寸图片幻灯片

图 5 自由面

Figure 5. Free surface

下载: 全尺寸图片幻灯片

图 6 坝体剖分网格Ⅱ

Figure 6. Diagram of dam body meshing Ⅱ

下载: 全尺寸图片幻灯片

图 7 闸基模型

Figure 7. Model for brake base

下载: 全尺寸图片幻灯片

图 8 闸基剖分网格图

Figure 8. Meshes of gate base

下载: 全尺寸图片幻灯片

图 9 防渗墙深5 m时水头等值线分布图

Figure 9. Distribution of water head contour line at depth of anti-seepage wall of 5 m

下载: 全尺寸图片幻灯片

图 10 防渗墙深10 m时水头等值线分布图

Figure 10. Distribution of water head contour line at depth of anti-seepage wall of 10 m

下载: 全尺寸图片幻灯片

图 11 防渗墙深15 m时水头等值线分布图

Figure 11. Distribution of water head contour line at depth of anti-seepage wall of 15 m

下载: 全尺寸图片幻灯片

表 1 剖分单元信息

Table 1 Information of elements

单元编号	单元类型	节点编号
1	3	3	2	4
2	3	3	1	2
3	4	3	4	5	6
4	5	3	6	7	8	1

下载: 导出CSV

表 2 节点坐标信息

Table 2 Information of node coordinates

节点	1	2	3	4	5	6	7	8
x	1	2	1	2	2	1	0	0
y	0	0	1	1	2	2	2	0

下载: 导出CSV

表 3 计算结果对比Ⅰ

Table 3 Comparison of calculated results Ⅰ

横坐标	解析解	FEM	误差1	WG	误差2
0	10.00	10.00	0.00	10.00	0.00
1	9.59	9.73	0.14	9.62	0.03
2	9.17	9.39	0.22	9.26	0.09
3	8.72	8.99	0.27	8.79	0.07
4	8.25	8.53	0.28	8.30	0.05
5	7.75	8.03	0.28	7.78	0.03
6	7.21	7.46	0.25	7.23	0.02
7	6.63	6.83	0.20	6.63	0.00
8	6.00	6.11	0.11	6.04	0.04
9	5.29	5.21	-0.08	5.34	0.05
10	4.47	4.50	0.03	4.50	0.03

下载: 导出CSV

表 4 计算结果对比Ⅱ

Table 4 Comparison of calculated results Ⅱ

横坐标	解析解	WG	误差
0	10.00	10.00	0.00
1	9.59	9.59	0.00
2	9.17	9.18	0.01
3	8.72	8.77	0.05
4	8.25	8.23	-0.02
5	7.75	7.74	-0.01
6	7.21	7.16	-0.05
7	6.63	6.61	-0.02
8	6.00	6.00	0.00
9	5.29	5.34	0.05
10	4.47	4.50	0.03

下载: 导出CSV

表 5 水力梯度值

Table 5 Values of hydraulic gradient

防渗墙深/m	防渗墙底端处水力梯度	闸基出口处水力梯度
5	0.323	0.332
10	0.333	0.262
15	0.358	0.190

下载: 导出CSV

参考文献(17)

[1]	BASSI F, REBAY S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[J]. Journal of Computational Physics, 1997, 131(2): 267-279. doi: 10.1006/jcph.1996.5572
[2]	COCKBURN B, SHU C W. The local discontinuous Galerkin method for time-dependent convection-diffusion systems[J]. SIAM Journal on Numerical Analysis, 1998, 35(6): 2440-2463. doi: 10.1137/S0036142997316712
[3]	COCKBURN B, GOPALAKRISHNAN J, LAZAROV R. Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[J]. SIAM Journal on Numerical Analysis, 2009, 47(2): 1319-1365. doi: 10.1137/070706616
[4]	何朝葵, 速宝玉, 盛金昌. 稳定渗流分析的局部间断伽辽金有限元法[J]. 河海大学学报(自然科学版), 2012, 40(2): 206-210. https://www.cnki.com.cn/Article/CJFDTOTAL-HHDX201202018.htm HE Zhaokui, SU Baoyu, SHENG Jinchang. Local discontinuous Galerkin finite element method for steady seepage analysis[J]. Journal of Hohai University (Natural Sciences), 2012, 40(2): 206-210. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HHDX201202018.htm
[5]	WANG J P, YE X. A weak Galerkin finite element method for second-order elliptic problems[J]. Journal of Computational and Applied Mathematics, 2013, 241: 103-115. doi: 10.1016/j.cam.2012.10.003
[6]	王军平, 叶秀, 张然. 弱有限元方法简论[J]. 计算数学, 2016, 38(3): 289-308. https://www.cnki.com.cn/Article/CJFDTOTAL-JSSX202001001.htm WANG Junping, YE Xiu, ZHANG Ran. Basics of weak Galerkin finite element methods[J]. Mathematica Numerica Sinica, 2016, 38(3): 289-308. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JSSX202001001.htm
[7]	朱弘泽, 林莉, 周晨光, 等. 弱Galerkin有限元法的稳定性[J]. 吉林大学学报(理学版), 2018, 56(6): 1427-1430. https://www.cnki.com.cn/Article/CJFDTOTAL-JLDX201806025.htm ZHU Hongze, LIN Li, ZHOU Chenguang, et al. Stability of weak Galerkin finite element method[J]. Journal of Jilin University (Science Edition), 2018, 56(6): 1427-1430. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JLDX201806025.htm
[8]	WANG J P, YE X. A weak Galerkin finite element method for the stokes equations[J]. Advances in Computational Mathematics, 2016, 42(1): 155-174. doi: 10.1007/s10444-015-9415-2
[9]	MU L, WANG J P, YE X, et al. A numerical study on the weak Galerkin method for the Helmholtz equation[J]. Communications in Computational Physics, 2014, 15(5): 1461-1479. doi: 10.4208/cicp.251112.211013a
[10]	MU L, WANG J P, YE X, et al. A weak Galerkin finite element method for the maxwell equations[J]. Journal of Scientific Computing, 2015, 65(1): 363-386. doi: 10.1007/s10915-014-9964-4
[11]	张然. 弱有限元方法在线弹性问题中的应用[J]. 计算数学, 2020, 42(1): 1-17. https://www.cnki.com.cn/Article/CJFDTOTAL-JSSX202001001.htm ZHANG Ran. Weak Galerkin finite element method for linear elasticity problems[J]. Mathematica Numerica Sinica, 2020, 42(1): 1-17. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JSSX202001001.htm
[12]	张然, 翟起龙. 偏微分方程特征值问题的弱有限元方法[J]. 中国科学: 数学, 2019, 49(12): 1979-1994. https://www.cnki.com.cn/Article/CJFDTOTAL-JAXK201912015.htm ZHANG Ran, ZHAI Qilong. Weak Galerkin finite element methods for PDE eigenvalue problems[J]. Scientia Sinica (Mathematica), 2019, 49(12): 1979-1994. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JAXK201912015.htm
[13]	吴梦喜, 张学勤. 有自由面渗流分析的虚单元法[J]. 水利学报, 1994, 25(8): 67-71. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB408.009.htm WU Mengxi, ZHANG Xueqin. Imaginary element method for numerical analysis of seepage with free surface[J]. Journal of Hydraulic Engineering, 1994, 25(8): 67-71. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB408.009.htm
[14]	速宝玉, 沈振中, 赵坚. 用变分不等式理论求解渗流问题的截止负压法[J]. 水利学报, 1996, 27(3): 22-29, 35. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB603.003.htm SU Baoyu, SHEN Zhenzhong, ZHAO Jian. The cut-off negative pressure method for soling filtration problems based on the theory of variational inequalities[J]. Journal of Hydraulic Engineering, 1996, 27(3): 22-29, 35. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB603.003.htm
[15]	王贤能, 黄润秋. 有自由面渗流分析的高斯点法[J]. 水文地质工程地质, 1997, 24(6): 1-4. https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG706.000.htm WANG Xianneng, HUANG Runqiu. Gavss point method for seepage analysis with free surface[J]. Hydrogeology and Engineering Geology, 1997, 24(6): 1-4. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG706.000.htm
[16]	党发宁, 王晓章, 郑忠安, 等. 有自由面渗流分析的变单元渗透系数法[J]. 西北水力发电, 2004, 20(1): 1-3. https://www.cnki.com.cn/Article/CJFDTOTAL-SXFD200401001.htm DANG Faning, WANG Xiaozhang, ZHENG Zhongan, et al. Variable element seepage coefficient method for seepage numerical analysis with free surface[J]. Journal of NOTHUCEST RYDOOELECTNICPOWER, 2004, 20(1): 1-3. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SXFD200401001.htm
[17]	孙伟建, 侯兴民, 李远东, 等. 基于虚单元法求解渗流自由面的曲线拟合法[J]. 水力发电, 2016, 42(11): 50-53, 67. https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD201611013.htm SUN Weijian, HOU Xingmin, LI Yuandong, et al. A curve fitting method based on virtual element method for solving seepage free surface[J]. Water Power, 2016, 42(11): 50-53, 67. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLFD201611013.htm

施引文献(5)

期刊类型引用(3)

1.	宋伟豪，朱晗玥，毛佳，赵兰浩. 三维流体非协调网格隐式插值方法. 水动力学研究与进展A辑. 2024(02): 196-203 . 百度学术
2.	马浩然，林涛，雷万能，荣云，李艺博，刘威. 智能扦检天轨桁架风载安全性分析. 粮油仓储科技通讯. 2024(04): 15-18 . 百度学术
3.	何朝葵，朱永忠，杨凤莲. 一维椭圆型方程界面问题的弱有限元方法. 数学的实践与认识. 2023(08): 237-247 . 百度学术