Modeling Transient Flows in Heterogeneous Layered Porous Media Using the Space–Time Trefftz Method

Abstract

In this study, we developed a novel boundary-type meshless approach for dealing with two-dimensional transient flows in heterogeneous layered porous media. The novelty of the proposed method is that we derived the Trefftz space–time basis function for the two-dimensional diffusion equation in layered porous media in the space–time domain. The continuity conditions at the interface of the subdomains were satisfied in terms of the domain decomposition method. Numerical solutions were approximated based on the superposition principle utilizing the space–time basis functions of the governing equation. Using the space–time collocation scheme, the numerical solutions of the problem were solved with boundary and initial data assigned on the space–time boundaries, which combined spatial and temporal discretizations in the space–time manifold. Accordingly, the transient flows through the heterogeneous layered porous media in the space–time domain could be solved without using a time-marching scheme. Numerical examples and a convergence analysis were carried out to validate the accuracy and the stability of the method. The results illustrate that an excellent agreement with the analytical solution was obtained. Additionally, the proposed method was relatively simple because we only needed to deal with the boundary data, even for the problems in the heterogeneous layered porous media. Finally, when compared with the conventional time-marching scheme, highly accurate solutions were obtained and the error accumulation from the time-marching scheme was avoided.

1. Introduction

The transport phenomena of flows through porous media have long been of research interest in numerous branches of engineering and science [1,2,3,4,5]. As heterogeneous porous media are more common than homogeneous soil in nature, flows in heterogeneous layered porous media have been widely explored [6,7,8,9], particularly utilizing mathematical modeling, numerical simulations, and experimental investigations [10,11,12].

A variety of mesh-generated techniques, such as the boundary element method [13], finite element method [14,15], finite volume method [16,17], and finite difference method [18,19], have been proposed to analyze heterogeneous layered porous media. In contrast with conventional mesh-generated techniques, meshfree approaches have attracted much attention in recent years because of their characteristics of simplicity, meshfree, and the capability to deal with engineering problems with complex geometry [20,21,22,23,24,25,26]. Of the wide variety of meshfree approaches, the Trefftz method is one of the widespread boundary-type meshless methods for dealing with steady-state Laplace-type problems, where computed results are approximated as a series of basis functions, completely satisfying the governing Laplace-type equations [27,28,29].

A comprehensive comparison of the collocation Trefftz method was proposed by Li et al. [30]. It has been concluded that the collocation Trefftz method is the simplest algorithm and provides the most accurate approximations with an optimal numerical stability. As the system of linear equations obtained from the collocation Trefftz method is ill-posed, the application of the collocation Trefftz method is less widespread [31]. Prior studies have demonstrated that applications of the collocation Trefftz method may be limited to linear and stationary problems. Recently, a study on modeling the subsurface flow problems with transient moving boundaries utilizing the Trefftz method was developed [12]. Moreover, the space–time collocation Trefftz method was developed to deal with the transient groundwater inverse problems [24]. However, the engineering applications of the Trefftz method with complete Trefftz functions for dealing with transient fluid flow through heterogeneous porous media have been less studied, which is what initiated this research.

The pore network model and the self-organized percolation (SOP) model may perform well compared with computational fluid dynamics modeling and the lattice Boltzmann method [32,33]. In contrast with other numerical methods, we developed a meshless method using the Trefftz space–time basis function. The proposed method is rooted in the Trefftz method, which is more like the semi-analytical method. Consequently, the proposed method is advantageous for dealing with problems with complex geometry because of its characteristics of simplicity and being meshfree. The proposed method may also be used to integrate regional-scale hydrological models [34].

In this study, a novel meshless approach for modeling transient flows in heterogeneous layered porous media is developed. In contrast to the conventional meshless method, the proposed approach is categorized into the boundary-type meshless method, such that the collocation points are placed only on the boundaries of the domain. Furthermore, we developed a boundary-type meshless method combining the conventional Trefftz method with the space–time collocation scheme. The proposed method is advantageous for dealing with engineering problems with a complex geometry because of its characteristics of simplicity and being meshfree. Based on the superposition principle, numerical approximations were approximated using the space–time basis functions of the governing equation. The numerical approximations of the problem were solved with boundary and initial data assigned on the space–time boundaries, which combined spatial and temporal discretizations in the space–time domain. Accordingly, the transient flows through the heterogeneous layered porous media in the space–time domain are solved without adopting the time-marching technique. The remainder of this research is organized as follows: the methodology is introduced in Section 2. In Section 3, a convergence analysis is carried out to investigate the accuracy and the stability of the proposed method. Several numerical examples are presented in Section 4. Finally, specific discussion of this research is provided in Section 5.

2. The Methodology

2.1. Governomg Equations

The two-dimensional diffusion equation used to describe transient flows is expressed as follows:

$$\frac{{\partial }^{2}h(r,\theta ,t)}{\partial r^2}+\frac{\partial h(r,\theta ,t)}{r\partial r}+\frac{{\partial }^{2}h(r,\theta ,t)}{r^2\partial {\theta }^2}=\frac{S_S}{k}\frac{\partial h(r,\theta ,t)}{\partial t}, (r,\theta ,t)\in {\Omega }_t,$$

(1)

where h is the total head, r is the radius, t is the time, $\theta$ is the polar angle, S_s is the volumetric specific storage, k is the hydraulic conductivity, and ${\Omega }_t$ is the space–time domain. For modeling transient flows in a porous media, the initial and boundary conditions are required, as follows:

$$h=H_0(r, \theta , t=0), (r,\theta ,t)\in {\Omega }_t,$$

(2)

$$h=H_D(r, \theta , t), (r,\theta ,t)\in \partial {\Omega }_t,$$

(3)

$$h_n=\frac{\partial h}{\partial n}=H_N(r, \theta , t), (r,\theta ,t)\in \partial {\Omega }_t,$$

(4)

where $H_0$ is the initial total head, $H_D$ is the Dirichlet boundary data, and $H_N$ is the Neumann boundary data.

2.2. The Space–Time Trefftz Method

The space–time Trefftz method utilizing Trefftz functions originates from the Trefftz method [12,21,24]. For a simply connected domain, one usually locates the source point in the space–time domain; the number of source points is only one for the proposed space–time Trefftz method. The approximation of Equation (1) can be expressed as follows:

$$h(r, \theta , t)\approx {\displaystyle \sum _{p=1}^m{\displaystyle \sum _{q=1}^s{A}_{pq}}}{T}_{pq}(r, \theta , t),$$

(5)

where $A_{pq}$ is a unknown coefficient, $T_{pq}(r, \theta , t)$ is the basis functions [24], and m and s denote the order of the basis functions. Considering the positive basis functions, we have the following:

$$T_{pq}(r, \theta , t)=\left\{r^q\cos (q\theta ),r^q\sin (q\theta ),e^{-p^{2}Dt}{J}_0(pr),e^{-p^{2}Dt}{J}_q(pr)\cos (q\theta ),e^{-p^{2}Dt}{J}_q(pr)\sin (q\theta )\right\},$$

(6)

where $D$ is the hydraulic diffusivity expressed as $D=k/S_S$, $J_0$ is the Bessel function of the first kind, and $J_q$ is the Bessel function of the first kind of the q order. Considering transient flows in heterogeneous layered porous media, we may write the following:

$$p=g\sqrt{\frac{D_1}{D_i}},$$

(7)

where g ranges from 1 to m; $D_1$ is the hydraulic diffusivity of first subdomain; and i denotes the ith subdomain, which means that the proposed method is applicable to transient flows in porous media including many layers. Inserting Equation (7) into Equation (6) obtains the following:

$$T_{pqi}(r, \theta , t)=\left\{r^q\cos (q\theta ),r^q\sin (q\theta ),e^{-g^2{D}_{1}t}{J}_0(g\sqrt{\frac{D_1}{D_i}}r),e^{-g^2{D}_{1}t}{J}_q(g\sqrt{\frac{D_1}{D_i}}r)\cos (q\theta ),e^{-g^2{D}_{1}t}{J}_q(g\sqrt{\frac{D_1}{D_i}}r)\sin (q\theta )\right\}.$$

(8)

Applying the Dirichlet boundary condition obtains the following:

$$\begin{array}{l}{h}_i(r, \theta , t)\approx A_1+{\displaystyle \sum _{q=1}^s{A}_2{}_q{\widehat{r}}^q\cos (q\theta )+}{\displaystyle \sum _{q=1}^s{A}_3{}_q}{\widehat{r}}^q\sin (q\theta )+{\displaystyle \sum _{g=1}^m{A}_4{}_g{e}^{-g^2\widehat{t}}{J}_0(g\sqrt{\frac{D_1}{D_i}}\widehat{r})}+\\ {\displaystyle \sum _{g=1}^m{\displaystyle \sum _{q=1}^s{J}_q(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})[A_5{}_{gq}\cos (q\theta )+A_6{}_{gq}\sin (q\theta )]e^{-g^2\widehat{t}}}},\end{array}$$

(9)

where $A_1,A_{2q},\cdots ,A_{6gq}$ denote arbitrary constants that have to be evaluated, $\widehat{r}$ is the dimensionless radius described as $\widehat{r}=r/R_0$, $R_0$ is the characteristic length, and $\widehat{t}$ is the dimensionless time written as $\widehat{t}=t{D}_1/R_0^2$. Applying the Neumann boundary condition obtains the following:

$$\frac{\partial h_i(r, \theta , t)}{\partial n}\approx {\displaystyle \sum _{q=1}^s{A}_{2q}{T_2}^{\prime }}+{\displaystyle \sum _{q=1}^s{A}_{3q}{T_3}^{\prime }}+{\displaystyle \sum _{g=1}^m{A}_{4g}{T_4}^{\prime }}+{\displaystyle \sum _{g=1}^m{\displaystyle \sum _{q=1}^s{A}_{5gq}{T_5}^{\prime }}}+{\displaystyle \sum _{g=1}^m{\displaystyle \sum _{q=1}^s{A}_{6gq}{T_6}^{\prime }}}.$$

(10)

where $n$ is the outward normal direction. The notation in Equation (10) is listed in

Table 1. The notation in Equation (10).

${T_2}^{\prime }$	$\frac{q}{R_0}\left[\begin{array}{l}\left({\widehat{r}}^{q-1}\cos (q\theta )\cos \theta +{\widehat{r}}^{q-1}\sin (q\theta )\sin \theta \right)n_x+\\ \left({\widehat{r}}^{q-1}\cos (q\theta )\sin \theta -{\widehat{r}}^{q-1}\sin (q\theta )\cos \theta \right)n_y\end{array}\right]$
${T_3}^{\prime }$	$\frac{q}{R_0}\left[\begin{array}{l}\left({\widehat{r}}^{q-1}\sin (q\theta )\cos \theta -{\widehat{r}}^{q-1}\cos (q\theta )\sin \theta \right)n_x+\\ \left({\widehat{r}}^{q-1}\sin (q\theta )\sin \theta +{\widehat{r}}^{q-1}\cos (q\theta )\cos \theta \right)n_y\end{array}\right]$
${T_4}^{\prime }$	$-\frac{g}{R_0}\sqrt{\frac{D_1}{D_i}}\left(J_1(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\cos \theta n_x+J_1(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\sin \theta n_y\right)e^{-g^2\widehat{t}}$
${T_5}^{\prime }$	$e^{-g^2\widehat{t}}\left[\begin{array}{l}\left(\frac{g}{2R_0}\sqrt{\frac{D_1}{D_i}}(J_{q-1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})-J_{q+1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r}))\cos (q\theta )\cos \theta +q{J}_{\mathrm{q}}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\sin (q\theta )\frac{\sin \theta }{r}\right)n_x+\\ \left(\frac{g}{2R_0}\sqrt{\frac{D_1}{D_i}}(J_{q-1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})-J_{q+1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r}))\cos (q\theta )\sin \theta -q{J}_{\mathrm{q}}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\sin (q\theta )\frac{\cos \theta }{r}\right)n_y\end{array}\right]$
${T_6}^{\prime }$	$e^{-g^2\widehat{t}}\left[\begin{array}{l}\left(\frac{g}{2R_0}\sqrt{\frac{D_1}{D_i}}(J_{q-1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})-J_{q+1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r}))\sin (q\theta )\cos \theta -q{J}_q(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\cos (q\theta )\frac{\sin \theta }{r}\right)n_x+\\ \left(\frac{g}{2R_0}\sqrt{\frac{D_1}{D_i}}(J_{q-1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})-J_{q+1}(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r}))\sin (q\theta )\sin \theta +q{J}_q(\mathrm{g}\sqrt{\frac{D_1}{D_i}}\widehat{r})\cos (q\theta )\frac{\cos \theta }{r}\right)n_y\end{array}\right]$

.

2.3. The Continuity Conditions of the Interface

To model transient flows in heterogeneous layered porous media, the domain decomposition method (DDM) was utilized [35]. The continuity of the head and the flux conservation at the interface between two layered porous media were considered.

An irregular domain, which was split into two subdomains, ${\Omega }_1$ and ${\Omega }_2$, as displayed in Figure 1a, was considered. To model transient flows in heterogeneous layered porous media, the irregular domain was divided into several sub-boundaries. The sub-boundaries at the interface must satisfy the following boundary conditions:

$$h|{}_{{\Gamma }_{12}}=h|{}_{{\Gamma }_{22}},k_1\frac{\partial h}{\partial n}|{}_{{\Gamma }_{12}}=k_2\frac{\partial h}{\partial n}|{}_{{\Gamma }_{22}},$$

(11)

where ${\Gamma }_{12}$ and ${\Gamma }_{22}$ are sub-boundaries at the interface shown in Figure 1a. Applying the boundary data on boundary points, the following simultaneous linear equations are obtained as follows:

$$A_W{c}_W=B_W,$$

(12)

$${\mathrm{A}}_W=\left[\begin{array}{cc}{\mathrm{A}}_{{\mathsf{\Omega }}_1}& 0_{{\mathsf{\Omega }}_2}\\ 0_{{\mathsf{\Omega }}_1}& {\mathrm{A}}_{{\mathsf{\Omega }}_2}\\ {\mathrm{A}}_{\mathrm{I}}|{}_{{\mathsf{\Gamma }}_{12}}& {\mathrm{A}}_{\mathrm{I}}|{}_{{\mathsf{\Gamma }}_{22}}\end{array}\right],{ \mathrm{c}}_W=\left[\begin{array}{c}{\mathrm{c}}_{{\Omega }_1}\\ {\mathrm{c}}_{\Omega }{}_{{}_2}\end{array}\right],{ \mathrm{B}}_W=\left[\begin{array}{c}{\mathrm{B}}_{{\mathsf{\Omega }}_1}\\ {\mathrm{B}}_{{\mathsf{\Omega }}_2}\\ {\mathrm{B}}_{\mathrm{I}}\end{array}\right],$$

(13)

$${\mathrm{A}}_{{\Omega }_i}=\left[\begin{array}{ccccc}1& {\widehat{r}}_1\cos ({\theta }_1)& {\widehat{r}}_1{}^2\cos (2{\theta }_1)& \cdots & J_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_1)\sin (q{\theta }_1)e^{-g^2{\widehat{t}}_1}\\ 1& {\widehat{r}}_2\cos ({\theta }_2)& {\widehat{r}}_2{}^2\cos (2{\theta }_2)& \cdots & J_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_2)\sin (q{\theta }_2)e^{-g^2{\widehat{t}}_2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\widehat{r}}_{\nu }\cos ({\theta }_{\nu })& {\widehat{r}}_{\nu }{}^2\cos (2{\theta }_{\nu })& \cdots & J_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_{\nu })\sin (q{\theta }_{\nu })e^{-g^2{\widehat{t}}_{\nu }}\end{array}\right],$$

(14)

$${\mathrm{A}}_I=\left[\begin{array}{ccccc}{(-1)}^i& {(-1)}^i{\left({\widehat{r}}_1\right)}^q\cos (q{\theta }_1)& {(-1)}^i{\left({\widehat{r}}_1\right)}^q\sin (q{\theta }_1)& \cdots & {(-1)}^i{J}_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_1)\sin (q{\theta }_1)e^{-g^2{\widehat{t}}_1}\\ {(-1)}^i& {(-1)}^i{\left({\widehat{r}}_2\right)}^q\cos (q{\theta }_2)& {(-1)}^i{\left({\widehat{r}}_2\right)}^q\sin (q{\theta }_2)& \cdots & {(-1)}^i{J}_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_2)\sin (q{\theta }_2)e^{-g^2{\widehat{t}}_2}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {(-1)}^i& {(-1)}^i{\left({\widehat{r}}_v\right)}^q\cos (q{\theta }_v)& {(-1)}^i{\left({\widehat{r}}_v\right)}^q\sin (q{\theta }_v)& \cdots & {(-1)}^i{J}_q(\sqrt{\frac{D_1}{D_i}}g{\widehat{r}}_v)\sin (q{\theta }_v)e^{-g^2{\widehat{t}}_v}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& k_i{{T_2}^{\prime }}_{(v)}& k_i{{T_3}^{\prime }}_{(v)}& \cdots & k_i{{T_6}^{\prime }}_{(v)}\end{array}\right],$$

(15)

where $A_{{\mathsf{\Omega }}_1}$ and $A_{{\mathsf{\Omega }}_2}$ are the ${\mathrm{A}}_{\mathrm{w}}$ matrix depicted in Equation (13) for ${\Omega }_1$ and ${\Omega }_2$, respectively. ${\mathrm{A}}_{\mathrm{I}}|{}_{{\mathsf{\Gamma }}_{12}}$ of the boundary ${\Gamma }_{12}$ and $A_{\mathrm{I}}|{}_{{\mathsf{\Gamma }}_{22}}$ of the boundary ${\Gamma }_{22}$ are the ${\mathrm{A}}_{\mathrm{I}}$ matrices at the interface. $0_{{\mathsf{\Omega }}_1}$ and $0_{{\mathsf{\Omega }}_2}$ are the zero matrices. $c_{{\mathsf{\Omega }}_1}$ and $c_{{\mathsf{\Omega }}_2}$ are vectors of the undetermined coefficients of ${\Omega }_1$ and ${\Omega }_2$, respectively. $B_{{\mathsf{\Omega }}_1}$ and $B_{{\mathsf{\Omega }}_2}$ are vectors of the boundary values of ${\Omega }_1$ and ${\Omega }_2$, respectively. $\nu$ and $\upsilon$ are the numbers of the subdomain and interface points, $B_I={\left[\begin{array}{cc}{0}_{{\Gamma }_{12}}& 0_{{\Gamma }_{22}}\end{array}\right]}^{\mathrm{T}}$, and $0_{{\Gamma }_{12}}$ and $0_{{\Gamma }_{22}}$ are zero vectors with a size of $2\upsilon \times 1$ for the Dirichlet and Neumann boundary conditions at the interface, satisfying flux conservation and continuity of the head for the heterogeneous layered porous media problem.

By solving Equation (12), we can obtain two sets of unknown coefficients, $c_{{\mathsf{\Omega }}_1}$ and $c_{{\mathsf{\Omega }}_2}$, for subdomains ${\Omega }_1$ and ${\Omega }_2$, respectively. The computed head at the inner points can then be computed by Equation (5), utilizing $A_{{\mathsf{\Omega }}_1}$ and $B_{{\mathsf{\Omega }}_1}$ for ${\Omega }_1$, and $A_{{\mathsf{\Omega }}_2}$ and $B_{{\mathsf{\Omega }}_2}$ for ${\Omega }_2$. In this study, the commercial program MATLAB backslash operator was utilized to solve Equation (12). The schematic diagram for modeling the transient flows in the heterogeneous layered porous media is displayed in Figure 2.

The continuity conditions at the interface of the subdomains are satisfied in terms of the DDM, such that the proposed method is capable of dealing with transient flows in heterogeneous layered porous media. Numerical solutions are approximated based on the superposition theorem using the space–time basis functions of the governing equation. Using the space–time collocation scheme, the numerical solutions of the problem are solved with boundary and initial data assigned on the space–time boundaries, as shown in Figure 1b. The spatial and temporal discretizations are then combined in the space–time domain. Accordingly, the transient flows through the heterogeneous layered porous media in the space–time domain can be solved without using the time-marching scheme.

3. Convergence Analysis of Transient Flows in Heterogeneous Layered Porous Media

In this section, convergence analysis of the transient flows in heterogeneous layered porous media is carried out. To illustrate the accuracy of the approximations, the error measures are defined as

$$\mathrm{RMSE}=\sqrt{\frac{1}{N_I}{\displaystyle \sum _{i=1}^{N_I}{[h(r_i,{\theta }_i,t_i)-\widehat{h}(r_i,{\theta }_i,t_i\left)\right]}^2}},$$

(16)

$$\mathrm{MAE}=\max \{\left|\widehat{h}(r_i,{\theta }_i,t_i)\right|-\left|h(r_i,{\theta }_i,t_i)\right|\}, 1\le i\le N_I,$$

(17)

where $h(r_i,{\theta }_i,t_i)$ is the exact solution, $\widehat{h}(r_i,{\theta }_i,t_i)$ is the approximate solution, RMSE denotes the root mean square error, MAE denotes the maximum absolute error, and $N_I$ is the number of scattered inner points. Here, the scattered inner points were uniformly distributed in the space–time domain.

Transient flows in two-layered porous media satisfy the diffusion equation depicted in Equation (1) in two dimensions over an irregular domain. The geometry for this problem is

$$\partial {\Omega }_t=\left\{(r, \theta , t)|\begin{array}{l}r=2(1+0.1\tanh (10\sin (10\theta +\frac{\pi }{3})),\\ 0\le \theta \le 2\pi ,0\le t\le 1\end{array}\right\}.$$

(18)

The initial data for subdomains ${\Omega }_1$ and ${\Omega }_2$ are imposed as follows:

$$h_1(r, \theta , t=0)=\sin (r\sin \theta )r\cos \theta ,(r, \theta , t=0)\in {\Omega }_1$$

(19)

$$h_2(r, \theta , t=0)=\sin (\sqrt{\frac{D_1}{D_2}}r\sin \theta )r\cos \theta ,(r, \theta , t=0)\in {\Omega }_2,$$

(20)

where $h_1(r, \theta , t=0)$ and $h_2(r, \theta , t=0)$ are the initial data for subdomains ${\Omega }_1$ and ${\Omega }_2$, respectively. $D_1$ and $D_2$ are the hydraulic diffusivity in layer 1 and 2, respectively. The boundary conditions for subdomain ${\Omega }_1$ and ${\Omega }_2$ are prescribed based on the following functions:

$$h_1(r, \theta , t)=e^{-2D_{1}t}\sin (r\sin \theta )r\cos \theta ,(r,\theta ,t)\in \partial {\Omega }_1,$$

(21)

$$h_2(r, \theta , t)=e^{-2D_{1}t}\sin (\sqrt{\frac{D_1}{D_2}}r\sin \theta )r\cos \theta ,(r,\theta ,t)\in \partial {\Omega }_2.$$

(22)

The exact solution is the function as follows:

$$h=e^{-2D_{1}t}\sin (\sqrt{\frac{D_1}{D_i}}r\sin \theta )r\cos \theta ,$$

(23)

where $D_i$ is the hydraulic diffusivity of the ith subdomain. The hydraulic diffusivity and hydraulic conductivity in two subdomains were set to be $D_1={1 \mathrm{m}}^2{\mathrm{h}}^{-1}$ and $D_2=0{.25 \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={2 \mathrm{m}}^1{\mathrm{h}}^{-1}$ and $k_2={1 \mathrm{m}}^1{\mathrm{h}}^{-1}$. The final elapsed time is 1 h. The boundary points are collocated on the boundary, as illustrated in Figure 1b. For simulating transient flows in two-layered porous media, DDM is utilized [35]. The domain boundary can be divided into several subdomains.

The Dirichlet boundary conditions are based on Equations (21) and (22). The sub-boundaries at the interface, including ${\Gamma }_{12}$ and ${\Gamma }_{22}$, must satisfy the continuity of flows and the flux conservation based on Equation (11).

Figure 3a shows the error versus the order of the basis functions. It appears that approximations with a high accuracy are acquired then the order is superior to 14. Additionally, promising results may be obtained when the order of the basis function ranges from 14 to 26. It seems that the optimal order may depend on the nature of the problems. The relationship between the absoluter error and the number of boundary points is displayed in Figure 3b. The results obtained illustrate that approximations with a high accuracy are achieved when the number of boundary points is greater than 1800. From the results of this convergence analysis, the number of the boundary points and the order of the basis functions were determined to be 2000 and 15, respectively.

4. Numerical Examples

4.1. Transient Flows in Two-Layered Porous Media

Transient flows in two-layered porous media satisfy the diffusion equation depicted in Equation (1) in two dimensions over an irregular domain. The geometry for this problem, as shown in Figure 4a, is

$$\partial {\Omega }_t=\left\{(r, \theta , t)|\begin{array}{l}r={|\sec (3\theta )|}^{\sin (6\theta )} ,\\ 0\le \theta \le 2\pi , 0\le t\le 1\end{array}\right\}.$$

(24)

The initial data for subdomains ${\Omega }_1$ and ${\Omega }_2$ are imposed as follows.

$$h_1(r,\theta ,t=0)=\sin \mathrm{h}(r\sin \theta )r\cos \theta ,(r, \theta , t=0)\in {\Omega }_1$$

(25)

$$h_2(r,\theta ,t=0)=\sin \mathrm{h}(\sqrt{\frac{D_1}{D_2}}r\sin \theta )r\cos \theta ,(r, \theta , t=0)\in {\Omega }_2$$

(26)

The boundary conditions for subdomains ${\Omega }_1$ and ${\Omega }_2$ are prescribed based on the following functions.

$$h_1(r,\theta ,t)=e^{2D_{1}t}\sin \mathrm{h}(r\sin \theta )r\cos \theta ,(r,\theta ,t)\in \partial {\Omega }_1$$

(27)

$$h_2(r,\theta ,t)=e^{2D_{1}t}\sin \mathrm{h}(\sqrt{\frac{D_1}{D_2}}r\sin \theta )r\cos \theta ,(r,\theta ,t)\in \partial {\Omega }_2$$

(28)

The exact solution is the following function:

$$h=e^{2D_{1}t}\sin \mathrm{h}(\sqrt{\frac{D_1}{D_i}}r\sin \theta )r\cos \theta .$$

(29)

In this example, the hydraulic diffusivity and hydraulic conductivity in the two subdomains were set to be $D_1=\frac{1}{9}{ \mathrm{m}}^2{\mathrm{h}}^{-1}$ and $D_2=1{ \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={2 \mathrm{m}}^1{\mathrm{h}}^{-1}$ and $k_2={6 \mathrm{m}}^1{\mathrm{h}}^{-1}$. The final elapsed time is 1 h. The order of the basis functions is 17. There are 2026 boundary points for this domain. This example involves two-dimensional space and one-dimensional time. The space–time domain is therefore transformed into a three-dimensional object domain, as shown in Figure 4b. The boundary points are placed on the boundary, as illustrated in Figure 4b. For simulating transient flows in two–layered porous media, DDM is utilized [35]. The domain boundary can then be split into two subdomains. The Dirichlet boundary conditions are prescribed based on Equations (27) and (28). Based on Equation (11), the sub-boundaries at the interface, including ${\Gamma }_{12}$ and ${\Gamma }_{22}$, must satisfy the continuity of flows and the flux conservation.

To compare the accuracy of our method with the analytical solution, the profiles of the results at different times are selected, as presented in Figure 5. Figure 5 shows that the computed results adopting our approach are completely consistent with the analytical solution.

To examine the accuracy of our method, the absolute error of the approximate solutions with the exact solution is evaluated. Figure 6 displays the absolute error at different simulation times. From Figure 6, the MAE is in the order of 10⁻¹⁰, which indicates that our approach may acquire results with a high accuracy.

To show the accuracy of the computed result from the proposed method and that from conventional mesh-generated techniques, we conducted a comparison example where transient flows in two-layered porous media satisfy the diffusion equation depicted in Equation (1) in two dimensions over a rectangular domain. The length and the width of the space domain are 4 m. The initial data for subdomains ${\Omega }_1$ and ${\Omega }_2$ are imposed as Equations (25) and (26). The boundary conditions for subdomains ${\Omega }_1$ and ${\Omega }_2$ are prescribed based on Equations (27) and (28). The exact solution for solving this example is as Equation (29).

In this example, the hydraulic diffusivity and hydraulic conductivity in the two subdomains were set to be $D_1=\frac{1}{9}{ \mathrm{m}}^2{\mathrm{h}}^{-1}$ and $D_2=1{ \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={2 \mathrm{m}}^1{\mathrm{h}}^{-1}$ and $k_2={6 \mathrm{m}}^1{\mathrm{h}}^{-1}$. The final elapsed time is 4 h. The order of the basis functions is 16. There are 2600 boundary points for this domain.

This example was solved using both the finite difference method (FDM) and the proposed method. For the FDM analysis, the spatial and temporal discretizations for the problem must be considered separately. The central difference approximation was adopted for the spatial discretization, and the implicit scheme was adopted for time discretization. For FDM, we considered each length of Δx along the x axis, each length of Δy along the y axis, and each length of Δt along the t axis to be 0.1 m, 0.1 m, and 0.001 h, respectively. We then compared the approximations of our method with those of the FDM.

Figure 7 presents the absolute error versus the simulation time. From Figure 7, the magnitude of the absolute error increases with the simulation time in the FDM. However, our approach may avoid error propagation. The results demonstrate that highly accurate numerical solutions can be obtained using our approach. Moreover, the accuracy of the error is of the order of 10⁻⁷ to 10⁻⁸.

4.2. Transient Flows in Two-Layered Porous Media

Transient flows in two-layered porous media satisfy the diffusion equation depicted in Equation (1) in two dimensions over a rectangular domain, as depicted in Figure 8a. Two cases with different combinations of composite porous media are considered. The geometry for case I, as shown in Figure 8a, is

$$\partial {\Omega }_{a1}=\left\{\left(x, y, t\right)|-5 \le x\le 0,-5 \le y\le 5, 0\le t\le 3\right\},$$

(30)

$$\partial {\Omega }_{a2}=\left\{\left(x, y, t\right)|0\le x\le 5, -5\le y\le 5,0\le t\le 3\right\}.$$

(31)

The hydraulic diffusivity and hydraulic conductivity in the two subdomains were set to be $D_1={100 \mathrm{m}}^2{\mathrm{h}}^{-1}$ and $D_2={10 \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={100 \mathrm{mh}}^{-1}$ and $k_2={10 \mathrm{mh}}^{-1}$. The initial condition is

$$h_a(x,y,t=0)=10 \mathrm{m}.$$

(32)

The boundary conditions are

$$h_a(-5,y,t)=10 \mathrm{m},$$

(33)

$$h_a(5,y,t)=4 \mathrm{m},$$

(34)

$$\frac{\partial h_a(x,-5,t)}{\partial n}=\frac{\partial h_a(x,5,t)}{\partial n}=0,$$

(35)

where $h_a$ is the total head for case I.

The geometry for case II, as shown in Figure 8b, is

$$\partial {\Omega }_{b1}=\left\{\left(x,y,t\right)|-5\le x\le -2, -5\le y\le 5,0\le t\le 3\right\},$$

(36)

$$\partial {\Omega }_{b2}=\left\{\left(x,y,t\right)|-2\le x\le 5, -5\le y\le 5,0\le t\le 3\right\}.$$

(37)

The hydraulic diffusivity and hydraulic conductivity in the two subdomains are set to be $D_1={10 \mathrm{m}}^2{\mathrm{h}}^{-1}$ and $D_2={100 \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={10 \mathrm{mh}}^{-1}$ and $k_2={100 \mathrm{mh}}^{-1}$. The initial condition is

$$h_b(x,y,t=0)=4 \mathrm{m}.$$

(38)

The boundary conditions are

$$h_b(-5,y,t)=10 \mathrm{m},$$

(39)

$$h_b(5,y,t)=4 \mathrm{m},$$

(40)

$$\frac{\partial h_b(x,-5,t)}{\partial n}=\frac{\partial h_b(x,5,t)}{\partial n}=0,$$

(41)

where $h_b$ is the total head for case II.

As this example may not have an exact solution to investigate the accuracy, we compare the approximations at the final time with the steady-state solution. The steady-state solutions for this two-layered porous media are prescribed based on the following equations:

$$h|{}_{{\Gamma }_{12}}=h|{}_{{\Gamma }_{22}}=\frac{L_2{k}_{1}h|{}_{{\Gamma }_{14}}+L_1{k}_{2}h|{}_{{\Gamma }_{24}}}{L_2{k}_1+L_1{k}_2},$$

(42)

$$h_1=\frac{h|{}_{{\Gamma }_{12}}-h|{}_{{\Gamma }_{14}}}{L_1}(x+5)+h|{}_{{\Gamma }_{14}},$$

(43)

$$h_2=\frac{h|{}_{{\Gamma }_{24}}-h|{}_{{\Gamma }_{22}}}{L_2}(x+5)+h|{}_{{\Gamma }_{22}}-(h|{}_{{\Gamma }_{24}}-h|{}_{{\Gamma }_{22}})\frac{L_1}{L_2}.$$

(44)

In case I and case II, the final time is 3 h. The order of the basis functions is 15. There exists 7200 boundary points for domain, as illustrated in Figure 9. For modeling transient flows in two-layered porous media, the DDM was utilized [35]. The domain boundary can then be divided into several subdomains. The Dirichlet boundary conditions are prescribed. Based on Equation (11), the sub-boundaries at the interface, including ${\Gamma }_{12}$ and ${\Gamma }_{22}$, must satisfy the continuity of flows and the flux conservation.

The profiles of total head at the final time were selected to express the approximations, as shown in Figure 10. The approximations at the final time from our approach were further compared with the steady-state solution. It is found that the approximations utilizing our approach may agree well with the steady-state solutions. The absolute difference at the final time is depicted in Figure 11. The absolute difference for case I and case II is in the order of 10⁻⁶ and 10⁻⁴, as given in Figure 11. Results demonstrate that our approach may obtain reasonable results for modeling transient flows in two-layered porous media.

4.3. Transient Flows in Four-Layered Porous Media

The last example is modeling of transient flows in four-layered porous media, as depicted in Figure 12. The diffusion equation depicted in Equation (1) in two dimensions over a rectangular domain is considered. The geometry for this example is defined as

$$\partial {\Omega }_1=\left\{\left(x, y,t\right)|-4.5\le x\le 0, -4.5\le y\le 0,0\le t\le 3\right\},$$

(45)

$$\partial {\Omega }_2=\left\{\left(x, y,t\right)|0\le x\le 4.5, -4.5\le y\le 0,0\le t\le 3\right\},$$

(46)

$$\partial {\Omega }_3=\left\{\left(x, y,t\right)|-4.5\le x\le 0, 0\le y\le 4.5, 0\le t\le 3\right\},$$

(47)

$$\partial {\Omega }_4=\left\{\left(x, y,t\right)|0\le x\le 4.5, 0\le y\le 4.5,0\le t\le 3\right\}.$$

(48)

The hydraulic diffusivity and hydraulic conductivity in the two subdomains are set to be $D_1={4 \mathrm{m}}^2{\mathrm{h}}^{-1}$, $D_2={3 \mathrm{m}}^2{\mathrm{h}}^{-1}$, $D_3=1{.5 \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $D_4={1 \mathrm{m}}^2{\mathrm{h}}^{-1}$, and $k_1={4 \mathrm{mh}}^{-1}$, $k_2={3 \mathrm{mh}}^{-1}$, $k_3=1.5{ \mathrm{mh}}^{-1}$, and $k_4={1 \mathrm{mh}}^{-1}$. The initial condition is

$$h(x,y,t=0)=5 \mathrm{m}.$$

(49)

The boundary condition is

$$h(x,y,t)=\sqrt{253-\frac{76x}{3}-\frac{76y}{3}} \mathrm{m}.$$

(50)

As this example may not have exact solutions to investigate the accuracy, we compare the approximations at different simulation times with the numerical solutions using the FDM. In this example, the final elapsed time is 3 h. The order of the basis functions needs to be 20. There are 6400 boundary points that exist for each domain, as illustrated in Figure 13. For simulating transient flows in four-layered porous media, DDM is utilized [35]. The domain boundary can then be divided into several subdomains. The Dirichlet boundary conditions are prescribed. At the interface, the sub-boundaries, including ${\Gamma }_{12}$, ${\Gamma }_{22}$, ${\Gamma }_{13}$, ${\Gamma }_{41}$, ${\Gamma }_{23}$, ${\Gamma }_{31}$, ${\Gamma }_{32}$, and ${\Gamma }_{42}$, must satisfy the continuity of flows and the flux conservation based on the following equations:

$$h|{}_{{\Gamma }_{12}}=h|{}_{{\Gamma }_{22}},k_1\frac{\partial h}{\partial n}|{}_{{\Gamma }_{12}}=k_2\frac{\partial h}{\partial n}|{}_{{\Gamma }_{22}},$$

(51)

$$h|{}_{{\Gamma }_{13}}=h|{}_{{\Gamma }_{41}},k_1\frac{\partial h}{\partial n}|{}_{{\Gamma }_{13}}=k_4\frac{\partial h}{\partial n}|{}_{{\Gamma }_{41}},$$

(52)

$$h|{}_{{\Gamma }_{23}}=h|{}_{{\Gamma }_{31}},k_2\frac{\partial h}{\partial n}|{}_{{\Gamma }_{23}}=k_3\frac{\partial h}{\partial n}|{}_{{\Gamma }_{31}},$$

(53)

$$h|{}_{{\Gamma }_{32}}=h|{}_{{\Gamma }_{42}},k_3\frac{\partial h}{\partial n}|{}_{{\Gamma }_{32}}=k_4\frac{\partial h}{\partial n}|{}_{{\Gamma }_{42}}.$$

(54)

The approximations at different simulation times from our approach are further compared with those of FDM. To apply FDM, discretization in the spatial and temporal domains must be considered. For the spatial discretization, the domain is divided into sections, and second derivatives of the diffusion equation for each grid point are approximated using central difference formulas. For the time discretization, the implicit scheme is adopted. For FDM, we consider each length of Δx along the x axis, each length of Δy along the y axis, and each length of Δt along the t axis to be 0.225 m, 0.225 m, and 0.01 h, respectively. Then, we compare the approximations of our method with those of FDM.

The profiles of the total head at different simulation times are selected to express the approximations, as shown in Figure 14. From Figure 14, it is found that the approximations utilizing our approach may be consistent with the FDM results. The results illustrate that our approach may obtain reasonable results for modeling transient flows in four-layered porous media.

5. Discussion

This study is rooted in the Trefftz method and demonstrates a promising numerical solution for dealing with two-dimensional transient flows in heterogeneous layered porous media. The comparison of the proposed method with previous studies depicts the following advantages.

Various mesh-based numerical approaches have been utilized for dealing with the transient and heterogeneous groundwater flow equation. In contrast to conventional mesh-generated techniques, we developed the meshless method using the Trefftz space–time basis function. The proposed method is advantageous for dealing with engineering problems with a complex geometry because of the characteristics of simplicity and being meshfree.

Contrary to the conventional meshless method, the proposed method is categorized into the boundary-type meshless method, such that the collocation points are only placed on the boundaries of the domain. Furthermore, the proposed method combines the conventional Trefftz method with the space–time collocation scheme, such that no inner points are required in the analysis and all collocation points needed to place the space–time boundaries are presented for the modeling of the two-dimensional transient flows in the heterogeneous layered porous media. As a result, both the initial and boundary data are regarded as boundary conditions on the space–time boundary, which is especially advantageous for an irregular domain shape.

Using the space–time collocation, the initial value problem for simulating transient flows in heterogeneous layered porous media is considered to be a problem of the inverse boundary value, where the time marching for the initial value problem can be avoided. Accordingly, the proposed method is advantageous for dealing with inverse problems as well.

On the other hand, the limitations of the proposed method may be raised, as we need to derive the Trefftz space–time basis function for the two-dimensional diffusion equation in layered porous media in the space–time domain, in which the Trefftz basis functions are obtained from the general solutions using the separation of variables. The solutions of the governing equation are then approximated numerically based on the superposition principle utilizing the space–time basis functions of the governing equation. Consequently, the proposed method is limited to the linear governing equation with accessible general solutions.

6. Conclusions

A novel meshless method for dealing with two-dimensional transient flows in heterogeneous layered porous media is presented in this article. The following key findings are given.

To deal with transient flows in heterogeneous layered porous media, the continuity conditions at the interface of the subdomains are satisfied in terms of the domain decomposition method. Numerical solutions are approximated based on the superposition principle adopting the space–time basis functions of the diffusion equation. Utilizing the space–time collocation scheme, numerical approximations are solved with boundary and initial data assigned on the space–time boundaries, which combine spatial and temporal discretizations in the space–time domain. The transient flows through heterogeneous layered porous media in the space–time domain may therefore be solved without using the time-marching scheme, which is the novelty of this study.

The results obtained show that an excellent agreement with the exact solution can be acquired. Additionally, the proposed approach may be relatively simple, because we only have to deal with the boundary data even with problems in heterogeneous layered porous media. Finally, compared with the conventional time-marching scheme, highly accurate solutions can be obtained and the error accumulation from the time-marching scheme can be avoided.