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.

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 [58]. 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 [1012], heat conduction [13, 14], fracture [1518], 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 Section 2, 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] where is a specific storage coefficient, is the source per unit volume, is the total head, is the derivative of the total head with respect to time, is the permeability matrix, and is the gradient operator, where . Applying the Fourier transform to the governing equation [32, 38] transforms it into the frequency domain as follows: where and are the Fourier transforms of and , respectively, and is the frequency.

As illustrated in Figure 1, the SBFEM presents a local coordinate system . The coordinates of any point along the radial line and inside the domain can be written as [9] where is the shape function matrix.

The differential operator in the Cartesian coordinate system can be transformed to the scaled boundary coordinate system as follows [9] where and the Jacobian matrix at the boundary can be written as

The head function at any point can be expressed as where is the nodal head vector and is the shape function matrix.

Using Equations (4) and (7), the flux can be written as where

Applying the weighted residual method and Green’s theorem and introducing the boundary conditions yield the following equations [9, 23] where

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, on the right-hand side of Equation (10) satisfies . Set in Equation (10). Then, the SBFEM equation for the steady-state seepage field can be written as

Introducing the variable that consists of the nodal water head functions and flux functions yields

The equation can be transformed into a first-order ODE: where the coefficient matrix is a Hamiltonian matrix. The solution for the bounded domain is obtained using the positive eigenvalues of . Hence, can be expressed as

The solution of Equation (18) can be obtained by computing the eigenvalue and eigenvector of the matrix, which yields where the real components of eigenvalues and are negative and positive, respectively.

The general solution of Equation (18) can be obtained as follows:

To obtain a finite solution at the scaling center (), must be equal to zero. The solution in the bounded domain can be written as

The relationship between and is expressed as where the steady-state stiffness matrix of the -element can be expressed as

3.2. Mass Matrix and Transient Solution

To determine the mass matrix of the SBFEM, the dynamic-stiffness matrix at is introduced:

For the bounded domain, the dynamic-stiffness matrix on the boundary formulated in the frequency domain is written as

To obtain the mass matrix of the bounded domain, the low-frequency case in the SBFEM is addressed, in which the dynamic-stiffness of a bounded domain is assumed to be

The steady-state stiffness matrix is computed using Equation (23). Substituting Equation (26) into Equation (25) leads to a constant term independent of , a term in , and higher-order term in , which are neglected. Additionally, the constant term vanishes. The coefficient matrix of can be expressed as

This is a linear equation used to solve the mass matrix . Using the eigenvalues and eigenvectors of matrix , Equation (27) can be written as where

After solving matrix in Equation (29), the mass matrix is obtained by

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: where the nodal water head 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 time , the water head change rate can be expressed as

Equation (35) is substituted into Equation (34), and the equation at time step can be obtained as follows:

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: where is the water head vector.

For transient seepage, AMATRX and RHS are defined as follows:

Figure 2 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 , , , and , which are used to construct the Hamilton matrix using Equation (18). The two eigenvector matrices (,) are constructed using eigenvalue decomposition. Finally, the stiffness matrix and mass matrix of the PSBFEM elements can be obtained.

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. Figure 3 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.

1USER ELEMENT, NODES=5, TYPE=U5, PROPERTIES=2, COORDINATES=2
2 8
3ELEMENT, TYPE=U5, ELSET=E5
4 3,2,3,4,8,7
5UEL PROPERTY, ELEST=E5
6 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: where is the numerical solution and 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 Figure 4(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  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 -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.

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 Figure 7(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  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.

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.

Figure 8 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.

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 () with a Stanford bunny cavity [41, 42] is considered, as shown in Figure 10. The coefficients of permeability in the and directions are considered: . The value of 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 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 12 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.

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 Figure 14(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 , , and 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.

A monitor point is chosen at the bottom dam to compare the PSBFEM and FEM results, as shown in Figure 14(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.

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.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Disclosure

The manuscript of an earlier version (https://arxiv.org/abs/2108.08434) is a preprint [43].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the Yunnan Postdoctoral Orientation Project. We thank Maxine Garcia, PhD, from Liwen Bianji (Edanz) (http://www.liwenbianji.cn) for editing the English text of a draft of this manuscript.