Abstract
A novel boundary-type meshless method for modeling geofluid flow in heterogeneous geological media was developed. The numerical solutions of geofluid flow are approximated by a set of particular solutions of the subsurface flow equation which are expressed in terms of sources located outside the domain of the problem. This pioneering study is based on the collocation Trefftz method and provides a promising solution which integrates the T-Trefftz method and F-Trefftz method. To deal with the subsurface flow problems of heterogeneous geological media, the domain decomposition method was adopted so that flux conservation and the continuity of pressure potential at the interface between two consecutive layers can be considered in the numerical model. The validity of the model is established for a number of test problems. Application examples of subsurface flow problems with free surface in homogenous and layered heterogeneous geological media were also carried out. Numerical results demonstrate that the proposed method is highly accurate and computationally efficient. The results also reveal that it has great numerical stability for solving subsurface flow with nonlinear free surface in layered heterogeneous geological media even with large contrasts in the hydraulic conductivity.
1. Introduction
Numerical approaches to the simulation of various subsurface flow phenomena using the mesh-based methods such as the finite difference method or the finite element method are well documented in the past [1–5]. Differing from conventional mesh-based methods, the meshless method has the advantages that it does not need the mesh generation. The meshless method has attracted considerable attention in recent years because of its flexibility in solving practical problems involving complex geometry in subsurface flow problems [6–9]. Chen et al. [10] conducted a comprehensive review of mesh-free methods and addressed that mesh-free methods have emerged into a new class of computational methods with considerable success. Subsurface flow problems are usually governed by second-order partial differential equations. Problems involving regions of irregular geometry are generally intractable analytically. For such problems, the use of numerical methods, especially the boundary-type meshless method, to obtain approximate solutions is advantageous.
Several meshless methods have been reported, such as the Trefftz method [11–16], the method of fundamental solutions [7, 17–19], the element-free Galerkin method [20], the reproducing kernel particle method [21, 22], the meshless local boundary integral equation method [23, 24], and the meshless local Petrov-Galerkin approach [25]. Proposed by Trefftz in 1926 [16], the Trefftz method is probably one of the most popular boundary-type meshless methods for solving boundary-value problems where approximate solutions are expressed as a linear combination of functions automatically satisfying governing equations. According to Kita and Kamiya [12], Trefftz methods are classified as either direct or indirect formulations. Unknown coefficients are determined by matching boundary conditions. Li et al. [14] provided a comprehensive comparison of the Trefftz method, collocation, and other boundary methods, concluding that the collocation Trefftz method (CTM) is the simplest algorithm and provides the most accurate solutions with optimal numerical stability.
To solve subsurface flow problems with the layered soil in heterogeneous porous media, the domain decomposition method (DDM) [26] is adopted because the DDM is natural from the physics of the problem to deal with different values of hydraulic conductivity in subdomains. The DDM can be divided into overlapping domain decomposition and nonoverlapping domain decomposition methods. In overlapping domain decomposition methods, the subdomains overlap by more than the interface. In nonoverlapping methods, the subdomains intersect only on their interface. One may need to use the DDM which decomposes the problem domain into several simply connected subdomains and to use the numerical method in each one. In this study, we adopted the nonoverlapping method to deal with the seepage problems of layered soil profiles. The problems on the subdomains are independent, which makes the DDM suitable for describing the layered soil in heterogeneous porous media.
The subsurface flow problem with a free surface is a nonlinear problem in which nonlinearities arise from the nonlinear boundary characteristics [27]. Such nonlinearities are handled in the numerical modeling using iterative schemes [28]. Techniques for solving problems with nonlinear boundary conditions have been investigated. Typically, the methods, such as the Picard method or Newton’s method, are iterative in that they approach the solution through a series of steps. Since the computation of the subsurface flow problem with a free surface has to be solved iteratively, the location of the boundary collocation points and the source points must be updated simultaneously with the moving boundary. Solving subsurface flow with a nonlinear free surface in layered heterogeneous soil is generally much more challenging. In addition, the convergence problems often arise from nonlinear phenomena. A previous study [28] indicates that the Picard scheme is a simple and effective method for the solution of nonlinear and saturated groundwater flow problems. Therefore, we adopted the Picard scheme to find the solution of the nonlinear free surface.
In this paper, we proposed a novel boundary-type meshless method. This pioneering study is based on the collocation Trefftz method and provides a promising solution which integrates the T-Trefftz method and F-Trefftz method for constructing its basis function using one of the particular solutions which satisfies the governing equation and allows many source points outside the domain of interest. To the best of the authors’ knowledge, the pioneering work has not been reported in previous studies and requires further research. Two important phenomena in subsurface flow modeling were explored in this study using the proposed method. We first adopted the domain decomposition method integrated with the proposed boundary-type meshless method to deal with the subsurface flow problems of heterogeneous geological media. The flux conservation and the continuity of pressure potential at the interface between two consecutive layers can be considered in the numerical model. Then, we attempted to utilize the proposed method to solve the geofluid flow with free surface in heterogeneous geological media.
The validity of the model is established for a number of test problems, including the investigation of the basis function using two possible particular solutions and the comparison of the numerical solutions using different particular solutions and the method of fundamental solutions. Application examples of subsurface flow problems with free surface were also carried out.
2. Solutions to the Subsurface Flow Equation in Cylindrical Coordinates
Consider a three-dimensional domain enclosed by a boundary . The steady-state subsurface flow equation can be expressed as withwhere denotes the outward normal direction, denotes the boundary where the Dirichlet boundary condition is given, and denotes the boundary where the Neumann boundary condition is given. Equation (1) is also known as the Laplace equation. In this study, we adopted the cylindrical coordinate system. In the cylindrical coordinate system, the Laplace governing equation can be written aswhere , , and are the radius, polar angle, and altitude in the three-dimensional cylindrical coordinate system. is the unknown function to be solved. Considering a two-dimensional domain in the polar coordinate, the Laplace governing equation can be written aswhere and are the radius and polar angle in the two-dimensional polar coordinate system. For the Laplace equation, the particular solutions can be obtained using the method of separation of variables. The particular solutions of (4) include the following basis functions:The definition of the particular solution in this study is in a wide sense which is to satisfy the homogenous or the nonhomogenous differential equations with or without part of boundary conditions. If we adopt the solution of a boundary value problem and enforce it to exactly satisfy the partial differential equation with the boundary conditions at a set of points, this leads to the CTM.
The CTM belongs to the boundary-type meshless method which can be categorized into the T-Trefftz method and F-Trefftz method. The T-Trefftz method introduces the T-complete functions where the solutions can be expressed as a linear combination of the T-complete functions automatically satisfying governing equations. On the other hand, the F-Trefftz method constructs its basis function space by allowing many source points outside the domain of interest. The solutions are approximated by a set of fundamental solutions which are expressed in terms of sources located outside the domain of the problem. The T-Trefftz method and the F-Trefftz method both required the evaluation of a coefficient for each term in the series. The evaluation of coefficients may be obtained by solving the unknown coefficients in the linear combination of the solutions which are accomplished by collocation imposing the boundary condition at a finite number of points.
The CTM begins with the consideration of T-complete functions. For indirect Trefftz formulation, the approximated solution at the boundary collocation point can be written as a linear combination of the basis functions. For a simply connected domain or infinite domain with a cavity, as illustrated in Figures 1(a) and 1(b), one usually locates the source point inside the domain or the cavity and the number of source points is only one for in the CTM [29].
(a) A simply connected domain
(b) An infinite domain with a cavity
(c) A doubly connected domain
(d) A multiply connected domain
For the doubly and multiply connected domains with genus greater than one, as illustrated in Figures 1(c) and 1(d), one may locate many source points in the domain. Usually, at least one source point inside the cavity is required. If we considered a simply connected domain, the T-complete basis functions can be expressed aswhere and . and is the order of the T-complete function for approximating the solution. For an infinite domain with a cavity as illustrated in Figure 1(b), one usually locates the source point inside the cavity, and the T-complete functions (negative T-complete set) includeThe accuracy of the solution for the CTM depends on the order of the basis functions. Usually, one may need to increase the value to obtain better accuracy. However, the ill-posed behavior also grows up with the value.
On the other hand, there is another type of the Trefftz method, namely, the F-Trefftz method, or the so-called method of the fundamental solutions (MFS) [14]. Instead of using only one source point and increasing the order of basis function, the MFS allows many source points outside the domain of interest. The solutions are approximated by a set of fundamental solutions which are expressed in terms of sources located outside the domain of the problem. Figures 2(a), 2(b), 2(c), and 2(d) illustrate the collocation of the boundary and the source points for a simply connected domain, an infinite domain with a cavity, doubly connected domains, and a multiply connected domain, respectively.
(a) A simply connected domain
(b) An infinite domain with a cavity
(c) A doubly connected domain
(d) A multiply connected domain
The unknown coefficients in the linear combination of the fundamental solutions which are accomplished by collocation imposing the boundary condition at a finite number of points can then be solved. Due to its singular free and meshless merits, the indirect type F-Trefftz method is commonly used. An approximation solution of the two-dimensional Laplace equation using the MFS can also be obtained aswhere and and is the number of source points which are placed outside the domain. The fundamental solution of Laplace equation can be expressed as is defined as the distance between the boundary point and source point, and . Then we selected a finite number of collocation points over the boundary and imposed the boundary condition at boundary collocation points to determine the coefficients of and for the CTM and the MFS, respectively.
For the conventional Trefftz method, the number of source points is only one. Theoretically, one may increase the accuracy by using a larger order of the basis functions [30]. Instead of using only one source point and increasing the order of basis functions, the MFS allows many source points but uses only one basis function, that is, the fundamental solution of the differential operator. One may be interested to investigate a method similar to the MFS which allows many source points but uses other basis functions.
In the following, we proposed a novel boundary-type meshless method. This pioneering study is based on the collocation Trefftz method and provides a promising solution which integrates the T-Trefftz method and F-Trefftz method for constructing its basis function using one of the particular solutions which satisfies the governing equation and allows many source points outside the domain of interest. Differing from the CTM and the MFS, the numerical solutions of the proposed method are approximated by a set of basis functions which are expressed in terms of source points located outside the domain. An approximation solution of the two-dimensional steady-state subsurface flow equation using the proposed method can be obtained aswhere is the spatial coordinate which is collocated on the boundary, is the source point, and is the number of source points which are placed outside the domain. The unknown coefficients can be expressed as . is the particular solution of Laplace equation. In this study, two different particular solutions of Laplace equation were adopted as the basis functions. Two possible particular solutions of Laplace equation can be expressed asThe determination of the unknown coefficients for the proposed method is exactly the same with those in the MFS as described in previous section. We first selected a finite number of collocation points over the boundary such thatwhere are the constant coefficients to be solved, and is the boundary condition imposed at boundary collocation points. Considering the boundary conditions, we havewhere represents the Dirichlet boundary condition; represents the Neumann boundary condition. Applying Dirichlet and Neumann boundary conditions, we obtainedwhere and . is the Neumann boundary condition imposed at boundary collocation points. The source points are on the artificial fictitious boundary, which are placed outside the domain to avoid the singularity of the solution at origin. The artificial fictitious boundary is often chosen as a circle with a radius. However, the position of source and collocation points may affect the accuracy. In order to determine the unknowns , collocating the numerical expansion of (12) at boundary conditions of (14) at boundary collocation points yields the following equations:where is a matrix which takes values of the solutions at the corresponding collocation points and source points, is a vector of unknown coefficients, and is a vector of function values at collocation points.
3. Validation of the Proposed Method
3.1. Investigation of the Basis Function
In this example, we adopted two possible particular solutions of Laplace equation as the basis functions. They are and where and , respectively. In this example, we verified the accuracy of the proposed method and also compared the numerical solution with the MFS. To compare the results with the analytical solution, we considered the subsurface flow problem with an exact solution.
For a two-dimensional simply connected domain enclosed by a boundary, the subsurface flow equation can be expressed as Laplace governing equation which can be expressed asThe two-dimensional object boundary under consideration is defined aswhere , .
The analytical solution can be found asThe Dirichlet boundary condition is imposed on the amoeba-like boundary by using the analytical solution as shown in (18) for the problem.
Figure 3 shows the collocation point for the boundary and the source points. To obtain a promising result of the location of the source points for the proposed method in this study, a sensitivity study was first carried out. An algorithm similar to the study conducted by Chen et al. [31] was adopted with scaling of the artificial boundary with the domain size. Assuming the boundary collocation points can be described as a known parametric representation as follows:The source points can also be described as a known parametric representation from the above equation:where is the dilation parameter and is greater than one. and are the angles of the discretization of the boundary for boundary and source points, respectively. and are the radiuses which represent the scale of the domain size for boundary and source points, respectively. The sensitivity example under investigation is in a simply connected domain. In this example, we investigated the accuracy by choosing locations of the source points through different values using the MFS. Figure 4 shows that could be the satisfactory location of the source points.
Using , we conducted an example to clarify the approximate number of boundary collocation and the source points. For simplicity, we took the same number of the boundary points. To examine the accuracy, we collocated 1074 uniformed distributed inner points inside the domain, as shown in Figure 5. The maximum absolute error can then be found by evaluating the absolute error for each inner point.
Figure 5 depicted the computed results of the maximum absolute error versus the number of source points. It is well known that the linear algebraic equation systems may be ill-conditioned while the global basis functions were adopted. To clarify this issue, we investigated the condition number versus the number of source points. Figure 6 shows that the relationship of the condition number versus the number of source points for the proposed method and the MFS. For simplicity, we adopted the commercial program MATLAB backslash operator to solve the linear algebraic equation systems. It is found that the proposed method remains relatively high accuracy compared to the MFS in this example. The best accuracy can reach the order of while the number of source points is greater than 180. On the other hand, the best accuracy of the MFS can reach only about in the same example.
3.2. Comparison of the Numerical Results
Similar to the previous example, we verified the accuracy of the proposed method with the consideration of a complex star-like boundary. For a two-dimensional simply connected domain enclosed by a boundary, the subsurface flow equation can be expressed as Laplace governing equation which can be expressed asThe two-dimensional object boundary under consideration is defined aswhere , .
The analytical solution can be found asThe Dirichlet boundary condition is imposed on the boundary by using the analytical solution as shown in (23) for the problem. Figure 7 shows the collocation point for the boundary and the source points. A sensitivity study using the MFS was first carried out and could be the satisfactory location of the source points, as shown in Figure 8. Also, to examine the accuracy, we collocated 3250 uniformed distributed inner points inside the domain, as shown in Figure 9. The maximum absolute error can then be found by evaluating the absolute error for each inner point.
Figure 9 depicted the computed results of the maximum absolute error versus the number of source points. Figure 10 shows the relationship of the condition number versus the number of source points for the proposed method and the MFS. It is found that the proposed method remains relatively high accuracy compared to the MFS in this example. The best accuracy can reach the order of while the number of source points is greater than 150. On the other hand, the best accuracy of the MFS can reach only about in the same example.
4. Application of the Proposed Method
4.1. Modeling of Subsurface Flow with Free Surface
The first application under investigation is a free surface seepage problem of a rectangular dam as depicted in Figure 11. The subsurface flow equation is the Laplace equation. The example with the upstream hydraulic head is 24 m, the downstream hydraulic head is 4 m, and the width of the rectangular dam is 16 m. The boundary conditions includes , , , , and , as depicted in Figure 11. In and , the Dirichlet boundary conditions are given as Based on the Bernoulli equation, we neglected the velocity head and the total head or the potential can be written aswhere is the elevation head, is the pressure head, and is the unit weight of fluid. In and , the free surface boundaries are given as overspecified boundary conditions asIn , the no-flow Neumann boundary condition to simulate the imperious boundary is given asSince is unknown a priori which needs to be determined iteratively after the initial guess of the free surface, the proposed method adopted to find the location of free boundary is expressed in the following section.
The subsurface flow with a free surface is a nonlinear problem in which nonlinearities arise from the nonlinear boundary characteristics. Such nonlinearities are handled in the numerical modeling using iterative schemes. Typically, the methods, such as the Picard method or Newton’s method, are iterative in that they approach the solution through a series of steps. In this study, the Picard method is adopted.
There are 16 boundary collocation nodes uniformly distributed in the initial guess of the moving boundary with the spacing of 1 m as shown in Figure 12. Figure 12 shows the computed results using the proposed method. There are 132 iterations to reach the stopping criterion using the Picard scheme. The numerical solutions of free surface were then compared with those obtained from Aitchison et al. [32, 33]. The separation point is the intersection of the free surface and the seepage face. The location of the separation point computed by this study is 13.19 m. It is found that the computed results agree closely with those from other methods.
4.2. Free Surface Seepage Flow through Layered Heterogeneous Geological Media
The previous examples have demonstrated that the proposed method can be used to deal with the subsurface flow with a free surface. Since the appearance of layered soil in heterogeneous geological media is much more common than homogeneous soil in nature, we further adopted the proposed method to deal with the subsurface flow problems of layered heterogeneous geological media using the DDM.
This example under investigation is a rectangular dam in layered soil as depicted in Figure 13. We considered the problem where the upstream hydraulic head is 10 m, the downstream hydraulic head is 2 m, and the height and the width of the rectangular dam are 10 m and 5 m, respectively. The boundary conditions including are shown in Figure 13. At and , the Dirichlet boundary conditions are given as At , , , and , the free surface boundaries are given as overspecified boundary conditions asAt and , the no-flow Neumann boundary condition to simulate the imperious boundary is given asTo deal with the geofluid flow through layered heterogeneous geological media, the domain decomposition method was adopted. The solution continuity or compatibility between different subdomains was assured by remaining equal values of the pressure potential and the flux at the interface between subdomains. For instance, the free surface seepage flow through layered heterogeneous geological media as at and , the flux conservation, and the continuity of pressure potential at the interface between two consecutive layers have to ensure the solution continuity. Accordingly, the following additional boundary conditions must be given:There are two soil layers in this example. The values of the hydraulic conductivity for layer 1 and layer 2 are and , respectively, and and .
In this study, we adopted the nonoverlapping method to deal with the subsurface flow problems of layered soil profiles. The problems on the subdomains are independent, which makes the DDM suitable for describing the layered soil in heterogeneous porous media.
For the modeling of the layered soil, we split the domain into smaller subdomains in which subdomains were intersected only on the interface between soil layers, as shown in Figure 13. For example, there is a problem with two soil layers as shown in Figure 13. The hydraulic conductivities are and for soil layer 1 and soil layer 2, respectively. The boundary and source points were collocated in each subdomain. At the interface, the boundary collocation points on left and right sides coincide with each other. The proposed method was then adopted to ensure that flux conservation and the continuity of pressure potential at the interface between two consecutive layers remain the same.
For the first subdomain, there are a total of 250 boundary collocation nodes where 50 boundary collocation nodes are uniformly distributed in the initial guess of the moving boundary. For the second subdomain, there are also a total of 250 boundary collocation nodes where 50 boundary collocation nodes are uniformly distributed in the initial guess of the moving boundary.
Figure 14 shows the computed results using the proposed method. There are 14 iterations to reach the stopping criterion using the Picard scheme. The numerical solutions of free surface were then compared with those obtained from previous studies [27, 34]. It is found that the computed results agree well with those from other methods.
4.3. Modeling of Three-Dimensional Subsurface Flow Problem
Because the basis function, , is also the particular solution of the Laplace equation in three-dimensional cylindrical coordinate system, it implies that the basis function proposed in this study can also be used to solve the three-dimensional subsurface flow problems. Accordingly, the last example under investigation is a three-dimensional homogenous isotropic steady-state subsurface flow problem. For a three-dimensional simply connected domain enclosed by a boundary as shown in Figure 15, the governing equation is expressed asThe boundary is defined aswhere , , and .
The analytical solution of the problem is given asThe Dirichlet boundary condition is imposed on the boundary by using the analytical solution as shown in (34) for the problem. Figure 15 shows the boundary collocation and the three-dimensional shape of the problem. A sensitivity study was first carried out and could be the satisfactory location of the source points, as shown in Figure 16. Also, to examine the accuracy, we collocated 889 uniformed distributed inner points inside the domain. The maximum absolute error can then be found by evaluating the absolute error for each inner point.
Figure 17 depicted the computed results of the maximum absolute error versus the number of source points. The best accuracy of the proposed method can reach the order of while the number of source points is greater than 350.
5. Conclusions
This study has proposed a novel boundary-type meshless method for modeling geofluid flow in heterogeneous geological media. The numerical solutions of geofluid flow are approximated by a set of particular solutions of the subsurface flow equation which are expressed in terms of sources located outside the domain of the problem. To deal with the subsurface flow problems of heterogeneous geological media, the domain decomposition method was adopted. The validity of the model is established for a number of test problems. Application examples of subsurface flow problems with free surface were also carried out. The fundamental concepts and the construct of the proposed method are addressed in detail. The findings are addressed as follows.
In this study, a pioneering study is based on the collocation Trefftz method and provides a promising solution which integrates the T-Trefftz method and F-Trefftz method for constructing its basis function using one of the negative particular solutions which satisfies the governing equation and allows many source points outside the domain of interest. The proposed method uses the same concept of the source points in the MFS, but the fundamental solutions can be replaced by the negative Trefftz functions. It may release one of the limitations of the MFS in which the fundamental solutions may be difficult to find.
It is well known that the system of linear equations obtained from the Trefftz method may also become an ill-posed system with the higher order of the terms. In this study, the proposed method integrates the collocation Trefftz method and the MFS which approximates the numerical solutions by superpositioning of the negative particular solutions as basis functions expressed in terms of many source points. As a result, only two Trefftz terms were adopted because many source points are allowed for approximating the solution. Meanwhile, the ill-posedness from adopting the higher order terms for the solution with only one source point in the collocation Trefftz method can be mitigated. In addition, results from the validation examples demonstrate that the proposed method may obtain better accuracy than the MFS.
The validity of the model is established for a number of test problems, including the investigation of the basis function using two possible particular solutions and the comparison of the numerical solutions using different particular solutions and the method of fundamental solutions. Application examples of subsurface flow problems with free surface were also carried out. Numerical results demonstrate that the proposed method is highly accurate and computationally efficient. This pioneering study demonstrates that the proposed boundary-type meshless method may be the first successful attempt for solving the subsurface flow with nonlinear free surface in layered heterogeneous geological media which has not been reported in previous studies. Moreover, the application example depicted that the proposed method can be easily applied to the three-dimensional problems.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was partially supported by the Ministry of Science and Technology, Taiwan. The authors thank the Ministry of Science and Technology for the generous financial support.