Free surface groundwater flow solution using boundary collocation methods

In this paper, two meshless numerical algorithms are developed for the solution of two-dimensional steady-state diffusion equation that describes the stationary groundwater flow. The proposed numerical methods, which are truly meshless, quadrature-free and boundary only, are based on the method of fundamental solutions and singular boundary method respectively. The diffusion equation is transformed into a Poissontype equation with a known fundamental solution. Numerical examples with moving boundary are presented and compared to the solutions obtained by the finite element method.


Introduction
Groundwater is one of the important aspects of the geotechnical design that the engineers have to deal with.That is the reason why the investigation of the numerical models for groundwater problems is an important task for geotechnical engineers and mathematicians.The groundwater flow models are solved by numerical techniques which have been developed to solve partial differential equations.One of the widespread numerical models -MODFLOW is using finite difference method (FDM) under the hood.The other wellknown model FEFLOW is based on the most widely used numerical method in engineering practice -finite element method (FEM).
Corresponding to these traditional methods, boundary-type methods have attracted considerable attention in recent decades, since methods of this type avoid the tedious domain mesh generation and reduce the geometrical dimension of the model by one.The flagship of a fore above-mentioned family of numerical methods is the boundary element method (BEM).Main drawbacks of these methods are the need for a fundamental solution of the differential operator to know, integration of the fundamental solution which contains the singularity around the source point and from the computational point of view, a solution of resulting fully populated, non-symmetric characteristic matrix.The boundary-type methods based on the fundamental solution are superior to the classical FEM or FDM formulations in the area of groundwater numerical models even in the case of the highly non-homogeneous domain of interest.The recent development of boundary-type methods is focused on the characteristic matrix assembly and solution scheme [1] and meshless approach to avoid creation of the boundary mesh for interpolation or integration purposes.
The meshless boundary type methods are based on the collocation of the fundamental solution and superposition principle to express the final solution as the linear combination of the fundamental solutions similarly to indirect BEM formulations.The promising advantage of this method is very easy programming and the method can be easily developed using few lines of code in Matlab or Octave, while it still maintains the advantage of boundary only formulation.Advantages of such approaches are even magnified when 3D problems come into play with complex boundaries; the creation of just boundary nodes is much more flexible and adaptive than creating the mesh of elements that have to obey some sort of tessellation.The problem for the developer of the method comes into play when he needs to deal with the proximity of fundamental solution singularity.In the classical BEM, thanks to the weak formulation, this singularity is treated using one of the strategies for integrating improper integrals.The origin intensity factor, can't be directly evaluated using fundamental solution so other techniques like moving source nodes outside the computational domain -method of fundamental solution [2], empirical formulas based on the mean value of boundary integral in the vicinity of source point -singular boundary method [3] or simply use of the known regular general solution for problem close enoughboundary knot method [4].
In this article, the MFS and SBM numerical schemes are presented for the solution of free-surface flow problems and the results are compared with the solutions of the same problems solved using FEM from the view computational resources and precision of results which are presented.

Method of fundamental solutions -MFS
The MFS implementation, as was mentioned above, is very easy and only a little boundary data preparation is required in terms of boundary collocation or boundary fitting which is simpler than the domain-type methods Although the MFS has been successfully applied to a variety of complicated problems [2,5] there have been very little applications of this method investigation of groundwater flow problems.
The ideas behind the MFS have been around for a few decades and were developed primarily by V. D. Kupradze and M. A. Alexidze in the late 1950s and early 1960s [2].
The MFS formulation is quite straightforward; consider a 2D Laplace partial differential equation we seek for its solution over domain Ω with Dirichlet and Neumann boundary conditions specified over part of the boundary of the computational domain ( , ) ( , ) Substituting the approximation into boundary condition yields the following matrix equation The unknown coefficients α can uniquely be determined by the above algebraic equation and then we can evaluate numerical solution at any location in the physical domain.
The 2D groundwater flow governing equation is expressed as where we denote the hydraulic conductivity, in the direction of the Cartesian axis, as K x and K y respectively.The h represents the hydraulic head and S is storativity.Because all examples in this article represent steady-state flow, the unstable term ∂h/∂t=0 making equation ( 1) homogeneous.To transform (4) to Laplace equation, we use following the transformation of coordinates and assuming constant hydraulic conductivities and we obtain The fundamental solution of the two-dimensional Laplace equation is when the source (s) and response points (x) are coincident the fundamental solution becomes singular, the diagonal terms in (3) can't be evaluated, and there is a need for a special strategy to deal with singularities.The singularity of the fundamental solution in the MFS is tackled using very simple approach -the source point is moved outside the computational domain -the fictitious boundary that doesn't need to have the same shape as the computational domain is created (See Fig. 1).Based on the nature of the fundamental solution the close proximity of fictitious boundary brings the influence of fundamental solution singularity and the source points placed too far from computational domain leads to the ill-conditioned characteristic matrix.
The fictitious boundary in the MFS is problematic because its position affects the results and there is no rigorous way how to choose the position for this boundary.Various empirical approaches have been developed [3], but none of them have general applicability.Some new techniques have recently been developed to cure this problem of MFS, such as the boundary knot method, singular boundary method, and regularized meshless method.

Singular boundary method -SBM
The singular boundary method (SBM) is based on the same principle of the method of fundamental solutions (MFS) but SBM employs the same set of collocation and source points placed on the boundary Γ.As in MFS, the solution of homogeneous Laplace equation is approximated by the linear combination of fundamental solutions.Then the singular term arises when the source and collocation points are coincident [3]  where  j are unknown coefficients, G(x,s) is the fundamental solution of Laplace equation ( 7) and u ii , q ii are defined as origin intensity factors (OIFs).Origin intensity factor q ii for the Neumann boundary condition can be calculated using a subtracting and adding-back technique where L j is the length of the appropriate part of the boundary around point x (see Figure 2).In the case of the Dirichlet boundary condition, there is no rigorous way h ofow to evaluate corresponding OIF.One of the best possibilities how to solve 2D Laplace equation is the empirical formula [3] 4 Numerical examplesfree surface flow

Analysis of plane free surface flow
In this example, we consider a boundary value problem frequently met in geomechanics in relation to the flow of water through dams.Free surface problems involve an upper boundary, the location of which is not known a priori so an iterative procedure is required to find it.The analysis starts by assuming an initial position for the free surface.A solution of Laplace's equation gives values of the total head along the free surface nodes which will not in general equal the elevation of the upper surface of the model.The elevations of the nodes along the upper surface are therefore adjusted to equal the total head values just calculated at those locations.The boundary condition for this BVP is presented on Figure 3 and Figure 4.The problem was solved using the finite element program [6] and MFS.
The FEM solution is obtained using 4 node bilinear elements.The corresponding solution obtained using MFS geometrically similar to the boundary of FEM mesh.The solution obtained using boundary method shows relatively good agreement with the FEM solution.
The iterative solution process performed by MFS showed instability and precision problems because the results were affected by the problematic position of the fictitious boundary which needs to be moving with deformation of the computational boundary in some way.
The new position of fictitious boundary nodes brings the instability into the solution and shows the problem of MFS in the case of complex geometrical conditions (Figure 4).Finally, the overall relative difference between MFS and FEM solution was 4.51×10 -2 .

Free-surface flow through an earth dam
A second example of an earth dam with sloping sides is presented in Figure 5 with a description of boundary conditions.The initial mesh is trapezoidal and starts with a horizontal free surface set at an elevation of 37.5 m, which is also the height of the starting mesh.The nodes on the upstream face of the dam are also set at a total head of 37.5 m, while the bottom two nodes on the downstream side are fixed at a total head of 7.5 m.The initial mesh is defined by the x-coordinates of the nodes at the base and the top.The solution obtained by SBM (Figure 6 bottom) agrees with the FEM solution (Figure 6 top) and the overall relative difference between free surface computed by FEM and SBM is 5.65×10 -3 .The SBM characteristic matrix is fully populated but only 34 equations have to be solved in contrast to FEM solution where the coarser mesh was used (6-noded elements) and 72 equations have to be solved.The SBM method showed much more stable behavior than MFS which doesn't show properties suitable for solving moving boundary problems.

Conclusions
This paper documents the attempt to show boundary collocation numerical methods MFS and SBM applied to free-surface flow problem.The two types of free-surface boundary value problems have been successfully tested.This type of methods seems to be a quite useful alternative to the solutions using the BEM and FEM, however the MFS needs to handle the two moving boundaries (domain boundary and fictitious boundary) to maintain the stability of the solution and this drawback with no rigorous way how to position the fictitious boundary, makes the MFS unsuitable for the problems with moving boundariesfree surface flow.The SBM, even though is similar to the MFS, is much more compact and stable and much more suitable for the free-surface flows, but the problem of Dirichlet OIF determination is handled using an empirical formula which can't be simply extended for the 3D problems.In the post-processing phase, the MFS's absence of source points in the computational domain makes the result smooth and without any influence of singularity.However the SBM have source points positioned on the boundary of the computational domain so results near boundary are affected by singularity, but this problem can be overcome using IIT technique to obtain results even near boundaries.The applications of SBM or MFS on a particle-driven flow and on a debris flow could be also very interesting.
This contribution is the result of the project funded by the Scientific Grant Agency of Slovak Republic(VEGA) No. 1-0716-17.

Fig. 1 .
Fig. 1.Distribution of boundary points and source points for MFS.

Fig. 2 .
Fig. 2. Part of boundary L around point x .

Fig. 4 .
Fig. 4. Boundary conditions, geometrical configuration and converged solution mesh obtained using MFS (The dotted lines represent boundary normals and plus signs are positions of source nodes).

Fig. 6 .
Fig. 6.Converged solution and final deformation of the computational mesh computed by FEM(top) and SBM (bottom).