Geofluids
Volume 2018 (2018), Article ID 9804291, 13 pages
https://doi.org/10.1155/2018/9804291

Research Article

A Novel Boundary-Type Meshless Method for Modeling Geofluid Flow in Heterogeneous Geological Media

Jing-En Xiao, Cheng-Yu Ku, Chih-Yu Liu, and Wei-Chung Yeih

Department of Harbor and River Engineering, National Taiwan Ocean University, Keelung, Taiwan

Correspondence should be addressed to Cheng-Yu Ku; chkst26@email.ntou.edu.tw

Received 3 July 2017; Accepted 18 December 2017; Published 16 January 2018

Academic Editor: Shujun Ye

Copyright © 2018 Jing-En Xiao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

with

where 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 as

where , , 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 as

where 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

Figure 1: Illustration of four different types of domain in the CTM.

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 as

where 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) include

The 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

Figure 2: Illustration of four different types of domain in the MFS.

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 as

where 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 as

where 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 as

The 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 that

where are the constant coefficients to be solved, and is the boundary condition imposed at boundary collocation points. Considering the boundary conditions, we have

where represents the Dirichlet boundary condition; represents the Neumann boundary condition. Applying Dirichlet and Neumann boundary conditions, we obtained

where 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 as

The two-dimensional object boundary under consideration is defined as

where , .

The analytical solution can be found as

The 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.

Figure 3: The collocation of boundary and source points.

Figure 4: The accuracy of the maximum absolute error versus .

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: The accuracy of the maximum absolute error versus the number of source points.

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.

Figure 6: The condition number versus the number of source points.

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 as

The two-dimensional object boundary under consideration is defined as

where , .

The analytical solution can be found as

The 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 7: The collocation of boundary and source points.

Figure 8: The accuracy of the maximum absolute error versus .

Figure 9: The accuracy of the maximum absolute error versus the number of source points.

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.

Figure 10: The condition number versus the number of source points.

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 as

where 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 as

In , the no-flow Neumann boundary condition to simulate the imperious boundary is given as

Since 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.

Figure 11: Subsurface flow with free surface through a rectangular dam.

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.

Figure 12: Result comparison of the computed free surface for a rectangular dam.

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 as

At and , the no-flow Neumann boundary condition to simulate the imperious boundary is given as

To 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 .

Figure 13: The collocation of boundary, source points (a) and configuration of boundary condition (b).

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.

Figure 14: Comparison of free surface for a rectangular dam in layered heterogeneous geological media.

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 as

The boundary is defined as

where , , and .

Figure 15: The boundary collocation points of three-dimensional subsurface flow problem.

The analytical solution of the problem is given as

The 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 16: The accuracy of the maximum absolute error versus .

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.

Figure 17: The accuracy of the maximum absolute error versus the number of source points.

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.

References

R. Allen and S. Sun, “Computing and comparing effective properties for flow and transport in computer-generated porous media,” Geofluids, vol. 2017, Article ID 4517259, 2017.

M. H. Bazyar and A. Graili, “A practical and efficient numerical scheme for the analysis of steady state unconfined seepage flows,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 36, no. 16, pp. 1793–1812, 2012.

J. T. Chen, H.-K. Hong, and S. W. Chyuan, “Boundary element analysis and design in seepage problems using dual integral formulation,” Finite Elements in Analysis and Design, vol. 17, no. 1, pp. 1–20, 1994.

J. Li and D. Brown, “Upscaled lattice boltzmann method for simulations of flows in heterogeneous porous media,” Geofluids, vol. 2017, Article ID 1740693, 2017.

J. A. Liggett and L. F. Liu, The Boundary Integral Equation Method for Porous Media Flow, George Allen & Unwin, London, UK, 1983.

C.-Y. Ku, C.-L. Kuo, C.-M. Fan, C.-S. Liu, and P.-C. Guan, “Numerical solution of three-dimensional Laplacian problems using the multiple scale Trefftz method,” Engineering Analysis with Boundary Elements, vol. 50, pp. 157–168, 2015.

W. Yeih, I.-Y. Chan, C.-Y. Ku, and C.-M. Fan, “Solving the inverse cauchy problem of the laplace equation using the method of fundamental solutions and the exponentially convergent scalar homotopy algorithm (ECSHA),” Journal of Marine Science and Technology (Taiwan), vol. 23, no. 2, pp. 162–171, 2015.

C.-Y. Ku, “On solving three-dimensional laplacian problems in a multiply connected domain using the multiple scale trefftz method,” CMES: Computer Modeling in Engineering & Sciences, vol. 98, no. 5, pp. 509–541, 2014.

C.-Y. Ku, “A novel method for solving ill-conditioned systems of linear equations with extreme physical property contrasts,” CMES. Computer Modeling in Engineering & Sciences, vol. 96, no. 6, pp. 409–434, 2013.

J.-S. Chen, M. Hillman, and S.-W. Chi, “Meshfree methods: Progress made after 20 years,” Journal of Engineering Mechanics, vol. 143, no. 4, Article ID 04017001, 2017.

C.-M. Fan, H.-F. Chan, C.-L. Kuo, and W. Yeih, “Numerical solutions of boundary detection problems using modified collocation Trefftz method and exponentially convergent scalar homotopy algorithm,” Engineering Analysis with Boundary Elements, vol. 36, no. 1, pp. 2–8, 2012.

E. Kita and N. Kamiya, “Trefftz method: an overview,” Advances in Engineering Software, vol. 24, no. 1-3, pp. 3–12, 1995.

E. Kita, N. Kamiya, and T. Iio, “Application of a direct Trefftz method with domain decomposition to 2D potential problems,” Engineering Analysis with Boundary Elements, vol. 23, no. 7, pp. 539–548, 1999.

Z.-C. Li, T.-T. Lu, H.-Y. Hu, and A. H.-D. Cheng, Trefftz and Collocation Methods, WIT Press, Southampton, UK, 2008.

C.-S. Liu, “A modified Trefftz method for two-dimensional Laplace equation considering the domain's characteristic length,” CMES. Computer Modeling in Engineering & Sciences, vol. 21, no. 1, pp. 53–65, 2007.

E. Trefftz, “Ein Gegenstück zum Ritzschen Verfahren,” in Proceedings of the 2nd International Congress for Applied Mechanics, pp. 131–137, 1926.

V. D. Kupradze and M. A. Aleksidze, “The method of functional equations for the approximate solution of certain boundary value problems,” USSR Computational Mathematics and Mathematical Physics, vol. 4, no. 4, pp. 82–126, 1964.

M. Katsurada and H. Okamoto, “The collocation points of the fundamental solution method for the potential problem,” Computers & Mathematics with Applications, vol. 31, no. 1, pp. 123–137, 1996.

J. T. Chen, C. S. Wu, Y. T. Lee, and K. H. Chen, “On the equivalence of the Trefftz method and method of fundamental solutions for Laplace and biharmonic equations,” Computers & Mathematics with Applications. An International Journal, vol. 53, no. 6, pp. 851–879, 2007.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

W. K. Liu, S. Jun, and Y. F. Zhang, “Reproducing kernel particle methods,” International Journal for Numerical Methods in Fluids, vol. 20, no. 8-9, pp. 1081–1106, 1995.

J.-S. Chen, C. Pan, C.-T. Wu, and W. K. Liu, “Reproducing kernel particle methods for large deformation analysis of non-linear structures,” Computer Methods Applied Mechanics and Engineering, vol. 139, no. 1-4, pp. 195–227, 1996.

T. Zhu, J. Zhang, and S. N. Atluri, “Meshless numerical method based on the local boundary integral equation (LBIE) to solve linear and non-linear boundary value problems,” Engineering Analysis with Boundary Elements, vol. 23, no. 5, pp. 375–389, 1999.

J. Sladek, V. Sladek, and S. N. Atluri, “Local boundary integral equation (LBIE) method for solving problems of elasticity with nonhomogeneous material properties,” Computational Mechanics, vol. 24, no. 6, pp. 456–462, 2000.

H. Lin and S. N. Atluri, “Meshless local Petrov-Galerkin (MLPG) method for convection-diffusion problems,” CMES-Computer Modelingin Engineering & Sciences, vol. 1, no. 2, pp. 45–60, 2000.

D.-L. Young, C.-M. Fan, C.-C. Tsai, and C.-W. Chen, “The method of fundamental solutions and domain decomposition method for degenerate seepage flownet problems,” Journal of the Chinese Institute of Engineers, Transactions of the Chinese Institute of Engineers,Series A/Chung-kuo Kung Ch'eng Hsuch K'an, vol. 29, no. 1, pp. 63–73, 2006.

S. J. Lacy and J. H. Prevost, “Flow through porous media: A procedure for locating the free surface,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 11, no. 6, pp. 585–601, 1987.

S. Mehl, “Use of Picard and Newton iteration for solving nonlinear ground water flow equations,” Groundwater, vol. 44, no. 4, pp. 583–594, 2006.

W. Yeih, C.-S. Liu, C.-L. Kuo, and S. N. Atluri, “On solving the direct/inverse cauchy problems of laplace equation in a multiply connected domain, using the generalized multiple-source-point boundary-collocation Trefftz method & characteristic lengths,” CMC: Computer, Materials, & Continua, vol. 17, no. 3, pp. 275–302, 2010.

C.-Y. Ku, J.-E. Xiao, C.-Y. Liu, and W. Yeih, “On the accuracy of the collocation Trefftz method for solving two- and three-dimensional heat equations,” Numerical Heat Transfer, Part B: Fundamentals, vol. 69, no. 4, pp. 334–350, 2016.

C. S. Chen, A. Karageorghis, and Y. Li, “On choosing the location of the sources in the MFS,” Numerical Algorithms, vol. 72, no. 1, pp. 107–130, 2016.

J. Aitchison, “Numerical treatment of a singularity in a free boundary problem,” Proceedings of the Royal Society A Mathematical, Physical and Engineering Sciences, vol. 330, pp. 573–580, 1972.

J.-T. Chen, C.-C. Hsiao, Y.-P. Chiu, and Y.-T. Lee, “Study of free-surface seepage problems using hypersingular equations,” Communications in Numerical Methods in Engineering, vol. 23, no. 8, pp. 755–769, 2007.

M. X. Wu, L. Z. Yang, and T. Yu, “Simulation procedure of unconfined seepage with an inner seepage face in a heterogeneous field,” Science China Physics, Mechanics & Astronomy, vol. 56, no. 6, pp. 1139–1147, 2013.