A Spacetime Meshless Method for Modeling Subsurface Flow with a Transient Moving Boundary

Abstract

In this paper, a spacetime meshless method utilizing Trefftz functions for modeling subsurface flow problems with a transient moving boundary is proposed. The subsurface flow problem with a transient moving boundary is governed by the two-dimensional diffusion equation, where the position of the moving boundary is previously unknown. We solve the subsurface flow problems based on the Trefftz method, 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 by the superposition theorem using the basis functions, which match the data at the spacetime boundary collocation points. Because the proposed basis functions fully satisfy the diffusion equation, arbitrary nodes are collocated only on the spacetime boundaries for the discretization of the domain. The iterative scheme has to be used for solving the moving boundaries because the transient moving boundary problems exhibit nonlinear characteristics. Numerical examples, including harmonic and non-harmonic boundary conditions, are carried out to validate the method. Results illustrate that our method may acquire field solutions with high accuracy. It is also found that the method is advantageous for solving inverse problems as well. Finally, comparing with those obtained from the method of fundamental solutions, we may obtain the accurate location of the nonlinear moving boundary for transient problems using the spacetime meshless method with the iterative scheme.

1. Introduction

The free surface flow in soils can be defined as moving boundary problems because the location of one or more of the domain boundaries is unknown [1,2,3,4]. The phreatic line is located between the fluid phase and the air phase of the soil. It is sometimes regarded as the phase change problem [5,6,7]. Phase change problems are often encountered in engineering, industry, and problems such as the design of roadways in cold regions [8,9,10]. These problems are usually non-stationary as well as nonlinear due to the phase change depending on time and complicated material properties [11]. Therefore, great challenges may be raised for solving the problems using analytical solutions.

Various numerical approaches [12], such as the boundary element method [13], the interpolation finite difference method [14], the finite element method [15], the finite volume method [16], the local radial basis function collocation method [17], the method of approximate particular solutions [18], and the method of fundamental solutions (MFS) [19,20], have been utilized for dealing with moving boundary problems. The collocation method can be categorized into one of the meshless methods [21,22]. The discretization of the domain for the collocation approaches is relatively simple because the arbitrary points are assigned only on the boundary if we find the basis functions, which must satisfy the governing equation [23,24]. This idea was developed by Erich Trefftz and known as the Trefftz method [25]. The Trefftz method is widely adopted for dealing with the boundary value problem (BVP) [26,27,28,29]. This method is often found to solve the Laplace-type equations. The reason is that the derivation of the Trefftz functions for other partial differential equations may be very challenging [30]. Previous studies have demonstrated that applications of the Trefftz method may be limited to linear as well as stationary problems. Recently, a study on solving subsurface flow problems with free and moving boundaries governed by the Laplace governing equation adopting the Trefftz method has been developed [23]. However, the engineering applications of the Trefftz method with complete Trefftz functions for dealing with transient problems are still hardly found, where solving transient moving boundary problems adopting the Trefftz method rarely exist.

In this study, we propose the spacetime meshless approach using Trefftz functions for solving subsurface flow problems with a transient moving boundary. The subsurface flow problem with a transient moving boundary is governed by the two-dimensional diffusion equation, where the position of the nonlinear moving boundary is previously unknown. We solve the subsurface flow problems utilizing the Trefftz method, in which the Trefftz basis functions can be obtained from the general solutions using the separation of variables. We propose the spacetime collocation scheme such that the solutions of the governing equation are approximated numerically by the superposition theorem using the proposed basis functions, which match the data at spacetime boundary collocation points. Because the proposed basis functions fully satisfy the diffusion equation, arbitrary nodes are collocated only on the spacetime boundaries for the discretization of the domain. Because the transient moving boundary problems exhibit nonlinear characteristics, the iterative scheme has to be used for solving the moving boundaries. Numerical examples, including harmonic and non-harmonic boundary conditions, were carried out to validate the method. The derivation of the spacetime collocation scheme utilizing Trefftz functions is depicted in the following section.

2. The Governing Equation

The transient moving boundary phenomenon is governed by the two-dimensional diffusion governing equation in the dimensionless form in polar coordinates as follows,

$$\frac{{\partial }^{2}h(r,\theta ,t)}{\partial r^2}+\frac{1}{r}\frac{\partial h(r,\theta ,t)}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}h(r,\theta ,t)}{\partial {\theta }^2}=\frac{\partial h(r,\theta ,t)}{\partial t} \mathrm{in} {\Re }_t,$$

(1)

where ${\Re }_t$ is the spacetime domain of the transient moving boundary problems, h is the total head, $\theta$ denotes the polar angle, r is the dimensionless variable, which is expressed as $r=\widehat{r}/R_0$, $\widehat{r}$ is the radius, $R_0$ is the maximum dimension of the problem, which is also named as the length of characteristic, t is the dimensionless variable, which is expressed as $t=\widehat{t}T/R_0^{2}S$, S is the storage coefficient, $\widehat{t}$ is the time, and T is the transmissibility coefficient. While $\widehat{r}$ is in the range of $0<\widehat{r}<R_0$, r is in the range of $0<r<1$. The initial condition for solving Equation (1) is as follows,

$$h(r,\theta ,t=0)=h_0(r,\theta ,t=0),$$

(2)

where $h_0$ is the initial total head. To solve Equation (1), the boundary data are expressed as

$$h(r,\theta ,t)=D(r,\theta ,t) \mathrm{on} {\mathsf{\Gamma }}_D^t,$$

(3)

$$h_n(r,\theta ,t)=\frac{\partial h(r,\theta ,t)}{\partial n}=N(r,\theta ,t) \mathrm{on} {\mathsf{\Gamma }}_N^t,$$

(4)

where ${\mathsf{\Gamma }}_D^t$ is the spacetime Dirichlet boundary condition, ${\mathsf{\Gamma }}_N^t$ is the spacetime Neumann boundary condition, the subscript D is the Dirichlet boundary data, the subscript N denotes the Neumann boundary data, the subscript n is the outward normal direction, $D(r,\theta ,t)$ and $N(r,\theta ,t)$ are the Dirichlet and Neumann boundary data of the transient moving boundary problems in the spacetime domain, respectively.

3. The Spacetime Meshless Method Using Trefftz Functions

3.1. Trefftz Functions for Transient Moving Boundary Problems

The spacetime meshless method using Trefftz functions is rooted in the Trefftz method. Thus, it is necessary for the nonlinear moving boundary problems to formulate the general solutions. To formulate the transient Trefftz functions for the nonlinear moving boundary problems, the separation of variables is adopted.

$$h(r,\theta ,t)=\phi (r,\theta )\Omega (t),$$

(5)

where $\phi (r,\theta )$ and $\Omega (t)$ are functions. The total head $h(r,\theta ,t)$ is a product of two functions. For simplicity, the following notations are considered.

$${\phi }_r=\frac{d\phi (r,\theta )}{dr}, {\phi }_{rr}=\frac{d^2\phi (r,\theta )}{d{r}^2}, {\phi }_{\theta \theta }=\frac{d^2\phi (r,\theta )}{d{\theta }^2} \mathrm{and} {\Omega }^{\prime }=\frac{d\Omega (t)}{dt},$$

(6)

where the subscript r denotes the first derivative with respect to r, the subscript rr denotes the second derivative with respect to r, the subscript $\theta \theta$ denotes the second derivative with respect to $\theta$. Inserting Equation (5) into Equation (1), by taking into account notation Equation (6), we have

$$({\phi }_{rr}+\frac{1}{r}{\phi }_r+\frac{1}{r^2}{\phi }_{\theta \theta })\Omega =\phi {\Omega }^{\prime }.$$

(7)

We further consider the following equation

$$\phi (r,\theta )=R(r)W(\theta ),$$

(8)

where $R(r)$ and $W(\theta )$ are functions. The function $\phi (r,\theta )$ is a product of two functions, including $R(r)$ and $W(\theta )$. Each function depends only on one of the variables r or $\theta$. Inserting Equation (8) into Equation (7), we obtain

$$R^{″}W\Omega +\frac{1}{r}{R}^{\prime }W\Omega +\frac{1}{r^2}R{W}^{″}\Omega =RW{\Omega }^{\prime },$$

(9)

where $R^{\prime }=\frac{dR(r)}{dr}$, $R^{″}=\frac{d^{2}R(r)}{d{r}^2}$, $W^{″}=\frac{d^{2}W(\theta )}{d{\theta }^2}$, and ${\Omega }^{\prime }=\frac{d\Omega (t)}{dt}$.

Dividing by $R(r)W(\theta )\Omega (t)$ on both sides in Equation (9), we can then obtain

$$\{\begin{array}{l}\frac{{\Omega }^{\prime }}{\Omega }=\lambda ,\\ \frac{r^2{R}^{″}+r{R}^{\prime }-r^{2}R\lambda }{R}=\chi \\ -\frac{W^{″}}{W}=\chi ,\end{array},$$

(10)

where $\lambda$ and $\chi$ are separation constants. We introduce p and q and assume that $\lambda =p^2$ or $\lambda =-p^2$ and $\chi =q^2$ or $\chi =-q^2$ to ensure $\lambda$ and $\chi$ to be positive or negative value, respectively. The formulation of Trefftz functions for transient moving boundary problems are expressed in the following description. Considering the combination of positive or negative values for $\lambda$ and $\chi$, there are six possible scenarios. If we consider the first scenario, $\lambda =0$ and $\chi =q^2$, we may obtain

$$\{\begin{array}{l}\Omega (t)=A_1,\\ R(r)=A_2{r}^q+A_3{r}^{-q},\\ W(\theta )=A_4\cos (q\theta )+A_5\sin (q\theta ),\end{array}$$

(11)

where $A_1$, $A_2$, $A_3$, $A_4$, and $A_5$ are arbitrary constants that have to be evaluated. Inserting Equation (11) into Equation (5) may yield

$$h(r,\theta ,t)={\overline{A}}_1{r}^q\cos (q\theta )+{\overline{A}}_2{r}^q\sin (q\theta )+{\overline{A}}_3{r}^{-q}\cos (q\theta )+{\overline{A}}_4{r}^{-q}\sin (q\theta ),$$

(12)

where ${\overline{A}}_1$, ${\overline{A}}_2$, ${\overline{A}}_3$, and ${\overline{A}}_4$ are arbitrary constants that have to be evaluated. We may find solutions for five other scenarios including $\lambda =0$ and $\chi =0$, $\lambda =p^2$ and $\chi =q^2$, $\lambda =p^2$ and $\chi =0$, $\lambda =-p^2$ and $\chi =q^2$, and $\lambda =-p^2$ and $\chi =0$ using the same procedure, as listed in Appendix A. As a result, we may obtain the complete Trefftz functions described as follows,

$$T=\left\{{\overline{\mathrm{T}}}_1,{\overline{\mathrm{T}}}_2,{\overline{\mathrm{T}}}_3, \dots ,{\overline{\mathrm{T}}}_{18}\right\},$$

(13)

where T denotes the Trefftz basis functions, and ${\overline{\mathrm{T}}}_1,{\overline{\mathrm{T}}}_2,{\overline{\mathrm{T}}}_3,\dots ,{\overline{\mathrm{T}}}_{18}$ denotes the functions, as listed in Appendix A. The transient numerical solution for the two-dimensional subsurface flow problem with a transient moving boundary is expressed by the series expansion as follows,

$$\begin{array}{cc}\hfill h(r,\theta ,t)=& \overline{a}+\overline{b}\ln r+{\displaystyle \sum _{q=1}^{\upsilon }\{\begin{array}{l}{\overline{c}}_{1q}{r}^q\cos (q\theta )+{\overline{c}}_{2q}{r}^q\sin (q\theta )+{\overline{c}}_{3q}{e}^{p^{2}t}{I}_0(qr)+{\overline{c}}_{4q}{e}^{-p^{2}t}{J}_0(qr)\\ +{\overline{c}}_{5q}{r}^{-q}\cos (q\theta )+{\overline{c}}_{6q}{r}^{-q}\sin (q\theta )+{\overline{c}}_{7q}{e}^{p^{2}t}{K}_0(qr)+{\overline{c}}_{8q}{e}^{-p^{2}t}{Y}_0(qr)\end{array}}\hfill \\ \hfill & +{\displaystyle \sum _{p=1}^{\upsilon }\left\{\begin{array}{l}{\overline{d}}_{1qp}{e}^{p^{2}t}{I}_q(pr)\cos (q\theta )+{\overline{d}}_{2qp}{e}^{p^{2}t}{I}_q(pr)\sin (q\theta )\\ +{\overline{d}}_{3qp}{e}^{-p^{2}t}{J}_q(pr)\cos (q\theta )+{\overline{d}}_{4qp}{e}^{-p^{2}t}{J}_q(pr)\sin (q\theta )\\ +{\overline{d}}_{5qp}{e}^{p^{2}t}{K}_q(pr)\cos (q\theta )+{\overline{d}}_{6qp}{e}^{p^{2}t}{K}_q(pr)\sin (q\theta )\\ +{\overline{d}}_{7qp}{e}^{-p^{2}t}{Y}_q(pr)\cos (q\theta )+{\overline{d}}_{8qp}{e}^{-p^{2}t}{Y}_q(pr)\sin (q\theta )\end{array}\right\}\begin{array}{l}\\ \\ \\ \\ \end{array}\},}\hfill \end{array}$$

(14)

where $\upsilon$ denotes the order of the Trefftz functions, and $\overline{a}$, $\overline{b}$, ${\overline{c}}_{1q}$ …, ${\overline{d}}_{8qp}$ denote unknown coefficients, $I_0$ and $I_q$ denote the modified Bessel functions of the first kind of zero order and of q order, respectively. $J_0$ and $J_q$ denote the Bessel functions of the first kind of zero order and of q order, respectively. $K_0$ and $K_q$ denote the modified Bessel functions of the second kind of zero order and of q order, respectively. $Y_0$ and $Y_q$ denote the Bessel functions of the second kind of zero order and of q order, respectively.

For the infinite domain or domain with cavities, the Trefftz basis functions are described as

$$\begin{array}{cc}h(r,\theta ,t)=\hfill & \overline{b}\ln r+{\displaystyle \sum _{q=1}^{\upsilon }\{ {\overline{c}}_{5q}{r}^{-q}\cos (q\theta )+{\overline{c}}_{6q}{r}^{-q}\sin (q\theta )+{\overline{c}}_{7q}{e}^{p^{2}t}{K}_0(pr)+{\overline{c}}_{8q}{e}^{-p^{2}t}{Y}_0(pr)}\hfill \\ \hfill & +{\displaystyle \sum _{p=1}^{\upsilon }\left\{\begin{array}{l}{\overline{d}}_{5qp}{e}^{p^{2}t}{K}_q(pr)\cos (q\theta )+{\overline{d}}_{6qp}{e}^{p^{2}t}{K}_q(pr)\sin (q\theta )\\ +{\overline{d}}_{7qp}{e}^{-p^{2}t}{Y}_q(pr)\cos (q\theta )+{\overline{d}}_{8qp}{e}^{-p^{2}t}{Y}_q(pr)\sin (q\theta )\end{array}\right\}\begin{array}{l}\\ \end{array}\}}.\hfill \end{array}$$

(15)

When the domain is simply connected, we may consider only positive basis functions. Consequently, the above equation is simplified as the following equation.

$$\begin{array}{cc}h(r,\theta ,t)=\hfill & \overline{a}+{\displaystyle \sum _{q=1}^{\upsilon }\{ {\overline{c}}_{1q}{r}^q\cos (q\theta )+ {\overline{c}}_{2q}{r}^q\sin (q\theta )+{\overline{c}}_{3q}{e}^{p^{2}t}{I}_0(qr)+{\overline{c}}_{4q}{e}^{-p^{2}t}{J}_0(qr)}\hfill \\ \hfill & +{\displaystyle \sum _{p=1}^{\upsilon }\left\{\begin{array}{l}{\overline{d}}_{1qp}{e}^{p^{2}t}{I}_q(pr)\cos (q\theta )+{\overline{d}}_{2qp}{e}^{p^{2}t}{I}_q(pr)\sin (q\theta )\\ +{\overline{d}}_{3qp}{e}^{-p^{2}t}{J}_q(pr)\cos (q\theta )+{\overline{d}}_{4qp}{e}^{-p^{2}t}{J}_q(pr)\sin (q\theta )\end{array}\right\}\begin{array}{l}\\ \end{array}\}}.\hfill \end{array}$$

(16)

To evaluate the unknown coefficients of $\overline{a}$, ${\overline{c}}_{1q}$, …, ${\overline{d}}_{4qp}$ in Equation (16), the spacetime collocation scheme must be utilized. Using the spacetime collocation scheme and applying the Dirichlet boundary data in Equation (16), a system of equations may then be yielded.

$$\left[\begin{array}{cccccc}1& r_1{}^q\cos (q{\theta }_1)& r_1{}^q\sin (q{\theta }_1)& & \cdots & e^{-p^2{t}_1}{J}_q(p{r}_1)\sin (q{\theta }_1)\\ 1& r_2{}^q\cos (q{\theta }_2)& r_2{}^q\sin (q{\theta }_2)& & \cdots & e^{-p^2{t}_2}{J}_q(p{r}_2)\sin (q{\theta }_2)\\ ⋮& ⋮& ⋮& & \ddots & ⋮\\ 1& r_s{}^q\cos (q{\theta }_s)& r_s{}^q\sin (q{\theta }_s)& & \cdots & e^{-p^2{t}_{\mathrm{s}}}{J}_q(p{r}_s)\sin (q{\theta }_s)\end{array}\right]\left[\begin{array}{c}\overline{a}\\ {\overline{c}}_{1q}\\ ⋮\\ {\overline{d}}_{4qp}\end{array}\right]=\left[\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_s\end{array}\right],$$

(17)

where $t_1, t_2,\cdots ,t_s$ are time in dimensionless form, $r_1, r_2,\cdots ,r_s$ are radiuses in dimensionless form, ${\theta }_1, {\theta }_2,\cdots ,{\theta }_s$ are polar angles in dimensionless form, $h_1, h_2,\cdots ,h_s$ are Dirichlet boundary data, the subscript s denotes the number of boundary points, and $\overline{a}, {\overline{c}}_{1q},\cdots ,{\overline{d}}_{4qp}$ denote the unknown coefficients. Equation (17) is expressed as

$$\mathbf{H}\mathbf{y}=\mathbf{Z},$$

(18)

where $\mathbf{H}$ denotes a matrix of the Trefftz basis functions with the size of $s\times w$, $\mathbf{y}$ denotes a vector of the unknown coefficients with the size of $w\times 1$, $\mathbf{Z}$ denotes a vector of accessible boundary value at boundary collocation points with the size of $s\times 1$, $s$ denotes the number of boundary points, $w$ denotes the term related to the order of the Trefftz basis function. Solving Equation (18), we may acquire the coefficients that are unknown for the spacetime domain. In addition, the Neumann boundary conditions are also considered in this study.

$$h_n(r,\theta ,t)=\frac{\partial h(r,\theta ,t)}{\partial n}=\nabla h(r,\theta ,t)\cdot \stackrel{\rightharpoonup }{n},$$

(19)

where $\stackrel{\rightharpoonup }{n}=(n_x,n_y)$ denotes the outward normal vector, $n_x$ and $n_y$ are the outward normal direction of the x and y axis, respectively. Adopting the chain rule, we may yield the formulations of $h_n$, $h_x$, and $h_y$, as listed in Appendix B.

3.2. The Spacetime Collocation Scheme

Instead of utilizing the original Euclidean space, we adopt a spacetime collocation scheme based on the Minkowski spacetime to perform the transient modeling of this problem. A two-dimensional transient nonlinear moving boundary problem is two-dimensional in space as well as one-dimensional in time, as displayed in Figure 1a. The spacetime region then becomes a domain in three dimensions, as shown in Figure 1b. To calculate the polar angle as well as the radius, we placed the source point as a reference point within the domain, as shown in Figure 1b. As a result, both the initial and boundary condition can be provided on the boundary of spacetime. Since the final time boundary data are unknown, the spacetime collocation scheme transforms a transient nonlinear moving boundary problem into an inverse BVP.

3.3. The Iterative Scheme for Modeling Transient Moving Boundary

For each collocation point on the moving surface, the total head is expressed as

$$h_{\phi }(r_j,{\theta }_j,t_j)=Y_j+\frac{p_j}{\gamma },$$

(20)

where $Y_j$ is the height above the sea level, $\gamma$ denotes the unit weight of water, $p_j$ denotes the pore water pressure, $h_{\phi }(r_j,{\theta }_j,t_j)$ is the total head, and the subscript j denotes the index of the points on the transient moving boundary to be renewed. The over-specified moving boundary conditions, including the no-flux and the zero pressure head, are described as

$$\frac{\partial h_{\phi }(r_j,{\theta }_j,t_j)}{\partial n}=0, h_{\phi }(r_j,{\theta }_j,t_j)=Y_j.$$

(21)

For each collocation point on the seepage face, the total head is expressed as

$$h_{\phi }(r_j,{\theta }_j,t_j)=Y_j.$$

(22)

As demonstrated in Equations (16) and (A11), the complete mathematical expressions of the Dirichlet and Neumann boundary conditions have been derived. Applying the Dirichlet and Neumann boundary values for boundary points on the moving surface may acquire

$$h_{\phi }(r_j,{\theta }_j,t_j)\approx {\displaystyle \sum _{q=1}^{\upsilon }{\displaystyle \sum _{p=1}^{\upsilon }{\mathbf{y}}_{qp}{\mathbf{H}}_{qp}(r_j,{\theta }_j,t_j)}},$$

(23)

$$\frac{\partial h_{\phi }(r_j,{\theta }_j,t_j)}{\partial n}\approx {\displaystyle \sum _{q=1}^{\upsilon }{\displaystyle \sum _{p=1}^{\upsilon }{\mathbf{y}}_{qp}\frac{\partial {\mathbf{H}}_{qp}(r_j,{\theta }_j,t_j)}{\partial n}}}.$$

(24)

From the above equations, the given boundary data are over-specified on the moving boundary. On the moving boundary, the location of the moving boundary is unknown. It can be referred to as the inverse geometric problem. For example, considering the no-flux and the zero pressure head boundary conditions, the unknowns are the coordinates of collocation points. From Equations (23) and (24), it is found that we may solve a nonlinear system of equations to obtain the coordinates of collocation points for the given time. The moving boundary problem may, therefore, exhibit the nonlinear characteristic. The inverse geometric problems are usually difficult to deal with because of the nonlinearity. For solving the inverse geometric problem, such as the moving surface flow problem, the iterative scheme is required. Previous studies have found it difficult to calculate the Jacobian matrix using Newton’s method. Thus, the Picard iterative method is used in this study [5]. The Picard iteration first begins from the initial guess of the location for the moving boundary. The iteration may be achieved by applying Equations (23) and (24).

$${\displaystyle \sum _{q=1}^{\upsilon }{\displaystyle \sum _{p=1}^{\upsilon }{\mathbf{y}}_{qp}^{{}^i}{\mathbf{H}}_{qp}^{}(r_j^i,{\theta }_j^i,t_j^i)}}=h^i(r_j,{\theta }_j,t_j),$$

(25)

$${\displaystyle \sum _{q=1}^{\upsilon }{\displaystyle \sum _{p=1}^{\upsilon }{y}_{qp}^{{}^i}\frac{\partial H_{qp}(r_j^i,{\theta }_j^i,t_j^i)}{\partial n}}}=\frac{\partial h^i(r_j,{\theta }_j,t_j)}{\partial n},$$

(26)

where $h^i(r_j,{\theta }_j,t_j)=Y_j^i$ and the superscripts i is the number of iteration steps. The iterative equation is depicted as

$$h^i(r_j,{\theta }_j,t_j)=h^{i-1}(r_j,{\theta }_j,t_j)+\epsilon (h^i(r_j,{\theta }_j,t_j)-h^{i-1}(r_j,{\theta }_j,t_j)),$$

(27)

where $h^i(r_j,{\theta }_j,t_j)$ is the total head to be renewed, and $\epsilon$ is the parameter of under-relaxation. The $\epsilon$ value is in the range of zero to one. The numerical procedure of the iteration starts by giving an initial value for the nonlinear moving boundary and ends while the stopping condition is achieved.

$$\frac{\sqrt{{\displaystyle \sum _{j=1}^J{\left[h^i(r_j,{\theta }_j,t_j)-h^{i-1}(r_j,{\theta }_j,t_j)\right]}^2}}}{\sqrt{{\displaystyle \sum _{j=1}^J{\left[h^{i-1}(r_j,{\theta }_j,t_j)\right]}^2}}}\le \omega ,$$

(28)

where $\omega$ is the stopping criteria, and J is the collocation point number on the moving boundary. In this study, we consider the stopping criteria to be $\omega ={10}^{-4}$.

4. Numerical Examples

4.1. Numerical Example 1

An example of the shape of a Cassini oval, as demonstrated in Figure 1a, is expressed as

$$\Re \in \left\{(x,y)|x=\widehat{r} \cos \theta ,y=\widehat{r}\sin \theta \right\},$$

(29)

where $\widehat{r}(\theta )=\sqrt[3]{\cos (3\theta )+\sqrt{2-\sin {(3\theta )}^2}} , 0\le \theta \le 2\pi$. The governing equation can be expressed as given in Equation (1). The harmonic data at the initial time is assumed as

$$h(x,y,\widehat{t}=0)=x^2+y^2.$$

(30)

The Neumann data are applied on the domain boundary $\mathsf{\Gamma }$, as displayed in Figure 1b. The Neumann boundary condition is considered as

$$\frac{\partial h(x,y,\widehat{t})}{\partial n}=2x\cdot n_x+2y\cdot n_y \mathrm{on} \mathsf{\Gamma }.$$

(31)

The following exact solution is adopted to validate the proposed method.

$$h(x,y,\widehat{t})=x^2+y^2+\frac{1}{4} \frac{T}{S}\widehat{t}$$

(32)

In this example, the storage coefficient S is ${10}^{-4}$, the transmissibility coefficient T is ${10}^{-5}$ m²/s, and final elapsed time ${\widehat{t}}_f^{}$ is 3 s. This example demonstrates space collocation points in two dimensions and time collocation points in one dimension. The spacetime collocation points can be regarded as a spacetime domain in three dimensions, as shown in Figure 1b. Due to the inaccessible final time boundary data, the two-dimensional initial value problem is transformed into a three-dimensional inverse BVP. The initial, as well as boundary data, are assigned on the circumferential and bottom sides of the spacetime domain in three dimensions, respectively. In this example, the source point is collocated on the origin as a reference point for calculating the polar angle and radius, as shown in Figure 1b.

The accuracy of the solution of our approach may be affected by the boundary collocation points number as well as the order of the Trefftz basis functions. A sensitivity analysis of the boundary collocation points number, and the order of the Trefftz basis functions is then carried out. To verify the stability of the proposed method, the accuracy of the numerical solution is measured by the following maximum absolute error (MAE).

$$\mathrm{MAE}=\max \left|u_E-u_N\right|,$$

(33)

where $u_E$ is the exact solution, and $u_N$ is the numerical solution.

Figure 2a demonstrates the relationship between the MAE and the order of the basis functions. It seems that accurate solutions are yielded when the order of the basis functions is greater than 10. Figure 2b depicts the number of boundary collocation points versus the MAE. It seems that accurate solutions may be achieved when the boundary collocation points number is greater than 700. Hence, the order of the basis functions and the boundary collocation points number are considered to be 11 and 918, respectively.

An example of a two-dimensional transient subsurface flow problem with harmonic initial and boundary conditions is then carried out to verify the computed result. To yield the computed total head and examine the accuracy of the proposed method, 3158 inner collocation points are uniformly collocated within the domain. The profiles of the computed total head are then chosen to compare with the analytical solution. Figure 3 depicts the exact solution as well as the field solution of the total head from the proposed method. It seems that by utilizing our method, the computed results are entirely consistent with the analytical solution. Figure 4 depicts the MAE of our method at different times. The MAE of our method is in the order of 10⁻¹³, as displayed in Figure 4. It is clear that our method may yield accurate results.

4.2. Numerical Example 2

This example is the forward and backward analyses of a two-dimensional transient problem. The governing equation is expressed in Equation (1). The non-harmonic initial and boundary data are considered. In this example, we assume the final elapsed time to be 1 min, the storage coefficient to be ${10}^{-4}$, the transmissibility coefficient to be ${10}^{-6}$ m²/s, the length to be 10 m, and the width to be 6 m. We apply the non-harmonic initial and boundary data as

$$h(x,y,\widehat{t}=0)=100\sin x\sin y,$$

(34)

$$h(x,y,\widehat{t})=100e^{-\frac{3T}{S}\widehat{t}}\sin x\sin y.$$

(35)

Because numerical example 2 may not exist an analytical solution to examine the accuracy, we conduct the forward modeling of the two-dimensional transient problem to compute the final time field solution of total head. The non-harmonic boundary data can be provided on vertical sides of the spacetime domain, as displayed in Figure 5a. A backward analysis using the final time results from the forward analysis, as displayed in Figure 5b, is then carried out to compute the field solution of initial head at the bottom of the spacetime domain. To verify the correctness of the field solution, the assigned non-harmonic initial data is compared with the computed initial head from the backward analysis of this problem.

In this example, there exists one source point collocated on the origin and 600 boundary points uniformly collocated on the boundary. The order of the Trefftz function is 21. To yield the computed total head and examine the accuracy of the proposed method, 1000 inner collocation points are uniformly collocated within the domain. Figure 6 shows the computed initial head with the assigned non-harmonic initial data at different times. It is found that the computed initial data from the backward analysis are consistent with the assigned non-harmonic data at time zero. Moreover, the absolute difference is in the order of ${10}^{-7}$, as given in Figure 7.

4.3. Numerical Example 3

The modeling of a two-dimensional transient moving boundary problem through a rectangular dam is considered, as depicted in Figure 8. The objective of this two-dimensional problem is to evaluate the time-dependent location of the phreatic surface. The governing equation is expressed in Equation (1). The rectangular dam, as shown in Figure 8, is constituted of five boundary lines, including ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, ${\mathsf{\Gamma }}_3$, ${\mathsf{\Gamma }}_4$, and ${\mathsf{\Gamma }}_5$. We assume the downstream head, upstream head, and the width to be 4, 24, and 16 m, respectively. The initial condition is expressed as

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

(36)

The Dirichlet data are applied on the domain boundary, including ${\mathsf{\Gamma }}_2$, ${\mathsf{\Gamma }}_3$, ${\mathsf{\Gamma }}_4$, and ${\mathsf{\Gamma }}_5$, as depicted in Figure 8. The boundary conditions are expressed as

$$h(x,y,\widehat{t})=4 \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_2,$$

(37)

$$h(x,y,\widehat{t})=24 \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_5,$$

(38)

$$h(x,y,\widehat{t})=y \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_3 \mathrm{and} {\mathsf{\Gamma }}_4.$$

(39)

On ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_4$, no-flow Neumann boundary data can be given as

$$\frac{\partial h(x,y,\widehat{t})}{\partial n}=0 \mathrm{on} {\mathsf{\Gamma }}_1 \mathrm{and} {\mathsf{\Gamma }}_4.$$

(40)

In this example, we assume that the final elapsed time is 700 days, storage coefficient is ${10}^{-3}$ and transmissibility coefficient is ${10}^{-6}$ m²/s. There exists one source point, where the location of the source point is (0,12). The order of the Treffttz basis functions and boundary collocation points number on a nonlinear moving boundary are 7 and 25,600, respectively. Since the numerical procedure for evaluating the location of the phreatic surface is considered as an inverse problem, the location of the separation point has to be investigated. To obtain the results, the number of iterations is 43 using the proposed method. Figure 9 shows the numerical solutions of the nonlinear moving boundary at different times. It seems that the transient nonlinear moving boundary can be obtained by utilizing our method.

Since several numerical approaches have been applied to solve this problem, we further compare the final time solutions of our method with those of Aitchison (1972) [31], Chen et al. (2007) [9], and Ku et al. (2019) [23], as depicted in Figure 10. It is found that the location of the nonlinear moving boundary using the proposed method agrees very well with the results from the previous studies at the final steady-state time.

4.4. Numerical Example 4

The final example is the modeling of a two-dimensional transient moving boundary problem through a trapezoidal dam, as displayed in Figure 11. The objective of this example is to evaluate the position of the transient moving boundary. The governing equation is expressed in Equation (1). The initial data is given as

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

(41)

The Dirichlet data are applied on the domain boundary, including ${\mathsf{\Gamma }}_2$, ${\mathsf{\Gamma }}_3$, ${\mathsf{\Gamma }}_4$, and ${\mathsf{\Gamma }}_5$, as depicted in Figure 11. The boundary conditions are described as

$$h(x,y,\widehat{t})=3 \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_2,$$

(42)

$$h(x,y,\widehat{t})=4 \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_5,$$

(43)

$$h(x,y,\widehat{t})=y \mathrm{m} \mathrm{on} {\mathsf{\Gamma }}_3 \mathrm{and} {\mathsf{\Gamma }}_4.$$

(44)

On ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_4$, the no-flow Neumann boundary data are given as

$$\frac{\partial h(x,y,\widehat{t})}{\partial n}=0 \mathrm{on} {\mathsf{\Gamma }}_1 \mathrm{and} {\mathsf{\Gamma }}_4.$$

(45)

In numerical example 4, the final elapsed time is 30 min, storage coefficient is ${10}^{-4}$, and transmissibility coefficient is 10⁻⁶ m²/s. The order of the Treffttz basis functions and number of boundary collocation points on the transient moving boundary are set to be 10 and 9663, respectively. There exists one source point, where the location of the source point is (1.2,2). Figure 12 shows the proposed spacetime collocation scheme of the two-dimensional transient moving boundary flow through a trapezoidal dam.

The profiles of the field solutions at different simulation times are chosen to clearly view the computed transient moving boundary. Figure 13 shows the numerical solutions of the transient moving boundary. The location of the moving surface at the final time from the proposed method is further compared with the steady-state solution by adopting the MFS to examine the accuracy of the proposed approach, as depicted in Figure 13. It is found that the location of the nonlinear moving surface using the proposed method agrees well with the steady-state solution by using the MFS.

5. Conclusions

This study is rooted in the Trefftz method and gives a promising numerical solution for the transient nonlinear moving boundary problems. To verify the proposed spacetime collocation scheme using Trefftz functions, we carried out several numerical problems. The key contributions of this study are as follows.

Previous studies have demonstrated that the engineering application of the Trefftz method with complete Trefftz functions for dealing with transient problems is still hardly found, where solving transient moving boundary problems using the Trefftz method rarely even exists. In this study, a pioneering attempt reveals that the transient nonlinear moving boundary problems governed by the two-dimensional diffusion equation are solved using the spacetime collocation scheme with complete Trefftz basis functions.

The significance of the proposed method rooted in the conventional Trefftz method is that the collocation points in our method are placed in the Minkowski spacetime rather than the Euclidean space. As a result, we may construct a spacetime domain in three dimensions, where both the boundary and initial data are given on the boundary of spacetime, which can be regarded as a BVP. Accordingly, the transient nonlinear moving boundary problems can be easily solved.

Because our method is a boundary-type meshless approach, the domain boundary has to be discretized by the boundary points. It depicts the simplicity of using the proposed method for dealing with problems of the transient moving boundary during the iterative process for searching the location of the free surface.