# A Novel Solution for Seepage Problems Implemented in the Abaqus UEL Based on the Polygonal Scaled Boundary Finite Element Method

Authors: Yang Yang; Zongliang Zhang; Yelin Feng; Kun Wang
Journal: Geofluids
 (2022)

---

## Abstract

The scaled boundary finite element method (SBFEM) is a semianalytical computational scheme based on the characteristics of the finite element method (FEM) and boundary element method that combines their respective advantages. In this paper, the SBFEM and polygonal mesh technique are integrated into a new approach to solve steady-state and transient seepage problems. The proposed method is implemented in Abaqus employing a user-defined element (UEL). A detailed implementation of the procedure is presented in which the UEL element is defined, the internal variables RHS and AMATRX are updated, and the stiffness/mass matrix is solved using eigenvalue decomposition. Several benchmark problems are solved to verify the proposed implementation. The results show that the polygonal element of the polygonal SBFEM (PSBFEM) is more accurate than the standard FEM element of the same element size. For transient problems, the results for the PSBFEM and FEM are in excellent agreement. Hence, the proposed method is robust and accurate for solving steady-state and transient seepage problems. The developed UEL source code and the associated input files can be downloaded from GitHub.


---

## Body

## 1. Introduction

Seepage analysis is an essential topic in civil engineering. Changes in soil pore water pressure may significantly affect the stability of structures, such as such as the slope [1, 2], tunnel [3], and earth-rock dam [4]. The finite element method (FEM) is the dominant method for seepage problems [5–8]. However, this method is cumbersome when dealing with singularities. Thus, alternative approaches have been proposed. The scaled boundary FEM (SBFEM) was developed in the 1990s. It is a semianalytical method that attempts to combine the advantages and characteristics of the FEM and the boundary element method into one new approach. In the SBFEM, only the boundaries of the domain are discretized in the circumferential direction. Then, in the radial direction, the partial differential equation is transformed to an ordinary differential equation (ODE), which can be solved analytically [9]. The SBFEM has been applied to many engineering problems, such as wave propagation [10–12], heat conduction [13, 14], fracture [15–18], acoustics [19], seepage [20, 21], elastoplastics [22], and fluids [23, 24].The polygonal SBFEM (PSBFEM) is a novel method that integrates the standard SBFEM and the polygonal mesh technique [22, 25]. Compared with the standard FEM element, a polygon element with more than four edges involves more nodes in the domain and is usually more accurate [22]. Polygons can discretize complex geometry flexibly. Furthermore, polygons have high geometric isotropy and eliminate the mesh dependence caused by the discretization of fixed meshes with standard triangles or quadrangles [26]. These advantages further motivate the choice of polygonal finite elements as an alternative to standard finite elements that use triangles or quadrangles.Recently, an alternative mesh technique has been widely used in geometric discretization. The quadtree algorithm is fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions [27]. As a result of hanging nodes between two adjacent elements of different sizes, the adaptation of quadtree meshes for direct computational analyses within the standard FEM is not widespread [27]. The presence of hanging nodes destroys the displacement compatibility between the adjacent elements. Several methods have been developed in the literature to resolve the displacement incompatibly introduced by hanging nodes for the FEM. Provatidis [28] developed arbitrary-noded large finite elements, which are based on the Coons–Gordon interpolation theory in conjunction with piecewise-linear, cubic B-splines, and Lagrangian interpolation of the potential between the nodal points arranged along the boundary of the problem domain. Duczek et al. [29] proposed a compatible transition element—xNy-element—which provides the capability of coupling different element types. Gupta [30] presented the formulation of the finite element to match one element with two elements side-by-side. The PSBFEM only discretizes in the geometric boundary. Hence, each element in a quadtree mesh is treated as a generic polygon regardless of the hanging nodes. This enables the structure of the quadtree to be exploited for efficient computation. The ability to assume any number of sides also enables the SBFEM to discretize complex curved boundaries.To date, few studies on SBFEM seepage analysis have been reported. Li and Tu [31] used the SBFEM to solve steady-state seepage problems with multimaterial regions. Bazyar and Talebi [32] simulated transient seepage problems in zoned anisotropic soils. Prempramote [33] developed a high-frequency open boundary for the transient seepage analyses of semi-infinite layers with a constant depth. Liu et al. [20] presented an iso-geometric SBFEM using nonuniform rational B-splines for the numerical solution of seepage problems in the unbounded domain. These studies demonstrate that the SBFEM has excellent accuracy, efficiency, and convergence rate. However, the PSBFEM has not been used to solve the seepage problem.Although the SBFEM has matured, these studies only exist as independent code, and the SBFEM is not available in commercial software. Therefore, engineers find it difficult to use this method to solve engineering problems. The commercial software Abaqus has powerful linear or nonlinear, static, or dynamic analysis capabilities [34]. Abaqus/standard analysis also provides a user-defined element (UEL) to define an element with an available option to interface with the code. Several researchers have focused on the implementation of the SBFEM in Abaqus. Yang et al. [35] developed UEL subroutines for steady-state and transient heat conduction analysis using the PSBFEM. Ya et al. [36] implemented an open-source polyhedral SBFEM element for three-dimensional and nonlinear problems through the Abaqus UEL. Yang et al. [37] implemented the SBFEM in Abaqus in linear elastic stress analyses. However, there are no subroutines available that enable the SBFEM to solve the seepage problem in Abaqus at the present time.The primary goal in this study is to implement a novel semianalytical approach by integrating the SBFEM and polygonal mesh technique to solve the seepage problem. This paper is divided into six sections: In Section2, the polygon seepage SBFEM concept is described. In Section 3, the solution procedures are outlined. In Section 4, the implementation of the PSBFEM for seepage problems using the Abaqus UEL subroutine is described. Then, several benchmark examples are presented in Section 5. Finally, the main concluding remarks for this study are presented in Section 6.


## 2. The PSBFEM for the Transient Seepage Problem

The governing equations related to two-dimensional transient seepage flow can be written as [21]
(1)∇⋅k∇h+p−Ssḣ=0,where Ss is a specific storage coefficient, p is the source per unit volume, h is the total head, ḣ is the derivative of the total head with respect to time, k is the permeability matrix, and ∇ is the gradient operator, where ∇=∂/∂xi⟶+∂/∂yj⟶. Applying the Fourier transform to the governing equation [32, 38] transforms it into the frequency domain as follows:
(2)∇⋅k∇h∼+p∼−iωSsh∼=0,where p∼ and h∼ are the Fourier transforms of p and h, respectively, and ω is the frequency.As illustrated in Figure1, the SBFEM presents a local coordinate system ξ,η. The coordinates of any point x∼,y∼ along the radial line and inside the domain can be written as [9]
(3)x∼=ξNηx,y∼=ξNηy,where Nη is the shape function matrix.Figure 1
Coordinate system of the SBFEM: (a)S-element; (b) local coordinate system of the SBFEM.
(a)(b)The differential operator in the Cartesian coordinate system can be transformed to the scaled boundary coordinate system as follows [9]
(4)∇=b1∂∂ξ+1ξb2∂∂η,where
(5)b1=1Jbyb,η00−xb,η,b2=1Jb‐yb00xb,and the Jacobian matrix at the boundary can be written as
(6)Jb=xbybxb,ηyb,η=xbyb,η−ybxb,η.The head function at any point can be expressed as(7)h∼ξ,η=Nuηh∼ξ,where h∼ξ is the nodal head vector and Nuη is the shape function matrix.Using Equations (4) and (7), the flux Q∼ξ,η can be written as
(8)Q∼ξ,η=−kB1ηh∼ξ,ξ+1ξB2ηh∼ξ,where
(9)B1η=b1ηNη,B2η=b2ηNη.Applying the weighted residual method and Green’s theorem and introducing the boundary conditions yield the following equations [9, 23]
(10)E0ξ2h∼ξ,ξζ+E0−E1+E1Tξh∼ξ,ξ−E2−iωM0ξ2h∼ξ=ξFξ,where
(11)E0=∫SB1ηTkB1ηJbdη,(12)E1=∫SB2ηTkB1ηJbdη,(13)E2=∫SB2ηTkB2ηJbdη,(14)M0=∫SNηTSsNηJbdη.


## 3. Solution Procedure for the PSBFEM Equation

### 3.1. Steady-State Solution

Equation (10) is a second-order Euler–Cauchy equation. In this paper, considering that the side of domains is adiabatic or the domain is closed, Fξ on the right-hand side of Equation (10) satisfies Fξ=0. Set ω=0in Equation (10). Then, the SBFEM equation for the steady-state seepage field can be written as
(15)E0ξ2h∼ξ,ξξ+E0+E1T−E1ξh∼ξ,ξ−E2h∼ξ=0.Introducing the variableXξ that consists of the nodal water head functions h∼ξ and flux functions Q∼ξ yields
(16)Xξ=h∼ξQ∼ξ.The equation can be transformed into a first-order ODE:(17)ξXξ,ξ−ZpXξ=0,where the coefficient matrix Zp is a Hamiltonian matrix. The solution for the bounded domain is obtained using the positive eigenvalues of Zp. Hence, Zpcan be expressed as
(18)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of Equation (18) can be obtained by computing the eigenvalue and eigenvector of the matrixZp, which yields
(19)Zpψ11ψ12ψ21ψ22=ψ11ψ12ψ21ψ22λnλp,where the real components of eigenvalues λn and λp are negative and positive, respectively.The general solution of Equation (18) can be obtained as follows:
(20)h∼ξ=ψ11ξ−λ1c1+ψ12ξ−λ2c2,Q∼ξ=ψ21ξ−λ1c1+ψ22ξ−λ2c2.To obtain a finite solution at the scaling center (ξ=0), c2must be equal to zero. The solution in the bounded domain can be written as
(21)h∼ξ=ψ11ξ−λ1c1,Q∼ξ=ψ21ξ−λ1c1.The relationship betweenh∼ξ and Q∼ξ is expressed as
(22)Ksth∼ξ=Q∼ξ,where the steady-state stiffness matrix of the S-element can be expressed as
(23)Kst=ψ21ψ11−1.


### 3.2. Mass Matrix and Transient Solution

To determine the mass matrixM of the SBFEM, the dynamic-stiffness matrix Kξ,ω at ξis introduced:
(24)Kξ,ωh∼ξ=Q~ξ.For the bounded domain, the dynamic-stiffness matrixKω on the boundary ξ=1formulated in the frequency domain is written as
(25)Kω−E1E0−1Kω−E1T−E2+2ωKω,ω−iωM0=0.To obtain the mass matrixMof the bounded domain, the low-frequency case in the SBFEM is addressed, in which the dynamic-stiffness Kω of a bounded domain is assumed to be
(26)Kω=Kst+iωM.The steady-state stiffness matrixKst is computed using Equation (23). Substituting Equation (26) into Equation (25) leads to a constant term independent of iω, a term in iω, and higher-order term in iω, which are neglected. Additionally, the constant term vanishes. The coefficient matrix of iω can be expressed as
(27)−Kst+E1E0−1−IM+ME0−1−Kst+E1T−I+M0=0.This is a linear equation used to solve the mass matrixM. Using the eigenvalues and eigenvectors of matrix ZP, Equation (27) can be written as
(28)I+λbm+mI+λb=ΦbT~TM0ΦbT~,where
(29)m=ΦbT~TMΦbT~.After solving matrixm in Equation (29), the mass matrix M is obtained by
(30)M=ΦbT~−1TmΦbT~−1.


### 3.3. Time Discretization

Equation (26) is substituted into Equation (24); then, the inverse Fourier transform is applied, and the nodal water head relationship for a bounded domain is expressed as a standard time domain equation using steady-state stiffness and mass matrices as follows:
(31)Kstht+Mh·t=Qt,where the nodal water head ht is a continuous derivative of time. Generally, it is difficult to obtain the function solution in the time domain. In this paper, the backward difference method [39, 40] is adopted to solve Equation (31). The time domain is divided into several time units, and the solution for the time node is obtained step by step from the initial conditions, and the nodal water head at any time is obtained by interpolation.At timet,t+Δt, the water head change rate h·t can be expressed as
(32)h·t=ΔhΔt=htt+Δt−httΔt.Equation (35) is substituted into Equation (34), and the equation at time step t+Δtcan be obtained as follows:
(33)Kstt+Δt+Mt+ΔtΔthtt+Δt=Qtt+Δt+Mt+ΔtΔthtt.


## 3.1. Steady-State Solution

Equation (10) is a second-order Euler–Cauchy equation. In this paper, considering that the side of domains is adiabatic or the domain is closed, Fξ on the right-hand side of Equation (10) satisfies Fξ=0. Set ω=0in Equation (10). Then, the SBFEM equation for the steady-state seepage field can be written as
(15)E0ξ2h∼ξ,ξξ+E0+E1T−E1ξh∼ξ,ξ−E2h∼ξ=0.Introducing the variableXξ that consists of the nodal water head functions h∼ξ and flux functions Q∼ξ yields
(16)Xξ=h∼ξQ∼ξ.The equation can be transformed into a first-order ODE:(17)ξXξ,ξ−ZpXξ=0,where the coefficient matrix Zp is a Hamiltonian matrix. The solution for the bounded domain is obtained using the positive eigenvalues of Zp. Hence, Zpcan be expressed as
(18)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of Equation (18) can be obtained by computing the eigenvalue and eigenvector of the matrixZp, which yields
(19)Zpψ11ψ12ψ21ψ22=ψ11ψ12ψ21ψ22λnλp,where the real components of eigenvalues λn and λp are negative and positive, respectively.The general solution of Equation (18) can be obtained as follows:
(20)h∼ξ=ψ11ξ−λ1c1+ψ12ξ−λ2c2,Q∼ξ=ψ21ξ−λ1c1+ψ22ξ−λ2c2.To obtain a finite solution at the scaling center (ξ=0), c2must be equal to zero. The solution in the bounded domain can be written as
(21)h∼ξ=ψ11ξ−λ1c1,Q∼ξ=ψ21ξ−λ1c1.The relationship betweenh∼ξ and Q∼ξ is expressed as
(22)Ksth∼ξ=Q∼ξ,where the steady-state stiffness matrix of the S-element can be expressed as
(23)Kst=ψ21ψ11−1.


## 3.2. Mass Matrix and Transient Solution

To determine the mass matrixM of the SBFEM, the dynamic-stiffness matrix Kξ,ω at ξis introduced:
(24)Kξ,ωh∼ξ=Q~ξ.For the bounded domain, the dynamic-stiffness matrixKω on the boundary ξ=1formulated in the frequency domain is written as
(25)Kω−E1E0−1Kω−E1T−E2+2ωKω,ω−iωM0=0.To obtain the mass matrixMof the bounded domain, the low-frequency case in the SBFEM is addressed, in which the dynamic-stiffness Kω of a bounded domain is assumed to be
(26)Kω=Kst+iωM.The steady-state stiffness matrixKst is computed using Equation (23). Substituting Equation (26) into Equation (25) leads to a constant term independent of iω, a term in iω, and higher-order term in iω, which are neglected. Additionally, the constant term vanishes. The coefficient matrix of iω can be expressed as
(27)−Kst+E1E0−1−IM+ME0−1−Kst+E1T−I+M0=0.This is a linear equation used to solve the mass matrixM. Using the eigenvalues and eigenvectors of matrix ZP, Equation (27) can be written as
(28)I+λbm+mI+λb=ΦbT~TM0ΦbT~,where
(29)m=ΦbT~TMΦbT~.After solving matrixm in Equation (29), the mass matrix M is obtained by
(30)M=ΦbT~−1TmΦbT~−1.


## 3.3. Time Discretization

Equation (26) is substituted into Equation (24); then, the inverse Fourier transform is applied, and the nodal water head relationship for a bounded domain is expressed as a standard time domain equation using steady-state stiffness and mass matrices as follows:
(31)Kstht+Mh·t=Qt,where the nodal water head ht is a continuous derivative of time. Generally, it is difficult to obtain the function solution in the time domain. In this paper, the backward difference method [39, 40] is adopted to solve Equation (31). The time domain is divided into several time units, and the solution for the time node is obtained step by step from the initial conditions, and the nodal water head at any time is obtained by interpolation.At timet,t+Δt, the water head change rate h·t can be expressed as
(32)h·t=ΔhΔt=htt+Δt−httΔt.Equation (35) is substituted into Equation (34), and the equation at time step t+Δtcan be obtained as follows:
(33)Kstt+Δt+Mt+ΔtΔthtt+Δt=Qtt+Δt+Mt+ΔtΔthtt.


## 4. Procedure Implementation

### 4.1. Implementation of the PSBFEM in Abaqus

The most critical task of UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX in the user subroutine interface provided in Abaqus [40]. For the steady-state seepage analysis, AMATRX and RHS are defined as follows:
(34)AMATRX=Kst,(35)RHS=−Ksth,where h is the water head vector.For transient seepage, AMATRX and RHS are defined as follows:(36)AMATRX=Kstt+Δt+Mt+ΔtΔt,RHS=−Kstt+Δtht+Δt−Mt+ΔtΔtht+Δt−ht.Figure2 shows a flow chart of the UEL subroutine for the PSBFEM. Based on the input file’s connectivity information, the UEL computes the scaling centers and transforms the global coordinate into the local coordinate. Equations (11)–(14) compute the coefficient matrices E0, E1, E2, and M0, which are used to construct the Hamilton matrix Zp using Equation (18). The two eigenvector matrices (Φqn,Φun) are constructed using eigenvalue decomposition. Finally, the stiffness matrix K and mass matrix M of the PSBFEM elements can be obtained.Figure 2
Flow chart of the UEL subroutine for the PSBFEM.


### 4.2. Defining the Element of the PSBFEM

The input file of Abaqus usually contains a complete description of the numerical model, such as nodes, elements, degrees of freedom, and materials. This information needs to be defined by the user in the “inp” file. Figure3 shows a simple polygon mesh of the PSBFEM to demonstrate the definition of elements in the UEL. The mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and pentagonal element (U5). As shown in Listing 1, the pentagonal element (U5) is defined as follows: 1–6 are the line numbers. Lines 1–6 are used to define the pentagonal element (U5): Line 1 assigns the element type, number of nodes, number of element properties, and number of degrees of freedom for each node; line 2 sets the active degrees of freedom for the pore pressure; lines 3–4 define the element sets E5; and lines 5–6 set the permeability coefficient of E5.Figure 3
Schematic diagram of the polygon meshes.Listing 1: Input file of the polygon element in Abaqus (cf. Figure 4).
1∗USER ELEMENT, NODES=5, TYPE=U5, PROPERTIES=2, COORDINATES=22 83∗ELEMENT, TYPE=U5, ELSET=E54 3,2,3,4,8,75∗UEL PROPERTY, ELEST=E56 0.003,0.003


## 4.1. Implementation of the PSBFEM in Abaqus

The most critical task of UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX in the user subroutine interface provided in Abaqus [40]. For the steady-state seepage analysis, AMATRX and RHS are defined as follows:
(34)AMATRX=Kst,(35)RHS=−Ksth,where h is the water head vector.For transient seepage, AMATRX and RHS are defined as follows:(36)AMATRX=Kstt+Δt+Mt+ΔtΔt,RHS=−Kstt+Δtht+Δt−Mt+ΔtΔtht+Δt−ht.Figure2 shows a flow chart of the UEL subroutine for the PSBFEM. Based on the input file’s connectivity information, the UEL computes the scaling centers and transforms the global coordinate into the local coordinate. Equations (11)–(14) compute the coefficient matrices E0, E1, E2, and M0, which are used to construct the Hamilton matrix Zp using Equation (18). The two eigenvector matrices (Φqn,Φun) are constructed using eigenvalue decomposition. Finally, the stiffness matrix K and mass matrix M of the PSBFEM elements can be obtained.Figure 2
Flow chart of the UEL subroutine for the PSBFEM.


## 4.2. Defining the Element of the PSBFEM

The input file of Abaqus usually contains a complete description of the numerical model, such as nodes, elements, degrees of freedom, and materials. This information needs to be defined by the user in the “inp” file. Figure3 shows a simple polygon mesh of the PSBFEM to demonstrate the definition of elements in the UEL. The mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and pentagonal element (U5). As shown in Listing 1, the pentagonal element (U5) is defined as follows: 1–6 are the line numbers. Lines 1–6 are used to define the pentagonal element (U5): Line 1 assigns the element type, number of nodes, number of element properties, and number of degrees of freedom for each node; line 2 sets the active degrees of freedom for the pore pressure; lines 3–4 define the element sets E5; and lines 5–6 set the permeability coefficient of E5.Figure 3
Schematic diagram of the polygon meshes.Listing 1: Input file of the polygon element in Abaqus (cf. Figure 4).
1∗USER ELEMENT, NODES=5, TYPE=U5, PROPERTIES=2, COORDINATES=22 83∗ELEMENT, TYPE=U5, ELSET=E54 3,2,3,4,8,75∗UEL PROPERTY, ELEST=E56 0.003,0.003


## 5. Numerical Examples

Four numerical examples are used to demonstrate the convergence and accuracy of the PSBFEM. Additionally, the results of the PSBFEM are compared with those of the standard FEM. The FEM analysis uses the commercial finite element software Abaqus. For validation, the relative errors in the water head are investigated as follows:(37)eL2=h‐hnL2Ω=∫Ωh−hnTh−hndΩ∫ΩhThdΩ,where h is the numerical solution and hn is the analytical or reference solution.


### 5.1. Steady-State Seepage Problem in a Concrete Dam

In the first example, a standard concrete dam foundation steady-state seepage problem is considered. The geometric model and boundary conditions are shown in Figure4(a). The dam is assumed to be impervious. The boundaries BC, AE, EF, and DF are defined as impermeable boundaries. The hydraulic head of AB is 80 m, and the hydraulic head of CD is 20 m. To verify the accuracy of the proposed method, three monitor points 1 (100, 80), 2 (120, 80), and 3 (140, 80) are chosen, as shown in Figure 4(a). The permeability coefficient of the dam foundation is kx=ky=1×10−5 cm/s. The FEM analysis uses the Abaqus CPE4P element in this study, as shown in Figure 4(b). The PSBFEM uses the polygonal element, as shown in Figure 4(c). The water head of the monitor points is presented in Table 1. The relative errors of Abaqus CPE4P and the PSBFEM are 1.38% and 0.90%, respectively. Hence, the PSBFEM is more accurate than the FEM for the same element size. Figure 5 shows that the results are virtually the same for the FEM and PSBFEM. Additionally, a convergence study is performed using mesh h-refinement. The meshes are refined successively following the sequence 20 m, 10 m, 5 m, 2.5 m, and 1.25 m. Figure 6 illustrates the water head of the monitor points at different degrees of freedom. It is noted that the convergence rates of the PSBFEM and Abaqus CPE4P are the same. Moreover, the PSBFEM is more accurate than Abaqus CPE4P for the same degrees of freedom.Figure 4
Steady-state seepage problem in a concrete dam: (a) geometric model and boundary conditions; (b) Abaqus CPE4P element; (c) PSBFEM polygonal element.
(a)(b)(c)Table 1
Comparison of the water head using different methods (element size 20 m).
MethodMonitor pointRelative erroreL2 (%)123Analytical solution (m)605040—Abaqus CPE4P (m)60.865039.141.38PSBFEM (m)60.4050.0639.320.90Figure 5
Water head distribution of the dam foundation: (a) Abaqus CPE4P element; (b) PSBFEM element.
(a)(b)Figure 6
Comparison of the convergence rates in the water head.


### 5.2. Steady-State Seepage Analysis in Permeable Materials

To demonstrate the flexibility of the PSBFEM using the quadtree mesh, a steady-state seepage problem for permeable materials is solved. The geometry is shown in Figure7(a), in which the interior of the permeable material contains an impermeable material. The length and width of the permeable material are both 1 m. The coefficient of permeability is kx=ky=5×10−4 cm/s. The quadtree mesh is shown in Figure 7(b). The mesh size for the quadtree mesh is the same for each material at the junction, and the mesh transition area can be effectively processed without further manual intervention.Figure 7
Permeable material’s geometric model and quadtree mesh: (a) geometric model; (b) Abaqus CPE4P mesh; (c) PSBFEM quadtree mesh.
(a)(b)(c)Moreover, the impermeable material does not divide the mesh, and only the impermeable boundary is set at the junction with the permeable material. To verify the accuracy of the quadtree mesh calculation, the results of the Abaqus CP4EP element with similar degrees of freedom are compared. The degrees of freedom for the quadtree and CPE4P elements are 11,447 and 11,749, respectively.Figure8 shows a comparison of the PSBFEM quadtree element and Abaqus CPE4P element in the water head. The relative errors of the left edge and right edge are 0.28% and 0.32%, respectively. Furthermore, the distributions of the water head obtained from the PSBFEM and FEM are illustrated in Figure 9. It can be observed that the contour plots present good agreement. Therefore, these results demonstrate the accuracy and reliability of the PSBFEM for the quadtree mesh.Figure 8
Comparison of the PSBFEM and FEM in the water head: (a) left edge; (b) right edge.
(a)(b)Figure 9
Water head distribution: (a) Abaqus CPE4P; (b) PSBFEM quadtree element.
(a)(b)


### 5.3. Transient Seepage Analysis for Complex Geometry

To demonstrate the PSBFEM’s ability to solve complex geometry in the transient seepage problem, a square plate (L=4.0m) with a Stanford bunny cavity [41, 42] is considered, as shown in Figure 10. The coefficients of permeability in the x and y directions are considered: kx=ky=5×10−6m/s. The value of Ss is 0.001 m-1. As shown in Figure 10, four monitor points A, B, C, and D are chosen to compare the results between the FEM and PSBFEM. The water heads at the top and bottom boundaries are specified as 10 m and 3 m, respectively. The total time 2000 s and time step Δt=1s are used for the PSBFEM and FEM. The PSBFEM and FEM are modeled using the quadtree element and CPE4P element, respectively. As shown in Figure 11, both approaches use the same element size.Figure 10
Geometry and boundary conditions for a square plate with a Stanford bunny cavity.Figure 11
Meshes of a square plate with a Stanford bunny cavity: (a) Abaqus mesh; (b) quadtree mesh.
(a)(b)Figure12 illustrates the history of the water head of the four monitor points. The solutions obtained by the PSBFEM are in excellent agreement with those obtained by the FEM for all points. When the time is greater than 1500 s, the water head of all points becomes stable. Additionally, the water head of the four monitor points at different times is presented in Table 2, which shows that the relative error of the four nodes is less than 1.6%. It is noted that the relative error reduces as the increment of time decreases. Furthermore, Figure 13 shows the distribution of the water head at different times. The water head distributions are virtually the same for the FEM and PSBFEM. Therefore, the PSBFEM with quadtree meshes demonstrates a good effect for solving complex geometry in the transient seepage problem.Figure 12
Comparison of the water head history of the four monitor points for the FEM and PSBFEM.Table 2
Water head of the monitor points at different times.
Time (s)ABCDRelative error (%)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)5007.5177.5714.6854.7833.2953.3484.8454.9401.5710008.4038.4286.0286.0704.2054.2286.2036.2440.5515008.5788.5886.2976.3174.3884.3996.4726.4900.2420008.6128.6186.3516.3644.4254.4316.5256.5370.15Figure 13
Water head distribution at different times using the FEM and PSBFEM.


### 5.4. Transient Seepage Analysis in a Concrete Dam with an Orthotropic Foundation

In the final example, the PSBFEM is applied to simulate a concrete dam with an orthotropic foundation. The geometry and boundary conditions are shown in Figure14(a). The initial water levels upstream and downstream are 10 m and 5 m, respectively. Moreover, the change of the water level upstream of the reservoir is illustrated in Figure 15. The material properties of kx, ky, and Ss are 0.001 m/day, 0.0005 m/day, and 0.001 m-1, respectively. Quadrilateral and polygonal meshes are used in Abaqus and the PSBFEM, respectively, as shown in Figures 14(b) and 14(c). The mesh size is 5 m.Figure 14
Transient seepage under a concrete dam constructed on anisotropic soil: (a) geometry and boundary conditions; (b) FEM mesh; (c) PSBFEM mesh.
(a)(b)(c)Figure 15
Variation of the water head with time upstream of a concrete dam constructed on anisotropic soil.A monitor point is chosen at the bottom dam to compare the PSBFEM and FEM results, as shown in Figure14(a). Table 3 shows the water head of the monitor points at six times. When the time is 500 days, the relative error of the water head is 3.958%. It can be observed that the relative error is 0.092% when the time is 3000 days. The history of the water head is shown in Figure 16, where the results obtained by the two methods correspond well. It is noted that the water head becomes stable when the time is more than 2000 days. Moreover, Figure 17 illustrates the distribution of the water head at different times using the PSBFEM and FEM. The results for the two methods are in excellent agreement.Table 3
Water head of the monitor points at different times.
Time (day)Monitor pointRelative error (%)FEM (m)PSBFEM (m)5005.8336.0643.958100012.81912.9961.379150015.77015.8790.695200016.87916.9370.341250017.28017.3090.170300017.42217.4380.092Figure 16
Comparison of the FEM and PSBFEM for the water head history of the monitor point.Figure 17
Comparison of the FEM and PSBFEM for the water head distribution at different times.


## 5.1. Steady-State Seepage Problem in a Concrete Dam

In the first example, a standard concrete dam foundation steady-state seepage problem is considered. The geometric model and boundary conditions are shown in Figure4(a). The dam is assumed to be impervious. The boundaries BC, AE, EF, and DF are defined as impermeable boundaries. The hydraulic head of AB is 80 m, and the hydraulic head of CD is 20 m. To verify the accuracy of the proposed method, three monitor points 1 (100, 80), 2 (120, 80), and 3 (140, 80) are chosen, as shown in Figure 4(a). The permeability coefficient of the dam foundation is kx=ky=1×10−5 cm/s. The FEM analysis uses the Abaqus CPE4P element in this study, as shown in Figure 4(b). The PSBFEM uses the polygonal element, as shown in Figure 4(c). The water head of the monitor points is presented in Table 1. The relative errors of Abaqus CPE4P and the PSBFEM are 1.38% and 0.90%, respectively. Hence, the PSBFEM is more accurate than the FEM for the same element size. Figure 5 shows that the results are virtually the same for the FEM and PSBFEM. Additionally, a convergence study is performed using mesh h-refinement. The meshes are refined successively following the sequence 20 m, 10 m, 5 m, 2.5 m, and 1.25 m. Figure 6 illustrates the water head of the monitor points at different degrees of freedom. It is noted that the convergence rates of the PSBFEM and Abaqus CPE4P are the same. Moreover, the PSBFEM is more accurate than Abaqus CPE4P for the same degrees of freedom.Figure 4
Steady-state seepage problem in a concrete dam: (a) geometric model and boundary conditions; (b) Abaqus CPE4P element; (c) PSBFEM polygonal element.
(a)(b)(c)Table 1
Comparison of the water head using different methods (element size 20 m).
MethodMonitor pointRelative erroreL2 (%)123Analytical solution (m)605040—Abaqus CPE4P (m)60.865039.141.38PSBFEM (m)60.4050.0639.320.90Figure 5
Water head distribution of the dam foundation: (a) Abaqus CPE4P element; (b) PSBFEM element.
(a)(b)Figure 6
Comparison of the convergence rates in the water head.


## 5.2. Steady-State Seepage Analysis in Permeable Materials

To demonstrate the flexibility of the PSBFEM using the quadtree mesh, a steady-state seepage problem for permeable materials is solved. The geometry is shown in Figure7(a), in which the interior of the permeable material contains an impermeable material. The length and width of the permeable material are both 1 m. The coefficient of permeability is kx=ky=5×10−4 cm/s. The quadtree mesh is shown in Figure 7(b). The mesh size for the quadtree mesh is the same for each material at the junction, and the mesh transition area can be effectively processed without further manual intervention.Figure 7
Permeable material’s geometric model and quadtree mesh: (a) geometric model; (b) Abaqus CPE4P mesh; (c) PSBFEM quadtree mesh.
(a)(b)(c)Moreover, the impermeable material does not divide the mesh, and only the impermeable boundary is set at the junction with the permeable material. To verify the accuracy of the quadtree mesh calculation, the results of the Abaqus CP4EP element with similar degrees of freedom are compared. The degrees of freedom for the quadtree and CPE4P elements are 11,447 and 11,749, respectively.Figure8 shows a comparison of the PSBFEM quadtree element and Abaqus CPE4P element in the water head. The relative errors of the left edge and right edge are 0.28% and 0.32%, respectively. Furthermore, the distributions of the water head obtained from the PSBFEM and FEM are illustrated in Figure 9. It can be observed that the contour plots present good agreement. Therefore, these results demonstrate the accuracy and reliability of the PSBFEM for the quadtree mesh.Figure 8
Comparison of the PSBFEM and FEM in the water head: (a) left edge; (b) right edge.
(a)(b)Figure 9
Water head distribution: (a) Abaqus CPE4P; (b) PSBFEM quadtree element.
(a)(b)


## 5.3. Transient Seepage Analysis for Complex Geometry

To demonstrate the PSBFEM’s ability to solve complex geometry in the transient seepage problem, a square plate (L=4.0m) with a Stanford bunny cavity [41, 42] is considered, as shown in Figure 10. The coefficients of permeability in the x and y directions are considered: kx=ky=5×10−6m/s. The value of Ss is 0.001 m-1. As shown in Figure 10, four monitor points A, B, C, and D are chosen to compare the results between the FEM and PSBFEM. The water heads at the top and bottom boundaries are specified as 10 m and 3 m, respectively. The total time 2000 s and time step Δt=1s are used for the PSBFEM and FEM. The PSBFEM and FEM are modeled using the quadtree element and CPE4P element, respectively. As shown in Figure 11, both approaches use the same element size.Figure 10
Geometry and boundary conditions for a square plate with a Stanford bunny cavity.Figure 11
Meshes of a square plate with a Stanford bunny cavity: (a) Abaqus mesh; (b) quadtree mesh.
(a)(b)Figure12 illustrates the history of the water head of the four monitor points. The solutions obtained by the PSBFEM are in excellent agreement with those obtained by the FEM for all points. When the time is greater than 1500 s, the water head of all points becomes stable. Additionally, the water head of the four monitor points at different times is presented in Table 2, which shows that the relative error of the four nodes is less than 1.6%. It is noted that the relative error reduces as the increment of time decreases. Furthermore, Figure 13 shows the distribution of the water head at different times. The water head distributions are virtually the same for the FEM and PSBFEM. Therefore, the PSBFEM with quadtree meshes demonstrates a good effect for solving complex geometry in the transient seepage problem.Figure 12
Comparison of the water head history of the four monitor points for the FEM and PSBFEM.Table 2
Water head of the monitor points at different times.
Time (s)ABCDRelative error (%)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)FEM (m)PSBFEM (m)5007.5177.5714.6854.7833.2953.3484.8454.9401.5710008.4038.4286.0286.0704.2054.2286.2036.2440.5515008.5788.5886.2976.3174.3884.3996.4726.4900.2420008.6128.6186.3516.3644.4254.4316.5256.5370.15Figure 13
Water head distribution at different times using the FEM and PSBFEM.


## 5.4. Transient Seepage Analysis in a Concrete Dam with an Orthotropic Foundation

In the final example, the PSBFEM is applied to simulate a concrete dam with an orthotropic foundation. The geometry and boundary conditions are shown in Figure14(a). The initial water levels upstream and downstream are 10 m and 5 m, respectively. Moreover, the change of the water level upstream of the reservoir is illustrated in Figure 15. The material properties of kx, ky, and Ss are 0.001 m/day, 0.0005 m/day, and 0.001 m-1, respectively. Quadrilateral and polygonal meshes are used in Abaqus and the PSBFEM, respectively, as shown in Figures 14(b) and 14(c). The mesh size is 5 m.Figure 14
Transient seepage under a concrete dam constructed on anisotropic soil: (a) geometry and boundary conditions; (b) FEM mesh; (c) PSBFEM mesh.
(a)(b)(c)Figure 15
Variation of the water head with time upstream of a concrete dam constructed on anisotropic soil.A monitor point is chosen at the bottom dam to compare the PSBFEM and FEM results, as shown in Figure14(a). Table 3 shows the water head of the monitor points at six times. When the time is 500 days, the relative error of the water head is 3.958%. It can be observed that the relative error is 0.092% when the time is 3000 days. The history of the water head is shown in Figure 16, where the results obtained by the two methods correspond well. It is noted that the water head becomes stable when the time is more than 2000 days. Moreover, Figure 17 illustrates the distribution of the water head at different times using the PSBFEM and FEM. The results for the two methods are in excellent agreement.Table 3
Water head of the monitor points at different times.
Time (day)Monitor pointRelative error (%)FEM (m)PSBFEM (m)5005.8336.0643.958100012.81912.9961.379150015.77015.8790.695200016.87916.9370.341250017.28017.3090.170300017.42217.4380.092Figure 16
Comparison of the FEM and PSBFEM for the water head history of the monitor point.Figure 17
Comparison of the FEM and PSBFEM for the water head distribution at different times.


## 6. Conclusions

In this study, the PSBFEM was proposed, which integrates the SBFEM and polygonal mesh technique to solve seepage problems. Several benchmark problems were solved to validate the implementation of the PSBFEM against the FEM.For steady-state problems, the accuracy rate of the polygonal element of the PSBFEM was higher than that of the standard FEM element for the same element size. The PSBFEM converged to an analytical solution with an optimal convergence rate. For transient problems, the results for the PSBFEM the FEM were in excellent agreement. Furthermore, the PSBFEM with quadtree meshes demonstrated a good effect for solving complex geometry in the seepage problem. Hence, the proposed method is robust and accurate for solving steady-state and transient seepage problems.



---

*Source: 5797014-2022-03-10.xml*