Next Article in Journal
Avoiding Conflicts between Future Freshwater Algae Production and Water Scarcity in the United States at the Energy-Water Nexus
Next Article in Special Issue
Assess Effectiveness of Salt Removal by a Subsurface Drainage with Bundled Crop Straws in Coastal Saline Soil Using HYDRUS-3D
Previous Article in Journal
Influence of Check Dams on Flood and Erosion Dynamic Processes of a Small Watershed in the Loss Plateau
Previous Article in Special Issue
Improved Solutions to the Linearized Boussinesq Equation with Temporally Varied Rainfall Recharge for a Sloping Aquifer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On Solving Nonlinear Moving Boundary Problems with Heterogeneity Using the Collocation Meshless Method

1
Department of Harbor and River Engineering, National Taiwan Ocean University, Keelung 20224, Taiwan
2
Center of Excellence for Ocean Engineering, National Taiwan Ocean University, Keelung 20224, Taiwan
*
Author to whom correspondence should be addressed.
Water 2019, 11(4), 835; https://doi.org/10.3390/w11040835
Submission received: 22 March 2019 / Revised: 16 April 2019 / Accepted: 17 April 2019 / Published: 20 April 2019

Abstract

:
In this article, a solution to nonlinear moving boundary problems in heterogeneous geological media using the meshless method is proposed. The free surface flow is a moving boundary problem governed by Laplace equation but has nonlinear boundary conditions. We adopt the collocation Trefftz method (CTM) to approximate the solution using Trefftz base functions, satisfying the Laplace equation. An iterative scheme in conjunction with the CTM for finding the phreatic line with over–specified nonlinear moving boundary conditions is developed. To deal with flow in the layered heterogeneous soil, the domain decomposition method is used so that the hydraulic conductivity in each subdomain can be different. The method proposed in this study is verified by several numerical examples. The results indicate the advantages of the collocation meshless method such as high accuracy and that only the surface of the problem domain needs to be discretized. Moreover, it is advantageous for solving nonlinear moving boundary problems with heterogeneity with extreme contrasts in the permeability coefficient.

1. Introduction

The free surface flow is a moving boundary problem governed by the Laplace equation but has nonlinear boundary conditions. The study of free surface seepage problem plays a crucial role in the analysis of hydraulic engineering. In the design of embankment, earth dams and rock–fill dams, finding the position of the moving boundary is of importance [1,2]. The solution of the Laplace governing equation may be carried out by solving the boundary value problem. Because the solution of the free surface seepage flow is nonlinear, iterative techniques are often required in the solution process for matching the over–specified boundary conditions. Various mesh–based numerical methods [3,4,5,6,7,8,9,10,11,12,13] have been used for the analysis of free surface flow. For the mesh–based methods, automatic grid regeneration [6] is commonly used to solve the free surface problems in mesh–based approaches. However, convergence problems are often raised due to the changing of element shapes and types in the process of the mesh generation. While the complexity of the boundary conditions is considered, mesh–based methods may become unstable since the automatic grid regeneration is likely to generate distorted meshes.
Recently, meshless methods have attracted much attention to solve free surface seepage problems [14]. Compared to mesh–based methods, the discretization of the domain for meshless methods is relatively simple because only arbitrary collocation points need to be placed in the physical domain without using any elements. If the basic functions satisfy the governing equation, the collocation points may be conducted only on the boundary. Accordingly, the meshless method has advantages with problems involving complex geometry [15]. Among several meshless methods, the collocation Trefftz method (CTM) may be regarded as an attractive boundary–type meshless method [16]. In the past, study of the Trefftz method was less widespread because the ill-posedness of the CTM limits the applications of the method. However, using the characteristic length [17,18,19], the CTM has been adopted to obtain accurate solutions for solving the Laplace governing equation in three–dimensions enclosed by simply and multiply connected domains [20].
The Trefftz method is first proposed by the German mathematician Erich Trefftz. Later, the CTM [21,22,23,24,25,26,27,28,29,30] is commonly used for solving partial differential equations. Since the CTM is categorized into the boundary–type meshless method, it approximates the solutions of the governing equation using the Trefftz basis functions where the solutions are described as the assembly of the Trefftz functions [31]. The CTM requires the evaluation of the coefficients in which they may be obtained by solving the linear simultaneous equations assembled by using the boundary conditions at a number of collocation points. Applications of the CTM such as Laplace and modified Helmholtz equations [32,33] and the problem of boundary detection [34] has been studied. Due to the complexity, applications of the CTM are most limited to the homogeneous problems. In addition, the study of nonlinear moving boundary problems with heterogeneity using the CTM has not been reported yet.
This paper presents the study on solving nonlinear moving boundary problems in heterogeneous geological media using the CTM. The free surface flow is a moving boundary problem governed by Laplace equation but has nonlinear boundary conditions. We adopt the CTM to approximate the solution using Trefftz base functions satisfying the Laplace equation. An iterative scheme in conjunction with the CTM for finding the phreatic line with over–specified nonlinear moving boundary conditions is developed. To deal with flow in the layered heterogeneous soil, the domain decomposition method is used so that the hydraulic conductivity in each subdomain can be different. The method proposed in this study is verified by several numerical examples. The formulation of the proposed method is described as follows.

2. Governing Equation and Boundary Conditions

The two-dimensional Laplace equation used to represent flow through a homogenous rectangular dam is expressed as
Δ φ = 0   in   Ω
and
φ = g   on   Γ D ,
φ n = f   on   Γ N ,
where φ is the head, Δ is the Laplacian, Ω represents the domain boundary of the problem, g and f denote the Dirichlet and Neumann boundary conditions, respectively. n denotes the normal vector. Γ D and Γ N denote the Dirichlet and Neumann boundary conditions.
The boundary conditions of the rectangular dam with a moving boundary, as depicted in Figure 1a, can be presented by Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 . The Dirichlet boundary conditions are imposed on the Γ 2 and Γ 5 , respectively.
φ = H 2   on   Γ 2 ,
φ = H 1   on   Γ 5 ,
where H 2 denotes the downstream elevation, H 1 denotes the upstream elevation. Neglecting the velocity head, the head is expressed as
φ = Y ( x ) + p γ ,
where Y ( x ) denotes the height above the sea level, p denotes the pore water pressure and γ denotes the water unit weight. On Γ 4 , the over–specified moving boundary conditions are given as
φ n = 0 ,   φ = Y ( x )   on   Γ 4 .
On Γ 3 , the seepage face boundary condition is depicted as
φ = Y ( x )   on   Γ 3 .
φ = Y ( x ) is unknown which can be solved by using the iterative scheme. On Γ 1 , the no–flow condition is given as
φ n = 0   on   Γ 1 .

3. The Collocation Trefftz Method

The Laplace equation in polar coordinate system is expressed as
2 φ r 2 + 1 r φ r + 1 r 2 2 φ θ 2 = 0 ,
where the radial coordinate is denoted by r and the angular coordinate is denoted by θ . The solution of the Laplace governing equation is approximated by using the Trefftz basis functions satisfying the governing equation, as shown in Equation (10). The Trefftz basis functions are obtained by finding the general solutions using the separation of variables method [35]. The Trefftz basis functions can be found to solve problems in a simply connected domain, as shown in Figure 2.

3.1. Formulation of T-Complete Basis Functions

We may apply the separation of variables [36]. The solution may be in the following form:
φ ( r , θ ) = U ( r ) V ( θ ) .
For simplicity, we let
U = d U ( r ) d r ,   U = d 2 U ( r ) d r 2   and   V = d 2 V ( θ ) d θ 2 .
Then, Equation (10) can be rewritten as follows.
U V + 1 r U V + 1 r 2 U V = 0 .
We divide U ( r ) V ( θ ) on both sides in the above equation and the equation can be rewritten as two differential equations as follows.
r 2 U U + r U U = λ r 2 ,
V V = λ .
Using the constant v to ensure positive or negative constants, we have λ = 0 , λ = v 2 and λ = v 2 . Considering the first scenario λ = 0 , we obtain the solutions as follows.
V = D 1 θ + D 2 ,
U = D 3 ln r + D 4 ,
where D 1 , D 2 , D 3 and D 4 are constants. Using the boundary conditions of V ( r , 0 ) = V ( r , 2 π ) , we may find that D 1 = 0 . Substituting Equations (16) and (17) into Equation (11), we have
φ = a 0 + b 0 ln r ,
where a 0 and b 0 denote the coefficients. Considering the second scenario, λ = v 2 , we obtain the following solutions.
V = D 5 cos ( ν θ ) + D 6 sin ( ν θ ) ,
U = D 7 r ν + D 8 r ν ,
where D 5 , D 6 , D 7 and D 8 are the coefficients. Inserting the above equations into Equation (11), we obtain
φ = a r ν cos ( ν θ ) + b r ν sin ( ν θ ) + c r ν cos ( ν θ ) + d r ν sin ( ν θ ) ,
where a , b , c and d denote the coefficients. Then, we may consider the last scenario λ = v 2 . Since there is not able to find any non–zero periodic solutions of differential system for U ( r ) , we may only find V ( θ ) = 0 . Collecting all the solutions from the above results, the linearly independent solutions to Laplace equation can be obtained as follows.
{ 1 , ln r , r ν cos ( ν θ ) , r ν sin ( ν θ ) , r ν cos ( ν θ ) , r ν sin ( ν θ ) } .
The Trefftz basis functions in a simply connected domain are as follows.
{ 1 r ν cos ( ν θ ) , r ν sin ( ν θ ) } .
In the numerical analysis, we approach the general solution in the form of infinite series of the Laplace equation in a simply connected domain by using a finite number of m . As a result, Equation (23) can be rewritten as
φ = a 0 + v = 1 m { a v r v cos ( v θ ) + c v r v sin ( v θ ) } ,
where m represents the terms of the Trefftz order. The above Equation (24) can be used to match the Dirichlet boundary condition. We may also need to consider the Neumann boundary conditions as follows.
φ n = φ n   on   Γ N .
Equation (25) can be rewritten as
φ n = φ n ,
where is the gradient and n = ( n x , n y ) denotes the normal vector. Equation (25) can then be written as
φ n = φ x n x + φ y n y ,
where
φ x = φ r r x + φ θ θ x ,   φ y = φ r r y + φ θ θ y .
Using Equation (24), we may find the derivatives of φ / r and φ / θ as follows.
φ r = ν = 1 m a ν ν r ν 1 cos ( ν θ ) + c ν ν r ν 1 sin ( ν θ ) ,
φ θ = ν = 1 m a ν ( ν ) r ν sin ( ν θ ) + c ν ν r ν cos ( ν θ ) .
Using Equations (29) and (30), Equation (27) leads to
φ n = v = 1 m { a v [ ( v r v 1 cos ( v θ ) cos θ + ( v ) r v sin ( v θ ) sin θ r ) n x + ( v r v 1 cos ( v θ ) sin θ + ( v ) r v sin ( v θ ) cos θ r ) n y ] + c v [ ( v r v 1 sin ( v θ ) cos θ + v r v cos ( v θ ) sin θ r ) n x + ( v r v 1 sin ( v θ ) sin θ + v r v cos ( v θ ) cos θ r ) n y ] }

3.2. The Characteristic Length

The characteristic length plays a crucial role in controlling the proposed numerical approach in a stable way. Because the matrix assembled with Trefftz trial functions is a full matrix, the resultant system of linear equations may be ill–posed [17,18]. The accuracy of the results from the Trefftz method depends sensitively on the order of the T-complete basis functions. Besides, the numerical solution may be unstable. Related to the CTM for solving two-dimensional Laplacian problems, Liu [17] proposed a characteristic length to mitigate the problems of the ill-posedness for the system of linear equations. Applying Dirichlet boundary condition, we obtain
φ = a 0 + v = 1 m [ a v ( r R ) v cos ( v θ ) + c v ( r R ) v sin ( v θ ) ] .
Using the CTM, we obtain the approximation solution of the Laplace equation as follows.
φ ( x ) = v = 1 m b v T v ( x ) ,
where b v = [ a 0 a v c v ] , T v = [ 1 ( r / R ) v cos ( v θ ) ( r / R ) v sin ( v θ ) ] , x is the coordinate of the collocation points and x Ω . Applying the Neumann boundary condition, we may obtain the following equations for simply connected domain using the characteristic length.
φ n = v = 1 m { a v [ ( v ( 1 R ) v r v 1 cos ( v θ ) cos θ + ( v ) ( 1 R ) v r v sin ( v θ ) sin θ r ) n x + ( v ( 1 R ) v r v 1 cos ( v θ ) sin θ + ( v ) ( 1 R ) v r v sin ( v θ ) cos θ r ) n y ] + c v [ ( v ( 1 R ) v r v 1 sin ( v θ ) cos θ + v ( 1 R ) v r v cos ( v θ ) sin θ r ) n x + ( v ( 1 R ) v r v 1 sin ( v θ ) sin θ + v ( 1 R ) v r v cos ( v θ ) cos θ r ) n y ] }
To mitigate the ill-posedness, the characteristic length [19], R , is adopted and is expressed as
R = 1.5 × m a x i m u m ( r ) ,
where m a x i m u m ( r ) denotes the maximum radial distance in the problem domain. After adopting the characteristic length in our numerical model, the ill-posed phenomenon is greatly reduced, and the accurate numerical solutions can be obtained. Collocating the numerical expansion from Equations (32) and (34) at boundary collocation points to match the given boundary conditions, we may obtain the following equation.
T b = B ,
T = [ 1 ( r 1 / R ) cos ( θ 1 ) ( r 1 / R ) sin ( θ 1 ) ( r 1 / R ) m cos ( m θ 1 ) ( r 1 / R ) m sin ( m θ 1 ) 1 ( r 2 / R ) cos ( θ 2 ) ( r 2 / R ) sin ( θ 2 ) ( r 2 / R ) m cos ( m θ 2 ) ( r 2 / R ) m sin ( m θ 2 ) 1 ( r i / R ) cos ( θ i ) ( r i / R ) sin ( θ i ) ( r i / R ) m cos ( m θ i ) ( r i / R ) m sin ( m θ i ) 0 N 1 , ν = 1 a 1 N 1 , ν = 1 c 1 N 1 , ν = m a ν N 1 , ν = m c ν 0 N 2 , ν = 1 a 1 N 2 , ν = 1 c 1 N 2 , ν = m a ν N 2 , ν = m c ν 0 N j , ν = 1 a 1 N j , ν = 1 c 1 N j , ν = m a ν N j , ν = m c ν ] , b = [ a 0 a 1 c 1 a m c m ] ,   B = [ g 1 g 2 g i f 1 f 2 f j ] ,
N j , ν a ν = v = 1 m [ ( v ( 1 R ) v r j v 1 cos ( v θ ) cos θ + ( v ) ( 1 R ) v r j v sin ( v θ ) sin θ r j ) n x + ( v ( 1 R ) v r j v 1 cos ( v θ ) sin θ + ( v ) ( 1 R ) v r j v sin ( v θ ) cos θ r j ) n y ]
N j , ν c ν = v = 1 m [ ( v ( 1 R ) v r j v 1 sin ( v θ ) cos θ + v ( 1 R ) v r j v cos ( v θ ) sin θ r j ) n x + ( v ( 1 R ) v r j v 1 sin ( v θ ) sin θ + v ( 1 R ) v r j v cos ( v θ ) cos θ r j ) n y ]
where T represents a l × M matrix, M = 2 m + 1 , b represents a M × 1 vector of unknown coefficients, B denotes a vector (size of l × 1 ) of given functions at boundary points, l represents the number of the boundary points and M represents the terms of the Trefftz order. i l and j l in which i and j are the number of boundary points for Dirichlet and Neumann boundary conditions, respectively. g 1 ,   g 2 ,   ,   g i and f 1 ,   f 2 ,   ,   f j denote the boundary values for Dirichlet and Neumann boundary conditions, respectively.
In this article, we adopt the domain decomposition method (DDM) [37,38] to solve the nonlinear moving boundary problems in heterogeneous geological media. The DDM is commonly used to solve the problem with different physical characteristics in each subdomain. We first split the domain into two subdomains which are intersected only at the interface. Hence, each subdomain can be regarded as an independent soil layer with its own hydraulic conductivity. At the interface, the flux and the head must satisfy the continuity condition. For instance, we consider a rectangular domain, Ω , which can be split into two intersected subdomains, Ω 1 and Ω 2 . Figure 3 shows that the rectangular domain is divided into Γ 1 , Γ 2 , …, Γ 8 sub boundaries where Γ 1 , Γ 2 , …, Γ 4 and Γ 5 , Γ 6 , …, Γ 8 are sub boundaries of subdomains Ω1. and Ω 2 , respectively. At the interface, the sub boundaries, Γ 2 and Γ 6 , are overlapped at the same location. Therefore, additional boundary conditions are imposed on the boundary points to ensure the flux and the head at the interface must be the same.
φ | Γ 2 = φ | Γ 6 ,   φ n | Γ 2 = φ n | Γ 6 .
Matching all given boundary conditions, we may obtain a system of linear equations as
T D b D = B D ,
T D = [ T Ω 1 0 Ω 2 T I | Γ 2 T I | Γ 6 0 Ω 1 T Ω 2 ] ,   b D = [ b Ω 1 b Ω 2 ] ,   B D = [ B Ω 1 B I B Ω 2 ] ,
where T Ω 1 with the size of l 1 × M 1 and T Ω 2 with the size of l 2 × M 2 are the T matrix shown in Equation (37) for Ω 1 and Ω 2 , respectively. l 1 and l 2 are the number of boundary points; M 1 and M 2 are the T-complete basis function order for Ω 1 and Ω 2 , respectively. T I | Γ 2 of the boundary Γ 2 with the size of l I × M 1 and T I | Γ 6 of the boundary Γ 6 with the size of l I × M 2 are matrices at the interface. l I represents the boundary point number at the interface, 0 Ω 1 and 0 Ω 2 are matrices which all values are zero with the size of l 2 × M 1 and l 1 × M 2 , respectively. b Ω 1 denotes a M 1 × 1 vector of unknown coefficients of Ω 1 , b Ω 2 denotes a M 2 × 1 vector of unknown coefficients of Ω 2 . B Ω 1 and B Ω 2 denote vectors of given functions at boundary points of Ω 1 and Ω 2 , respectively. B I = [ 0 g 0 f ] T , 0 g and 0 f are vectors which all values are zero with the size of l I × 1 . The total head can be determined by collocating the inner points within subdomains, Ω 1 and Ω 2 . Consequently, the value of the total head, φ , can then be approximated by using Equation (33).

3.3. The Iterative Scheme for Solving Free Surface

The nonlinearity of the moving surface flow is caused by the nonlinear boundary conditions. For solving free surface problems with nonlinear boundary conditions, the iterative scheme must be used in the numerical modeling. Due to the difficulty of the computation of the Jacobian matrix for Newton’s method, the Picard scheme is adopted in this study. Applying boundary conditions, we obtain
φ ( x k ) v = 1 m b v T v ( x k ) = g ( x k ) ,
φ ( x k ) n v = 1 m b v n T v ( x k ) = f ( x k ) ,
where k denotes the index of the collocation points on the free surface to be updated. The head at J t h iteration is given as
φ J ( x k ) v = 1 m b v J T v ( x k J ) ,
φ J ( x k ) n = φ J n v = 1 m b v J n T v ( x k J ) ,
where J denotes the number of iteration steps. We may obtain the following iterative equation [39].
φ J ( x k ) = φ J 1 ( x k ) + β ( φ J ( x k ) φ J 1 ( x k ) ) ,
where β is the under–relaxation factor and φ J ( x k ) is the head to be updated. The value of β is ranging from 0 to 1. The iterative process begins from the first guess values for the moving surface and ceases while the stopping condition expressed in the following equation is achieved.
ε = k = 1 n i ( φ J ( x k ) φ J 1 ( x k ) ) 2 k = 1 n i ( φ J 1 ( x k ) ) 2 10 4 ,
where n i is the collocation point number on the free surface.

4. Validation Examples

4.1. Laminar Flow around a Cylinder

The first example under consideration is a laminar flow around a cylinder. The dimensions of the problem are depicted in Figure 4a. The radius of the cylinder at the center is 1 m. Because the geometry of the problem is clearly symmetric, we considered the half symmetry model. For a two-dimensional domain, Ω , enclosed by a boundary, the Laplace equation is expressed as
Δ φ = 0   in   Ω
The Dirichlet and Neumann data are applied on the domain boundary, Γ , including Γ 1 , Γ 2 ,…, Γ 6 , as shown in Figure 4a. On Γ 2 and Γ 4 , the Dirichlet data are ϕ = 4 m on Γ 2 and φ = 6 m on Γ 4 . On Γ 1 , Γ 3 , Γ 5 and Γ 6 , the no–flow Neumann boundary data are given as follows.
φ n = 0 .
The order of the Trefftz basis functions, m , is 150. Totally, 900 collocation points including boundary points and sources are uniformly placed on the domain boundary, as shown in Figure 4a. In order to examine accuracy of the proposed method, 7786 inner points are collocated within the domain boundary. The computed results are validated with SVFLUX [40] which is a finite element seepage analysis program. Figure 4b shows the finite element mesh of SVFLUX. The results of the computed head are compared with the exact solution, as shown in Figure 5. It is found that the numerical solutions agree very well with those of the SVFLUX.

4.2. Nonlinear Moving Surface through a Rectangular Dam

The second example is a nonlinear moving surface through a rectangular dam, as shown in Figure 1a. This problem can be viewed as a benchmark problem in geotechnical engineering for finding the position of the moving boundary [3,11,14]. The upstream and downstream heads are 24 m and 4 m, respectively. The dimensions of the problem are depicted in Figure 6. This rectangular dam is assembled with five boundary lines, including Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 , as depicted in Figure 1b.
The order of the Trefftz basis functions, m , is 75. We collocate 120 collocation points on moving boundary, as depicted in Figure 6. Since the process for finding the position of the nonlinear free surface is regarded as an inverse problem, the location of the separation point may also need to be examined. In this study, the initial guess of the separation point is at 14.2 m.
The numerical solutions of the free surface are shown in Figure 6 and the results are compared with those from other methods, such as Aitchison [3], Chen et al. [11] and Xiao et al. [14]. Figure 7 shows that for solving the free surface the number of iteration step is 112 to reach the stopping criterion using the proposed iteration scheme. To further explore the accuracy of the computed results, we compare the computed location of the separation point with those from other numerical methods [3,11,14]. As depicted in Table 1, the position of the free surface is almost identical with other methods. The results of the separation point calculated by three different methods [3,11,14] are 12.79, 12.68 and 12.84 m, respectively. It is found that the location of the separation point by using the CTM is 12.85 m which is close to that from previous studies.

4.3. Nonlinear Moving Surface through an Earth Dam

The third example is a nonlinear moving surface through an earth dam, as shown in Figure 8. The upstream and downstream hear are 18 m and 8 m, respectively. The dimensions of the problem are depicted in Figure 9a. The boundaries of the earth dam include Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 in which Γ 2 and Γ 5 are the downstream and upstream boundaries, as depicted in Figure 8. The boundary values are given as
H 1 = 18   m   on   Γ 5   and   H 2 = 8   m   on   Γ 2 .
Γ 1 is the bottom of the dam where the no–flow condition is given as
φ n = 0   on   Γ 1 .
On Γ 3 , the seepage face boundary condition is as follows.
φ = Y ( x )   on   Γ 3 .
On Γ 4 , the over–specified moving boundary conditions are as follows.
φ n = 0 ,   φ = Y ( x )   on   Γ 4 .
For this example, the Trefftz basis function order, m , is 150. 750 points are collocated on the boundary. The first guess of the free surface is depicted in Figure 9a. The numerical solutions of the free surface are compared with those from other methods [8,14], as shown in Figure 9b. The results illustrate that the computed results are almost identical with other methods.

4.4. Flow through Two Layered Soils

The modeling of two-dimensional heterogeneous isotropic subsurface flow in two layered soils is depicted in Figure 3. The coefficients of the permeability with extreme contrasts for two different soils, k 1 and k 2 , are adopted in which the k 1 = 10 1 and k 2 = 10 15 cm/s. The analytical solution of this example as shown in follows.
φ = φ 2 φ 1 L 1 x + φ 1 ,   0 x L 1 ,
φ = φ 3 φ 2 L 2 x + ( φ 2 φ 3 φ 2 L 2 L 1 ) ,   L 1 x L 1 + L 2 ,
φ 2 = k 1 φ 1 L 2 + k 2 φ 3 L 1 k 1 L 2 + k 2 L 1 ,
where L 1 is the width of the layer 1, L 2 is the width of the layer 2 and φ 2 is the head at the interface. The domain is split into two sub-domains which are simply connected. For each sub-domain, 112 boundary points are uniformly collocated. Figure 3 depicts that the rectangular domain boundary is split into eight sub–boundaries: Γ 1 , Γ 2 ,…, Γ 8 . At the interface, we have the following additional boundary conditions.
φ | Γ 2 = φ | Γ 6 ,   φ n | Γ 2 = φ n | Γ 6 .
Totally, 1520 interior points are collocated within the domain to approximate the head for examining the accuracy. The results of the computed head are compared with the exact solution, as shown in Figure 10. In addition, the accuracy of the results can reach to the order of 10 11 , as shown in Figure 11.

4.5. Nonlinear Moving Surface through a Zoned–Earth Dam

The last case is a nonlinear moving surface through a zoned–earth dam, as shown in Figure 12. For the zoned–earth dam, the upstream and downstream heads are 18 m and 2 m, respectively. The dimensions of the problem are depicted in Figure 12. The values of the permeability for the upstream shell, the core and the downstream shell are 1.43 × 10 4 , 1.43 × 10 5 and 1.43 × 10 4 cm/s, respectively.
For the zoned–earth dam, we divide the domain into three smaller sub-regions, as shown in Figure 13. For each sub-region, we collocate 400, 216 and 191 points on the sub-boundaries for the first, second and third sub-regions, respectively. Besides, we place 50, 66 and 66 collocation points on the moving boundaries, respectively. To validate the results, we compare the computed free surface with that from the finite difference seepage analysis program SEEP2D [41], as shown in Figure 14. It is found that the numerical results agree well with those obtained from the SEEP2D.

5. Discussion

In this article, the CTM is adopted to solve the nonlinear moving boundary problems in layered heterogeneous media. Because of the characteristics of the non-linearity, solving nonlinear moving boundary problems with a moving surface remains a challenge. For modeling moving surface problems with nonlinear boundary conditions, the iterative scheme must be used. The sophisticated automatic mesh generation may be required using conventional mesh-based approaches. In addition, the complicated remeshing process in the iterative scheme may lead to problems of the convergence. In this study, we just need to place the boundary points on the domain boundary. Furthermore, only boundary nodes are required to adjust during the iteration process for find the moving boundary. Comparing with the traditional numerical methods, the proposed method is highly accurate. Therefore, the proposed method is advantageous for the nonlinear moving boundary analysis with a phreatic line.
To solve the flow through the layered soils, we adopt the CTM in conjunction with the DDM so that the hydraulic conductivity in each subdomain can be different. To verify the proposed method, numerical examples with nonlinear moving boundary are conducted. Besides, free surface flow through a zoned–earth dam is also carried out. The comparison of the results using the DDM with the exact solution depicts that the CTM with the use of the DDM can reach the accuracy in the order of 10 11 . Although the CTM may be regarded as an attractive boundary–type meshless method, limitations for the applications may still remain due to the ill-posedness of the method.

6. Conclusions

This paper presents a study on solving nonlinear moving boundary problems in heterogeneous soils using the CTM. The method is verified by several numerical examples. Application examples are also carried out. The findings are concluded.
The appearance of heterogeneous soils is often found in free surface flow problems. In the past, the CTM is usually limited to linear, homogeneous problems. In this article, we propose a novel idea for solving nonlinear moving boundary problems in layered heterogeneous soils using the collocation meshless method. The results show that the proposed method can be used to deal with nonlinear moving boundary problems in heterogeneity with large permeability contrasts. The method proposed in this study is capable of solving nonlinear moving boundary problems in layered heterogeneous media. However, it is still limited to the layered or zoned porous media in which the porous medium is still homogenous in each zoned domain. In addition, the proposed method can only be applied in saturated soils.

Author Contributions

Conceptualization, C.-Y.K.; Data curation, C.-Y.L.; Methodology, J.-E.X.; Validation, J.-E.X.; Writing—original draft, C.-Y.K.; Writing—review & editing, C.-Y.K. and J.-E.X.

Funding

This research was funded by Ministry of Science and Technology of the Republic of China under grant MOST 108-2119-M-019-001.

Acknowledgments

The authors thank the Ministry of Science and Technology for the generous support. The first author would like to express his gratitude to former student, Yi-Wun Chen, for her assistance of this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

References

  1. Zheng, H.; Liu, F.; Li, C.G. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method. Appl. Math. Model. 2015, 39, 794–808. [Google Scholar] [CrossRef]
  2. Wang, Y.; Hu, M.; Zhou, Q.; Rutqvist, J. A new second-order numerical manifold method model with an efficient scheme for analyzing free surface flow with inner drains. Appl. Math. Model. 2015, 40, 1427–1445. [Google Scholar] [CrossRef]
  3. Aitchison, J. Numerical treatment of a singularity in a free boundary problem. Proc. R. Soc. Lond. 1972, 330, 573–580. [Google Scholar] [CrossRef]
  4. Fukuchi, T. Numerical analyses of steady-state seepage problems using the interpolation finite difference method. Soils Found. 2016, 56, 608–626. [Google Scholar] [CrossRef]
  5. Darbandi, M.; Torabi, S.O.; Saadat, M.; Daghighi, Y.; Jarrahbashi, D. A moving-mesh finite-volume method to solve free-surface seepage problem in arbitrary geometries. Int. J. Numer. Anal. Methods Geomech. 2007, 31, 1609–1629. [Google Scholar] [CrossRef]
  6. Bresciani, E.; Davy, P.; Dreuzy, J.R. A finite volume approach with local adaptation scheme for the simulation of free surface flow in porous media. Int. J. Numer. Anal. Methods Geomech. 2011, 36, 1574–1591. [Google Scholar] [CrossRef]
  7. Sauerland, H.; Fries, T.P. The extended finite element method for two-phase and free-surface flows: A systematic study. J. Comput. Phys. 2011, 230, 3369–3390. [Google Scholar] [CrossRef]
  8. Bazyar, M.H.; Talebi, A. Locating the free surface flow in porous media using the scaled boundary finite-element method. Int. J. Chem. Eng. Appl. 2014, 5, 155–160. [Google Scholar] [CrossRef]
  9. Bathe, K.J.; Khoshgoftaar, M.R. Finite element free surface seepage analysis without mesh iteration. Int. J. Numer. Anal. Methods Geomech. 1979, 3, 13–22. [Google Scholar] [CrossRef] [Green Version]
  10. Neuman, S.P.; Witherspoon, P.A. Finite element method of analyzing steady seepage with a free surface. Water Resour. Res. 1970, 6, 889–897. [Google Scholar] [CrossRef]
  11. Chen, J.T.; Hsiao, C.C.; Chen, K.H. Study of free surface seepage problems using hypersingular equations. Commun. Numer. Methods Eng. 2007, 23, 755–769. [Google Scholar] [CrossRef]
  12. Rafiezadeh, K.; Ataie-Ashtiani, B. Transient free-surface seepage in three-dimensional general anisotropic media by BEM. Eng. Anal. Bound. Elem. 2014, 46, 51–66. [Google Scholar] [CrossRef]
  13. Leontiev, A.; Huacasi, W. Mathematical programming approach for unconfined seepage flow problem. Eng. Anal. Bound. Elem. 2016, 25, 49–56. [Google Scholar] [CrossRef]
  14. Xiao, J.E.; Ku, C.Y.; Liu, C.Y.; Fan, C.M.; Yeih, W. On solving free surface problems in layered soil using the method of fundamental solutions. Eng. Anal. Bound. Elem. 2017, 83, 96–106. [Google Scholar] [CrossRef]
  15. Pepper, D.W.; Waters, J. A local meshless method for approximating 3D wind fields. J. Appl. Meteorol. Clim. 2016, 55, 163–172. [Google Scholar] [CrossRef]
  16. Li, Z.C.; Lu, T.T.; Hu, H.Y.; Cheng, A.H.D. Trefftz and Collocation Methods; WIT Press: Southampton, UK; Boston, MA, USA, 2008. [Google Scholar]
  17. Liu, C.S. A modified Trefftz method for two-dimensional Laplace equation considering the domain’s characteristic length. CMES-Comput. Model. Eng. 2007, 21, 53–65. [Google Scholar]
  18. Liu, C.S. A highly accurate collocation Trefftz method for solving the Laplace equation in the doubly connected domains. Numer. Meth. Part. Differ. Equ. 2008, 24, 179–192. [Google Scholar] [CrossRef]
  19. Ku, C.Y.; Xiao, J.E.; Liu, C.Y.; Yeih, W. On the accuracy of the collocation Trefftz method for solving two- and three-dimensional heat equations. Numer. Heat Transf. B-Fundam. 2016, 69, 334–350. [Google Scholar] [CrossRef]
  20. Ku, C.Y. On solving three-dimensional Laplacian problems in a multiply connected domain using the multiple scale Trefftz method. CMES-Comput. Model. Eng. 2014, 98, 509–541. [Google Scholar]
  21. Huang, H.T.; Li, Z.C. Global superconvergence of Adini’s elements coupled with the Trefftz method for singular problems. Eng. Anal. Bound. Elem. 2003, 27, 227–240. [Google Scholar] [CrossRef]
  22. Lee, M.G.; Young, L.J.; Li, Z.C.; Chun, P.C. Combined Trefftz methods of particular and fundamental solutions for corner and crack singularity of linear elastostatics. Eng. Anal. Bound. Elem. 2010, 34, 632–654. [Google Scholar] [CrossRef]
  23. Kita, E.; Kamiya, N. Trefftz method: An overview. Adv. Eng. Softw. 1995, 24, 3–12. [Google Scholar] [CrossRef]
  24. Li, Z.C.; Lee, M.G.; Chiang, J.Y.; Liu, Y.P. The Trefftz method using fundamental solutions for biharmonic equations. J. Comput. Appl. Math. 2011, 235, 4350–4367. [Google Scholar] [CrossRef] [Green Version]
  25. Fu, Z.J.; Qin, Q.H.; Chen, W. Hybrid-Trefftz finite element method for heat conduction in nonlinear functionally graded materials. Eng. Comput. 2011, 28, 578–599. [Google Scholar] [CrossRef]
  26. Borkowski, M.; Hao, F.; Wang, Y.; Chen, C.S. The scalability of the matrices in direct Trefftz method in 2D Laplace problem. Eng. Anal. Bound. Elem. 2016, 64, 12–18. [Google Scholar] [CrossRef]
  27. Li, Z.C.; Lee, M.G.; Chiang, J.Y. Collocation Trefftz methods for the Stokes equations with singularity. Numer. Methods Part. Differ. Equ. 2012, 29, 361–395. [Google Scholar] [CrossRef]
  28. Li, Z.C.; Huang, H.T.; Huang, J.; Ling, L. Stability analysis for the penalty plus hybrid and the direct Trefftz methods for singularity problems. Eng. Anal. Bound. Elem. 2007, 31, 163–175. [Google Scholar] [CrossRef]
  29. Kong, W.; Wu, X. The Laplace transform and polynomial Trefftz method for a class of time dependent PDEs. Appl. Math. Model. 2009, 33, 2226–2233. [Google Scholar] [CrossRef]
  30. Ku, C.Y.; Liu, C.Y.; Xiao, J.E.; Yeih, W. Transient modeling of flow in unsaturated soils using a novel collocation meshless method. Water 2017, 9, 954. [Google Scholar] [CrossRef]
  31. Ku, C.Y.; Kuo, C.L.; Fan, C.M.; Liu, C.S.; Guan, P.C. Numerical solution of three-dimensional Laplacian problems using the multiple scale Trefftz method. Eng. Anal. Bound. Elem. 2015, 50, 157–168. [Google Scholar] [CrossRef]
  32. Lv, H.; Hao, F.; Wang, Y.; Chen, C.S. The MFS versus the Trefftz method for the Laplace equation in 3D. Eng. Anal. Bound. Elem. 2017, 83, 133–140. [Google Scholar] [CrossRef]
  33. Fan, C.M.; Liu, Y.C.; Chan, H.F.; Hsiao, S.S. Solutions of boundary detection problem for modified Helmholtz equation by Trefftz method. Inverse Probl. Sci. Eng. 2014, 22, 267–281. [Google Scholar] [CrossRef]
  34. Fan, C.M.; Chan, H.F.; Kuo, C.L.; Yeih, W. Numerical solutions of boundary detection problems using modified collocation Trefftz method and exponentially convergent scalar homotopy algoithm. Eng. Anal. Bound. Elem. 2012, 36, 2–8. [Google Scholar] [CrossRef]
  35. Yeih, W.; Liu, C.S.; Kuo, C.L.; Atluri, S.N. 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-Comput. Mater. Contin. 2010, 17, 275–302. [Google Scholar]
  36. Huang, F.K.; Chuang, M.H.; Wang, G.S.; Yeh, H.D. Tide-induced groundwater level fluctuation in a U-shaped coastal aquifer. J. Hydrol. 2015, 530, 291–305. [Google Scholar] [CrossRef]
  37. Xiao, J.E.; Ku, C.Y.; Liu, C.Y.; Yeih, W. A novel boundary-type meshless method for modeling geofluid flow in heterogeneous geological media. Geofluids 2018, 2018, 9804291. [Google Scholar] [CrossRef]
  38. Young, D.L.; Fan, C.M.; Tsai, C.C.; Chen, C.W. The method of fundamental solutions and domain decomposition method for degenerate seepage flow net problems. J. Chin. Inst. Eng. 2006, 29, 63–73. [Google Scholar] [CrossRef]
  39. Mehl, S. Use of Picard and Newton iteration for solving nonlinear groundwater flow equations. Groundwater 2006, 44, 583–594. [Google Scholar] [CrossRef] [PubMed]
  40. SoilVision Systems Ltd. SVFLUX Verification Manual; SoilVision Systems Ltd.: Saskatoon, SK, Canada, 2018. [Google Scholar]
  41. AQUAVEO. SEEP2D Earth Dam; AQUAVEO: Provo, UT, USA, 2016. [Google Scholar]
Figure 1. Nonlinear moving surface through a rectangular dam. (a) The cross section and boundary conditions and (b) collocation points on the boundary.
Figure 1. Nonlinear moving surface through a rectangular dam. (a) The cross section and boundary conditions and (b) collocation points on the boundary.
Water 11 00835 g001
Figure 2. Illustration of the collocation scheme in the CTM.
Figure 2. Illustration of the collocation scheme in the CTM.
Water 11 00835 g002
Figure 3. The cross section and the collocation scheme of the two layered soil for the analysis.
Figure 3. The cross section and the collocation scheme of the two layered soil for the analysis.
Water 11 00835 g003
Figure 4. Comparison of the domain discretization. (a) The CTM and (b) the finite element method.
Figure 4. Comparison of the domain discretization. (a) The CTM and (b) the finite element method.
Water 11 00835 g004
Figure 5. Comparison of computed results of the CTM with those from SVFLUX. (a) Computed results of the CTM and (b) SVFLUX results.
Figure 5. Comparison of computed results of the CTM with those from SVFLUX. (a) Computed results of the CTM and (b) SVFLUX results.
Water 11 00835 g005
Figure 6. Comparison of moving boundary results with other numerical methods.
Figure 6. Comparison of moving boundary results with other numerical methods.
Water 11 00835 g006
Figure 7. The convergence of the iteration step for finding the free surface.
Figure 7. The convergence of the iteration step for finding the free surface.
Water 11 00835 g007
Figure 8. The cross section of the homogeneous earth dam for the analysis.
Figure 8. The cross section of the homogeneous earth dam for the analysis.
Water 11 00835 g008
Figure 9. Free surface flow through a homogeneous dam. (a) Initial guess and the computed free surface and (b) comparison of the results with those from other methods.
Figure 9. Free surface flow through a homogeneous dam. (a) Initial guess and the computed free surface and (b) comparison of the results with those from other methods.
Water 11 00835 g009
Figure 10. The computed head and the analytical solution in two-layered soils.
Figure 10. The computed head and the analytical solution in two-layered soils.
Water 11 00835 g010
Figure 11. Absolute error of the example in two-layered soils.
Figure 11. Absolute error of the example in two-layered soils.
Water 11 00835 g011
Figure 12. The cross section and soil layer configuration of the zoned–earth dam for the analysis.
Figure 12. The cross section and soil layer configuration of the zoned–earth dam for the analysis.
Water 11 00835 g012
Figure 13. The collocation scheme of the zoned-earth dam for the DDM.
Figure 13. The collocation scheme of the zoned-earth dam for the DDM.
Water 11 00835 g013
Figure 14. Comparison of the computed free surface with SEEP2D.
Figure 14. Comparison of the computed free surface with SEEP2D.
Water 11 00835 g014
Table 1. Comparison of computed result of the separation point with those from references.
Table 1. Comparison of computed result of the separation point with those from references.
ReferenceHeight (m)
This study12.85
Aitchison [3]12.79
Chen, Hsiao, Chiu and Lee [11]12.68
Xiao, Ku, Liu, Fan and Yeih [14]12.84

Share and Cite

MDPI and ACS Style

Ku, C.-Y.; Xiao, J.-E.; Liu, C.-Y. On Solving Nonlinear Moving Boundary Problems with Heterogeneity Using the Collocation Meshless Method. Water 2019, 11, 835. https://doi.org/10.3390/w11040835

AMA Style

Ku C-Y, Xiao J-E, Liu C-Y. On Solving Nonlinear Moving Boundary Problems with Heterogeneity Using the Collocation Meshless Method. Water. 2019; 11(4):835. https://doi.org/10.3390/w11040835

Chicago/Turabian Style

Ku, Cheng-Yu, Jing-En Xiao, and Chih-Yu Liu. 2019. "On Solving Nonlinear Moving Boundary Problems with Heterogeneity Using the Collocation Meshless Method" Water 11, no. 4: 835. https://doi.org/10.3390/w11040835

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop