A numerical model to investigate the effects of propagating liquefied soils on structures

Computers and Geotechnics 33 (2006) 108–120 www.elsevier.com/locate/compgeo

A numerical model to investigate the effects of propagating liquefied soils on structures Sami Montassar, Patrick de Buhan

*

Laboratoire des Mate´riaux et Structures du Ge´nie Civil (LCPC, ENPC, CNRS UMR 113), Ecole Nationale des Ponts et Chaussees, 6 et 8 av. Blaise Pascal, Cite Descartes, Champs/Marne, 77455 Marne-la-Valle´e Cedex 02, France Received 23 May 2005; received in revised form 3 February 2006; accepted 8 February 2006 Available online 18 April 2006

Abstract This contribution is concerned with the development of a computational method for simulating the propagation of a liquefied soil, modelled as a Bingham material, while evaluating its effects on adjacent structures. The method is based upon the existence of a minimum principle for the velocity field governing the instantaneous evolution of the liquefied soil mass, which is solved by resorting to a specific numerical algorithm implemented within the context of a finite element formulation. Making use of simplifying assumptions, the calculation of the drag force induced by the moving soil mass on structures schematized as rigid obstacles, can be performed independently, exploiting the same numerical tool. The whole procedure is illustrated on the typical example of a liquefied soil embankment interacting with a single row of regularly spaced foundation piles placed across the propagating soil mass.  2006 Elsevier Ltd. All rights reserved. Keywords: Liquefied soil; Residual shear strength; Viscosity; Bingham model; Drag force; Minimum principle; Finite element method; Decomposition– co-ordination method

1. Introduction Initial flow failure and subsequent propagation of liquefied soil masses triggered by the well-known seismically induced liquefaction phenomenon of saturated granular soils, may have very serious consequences on the stability of civil engineering structures, such as bridge piers, foundation piles or lifelines facilities. It is therefore essential to devise rational design methods capable not only to predict how such liquefied masses actually propagate, but also to quantitatively assess the potential damage they may cause to structures. A relatively abundant literature has been devoted to dealing with the first aspect of this question. In an attempt to estimate the maximum permanent ground displacement due to an earthquake, the following empirical formula has been for instance proposed by Hamada et al. [17] from exploiting *

Corresponding author. E-mail address: [email protected] (P. de Buhan).

0266-352X/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2006.02.003

extensive measurements performed on the site of the 1964 Niigita earthquake and 1983 Nihonkai-Chubu earthquake: ffiffiffi pffiffiffiffip 3 D ¼ 0:75 H h ð1Þ where D (m) is the maximum horizontal displacement, H (m) the depth of the liquefied soil layer and h (%) the greater of the slope of ground surface and lower boundary face of liquefied soil. Such a formula also incorporates the observations made on the San Fernando earthquake [9]. Likewise, Bartlett and Youd [4] have derived the same kind relationship based upon data gathered on six earthquakes in the United States, including the one having occurred in Alaska in 1964. Even though they remain empirical, that is based upon the fitting of observed field data, such relationships provide a first attempt to devise predictive methods for estimating the amount of permanent ground displacement induced by seismic liquefaction. In order to avoid the use of complex finite element non linear simulations, simplified analytical models and related calculation methods incorporating the rheological

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

characteristics of the liquefied soil, have been developed by Towhata et al. [38–40], in the case of a liquefied soil layer topped by a non liquefiable soil, in which lateral flow is predominant. The method, which is based on the application of the minimum principle for the potential energy, is focused on the evaluation of the maximum displacements experienced by the liquefied ground under the action of gravity. It relies upon the following assumptions, some of them being supported by experimental evidence.  The lateral displacement is supposed to vary as a sinusoidal function of depth, as suggested for instance by the observations of centrifuge tests, leading to a much easier numerical treatment of the problem than in the case of finite element simulations. This makes it possible to generalize the method to three dimensional situations, as it has been done by Orense and Towhata [28].  The liquefied soil behaves as a Bingham or Newtonian incompressible fluid, which may exhibit viscous as well as residual shear strength properties. A constitutive relationship of the following form is notably proposed by the authors [38,39] s¼G

ou þ sr oz

ð2Þ

where u is the lateral permanent displacement, function of the depth z, s (respectively sr) is the shear stress (respectively residual shear strength) and G the liquefied soil shear modulus. While the first assumption seems quite legitimate, the adoption of a constitutive law such as (2) for modelling the constitutive behaviour of liquefied soil may appear somewhat questionable (apart from the particular situation where G = sr = 0), unless the displacement u is replaced by the velocity (and the corresponding shear strain by the shear strain rate), whereas G would denote the viscosity coefficient instead of the elastic shear modulus (see Eq. (3) below). It is widely acknowledged today that, in accordance with the terminology used, a ‘‘liquefied’’ soil behaves as an incompressible fluid, displaying both viscous and residual strength properties. Ample experimental evidence of the viscosity of a liquefied soil, evaluated by different measurement techniques (viscometer, lateral spreading, load acting on a sphere, etc.) has been provided, the synthesis of which may be found in Uzuoka et al. [42], Towhata et al. [38] or Hadush et al. [15,16]. The same kind of experiments makes it possible to evaluate the undrained residual shear strength of a liquefied soil [35,21,11]. This clearly indicates that the rheology of a liquefied soil under pure shear conditions, may be appropriately described by the following constitutive equation: s ¼ l_c þ k

ð3Þ

where k is the residual shear strength, c_ the shear strain rate and l the viscosity coefficient. This is the characteristic

109

equation of a ‘‘yield stress’’ fluid or Bingham material, adopted by Uzuoka et al. [42] for simulating the post liquefaction behaviour of liquefied ground, or for dealing with closely related problems, such as flowslides [29], slurry flows [6], or mudflows [7,18]. Several simulation techniques have been developed using such a constitutive equation, namely the so-called ‘‘shallow water approximation’’ [22] or the ‘‘depth integrated method’’ [6,29], based upon the simplifying assumption that the flow is predominantly parallel to the slope, making it possible to account for very large displacements of the laterally spreading soil mass, inducing considerable changes of geometry. In the more general situation where such a simplification cannot be made, specific simulation techniques should be developed such as fluid dynamics based numerical methods [42], finite element [37] or cubic interpolated pseudoparticle (CIP) method [16]. It should be kept in mind however, that the adoption of a viscoplastic Bingham model, expressed by Eq. (3), is only justified in the situation when the consolidation time which governs the pore pressure dissipation, may be considered large enough in comparison with the time of propagation. Indeed, in the case when these two times are of the same order, the excess pore pressure which is at the origin of the liquefaction phenomenon, may significantly dissipate during the propagation, thus allowing the granular soil to recover its frictional properties, bringing the moving soil mass to stabilize again. The fundamental assumption upon which all the previously described approaches, as well as the numerical scheme developed in this paper, are based, is that the rheological properties of the liquefied soil, which may be conveniently captured by means of a Bingham model, remain constant throughout the whole propagation process. Analysing the potentially damaging effects of a propagating liquefied soil on the neighbouring structures, requires being able to evaluate the corresponding additional loads applied to these structures by the moving soil mass. Although being essential for engineering design purposes, this second aspect has almost exclusively been treated from an experimental point of view. One may quote for instance the experimental work by Towhata and Hussaini [38], who performed tests aimed at evaluating the lateral loading exerted by a submarine mudflow on offshore structures, as also the contributions of Yasuda et al. [44] and Towhata et al. [41] concerning shaking table tests carried out on buried pipelines. Several experimental studies have also been undertaken on reduced scale centrifuged models, in order to investigate the lateral resistance of foundation piles subject to liquefaction [5,1,19,43]. To the authors’ knowledge, apart from the works of Dobry et al. [8] and Ishihara and Cubrinovski [20], who have proposed simplified design analyses of piles subject to liquefaction induced lateral loading, no other significant contributions can be quoted as regards the numerical modelling of this kind of problem. The present contribution is concerned with the presentation of a fully integrated numerical tool allowing to

110

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

simultaneously analyse the post liquefaction behaviour of a water-saturated soil mass and its consequences in terms of additional loadings exerted on the surrounding structures. The whole procedure, which relies on the application of a minimum principle for the velocities, analogous to that employed by Towhata et al. [38–40], is developed in the context of a liquefied soil modelled as a viscoplastic Bingham material. It is worth emphasizing that the ability of such a numerical method to produce reliable predictions, is subject to the condition that the extent and depth of the liquefied zone, along with its viscous and residual strength characteristics, are being known quantities, incorporated as input data in the analysis. 2. Problem statement Analysing the post-liquefaction behaviour of a saturated soil mass entails two different aspects, which are sketched in Fig. 1.

Fig. 1. (a) Evolution problem of a liquefied soil mass; and (b) evaluation of its action on structures.

(a) The liquefied soil, whose strength has been considerably reduced as a consequence of the liquefaction process, is first set into motion under the action of gravity, undergoing very large displacements as it propagates downhill. It is therefore essential to devise appropriate methods for tracking the geometry changes of the moving soil mass as a function of time: this is the evolution problem. (b) Assuming that the evolution problem is being solved, which means that the successive configurations of the liquefied soil mass are known at any time along with the attached velocity field, it is essential to evaluate the efforts (drag forces) exerted by the flowing soil mass on neighbouring civil engineering structures, such as bridge piers or deep foundations piles. In the general case, the comprehensive treatment to this twofold problem requires a fully coupled complex analysis, which would lead us far beyond the scope of engineering design methods. The present contribution is devoted to the development of a simplified evaluation procedure aimed at dealing with the particular configuration of a single row of vertical cylindrical obstacles placed across the flow of a liquefied soil, as shown in Fig. 2. The following simplifying assumption are made throughout the paper.  The obstacles are fixed, vertical and remain undeformed under the action of the liquefied soil. This hypothesis should be reconsidered in the case of obstacles such as for example flexible foundation piles, the deflection of which would of the same order of magnitude as the liquefied soil displacements.  The presence of such obstacles does not significantly affect the propagation of the liquefied soil. This means from a practical point of view that the spacing 2H

Fig. 2. Different views of a row of fixed cylindrical obstacles placed across a propagating liquefied soil mass.

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

111

between two adjacent obstacles may be considered as significantly larger than their diameter D = 2R. Such an hypothesis implies that the problem of predicting the evolution of the liquefied ground can be solved independently from the calculation of the efforts on the obstacles.  Furthermore, let F(y) denote the density per unit vertical length of the drag force horizontal component exerted by the liquefied soil on each individual obstacle. This density may be evaluated from the solution to the twodimensional stationary flow problem sketched in Fig. 2c, where uniform velocity boundary conditions are imposed at sufficient upstream and downstream distances from the obstacle. Its value, considered as an input parameter for the stationary flow problem, is equal to the horizontal component Ux(y) of the velocity of the liquefied ground at the point where the obstacle is located, which has been previously calculated from the evolution problem.

uefied soil mass, subject to the sole action of gravity, is characterised by the velocity field U(x,t) defined at any point x on its current configuration Xt at time t P 0. As sketched in Fig. 3 the external boundary oXt of Xt is subdivided into two complementary surfaces oXrt and oXUt . oXrt remains stress-free, that is:

Solving the evolution problem of the liquefied soil mass as well as the stationary flow problem around the obstacle in order to calculate the drag force, requires to specify the rheological behaviour of the liquefied soil. As previously discussed in the introduction, the Bingham model appears to provide a realistic description of such a post-liquefaction behaviour, both in terms of residual shear strength and viscosity exhibited by such kind of material. Expressed in full tensorial form as a stress–strain rate relationship, such a constitutive behaviour may be written as (see [25], for more details):

8t; 8x 2 Rt ;

pffiffiffi d r ¼ pI þ 2ld þ k 2 kdk

if

d 6¼ 0

ksk f ðrÞ ¼ pffiffiffi  k 6 0 2 div U ¼ tr d ¼ 0

d¼0

if

ð4Þ

8t;

8x 2 @Xrt ;

rðx; tÞ  nðx; tÞ ¼ 0

ð7Þ oXrt ,

while slipwhere n(x,t) is the external unit normal to page boundary conditions proposed by Fortin et al. [13], are imposed on Rt which denotes the interface between the liquefied soil and the rigid substratum, located within the liquefied soil mass along the boundary surface oXUt (Fig. 3). The velocity V of any point belonging to the liquefied soil mass at the interface Rt, which is equal to the velocity jump across Rt following the unit normal m, since the velocity of any point of the rigid substratum is zero, must remain tangential V ðx; tÞ  mðx; tÞ ¼ 0

ð8Þ

and the following slip boundary conditions are applied on Rt (Fig. 4): s ¼ l V þ k 

V if V 6¼ 0 jV j

V ¼ 0 if jsj 6 k 

ð9Þ ð10Þ

where s = T  (T Æ m) m is the shear stress vector exerted by the liquefied soil on the substratum, k* is the interface shear strength below which no slippage occurs, and l* is a viscosity coefficient. The particular case of the perfect bonding condition is recovered when taking k* ! 1. The similarity between the above boundary conditions and the constitutive Eqs. (4)–(6) of the Bingham model is to be noticed,

ð5Þ ð6Þ

In these equations, r is the Cauchy stress tensor, s its deviatoric part, p =  1/3tr r is the mean pressure, which can be interpreted as the Lagrange multiplier associated with the incompressibility condition (6), where d is the strain rate tensor, and f(r) is the von Mises function which involves the pffiffiffiffiffiffiffiffi norm of the stress deviator ksk ¼ s : s. It is worth noting that the purely viscous Newtonian model is recovered from the Bingham model when the shear strength k is reduced to zero, while the so-called ‘‘rigid-plastic’’ model is obtained when the viscosity coefficient l tends to zero.

Fig. 3. Initial and current configurations of a liquefied soil mass subject to gravitational forces.

3. A velocity minimum principle for solving the evolution problem The liquefied soil mass is described within the framework of three-dimensional continuum mechanics as occupying a geometrical domain X0 in its initial configuration at time t = 0, immediately after the liquefaction process has been triggered. The subsequent propagation of this liq-

Fig. 4. Imposition of a non-linear friction boundary condition allowing for possible slippage at the liquefied soil/substratum interface.

112

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

because it will be helpful in the numerical treatment of the problem. Let us now introduce the sets S(t) of statically admissible stress fields r 0 , and K(t) of kinematically admissible incompressible velocity fields U 0 defined, respectively as: ( div r0 þ qg ¼ 0 8x 2 Xt r0 2 SðtÞ () ð11Þ r0  n ¼ 0 8x 2 @Xrt  div U 0 ¼ 0 on Xt 0 U 2 KðtÞ () ð12Þ V 0  m ¼ 0 on Rt A velocity field U 2 K(t) is a solution for the evolution problem if it is associated with a stress field r 2 S(t) through the constitutive Eqs. (4), (5), (9) and (10). It has been proved in [24,25] that such a solution U realises the minimum of the following functional defined on K(t): F ðU Þ ¼ MinfF ðU 0 Þ; U 0 2 KðtÞg with F ðU 0 Þ ¼

Z

Z



l 2 2 jV 0 j dRt  lkd 0 k dXt þ Rt 2 Z pffiffiffi Z þ k 2kd 0 kdXt þ k  jV 0 jdRt Xt

Xt

ð13Þ Z

qg  U 0 dXt Xt

rate field and velocity jump distribution across the soil/substratum interface, as minimization variables denoted by d 0 and v 0 , respectively, linked to the velocity field by the following constraints: d0 ¼ dðU 0 Þ ¼ 1=2ðgrad U 0 þT grad U 0 Þ on Xt ð15Þ and v0 ¼ V 0 ; V 0  m ¼ 0; on Rt ð16Þ 0 0 Denoting by k and w the Lagrange multipliers associated with the above constraints (15) and (16), the minimisation problem defined by (10) and (11) leads to a generalised saddle-point problem of the form: Lr ðU  ; p ; d ; k ; v ; w Þ max min max min max Lr ðU 0 ; p0 ; d0 ; k0 ; v0 ; w0 Þ ¼ min 0 0 0 0 0 0 U

p

d

ð14Þ

The two first integrals appearing on the upper line of (14) are equal to half the viscous dissipation, while the third integral corresponds to the work of gravity, which is the unique external force. The lower line of integrals is equal to the plastic dissipation. This minimum principle is analogous to the minimum of potential energy invoked by Towhata et al. [38–40], which has been briefly described in Section 1. In the particular case when the residual shear strength is equal to zero, so that the liquefied soil is modelled as a Newtonian fluid, F reduces to a quadratic functional, and the principle is analogous to the minimum of potential energy established in the context of linear elasticity (see for instance [33]), except that the velocity field is substituted to the displacement field. In the general case, this functional is non-differentiable due to the presence of the plastic dissipation terms. The resolution of the minimisation problem (13) then requires a special treatment. The decomposition–co-ordination method by augmented Lagrangian, proposed by Fortin and Glowinski [12], has proved its efficiency for dealing with such a numerical problem (see [24,25] for a comparison with other numerical methods). 4. Numerical treatment: the decomposition–coordination method by augmented Lagrangian The present section is devoted to a general outline of the numerical method adopted for performing the above minimization procedure, without going into details of any mathematical and theoretical justification, which may be found for instance in [12] Glowinski and Le Tallec [14], Fortin et al. [13] or more recently Roquet and Saramito [31]. Its principle is based on the introduction of the strain

w

ð17Þ where Lr ðU 0 ;p0 ;d0 ;k0 ;v0 ; w0 Þ ¼

Z

2

lkd0 k dXt þ

Z

X

Zt

pffiffiffi k 2kd0 kdXt

Xt 0

Z

p0 divU 0 dXt X Xt Z t Z  l 02 0 0 0 jv j dRt þ k : ðd  d ÞdXt þ Xt Rt 2 Z Z þ k  jv0 jdRt þ w0  ðV 0  v0 ÞdRt R Rt Z t r1 2 ðdivU 0 Þ dXt þ Xt 2 Z r2 0 2 kd  d0 k dXt þ Xt 2 Z r3 0 2 jV  v0 j dRt þ ð18Þ 2 R is the augmented Lagrangian functional defined for all r1 > 0, r2 > 0 and r3 > 0, whereas p 0 is the Lagrange multiplier associated with the incompressibility condition. It is to be noticed that, for any choice of parameters ~ ; ~p; ~d; ~k; ~v; w ~ Þ to the above sad(r1, r2, r3), the solution ðU dle-point problem (17) verifies in particular: ~ ¼U U ð19Þ 

Rt

v

k

qg  U dXt 

Such a solution can be calculated using the following generalised Uzawa’s algorithm proposed by Fortin and Glowinski [12] and adapted to our problem [24,25]  Step 1: Initialisation at iteration i = 0; p(0), k(0), w(0) and U(0) are arbitrarily specified.  Step 2: for any i P 1 : – Calculate d(i) as a function of k(i1) and U(i1) using the following relationship: * pffiffiffi + 1 k 2 ðiÞ d ¼ 1  ði1Þ bði1Þ with 2l þ r2 kb k bði1Þ ¼ kði1Þ þ r2 dðU ði1Þ Þ

and

hi ¼ supf; 0g ð20Þ

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

– Calculate v(i) as a function of w(i1) and U(i1) using :   1 k ðiÞ v ¼  1  ði1Þ xði1Þ with l þ r3 jx j xði1Þ ¼ wði1Þ þ r3 V ði1Þ

ð21Þ

– Update the Lagrange multipliers associated with constraints (15) and (16) : kðiÞ ¼ kði1Þ þ r2 ðdðU ði1Þ Þ  dðiÞ Þ ðiÞ

w

ði1Þ

¼w

þ r3 ðV

ði1Þ

v Þ

ð23Þ

– Find U(i) as a function of p(i1), d(i), k(i), v(i) and w(i) by minimising the following functional, which is quadratic with respect to U 0 : GðU ðiÞ ; pði1Þ ; dðiÞ ; kðiÞ ;vðiÞ ; wðiÞ Þ 8R R ði1Þ r1 0 2 0 > t > > Xt h2ðdivU Þ dXt  Xt p divU dX < i R r2 0 ðiÞ ðiÞ 0 ¼ Min ðdðU Þ  2d Þ þ k : dðU Þ dXt Xt 2 U0 > > > : þ R r3 ðV 0  2vðiÞ Þ þ wðiÞ  V 0  dR  R qg  U 0 dX Rt

2

terms have to be taken into account, notably in the initial phase of propagation, since such a principle still applies, provided that the acceleration field is being known. In such a case the linear system (26) becomes: ½DfU g ¼ fF g  fAg

t

Xt

9 > > > = t

> > > ;

ð24Þ – Update the Lagrange multiplier associated with the incompressibility condition: pðiÞ ¼ pði1Þ  r1 div U ðiÞ

ð25Þ

 Step 3: Application of convergence tests. According to Fortin and Glowinski [12], the explicit formula (20) (respectively (21)) represents the solution to the minimisation of Lr with respect to d (respectively v), with U, p, k, v and w (resp. U, p, d, k and w) being kept fixed. Owing to the decomposition procedure, the above described algorithm transforms the initial non-differentiable problem (13) into a family of classical problems (explicit computations for (20) and (21) and standard minimisation of a quadratic functional for (24)), co-ordinated via the Lagrange multipliers, updated at each iteration through Eqs. (22), (23) and (25). The minimisation procedure of the quadratic form (24) is carried out numerically in the context of the finite element method. The evolution problem being treated as a plane-strain two-dimensional problem, the volume of liquefied soil is discretized into a finite element mesh made of four-noded quadrilateral elements with linearly varying velocities, so that performing the minimization (24) reduces to solving a linear system of the general form: ½DfU g ¼ fF g ð26Þ where {U} is the column matrix formed by all the nodal velocities, and {F} represents the external forces (gravity).The choice of this kind of elements is primarily justified by the fact that it lends itself very easily to the simulation of the very large deformations experienced by the liquefied soil during it propagation [24]. The minimum principle and related numerical algorithm may be readily extended to the case when acceleration

ð27Þ

where {A} represents the acceleration terms, which may be approximated as follows at current time t:

ð22Þ

ðiÞ

113

fAgðtÞ ¼ ½M

fU gðtÞ  fU gðt  DtÞ Dt

ð28Þ

where [M] is the matrix of masses, and {U}(t) (respectively {U}(t  Dt)) is the vector of nodal velocities at time t (respectively at the previous time step t  Dt), so that the linear system to be solved becomes:

½M ½M ½D þ fU gðt  DtÞ ð29Þ fU gðtÞ ¼ fF g þ Dt Dt Finally, the simulation of the evolution problem is carried out in an incremental way, by updating the geometry of the propagating soil mass from time t to time t + Dt, according to the following explicit relationship: fxgðt þ DtÞ ffi fxgðtÞ þ DtfU gðtÞ

ð30Þ

where {x} denotes the set of co-ordinates of all the nodes of the finite element mesh, considered as material points. The whole above described numerical procedure for simulating the lateral flow of liquefied soil has been favourably compared [24] with the reference analytical solution obtained in the particular situation of a uniform layer of liquefied soil moving on an inclined substratum [24,26]. 5. Evaluation of the efforts exerted on structures The previously elaborated numerical method, based on the minimum principle for the velocities, is now applied to calculate the drag force exerted by the propagating liquefied soil on structures, in the particular situation shown in Fig. 2c, where it reduces to solving a quasistatic steady flow problem of a Bingham material around a circular obstacle in the (x, z)-plane. On account of the different symmetries exhibited by such a problem, it may be easily shown that the flow domain to be considered in the analysis may be restricted to the region of the plane defined by the following conditions: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 6 x 6 L; 0 6 z 6 H ; r ¼ x2 þ z2 P R ð31Þ where the horizontal extension L is significantly larger than H. The corresponding boundary conditions are (Fig. 5):  Uniform velocity imposed on the right-hand side x = L: Ux ¼ U; Uz ¼ 0

ð32Þ

 Left-hand side x = 0: T x ¼ 0;

Uz ¼ 0

ð33Þ

 Upper and lower sides z = 0, H: T x ¼ 0;

Uz ¼ 0

ð34Þ

114

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120 100

80

present simulations

60 ∗ (N

F

Mitsoulis (2004)

)

40

20

Fig. 5. Boundary conditions for analysing the stationary flow problem relative to the evaluation of the drag force. 0 0

 Contact surface with the obstacle r = R: U ¼0

The numerical simulations, have been carried out for L = 6H, varying the ratio H/R between 2 and 50. Fig. 6 shows a typical finite element mesh of the flow domain (H/R = 3). The drag force per unit transversal length along the y-direction may be calculated as follows: Z p=2 F ¼4 ex  rðr ¼ R; hÞ  er R dh ð36Þ 0

since r(R,h) Æer is the stress vector applied by the liquefied soil in each point (R, h) of the obstacle. 5.1. Newtonian model of liquefied soil (k = 0) It may be easily shown from dimensional analysis arguments that, in the case of a purely viscous Newtonian fluid, characterized by its sole viscosity coefficient l, the drag force density may be put in the following form: ð37Þ

where F ðN Þ is a dimensionless factor which is a function of the sole relative spacing H/R. Fig. 7 displays the evaluation of this factor (triangles) obtained from the finite element simulations. Following the suggestion made by Faxe´n [10] these numerical points

15

10

20

25

30

H /R

ð35Þ

Note that in conditions (33) and (34), Tx represents the horizontal component of the stress vector acting upon the corresponding boundary.

F ðN Þ ¼ lUF ðN Þ ½H =R

5

Fig. 7. Evaluations of the non dimensional drag force as a function of the relative spacing between obstacles: case of Newtonian fluid.

may be fitted by a curve, drawn in the same figure (solid curve), defined by the following equation: F ðN Þ ¼

4p 2

lnðH =RÞ þ a1 þ a2 ðR=H Þ þ a3 ðR=H Þ4 þ a4 ðR=H Þ6 þ a5 ðR=H Þ8

ð38Þ where the coefficients ai, i = 1, . . ., 5, involved in the above formula, have been calculated using the least square method: a1 ¼ 0:661;

a2 ¼ 1:985;

a4 ¼ 193:823;

a5 ¼ 443:505.

a3 ¼ 25:332; ð39Þ

Referring to the slightly different problem of the flow around a cylinder kept between parallel plates, Faxe´n [10], and quite recently Mitsoulis [23], have produced the same kind of evaluation reported in the same figure, where the squares correspond to the numerical simulations performed by Mitsoulis, while the dashed curve is obtained from fitting formula (38). As could be expected, the numerical simulation of the latter problem leads to slightly overestimating the non dimensional drag force, both estimates converging to zero as the relative spacing tends to infinity. 5.2. Viscoplastic Bingham model This kind of problem has already been investigated in the case of a non-viscous soil (‘‘rigid-plastic’’ model:

Fig. 6. Finite element mesh of the flow domain adopted for H/R = 3.

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

l = 0) by several authors such as Randolph and Houlsby [30], Murff et al. [27], or more recently Salenc¸on [34] in the context of the limit analysis or yield design theory. It has been extended to the case of a Bingham material by Adachi and Yoshioka [2] making use of variational principles. Let F(B) denote the drag force per unit length of the obstacle exerted by the liquefied soil modelled as a viscoplastic Bingham medium, characterized by its residual shear strength k along with its viscosity coefficient l. Resorting once again to dimensional analysis, such a drag force density may be written as: F ðBÞ ¼ lUF ðBÞ ½H =R; Bn 

ð40Þ

where F ðBÞ is the non dimensional drag force coefficient depending not only of the relative spacing H/R, but also on the dimensionless Bingham number defined as: Bn ¼

kR lU

ð41Þ

which compares the contribution of the shear strength to that of the viscosity of the Bingham material. The case of Newtonian fluid is obviously recovered when Bn ! 0. The numerical simulations where carried out for different values of the ratio H/R, varying the residual shear strength k, while keeping all the other parameters of the problem, namely R, l and U, constant and equal to unity. The results of such simulations are shown in Fig. 8, in the form of a series of curves giving the non dimensional drag force coefficient F ðBÞ as a function of H/R and Bn. As suggested by Zisis and Mitsoulis [45] who carried out the same kind of FEM simulations using a regularized Bingham model, the numerical points reported in this figure have been fitted by solid curves obeying the following analytical expression: b

F ðBÞ ðH =R; Bn Þ ¼ F ðN Þ ðH =RÞ½1 þ aBn 

ð42Þ

F ðBÞ

is given by Eq. (38), whereas coefficients a and b where are both functions of H/R only.

100000

115

Table 1 Coefficients a and b involved in Eq. (42) evaluated for different relative spacings of the obstacles H/R

F ðNÞ

a

b

2/1 3/1 4/1 10/1 15/1 20/1 25/1 30/1 35/1 40/1 50/1

52.059 23.117 15.886 7.58 6.318 5.701 5.437 5.212 5.056 4.936 4.774

0.531 1.197 2.127 7.345 9.835 15.021 18.88 19.457 25.89 25.23 29.486

0.963 0.958 0.881 0.801 0.855 0.796 0.788 0.801 0.762 0.784 0.764

All the results derived from these numerical simulations are summed up in Table 1 below. This table, combined with formula (38) along with (39), makes it possible to compute in a very straightforward manner, the drag force per unit length exerted by a propagating liquefied soil on a row of evenly spaced cylindrical obstacles. 6. An illustrative example The whole numerical procedure will now be illustrated on the example sketched in Fig. 9. The liquefied soil mass under consideration occupies in its initial configuration at time t = 0 a quadrilateral domain of height 1 m and horizontal extension equal to 3 m along its bottom, where it relies upon a horizontal substratum. Its left hand side also remains in contact with such a substratum, while no surcharge is applied to the upper and right hand surfaces. The only driving force is the liquefied soil specific weight equal to 18 kN/m2. A row of regularly spaced vertical cylindrical obstacles is placed at a horizontal distance of x = 2 m from the lower left hand corner of the embankment, taken as the origin of the axes. Their diameter is

F(B∗ ) 2m

0.2m

F (y)

10000

1000

y

H /R=2

100

2m

3 10

ρg = 18 kN/m3

10 50

Bn

1 0.001

00.1

0.1

1

F (y)

10 10

100

1m

1000

Fig. 8. Numerical evaluations of the non dimensional drag force coefficient as a function of the Bingham number for different values of the relative spacing between obstacles.

x 1

3m

Fig. 9. Liquefied soil embankment in its initial configuration.

116

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

equal to 0.2 m, while the interval between two neighbouring obstacles is taken equal to 2 m, so that H/R = 10. The different simulations are performed for a coefficient of viscosity equal to l = 1 kPa s, the undrained shear strength being progressively increased from 0 to 4 kPa. Referring to the values reported in Table 1, it follows that the drag force density applied by the propagating liquefied soil characterized locally by its horizontal velocity Ux (x = 2 m, y) at depth y is:  0:801 0:1k ð43Þ F ðyÞ ¼ 7:58U x ð2; yÞ 1 þ 7:345 U x ð2; yÞ where the value of F(y) is expressed in kN/m, while the undrained shear strength k and the velocity Ux(y) are expressed in kPa and m/s, respectively. A relatively coarse and regular mesh comprising 10 · 30 quadrilateral elements is used throughout the different simulations. 6.1. Influence of the undrained shear strength A first series of simulations is performed, assuming that the liquefied soil mass is perfectly bonded to the substratum along the lower and left-hand sides, which means that no velocity discontinuity is allowed along these boundaries. Adopting a time increment equal to Dt = 0.02 s, the evolution of the liquefied embankment is simulated up to t = 1 s. The corresponding deformed configurations of the embankment for increasing values of the undrained shear strength are pictured in Fig. 10. This figure clearly shows that the adoption of a Bingham instead of a Newton model for describing the rheological behaviour of the liquefied soil, that is the introduction of an undrained shear strength in addition to the viscosity, results in considerably limiting the amplitude of the geometry changes experienced by the liquefied soil. This is con-

firmed by the observation of the instantaneous velocity field represented by a distribution of arrows attached to the nodal points of the mesh: the volume of liquefied soil involved in the flow failure of the embankment, that is experiencing large deformations, is becoming less extensive as the liquefied soil strength increases, almost disappearing for k = 4 kPa. Furthermore, it should be noted that, contrary to what is observed, either in centrifuge experiments or by means of the present numerical procedure [24], in the case of an infinite layer of liquefied soil flowing down an inclined slope, the point of maximum velocity (and hence displacement) is not located on the ground surface but, as shown in Fig. 10 tends to approach the toe of the embankment, as the residual shear strength increases. This is of course to be attributed to the particular configuration of the liquefied soil mass, where the vertical facing remains stress free. This is corroborated by referring to a limit analysis (or yield design) point of view. Indeed, it may be shown [18,26], that viscoplastic flow failure of any liquefied soil structure will occur as soon as the latter is proved to be unstable in the sense of limit analysis. In the present case, a classical result states (see for instance [32]) that the stability of a vertical embankment of height h, as that considered in our simulations, is ensured as far as its stability factor remains lower than a critical value, which is close to 3.8: qgh 6 K ffi 3:8 k

ð44Þ

thus indicating that the liquefied soil embankment will not propagate as soon as the undrained shear strength exceeds the following value: kP

qgh 18  1 ffi ffi 4:73 kPa K 3:8

ð45Þ

Fig. 10. Deformed configurations at time t = 1 s of the liquefied soil embankment for different values of the undrained shear strength, and corresponding drag force profiles on the obstacles.

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

The distributions of the drag force F(y) exerted by the propagating soil on the row of vertical piles are pictured in the same figure. They are derived from the profiles of horizontal velocities of the liquefied soil along the vertical line where the obstacles are located, previously deduced from the evolution problem, by using formula (43). These profiles are decreasing from a maximum value reached at the top of the embankment to zero at the bottom. Moreover, it can be observed that, every other parameter being held constant, this profile increases in a first stage with the undrained shear strength, passing through a maximum for k = 1–2 kPa, then decreases again down to almost vanishing for k = 4 kPa. This could be tentatively explained from examining formula (43), which shows that the drag force is an increasing function of k, while at the same time the horizontal velocity component involved in this equation is a decreasing function of k. 6.2. Introduction of a friction boundary condition The perfect bonding condition imposed at the interface between the liquefied soil and the substratum in the previous set of simulations, may lead to somewhat unrealistic results as regards for instance the deformations experienced by the liquefied soil in the vicinity of the toe of the embankment, where excessive mesh distortion appears as time passes, leading to an ill-conditioned problem. This is why such a condition should be relaxed by the introduction of a friction boundary condition, allowing for possible slippage beyond a shear stress value. Therefore the same kind of simulation has been performed for the following set of rheological parameters:

l ¼ 1 kPa s; k  ¼ 1 kPa

k ¼ 2 kPa;

117

l ¼ 1 kPa s m; ð46Þ

It should be noted in particular, that the selected value of the soil-substratum strength parameter k* should be lower than that of the liquefied soil strength k, otherwise the introduction of such a friction boundary condition will have no influence of the liquefied soil propagation. Fig. 11 pictures a selected set of successive configurations of the embankment at regular time intervals between t = 0 s and t = 1.5 s. Comparing the configuration obtained for t = 1 s with that of Fig. 10 for k = 2 kPa, gives clear evidence of the influence of the friction boundary condition both on the amplitude of the liquefied soil propagation and on the drag force profile. As regards the first aspect, Fig. 11 shows that the liquefied soil is no more adherent to the substratum, the slippage condition being activated along the bottom horizontal substratum. This leads to much larger values for the velocities, and hence for the cumulated displacements, notably near the toe of the embankment, while the drag force profile is significantly different from that observed in the case of the perfect bonding condition. Furthermore, the extension of the volume of propagating liquefied soil is larger in this second simulation. It is to be noted that selecting appropriate values of parameters l* and k* to be incorporated into the analysis, remains a difficult task. This is due to the complete lack of available experimental data in the field of geotechnical engineering, whereas the introduction of such a slip boundary condition is of primary importance for the simulation of non-Newtonian fluids in the material process industries [13]. It could provide however, a relevant description for a

Fig. 11. Successive configurations of the liquefied soil embankment and associated drag force profiles as a function of time, taking a friction boundary condition into account.

118

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

40

quasi - static evolution 30

R (k N )

20

dynamic evolution

10

t(s)

0 0

0.5

1

1.5

2

2.5

3

Fig. 12. Evolution of the resultant drag force as a function of time: quasi-static vs. dynamic simulations.

thin layer of liquefied soil, modelled as a Bingham fluid, across which the velocity increases from zero to a finite value. Finally, in order to assess the role played by the acceleration forces in the propagation of the liquefied soil mass and consequently on the evaluation of the drag force exerted on the structures, curves giving the resultant drag force defined as: Z y¼lðtÞ RðtÞ ¼ F ðy; tÞ dy ð47Þ y¼0

where l(t) denotes the length of the obstacle remaining in contact with the liquefied soil embankment in its current configuration (l(t = 0) = 1 m), have been drawn in Fig. 12 as a function of time up to 3 s. The curve exhibiting a maximum value of R(t @ 0.4s) @ 21.6 kN corresponds to the dynamic simulation, where inertia forces have been accounted for, while the curve which monotonically decreases from an initial maximum value of R(t = 0) @ 42 kN is associated with a quasi-static simulation, where the acceleration forces have been neglected. The fact that, in the latter simulation, the initial value of the drag force is not equal to zero, comes from the quasi-static approximation, since in such a case, the velocity field attached to the initial configuration (and hence the drag force distribution: see Fig. 11) is also nonzero. Fig. 12 shows that, even though both curves become almost coincident after an elapsed time of 1.5 s, they are quite significantly different in the first stage of propagation. It turns out in particular that, owing to such acceleration forces, the maximum peak value of the resultant drag force is significantly reduced (almost cut in half in the present example). 7. Concluding remarks A general procedure based upon the numerical implementation of a variational principle for the velocities, combined with a time-step integration algorithm, has been proposed in order to predict the propagation of a liquefied soil mass, modelled as a Bingham material, even in the case when large displacements are involved. The effect of such a moving liquefied soil on structures is then investigated

using the same numerical approach (decomposition– co-ordination method), making it thus possible for instance to evaluate the drag force applied by the liquefied soil on a pile foundation modelled as a fixed obstacle. A fully integrated computational tool dealing with both aspects simultaneously has been set up in the simplified situation when the evolution problem can be solved independently from the evaluation of the drag forces applied to the structures. Its feasibility has been demonstrated on an illustrative example, and preliminary comments have been made on the relative influence of some key parameters such as the liquefied soil or soil-substratum shear strengths. Looking forward to devising an engineering design-oriented numerical code, further improvements may be considered in the near future:  The generalization of the method to three dimensional situations, which would probably require significant progress as regards finite element discretization techniques, although decisive simplifications, such as the ‘‘shallow water approximation’’, could be envisaged in the frequently encountered situation when the liquefied soil flow is predominantly parallel to a plane.  The use of remeshing techniques such as those developed in the quite different context of the simulation metal-forming processes (see for instance [3]) where large strains are also involved, which may induce excessive mesh distortion and thus affect the numerical accuracy of the computations.  Extending the approach to the case when the evolution problem can no longer be disconnected from the evaluation of the forces exerted on structures. Such a situation arises when the characteristic size of the latter is not sufficiently small with respect to the overall dimensions of the propagating liquefied soil mass. It should however be kept in mind that the validity of the whole procedure relies upon the assumption that both the viscosity and shear strength parameters governing the liquefied soil, which are captured through the adoption of a Bingham model, remain constant and unaffected by the propagation. A procedure could however be imagined in

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

which for instance, due to the dissipation of the excess pore pressures as time passes, the soil may be recovering its strength, which could be modelled by a corresponding increase with time of parameter k and possibly l in the Bingham model. The precise evaluation of these parameters by appropriate in situ or laboratory experimental tests remains a difficult task for the geotechnical engineer. A possible line of research, which is currently in progress, consists in implementing the numerical procedure developed in this paper for simulating the response of liquefied sand in the free field, and compare the velocity and displacement profiles soobtained with experimental results drawn from centrifuge model tests on earthquake-induced lateral spreading using a laminar box [36]. This would serve the double purpose of validating the proposed numerical scheme and related computer code on the one hand, assessing realistic values to the key constitutive parameters of liquefied soil (k and l) from a back calculation procedure, on the other hand. References [1] Abdoun T, Dobry R, O’Rourke TD, Goh SH. Pile response to lateral spreads: centrifuge modeling. J Geotech Geoenviron Eng 2003;129(10):869–78. [2] Adachi K, Yoshioka N. On creeping flow of a visco-plastic fluid past a circular cylinder. Chem Eng Sci 1973;28:215–26. [3] Aymone JLF, Bittencourt E, Creus GJ. Simulation of 3D metalforming using an arbitrary Lagrangian–Eulerian finite element method. J Mater Process Technol 2001;110:218–32. [4] Bartlet SF, Youd TL, Empirical prediction of lateral displacement. In: Proceedings of the 4th US–Japan workshop on earthquake resistant design of lifeline facilities and countermeasures against soil liquefaction, vol. 1, Technical Report, NCEER-92-0019, Honolulu, USA; 1992. p. 65–351. [5] Boulanger RW, Wilson WW, Kutter BL, Soil–pipe-superstructure interaction in liquefiable sand, Transportation Research Record 1569, Transportation Research Board, Washington (DC); 1997. [6] Chen H, Lee CF. Runout analysis of slurry flows with Bingham model. J Geotech Geoenviron Eng 2002;128(12):1032–42. [7] Coussot PH. Steady, laminar flow of concentrated mud suspension in open channel. J Hydraul Res 1994;32(4):535–59. [8] Dobry R, Abdoun T, O’Rourke TD, Goh SH. Single piles in lateral spreads: field bending moment evaluation. J Geotech Geoenviron Eng 2003;129(10):879–89. [9] Fallgren RB, Smith JL , Ground displacement at San Fernando valley juvenile hall during San Fernando earthquake. In: San Fernando California earthquake of February 9, 1971, vol. 3, US Dept. of Commerce, NOAA, Washington (DC); 1973. p. 189–196. [10] Faxe`n OH. Forces exerted on a rigid cylinder in a viscous fluid between two parallel fixed planes. Proc Roy Swedish Acad Eng Sci 1946;187:1–13. [11] Finn WDL. State-of-the-art of geotechnical earthquake engineering practice. Soil Dynam Earthquake Eng 2000;20:1–15. [12] Fortin M, Glowinski R. Augmented Lagrangian methods. Amsterdam: North-Holland; 1983. [13] Fortin A, Coˆte´ D, Tanguy PA. On the imposition of friction boundary conditions for the numerical simulations of Bingham fluid flows. Comput Meth Appl Mech Eng 1991;88:97–109. [14] Glowinski R, Le Tallec P. Augmented Lagrangian and operator splitting method in non-linear mechanics. SIAM Stud Appl Math 1989. [15] Hadush S, Yashima A, Uzuoka R. Importance of viscous fluid characteristics in liquefaction induced lateral spreading analysis. Comput Geotech 2000;27:199–224.

119

[16] Hadush S, Yashima A, Uzuoka R, Moriguchi S, Sawada K. Liquefaction induced lateral spread analysis using the CIP method. Comput Geotech 2001;28:549–74. [17] Hamada M, Towhata I, Yasuda S, Isoyama R. Study on permanent ground displacement induced by seismic liquefaction. Comput Geotech 1987;4:197–220. [18] Hild P, Ionescu IR, Lachand-Robert T, Rosca I. The blocking of an inhomogeneous Bingham fluid. Applications to landslides. Math Modell Numer Anal 2002;36(6):1013–26. [19] Imamura S, Hagiwara T, Tsukamoto Y, Ishihara K. Response of pile groups against seismically induced lateral flow in centrifuge model tests. Soils Foundat 2004;44(3):39–55. [20] Ishihara K, Cubrinovski M. Soil–pile interaction in liquefied deposits undergoing lateral spreading. In: Maric B et al., editors. Geotechnical hazards. Rotterdam: Balkema; 1998. p. 51–64. [21] Kawakami T, Suemasa N, Hamada M, Sato H, Takada T, Experimental study on mechanical properties of liquefied sand. In: Proceedings of the 5th US–Japan workshop on earthquake resistant design of lifeline facilities and countermeasures against soil liquefaction, vol. 1, Technical Report, NCEER-94-0026, Salt Lake City, USA; 1994. p. 99–285. [22] Liu KF, Mei CC. Slow spreading of a sheet of Bingham fluid on an inclined plane. J Fluid Mech 1989;207:505–29. [23] Mitsoulis E. On creeping drag flow of a viscoplastic fluid past a circular cylinder: wall effects. Chem Eng Sci 2004;59:789–800. [24] Montassar S, Contribution to the numerical simulation of liquefied ground flow and its effects on structures. PhD thesis (in French), ENPC, Paris; 2005. [25] Montassar S, de Buhan P. Minimum principle and related numerical scheme for simulating initial flow and subsequent propagation of liquefied ground. Int J Numer Anal Meth Geomech 2005;29:1065–86. [26] Montassar S, de Buhan P. Some general results on the stability and flow failure of rigid viscoplastic structures. Mech Res Commun 2006;33:63–71. [27] Murff JD, Wagner DA, Randolph MF. Pipe penetration in cohesive soil. Ge´otechnique 1989;39(2):213–29. [28] Orense RP, Towhatal I. Three dimensional analysis on lateral displacement of liquefied subsoil. Soils Foundat 1998;38(4):1–15. [29] Pastor M, Quecedo M, Fernandez Merodo JA, Herreros MI, Gonzalez E, Mira P. Modelling tailing dams and mine waste dumps failures. Ge´otechnique 2002;52(8):579–91. [30] Randolph MF, Houlsby GT. The limiting pressure on a circular pipe loaded laterally in cohesive soil. Ge´otechnique 1984;34(4):613–23. [31] Roquet N, Saramito P. An adaptive finite element method for Bingham fluid flows around a cylinder. Comput Meth Appl Mech Eng 2003;192:3317–41. [32] Salenc¸on J. An introduction to the yield design theory and its applications to soil mechanics. Eur J Mech, A/Solids 1990;9(5): 477–500. [33] Salenc¸on J. Handbook of continuum mechanics. Berlin: Springer; 2000. [34] Salenc¸on J, Indentation of a cohesive soil by a circular pipe. In: Proceedings 15th international congress on soil mechanics and foundation engineering, Istanbul; 2001. p. 1311–4. [35] Seed HB, Harder J, SPT-based analysis of cyclic pore pressure generation and undrained residual strength. In: Duncan JM, editor. Proceedings H. Bolton Seed memorial symposium, vol. 2, University of California, Berkeley, May 9–11; 1990. p. 351–76. [36] Taboada V, Dobry R. Centrifuge modelling of earthquake-induced lateral spreading in sand. J Geotech Geoenviron Eng 1998;124(12): 1195–206. [37] Tamate S, Towhata I. Numerical simulation of ground flow caused by seismic liquefaction. Soil Dynam Earthquake Eng 1999;18: 473–85. [38] Towhata I, Al-Hussaini TM. Lateral loads on offshore structures exerted by submarine mudflows. Soils Foundat 1988;32:26–34. [39] Towhata I, Orense RP, Toyota H. Mathematical principles in prediction of lateral ground displacement induced by seismic liquefaction. Soils Foundat 1999;39(2):1–19.

120

S. Montassar, P. de Buhan / Computers and Geotechnics 33 (2006) 108–120

[40] Towhata I, Sasaki Y, Tokida KI, Matsumoto H, Tamari Y, Yamada K. Mechanism of permanent displacement of ground caused by seismic liquefaction. Soils Foundat 1992;32(3):79–96. [41] Towhata I, Toyota H, Vargas-Monge W, Dynamics in lateral flow of liquefied ground. In: Proceedings of the 10th Asian regional conference on SFME, vol. 1, Beijing, 1995. p. 497–500. [42] Uzuoka R, Yashima A, Kawakami T, Konrad J-M. Fluid dynamics based prediction of liquefaction induced lateral spreading. Comput Geotech 1998;32(3–4):243–82.

[43] Wilson DW, Boulanger RW, Kutter BL. Observed seismic lateral resistance of liquefying sand. J Geotech Geoenviron Eng 2000;126(12):898–906. [44] Yasuda S, Kiku S, Yosida T, Shaking table tests on the friction force between the buried pipelines and the ground during liquefaction. In: Proceedings of the 24th national conference JSSMFE (in Japanese); 1989. p. 759–60. [45] Zisis TH, Mitsoulis E. Viscoplastic flow around a cylinder kept between parallel plates. J Non-Newtonian Fluid Mech 2002;105:1–20.