Abstract
This paper presents a PSBFEM approach that integrates the quadtree mesh generation technique based on digital images for solving seepage problems. The quantitative representation of the distribution of geometrical information and material parameters is achieved by utilizing the color intensity of each pixel, which can then be applied to seepage analysis. The presented method addresses the issue of hanging nodes by treating them as nodes of a polygonal element. We validate the proposed technique by solving three benchmark seepage problems. Results show that the image-based domain can be automatically discretized using a quadtree decomposition of the images. Furthermore, the computational efficiency and precision of the PSBFEM surpass that of the standard FEM. Therefore, the proposed technique allows for the convenient automatic discretization of the domain using pixel meshes to solve seepage problems in engineering applications.
1. Introduction
In geotechnical engineering, the seepage problem is a critical aspect of a structure’s safe [1–4]. Seepage directly affects the stability of structures due to the increase in the pore water pressure. The application of the finite element method (FEM) is prevalent in the field of seepage problems [5–8]. A high-quality mesh is essential for accurate numerical solutions in the FEM. It helps to capture the complex geometry, boundaries, and variations in the physical properties of the problem domain. By properly discretizing the domain, the FEM can provide more precise results. However, it is time-consuming work to discrete a high-quality mesh. Sandia National Laboratories estimates that over 80% of the total analysis time is allocated to preprocessing [9]. To alleviate the preprocessing burden, researchers have proposed several alternative approaches, such as the boundary element method (BEM) [10], the isogeometric analysis (IGA) [11], the meshfree method [12–14], and deep neural networks (DNNs).
Due to the complexity of soil behavior under varying hydraulic heads, numerical methods are frequently utilized to obtain accurate results. In the FEM analysis, we need to generate a computational mesh. However, it should be emphasized that mesh generation would require a substantial increase in preprocessing time and the level of human intervention, particularly in instances where the stratigraphy comprises a complex array of intersecting material layers [15]. In seepage analyses, it is commonly assumed that soil layers exhibit homogeneity [16, 17]. However, implementing such a method would require a substantial increase in both the preprocessing time and the level of human intervention, particularly in instances where the stratigraphy comprises a complex array of intersecting material layers.
In order to efficiently and automatically obtain high-quality meshes, researchers proposed a technique for automatic mesh generation based on quadtree data structures [18]. The meshing technique especially is well-suited for an image-driven automated analysis [19]. The color intensity within digital images offers information into geometrical characteristics and material distribution, enabling their effective utilization for numerical analysis [19]. Hanging nodes are generated at the interface of neighboring elements with varying sizes in the quadtree. Hanging nodes can compromise displacement compatibility among neighboring elements, leading to disturbances in the accuracy and convergence of numerical simulations. Consequently, using quadtree meshes for direct numerical simulation in the conventional FEM is not widely adopted [20].
To address the issue of hanging nodes, the researcher proposed the polygonal technique [21–24]. Huynh et al. [24] proposed an extended polygonal finite element method to model the fracture behavior of hyperelastic materials experiencing significant deformation. Biabanaki et al. [25] constructed polygonal elements using various polygonal shape functions to solve frictionless contact-impact problems. Moreover, Li and Cui [26] introduced the N-sided polygonal smoothed finite element method for phase-field fracture simulation, incorporating nonconforming meshes for local refinement. The core idea of these studies is to construct polygonal elements and treat hanging nodes of nonmatching elements as nodes of polygonal elements. Overall, these methods adopt a weak form in each direction of discretization, resulting in relatively low computational accuracy.
The scaled boundary FEM (SBFEM) is a classical approach to handling the presence of hanging nodes on quadtree meshes [17, 19]. The SBFEM integrates the beneficial aspects and properties of both the FEM and BEM, resulting in a unique semianalytical approach. The SBFEM discretizes only in the circumferential direction. It should be emphasized that the governing equations assume a weak form in the circumferential direction while retaining a strong form in the radial direction [27]. Given these benefits, the SBFEM has found widespread application in various engineering scenarios [28–36].
Because the SBFEM only discretizes in the circumferential direction, we could develop a polygonal SBFEM (PSBFEM) to address the hanging nodes, which treats each node in a nonmatching mesh as the node of a generic polygon. Therefore, the PSBFEM could use the quadtree mesh to solve the engineering problem of complex geometry efficiently. Few researchers have recently begun using the SBFEM to solve seepage problems [17, 37–39]. These works demonstrate that the SBFEM exhibits better accuracy and efficiency than the conventional FEM in solving seepage problems. However, currently, the technology of pixel quadtree mesh has not been employed in the seepage analysis utilizing the SBFEM.
This work addresses the challenges of utilizing the SBFEM for seepage analysis in the presence of complex geometries, stratigraphy, and material by the pixel quadtree mesh approach. The proposed method integrates the PSBFEM and pixel quadtree decomposition for seepage analysis. The rest of this study is structured into five sections. Section 2 presents the formulation of the SBFEM for seepage analysis. Section 3 outlines the technique used for generating the pixel quadtree mesh. Numerical examples are provided in Section 4. Finally, Section 5 summarizes the main conclusions of this work.
2. Solution Procedure for the SBFEM
In this study, we explore a seepage problem in a two-dimensional domain, allowing us to derive the governing equations in the following manner [17, 40]: where is a storage coefficient, is the volumetric source term, is the total hydraulic head, denotes the rate of change of the total hydraulic head with respect to time, is the hydraulic conductivity matrix, and is the differential operator.
By subjecting the governing equation to the Fourier transform [39, 41], we can express it in the frequency space as presented below: where the Fourier transforms of and are represented by and , respectively, and denotes the frequency.
Figure 1 illustrates that the SBFEM utilizes a definition of coordinate . It is possible to express the coordinates of any arbitrary point located along the radial direction within the domain using the methodology described in [27]. where and denote the radial and circumferential coordinates, respectively, and denote the shape function matrix.
The transformation from the Cartesian coordinate system to the coordinate system of SBFEM results in the conversion of the differential operator, which can be expressed in the following manner [27]: with and is the Jacobian matrix, which is written as
The hydraulic head field at a point can be written as where denotes the nodal hydraulic head vector and denotes the matrix of shape function.
The flux can be represented in the coordinate system of the SBFEM as follows: with
The radial hydraulic head function emerges as the solution to the SBFEM equation governing the hydraulic head, as presented in [27, 35]: with
The steady-state seepage field in the SBFEM formulation can be derived by substituting in Equation (10), resulting in the following:
Through incorporating the variable ,
Equation (12) can be reformulated as follows: where the coefficient matrix is expressed as
The eigenvalue decomposition of can be written as where and are negative and positive of the real components of eigenvalues, respectively.
By utilizing the following approach, it is possible to derive the overall solution for hydraulic head functions and nodal fluxes:
The solution at the scaling center () can be solved when is equal to zero. Consequently, the expression for the solution within the bounded domain can be represented as
Relation between and can be expressed as
The stiffness matrix for the steady state can be formulated as follows:
We introduce the dynamic-stiffness matrix to obtain the mass matrix in the SBFEM. The following expression can represent the dynamic-stiffness matrix associated with the boundary :
The dynamic-stiffness formulation of the SBFEM equation is given by
In the low-frequency scenario, the dynamic stiffness of the bounded domain is assumed to have the following properties:
By solving Equation (20), the steady-state stiffness matrix is obtained. Upon substituting Equation (23) into Equation (22), a constant term emerges, which is independent of . The higher-order terms can be disregarded, while the constant term can be set to zero. Consequently, the coefficient matrix of can be expressed as follows:
Through the utilization of the eigenvalues and eigenvectors of matrix , Equation (24) can be reformulated as follows: with the transformation
By evaluating the matrix in Equation (26), the mass matrix can be derived, resulting in
By substituting Equation (23) into Equation (21) and applying the inverse Fourier transform, the relationship between nodal flux and hydraulic head within a confined domain can be expressed as a conventional equation in the time domain. This equation utilizes the static stiffness and mass matrices and is illustrated below:
3. Pixel Quadtree Mesh Generation
We introduce the technique of pixel quadtree mesh in this section. A basic illustration of a quadtree discretization is presented in Figure 2. The quadtree discretization of images is aimed at partitioning the image into square elements of approximately equal colors [27].
Moreover, to show the process of quadtree decomposition, we show an pixel digital gray-scale image for employing as an example in Figure 3(a). The image matrix of the gray-scale digital image is shown in Figure 3(b). The color of a gray-scale image varies from 0 to 255, with 0 representing black and 255 representing white. The numbers closer to zero indicate darker shades, whereas those closer to 255 indicate lighter shades or white. Digital images can therefore be described as storing information about the feature of the material, characterized by the gray intensity of the pixels.
(a)
(b)
The quadtree discretization of an image can be achieved using MATLAB’s qtdecomp function. Figure 4 illustrates the quadtree discretization process applied to the digital matrix of the image presented in Figure 3(b). In this example, the color intensity threshold has been set to two. The smallest dimension of an element is one pixel, while the maximum dimension is four pixels. Firstly, the image is decomposed into four squares of the same size ( pixels), as shown in Figure 4(a). The element will be decomposed when there is a more significant color intensity difference than the threshold. Except for the element at the upper-left corner, all other elements are divided into four separate elements of pixels, as depicted in Figure 4(b). The discretization achieved after the second iteration by repeating the identical operations is shown in Figure 4(c). It is worth mentioning that in the final iteration, none of the elements exhibit a color intensity difference larger than the threshold.
(a) Initial decomposition
(b) First iteration
(c) Second iteration
Considering rotational symmetry, Figure 5 illustrates the six types of quadtree mesh, each representing a unique configuration. Therefore, we only need to calculate the stiffness matrix and mass matrix for these six types of parent elements once in the simulation. Other elements’ stiffness and mass matrices can be matched to these six parent elements. Hence, this approach will significantly reduce computational expenses.
4. Numerical Examples
We verify the validity of the proposed method by performing calculations on three numerical examples. Subsequently, we compare the convergence and accuracy of our method with the conventional FEM, where the FEM is implemented using the commercial software ABAQUS. The relative error norm in the results is conducted based on where is the numerical result and is the reference result.
4.1. Seepage Analysis in a Multiple Hole Plane
To validate the proposed approach for the seepage analysis, we solve the steady-state problem for multiple holes’ square plate (). Figure 6 shows the geometric model and the mesh. In order to verify the accuracy and convergence of the proposed approach, four monitor points are selected, namely, (1.5, 3.5), (0.5, 2.5), (2.5, 2.25), and (2.5, 0.75). The permeability coefficient is equal to cm/s. We set the hydraulic head boundary conditions at the bottom and top to 30 m and 80 m, respectively. Furthermore, we conduct a convergence study by progressively refining the meshes according to sequences , , , and . The quadtree mesh, as depicted in Figure 6(b), is utilized in the PSBFEM analysis. The transition of the mesh between holes of varying sizes is efficiently managed without the need for additional human intervention. The quadtree mesh contains 13674 elements. Among them, 2945 are polygonal elements, and 10729 are quadrilateral elements. In addition, the FEM mesh is shown in Figure 6(c).
(a)
(b)
(c)
Figure 7 illustrates the relative error in the hydraulic head as the mesh is refined. Notably, the PSBFEM exhibits superior convergence compared to the FEM as the mesh undergoes refinement. Moreover, Table 1 illustrates the relative error in the hydraulic head at various mesh sizes. A noticeable trend is observed where the relative errors decrease with increasing mesh refinement. The relative errors of the PSBFEM are lower than those of the FEM when considering equivalent mesh sizes. Figure 8 depicts the hydraulic head contours of both the PSBFEM and the FEM. This figure demonstrates a high degree of similarity between the results obtained from the two approaches. Furthermore, the comparative computational costs between the PSBFEM and FEM are illustrated in Figure 9. It is evident that the computational cost of PSBFEM is less than half of FEM. Hence, when we use the six kinds of parent elements to match the model, the calculation cost of PSBFEM will be lower than that of standard FEM.
(a)
(b)
4.2. Transient Seepage Analysis in a Square Plate with a Rabbit Cavity
To showcase the computational capabilities of the proposed method in handling intricate geometries, we analyzed a square plate featuring a cavity modeled after the rabbit [41,42], as illustrated in Figure 10. The square plate is partitioned into four sections based on four distinct grayscale regions. The coefficients of permeability are illustrated in Table 2. The storage coefficient is 0.001 m-1. The selected monitoring points (, , , and ) for assessing and contrasting the outcomes of both PSBFEM and FEM are shown in Figure 10. The hydraulic head boundary conditions were set to 7 m and 2 m for the upper and lower boundaries, respectively. We set the total time and time interval as 200 s and 5 s for the analysis. Moreover, the quadtree mesh is divided into four element sets by color intensity, as illustrated in Figure 11.
(a)
(b)
The hydraulic head history of four monitoring points is depicted in Figure 12, showcasing the congruence of the solutions obtained via PSBFEM with the reference solution for all points. This finding underscores the efficacy of the PSBFEM method in accurately predicting hydraulic head dynamics.
Table 3 presents the relative error of the hydraulic head for the monitor points. The FEM exhibits a relative error of 0.29%, while the PSBFEM demonstrates a lower relative error of 0.12%. Thus, the accuracy of the PSBFEM surpasses that of the FEM. Additionally, the hydraulic head distribution at various times is illustrated in Figure 13, indicating a nearly identical hydraulic head distribution between the two methods. Thus, employing PSBFEM utilizing quadtree meshes exhibits a favorable outcome in addressing intricate geometries related to transient seepage concerns.
4.3. Seepage Simulation considering a Heterogeneous Foundation of Dam
We consider a concrete dam with a heterogeneous foundation to show the capability of pixel quadtree mesh—the ability to handle complex foundations. Figure 14(a) presents the geometry and boundary conditions. The upstream and downstream initially have hydraulic heads of 10 m and 8 m, respectively. Furthermore, Figure 15 illustrates the variation in the upstream hydraulic head. The color intensity of the soil layer is shown in Table 4. Figure 14(b) displays the quadtree, while Figure 14(c) depicts the FEM mesh. It is noted that the soil layers are discretized as five mesh zones by the pixel quadtree mesh generation. The transition of the element between the boundaries of different soils is effectively managed. The mesh size is 2 m. The quadtree mesh contains 6864 polygonal elements and 11872 quadrilateral elements, respectively. The permeability coefficients in the foundation of the dam are shown in Table 4.
(a)
(b)
(c)
To compare the results of the PSBFEM and FEM, Figure 14(a) illustrates the selected monitor point located at the bottom dam. Table 5 shows the hydraulic head at the monitoring point at six different time intervals. For time intervals of 500 and 1000 days, the relative error of the FEM is lower than that of the PSBFEM. However, once the time exceeds 1500 days, the relative error of the PSBFEM becomes smaller than that of the FEM. Thus, it is evident that the relative errors of the PSBFEM are lower than those of the FEM as time increases. Furthermore, Figure 16 presents the hydraulic head history, demonstrating a good agreement between the results obtained by both methods. Additionally, Figure 17 displays the variation of the hydraulic head at different time points utilizing the FEM and PSBFEM methods. The outcomes from these two approaches exhibit remarkable concurrence.
5. Conclusions
This study introduces a novel PSBFEM approach incorporating the pixel quadtree mesh generation technique for solving seepage problems. We demonstrate the ability of the proposed method to solve the seepage problems considering the complex geometries, stratigraphy, and material by simulating three numerical examples.
The results show that digital images can be described as storing information about the material’s geometry and distribution through the pixels’ gray intensity. The utilization of polygonal elements, implemented through the proposed approach, effectively resolves the challenge of mesh incompatibility often caused by hanging nodes between neighboring meshes of different sizes within the quadtree structure. Moreover, when mesh refinement is employed, the PSBFEM demonstrates superior convergence and accuracy compared to the FEM. By employing six types of parent elements to match the model, the PSBFEM achieves a reduced computational cost compared to standard FEM. Hence, the proposed approach proves valuable in accurately and efficiently simulating various complexities and challenges typically encountered in seepage analysis.
Data Availability
The data used to support the findings of this study are available from the corresponding authors upon request.
Conflicts of Interest
The authors affirm that there are no conflicts of interest.
Acknowledgments
The Yunnan Postdoctoral Orientation Project and the Research on Key Technologies of Green Energy Security and Sustainable Development in Yunnan Province (2021KXJGZS01) provided support for this study.