**2. Description of the numerical method**

Sugano 2000; Ichii et al. 2000; Inoue et al. 2003; Nozu et al. 2004; Mostafavi Moghadam et al. 2009 and 2011). A new design methodology, named performance-based design, has born from lessons learned caused by earthquakes in 1990's to overcome the limitations of conventional seismic design (PIANC 2001). In this framework, lateral spreading of the saturated backfill and foundation soils along with the effect of quay wall as the supporting structure (saturated

Predicting the response of a structure retained a liquefiable soil during an earthquake is highly dependent on adequately accounting for the effects of pore water pressure development, stress-strain softening and strength reduction in the soil on the system behavior. Thus, it is required to perform dynamic analyses that account for the saturated soil-structure interaction effects using numerical modeling techniques. Several researchers have reported the use of numerical analysis for predicting the behavior of liquefiable soil measured in laboratory tests or field case histories. Iai et al. (1998) reported that FLIP code can successfully predict the seismic behavior of port structures. Yang and Parra developed CYCLIC code and reported successful predictions of the seismic behavior of gravity quay wall placed on liquefiable sand and calibrated the numerical results with centrifuge tests (Parra 1996; Yang 2000). Both the codes DYSAC2 and DYNAFLOW are reported by Arulanandan (1996) as having adequately predicted the response of a submerged embankment subjected to dynamic loading in a centrifuge test. According to Madabhushi and Zeng (1998), the code SWANDYNE successfully predicted the seismic response of a gravity quay wall rested on liquefiable sand modeled in the centrifuge. Finn reported the successful validation of TARA-3 using centrifuge tests results (Finn et al. 1991). The successful use of FLAC package for prediction of the behavior of caisson retaining walls in liquefiable soils was reported by Dickenson and Yang (1998). They used a nonlinear effective stress analysis method based on the Mohr-Coulomb constitutive model and a pore water pressure increment scheme based on the work of Seed and his co-workers (e.g. Martin et al. 1975; Seed and DeAlba 1986). Likewise, McCullough and Dickenson, who used the same analysis method and soil model in FLAC, reported fairly good agreement between predicted and measured permanent horizontal displacements at top of five anchored sheet pile bulkhead walls in liquefiable soils subjected to different earthquakes in Japan between 1987 and 1993 (McCullough and Dickenson 1998). It should be noted that the assumption of the Yang and McCullough was that the foundation soil is non-liquefiable. A hyperbolic type stress-strain formulation developed by Pyke (1979) along with the pore water pressure build up model proposed by Byrne (1991) was implemented in FLAC code by Cooke (2001). They concluded that Pyke-Byrne model over-predicts the horizontal displacement of gravity quay

soil-structure interaction) are taken into account as a more logical design.

258 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

walls modelled in centrifuge within a factor of approximately two.

Recently, Ebrahimian et al. (2009) have carried out a series of two dimensional fully coupled effective stress dynamic analyses in order to study the deformation of quay walls and the liquefaction potential of backfill soils. Additionally, several 1g shaking table tests have been executed to verify the obtained numerical results. They showed based on the lessons learned obtained from numerical results, a safer design of gravity-type quay walls can be developed. Correspondingly, Mostafavi Moghadam et al. (2009) conducted finite difference effectivestress analyses to investigate the seismic performance of caisson quay walls. Their obtained In this research, a two-dimensional (2D) reference model is developed to simulate the seismic performance of gravity-type quay walls in a rational way. Nonlinear time history dynamic analysis is conducted using computer program FLAC 2D (Itasca 2004). FLAC 2D is an explicit finite difference program for modeling soil-structure interaction analysis under static and seismic loading conditions. Here, numerical approach is based on a continuum finite difference discretization applying Lagrangian approach (Itasca 2004). Every derivative in the set of governing equations is directly replaced by algebraic expression written in terms of field variables (e.g., stress or displacement) at discrete point in space. Regarding dynamic analysis, explicit finite difference scheme is applied to solve the full equation of motion using the lumped grid point masses derived from the real density of surrounding zone. The calculation sequence first invokes the equations of motion for deriving new velocities and displacements from stresses and forces; then, strain rates are derived from velocities, and new stresses from strain rates. Every cycle around the loop corresponds to one time step. Each box updates all grid variables from known values which are fixed over the time step being executed (Figure 1).

**Figure 1.** Basic explicit calculation cycle (Itasca 2004)

The equation of motion, in the simplest form, relates the acceleration (*du*˙ / *dt*) of a mass (*m*) to the applied force (*F*) which may vary with time. Newton's law of motion for the mass-spring system is:

$$m\frac{d\dot{u}}{dt} = F\tag{1}$$

s

retained by quay wall rested on a very dense foundation.

**2.2. Soil constitutive model**

may or may not exist).

**2.1. Modeling procedure**

 sk

where, *M* = specific rule of behaviour; *κ*= history parameters (based on the specific rules which

Firstly, a static analysis considering the effect of gravity with the Mohr-Coulomb elastic perfectly plastic constitutive soil model is performed to establish the in-situ stresses before seismic loading. The water within the soil is modeled directly, and is allowed to flow during the static solutions. Static solutions are obtained by including damping terms that gradually remove kinetic energy from the system. Then, boundary pressures equivalent to fluid weight is applied along the bottom of the sea and the front side of the caisson. Once initial stress state is established in the model and the soil model is changed to a pore pressure generation constitutive model, the effective stress dynamic analysis is started. For dynamic analysis, the acceleration record is applied to the nodes along the bottom of numerical grid. During the dynamic solutions excess pore water pressures are allowed to generate and also the dissipation of these pore pressures is modelled. To this end, the Finn and Byrne model (Finn et al. 1977; Byrne 1991) is modified and used to carry out coupled dynamic groundwater flow calculations. This effective stress analysis which take into account the effects of seismically induced pore water pressures is used to investigate the degree and extend of liquefaction that occur in backfill. It should be noted that, all analyses are performed under fully drained condition. One key advantage of the coupled numerical approach is the ability to account for the interde‐ pendent effects of various mechanisms and phenomena on each other as the numerical computations proceed. For instance, when an effective stress formulation is used, the inclusion of pore water pressure generation in the simulation can impact the strength and stress-strain behavior of the soil during shaking in the same manner as in the field. For the sake of simplicity in defining the performance of a quay wall during an earthquake, the author has focused on the horizontal movement because most of the gravity-type quay walls exhibited a failure pattern predominantly in the form of excessive horizontal movement rather than other damages. The numerical grid for a typical quay wall and foundation soil used in the current study is illustrated in Figure 2. Numerical model represents a backfill of constant depth

In static analysis, the Mohr-Coulomb constitutive model is used to model the behaviour of sandy soil. The linear behaviour is defined by elastic shear and drained bulk modulus. The shear modulus of sandy soil is calculated with the formula given by Seed and Idriss (1970):

max 2 max 1000 *G k <sup>m</sup>* =

( ) 0.5

¢ (6)

s

*ij ij ij* = *M e* ( , , ) & (5)

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

http://dx.doi.org/10.5772/55027

261

In a continuous solid body, Eq. (1) is generalized as follows:

$$
\rho \frac{\partial \dot{u}\_i}{\partial t} = \frac{\partial \sigma\_{ij}}{\partial \mathbf{x}\_j} + \rho \mathbf{g}\_i \tag{2}
$$

where, *ρ* = mass density; *t* = time; *xj* = components of coordinate vector; *gi* = components of gravitational acceleration (body forces); *σij* = components of stress tensor; *i* = components in a Cartesian coordinate frame.

For problem analysis, the strain rate tensor and rotation rate tensor, having the velocity gradients, are calculated by the following equations:

$$e\_{ij} = \frac{1}{2} \left[ \frac{\partial \dot{u}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \dot{u}\_j}{\partial \mathbf{x}\_i} \right] \tag{3}$$

$$\alpha o\_{ij} = \frac{1}{2} \left[ \frac{\partial \dot{u}\_i}{\partial \mathbf{x}\_j} - \frac{\partial \dot{u}\_j}{\partial \mathbf{x}\_i} \right] \tag{4}$$

where, *eij*= components of strain rate; *ωij*= components of rotation rate; *u*˙*<sup>i</sup>* = components of velocity.

The specific mechanical relationship is used in order to obtain the stress tensor as below:

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls http://dx.doi.org/10.5772/55027 261

$$
\sigma\_{ij} = M \left( \sigma\_{i\rangle}, \dot{e}\_{i\rangle}, \kappa \right) \tag{5}
$$

where, *M* = specific rule of behaviour; *κ*= history parameters (based on the specific rules which may or may not exist).

#### **2.1. Modeling procedure**

**Figure 1.** Basic explicit calculation cycle (Itasca 2004)

260 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

where, *ρ* = mass density; *t* = time; *xj*

Cartesian coordinate frame.

velocity.

system is:

The equation of motion, in the simplest form, relates the acceleration (*du*˙ / *dt*) of a mass (*m*) to the applied force (*F*) which may vary with time. Newton's law of motion for the mass-spring

*dt* <sup>=</sup> & (1)

& (2)

& & (3)

& & (4)

= components of

= components of

*du m F*

*ij i*

*u*

¶ ¶ = + ¶ ¶

> 1 2

> > 1 2

*ij*

*ij*

where, *eij*= components of strain rate; *ωij*= components of rotation rate; *u*˙*<sup>i</sup>*

w

*e*

r

*j*

gravitational acceleration (body forces); *σij* = components of stress tensor; *i* = components in a

For problem analysis, the strain rate tensor and rotation rate tensor, having the velocity

*j i*

*j i u u*

*j i*

é ù ¶ ¶ = - ê ú ê ú ¶ ¶ ë û

The specific mechanical relationship is used in order to obtain the stress tensor as below:

*j i u u x x*

*x x* é ù ¶ ¶ = + ê ú ê ú ¶ ¶ ë û

*<sup>g</sup> t x* s

*i*

= components of coordinate vector; *gi*

 r

In a continuous solid body, Eq. (1) is generalized as follows:

gradients, are calculated by the following equations:

Firstly, a static analysis considering the effect of gravity with the Mohr-Coulomb elastic perfectly plastic constitutive soil model is performed to establish the in-situ stresses before seismic loading. The water within the soil is modeled directly, and is allowed to flow during the static solutions. Static solutions are obtained by including damping terms that gradually remove kinetic energy from the system. Then, boundary pressures equivalent to fluid weight is applied along the bottom of the sea and the front side of the caisson. Once initial stress state is established in the model and the soil model is changed to a pore pressure generation constitutive model, the effective stress dynamic analysis is started. For dynamic analysis, the acceleration record is applied to the nodes along the bottom of numerical grid. During the dynamic solutions excess pore water pressures are allowed to generate and also the dissipation of these pore pressures is modelled. To this end, the Finn and Byrne model (Finn et al. 1977; Byrne 1991) is modified and used to carry out coupled dynamic groundwater flow calculations. This effective stress analysis which take into account the effects of seismically induced pore water pressures is used to investigate the degree and extend of liquefaction that occur in backfill. It should be noted that, all analyses are performed under fully drained condition. One key advantage of the coupled numerical approach is the ability to account for the interde‐ pendent effects of various mechanisms and phenomena on each other as the numerical computations proceed. For instance, when an effective stress formulation is used, the inclusion of pore water pressure generation in the simulation can impact the strength and stress-strain behavior of the soil during shaking in the same manner as in the field. For the sake of simplicity in defining the performance of a quay wall during an earthquake, the author has focused on the horizontal movement because most of the gravity-type quay walls exhibited a failure pattern predominantly in the form of excessive horizontal movement rather than other damages. The numerical grid for a typical quay wall and foundation soil used in the current study is illustrated in Figure 2. Numerical model represents a backfill of constant depth retained by quay wall rested on a very dense foundation.

#### **2.2. Soil constitutive model**

In static analysis, the Mohr-Coulomb constitutive model is used to model the behaviour of sandy soil. The linear behaviour is defined by elastic shear and drained bulk modulus. The shear modulus of sandy soil is calculated with the formula given by Seed and Idriss (1970):

$$G\_{\text{max}} = 1000 \, k\_{2\,\text{max}} \left( \sigma\_m' \right)^{0.5} \tag{6}$$

**Figure 2.** Numerical grid constructed in FLAC

where, Gmax is the maximum (small strain) shear modulus in pounds per square foot, psf (it is later converted to kPa to be consistent with metric units being used), k2max is the shear modulus number (Seed and Idriss 1970), and *σ* ′ *<sup>m</sup>* is the mean effective confining stress in psf. The Poisson's ratio is taken as 0.35.

For pore water generation during dynamic analysis, the updated model proposed by Byrne (1991) is incorporated to account the development of pore water pressure build up as an effect of volumetric strain induced by the cyclic shear strain using the following formulation:

$$
\Delta \varepsilon\_v = \mathbb{C}\_1 \exp \left( \mathbb{C}\_2 \, \varepsilon\_v / \gamma \right) \tag{7}
$$

Modulus degradation curves imply a nonlinear stress-strain curve. An incremental constitu‐ tive relation can be derived from the degradation curve, described by τ/γ = Ms, where τ is the normalized shear stress, γ is the shear strain and Ms is the normalized secant modulus. The

> g

Equivalent linear analysis is the common method used for evaluating the seismic behaviour of earth structures. In this approach, first, the responses are linearly analyzed using the initial values of damping ratio and shear modulus. Then, the new values of damping ratio and shear modulus are estimated, using maximum value of shear strain and laboratory curves. These values are used for redoing the analysis. This procedure is repeated several times until the material properties show no variation. Therefore, no non-linear effect is directly captured by this method as it assumes linearity during the solution process. Strain-dependent modulus and damping functions are considered roughly in order to approximate some effects of non-

In the non-linear analysis, employed in this study, the non-linear stress-strain relationship is directly followed by each zone. Damping ratio and shear modulus of the materials are calculated automatically at different strain levels. The real behaviour of soil, under cyclic loading, is non-linear and hysteretic. Such behaviour can be simulated by Masing model (Masing 1926), which can model the dynamic behaviour of soil. In this model, the shear

> ( ) ( ) max max max <sup>1</sup> *bb G*

g

t

where, *Fbb(γ)* = backbone or skeleton function; *γ*= shear strain amplitude; *Gmax* = initial shear

Stress-strain curve follows the backbone curve in the first loading, as shown in Figure 3(a); however, for explaining the unload-reload process, the above equation should be modified. If load reversal occurs at the point (*τr*, *γr*), stress-strain curve follows the path given by the below

 g

<sup>=</sup> <sup>+</sup> (11)


*G*

2 2 *r r bb F*

 gg

tt

 g

= =+ (10)

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

, where G is the

http://dx.doi.org/10.5772/55027

263

, is descried as

t

g

The incremental shear modulus in a nonlinear simulation is then given by GMt

. *<sup>s</sup> t s <sup>d</sup> dM M M d d*

normalized tangent modulus, Mt

small-strain shear modulus of the material.

linearity (damping and material softening).

behaviour of soil may be explained by a backbone curve as:

*F*

modulus; *τmax* = maximum shear stress amplitude.

formula:

g

**2.3. Full non-linear dynamic analysis**

where, Δεv is the volumetric strain increment that occurs over the current cycle, εv is the accumulated volumetric strain for previous cycle, γ is the shear strain amplitude for the current cycle, and C1 and C2 are constants dependent on the volumetric strain behavior of sand. According to Byrne (1991), the constant C1 in equation 7 controls the amount of volumetric strain increment and C2 controls the shape of volumetric strain curve. These constants are estimated using:

$$C\_1 = 7600 \left(Dr\right)^{-2.5} \tag{8}$$

$$\mathbf{C\_2} = 0.4 \mathbf{\dot{/C\_1}}\tag{9}$$

where, Dr is the relative density of soil in percent. To provide constitutive model that can better fit the curves of shear modulus reduction and damping ratio derived from the experimental tests data, two different modifications to soil model are implemented as a part of this research to assess the potential for predicting liquefaction phenomenon and associated deformations. To represent the nonlinear stress-strain behavior of soil more accurately that follows the actual stress-strain path during cyclic loading, the masing behavior is implemented into FLAC which works with Byrne model by a FISH subroutine as a first modification.

Since, there is a need to accept directly the same degradation curves derived from the test data in fully nonlinear method to model the correct physics, the second modification is related to incorporate such cyclic data into a hysteretic damping model for FLAC.

Modulus degradation curves imply a nonlinear stress-strain curve. An incremental constitu‐ tive relation can be derived from the degradation curve, described by τ/γ = Ms, where τ is the normalized shear stress, γ is the shear strain and Ms is the normalized secant modulus. The normalized tangent modulus, Mt , is descried as

$$\mathbf{M}\_{\rm t} = \det\_{\mathbf{d}\chi} \Big/\_{\mathbf{d}\chi} = \mathbf{M}\_{\rm s} + \boldsymbol{\chi}. \quad \Big/\_{\mathbf{d}\chi} \Big/\_{\mathbf{d}\chi} \tag{10}$$

The incremental shear modulus in a nonlinear simulation is then given by GMt , where G is the small-strain shear modulus of the material.

#### **2.3. Full non-linear dynamic analysis**

**Figure 2.** Numerical grid constructed in FLAC

Poisson's ratio is taken as 0.35.

estimated using:

number (Seed and Idriss 1970), and *σ* ′

262 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

where, Gmax is the maximum (small strain) shear modulus in pounds per square foot, psf (it is later converted to kPa to be consistent with metric units being used), k2max is the shear modulus

For pore water generation during dynamic analysis, the updated model proposed by Byrne (1991) is incorporated to account the development of pore water pressure build up as an effect of volumetric strain induced by the cyclic shear strain using the following formulation:

> eg

where, Δεv is the volumetric strain increment that occurs over the current cycle, εv is the accumulated volumetric strain for previous cycle, γ is the shear strain amplitude for the current cycle, and C1 and C2 are constants dependent on the volumetric strain behavior of sand. According to Byrne (1991), the constant C1 in equation 7 controls the amount of volumetric strain increment and C2 controls the shape of volumetric strain curve. These constants are

( ) 2.5

where, Dr is the relative density of soil in percent. To provide constitutive model that can better fit the curves of shear modulus reduction and damping ratio derived from the experimental tests data, two different modifications to soil model are implemented as a part of this research to assess the potential for predicting liquefaction phenomenon and associated deformations. To represent the nonlinear stress-strain behavior of soil more accurately that follows the actual stress-strain path during cyclic loading, the masing behavior is implemented into FLAC which

Since, there is a need to accept directly the same degradation curves derived from the test data in fully nonlinear method to model the correct physics, the second modification is related to

( ) 1 2 exp D =

e

works with Byrne model by a FISH subroutine as a first modification.

incorporate such cyclic data into a hysteretic damping model for FLAC.

*<sup>m</sup>* is the mean effective confining stress in psf. The

*v v C C* (7)

<sup>1</sup> *C Dr* <sup>7600</sup> - <sup>=</sup> (8)

2 1 *C C* = 0.4 (9)

Equivalent linear analysis is the common method used for evaluating the seismic behaviour of earth structures. In this approach, first, the responses are linearly analyzed using the initial values of damping ratio and shear modulus. Then, the new values of damping ratio and shear modulus are estimated, using maximum value of shear strain and laboratory curves. These values are used for redoing the analysis. This procedure is repeated several times until the material properties show no variation. Therefore, no non-linear effect is directly captured by this method as it assumes linearity during the solution process. Strain-dependent modulus and damping functions are considered roughly in order to approximate some effects of nonlinearity (damping and material softening).

In the non-linear analysis, employed in this study, the non-linear stress-strain relationship is directly followed by each zone. Damping ratio and shear modulus of the materials are calculated automatically at different strain levels. The real behaviour of soil, under cyclic loading, is non-linear and hysteretic. Such behaviour can be simulated by Masing model (Masing 1926), which can model the dynamic behaviour of soil. In this model, the shear behaviour of soil may be explained by a backbone curve as:

$$F\_{kb}\left(\boldsymbol{\gamma}\right) = \frac{G\_{\text{max}}\boldsymbol{\gamma}}{1 + \left(G\_{\text{max}}/\tau\_{\text{max}}\right) \left|\boldsymbol{\gamma}\right|}\tag{11}$$

where, *Fbb(γ)* = backbone or skeleton function; *γ*= shear strain amplitude; *Gmax* = initial shear modulus; *τmax* = maximum shear stress amplitude.

Stress-strain curve follows the backbone curve in the first loading, as shown in Figure 3(a); however, for explaining the unload-reload process, the above equation should be modified. If load reversal occurs at the point (*τr*, *γr*), stress-strain curve follows the path given by the below formula:

$$\frac{\tau - \tau\_r}{2} = F\_{bb} \left[ \frac{\gamma - \gamma\_r}{2} \right] \tag{12}$$

In other words, the shapes of unload-reload curves are similar to that of backbone curve (with the origin shifted to the loading reversal point) except they are enlarged by a factor of 2, as shown in Figure 3(b). The Equations (9) and (10) describe the Masing behaviour (Masing 1926).

Masing rules seem not to be enough for precise explanation of soil response under general cyclic loading. Finn et al. (1977) developed modified rules to describe the irregular loading. They suggested that unloading and reloading curves follow the concerning two rules. If the new unloading or reloading curve exceeds the last maximum strain and cut the backbone curve, it will follow the backbone curve up to meeting the next returning point, as shown in Figure 3(c). If a new unloading or reloading curve passes through the previous one, it will follow the former stress-strain curve, as shown in Figure 3(d). According to this model, the tangent shear modulus can be defined at the points on the backbone and new reloadingunloading curves by the Formulas (11) and (12), respectively, as:

$$G\_t = G\_{\text{max}} \sqrt{\left[1 + \frac{G\_{\text{max}} \left| \mathcal{V} \right|}{\tau\_{\text{max}}} \right]^2} \tag{13}$$

where, *α* =3.3 , *β* =1.3 and *η* =10 are constant coefficients, *f* is the base acceleration frequency

(g<sup>r</sup> ,tr) *Gmax*

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

http://dx.doi.org/10.5772/55027

265

*Gmax*

*General Loading*

(a) Rule 1 (b) Rule 2

*Gmax*

*Local Loading*

*Gmax*

*(kN/m3) (MPa) (MPa) (degree)* Caisson quay wall 2400 83.33 111.11 - - - Backfill soil 1520 35 75.83 25 2.432 0.164 Foundation soil 1600 130 281.66 38 0.114 3.508 Rubble foundation 1800 155 206.66 40 - -

The structural element consists of a gravity quay wall supported on very dense foundation soil. The structural element is modeled within FLAC as elastic element. The soil-structure interaction between soil and wall including normal and shear springs. The FLAC manual (Itasca 2004) recommends as a rule of thumb that ks and kn be set to ten times the equivalent

(d) Rule 4

**γ***<sup>d</sup>* **G K Friction C1 C2**

(g<sup>r</sup> ,tr)

*Gmax*

(g<sup>r</sup> ,tr)

**Figure 3.** General patterns of loading, unloading and reloading paths in Masing model

stiffness of the stiffest neighboring zone (see equation (16) below).

*Gmax*

(c) Rule 3

*Gmax*

and *t* is the time.

**Table 1.** Soil properties

**2.5. Structural properties**

$$\mathbf{G}\_t = \mathbf{G}\_{\text{max}} \left\langle \left[ 1 + \frac{\mathbf{G}\_{\text{max}}}{2\tau\_{\text{max}}} \left| \boldsymbol{\nu} - \boldsymbol{\nu}\_r \right| \right]^2 \right. \tag{14}$$

Based on the results, obtained in this research, the shear stress decreases as the number of load cycles increases; it means that shear stress-strain curves are more inclined. In this study, Masing rules are implemented into FLAC via a series of FISH functions in order to simulate the non-linear stress-strain relationships.

#### **2.4. Material properties and input seismic loading**

For the gravity quay wall, the concrete caisson is modeled by linear elastic elements. The granular backfill is modeled as a purely frictional, elastic–plastic soil with a Mohr–Coulomb failure criterion. Strength properties for sandy soil at relative densities of 25 and 85 percent are obtained based on correlations in literature. Maximum shear moduli (i.e. modulus at a strain level of approximately 0.0001%) for Sandy soil are determined from formula developed by Seed and Idriss (1970), given in Equation (6). Hydraulic conductivities are derived from constant head permeability tests. Values for the constants C1 and C2 used in the volumetric strain equation (Equation (7)) are obtained from Equations 8 and 9. The material properties of soil layers are listed in Table 1.

After static equilibrium is achieved (end of static construction stage), the full width of the foundation is subjected to the variable-amplitude harmonic ground motion record illustrated in Figure 4. The mathematical expression for input acceleration is given by:

$$
\ddot{\mathbf{U}}(t) = \sqrt{\beta \mathbf{e}^{-\alpha t} \, t^{\eta}} \sin \left( 2 \pi f t \right) \tag{15}
$$

where, *α* =3.3 , *β* =1.3 and *η* =10 are constant coefficients, *f* is the base acceleration frequency and *t* is the time.

**Figure 3.** General patterns of loading, unloading and reloading paths in Masing model


**Table 1.** Soil properties

In other words, the shapes of unload-reload curves are similar to that of backbone curve (with the origin shifted to the loading reversal point) except they are enlarged by a factor of 2, as shown in Figure 3(b). The Equations (9) and (10) describe the Masing behaviour (Masing 1926).

Masing rules seem not to be enough for precise explanation of soil response under general cyclic loading. Finn et al. (1977) developed modified rules to describe the irregular loading. They suggested that unloading and reloading curves follow the concerning two rules. If the new unloading or reloading curve exceeds the last maximum strain and cut the backbone curve, it will follow the backbone curve up to meeting the next returning point, as shown in Figure 3(c). If a new unloading or reloading curve passes through the previous one, it will follow the former stress-strain curve, as shown in Figure 3(d). According to this model, the tangent shear modulus can be defined at the points on the backbone and new reloading-

2

2

(13)

(14)

max

*G*

= + ê ú

t

max

é ù

max

Based on the results, obtained in this research, the shear stress decreases as the number of load cycles increases; it means that shear stress-strain curves are more inclined. In this study, Masing rules are implemented into FLAC via a series of FISH functions in order to simulate

For the gravity quay wall, the concrete caisson is modeled by linear elastic elements. The granular backfill is modeled as a purely frictional, elastic–plastic soil with a Mohr–Coulomb failure criterion. Strength properties for sandy soil at relative densities of 25 and 85 percent are obtained based on correlations in literature. Maximum shear moduli (i.e. modulus at a strain level of approximately 0.0001%) for Sandy soil are determined from formula developed by Seed and Idriss (1970), given in Equation (6). Hydraulic conductivities are derived from constant head permeability tests. Values for the constants C1 and C2 used in the volumetric strain equation (Equation (7)) are obtained from Equations 8 and 9. The material properties of

After static equilibrium is achieved (end of static construction stage), the full width of the foundation is subjected to the variable-amplitude harmonic ground motion record illustrated

> p


( ) e sin 2( ) *<sup>t</sup> U t t ft* h

in Figure 4. The mathematical expression for input acceleration is given by:

b

ë û

t

= +- ê ú

ê ú ë û

é ù

max

g g

g

unloading curves by the Formulas (11) and (12), respectively, as:

264 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

*G G*

max

1 2 *t r*

max

the non-linear stress-strain relationships.

soil layers are listed in Table 1.

**2.4. Material properties and input seismic loading**

*<sup>G</sup> G G*

1 *<sup>t</sup>*

#### **2.5. Structural properties**

The structural element consists of a gravity quay wall supported on very dense foundation soil. The structural element is modeled within FLAC as elastic element. The soil-structure interaction between soil and wall including normal and shear springs. The FLAC manual (Itasca 2004) recommends as a rule of thumb that ks and kn be set to ten times the equivalent stiffness of the stiffest neighboring zone (see equation (16) below).

(16)

**2.7. Hydrodynamic effects**

parabolically with depth and is defined by

to move in the horizontal direction.

**2.8. Boundary conditions**

waves back into the model.

**2.9. Damping**

dynamic interaction of the wall and water.

Based on Westergaard's approach (1931), effects of hydrodynamic pressure can be approxi‐ mately considered as added masses acting with the quay wall. The added mass increases

Where m(y) is the variation of mass with depth y. Mw is the mass density of water and h is the overall depth of the water in front of the quay wall. The added mass can be reasonably modeled as a simple-supported beam element which possess the corresponding mass m(y) and is free

The water in front of the wall is modeled indirectly by including the resulting water pressures along the boundaries. This allows a simple modeling of the water, but does not account for the

Many geotechnical problems can be idealized by assuming that the regions remote from the area of interest extend to infinity. As the capability of computer is limited, the unbounded theoretical models have to be truncated to a manageable size by using artificial boundaries. In practice, the numerical model for the quay wall system should be extended to a sufficient depth below the ground level and to a sufficient width to consider local site effects and soil-structure interaction. During the static analysis, the bottom boundary is fixed in the both horizontal and vertical directions and the lateral boundaries are just fixed in the horizontal direction. In dynamic problems fixed boundary condition will cause the reflection of outward propagating

Thus, absorbing boundaries proposed by Lysmer and Kuhlemeyer (1969) are applied to the lateral boundaries in the model during dynamic analyses. It is based on the use of independent

Material damping in a soil is generally caused by its viscous properties, friction and the developmentofplasticity.Indeed,theroleofthedampinginthenumericalmodelsistoreproduce in magnitude and form the energy losses in the natural system when subject to a dynamic load. The dynamic damping in the model is provided by the Rayleigh damping option provided in FLAC. A damping percentage of 5 percent is used which is a typical value for geologic materi‐ als (Itasca 2004). The damping frequency is chosen by examining the undamped behavior of the numerical model. A damped frequency of 1 Hz is used for the present model. In each dynamic analysis,5percentRayleighdampingisincludedforthesoilelementsinadditiontothehysteretic

dashpots in the normal and shear directions at the lateral boundaries.

damping already incorporated in the nonlinear stress-strain model.

( ) <sup>7</sup> . *m y M hy* <sup>=</sup> <sup>8</sup> *<sup>W</sup>* (17)

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

http://dx.doi.org/10.5772/55027

267

**Figure 4.** Seismic excitation applied to the bottom of numerical model

where K and G are the bulk and shear moduli; and *Δz*min is the smallest width of an adjoining zone. The max [ ] notation indicates the maximum value over all zones adjacent to the interface.

In the case of a rough wall, modelling the interface between the soil and the wall is invariably an integral part of the analysis. In the case of soil–structure interaction, the interface is considered stiff compared to the surrounding soil, but it can slip and open in response to the loading. Joints with zero thickness are more suitable for simulating the frictional behaviour at the interface between the wall and the soil. The interface model is described by Coulomb law (Itasca 2004) to simulate the soil/wall contact. The interface element properties are summarized in Table 2.


**Table 2.** Interface element properties

#### **2.6. Fluid properties**

The modeled fluid is representative of seawater with a unit weight of 1027 kg/m3 and a bulk modulus of 2 GPa. An estimated pore pressure distribution is initialized in the model prior to the initial static stress state solution.

## **2.7. Hydrodynamic effects**

( )

ë û

*Z*

**0 5 10 15 TIME (s)**

where K and G are the bulk and shear moduli; and *Δz*min is the smallest width of an adjoining zone. The max [ ] notation indicates the maximum value over all zones adjacent to the interface.

In the case of a rough wall, modelling the interface between the soil and the wall is invariably an integral part of the analysis. In the case of soil–structure interaction, the interface is considered stiff compared to the surrounding soil, but it can slip and open in response to the loading. Joints with zero thickness are more suitable for simulating the frictional behaviour at the interface between the wall and the soil. The interface model is described by Coulomb law (Itasca 2004) to simulate the soil/wall contact. The interface element properties are summarized

> Quay wall with rubble 1e11 1e11 20 Quay wall with backfill soil 1e11 1e11 12

The modeled fluid is representative of seawater with a unit weight of 1027 kg/m3 and a bulk modulus of 2 GPa. An estimated pore pressure distribution is initialized in the model prior to

**kn ks Friction** *(Pa/m) (Pa/m) (degree)*

*K G*

<sup>3</sup> 10 max *n s*

é ù <sup>+</sup> ê ú ==´ ê ú <sup>D</sup> ê ú

*k k*

266 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

**-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5**

**Figure 4.** Seismic excitation applied to the bottom of numerical model

**Surface**

**Table 2.** Interface element properties

the initial static stress state solution.

**2.6. Fluid properties**

in Table 2.

**ACCELERATION (m/s**

**2**

**)**

4

min

(16)

Based on Westergaard's approach (1931), effects of hydrodynamic pressure can be approxi‐ mately considered as added masses acting with the quay wall. The added mass increases parabolically with depth and is defined by

$$
\Delta m(y) = \bigvee\_{8} M\_{\mathcal{W}} \sqrt{\hbar, y} \tag{17}
$$

Where m(y) is the variation of mass with depth y. Mw is the mass density of water and h is the overall depth of the water in front of the quay wall. The added mass can be reasonably modeled as a simple-supported beam element which possess the corresponding mass m(y) and is free to move in the horizontal direction.

The water in front of the wall is modeled indirectly by including the resulting water pressures along the boundaries. This allows a simple modeling of the water, but does not account for the dynamic interaction of the wall and water.

### **2.8. Boundary conditions**

Many geotechnical problems can be idealized by assuming that the regions remote from the area of interest extend to infinity. As the capability of computer is limited, the unbounded theoretical models have to be truncated to a manageable size by using artificial boundaries. In practice, the numerical model for the quay wall system should be extended to a sufficient depth below the ground level and to a sufficient width to consider local site effects and soil-structure interaction. During the static analysis, the bottom boundary is fixed in the both horizontal and vertical directions and the lateral boundaries are just fixed in the horizontal direction. In dynamic problems fixed boundary condition will cause the reflection of outward propagating waves back into the model.

Thus, absorbing boundaries proposed by Lysmer and Kuhlemeyer (1969) are applied to the lateral boundaries in the model during dynamic analyses. It is based on the use of independent dashpots in the normal and shear directions at the lateral boundaries.

### **2.9. Damping**

Material damping in a soil is generally caused by its viscous properties, friction and the developmentofplasticity.Indeed,theroleofthedampinginthenumericalmodelsistoreproduce in magnitude and form the energy losses in the natural system when subject to a dynamic load. The dynamic damping in the model is provided by the Rayleigh damping option provided in FLAC. A damping percentage of 5 percent is used which is a typical value for geologic materi‐ als (Itasca 2004). The damping frequency is chosen by examining the undamped behavior of the numerical model. A damped frequency of 1 Hz is used for the present model. In each dynamic analysis,5percentRayleighdampingisincludedforthesoilelementsinadditiontothehysteretic damping already incorporated in the nonlinear stress-strain model.

#### **2.10. Element size**

Toavoidthenumericaldistortionofthepropagatingwaveindynamicanalysisthespatialelement size, ΔL, must be smaller than approximately one-tenth to one-eighth of the wavelength associat‐ ed with the highest frequency component of the input wave (Kuhlemeyer and Lysmer 1973):

$$
\Delta \mathcal{L} = \bigvee\_{\Delta} \tag{18}
$$

The one-zone sample is modeled with FLAC that this consist of a sandy soil which is given a periodic motion at its base. Vertical loading is by gravity only. Equilibrium stresses and pore pressures are installed in the soil, and pore pressure and effective stress (mean total stress minus the pore pressure) are established within the soil. The modified Byrne model is applied

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

http://dx.doi.org/10.5772/55027

269

Figure 6 indicates the pore pressure build up in a single zone. It can be seen that the effective stress reaches zero after about 20 cycles of shaking (2 seconds, at 5 Hz). At this point, lique‐ faction can be said to occur. This test is strain-controlled in the shear direction. The stress/strain loops for the one-zone sample for several cycles are shown in Figure 7. It can be observed that shear modulus decreases with increasing shear strain. The hysteretic model seems to handle multiple nested loops in a reasonable manner. There is clearly energy dissipation and shear stiffness degradation during seismic loading. Due to the satisfactory modeling of the validation case, the numerical model is used to perform parametric studies on caisson quay wall, as

for the soil. The results based on Equation (7) are shown in Figures 6 and 7.

**Figure 6.** Pore pressure generation and effective stress time histories during seismic loading

**-1.5E-03 -1.0E-03 -5.0E-04 0.0E+00 5.0E-04 1.0E-03 Shear strain**

**-3.0E+05 -2.5E+05 -2.0E+05 -1.5E+05 -1.0E+05 -5.0E+04 0.0E+00 5.0E+04 1.0E+05 1.5E+05 2.0E+05 2.5E+05**

**Shear stress (Pa)**

**Figure 7.** Hyteresis loop in a one-zone sample element

described in the pervious sections.

In general, the cut-off frequency for geotechnical earthquake engineering problems should be no less than 10 Hz (ASCE 2000). Considering above criteria, element size is defined small enough to allow seismic wave propagation throughout the analysis. A finer mesh is used in sensitive areas such as below and near quay wall. A coarser mesh has been chosen for other areas in order to save computer analysis time.

#### **2.11. Time step**

To complete the numerical solution, it is necessary to integrate the governing equations in time in an incremental manner. The time step of the solution should be sufficiently small to accurately define the applied dynamic loads and to ensure stability and convergence of the solution. In the current FLAC model, the time step is approximately 10-6 second.

The convergence criterion for FLAC is the ratio defined to be the maximum unbalanced force magnitude for all the gridpoints in the model divided by the average applied force magnitude for all the gridpoints.

### **3. Verification**

To validate the implementation of the masing rule and hysteresis damping in FLAC program, the simulation of one-zone sample with modified Byrne model is done by using the unit cell as shown in Figure 5 incorporated with the implemented rules.

**Figure 5.** One-zone model in FLAC for simulating cyclic simple shear test

The one-zone sample is modeled with FLAC that this consist of a sandy soil which is given a periodic motion at its base. Vertical loading is by gravity only. Equilibrium stresses and pore pressures are installed in the soil, and pore pressure and effective stress (mean total stress minus the pore pressure) are established within the soil. The modified Byrne model is applied for the soil. The results based on Equation (7) are shown in Figures 6 and 7.

Figure 6 indicates the pore pressure build up in a single zone. It can be seen that the effective stress reaches zero after about 20 cycles of shaking (2 seconds, at 5 Hz). At this point, lique‐ faction can be said to occur. This test is strain-controlled in the shear direction. The stress/strain loops for the one-zone sample for several cycles are shown in Figure 7. It can be observed that shear modulus decreases with increasing shear strain. The hysteretic model seems to handle multiple nested loops in a reasonable manner. There is clearly energy dissipation and shear stiffness degradation during seismic loading. Due to the satisfactory modeling of the validation case, the numerical model is used to perform parametric studies on caisson quay wall, as described in the pervious sections.

**Figure 6.** Pore pressure generation and effective stress time histories during seismic loading

**Figure 7.** Hyteresis loop in a one-zone sample element

**2.10. Element size**

**2.11. Time step**

for all the gridpoints.

**3. Verification**

areas in order to save computer analysis time.

268 Engineering Seismology, Geotechnical and Structural Earthquake Engineering

Toavoidthenumericaldistortionofthepropagatingwaveindynamicanalysisthespatialelement size, ΔL, must be smaller than approximately one-tenth to one-eighth of the wavelength associat‐ ed with the highest frequency component of the input wave (Kuhlemeyer and Lysmer 1973):

In general, the cut-off frequency for geotechnical earthquake engineering problems should be no less than 10 Hz (ASCE 2000). Considering above criteria, element size is defined small enough to allow seismic wave propagation throughout the analysis. A finer mesh is used in sensitive areas such as below and near quay wall. A coarser mesh has been chosen for other

To complete the numerical solution, it is necessary to integrate the governing equations in time in an incremental manner. The time step of the solution should be sufficiently small to accurately define the applied dynamic loads and to ensure stability and convergence of the

The convergence criterion for FLAC is the ratio defined to be the maximum unbalanced force magnitude for all the gridpoints in the model divided by the average applied force magnitude

To validate the implementation of the masing rule and hysteresis damping in FLAC program, the simulation of one-zone sample with modified Byrne model is done by using the unit cell

solution. In the current FLAC model, the time step is approximately 10-6 second.

as shown in Figure 5 incorporated with the implemented rules.

**Figure 5.** One-zone model in FLAC for simulating cyclic simple shear test

(18)

<sup>9</sup> D = *<sup>L</sup>* l

A series of 1g shaking table tests have been executed in order to verify the obtained numerical results. It is attempt to create almost similar conditions between laboratory model test and numerical model. The liquefiable soil is modeled by loose sand and non-liquefiable soil is modeled by very dense sand. The seismic excitation is shown in Figure 8. The numerical results are presented and compared to those of corresponding shaking table test. Figure 9(a) shows the permanent deformation pattern of the numerical model after dynamic excitation. The nodal displacement vectors are presented in Figure 9(b). As may be expected, more ground surface settlement is observed in the backfill near the wall than at far field. A rigid body rotation of the wall (tilt) to the seaward direction is also clearly seen. The deformation pattern of model test at the end of seismic loading is presented in Figure 10. The trend of deformation behind quay wall and movement of the wall are in fairly good agreement with numerical results. Comparisons between the calculated and measured results (Calculated: numerical and Measured: experimental results) are made in Figure 11.

One might notice that the calculated results are rather close to measured values. This clearly demonstrates that the current numerical procedure captures very well the seismic behavior of

> **Calculated Measured**

**Calculated Measured**

(a) (b)

**Total pressure (Pa)**

**Horizontal Displacement (cm)**

> **0 5 10 15 20 Time (second)**

Numerical Modelling of the Seismic Behaviour of Gravity-Type Quay Walls

**0 5 10 15 20 TIME (s)**

**Calculated Measured**

http://dx.doi.org/10.5772/55027

271

**Calculated Measured**

(c) (d)

**Figure 11.** Recorded versus computed (a) acceleration, (b) horizontal displacement at top of the quay wall, and (c)

Results of nonlinear effective-stress dynamic analyses are presented in this section to investi‐ gate the effects of soil properties and input excitation characteristics on liquefaction potential, deformation of quay wall and failure mechanisms of soil-wall system during seismic loading. For gravity quay walls on firm foundations, typical failure modes during earthquakes are seaward displacements and tilting. Therefore, the horizontal displacement of quay wall head

Three additional sets of soil properties with different relative densities are selected for backfill material (Dr = 15%, 25% and 40%). Figure 12(a) depicts the computed lateral displacement of the quay wall's head. The horizontal deformation for all relative densities is greater than allowable value proposed by PIANC (2001) and the quay wall system goes toward failure. After the earthquake, the system reaches to equilibrium. The final permanent deformations

is selected as a key parameter to judge about the stability of quay wall.

**-1500 -1000 -500 0 500 1000 1500**

the gravity-type quay wall and surrounding soils.

**0 5 10 15 20 TIME (s)**

**0 5 10 15 20 TIME (s)**

excess pore water pressure, (d) total pressure behind the quay wall

**4. Numerical results and discussion**

**4.1. Influence of relative density of backfill soil**

**-3 -2 -1 0 1 2 3**

**-800 -600 -400 -200 0 200 400 600 800**

**Excess pore water pressure (Pa)**

**ACCELERATION (m/s**

**2**

**)**

(a) (b)

**Figure 9.** Computed post-earthquake (a) deformed shape, (b) nodal displacement vectors for the quay wal after the seismic loading


**Figure 10.** Measured post-earthquake deformed shape after the seismic loading

One might notice that the calculated results are rather close to measured values. This clearly demonstrates that the current numerical procedure captures very well the seismic behavior of the gravity-type quay wall and surrounding soils.

**Figure 11.** Recorded versus computed (a) acceleration, (b) horizontal displacement at top of the quay wall, and (c) excess pore water pressure, (d) total pressure behind the quay wall
