func-MESHLESS ANALYSIS OF AN EMBANKMENT USING LOCAL GALERKIN RADIAL POINT INTERPOLATION METHOD (LGRPIM) MESHLESS ANALYSIS OF AN EMBANKMENT USING LOCAL GALERKIN RADIAL POINT INTERPOLATION METHOD (LGRPIM)

The paper deals with use of the meshless method for slope stability analysis. There are many formulations of the meshless methods. The article presents the Local Galerkin Radial Point Interpolation method (LGRPIM) – local weak formulation of the equilibrium equations. The main difference between meshless methods and the conventional finite element method (FEM) is that meshless shape functions are constructed using ran-domly scattered set of points without any relation between points. The shape function construction is the crucial part of the meshless numerical analysis in the construction of shape functions. The article presents the radial point interpolation method (RPIM) for the shape functions construction. The numerical example of the slope stability was calculated using meshless computer code and compared with FEM results.


Introduction
The investigation of soil behavior through the field observations or laboratory modeling is usually expensive and time-consuming [1]. So it is profitable to develop numerical solutions which may be used to predict the load-deformation characteristics and the stability of the soil structures [2]. The most popular method for numerical analyses in geotechnical engineering is the finite element method [3]. Although finite element method -FEM is widely used, there are some deficiencies such as the discontinuity of stresses on element boundaries or low accuracy at large deformation analysis. The mentioned shortcomings are mainly related to the mesh definition. Therefore, the new families of numerical methods independent of mesh are applied to the geotechnical problems -the slope stability analysis in this article. These new families of numerical methods are widely named as meshless methods.
In this article, the elasto-plastic analysis of soil slope structure using the Local Galerkin Radial Point Interpolation method (LGRPIM) is presented. LGRPIM is the "true" meshless method based on the local weak-formulation of the equilibrium equations [4]. The Radial point interpolation method (RPIM) is used to construct shape functions. The system of equations has been derived and the computer code developed. The validity has been investigated by solving an example and the result comparison is provided at the end of the article.

Meshless shape function construction
The crucial part of meshless methods is the construction of the shape function used to approximate the unknown field function [3]. These shape functions are locally supported, because only set of field nodes in a small local domain are used in the construction of shape function. In contrast with FEM, where the shape func-tion is based on nodes with defined relationship (elements), the meshless shape function is based on the group of nodes arbitrarily scattered within the support domain [5]. There are many techniques used to construct meshless shape functions. The point interpolation method (PIM) and radial point interpolation method (RPIM) are presented.

Point interpolation method -PIM.
The polynomials have been used as basis functions for interpolation in various numerical methods such as FEM. In FEM the interpolation process is based on connected elements without gaps or overlaps. In PIM, the interpolation is based on small portion of nodes around the desired point, forming the local support ( Fig. 1) [5]. Support domains of different points can overlap each other. The continuous function u(X) (displacement function) can be expressed around point X as (1) where p i (X) is polynomial function of X ϭ [x,y] T , n is count of nodes in support domain of point X and a i is suitable coefficient of the basis function. Unknown coefficients ai in Eq. 1 can be obtained by setting the u(X) to be nodal displacement at n nodes of the support domain. This can be written using matrix notation as: where U S is the vector of nodal displacement, and a is the vector of unknown coefficients a ϭ [a 1 a 2 a 3 … a n ] T (4) and P m is polynomial moment matrix (5) considering m ϭ n one can assume that the inversion of moment matrix P m Ϫ1 exists and the unique solution for coefficients a can be obtained as a ϭ P m Ϫ1 U S (6) It has to be noted that the coefficients a are constants even if the point of interest X changes, as long as the same set of nodes is used in the interpolation, because the P m is matrix of constants for given set of nodes. f where Φ(X) is the vector of PIM shape functions defined by (8) The derivatives of the shape functions can be easily obtained because PIM function is of polynomial form. The I-th derivative of PIM shape functions can be written The shape functions constructed by PIM have the Kronecker delta function property, which allows the simple imposition of essential boundary conditions as in conventional FEM [6].

Radial point interpolation method -RPIM
The PIM is accurate and easy to use. However, an inappropriate choice of polynomial basis or the improper position of nodes inside the support domain near X results in singular matrix P m [6]. Several strategies have been proposed to overcome this problem. Using radial basis functions (RBF) is one of the best solutions to guarantee the invertibility of P m [5]. There are various classes of functions such as multi-quadratic, Gaussian or Logarithmic that can be used as RBF. Multi-quadratic is one of the most popular radial function [5] and it is defined as where r i is the distance between the desired point (X) and the field node i(X i ) defined simply as 2D Euclidean distance Constants ε and q in Eq.10 are constants that depend on the type of problem. For the solid mechanics the suggested values are 1.42 and 0.98 for ε and q, respectively [3].
The RPIM interpolation augmented with polynomials can be written as where R i (X) is the radial basis function (RBF), n is the number of RBFs, p j (X) is polynomial basis function, m is number of polynomial basis function, a i and b j are interpolation coefficients. In order to determine a i and b j a support domain is formed for the point of interest at X, and n field nodes are included in the support domain.
Interpolation coefficients can be determined by enforcing Eq. 12 to be satisfied at these n nodes surrounding the point of interest X. This leads to n linear equations, one for each node. The equation system in matrix form can be expressed as where U S is the vector of displacement function values, the RBF moment matrix is (14) and the polynomial moment matri is (15) the vector of RBF interpolation coefficients is defined as And the vector of polynomial interpolation coefficients is defined as However, there are n ϩ m variables in Eq. 13. The additional m equations can be added using the following m constraint conditions.
(18) Combining Eq. 13 and Eq. 18 yields the following set of equations in the matrix form expressed as Because the matrix R is symmetric, the composed matrix G will be also symmetric. Solving the equation system (19), we obtain (22) Equation (12) can be subsequently rewritten as f Using equation (22) we can obtain (24) and finally the RPIM shape functions corresponding to the nodal displacements vector Φ(X) are obtained as (25)   (26) Then the displacement function approximation can be written using the RPIM shape functions and nodal displacements (27) and the I-th derivatives of displacement function u(X) are easily obtained as (28) Note that R 0 Ϫ1 usually exists for arbitrarily scattered nodes. In addition, the order of polynomial used in RPIM shape functions is relatively low. Therefore, there is no singularity problem in the RPIM as a small number of nodes is used in the local support domain [7].

Weak formulation of the equilibrium equations
The LGRPIM uses the local weighted residual formulation of the equilibrium equations rather than the global energy principle to create the discretized equation system. The compatibility of the shape functions in whole domain is not required as long as the field approximation is continuous at any point in the local quadrature domain. In other words the LGRPIM method only requires the local compatibility in the local quadrature domain. The RPIM shape function formulated in previous chapter satisfies all these requirements, in addition to its delta function property.

Local weighted residual formulation of the equilibrium equations.
The general local weighted residual form defined over local quadrature domain Ω q bounded by Γ q has following matrix form [5] K I u ϭ f where K I is the matrix called nodal stiffness matrix for the I-th field node, which is computed using following formula and the f I is a nodal force vector with contribution from the body forces applied in the model domain, and the tractions applied on the natural boundary The terms W and V in the previous definitions is the weight function matrix and weight function derivatives matrix respectively [6], used to smear the numerical residuals over the local quadrature. In equation (30) D represents the matrix of elastic constants, B is the strain-displacement matrix, n is matrix containing the components of the unit outward vector on the domain boundary, t is the traction defined on the domain boundary and b is the body force.
The weight function W definition used in this article follows the rectangular shape of the quadrature domain. Cubic spline type weight function was defined as follows [5] where d sx and d sy is the size of the local quadrature domain. (34) Assembling the nodal stiffness matrices into a global stiffness matrix based on node numbering, the global equation system is obtained. The global stiffness matrix is asymmetric, which makes the computation more expensive than conventional FEM [4]. For numerical example in this article the Gaussian elimination solver was used to solve global equation system [8]. Because the analysis of the soil models using the elasto-plastic models needs the nonlinear solution scheme, the modified Newton-Raphson scheme was adopted for solving the models. The nodal loads are generated from the self-weight of the slope body and redistributed over the model nodes. Once the plastic flow occurs, the stresses are computed using elasto-plastic matrix (Eq. 36) and the unbalanced forces are redistributed over adjacent nodes. This process is repeated until there is no unbalanced force or maximum iteration count was achieved.

Elasto-plastic constitutive equations
In elasto-plastic analysis, the incremental stress-strain relation is as follows [5]: where D ep is elasto-plastic matrix. In this article, the Mohr-Coulomb yield function was used as the most popular yield criterion used in soil mechanics.
where c and ϕ means the cohesion and internal friction angle respectively. The quantities σ m , σ, θ represents the stress invariants in the principal stress space. The plastic potential equation is similar to the yield function (Eq. 36), except that instead of internal friction angle the dilatation angle -ψ is used.

Numerical example of the slope stability computation
The LGRPIM model has been verified with the computation of the slope stability problem. Because there is no exact analytical solution of this kind of problem, the model of the embankment slope was calculated using conventional finite element code for geotechnical problems -PLAXIS and the Janbu's limit equilibrium method (LEM) implemented in GeoStudio 2003. Figure 2 shows the geometry and dimensions of the slope model and parameters for the Mohr-Coulomb elastoplastic model are presented in Table 1. The model represents simple homogenous slope body with "roller" essential boundary condition on left and right boundary (u r x ϭ u l x = 0) and fixed (u b x ϭ u b y ϭ 0) movement at the model bottom.
The result of the slope stability analysis is FS -factor of safety and the shape of the failure in the slope body. Factor of safety (FS) in numerical analysis is defined as the strength reduction factor (SRF) which represents the boundary between stable and unstable state of the model. In the case of the limit equilibrium methods (Janbu), the factor of safety is defined as the ratio between active (moving) and passive (resisting) forces along the slip surface.
To perform the integrations for the local weak-form, local quadrature domains are needed. The local quadrature domain is simple rectangular domain which is easy to handle. The size of the quadrature domains used in LGRPIM implementation is r qx ϭ ϭ r qy ϭ 3.75. The quadrature domain needs to be divided into 4 subdomains to ensure accuracy of the solution. Each partition contains 16 Gauss points [5]. The support domain used to construct RPIM shape functions has circular shape with radius r s ϭ 3.75. The cubic spline function (34) is used as the test function for the LGRPIM.
The Mohr-Coulomb model was used in PLAXIS computation as well as in the analysis performed using LGRPIM. The Janbu's limit equilibrium method is also using the Mohr-Coulomb strength criterion for stability analysis. The result of the LGRPIM analysis showing the slope failure (Fig. 2) and plastic flow vectors (Fig. 3) at the FS ϭ 1.63.
The parametric study with varying slope inclination angle (β - Fig. 2) was also performed to show how the factor of safety (FS) will develop. Parametric study was performed using LGRPIM, FEM and Janbu's limit equilibrium method. The result (Fig. 5) indicates close agreement of LGRPIM with FEM and Janbus's method for presented slope stability problem.
The relative difference of the computation results performed by LGRPIM against the FEM code are shown in the Table 2. The safety factor was developing similary with maximum overall difference 9.64%. Soil parameters used in the numerical study [1]

Conclusion
In this paper the application of the radial basis point interpolation method (RPIM) as an local weak formulated meshless method (LGRPIM) for the elasto-plastic analysis of slope structure. There is no need of mesh in the traditional sense and the shape functions are based on nodes. This method uses rectangular local quadratures with 4 sub-domains and 16 Gauss points per subdomain for numerical integration. The results of the analyses with the LGRPIM, FEM and LEM show very good agreement with each other. Hence the result of the parametric study presented in Figure  5 shows very good agreement with the results of the FEM and LEM result. This can be a positive outlook for the application of meshless methods in slope stability analysis. However, as the present study is the first one in this category, it is too soon to conclude generally, and more theoretical and experimental research is needed.

Acknowledgements
This contribution is the result of the project implementation: "Support of Research and Development for Centre of Excellence in Transport Engineering" (ITMS: 26220120031) supported by the Research & Development Operational Programme funded by the ERDF.
This contribution is the result of the project implementation: Scientific Grant Agency of the Slovak Republic (VEGA) No. 1-0789-12.

Fig. 5 Factor of safety variation with slope inclination angle
Factor of safety variation with inclination angle Table 2 Slope inclination [°] Factor of safety FS [-] LGRPIM