Dual-phase coupled u–U analysis of wave propagation in saturated porous media using a commercial code

Computers and Geotechnics 55 (2014) 316–329

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Dual-phase coupled u–U analysis of wave propagation in saturated porous media using a commercial code Feijian Ye ⇑, Siang Huat Goh, Fook Hou Lee Department of Civil and Environmental Engineering, National University of Singapore, Block E1A, #07-03, No. 1 Engineering Drive 2, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 4 July 2013 Received in revised form 28 August 2013 Accepted 5 September 2013 Available online 18 October 2013 Keywords: Saturated Porous media Wave propagation Dual-phase Fully-coupled Finite element

a b s t r a c t This paper presents a method for incorporating dual-phase u–U wave propagation analysis for saturated porous media in commercial codes. Biot’s formulation is first re-formulated in terms of effective stress and pore pressure, followed by its finite element discretization. The method of incorporating the dualphase computation involves representing the solid and pore fluid phases by two overlapping meshes with collocated elements and nodal points. Viscous coupling between the solid and fluid phases are realized by connecting each pair of collocated nodal points using viscous Cartesian connectors. Volumetric compatibility between the solid and fluid phases is enforced by coding the behaviour of the two phases and their volumetric interaction within a user-defined material subroutine. Non-linear behaviour of the soil skeleton can also be prescribed within the same subroutine. Comparison with existing analytical or numerical solutions for three one-dimensional elastic problems shows remarkably good agreement. Finally, a three-dimensional elasto-plastic example involving impulsive surface loading on dry sand overlying saturated sand, which is akin to a dynamic compaction problem, is analyzed to demonstrate the potential application of the u–U analysis on a more realistic transient loading problem. The results of the analysis highlight the importance of the slow wave in dynamic compaction process. Although this approach was developed on the ABAQUS Explicit platform, it is, in principle, applicable to other existing codes with the right features. Its main advantage is that users can still access the computational functionalities and stability, as well as pre- and post-processing features that often come with commercial codes. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Many soil dynamic problems such as blast loading and heavy tamping in saturated soils involve wave propagation through solid and fluid phases. Where the loading rate is much higher than the rate of drainage, an undrained condition is often assumed to prevail, in which there is no movement of pore fluid relative to the soil skeleton, so that the two phases are considered as a single phase [1]. However, this assumption may not be valid for all loading and ground conditions. Biot [2,3] showed that, even when there is no drainage of pore fluid, relative oscillatory motion can still exist between soil skeleton and pore fluid, and different dilatational waves propagate through pore fluid and soil skeleton at different speeds. In partially drained conditions, relative movement between soil skeleton and water will invariably occur. Hence, the soil skeleton and pore fluid is more accurately considered as two coupled phases in highly transient loading events. Truesdell’s Mixture Theory [4,5] and Biot’s dual-phase formulations [2,3] are commonly used to describe solid, fluid interaction under dynamic

⇑ Corresponding author. E-mail address: [email protected] (F. Ye). 0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.09.002

conditions. Of these two, Biot’s formulation is the more widely used. Biot’s formulation considers the soil skeleton and pore fluid to be two phases coupled together by inertial, viscous and volumetric coupling. The primary variables of this u–U or u–U–p formulation are the displacement of the soil skeleton (u), displacement of the pore fluid (U) and partial fluid pressure (p). For problems involving lower loading rates, Zienkiewicz et al. [1,6] proposed that the inertial effect due to the relative acceleration between solid and pore fluid phases can be neglected. This leads to the u–p formulation, in which the primary variables are the displacement of the soil skeleton (u) and the pore pressure (p). The u–p formulation is sufficient for seismic problems; however, Zienkiewicz et al. [6] pointed out that it may be less accurate for high frequency, short duration problems. For such problems, a fully dynamic formulation, such as the u–U framework considered in this paper, is preferable. Many dynamic soil–pore water interaction problems involving transient loading are still not well-studied. An example of this is the propagation of ground waves, such as that induced by explosion [7–9] or dynamic compaction [10–12], through full or partially saturated soils. Another example is excess pore pressure

317

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

induced by breaking waves in rubble mound breakwaters [13–16] which may generate transient high wave loads on structures. Many of the numerical studies associated with explosion [7–9] used onephase concept whereby the effect of fluid flow was not considered. Only [12] carried out a fully-coupled analysis of dynamic compaction using u–p formulation. Existing analytical solutions for Biot’s dual-phase formulation [17–21] are restricted to highly idealized one-dimensional problems of wave propagation through homogeneous linear elastic medium. Published numerical solutions are also largely limited to linear elastic behaviour and one dimensional idealization [22– 27]. There are relatively few existing u–U codes and these, e.g. POROUS [17] and DYNAFLOW [28], tend to be codes developed in-house. Commercial codes, such as ABAQUS and LS-DYNA, are increasingly used in engineering analysis because they have well-developed 2- and 3-dimensional capabilities, portals for user-defined constitutive models, user interfaces and visualization tools. However, they do not have facilities for u–U analyses. This is possibly one reason for the dearth of u–U analyses. Instead of developing a purpose-built u–U code, this paper proposes that u–U analysis can be implemented by incorporating appropriate user-defined subroutines into a commercial code, in this case, ABAQUS. By so doing, fully-coupled u–U dual phase analysis can be realized while retaining the functionality, user-interfaces and visualization tools of ABAQUS. This can potentially bring u–U dual phase analysis within the reach of many more engineers and researchers. This dual-phase coupled (DPC) approach is implemented on ABAQUS, but it can, in principle, be implemented on other commercial codes with similar user-defined portals. In the discussion below, Biot’s u–U formulation [2] is first summarized in terms of effective stress and pore pressure, following which its implementation using ABAQUS Explicit is presented. This will be followed by three validation examples involving onedimensional wave propagation in saturated soil. Finally, one three-dimensional, non-linear, transient loading example is presented. 2. Formulations

where K0 and Ks are the bulk moduli of the soil skeleton and soil grain, respectively. In most geotechnical problems, K0  Ks, so that a  1. Substituting Eq. (2) into Eq. (1), and assuming q12 to be negligible (to be discussed later), leads to Biot’s u–U formulation in terms of effective stress and pore pressure, hereafter termed u–U– p formulation,

r0ij;j  ð1  nÞp;j þ cðU_ i  u_ i Þ  ð1  nÞqs u€i ¼ 0

ð4aÞ

€i ¼ 0 np;j  cðU_ i  u_ i Þ  nqf U

ð4bÞ

in which qs and qf are the mass densities of the solid grain and pore fluid respectively. The viscous damping coefficient c is related to the Darcy’s coefficient of permeability k and unit weight of the pore fluid cw by

c¼

n2 cw k

The mass coefficients are related to the mass density of fluid and solid grains by the relations

q11 þ q12 ¼ ð1  nÞqs

ð6aÞ

q22 þ q12 ¼ nqf

ð6bÞ

The coefficient q12 is an inertial coupling coefficient which relates the stresses of one phase to the acceleration of the other. We can alternatively express the inertial interactions between the soil skeleton and fluid in Eq. (1) as

rsij;j þ cðU_ i  u_ i Þ  ðq11 þ q12 Þu€i  q12 ðU€ i  u€i Þ ¼ 0

ð7aÞ

p;j  cðU_ i  u_ i Þ  ðq22 þ q12 ÞU€ i þ q12 ðU€ i  u€i Þ ¼ 0

ð7bÞ

If u = U, then two independent terms, namely (q11 + q12) and (q22 + q12), are sufficient to describe the inertial effects on the respective phases, and q12 on its own, which represents the advective inertial coupling arising from the relative fluid–solid motion, is not needed. The value of q12 is not readily determined and is often assumed to be zero (e.g. [17,18,23]). The pore pressure is related to the strain tensor via

p ¼ K f

2.1. Biot’s u–U/u–U–p formulation The original Biot’s u–U/u–U–p formulation can be written in the form

rsij;j þ cðU_ i  u_ i Þ  q11 u€i  q12 U€ i ¼ 0

ð1aÞ

p;j  cðU_ i  u_ i Þ  q12 u€i  q22 U€ i ¼ 0

ð1bÞ

where r and p are the partial solid stress and partial fluid pressure, respectively, with positive stresses denoting tension; u and U are the solid and fluid displacement at time t; the single and double overdots denoting first and second derivatives with respect to time, respectively; c is the viscous damping coefficient; q11, q12 and q22 are mass coefficients which relate the inertial forces to the acceleration in the solid and fluid phases. The partial solid stress and fluid pressure are the average stresses on the solid and fluid phases respectively, and can be related to the effective stress r0ij and pore pressure p via the relations s ij

ð5Þ

ða  nÞekk þ nfkk K

n þ ð1  nÞ K fs þ

ð8Þ

K0 Kf K 2s

in which Kf is the bulk modulus of the fluid, ekk is the volumetric strain of the solid phase, and fkk is the volumetric strain of the fluid phase, denoted as

fkk ¼ U i;i

ð9Þ 0

In most saturated soils, K  Kf  Ks, and Eq. (8) reduces to

p¼

Kf ½ekk  nðekk  fkk Þ n

ð10Þ

If we further consider a completely ‘‘undrained’’ situation wherein the solid phase displacement and strain are the same as that of the pore fluid phase, then ekk = fkk and Eq. (10) reduces to the standard form for pore pressure in an undrained situation, that is

p¼

Kf n

ekk

ð11Þ

rij ¼ r0ij  dij ap ¼ rsij þ dij p

ð2aÞ

2.2. Finite element discretization and time stepping

p ¼ p=n

ð2bÞ

Following [6], application of the Galerkin method to Eq. (4) leads to the finite element form

in which rij is the total stress, p the pore pressure (positive in compression), n the porosity, dij the Kronecker delta and

a¼1

K0 Ks

ð3Þ



Ks þ Kd Ksf Ksf

Kf

" _ #  " € # " ð1Þ #   Ms 0 u u f þ þ _ € ¼ ð2Þ H H 0 Mf U U U f

   u



H H

ð12Þ

318

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

 and U are the nodal displacements of the solid and fluid in which u phases and the matrices are given by

Z

Ks ¼

BT D0 BdX

ð13Þ

X

Kd ¼

ð1  nÞ2 Kf n Z

Kf ¼ nK f

Z

RT RdX

ð14Þ

X

RT RdX

ð15Þ

X

Ksf ¼ ð1  nÞK f

Z

RT RdX

ð16Þ

NT NdX

ð17Þ

X

Ms ¼ ð1  nÞqs

Mf ¼ nqf

H¼c

Z

Z

Z X

NT NdX

ð18Þ

X

NT NdX

ð19Þ

X

f

ð1Þ

¼

Z

Z

NT tdC þ ð1  nÞ C

f

ð2Þ

¼n

C

Z C

NT qdC þ nqf

Z

NT qdC þ ð1  nÞqs

Z

NT b s d X

ð20Þ

X

NT bf dX

ð21Þ

X

where N is the shape function matrix; B the strain–displacement 0 matrix; R = rTN; D the stress–strain matrix; bs and bf the body force vectors for solid and fluid phases, respectively; t and q the surface tractions on the solid and fluid phases respectively. Eq. (12) can be expressed in the form

€ þ Cu_ þ Ku ¼ F Mu

3.2. Viscous coupling Solid–fluid coupling terms are represented in Eq. (12) by the off-diagonal sub-matrices H and Ksf. The former term represents viscous coupling while the latter represents constitutive coupling. The viscous coupling between the solid and fluid meshes specified in Eq. (19) is implemented in ABAQUS using Cartesian connector elements between the collocated nodes. ABAQUS also has an axial connector element, but this is restricted to one-dimensional analysis only. The expansion of Eq. (19) indicates that viscous damping coefficient c has to be scaled according to the element size. The scaled viscous damping coefficient c of the connector elements can be computed using the following scaling law

c ¼ gXc

ð25Þ

where X is the element volume (area) for three- (two-)dimensional elements and g is a factor that prescribes how the total viscous coupling within the element is to be distributed to the various nodes. Two approaches for discretization the viscous coupling are possible, namely a consistent damping scheme and a lumped damping scheme. The consistent damping scheme is based on direct integration of Eq. (19). However, it is not readily implemented using connectors elements since the coupling between each pair of solid and fluid nodes depends on the relative velocity of all pairs with different distribution factors g. The lumped damping discretization is more heuristic and only considers the damping force between collocated solid and fluid nodes as illustrated in Fig. 1. This is the method used herein. The contribution of one node of one phase to all non-collocated nodes of the other phase is assumed to be zero. The total viscous damping is then distributed evenly amongst the nodal points, thereby leading to g-values as shown in Table 1. This is similar to the lumped mass discretization scheme commonly used in explicit time integration. For the four-noded quadrilateral, the lumped damping matrix H is given by

ð22Þ

where

 M¼

Ms 0

    H H Ks þ Kd ;C ¼ ;K ¼ Mf H H Ksf 0

" #  ð1Þ f ; F ¼ ð2Þ Kf f

Ksf

ð23Þ  u¼

( ) ( )  € _ u u €¼ € ; u_ ¼ _ ; u U U U  u

ð24Þ

which can be solved explicitly or implicitly. In ABAQUS Explicit, a central difference rule is used.

Fluid mesh

Solid mesh

Fluid node

Solid node 3. Implementation on ABAQUS

3D view

3.1. Overlapping meshes In purpose-built dual-phase finite element codes [22–27], each nodal point has six translational degrees-of-freedom, three for the solid phase and three for the pore fluid phase. As such special elements are unavailable in ABAQUS, the DPC method employs two identical, overlapping meshes of completely collocated elements, both of which are analyzed concurrently. One mesh represents the soil skeleton and the other represents the pore fluid. The deformation of the solid mesh reflects the deformation of the porous skeleton while the deformation of the pore fluid mesh reflects the pore fluid displacement.

Fig. 1. Illustration of how the connector elements work in the lumped damping scheme (2D four-node elements).

Table 1 The scaling law for typical elements corresponding to lumped damping scheme. Element type

X

g

2D 2D 3D 3D

A A V V

1/3 1/4 1/4 1/8

three-node triangle four-node quadrilateral four-node tetrahedral eight-node hexahedral

319

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

2

axis of symmetry

7 0 0 0 0 0 07 7 0 14 0 0 0 0 0 7 7 7 0 0 14 0 0 0 0 7 7 0 0 0 14 0 0 0 7 7 7 0 0 0 0 14 0 0 7 7 0 0 0 0 0 14 0 7 5 0 0 0 0 0 0 0 14

(z3,r3)

4

3 (r 0,z0)

1

1 4

For two-dimensional axisymmetric elements, the H-matrix of a four-noded quadrilateral can be represented using cylindrical coordinates as

2

H¼

Z

NT cNdX ¼ 2pc

ε kk

NT Nrdrdz

where r can be expressed by using the shape functions.

Fig. 2. The axisymmetric four-node quadrilateral element.

Common Variable

Z A

X

Given σ ijn and d ε ij Compute trial stress e dσ ijT = Dijkl d ε kl

σ ijT = σ ijn + dσ ijT

Solid Subroutine

If f (σ ijT ) ≥ 0

No

σ ijn +1 = σ ijT

Yes Stress return algorithm ep Compute Dijkl ( Dep )

ep dσ ijn +1 = Dijkl d ε kl

σ ijn +1 = σ ijn + dσ ijn +1 σ ijn +1

Given p n and d ζ ij Fluid Subroutine

ð26Þ

r

(z2,r2)

(z1,r1)

0 0 0 0 0 0 0

6 60 6 60 6 6 60 H¼Ac ~ 6 60 6 6 60 6 60 4

z (z4,r4)

1 4

3

Common Variable

ζ kk

Compute trial stress Q R dp = d ε kk + d ζ kk n n p n +1 = p n + dp

p n +1

Fig. 3. Flow chart showing the DPC implementation for a dual-phase elasto-plastic model. Dep represents the general elasto-plastic stress–strain matrix.

ð27Þ

320

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

Fig. 4. DPC numerical model for Example 1.

1 ½ðr 1 þ r 2 þ r 3 þ r 4 Þ þ 4an ¼ r 0 þ an 4

ð28Þ

1.2

Wave I

where ri (i = 1,2,3,4) represents the distance of ith node from the axis of symmetry, a the half-length of the element in the radial direction, n the local radial coordinate and r0 ¼ 14 ðr 1 þ r 2 þ r3 þ r 4 Þ is the average radial distance of the nodes from the axis of symmetry. For the four-noded quadrilateral element shown in Fig. 2, the axisymmetric lumped damping matrix can be written as 3 3r 0  a 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 7 3r 0  a 7 6 7 6 6 0 0 0 0 0 0 7 0 3r 0 þ a 7 6 6 0 0 0 0 7 0 0 3r 0 þ a pAc 6 0 7 H¼ 7 6 6 6 0 0 0 0 7 0 0 0 3r 0 þ a 7 6 6 0 0 0 7 0 0 0 0 3r 0 þ a 7 6 7 6 4 0 0 5 0 0 0 0 0 3r 0  a 0 0 0 0 0 0 0 3r 0  a ð29Þ 2

1

Velocity (cm/s)

r¼

0.8

Wave II

0.6

0.4

Analytical (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

0.2

0 0

20

40

where the distribution coefficients for the nearer nodes and the further nodes are, respectively:

g1;4 ¼

60

80

100

120

Time (μs)

(a) Solid 1.2

3r 0  a 6

ð30aÞ

Wave II

g2;3

ð30bÞ

3.3. Constitutive coupling As Eq. (12) shows, the combined stiffness matrix K is needed to evaluate the internal forces, where

 K¼

Ks þ Kd

Ksf

Ksf

Kf

Velocity (cm/s)

1

3r 0 þ a ¼ 6

0.8

0.6

Wave I

0.4

 ð31Þ

The sub-matrices Ks and Kf are the stiffness matrices of the solid skeleton and fluid phases, respectively, and can be evaluated in standard fashion. Kd represents the effect of pore pressure gradient on the soil skeleton, that is the term (1  n)p,j in Eq. (4a). Non-linear soil behaviour can be reflected through the D0 matrix in Eq. (13). The pore fluid phase is assumed to be linear elastic. The coupling sub-matrix Ksf accounts for the volumetric compatibility

Analytical (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

0.2

0 0

20

40

60

80

100

120

Time (μs)

(b) Fluid Fig. 5. Solid and fluid nodal velocity history at the monitoring point for low drag.

321

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

1.2

1.2

Wave I 1

0.8

Wave II

Velocity (cm /s)

Velocity (cm/s)

Analytical (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

1

0.6

0.4

Numerical Inv. (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

0.2

0.8

0.6

0.4

0.2

0

0 0

20

40

60

80

100

22

120

23

24

25

26

Time (μs)

(a) Solid

28

29

30

31

32

33

29

30

31

32

33

(a) Solid

1.2

1.2

1

1

Velocity (cm /s)

Velocity (cm/s)

27

Time (μs)

0.8

Wave I 0.6

Wave II

Analytical (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

0.8

0.6

0.4

0.4

Numerical Inv. (Garg et al. 1974) POROUS (Garg et al. 1974) Qiu & Fox (2008) DPC

0.2

0.2

0

0 0

20

40

60

80

100

22

120

23

24

25

26

27

28

Time (μs)

Time (μs)

(b) Fluid

(b) Fluid

Fig. 6. Solid and fluid nodal velocity history at the monitoring point for medium drag.

Fig. 7. Solid and fluid nodal velocity history at the monitoring point for high drag.

1.5

1

0.5

Amplitude

condition which must exist between the two phases. This condition can be illustrated by considering a simple example in which  the fluid nodal displacement U is zero, while the soil skeleton is compressed. The resulting volumetric compression of the soil skeleton will lead to a reduction in the volume of the pores, thereby compressing the pore fluid and causing a rise in pore pressure, even though there is no fluid nodal displacement. ABAQUS, as well as many commercial finite element softwares, has no built-in feature to reflect this volumetric coupling. For this reason, constitutive coupling is incorporated into ABAQUS in the form of a user-defined material through the user material interface VUMAT. ABAQUS Explicit only allows one user-defined material subroutine to be incorporated, which is invoked once for each integration point. ABAQUS Explicit passes the strain tensor for each integration point into VUMAT. VUMAT returns the solid and fluid stress tensors to ABAQUS Explicit, which uses it to compute the stiffness-induced force vector. VUMAT computes the stress for both the solid and fluid phases by calling two subroutines, one for each phase. A controlling routine within VUMAT directs the computation to the respective subroutines depending upon whether the element is solid or fluid. This is accomplished via

0 -10

-0.5

-1

-1.5

0

10 10

20 20

30 30

4040

5050

Normalized Time (τ = t/ρk) Step Sinusoidal Spike Fig. 8. The amplitude of three loadings in Example 2.

the material designation parameter ‘‘cmname’’. For each pair of collocated element, VUMAT is written to loop through the solid phase element first, followed by the fluid phase element. The volumetric strains are stored within VUMAT and its subroutines as global

322

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

Table 2 Summary of material properties for Example 2 (self-consistent unit). Porosity Young’s modulus of skeleton Poisson’s ratio Density of solid Density of fluid Density of bulk mixture Darcy’s coefficient of permeability Magnitude of stress loading Material constant Material constant

n E

t qs qf q k

r0 a M

0.333 3000 0.2 0.3102 0.2977 0.3060 0.004883 100 0.667 13848

variable through the COMMON declaration. The procedure is illustrated in Fig. 3 for elasto-plastic soil behaviour. 3.4. Mesh generation and data input The following procedures are used in the preparation of the finite element model: (a) The nodal and element information for the overlapping solid and fluid meshes are first generated. (b) A MATLAB-based program generates the connector damping elements. This program reads in the nodal and element information from the input file previously generated.

Fig. 10. Solid and fluid nodal displacement history at the surface for the sinusoidal loading.

(c) A new input file following the ABAQUS format is then generated which combines the mesh information from (a), together with the connector damping elements generated in (b). The scaled damping coefficients and the material parameters are also added to the file. (d) The numerical model is then sent into the solver, together with the corresponding user subroutine.

4. Verification examples In this section, the DPC method is compared with available analytical or numerical solutions for three transient dynamic problems.

4.1. Example 1: Fully saturated porous half space subject to step velocity loading

Fig. 9. Solid and fluid nodal displacement history at the surface for the step loading.

Garg et al. [17] presented the analytical and numerical solutions for the problem of a transient step velocity input of 1 cm/s applied uniformly across the entire surface of a fluid saturated, porous, elastic half-space. Qiu and Fox [25] repeated the same problem using their approach. The problem can be simplified into a laterally-constrained, infinitely long, plane strain, saturated column, Fig. 4. The length of the column is set long enough so that the reflected waves will not interfere with the capturing of the incident ones.

323

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

1.2

0 -2

1

Amplitude

DPC Analytical

-4 -6

0.8 0.6

Load I H I (t )

0.4

Load II H II (t )

-8 0.2

-10

0 -200

-12

0

200

400

600

800

1000

Time (μs) 0

10

20

30

40

50 Fig. 12. Load profiles for Example 3.

(a) Solid (u) 3.5 3

DPC Analytical

2.5 2 1.5 1 0.5 0 0

10

20

30

40

50

(b) Fluid respect to solid (w) Fig. 11. Solid and fluid nodal displacement history at the surface for the spike loading.

The material properties used are the same as those of Garg et al. and will not be repeated herein. Following Garg et al., three values of viscous drag are analyzed with c values of 2.19  104, 2.19  106 and 2.19  1010 kg/m3 s; termed hereafter as low, medium and high drag cases. Garg et al. assumed solid–fluid inertial coupling to be zero. It should be noted that Garg et al.’s parameters are not necessarily realistic for soils. They are only used here for purpose of verification. Both automatic time step and user-specified time step of 106 s is used for this example. The waves are monitored at a point located at a depth of 10 cm beneath the soil surface, and the results are plotted for the solid and fluid particle velocity at this location. Fig. 5 shows the solid and fluid particle velocity time histories at the monitoring point for the low drag case. The presence of two dilatational waves can be seen clearly in the fluid velocity response of Fig. 5b. The fast wave (wave I) arrives at about 28 ls, leading to an almost instantaneous increase in the velocity to 0.5 cm/s. The subsequent arrival of the slow wave (wave II) at about 78 ls doubles the fluid velocity to about 1 cm/s, which is the prescribed value at the boundary. The solid phase velocity, on the other hand, achieves a value of about 1 cm/s upon the arrival of wave I, while the influence of wave II is weakly discernable (Fig. 5a) and not as prominent as in the fluid phase. It should be pointed out that the solid and fluid velocity trends shown in Fig. 5 are significantly influenced by the material parameters adopted in the analysis, which, as was pointed out earlier, are not necessarily realistic for soils in this example.

The dual-wave structure is also evident in the medium-drag cases, Fig. 6. Comparing Figs. 6a and 5a, it is noted that there is no significant difference in the solid velocity response between the low and medium drag cases. However, Fig. 6b shows that the increased drag between the two phases causes a gradual increase of the fluid velocity in the interval between the arrivals of waves I and II, in contrast with the stepped response of Fig. 5b. For the high drag case of Fig. 7 (plotted for a smaller time interval), the drag forces are high enough to constrain the fluid and solid particles to move virtually in unison; hence the dual-wave structure is no longer apparent and Fig. 7a and b appears to be almost identical. Garg et al.’s numerical analysis was conducted using an inhouse finite difference code, POROUS. As Figs. 5–7 show, neither the DPC, POROUS nor Qiu & Fox’s results can completely capture the instantaneous increase in velocity that the analytical results show; with all showing some dispersion of the wavefront. This dispersion is influenced by the element sizes according to the element size study which is not discussed herein. However, it should be noted that, in all cases, the rise times are very short, the longest being no more than about 20 ls. The agreement between DPC, POROUS and Qiu and Fox’s results is remarkably good for the low- and medium-drag cases, with the DPC showing slightly more overshoot on the rise. This overshoot appears to be accentuated for the high-drag case. Qiu and Fox’s results, on the other hand, exhibit less overshoot than DPC but more drift on the rise than DPC and POROUS. Although the time step effect is not shown here, using automatic time step control predicts almost the same results as user-specified time step control. 4.2. Example 2: Fully saturated porous half space subject to dynamic stress loading Simon et al. [18] presented analytical solutions for the problem of an infinitely thick, saturated, elastic soil layer subjected to three types of uniformly distributed stress loadings on the soil surface, with load-time history as shown in Fig. 8. A 25-m long, plane strain soil column, discretized into 5 cm-thick finite elements is used to model the soil layer. The surface stress loading was only imposed on the solid phase. The material properties are presented in Table 2 for a consistent set of units. In this Table, a is given in Eq. (3) and M is a material constant defined by

M¼

n Kf

1 þ aKn s

ð32Þ

The results, shown in Figs. 9–11, are expressed in normalized displacements u0 and w0 and normalized time s, such that

324

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

Fig. 13. Solid and fluid nodal velocity history at the monitoring point for low drag (load profile I).

u0 ¼ u

Vc

qr0

Vc w0 ¼ w kr0

s¼

t

ð33Þ

ð34Þ

ð35Þ

qk

where q, k, r0 are defined in Table 2, w is the relative displacement of fluid (U) with respect to solid (u), denoted as

w ¼ nðU  uÞ

ð36Þ

Vc is the wave speed of the solid and fluid in a situation where the particle motions of the two phases are in complete unison and is given by

Vc ¼

sffiffiffiffiffi Ec

q

ð37Þ

Fig. 14. Solid and fluid nodal velocity history at the monitoring point for high drag (load profile I).

in which

4 Ec ¼ K 0 þ G þ M 3

ð38Þ

is the equivalent constrained modulus and q is the bulk density. If we assume that the solid particles are incompressible, then Ks = 1, and Eq. (38) reduces to

Kf 4 Ec ¼ K 0 þ G þ 3 n

ð39Þ

For dry soil, Kf = 0, and Ec reduces to the constrained modulus of the soil skeleton. The numerical approach seems to be in reasonable agreement with Simon et al.’s results, Figs. 9–11, although, for this problem, the DPC algorithm seems to exhibit some drift. For all three loadtime histories, the solid nodes move downwards whereas the fluid nodes move in the opposite direction. This is the effect of the volumetric coupling. As the solid nodes move downwards, the soil skeleton is compressed and the resulting expulsion of fluid is

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

325

Fig. 15. Solid and fluid nodal velocity history at the monitoring point for low drag (load profile II).

Fig. 16. Solid and fluid nodal velocity history at the monitoring point for high drag (load profile II).

reflected in the upward movement of the fluid nodes. The reverse occurs when the solid nodes move upwards. The surface displacement trends of Figs. 9–11 are consistent with the applied stress loadings. The step loading gives rise to a monotonically increasing displacement and the sinusoidal loading gives rise to a sinusoidal one. The spike loading shows a slight rebound in the soil skeleton, this effect being more significant in the fluid phase. The rebound phenomenon is related to drag effects associated with the soil permeability. Although not discussed here, an additional analysis shows that, when drag effects are negligible, the rebound effect is not discernable and the solid and fluid displacements are relatively constant with time.

of 0.3 m beneath the surface. Using Laplace Transform inversion, Hiremath et al. [24] obtained analytical solutions to this problem subjected to two velocity profiles as shown in Fig. 12 where profile I is the step function and profile II follows the relation

4.3. Example 3: Fully saturated porous layer with a finite thickness subjected to velocity loading This example is similar to Example 1 except that the saturated porous layer has a finite thickness of 0.5 m and the bottom boundary is rigid and impermeable. The monitoring point is set at a depth

HII ðtÞ ¼ 1  0:2et=T

ð40Þ

where T = 140.9 ls. Each velocity profile is imposed on the surface nodes of the fluid and solid mesh. The same material properties as Example 1 given by Garg et al. [17] are used. Automatic time step size is used in the DPC analyses. Figs. 13 and 14 show the computed velocity histories at the monitoring point obtained using load profile I. Due to multiple reflection effects at the fixed and free ends of the saturated porous layer, the velocity response appears as a series of finite pulses. As in Example 1, the effect of wave II is more pronounced in the fluid phase in the low drag case, Fig. 13b. For the high drag case, Fig. 14 shows that the solid and fluid phases appear to move in unison. Similar observations may be made for the computed responses of Figs. 15 and 16, corresponding to load profile II.

326

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

2m 2m

Cushion P (MPa)

D = 20m

1.0

(0.001, 1)

0 0 .0 0 1

0.01

T (s)

1

T (s)

T=0.01s P (MPa) 1.0

(0.1, 1)

0 0 .1

T=1.0s

20m

Fig. 17. Simulation of dynamic compaction using a load history curve.

surfaces are well-captured. For the high-drag cases, some overshoot is evident for the DPC results; this is also consistent with the trend in Example 1. The arrivals of incident wave I and incident/reflected wave II are illustrated respectively in Figs. 13b and 15b, which reveal that the velocity responses at the recorded station are interfered by the recurrent incident and reflected wave I & II. Examples 1 and 3 illustrate the influence of the drag interaction between the fluid and solid phases. Under low to medium drag conditions, the solid and fluid phases exhibit distinct velocity responses, as shown in Figs. 5, 13 and 15. On the other hand, the high drag condition results in strong viscous coupling between the two phases, thus suppressing relative movement and constraining the two phases to virtually move in tandem, as shown in Figs. 7, 14 and 16.

Table 3 Material properties for water saturated soils (application example). Skeleton and pore fluid qs (kg/m3) qf (kg/m3) 2670 1000 E (MPa) t 50 0.3

Ks (GPa) 34.5 c (kPa) 10

Kf (GPa) 2.0 0 / 30°

0°

Dry soils qdry (kg/m3) 1735.5

t 0.3

c (kPa) 10

/ 30°

Edry (MPa) 50.0

n 0.35

k (m/s) 1.0  106

w

0

w 0°

qs is the Density of soil particles, qf the density of fluid, qdry the Density of dry soils, Ks the bulk modulus of soil particles, Kf the bulk modulus of fluid, n the porosity, k the coefficient of permeability, E the Young’s modulus of solid skeleton of saturated soils, Edry the Young’s modulus of solid skeleton of dry soils, t the Poisson’s ratio, c the cohesion, /0 the effective angle of friction, and w is the angle of dilation.

As shown in Figs. 13–16, there is overall good agreement between the DPC predictions and Hiremath et al.’s results. Good agreement is also achieved in a comparison with Qiu and Fox’s [25] results for low-drag case with profile I, as shown in Fig. 13b. In all four cases, the effects of wave reflection at the free and fixed

5. Application example In this section, the DPC algorithm is applied to a three-dimensional, quadrant-symmetric problem of uniformly distributed impulsive loading over a square patch of the ground surface. The

Z/D

0

0

-0.1

-0.1

-0.02

-0.2

-0.2

-0.04

-0.3

-0.3

-0.06

-0.4

-0.4

-0.08

-0.5

-0.5

-0.1

-0.6

-0.6

-0.12

-0.7

-0.7

-0.14

-0.8

-0.8

-0.16

-0.9 -1

-0.9

(a) 0

0.1

0.2

0.3

0.4

0.5

X/D

0.6

0.7

0.8

0.9

1

-1

-0.18

(b) 0

0.1

0.2

0.3

0.4

0.5

X/D

0.6

0.7

0.8

Fig. 18. Normalized (a) vertical effective stress and (b) pore pressure contours at time t = 10 ms.

0.9

1

-0.2

327

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

Z/D

0

0

-0.1

-0.1

-0.02

-0.2

-0.2

-0.04

-0.3

-0.3

-0.06

-0.4

-0.4

-0.08

-0.5

-0.5

-0.1

-0.6

-0.6

-0.12

-0.7

-0.7

-0.14

-0.8

-0.8

-0.16

-0.9 -1

-0.9

(a) 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

X/D

-1

-0.18

(b) 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

-0.2

X/D

Fig. 19. Normalized (a) vertical effective stress and (b) pore pressure contours at t = 16 ms.

soil profile consists of a thin layer of dry sand overlying saturated sand, Fig. 17. A rigid cushion is prescribed over the loaded area to distribute the impulsive load. The top boundary of the fluid mesh coincides with the groundwater table; the dry sand layer is modelled only by the soil skeleton. The contact between the cushion and the sand is assumed to be smooth. The sand is modeled using a Drucker–Prager model with non-associate flow rule. The material properties used are shown in Table 3. The analysis is conducted for two different coefficients of permeability, viz. 104 m/s and 102 m/s. These two values are representative of the typical range of permeability for sandy soils. The impulsive loading is represented by a triangular load-time history, with a sharp rise and a gradual decay, Fig. 17. Two loading durations are analyzed, namely 0.01 s and 1 s. The peak load P0 is 1 MPa for both cases. The speed of the dilatational wave in the dry sand can be readily shown to be approximately 197 m/s; this wave will be termed hereafter as the ‘‘dry solid wave’’. Fig. 18 shows the contours of the vertical effective stress and pore pressure, both normalized by the load peak, 10 ms after commencement of loading, for permeability of 104 m/s and loading duration of 0.01 s. In the analysis, negative pore pressure and stresses represents compression whereas positive represents dilatation. At this point of time, the dry solid wave has just traversed the entire thickness of the dry sand and reached the groundwater table. Fig. 18b shows the high pore pressure generated just below the groundwater table, corresponding to the initiation of the fast dilatational wave. Fig. 19 shows the corresponding contours 16 ms after commencement of loading. The velocities of the fast and slow dilatational waves in this problem are approximately 1580 m/s and 1 m/s. The pore pressure bulb at the mid-depth of the saturated layer is consistent with the arrival of the fast wave. On the other hand, the slow wave, which has a much lower velocity and is highly attenuated, is not evident from the effective stress contours in Fig. 19a. The sharp decrease in effective stress beneath the loaded area at the groundwater table and the rise in the pore pressure indicates that, at the groundwater table, the impulse from the dry solid wave is strongly converted to excess pore pressure, which propagates downwards via the fast wave. As Fig. 20a shows, apart from a very thin layer of saturated sand just beneath the groundwater table, the impulsive loading is not strongly reflected in effective stress fluctuations at subsequent time instances. Reflection of the effective stress wave at the groundwater is also evident from Fig. 20,

owing to the poor impedance coupling between the dry soil skeleton and pore water. Fig. 21 shows the corresponding effective stress and pore pressure profiles for the same permeability but with the much longer loading duration of 1 s. As can be seen, the slow wave in the saturated sand is now much better formed. The effective stress wavefront in the saturated sand is also much better formed, albeit still highly attenuated at depths greater than 4 m (Z/D = 0.2). High excess pore pressure is developed over a significant depth in the saturated zone; this being consistent with the longer duration of the loading. The propagation of the pore pressure peak suggests a wave velocity of only about 8 m/s, which is consistent with that of the slow dilatational wave. Hence, even for the longer loading duration, much of the impulsive loading is still translated into pore pressure; only a small proportion appears as effective stress increase. Fig. 22 shows the corresponding profiles for a higher permeability of 102 m/s and loading duration of 1 s. As can be seen, the effective stress excursion in the saturated sand is now significantly higher and the excess pore pressure amplitude is correspondingly lower. The above discussion has practical implications for dynamic compaction of loose sand fill with a shallow groundwater table. In dynamic compaction of dry sand, an important compaction mechanism is transient effective stress excursion during compaction [29]. In coastal reclaimed land, loose hydraulic sand fill, with shallow water table, is often encountered [30,31]. Hence, dynamic compaction in such areas is akin to the application of impulsive loading on dry sand underlain by saturated sand. As the discussion shows, if the loading is too rapid or if the permeability of the sand is too low, the dry solid wave reaching the groundwater table is translated into a fast wave which propagates through the pore water; the slow wave is insignificant. Since the pore water has a much higher modulus than the soil skeleton, the latter will experience almost no excursion in effective stress. In other words, the fast wave has little or no compaction effect. As the loading rate is reduced, the slow wave becomes much better formed. The effect of the slow wave is to cause compression of the soil skeleton. For this to occur, pore water must be expelled from the voids. If the permeability of the sand is low, the viscous drag between soil skeleton and water will cause both to compress concurrently. Since the pore water is much less compressible than the soil skeleton, much of the loading will still be translated into excess pore pressure, Fig. 21. To increase effective stress excursion in the saturated sand,

328

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

0

-0.2

-0.2

-0.4

-0.4

Z/D

Z/D

0

-0.6

-0.6

t=10ms t=16ms t=32ms t=80ms t=200ms

-0.8

-0.8

(a)

-1 -1

-0.8

-0.6

t=10ms t=16ms t=32ms t=80ms t=200ms

-0.4

-0.2

0

(a)

-1 -1

0.2

0

0.2

σ 'v /P0

σ 'v /P0 0

0

-0.2

-0.2

-0.4

Z/D

Z/D

-0.4

-0.6

-0.6

t =10ms t =16ms t =32ms t =80ms t =200ms

-0.8

-0.8

(b)

-1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

u/P0 Fig. 20. Normalized (a) vertical effective stress and (b) pore pressure profiles beneath the centre of the rigid cushion at different instances, for loading duration of 0.01 s and k = 104 m/s, after removal of geostatic effective stress and hydrostatic pore pressure.

the loading duration has to be increased and the soil must be sufficiently permeable. The former condition may be achieved, for instance, by underlining the tamper with a soft cushion. 6. Conclusions The foregoing discussion demonstrates a readily implementable method for fully coupled u–U type analysis within the framework of a commercial code. The examples presented show good agreement with published analytical and numerical solutions. The DPC algorithm also shows a high degree of computational stability. The three-dimensional elasto-plastic example also demonstrates the potential application of DPC u–U analysis on more realistic transient loading problems, while allowing users to access the computational functionalities and stability, as well as pre- and post-processing features that come with such commercial softwares. Although the examples considered in this paper focus

t=10ms t=16ms t=32ms t=80ms t=200ms

(b)

-1 0

0.1

u/P0 Fig. 21. Normalized (a) vertical effective stress and (b) pore pressure profiles beneath the centre of the rigid cushion at different instances, for loading duration of 1 s and k = 104 m/s, after removal of geostatic effective stress and hydrostatic pore pressure.

mainly on the dilatational wave response of the fluid and solid phases, the proposed formulation is, in principle, applicable for problems involving S- and/or R-wave propagation in a 2D or 3D domain with appropriate boundary conditions. The practical example presented is that of an impulsive load over a square patch on the surface of a sandy ground with high water table. The results of the analysis highlight the importance role played by the slow dilatational wave in the compaction process. For the slow wave to be effective, the results show that the soil has to be sufficiently permeable and the loading duration sufficiently long. It is plausible that there is a relationship with permeability and loading duration. However, this is outside the scope of this paper which is to present the DPC methodology and to demonstrate its usage. The last example also demonstrates that the DPC algorithm is also applicable to problems with drainage, as long as the relative movement between soil skeleton and water is not too large to

F. Ye et al. / Computers and Geotechnics 55 (2014) 316–329

329

References

0

-0.2

Z/D

-0.4

-0.6

t=10ms t=16ms t=32ms t=80ms t=200ms

-0.8

(a)

-1 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

σ 'v /P0 0

-0.2

Z/D

-0.4

-0.6

t=10ms t=16ms t=32ms t=80ms t=200ms

-0.8

(b)

-1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

u/P0 Fig. 22. Normalized (a) vertical effective stress and (b) pore pressure profiles beneath the centre of the rigid cushion at different instances, for loading duration of 1 s and k = 102 m/s, after removal of geostatic effective stress and hydrostatic pore pressure.

violate the small deformation assumption. Although the DPC approach was developed on the ABAQUS platform, it is applicable, in principle, to any existing codes with the following features: (1) A Cartesian damper element which allows viscous drag forces to be related to the relative nodal velocity. (2) An interface for user-defined material which allows constitutive coupling between the two phases to be programmed.

Acknowledgement The research scholarship for the first author provided by the National University of Singapore is hereby gratefully acknowledged.

[1] Zienkiewicz OC, Chang CT, Bettess P. Drained, undrained, consolidating and dynamic behaviour assumptions in soils. Geotechnique 1980;30(4):385–95. [2] Biot MA. Theory of propagation of elastic waves in a fluid-saturated porous solid, I. Low-frequency range. J Acoust Soc Am 1956;28:168–78. [3] Biot MA. Theory of propagation of elastic waves in a fluid-saturated porous solid, II. Higher frequency range. J Acoust Soc Am 1956;28:179–91. [4] Truesdell C. Sulle basi della termomeccanica. Rend Lincei 1957;22. p. 33–8, 158–66. [5] Truesdell C, Toupin R. The classical field theories. In: Flugge S, editor. Handbuch der Physik, vol. III/1. Berlin: Springer-Verlag; 1960. [6] Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T. Equations governing the dynamic, soil-pore fluid, interaction. In: Computational geomechanics with special reference to earthquake engineering. New York: Wiley & Sons; 1998. p. 17–52. [7] Wang ZQ, Hao H, Lu Y. A three-phase soil model for simulating stress wave propagation due to blast loading. Int J Numer Anal Meth Geomech 2004;28:33–56. [8] Grujicic M, Pandurangan B, Huang Y, Cheeseman BA, Roy WN, Skaggs RR. Impulse loading resulting from shallow buried explosives in water-saturated sand. Proc Inst Mech Eng L: J Mater Des Appl 2006;220:1–15. [9] An J, Tuan CY, Cheeseman BA, Gazonas A. Simulation of soil behaviour under blast loading. Int J Geomech 2011;11:323–34. [10] Lu CW. The experimental study on the reclaimed soils subjected to various impact energy. Master thesis, National Cheng Kung University, Taiwan; 2008. [11] Thevanayagam S, Nashed R, Martin GR. Dynamic compaction of saturated sands and silty sands: theory. Ground Improve 2009;162(G12):57–68. [12] Ghassemi A, Pak A, Shahir H. Numerical study of the coupled hydromechanical effects in dynamic compaction of saturated granular soils. Comput Geotech 2010;37:10–24. [13] Oumeraci H, Partenscky HW. Wave-induced pore pressure in rubble mound breakwaters. In: Edge BL, editor. Proceedings of 22nd conference on coastal engineering. The Netherlands: ASCE; 1990. p. 11334–47. [14] Mase H, Sakai T, Sakamoto M. Wave-induced porewater pressures and effective stresses around breakwater. Ocean Eng 1994;21(4):361–79. [15] Troch P, De Rouch J, Burcharth HF. Experimental study and numerical modelling of wave induced pore pressure attenuation inside a rubble mound breakwater. In: Smith JM, editor. Proceedings of the 28th international conference on coastal engineering (ICCE). UK: World Scientific; 2003. p. 1607–9. [16] Oumeraci H, Kortenhaus A, Allsop W, De Groot M, Crouch R, Vrijling H, et al. Probabilistic design tools for vertical breakwaters. Rotterdam: A. A. Balkema Publishers; 2001. [17] Garg SK, Nayfeh AH, Good AJ. Compressional waves in fluid-saturated elastic porous media. J Appl Phys 1974;45(5):1968–74. [18] Simon BR, Zienkiewicz OC, Paul DK. An analytical solution for the transient responses of saturated porous elastic solids. Int J Numer Anal Meth Geomech 1984;8:381–98. [19] Gajo A, Mongiovi L. An analytical solution for the transient response of saturated linear elastic porous media. Int J Numer Anal Meth Geomech 1995;19:399–413. [20] Schanz M, Cheng AHD. Transient wave propagation in a one-dimensional poroelastic column. Acta Mech 2000;145:1–18. [21] Shan ZD, Ling DS, Ding HJ. Exact solutions for one-dimensional transient response of fluid-saturated porous media. Int J Numer Anal Meth Geomech 2011;35:461–79. [22] Ghaboussi J, Dickmen SU. Liquefaction analysis of horizontally layered sands. J Geotech Eng Div 1978;104(GT3):341–56. [23] Zienkiewicz OC, Shiomi T. Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution. Int J Numer Anal Meth Geomech 1984;8:71–96. [24] Hiremath MS, Sandhu RS, Morland LW, Wolfe WE. Analysis of onedimensional wave propagation in a fluid-saturated finite soil column. Int J Numer Anal Meth Geomech 1988;12:121–39. [25] Qiu T, Fox PJ. Numerical analysis of 1-D compression wave propagation in saturated poroelastic media. Int J Numer Anal Meth Geomech 2007;32:161–87. [26] Kim KJ, Blouin SE, Timian DA. Experimental and theoretical response of multiphase porous media to dynamic loads. Annual Report 2 F49620–85-C0102, US Air Force Office of Scientific Research; 1987. [27] Kim KJ, Blouin SE, Chitty DE, Merkle DH. Experimental and theoretical response of multiphase porous media to dynamic loads. Final Report F49620– 85-C-0102, US Air Force Office of Scientific Research; 1988. [28] Prevost JH. DYNAFLOW: a nonlinear transient finite element analysis programme. US: Princeton University; 1981. [29] Gu Q, Lee FH. Ground response to dynamic compaction of dry sand. Geotechnique 2002;52(7):481–93. [30] Na YM, Choa V, Teh CI, Chang MF. Geotechnical parameters of reclaimed sandfill from the cone penetration test. Can Geotech J 2011;42:91–109. [31] Bo MW, Chang MF, Arulrajah A, Choa V. Ground investigations for Changi east reclamation projects. Geotech Geol Eng 2012;30:45–62.