Green ’ s function based finite element formulations for isotropic seepage analysis with free surface 1

A solution procedure using the Green’s function based finite element method (FEM) is presented for two-dimensional nonlinear steadystate seepage analysis with the presence of free surface in isotropic dams. In the present algorithm, an iteration strategy is designed to convert the over-specified free surface problem to a regular partial differential equation problem. Then, at each iteration step, the Green’s function for isotropic linear seepage partial differential equation is employed to construct the element interior water head field, while the conventional shape functions are used for the independent element frame water head field. Then these two independent fields are connected by a double-variable hybrid functional to produce the final solving equation system. By means of the physical definition of Green’s function, all two-dimensional element domain integrals in the present algorithm can reduce to one-dimensional element boundary integrals, so that versatile multinode element is constructed to simplify mesh reconstruction during iteration. Finally, numerical results from the present Green’s function based FEM with isotropic Green’s function kernels are compared with other numerical results to verify and demonstrate the performance of the present method.

Green's function based finite element formulations for isotropic seepage analysis with free surface

INTRODUCTION
During the treatment of water hazard, one of practical measures is to construct dams.So the drawdown seepage analysis is fatal for the security of dam subjected to water level change.One of the main difficulties is that the prediction of the free surface, which is usually over-specified.
Due to the complex boundary conditions and geometrical domain, numerical methods including FEM (Neuman, 1973;Bathe and Khoshgoftaar, 1979;Bathe, 2006;Luo et al., 2008) and boundary element method (BEM) (Brebbia and Chang, 1979;Chen et al., 1994;Tsay et al., 1997) are often used to analyze seepage problems of the dam with free surface.For example, Bath and Khoshgoftaar developed a FE solution procedure for nonlinear free surface seepage analysis without mesh iteration (Bathe and Khoshgoftaar, 1979).Luo et al. applied the FEM to establish three dimensional computational model of transient seepage for deep foundation pit dewatering in the Yangtze River Delta (Luo et al., 2008).Brebbia and Chang established boundary element formulation for seepage problem in anisotropic soils (Brebbia and Chang, 1979).Chen et al. constructed boundary element formulation for seepage problems by introducing a dual integral equation with hypersingular integral (Chen et al., 1994).In addition to the FEM and BEM mentioned above, meshless/meshfree methods including the method of fundamental solutions (MFS) (Young et al., 2006;Chaiyo et al., 2011), the natural element method (NEM) (Jie et al., 2013), the hybrid boundary node method (Tan et al., 2010), the smoothed fixed grid finite element method (Kazemzadeh-Parsi and Daneshmand, 2012) and the scaled boundary finite-element method (Bazyar and Talebi, 2014) were recently developed for the numerical determination of the free surface in seepage problems.
Generally, the iteration procedure is necessary for the free surface problems.However, the mesh reconstruction during iteration increases extremely the difficulty of solution.To bypass this difficulty, development of super element (or multi-node large element) is desired.In this study, another numerical method different to the conventional FEM and BEM (Bathe, 2006;Qin, 1994;1995;2003;Qin and Mai, 2002), called the fundamental solution-based (or Green's function based) hybrid FEM (HFS-FEM), is formulated for solving such problems in two-dimensional isotropic dams and a multi-node element is developed to model the region close to the free surface for simplifying the mesh redefinition.The HFS-FEM was firstly presented by Wang and Qin for heat transfer analysis (Wang and Qin, 2009) and then was extended to analyze elastic stress field (Wang and Qin, 2010b;Wang and Qin, 2011a;Wang and Qin, 2012b), thermal properties of advanced functional/composite materials (Wang and Qin, 2011b;Wang et al., 2012;Wang et al., 2013) and bioheat transfer in biological tissues (Wang and Qin, 2010a;Wang and Qin, 2012a) with general/special elements to achieve the purpose of high accuracy and mesh reduction.It should be mentioned that convergence of the HFS-FEM was fully discussed in these works.Different to the conventional FEM, the HFS-FEM involves element boundary integrals only, and allows for constructing arbitrary n-sided polygonal elements and accurately calculating fields everywhere in the domain, due to the use of fundamental solutions.Besides, both the HFS-FEM and the BEM are dependent of fundamental solutions, however, the BEM seems more complicated for solving multi-material problems, because it requires separate boundary integral equation for each material domain and additional connective equations for satisfying the continuous conditions between adjacent material domains.In this study, the HFS-FEM was extended for solving two-dimensional nonlinear seepage problems with free surfaces.In the present hybrid FE algorithm, the linear combination of Green's function of isotropic seepage governing equation is employed to construct the element interior water head field by placing source points outside the physical boundary of the domain, as done in the classical MFS (Chen et al., 2008), while the conventional shape functions are used for the element frame water head field, which is independent of the element interior field.Then these two fields are connected by a double-variable hybrid functional.Due to physical definition of the Green's function, the element interior field can naturally satisfy the isotropic seepage governing equation and thus, all element domain integrals in the functional are converted Latin American Journal of Solids and Structures 12 (2015) 1991-2005 into nonsingular element boundary integrals, which have lower dimensions than domain integrals.Most importantly, one can flexibly construct large elements with more nodes and edges than the conventional plane elements for practical use.In the paper, a large multi-node element is designed to model the region nearby the free surface to alleviate difficulty of unavoidable mesh reconstruction during iteration procedures.
The outline of the present work is as follows: Details of the seepage boundary value problem are provided in Section 2, and the present Green's function based FEM is formulated in Section 3. Subsequently, numerical experiments are conducted by the present method in Section 4 and some conclusions are finally presented in Section 5.

Governing equation and boundary conditions
Generally, the governing partial differential equation of two-dimensional incompressible steadystate flow through isotropic dams can be described as: where  is a bounded domain in two-dimensional space 2  refering to the rectangular coordinate H is the unknown water head in the domain  , and 0 0 k  denotes isotropic hydraulic conductivity or permeability coefficient of materials.
According to the Darcy's law, the seepage velocity components ( 1, 2) i q i  can be written in terms of the hydraulic head H, that is, Besides, generalized boundary conditions are of the type where H and q are respectively the known boundary data on the boundary H  and q  .q is the velocity in the direction normal to the boundary, and ( 1, 2) i n i  are the direction cosines of the normal with respect to the global coordinate system.
For the convenience of following derivation, the normal velocity q is rewritten in matrix form as with Latin American Journal of Solids and Structures 12 (2015) 1991-2005 Figure 1 shows a two-dimensional seepage problem with free surface in a dam, in which the material is assumed to be isotropic.AE and DG are the boundaries with specified constant water heads 1 H and 2 H , respectively.The bottom surface AD of the dam is assumed to be impermeable, that is 0 q  .The curve EF is the free surface and F is the separation point.On the free surface EF, two boundary conditions should be satisfied simultaneously: and 0 q  , so EF is overspecified.On the seepage surface FG, there is a given water head in terms of the coordinate variable 2 x , i.e.

H
x  .Since the water outflows the dam on this boundary, it should also satisfy the condition of 0 q  on FG.In summary, the boundary sections AE, EF, FG, GD, and AD consist of the complete the boundary of the solution domain.For the sake of clarity, all boundary conditions related to the seepage shown in Figure 1 are given as

Iteration of the free surface
Since the free surface problem is a nonlinear problem, it is inevitable to perform iteration in determining the location of the free surface displayed in Figure 1.
Latin American Journal of Solids and Structures 12 (2015) 1991-2005 (1) Firstly, an initial guess of the location of free surface is assumed to start the computation.For example, a straight line of EF can be initially given by setting the separation point F as the midpoint of the line CG to determine the solution domain AEFDA (see Figure 2).
(2) During iteration, only the condition 0 q  is specified on the free surface EF and the linear water condition is given on the seepage surface FG.
(3) After each step of iteration, the water head at the nodes on the free surface can be obtained using the proposed numerical method.
(4) Examining the error norm defined by where 2i x are the nodal vertical coordinates on the free surface and i H the nodal water heads on the free surface.
(5) Updating the nodal vertical coordinates 2i x on the free surface by setting

HYBRID FINITE ELEMENT FORMULATION
In the present Green's function based FE formulation, two independent element fields are separately defined and then are linked by a hybrid functional, from which the stiffness equation in terms of unknown nodal potential and the relationship of these two independent element fields can be obtained.

Hydraulic head within the element
In the interior of an element, say element e , the hydraulic head field can be approximated by a linear combination of Green's function as where s N is the number of source points locating outside the element domain.
( 1 ) s s c s N   present unknown interpolation coefficients and * ( , ) s H P Q stands for the isotropic Green's function related to the field point P and the source point s Q , which locates outside the element to avoid the singularity of the Green's function, as done in the MFS (Chen et al. and Karageorghis et al., 2008;Chaiyo et al. and Rattanadecho et al., 2011).Compared to the MFS, which employs the boundary discretization of the entire domain, the present method divides the entire domain into several small elements and in each element, the linear combination of Green's function is employed to represent the interior hydraulic head.Thus, the present method has better interpolation stability than the MFS and the location of the source points can be flexibly arranged outside the element domain (Wang and Qin, 2009;Wang and Qin, 2010a;Wang et al., 2013).
In Eq. ( 9), the vectors N and e c are respectively and Differentiating Eq. ( 9) yields the following normal seepage flow where The Green's function for a two-dimensional isotropic seepage problem is defined as whose solution is is an Euclidian distance between the field point 1 2 ( , ) P x x and the source point and  is the Dirac delta function.
Furthermore, the derivatives of the Green's function in terms of coordinate variables can be written as Obviously, the assumed interior water head field (9) analytically satisfies the governing equation ( 1), and this feature will be used to simplify the hybrid functional.

Hydraulic head along the element boundary
In order to enforce conformity on inter-element boundary between two adjacent elements, the frame water field H  is independently prescribed in terms of element nodal hydraulic head e d by ( ) where N  is an interpolation shape function vector, in which the interpolation functions are chosen to be the conventional quadratic shape functions commonly used in the isoparametric quadratic line elements in the FEM (Bathe, 2006) and the BEM (Brebbia, 1982).

Double-variable hybrid functional
For the particular element e , which occupying a sub-domain e  with boundary e  , the hybrid functional in terms of H and H  is defined as (Qin, 2000)   Applying the Gauss theorem, one has can be produced.In Eq. ( 25), the symmetric and sparse coefficient matrix is written by 4 NUMERICAL EXPERIMENTS

Verification of the present method
To verify the accuracy of the present approach, let's consider a linear isotropic seepage problem in a rectangular dam by fixing the free surface.The rectangular dam is 24 m in height and 16 m in width.The upstream (left) water head is 24 m, and the downstream (right) water head is 4 m, as shown in Figure 3.The permeability coefficient is . The free surface is assumed to be an inclined straight line connecting the points (0, 24) and (16, 14), which is the initial guess of the free surface in the next example.
In the computation, the FE mesh shown in Figure 3 is used to model the rectangular dam.In the meshing scheme, three 8-node general elements is used to model the bottom region of the dam and one 22-node super element is employed in the region close to the free boundary to provide more efficient remeshing process than the conventional FEM during iteration by updating only the locations of nodes on the free surface and the seepage surface.Besides, the use of large element is also beneficial to decrease the computation time, in contrast to the conventional FEM.To demonstrate modeling efficiency and solution accuracy of the present method, numerical results from the conventional FEM implemented by ABAQUS are also provided.Figures 4 and 5 respectively show the variation of the water head on the impermeable bottom surface and the straight-line free surface.It is observed that there is a good agreement between the present HFS-FEM and the FEM although much less elements are used in the proposed element model.It can be also seen from Figure 5 that there is an sudden jump of the water head at those points close to the seepage surface, due to the restriction of natural boundary condition at the separation point, so that curve fitting technology should be used to obtain a smooth free surface and to update the location of the separation point.Finally, the contour plots from the conventional FEM and the present method for the water head are presented in Figure 6, and similar distribution can be found for both methods.Therefore, the present method is efficient and accurate to simulate the distribution of the water head in the domain.In this example, the rectangular dam in the first test is considered again to determine the final location of the free surface by mean of iteration approach.In order to start the iteration, the inclined straight line given in Figure 3 is selected as an initial guess.The error tolerance in the iteration process is taken as 0.01.
The iteration results are shown in Figure 7.After 12 iterations, the iteration converges.Figure 8 shows the converged configuration of the free surface from the present method.Also, the results from the BEM (Chen et al., 2007) and the MFS (Chaiyo et al., 2011) are provided in Figure 8 for comparison.We can see that there is no significant difference on the computation results between the present method and other methods.The small deviation at the separation point is mainly caused by different curve fitting techniques.In the developed algorithm, the cubic polynomial interpolation is employed.For the sake of clarity, the height of the separation point respectively predicted by Latin American Journal of Solids and Structures 12 (2015) 1991-2005 the present HFS-FEM, the BEM, the MFS, and the FDM (Aitchison, 1972) are listed in Table 1, and it's found that they agree very well.Finally, the converged shape of the dam and the corresponding contour plot for the water head are listed in Figure 9.As another example in the literature (Jie et al., 2013), a rectangular dam with 4 m in width and 6 m in height is considered and its geometrical size is shown in Figure 10.The water levels upstream (left) and downstream (right) are 6 m and 1 m, respectively.The permeability coefficient is 0 k  0.1 / m d .The error tolerance in the iteration process is taken as 0.01.The example is analyzed by using the present HFS-FEM.Figure 11 displays the initial computational mesh, in which three general 8-node hybrid elements and one 22-node hybrid element are employed to model the domain.After 11 iterations, the mesh of the computing domain is updated to a convergent configuration, which is displayed in Figure 11 too.Figure 12 plots the final location of the free surface predicted from the present method and the test (Jie et al., 2013), and it can be found that a good agreement is achieved between them.This dam has also been analyzed by Jie et.al.using the natural element method (NEM) (Jie et al., 2013).For the purpose of comparison, the heights of the separation point from different methods are tabulated in Table 2, from which it can be seen that the present method has smaller derivation to the test value.

CONCLUSIONS
As a kind of domain-type numerical methods, the present hybrid FEM with Green's function kernels needs only element boundary integrations and versatile element construction can be fulfilled in the HFS-FEM for your practical needs.Thus, the present HFS-FEM is more suitable for dealing with problems with moving boundaries than the conventional FEM.In this study, the HFS-FEM is formulated for the determination of the free surface and the height of the separation point in isotropic dams, and a large multi-node element is employed to model the region close to the free surface for the purpose of simplifying mesh reconstruction during iteration process, so that only nodes on the free surface and the seepage surface are updated during iterations.The numerical results demonstrate that the developed Green's function-based hybrid FEM is capable to calculate free surface with good accuracy and efficiency.

Figure 1 :
Figure 1: Two-dimensional seepage problems with free surface in an isotropic dam.
The polynomial interpolation of the vertical coordinates of nodal points on the free surface is used to determine the new location of the separation point F.(7) Updating the nodal coordinates and linear water head boundary conditions on the seepage surface.(8)The iteration process can be terminated if the error  or the difference between ( ) to or smaller than a given tolerance.

Figure 2 :
Figure 2: Initial guess of the free surface.
element boundary integrals are involved.The substitution of the interior field (9) and the frame field (19) into Eq.(21) yields system of equations in terms of the nodal hydraulic head vector e d

Figure 3 :
Figure 3: Rectangular dam with straight free surface and Mesh strategy in the HFS-FEM.

Figure 4 :Figure 5 :
Figure 4: Distribution of water head along the bottom surface of the dam.Latin American Journal of Solids and Structures 12 (2015) 1991-2005

Figure 6 :
Figure 6: Contour plots for the water head in the dam.

4. 2
Isotropic seepage in a rectangular dam with size 24 m x 24 m

Figure 7 :
Figure 7: The change of location of the free surface during iterations.

Figure 8 :Table 1 :Figure 9 :
Figure 8: The convergent location of the free surface.

Figure 10 :
Figure 10: Rectangular dam with a height of 6 m and a width of 4 m.Latin American Journal of Solids and Structures 12 (2015) 1991-2005

Figure 12 :
Figure 12: Final location of the free surface.

Table 2 :
The height of the separation point.