The 12th International Conference of

International Association for Computer Methods and Advances in Geomechanics (IACMAG)

1-6 October, 2008

Goa, India

A Stabilization Procedure for Soil-water Coupled Problems Using

the Mesh-free Method

T. Shibata

Dept. of Civil and Environmental Engineering, Matsue National College of Technology, Shimane, Japan

A. Murakami

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.

1 Introduction

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

64

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

as follows:

∫div ⋅ +∫( − )⋅ = 0 V Γ

t

u

S& δ v dV ρ v v δ v d S (1)

∫( ) ∫ ∫ ∫ ( ) Γ Γ

− + ⋅ − − − =

q h

tr hdV grad hdV q hdS pw pw hdS

V

w

V

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

pressure.

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 V w w w

p dV

tr hdV grad h dV q h dS p p h dS

q h

δ α

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

Order

Number

：Pressure

：Displacement

Figure 1. Rearrangement and allocation of background cells to cover the chang ing domain

65

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

30m

z = 0

z

x

= 0

∂

∂

z

pw

0

0

=

=

z

x

u

u

z = 30 = 0 w p

1m

q

0

0

=

=

x

z

u

u

= 0

∂

∂

x

pw

Figure 3. Effect of the dimensionless parameter

0 1 2 3 4

[×105]

0

10

20

30

=1.0× 10-7( )

Pore water pressure ( )

Depth ( )

k

Pa

m

cm/sec

without stabilization term

α =0.1

α =0.01

α =0.001

66

adopted.

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.

2.0

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

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

(a) Geometry and boundary conditions

λ = 0.11

κ = 0.04

M = 1.42

ν = 0.333

pw =0

10.0m

5.0 m

4.0m 2.0m 4.0m

e0 = 0.83

v0 = 1+ e0 =1.83

p0′ = 294 (kPa)

Loading plate

Figure 4. Sketch of the problem

(b) Initial collocation of the nodal points

- 50000

- 30000

- 10000

10000

30000

50000

70000

90000

[kPa]

(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

67

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

0

0.05

0.1

0.15

0.2

0.25

0.3

0 5

10

15

20

25

30

35

40

45

50

[kPa]

(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

0.1

0.2

Prandtl's solution

EFG solution Applied footing stress [

MPa]

Settlement [m]

68

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)

where dt

t = ∫

0

ε 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.

4 Conclusion

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.

5 References

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,

440-449.

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).

69

Wang JG, Karim MR, Lin PZ., 2007, Analysis of seabed instability using element free Galerkin method, Ocean Engineering,

34(2), 247-260.

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.

70

Home »
article »
A Stabilization Procedure for Soil-water Coupled Problems Using the Mesh-free Method

# A Stabilization Procedure for Soil-water Coupled Problems Using the Mesh-free Method

Subscribe to:
Post Comments (Atom)

## 0 comments:

## Post a Comment