The 12th International Conference of
International Association for Computer Methods and Advances in Geomechanics (IACMAG)
1-6 October, 2008
A Stabilization Procedure for Soil-water Coupled Problems Using
the Mesh-free Method
Dept. of Civil and Environmental Engineering, Matsue National College of Technology, Shimane, Japan
Graduate School of Environmental Science, Okayama University, Okayama, Japan
Keywords: Mesh-free method, soil-water coupled problem, stabilization procedure
ABSTRACT: The development of stability problems related to classical mixed methods has recently been
observed. In this study, a soil-water coupled boundary-value problem, one type of stability problem, is presented
using the Element-free Galerkin Method (EFG Method). In this soil-water coupled problem, anomalous behavior
appears in the pressure field unless stabilization techniques are used. The remedy to such numerical instability
has generally been to adopt a higher interpolation order for the displacements than for the pore pressure. As an
alternative, however, an added stabilization term is incorporated into the equilibrium equation. The advantages of
this stabilization procedure are as follows: (1) The interpolation order for the pore pressure is the same as that for
the displacements. Therefore, the interpolation functions in the pore pressure field do not reduce the accuracy of
the numerical results. (2) The stabilization term consists of first derivatives. The first derivatives of the
interpolation functions for the EFG Method are smooth, and therefore, the solutions for pore pressure are
accurate.In order to validate the above stabilization technique, some numerical results are given. It can be seen
from the results that a good convergence is obtained with this stabilization term.
The development of numerical computation technologies has enabled a variety of engineering problems to be
solved and has brought about remarkable progress in recent decades. Among the related findings, meshless
and/or mesh-free methods in particular have been applied to some problems for which the usual finite element
method is ineffective in dealing with significant mesh distortion brought about by large deformations, crack growth,
and moving discontinuities.
Various meshless and/or mesh-free methods have been used for geotechnical problems, instead of the finite
element method, to overcome the above-mentioned difficulties. Consolidation phenomena have been analyzed
by means of EFGM (Modaressi et al., 1998)(Nogami et al., 2004)(Murakami et al., 2005)(Wang et al., 2006), the
point/radial point interpolation method (PIM/RPIM) (Wang et al., 2001)(Wang et al., 2002), the local RPIM (Wang
et al., 2005), RKPM (Chen et al., 2001)(Zhang et al., 2005), and the natural neighbor method (Cai et al., 2005),
the transient response of saturated soil has been dealt with under cyclic loading by means of EFGM (Karim et al.,
2002)(Sato et al., 2006), wave-induced seabed response and instability have been examined by EFGM (Wang et
al., 2007) and RPIM (Wang et al., 2004), slip lines have been modeled by geological materials using EFGM
(Rabczuk et al., 2006), and a Bayesian inverse analysis has been carried out in conjunction with the meshless
local Petrov-Galerkin method (Sheu, 2006).
However, unless certain requirements are met in dealing with soil-water coupled problems for the finite element
computation, based on the coupled formulation becoming ill-conditioned, numerical instabilities will occur
(Chapelle et al., 1993). In order to overcome these weaknesses, several strategies have been proposed (Pastor
et al., 1999). For example, as a necessary condition for stability, the interpolation degree of the displacement
field must be higher than that of the pore pressure field. An alternative means of stabilization was also proposed
based on the Simo-Rifai enhanced strain method which even allows an equal order of interpolation degree for
both variables. However, these strategies are not directly applicable to meshless/mesh-free methods, because
all the nodal points simultaneously have the same degree of freedom for both the displacement field and the pore
pressure field, and no information between the element and the nodes can be utilized.
The purpose of this paper is to present a stabilization methodology for the mesh-free analysis of soil-water
coupled problems by incorporating the stabilizing term into the weak form. The following sections deal with
descriptions of the formulation, including the stabilization term. In Section 3, two applications of the strategy to
soil-water coupled problems are analyzed, one being the saturated soil column test appearing in Mira et al.(2003),
to demonstrate the effectiveness of the strategy, and the other being the foundation behavior under
displacement-controlled condition, for which the feasibility of theanalysis will be thoroughly discussed while
comparing the finite element solutions. Conclusions will follow in Section 4.
2 The Stabilization term
The governing equations for soil-water coupled problems with boundary conditions and initial conditions are given
∫div ⋅ +∫( − )⋅ = 0 V Γ
S& δ v dV ρ v v δ v d S (1)
∫( ) ∫ ∫ ∫ ( ) Γ Γ
− + ⋅ − − − =
tr hdV grad hdV q hdS pw pw hdS
Dδ v δ δ β δ 0 (2)
where S& t is the nominal stress rate, v is the velocity, v is the boundary value of the velocity, D is the stretching,
h is the total head, vw is the velocity of the pore water, q is the discharge per unit area with units of length per
time, β is Lagrange multiplier, pw is the pore water pressure, and pw is the boundary value of the pore water
As previously mentioned, mixed displacement-pressure formulations (e.g., finite element methods) produce
locking phenomena in the pressure field unless stabilization techniques are used. The remedy for such numerical
instability has generally been to adopt a higher interpolation order for the displacements than for the pore
pressure. As an alternative, however, an added stabilization term is incorporated.
It is shown in this study that the instability can be eliminated by adding the stabilization term to the continuity of
pore water (Equation (2)) which consists of the square of the pore water pressure of the first derivatives.
( ) ( )
( ) 0 + ′ 2 =
− + ⋅ − − −
∫ ∫ ∫ ∫ Γ Γ
V V w w w
tr hdV grad h dV q h dS p p h dS
Dδ v δ δ β δ (3)
where α is the stability parameter, p′ is the differentiate p with respect to z , and z is the vertical axis. In this
study, the soil column test created by Zienkiewicz et al. is performed in order to examine the numerical stability of
this procedure, and the values for the pore water pressure are illustrated along the vertical axis of the soil column.
In order to smooth the values for the pore water pressure along the vertical axis with the continuity
condition of the pore water, the stabilization term is differentiated with respect to z .
The two advantages of this stabilization procedure are as follows: (1) The interpolation order for the pore
pressure is the same as that for the displacements, namely, a lower interpolation order is not adopted for the pore
Table 1. Order of the interpolation function in FEM
0-1 1-1 1-2 2-2
Pressure 0 1 1 2
Displacement 1 1 2 2
Figure 1. Rearrangement and allocation of background cells to cover the chang ing domain
pressure. Therefore, the interpolation functions in the pore pressure field do not reduce the accuracy of the
numerical results. (2) Table 1 summarizes the order of the interpolation functions for the pressure field and for
the displacement field in FEM. The first derivatives of the interpolation functions in the pore pressure field are the
zero order or the first order, as shown Table 1. Therefore, accuracy in the numerical results cannot be obtained.
With EFGM, however, the interpolation functions are derived by the MLS approximant and a linearbased
polynomial is used, in other words, the resultant interpolation functions are smooth.
The weak forms are integrated using the MLS interpolation functions in space and employing the explicit time
scheme in time. In the manipulation of the stiffness matrix, a numerical integration is performed on the
background cells. The background cells, which intersect or make contact with the boundary, are divided into four
finer cells. The cells are then rearranged during the computation, according to changes in the domain, as
shown in Figure 1.
3 Numerical Examples
3.1 The saturated column soil test
For soil-water coupled problems, numerical instabilities are often encountered at the initial stage under undrained
conditions unless stabilization techniques are performed. In this chapter, the numerical stability of an
EFG computation is examined using the proposed stabilization term described in the last chapter.
Let’s consider the 1D problem analyzed by Zienkiewicz et al.(1986) in which the saturated soil column test in
Figure 2 uses the material parameters listed in Table 2. Here, a sand permeability of 1.0 × 10-7 cm/s is employed.
A model discretized by 62 nodal points is adopted. Linear-based polynomials are applied for the interpolation
functions of the EFG Method. The functions have the same order for both the displacements and the pore
pressure. The weight function is a quartic spline type of weight function and its radius of support is 1.0. Here, 5 ×
5 Gaussian points are used. The background cells are 1.0m within the domain and 0.5m outside of the domain.
The scale factor, which is defined as the magnification of the support diameter to the side length of the square
background cell, is 1.5. This study employs penalty methods to apply the boundary conditions, and the value of
the penalty factor is 1.0 × 106. In order to solve the stiffness equation, the forward difference is
Table 2. Material parameters
Compression index λ 0.11
Swelling index κ 0.04
Critical state parameter M 1.42
Poisson’s ratio ν 0.333
Initial void ratio e0 0.83
Initial volume ratio v0 =1+ e0 1.83
Initial consolidation stress (kPa) p0′ 294
Figure 2. Saturated column test
z = 0
z = 30 = 0 w p
Figure 3. Effect of the dimensionless parameter
0 1 2 3 4
=1.0× 10-7( )
Pore water pressure ( )
Depth ( )
without stabilization term
Figure 3 shows the effect of the proposed stabilization term and dimensionless parameter α on the numerical
profile of the pore pressure just after loading under undrained conditions, for which a time difference of Δt = 0.01
and a value of permeability of k = 1.0×10-7 cm/s are adopted. There are 62 nodal points and the weight function
is a quartic spline. From Figure 3, it can be concluded that the numerical solution has been improved in the case
where the value of the dimensionless parameter is 0.01. In subsequent examinations, a permeability of
1.0 × 10-7cm/s and a time difference of 0.01 are adopted.
Here, the numerical instability is caused from the large difference of the value in the stiffness matrix in between
displacements and pore water pressures. Generally, the values of stiffness matrix in pore water pressure are
smaller than in displacements. However, if the lower order interpolation function in pore water pressures is
adopted, the values of stiffness matrix in pore water pressures become large, so that the difference of the values
in stiffness matrix becomes small. Similarly, if the stabilization term is considered, the values of stiffness matrix
(a) Geometry and boundary conditions
λ = 0.11
κ = 0.04
M = 1.42
ν = 0.333
4.0m 2.0m 4.0m
e0 = 0.83
v0 = 1+ e0 =1.83
p0′ = 294 (kPa)
Figure 4. Sketch of the problem
(b) Initial collocation of the nodal points
(a) Contours of the pore pressure (b) Contours of the normalized strain
Figure 5. Numerical results without stabilization term
(c) Displacement distributions (d) Collocation of the nodal points
in pore water pressure also become large. Therefore, the stabilization term quells the anomalous pore pressure
behavior because of the small difference.
3.2 Foundation Problem Subjected to Strip Loading
In order to examine the numerical availability of the stabilization procedure to EFGM, a 2D soil-water coupled
problem is solved in relation to a soft soil foundation. In this problem, the p& in Equation (1) is the differentiate p
with respect to x and z , with x being the horizontal axis and z being the vertical axis. The geometry and
the boundary conditions are given in Figure 4(a), while the initial collocation of the nodal points is shown in Figure
4(b). Vertical displacements are applied at the top face to simulate loading on the foundation, while the
transverse direction is restrained. Thus, the boundary conditions under the loading surface are given by the
settlements, and the other locations on the upper boundary are stress free. The bottom boundary is fixed in all
directions, while the side boundaries are fixed in only the horizontal direction, such that vertical displacement is
allowed. Hydrostatic pressure is used here for the initial conditions in the pore pressure field. The model is
generated with 1,326 nodal points, in other words, 26 vertical nodal points and 51 horizontal nodal points. The
size of background cells is 0.2m, and the value of the dimensionless parameter α is same as the previous
analysis. Figures 5 and 6 present the numerical results without and with the stabilization term, respectively, in
which (a), (b), (c) and (d) show the contours of the pore water pressure, the contours of the normalized strain, the
displacement distributions and collocation of the nodal points, respectively. In these results, the settlements
(a) Contours of the pore pressure (b) Contours of the normalized strain
(c) Displacement distributions
Figure 6. Numerical results with stabilization term
(d) Collocation of the nodal points
Figure 7. Comparison EFG solution with Prandtl’s solution
00 0.1 0.2 0.3 0.4
EFG solution Applied footing stress [
under the loading surface are 0.04m in Figures 5 and 0.4m in Figures 6, respectively. The normalized strain
measure is given as
ε = tr(εε T ) (2)
t = ∫
ε D (Yatomi et al., 1989).
Spurious oscillations arise in the pore pressure field, as can be seen in Figure 5(a). In contrast to the results
using the stabilization procedure, high values are obtained for the pore pressure and the normalized strain in
Figures 5(a) and (b), respectively. In particular, anomalous behavior appears on left side in these figures.
Moreover, the directions of the displacement vectors are disorderly in Figure 5(c) because of the interaction
between the pore pressure field and the displacement field. If stabilization techniques are used, however, no
oscillations appear in the solution, as observed in Figure 6. Figure 6(a) accurately displays the rise in the pore
water pressure below the loading surface. In Figure 6(b), the prominently localized zones of the normalized strain
occur just beneath the edge of the loading surface. The deformed pattern in Figure 6(c) is very similar to the
classical slip line solution obtained by Prandtl. The shear bands are recognized as the localized deformation.
Figure 7 compares the EFG solution with Prandtl’s solution. The Prandtl’s solution q f is expressed as
q f =5.14 cu (3)
where cu is the undrained shear strength. From this figure, it is revealed that the EFG solution approaches
Prandtl’s solution. This result shows that the EFG method with stabilizing procedure is very efficient for solving
problems of computational geomechanics.
In this paper, we have proposed a stabilization method for soil-water coupled problems using the Element-free
Galerkin Method. A stabilization term has been presented by the addition of a continuity condition. The proposed
stabilization procedure has the following two characteristics: (1) The interpolation order for the pore pressure field
is the same as that for the displacements, in other words, a lower interpolation order for the pore pressure is not
adopted. (2) The stabilization term consists of first derivatives in which the interpolation functions are smooth
because of the MLS approximation. The saturated column test and the foundation loading problem have been
solved using the st a b i l izat ion procedure. Numerical examples have shown that the
stabilization method can indeed quell the anomalous pore pressure behavior.
Modaressi H, Aubert P., 1998, Element-free Galerkin method for deforming multiphase porous media, International Journal for
Numerical Methods Engineering, 42, 313-340.
Nogami T, Wang W, Wang JG., 2004, Numerical method for consolidation analysis of lumpy clay fillings with meshless method,
Soils and Foundations, 44(1), 125-142.
Murakami A, Setsuyasu T, Arimoto S., 2005, Mesh-free method for soil-water coupled problem within finite strain and its
numerical validity, Soils and Foundations, 45(2).
Wang ZL, Li YC., 2006, Analysis of factors influencing the solution of the consolidation problem by using an element-free
Galerkin method, Computers & Geosciences, 32, 624-631.
Wang JG, Liu GR, Wu YG., 2001, A point interpolation method for simulating dissipation process of consolidation, Computer
Methods in Applied Mechanics and Engineering, 190, 5907-5922.
Wang JG, Liu GR, Lin P., 2002, Numerical analysis of Biot’s consolidation process by radial point interpolation method,
International Journal of Solids and Structures, 39, 1557-1573.
Wang JG, Yan L, Liu GR., 2005, A local radial point interpolation method for dissipation process of excess pore water pressure,
International Journal of Numerical Methods for Heat & Fluid Flow, 15(6), 567-587.
Chen JS, Wu CT, Chi L, Huck F., 2001, A meshfree method for geotechnical materials, Journal of Engineering Mechanics, 127,
Zhang JF, Zhang WP, Zheng Y., 2005, A meshfree method and its applications to elasto-plastic problems, Journal of Zhejiang
University SCIENCE, 6A(2), 148-154.
Cai YC, Zhu HH, Xia CC., 2005, Meshless natural neighbor method and its application to complex geotechnical engineering,
Yanshilixue Yu Gongcheng Xuebao/Chinese Journal of Rock Mechanics and Engineering, 24(11), 1888-1894.
Karim MR, Nogami T, Wang JG., 2002, Analysis of transient response of saturated porous elastic soil under cyclic loading
using element-free Galerkin method, International Journal of Solids and Structures, 39, 6011-6033.
Sato T, Matsumaru T., 2006, Numerical simulation of liquefaction and flow process using mesh free method, Journal of
Geotechnical Engineering, 813, 89-101 (in Japanese).
Wang JG, Karim MR, Lin PZ., 2007, Analysis of seabed instability using element free Galerkin method, Ocean Engineering,
Wang JG, Nogami T, Dasari GR, Lin PZ., 2004, A weak coupling algorithm for seabed-wave interaction analysis, Computer
Methods in Applied Mechanics and Engineering, 193, 3935-3956.
Rabczuk T, Areias PMA., 2006, A new approach for modeling slip lines in geological materials with cohesive models,
International Journal for Numerical and Analytical Methods in Geomechanics, 30, 1159-1172.
Sheu GY., 2006, Direct back analysis by the meshless local Petrov-Galerkin method and Bayesian statistics, International
Journal for Numerical and Analytical Methods in Geomechanics, 30(8), 823-842.
Chapelle D, Bathe KJ., 1993, The inf-sup test, Computers & Structures, 47, 537-545.
Pastor M, Li T, Liu X, Zienkiewicz OC., 1999, Stabilized low-order finite elements for failure and localization problems in
undrained soils and foundations, Computer Methods in Applied Mechanics and Engineering, 174, 219-234.
Mira P, Pastor M, Li T, Liu X., 2003, A new stabilized enhanced strain element with equal order of interpolation for soil
consolidation problems, Computer Methods in Applied Mechanics and Engineering, 192, 4257-4277.
Zienkiewicz OC, Qu S, Taylor RL, Nakazawa S., 1986, The patch test for mixed formulations, International Journal for
Numerical Methods Engineering, 23, 1873-1883.
Yatomi C, Yamashita A, Iizuka A, Sano I., 1989, Shear bands formation numerically simulated by a non-coaxial Cam-Clay
model, Soils and Foundations, 29, 1-13.
The 12th International Conference of