Journal of Natural Disaster Science, Volume 26, Number 1, 2004, pp37-48

Adaptive Mesh Refinement and Error Estimate for 3-D Seismic
Analysis of Liquefiable Soil Considering Large Deformation
Xiaowei TANG and Tadanobu SATO
Disaster Prevention Research Institute, Kyoto University, Uji Kyoto 611-0011, Japan(Received for 28. Sep., 2004 and in revised from 27 Dec., 2004)

ABSTRACTThe focus of this paper is on 3-D adaptive analysis method of liquefiable soil improvement of the approximate
quality of nonlinear numerical simulation of the liquefaction process with large deformation. This adaptive technique was applied to the 3-D non-linear FE analysis of liquefiable soil considering large deformation, including
the liquefaction phenomenon. The cyclic elasto-plastic model and updated Lagrangian formulation were adopted
in the three-dimensional FE analysis. The fission procedure belonging to the h-refinement indicated by the error
estimator, a posteriori error estimate procedure depending on L2-norm of strain and superconvergent patch recovery method, was used. The convergence of this error estimate scheme and the effectiveness of h-adaptive meshrefinement were shown by simple examples. The adaptive FE method was applied to the three-dimensional practical seismic analysis of an embankment constructed on liquefiable soil, and its efficacy demonstrated in detail.
Application of the adaptive technique to the nonlinear analysis of saturated soil, including the liquefaction
process, is a valuable step in the three-dimensional seismic analysis of liquefaction.

1. INTRODUCTION
Such soil structures as river dikes, high way embankments
and earth dams frequently have been damaged during major earthquakes, often due to liquefaction of embankments and foundationsoils. In most cases, unacceptable, large, permanent deformation
has occurred due to liquefaction of the supporting, loose, cohesionless foundation soil. The seismic design of soil structures in theseliquefiable soils poses very difficult problems for analysis anddesign. The possibility of liquefaction necessitates further investigation not only experimentally but by numerical analysis.
In computational analysis, saturated soil is described as a typeof porous media. Biot’s two-phase mixture theory leads to theapplication of an effective nonlinear finite element method throughderivation of the governing equations expressed by u-p formulation; for example, a finite element and finite difference (FE-FD)
coupled scheme (Akai and Tamura, 1978). That scheme used toanalyze liquefaction with infinitesimal strain condition, reduces thetotal degree of freedom of the coupled equations. It avoids shear
locking under the undrained condition and involves less computational effort than the generalized finite element method. An appropriate constitutive model also is important in the nonlinear analysis
of saturated soil. Many effective constitutive models have beendeveloped in recent decades. The cyclic elasto-plastic model (Oka,
1992) is one of most effective. It was proposed as a numericalmethod to simulate the liquefaction of saturated soils and is basedon the non-linear kinematic hardening rule. Because deformation
caused by liquefaction usually is very large, the finite deformation
theory was introduced to the finite element method. Further development of the FE-FD coupled scheme was achieved in two- andthree- dimensional analysis by use of the updated Lagrangian for-mulation (Sato and Di, 2001; Tang and Sato, 2002). The cyclic
elasto-plastic model was adopted as the constitutive model of liquefiable soil in FE-FD coupled scheme used in the study reported
here. The material nonlinearity of saturated sand is simulated
effectively by the constitutive model. The updated Lagrangian formulation, which belongs to the finite deformation theory, also wasapplied to the three-dimensional FE-FD coupled scheme to deal
with the geometrical nonlinearity of liquefiable soil caused bylarge deformation.
As a numerical approximation method, errors are inevitable inanalytical results obtained by the finite element method. Solutions
do not always produce the desired accuracy, sometimes presenting
serious analysis problems, especially when considering large
deformation. Such error is caused by discretization in the FEM
process. In practical problems of liquefied soil flow caused byearthquakes, soil deformation is not uniform and where large
deformation is developed can not be predicted. With elements thathave large deformation, large numerical error develops. Evidently,
uniformly reducing element size during discretization minimizes
error, but the number of nodes and elements are increased at thesame time, producing a heavy calculation burden. In fact, a region
with small deformation does not require a finer mesh. It is better
to use a fine mesh in a region with large error and a normal orcoarse mesh in a low-error region. A method called the adaptive
technique or adaptive mesh refinement has been developed and utilized to reduce discretization error. It has been used successfullyin many fields, including solid and fluid mechanics, to linear andnonlinear problems and to solve static and transient behavior of
two-dimensional continua. Extension to three-dimensional,
dynamic applications, however, is in its infancy. It is still a major
challenge to apply adaptive mesh refinement to the seismic analy-

KEY WORDS: liquefiable soil, large deformation, updated Lagrangian formulation, error estimate, adaptive mesh refinement

38

X. TANG AND T. SATO

sis of liquefaction.
A posteriori error estimate based on the L2 norm of strain error
was adopted in the adaptive analysis used in our study. It effectively
estimates the element error after each calculation step in the nonlinear
FE analyses of soil. The superconvergent patch recovery technique(Zienkiewicz and Zhu, 1992) was used in the error estimates. The
convergence of this error estimate method was tested in two loadcases for a simple example with different quality meshes.
Convergence ratios obtained by our method agree with the predictedvalues. The fission procedure of h-refinement, indicated by the error
estimates of the elements in adaptive analysis, was used. The approximation was refined successively so as to satisfy the predetermined
standard of accuracy. The efficacy of this method was confirmed
for another simple example of consolidation of saturated soil.
Lastly, a practical example of seismic analysis of an embankment, based on the finite deformation theory, was determined withthe adaptive FE method. Adaptive analysis results for mesh deformation, error distribution, average error history, and displacement
responses were obtained and compared with results calculated
without the adaptive procedure. The efficacy of three-dimensional
adaptive analysis also was demonstrated by means of practicalexample.
2.

GOVERNING EQUATIONS USING THE UPDATEDLAGRANGIAN FORMULATION

2.1 Constitutive equation
In this research, the effective cyclic elasto-plastic constitutive
model was used to describe the non-linear behavior of saturated
soils. The deformation tensor rate is suitable for application to a
constitutive law because it vanishes when the body has rigid-body
motion. The stress rate used for a constitutive law therefore must
be invariable with respect to rigid-body rotation. The material
derivative ds ij / dt does not, in general, satisfy invariance due torigid-body rotation. To extend this model to a large deformation
problem it is necessary to use an invariant stress rate with respect
to rigid-body rotation for the constitutive relations. The Jaumanneffective stress rate, an objective measure of the stress rate, was;s˙ ij¢ = s˙ ij¢ - s ik¢ w jk - s kj¢ w ik
J(1)

where s ij¢ is the effective Cauchy stress, s˙ ij¢ the effective Cauchy
stress rate, and w ij the antisymmetric spin tensor.
Generally, large deformation cannot be expressed linearly interms of displacements, because the elastic and plastic parts of thedeformation rate are not summable. If each time step in the updated Lagrangian analysis is small enough, the total Lagrangian strain
rate Ėij can be decomposed
into the elastic component, Ėije , and.p
the plastic component, Eij.
The linear relationship between the objective stress and deformation rates (i.e., the effective stress-strain law) iss˙ = s˙ ij¢ - ˙pd ij = D l - ˙pd ij
J
ij

J

ep
ijkl kl(2)

epwhere ṗ is the rate of pore pressure, and Dijkl
the elasto-plastic
tensor of the cyclic elasto-plastic constitutive model.

2.2 Equilibrium equation
According to Biot’s two-phase mixture theory, the local equilibrium equation of motion for the total saturated porous media is∂s ji+ rbi - (1 - n)r s v˙ is - nr f v˙ if = 0
∂x j

(3)

where vis is the velocity of the soil skeleton, vif the velocity of thefluid phase, n the porosity, r the mixture density, r s the solid
phase density and r f the liquid phase density.
.
The average seepage velocity w i is(ẇi = n vif - vis)(4)

Substituting Eq.(4), Eq.(3) is rewritten
˙
∂s ji
Ê ∂w
ˆ
˙ kw
˙ i ,k ˜ = 0+ rbi - r˙˙ui - r f Á i + w
Ë ∂t
¯
∂x j

(5)

where s ij is the Cauchy total stress in the combined solid and fluid
mixture, and bi the body force acceleration.
For the pore fluid, the equation of motion is-∂(np)+ nr f bi + Ri - nr f v˙ if = 0
∂xi

(6)

where p is the pore pressure (taken as positive when compressive),and Ri the viscous drag force of the soil skeleton acting on the fluid
phase, which, according to Darcy’s law, can be written
Ri = - n

gf
˙
w
k i

(7)

where k is the Darcy permeability coefficient in isotropic condition, and g f the unit fluid weight.
Substituting Eqs.(4) and (7) in Eq.(6) gives
Ê
∂(np)
∂w˙
∂w˙ ˆ
- nr f bi + ng f kij-1w˙ i + r f Á nu˙˙i + ( i + w˙ k i )˜ = 0
∂xi
∂t
∂xk ¯
Ë

(8)

We assume that the acceleration of the fluid phase relative tothe soil skeleton can be neglected. The u-p approximation(Zienkiewicz et al., 1982 and 1984) is valid for a low frequency
dynamic analysis problem. The equations of motion for the total
mixture given by Eq.(5) can be simplified to∂s ij+ rbi - ru̇˙i = 0
∂x j

(9)

and the pore fluid motion equation (8) rewritten as∂(np)
˙ i + nr f ˙˙ui = 0- nr f bi + ng f k -1w
∂xi

(10)

Defining the excess pore pressure pE , by∂(np - npE )= nr f bi
∂xi

(11)

Eq.(10) becomes∂(npE )
˙ i + nr f ˙˙ui = 0+ ng f k -1w
∂xi

(12)

From Eq.(12), we obtain
˙i = wˆk Ê 1 ∂(npE )+ r f ˙˙ui ˜
g f ÁË n ∂xi¯(13)

2.3 Continuity equation
According to the law of mass conservation, the continuity

ADAPTIVE MESH REFINEMENT AND ERROR ESTIMATE FOR 3-D SEISMIC ANALYSIS OF LIQUEFIABLE SOIL CONSIDERING LARGE DEFORMATION

define a coupled set of equations for saturated soils, in which ui
and pE are the only unknown variables.

equation of the soil skeleton in the local form is() + ∂( r (1 - n)v ) = 0∂ r (1 - n)
s

∂t

s

s
i

(14)

∂xi

Similarly, for the fluid it is( ) + ∂(nr v ) = 0∂ nr f

f

∂t

∂xi

f

i

(15)

After manipulation (Oka, et al., 1994), Eqs.(14) and (15) give
˙i
∂w
r˙
r˙
+ lii + n f + (1 - n) s = 0
∂xi
r
r
f

s

(16)

where l i j is the symmetric rate of the deformation tensor. r s isconstant and ṙ s is zero because soil skeletons are assumed to be.
incompressible. r f is the material derivative of the fluid phase.
density related to the material derivative, pE , of the excess porepressure;

2.4 u-p equations in spatial discrete form
The equilibrium equation is satisfied at the end of each time
step, t + Dt . In the updated Lagrangian method, relevant quantities, such as stress and strain, are correlated with the reference configuration at time t, and the weak formulation of Eq.(9) istt + Dt)S˙ ij dt dE˙ ij d tV + Ú t v t s ij dE˙ ij d tV = t + Dt F

(24)

where s ij is the Cauchy stresstensor at time t, Ṡij the second.
Piola-Kirchhoff stress rate, E ij the Lagrangian strain rate tensor,
and t + Dt F the external virtual work done by the applied bodyforces and tractions. Let t + Dt b i be the force acceleration per unit
volume and t + Dt t i the traction, then t + Dt F istt + DtF = Ú t + Dt A t + Dt ti dvi d t + DtA + Ú t + Dt V r t + Dtbi dvi d t + DtV(17)
Ú t V r ˙lii d V - Ú t V
fwhere K f is the bulk modulus of the fluid phase.
Substituting Eq.(17) in Eq.(16) gives(25)

Ê ∂ 2 ( pE ) ˆ t
gf
lii d tV + Ú t V ÁdV
2 ˜kË ∂xi ¯
i

(26)

(18)

.
From Eq.(13) the sum of the partial derivatives of w i withrespect to the coordinate xi isÊ∂w˙ i
1 Ê ∂ 2 p ˆ ˆ k 1 ∂pE
1 ∂pE ∂r fkÊ= - Á l˙ii + f Á 2 E ˜ ˜ - Á f2r Ë ∂ xi ¯ i ¯ g ÁË nr ∂xi∂xi
gË
r f ∂xi ∂xi( )ˆp Ê ∂ 2n ˆpE ∂n ∂r f
p ∂n ∂n+ Ef Á 2 ˜ - 2Ef˜2nr Ë ∂ xi ¯ n r f ∂xi ∂xi n r ∂xi ∂xi ˜¯

( )
)

(19)

(20)

The porosity distribution is assumed to be smooth enough inthe soils, therefore satisfying
Ê ∂ 2 (ln n) ˆÁ ∂x 2 ˜ = 0Ë¯ii

By means of the FE-FE coupled method (Oka et al, 1991 and1994), the terms associated with excess pore pressure in Eq.(26)
are discretized by the finite difference method in the space domain
in a grid pattern that is the same as that of the mesh for the finite
element method. By the Gauss theorem, the third term on the left
side of Eq.(26) becomesÚ t V — pE d V = Ú t A
2

where g is the acceleration caused by gravity.
From Eq.(17),
f

t

ng f
˙pE d tV = 0
-Út V
kK f

˙i
∂w
n
+ lii + f P˙ E = 0
∂xi
K

(

(

Ú t v r˙˙ui dvi d V + Ú t v Út

Integrating Eq.(23) for element volume gives

˙pr˙ f = Ef r fK∂ ln r
∂pE
=Kf
∂xi
∂xi

39

(21)

t

∂pE t
ni d A∂xi(27)

where tA is the element surface and ni the normal direction vector
of tA.
The value of the partial derivative ∂pE / ∂xi on the common
boundary of two elements is approximated by a difference expression in terms of excess pore pressure values at the gravity centersof the element and adjacent elements. It is calculated by
Út A

6
6
∂pE t
ni d A ª pE Â n Ai ◊ n Ei Ai + Â pEi n Ai ◊ n Ei
i Ai
i
=
i
=
1
1
∂xi

(

)

(

)

where pE is the excess pore pressure value at the gravity center ofthe element; pE i that of the adjacent element; A i the area of the( )If the gradients of ln(n) and ln r
are so small that theirquadratic terms in the above expressions can be ignored, considering Eq.(21), then Eq.(19) is expressed byf

˙i∂w
1 Ê ∂2 p ˆ ˆkÊ= - Á ˙lii + f Á 2 E ˜ ˜
r Ë ∂ xi ¯ i ¯
∂xi
gË

(22)

Finally, substituting Eq.(22) in Eq.(18), the final form of thecontinuity equation isk˙k Ê ∂2 p ˆnlii - lii + f Á 2 E ˜ - f ˙pE = 0
g Ë ∂ xi ¯ i K
g(23)

Clearly, Eqs.(14) and (23), together with the constitutive law,(28)

Fig. 1

Dissipation pore water between elements

40

X. TANG AND T. SATO

joint surface between the element and the adjacent one i; n Ai thenormal direction vector of A i; n Ei the normal direction vector of
dissipation from the element to the adjacent one i; n Ei i the normal
direction vector of dissipation from the adjacent element i to theelement (Fig.1).

norm for the entire solution domain. Its definition is
ee

h=

+ eh

2

nel
2
eh = ÊÂ eh i ˆ
Ë i =1
¯

1/ 2

2

ee

(32)

where
3. A POSTERIORI ERROR ESTIMATE
Error estimate is the first and most important procedure in
adaptive FE analysis. It indicates the next step, the mesh refinement. A posteriori error estimate procedure based on evaluating
the L 2-norm of variable error and the superconvergent patch
recovery (SPR) procedure (Zienkiewicz and Zhu, 1992) was usedin this study. The formulations were derived for linear hexahedral
elements. A simple example of compression of saturated soil was
analyzed in order to discuss the convergence of this error estimator.
3.1 Definition and evaluation of error
Error is defined as the difference between the exact solution
and the finite element approximation value. Variables consideredin the error estimate are displacement, strain, and stress. For
example, an error in strain (Kelly, et al., 1983) is described as
ee = e * -e h(29)where e * represents exact solutions and e the values of finite
element approximation.
h
A one-dimensional linear approximation of strain e and e *is shown in Fig.2 to explain the approximate value of FEM andexact solution.
The direct definitions of error given by Eq. (29) are not convenient for use in the process of error estimation. Usually, scalarnorms, such as energy or L 2 norms, are used to measure error.
Scalar measure corresponds to the square root of the quadraticerror. In this study, the L 2 norm was used to measure error, as it
ise associated with errors of any quantity. For the strain in element
i, the L 2 norm of the error, ee i , is defined ash(ee i = ÚVi {ee }i {ee }i dv / Vi
T

)

1/ 2

(30)

The error of entire solution domain ee is calculated by summing
the errors of all the elements.
ee = Ê Â ee i ˆË i =1
¯
nel

2

1/ 2

The local error indicator for the i-th element is defined similarly.
Clearly, the energy and L 2 norms are related to the strain energy,
ee i

hi =

2

ee

+ eh

2

(34)
/ nel

3.2 Smoothing algorithm with superconvergent patch recovery
In the error estimate process, the relatively accurate values
rather than the exact solution are used to calculate errors because
the exact solution is not easy or impossible to obtain. The superconvergent patch recovery technique (Zienkiewicz and Zhu, 1992)
was used. It is a single, continuous polynomial expansion of thefunction describing the derivatives and used in an element patch
surrounding nodes at which recovery is desired. This expansion
can be made to fit the superconvergent points locally in a leastsquares manner or can simply be an L 2 projection of the consistent
finite element derivatives. In this technique, nodal values are
assumed to belong to a polynomial expansion of the same order asthat presented in the shape function. This polynomial expansion isvalid over an element patch surrounding the particular assembly
node considered. Such a “patch” represents a union of elements
containing the vertex node. A typical patch for three-dimensionalhexahedral elements is shown in Fig.3. This polynomial expansion
is used to obtain the nodal values of strain;e * = {P}{a}

T

(35)

where {P} contains the appropriate polynomial terms, and {a} is aset of unknown parameters. For the eight nodes the hexahedral
element in three-dimensions{P} = {1, x , y, z, xy, yz, zx , xyz}(36){a} = {a , a , a ,L, a }

(37)

1

(31)

(33)

2

3

8

where nel is the total number of elements.
Generally, a relative error is used in the practical adaptivity
process because it is more easily interpreted. The relative error of
the solution is estimated by use of the error norm and strain L 2-

Fig. 3Fig. 2

Approximate values and exact solution

Three-dimensional superconvergent patch recovery
for hexahedra elements

ADAPTIVE MESH REFINEMENT AND ERROR ESTIMATE FOR 3-D SEISMIC ANALYSIS OF LIQUEFIABLE SOIL CONSIDERING LARGE DEFORMATIONThe unknown parameter vector {a} of the expansion in Eq.
(35) is determined by ensuring a least squares fit of this expansion
to the existing set of superconvergent, or at least high accuracy,
sampling points in the patch considered if such points are available.

(

)

2

(

F (a ) = Â e h ( xi , yi ) - e * ( xi , yi ) = Â e h ( xi , yi ) - {P}{a}
n

i =1

n

i =1

T

)

2

(38)

where (xi , yi ) are co-ordinates of a group of sampling points and n
is the total number of sampling points. The minimization condition of F(a) implies that {a} satisfies

{

}

{

}

Â PT ( xi , yi )P( xi , yi ) {a} = Â e h ( xi , yi ) PT ( xi , yi )
n

i =1

T

n

i =1

(39)

This can be solved in matrix form as
a = A -1b

(40)

where
A = Â P T ( xi , yi )P( xi , yi ) and b = Â P T ( xi , yi )e h ( xi , yi ) ;
n

n

i =1

i =1

P( xi , yi ) = {P( xi , yi )}The discrete expression form for three-dimension hexahedra
elements is:

1×3 m2 beside a side. Three meshes with 27, 216, and 729 elements were analyzed. For the boundary constraints of these meshes, the bottom nodes were fixed, and the nodes on side surfaces
allowed to slide on the surfaces.
Figs.5-7 show the 27-, 216- and 729-element meshes and therespective error contour results for Load 1. Figs.8-10 show the 27-,
216- and 729-element meshes and their respective error contour
results for Load 2. The relative error of an element was calculated
at its center and the contour obtained by interpolation from the relative error values at the nodes. By comparing the average relative
error values of the meshes with different quality levels, the relative
error of the coarse mesh is readily shown to be larger than that of
the fine mesh. In addition to the point with the peak value of relative error, the estimation of relative error is affected by the size ofthe element. When the size is small, the distance between the center of the element and point with the peak value is short, and therelative error of this element approaches the peak value. In contrast, the relative error of this element in the coarse mesh does notaccurately approach the peak value. This is why the relative error
value of a fine mesh is larger than that of a coarse mesh at somepoints. This phenomenon occurs in comparisons between mesheswith small differences in mesh quality levels but evidently does
not affect the average relative value of the entire mesh.
In Fig.11, the lines show the relationship between nodes number and the average relative error value, eta_m, calculated by Eq.

È Â P (x ¢,h¢,z ¢)P (x ¢,h¢,z ¢) L Â P (x ¢,h¢,z ¢)P (x ¢,h¢,z ¢)˘
1
i
i
i
8
i
i
i
˙
Íi 1 i i i 1 i i i
i
˙
Í
M
O
M
˙
Í
ÍÎÂi P8 (xi¢,hi¢,z i¢)P1 (xi¢,hi¢,z i¢) L Âi P8 (xi¢,hi¢,z i¢)P8 (xi¢,hi¢,z i¢)˙˚hÏa1 ¸ Ï Â P1 (xi¢,hi¢,z i¢)e i ¸
Ô
ÔÔ
i
Ô Ô Ô
M
ÌM ˝=Ì˝Ôa Ô Ô Â P (x ¢,h¢,z ¢)e h ÔÓ 8 ˛ ÓÔ i 8 i i i i ˛Ô(41)

where (xi¢,hi¢,z i¢) are the normalized coordinates of a group of
sampling points in three-dimensions.

Fig. 5

Mesh, load, and error contour in Load 1 (27 elements)

Fig. 6

Mesh, load, and error contour in Load 1 (216 elements)

Fig .7

Mesh, load, and error contour in Load 1 (729 elements)

3.3 Convergence of error estimator in three-dimension
A simple example was calculated to test this error estimator inthree-dimensions. A 3×3×3m3 cubic block of saturated soil in acontainer had a vertical load applied to the top surface locally(Fig.4). The load was increased linearly to 90kN/m 2 until
t=0.3second. Only the upper surface was drained. Two load cases
are shown in Fig.4. In case 1, the load was applied to an area with
1×1 m2 in corner. In case 2, the load is applied to an area with

Fig. 4

Local compression of soil with Loads 1 and 2

41

42

X. TANG AND T. SATO

(32) for two load cases. The rate of convergence in the L 2-norm ofstrain is of the order O(h p). p is the degree of the polynomial utilized in the finite element approximation, for linear element p=1; hthe size of the element; and hndim is proportional to the volume,
where ndim is the number of spatial dimensions in the problem (inthree dimensions, ndim=3). The number of nodes or elements isinversely proportional to the volume. The rate of convergence
therefore is O(h1/3). The rates of convergent obtained by this error
estimator are almost same as the theoretical values.

Fig. 8

Mesh, load, and error contour in Load 2 (27 elements)

Fig. 9Mesh, load, and error contour in Load 2 (216 elements)

Fig .10

Mesh, load, and error contour in Load 2 (729 elements)

Fig. 11

Convergence of error estimator in three-dimensions

4. ADAPTIVE MESH REFINEMENT
A fission procedure belonging to h-refinement was applied toadaptive mesh refinement. Details of this procedure, including
implementation for quadrilateral elements in two-dimensions andfor hexahedral elements in three-dimensions, are presented here.
Its application to soil-pile interaction analysis also is shown. A
simple example was calculated to test the use of adaptive meshrefinement with the fission procedure.
4.1 3-D fission procedure for linear hexahedra elements
Linear hexahedral elements were used in the three-dimensional adaptive finite element analysis. When element error exceeded
an acceptable limit, the element was fissioned into eight smaller
elements, as shown in Fig.12. The circled numbers show element
numbers, the other numbers node numbers. The initial meshshown is mesh-a. After error estimation, the relative error for element 3 exceeded the error limit, and the element was fissioned into
eight smaller elements: 3, 8-11. Nineteen new nodes, 19-37, werecreated: in the middle of each of the twelve sides, at the centers of
the six surfaces and center of each element. The new refined meshshown is mesh-b. The parameters of element 3 in mesh-a weretransferred to elements 3, 5-11 in mesh-b, and the variables of the
new elements interpolated from the variables of element 3 inmesh-a. Pore pressure values of the new smaller elements are thesame as the value of the old parent element. The displacement,
velocities, and accelerations of the new nodes also were interpolated from values of the old nodes in mesh-a. The next calculation
step is based on mesh-b.

Fig. 12Fission procedure of hexahedral elements

ADAPTIVE MESH REFINEMENT AND ERROR ESTIMATE FOR 3-D SEISMIC ANALYSIS OF LIQUEFIABLE SOIL CONSIDERING LARGE DEFORMATION

When an element is fissioned next to an unfissioned one, slavenodes are created that are constrained by the compatibility condition of the master node. There are two kinds of slave nodes. One
is a slave node on the side of an element (Fig.13). Motion of slave
node 0 should be governed by master node 1 and node 2 as
V3 = T{V1

V2 }

T

(42)

where T is defined by [I/2,I/2], in which I is a unit matrix.
The equation of motion is not evaluated at the slave node.
Instead, nodal forces at the slave nodes are added to the forces atthe corresponding master nodes;{F

1F2 } = {F1*F2 } + TF3(43)

where F3 denotes nodal forces at the slave node and {F}* those atnodes 1 and 2 prior to the consideration of F3.
The secong type is a slave node on the surface of the element(Fig.14). Motion of slave node 0 should be governed by the master node 1 up to node 4;

Fig. 13

Slave node on the side of an element

Fig. 14

Slave node on the surface of an element

Table 1V5 = T{V1

V2

V3

V4 }

T

43

(44)

where T is defined by [I/4,I/4,I/4,I/4], in which I is a unit matrix.
The equation of motion is not evaluated at the slave node. Instead,nodal forces at slave nodes are added to the forces at the corresponding master nodes;

{F

1

F2

F3

F4 } = {F1

F2

F3

*

F4 } + TF5(45)

where F5 denotes nodal forces at the slave node and {F}* those atnodes 1 and 2 prior to the consideration of F5.
An array NABOR (NE, I, J) is used to define the relationshipof neighbor elements. It is composed of the twenty-four neighboring elements on the six joint surfaces with other elements. NE iselement number. I surface number and J neighbor element number
defined in Fig.15.
To clarify use of the array, consider mesh-a and mesh-b(Fig.12) and the fission of element 3. NABOR arrays for elements
1 and 3 are given in Table 1 for both meshes. Zero indicates theboundary at the surface of an element. The four identical numbers
on a surface indicate only one neighbor element on that surface.
When an element is beside fissioned elements, all the numbers of
the four neighboring elements on this surface are recorded in thearray.
4.2 Modifying permeance of pore water between hexahedral
elements
If an unfissioned element is beside a fissioned element, theusual permeance route crossing the intersurface should be modified, as shown in Fig.16. The place of permeance between two
unfissioned elements is taken by the permeance between an unfissioned element and four fissioned elements. The terms describing
the permeance of pore water in the continuity equation therefore
are modified easily.

Fig. 15

Definition of I and J for array NABOR (NE,I,J)

Nabor (JE,I,J) for mesh-a and mesh-b

44

X. TANG AND T. SATO

Fig. 16

Fig. 17

Fig. 18

Increase in node number

Fig. 19

Decrease in relative error

Modifying permeance of pore water in three dimensions

Refined meshes and relative error contours at different steps

4.3 3-D Soil block compression
A simple three-dimensional example was analyzed to test theadaptive mesh refinement method. The example model is same asthat used to check the error estimator in section 3.3. The loading
condition is Load 1. h-Adaptive FE analysis starts from an initialmesh with 27 elements, and the relative error limit is defined as
0.02. Mesh refinement starts at t=0.15 second. The time increment of finite element calculation is 0.01 second. Mesh refinement is implemented every five steps in the finite element calculation. The initial and refined meshes at three steps of adaptive meshrefinement with relative error contours are shown in Fig.17.
Elements whose relative errors are larger than the 0.02 limit are
refined in every step of adaptive refinement.
Fig.18 shows that the number of elements increases rapidly.
At the same time, relative error decreases, as shown in Fig.19 inwhich amr denotes results of the adaptive FE method. Because
deformation is very small, relative error does not appear toincrease evidently. Fig.20 shows the efficacy of the method. The
results obtained show that this adaptive scheme produces substantial improvements in accuracy with limited computation.

Fig. 20

Convergence of error

5. NUMERICAL EXAMPLE
Three-dimensional seismic analysis of an embankment by theh-adaptive FE method is introduced. The updated Lagrangian formulation was adopted for the governing equations of the finite element method to deal with large deformation of liquefied soil. An
effective cyclic elasto-plastic model was used to describe the material nonlinearity of saturated soil.
An embankment constructed on saturated sand is shaken by anearthquake in the x-direction. As shown in Fig. 21, the embankment is 8m wide at the top and 16m wide at the bottom. There is a
type of filling soil without pore water in this embankment composed of Edosaki sand 1. The foundation soil region assessed isADAPTIVE MESH REFINEMENT AND ERROR ESTIMATE FOR 3-D SEISMIC ANALYSIS OF LIQUEFIABLE SOIL CONSIDERING LARGE DEFORMATION

Fig. 2145

Embankment example of 3-D analysis considering large deformation

Table 2

Soil parameters of the embankment example

48m wide, 12m thick, and 16m deep. If an area larger than this isconsidered, deformation caused by liquefaction obviously will belarger but only by a small extent. There are two kinds of saturated
soil, Edosaki sand 2 and silica, in the soil layers. The upper layer
is composed of Edosaki sand 2 that is 12m deep. The bottom layerof silica is 4m deep. Dynamic parameters of these soils are given
in Table 2. Soil parameters were obtained experimently and by
evaluation. The parameter definitions are those introduced by Oka
et al. (1992). Initial stresses of the elements were calculated forgravity. A linear hexahedral element with eight nodes was used.
Displacement of the bottom boundary was fixed, the side boundary
was fixed horizontally but allowed to slide on the surface. Drainage
was allowed only on the top boundary surface of the mesh.
Input earthquake acceleration was a modified earthquake
acceleration history at 60% amplitude of the original record, Fig.
22. The original acceleration was recorded in the 1995 HyogokenNanbu earthquake. The maximum value is 433 gal. Fig. 23 gives
the time history of the EPWPR of soil 6m deep. Reduction of the

46

X. TANG AND T. SATO

effective soil stress due to the increased excess pore water pressure
ratio leads to significant loss of soil strength and stiffness. Whenthe excess pore water pressure ratio reaches 1.0, full liquefaction
occurs at about 6 second.

Fig. 22

Fig. 23

Input acceleration for a 3-D embankment

EPWPR response of the 3-D embankment example

Fig. 24 Coarse mesh and relative error results for the 3-D
embankment

The deformed meshes and their respective relative errors distributions at t=7.0, 11.0, and 15.0 seconds analyzed with a fixedcoarse mesh (168 elements) are given in Fig.24. (a), (b) and (c).
The L 2-norm of strain was used to evaluate the error. Deformationand strain in regions besides the foot of and under the embankment
are larger than in other area, and relative error values in those area
are increased.
The adaptive mesh refinement process starts at t=6.5 seconds
for the initial coarse mesh with 168 elements. The relative error
limit is 0.07. Adaptive mesh refinement was done once per second, and the time increment of finite element calculation was
0.005 second. Only two levels of refinement were allowed, whichmeans that an element is fissioned only twice. The refined meshes
analyzed by the h-adaptive FEM using the updated Lagrangian
method at times t=7.0, 11.0, and 15.0 seconds respectively are
shown in Fig.25. (a), (b) and (c).
Elements in regions with large errors were refined step bystep. The relative error value of the refined mesh decreased as theelements were fissioned. In three-dimensional adaptive analysis,as the mesh becomes fine, the number of elements increases rapidly. No relative error value larger than the 0.07 limit is found in themesh with an adaptive level lower than 2 in Fig. 25. (a), (b) and (c).

Fig. 25

Final mesh with adaptive mesh refinement

ADAPTIVE MESH REFINEMENT AND ERROR ESTIMATE FOR 3-D SEISMIC ANALYSIS OF LIQUEFIABLE SOIL CONSIDERING LARGE DEFORMATIONThis example also was analyzed with a fixed fine mesh of1344 elements and cubic, liquefiable soil elements with 2m-longsides. Meshes with relative error distributions at t=7.0, 11.0 and
15.0 seconds are shown respectively in Fig. 26. (a), (b) and (c).
The relative error value of the fine mesh is lower than that of the
fixed coarse mesh, and the top surface of the liquefiable soil muchsmoother. The high relative error value distribution also is concentrated in regions near the foot of or under the embankment where
the elements have undergone large deformation and strain. The
meshes of the adaptive FE analysis are similar to the meshes calculated with the fixed fine mesh. Clearly, the refinement of these elements improves the accuracy of embankment displacement results.
In Figs. 27 and 28, horizontal and vertical displacements on
the top of the embankment use for the adaptive FE analysis arecompared with those for finite element analysis with fixed coarse
and fixed fine meshes without the adaptive procedure consideringlarge deformation. Horizontal displacement on the top of theembankment, calculated with the fixed fine mesh, is larger thanwith the coarse mesh (Fig. 27). In contrast, vertical displacementon the top of the embankment calculated with the fixed fine meshis smaller than with the coarse mesh (Fig. 28).
As mesh refinement proceeds, horizontal displacement calculated by the adaptive FE method becomes large, and vertical dis-placement small and closer to the displacement values of the fixed
fine mesh. As mesh refinement was not carried out before t=6.5
seconds, it is reasonable that the adaptive results are not as accurate as fine mesh results. From a comparison of the horizontal andvertical displacement lines for the three cases, it concluded that theh-adaptive FE method, which considers large deformation, effectively improves accuracy.
In Fig. 29, accuracies of three different cases are compared byshowing the three the average relative errors curves. As the quality of the entire mesh affects seismic analysis results, the general
quality of the entire mesh is determined by evaluating the averagerelative error. The average relative error of the fixed coarse meshwith 168 elements is larger than that of the fixed fine mesh with
1344 elements. In the adaptive FE analysis, the average relative
error value decreased once the adaptive procedure started at t=6.5
seconds. As the adaptive process proceeds, the average relative
error value approaches the value for the fixed fine mesh.
Reduction of the average relative error confirms the efficacy of theFig. 27 Horizontal displacement on the top of the embankment for 3-D analysis

Fig. 28 Vertical displacement on the top of the embankmentfor 3-D analysis

Fig. 26

Fixed fine mesh and relative error results

47

Fig. 29 Average relative error of 3-D analysis of theembankment

48

X. TANG AND T. SATO

adaptive FE method.
By this example, the h-adaptive FE method was shown to beeffective for three-dimensional elasto-plastic analysis of anembankment constructed on liquefiable soil. Moreover it improves
the accuracy of the finite element analysis by reducing the size oflarge error elements.
6. CONCLUSION
The h-adaptive technique was applied to non-linear FE analysis of saturated soil considering large deformation that includes theliquefaction phenomenon. In previous research, adaptive techniques have been used to solve static or transient problems withmonotonous loads. In this study, the adaptive FE method wasapplied to seismic analysis in which the load was an earthquake
acceleration history that included a few thousands steps of timeincrements and cyclic vibration. The fission procedure belonging
to h-refinement, indicated by element error estimate in the threedimensional adaptive FE method was used. The approximation
was refined successively, satisfying the predetermined standard ofaccuracy, and the efficacy of this method in finite element analyses
was confirmed. This method is easily used to solve practical andengineering problems. In the transient analysis, the relative errors
limits used ranged from 0.01 to 0.05. Although there have beenfew studies on the use of adaptive mesh refinement in dynamics
analysis, the reported relative errors limits in dynamics analysis
have ranged from 0.02 to 0.07. Selection of an appropriate adaptive frequency also is recommended.
A posteriori error estimate based on the L 2 norm of strain
error was used in the adaptive FE method. It effectively estimatedelements error after each calculation step in the nonlinear FE
analyses of soil. The superconvergent patch recovery technique
was used for the error estimates. The convergence of this error
estimate method was tested with three-dimensional examples anddifferent quality meshes. The convergence ratios matched the theoretical values. This method is readily implemented for any code,and calculations based on it are very simple. Moreover, the advantage of saving computation time is evident, and the technique is areliable indicator in mesh refinement.
The h-adaptive FE method is well suited to the elasto-plastic
analysis of saturated soil, especially to problems that include liquefaction analysis. Refining elements with large deformation orlarge error due to liquefaction flow effectively improves the accuracy of nonlinear numerical analysis. The h-adaptive FE method
also is useful for three-dimensional dynamic analysis. The application of adaptivity to the three-dimensional nonlinear analysis ofsaturated soil is valuable for the development of adaptive techniques. By giving the limit of the adaptive level and selecting asuitable adaptive frequency, the calculation time is reduced andaccuracy is increased. All the results obtained in this study showthat this adaptive scheme produces substantial improvements inaccuracy under a limited computational effort.
REFERENCES
Akai, K. and Tamura, T. (1978). Study of two-dimensional consolidationaccompanied by an elastic plastic constitutive equation, Proceedingsof JSCE, No. 269, 98-104.
Belytschko, T. and Tabbara, M. (1993). H-adaptive finite element methods
for dynamic Problem, with Emphasis on Localization, International
Journal for Numerical Methods in Engineering, Vol. 36, 4245-4265.
Belytschko, T., Wong, B. L. and Plaskacz, E. J. (1989). Fission-fussion
adaptivity in finite elements for nonlinear dynamics of shells,
Computers & Structures, Vol. 33, 1307-1323.
Biot, M. A. (1962). Mechanics of deformation and acoustic propagation in
porous media, Journal of Applied Physics, Vol. 33, 1482-1498.
De, J. P., Gago, S. R. Kelly, D., W., Zienkiewicz, O., C. and Babuska, I.
(1983). A posteriori error analysis and adaptive processes in the finite
element method: Part II-Adaptive mesh refinement, International
Journal for Numerical Methods in Engineering, Vol. 19, 1621-1656.
Di, Y. and Sato, T. (2001). Large deformation dynamic analysis of porous
medium using FEM-FDM coupled method, Advances in
Computational Engineering & Sciences (Atluri, S.N.,Nishioka, T. andKikuchi, M.,Eds.), Chapter 12 Multiphysics & Multibody Dynamic,
Tech Science Press (CD-ROM).
Hamada, M., Sato, H. and Kawakami, T. (1994). A consideration of themechanism for liquefaction-related large ground displacement, A proceedings from the fifth U.S.-Japan workshop on earthquake resistantdesign of lifeline facilities and countermeasures against soil liquefaction, Technical Report NCEER-94-0026, 217-232.
Hicks, M., A. (2000). Coupled computations for an elastic-perfectly plastic
soil using adaptive mesh refinement, International Journal for
Numerical. Analysis Methods of Geomechnics, Vol. 24, 453-476.
Kelly, D.W., De, J. P., Gago, S., R. Zienkiewicz, O., C. and Babuska, I.
(1983). A posteriori error analysis and adaptive processes in ihe finite
element method: Part I: Error analysis, International Journal for
Numerical Methods in Engineering, Vol. 19, 1593-1619.
Oka, F., Yashima, A., Shibata, T., Kato, M. and Uzuoka, R. (1994). FEMFDM coupled liquefaction analysis of a porous soil using an elastoplastic model, Applied Scientific Research, Vol. 52, 209-245.
Oka, F., Yashima, A., Kato, M. and Sekiguchi, K. (1992). A constitutive
model for sand based on the non-linear kinematic hardening rule andits application, Proceedings of 10th World Conference on earthquake
Engineering, Barcelona, 2529-2534.
Tang, X. and Sato, T. (2004). A posteriori error estimate and h-adaptive FE
analysis of Liquefaction with large deformation, Journal of
Geotechnical Engineering, No. 788/III-69, 1-15.
Tang, X. and Sato, T. (2004). H-adaptivity applied to liquefiable soil innonlinear analysis of soil-pile interaction, Proceedings of the 11th
International Conference on Soil Dynamics & Earthquake
Engineering & The 3rd International Conference on Earthquake
Geotechnical Engineering, 760-767.
Zienkiewicz, O. C. and Zhu, J. Z. (1992).The supercovergent patch recovery and a posteriori error Estimates. Part I: The Recovery Technique,International Journal for Numerical Methods in Engineering, Vol. 33,
1331-1364.
Zienkiewicz, O. C. and Zhu, J. Z. (1992).The supercovergent patch recovery and a posteriori error estimates. Part II: Error estimates and adaptivity, International Journal for Numerical Methods in Engineering,Vol. 33, 1365-1382.
Zienkiewicz, O. C. and Zhu, J. Z. (1992).The superconvergent patch recovery(SPR) and adaptive finite element refinement, Computational Methods
and Application of Mechanics Engineering, Vol. 101, 207-224.