Geofluids

Volume 2018, Article ID 9804291, 13 pages

https://doi.org/10.1155/2018/9804291

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

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

Correspondence should be addressed to Cheng-Yu Ku; wt.ude.uotn.liame@62tskhc

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