Finite Element Modelling of Seismic Liquefaction in Soils

Finite Element Modelling of Seismic Liquefaction in Soils

'Research clepartment,Plaxisbv, Delft, the Netherlaticts

^Geo-Engineering section. Delft University of Technology, Delft, the Netherlands E-mail: v.galavi@plaxis.nl

A B S T R A C T : Numerical aspects of seismic liquefaction i n soils as implemented in the finite element code, PLAXIS, is described in this paper. After description of finite element equations of dynamic problems, three practical dynamic boundary conditions, namely viscous boundary tractions, tied degrees of freedom and free field elements are reviewed. Possibilities and limitation of each type of boundary condition is highlighted. The formulation of a constitutive model, called as UBC3D-PLM, which describes the mechanical behaviour of soils under cycUc loading is also presented. The model is an extension of the two dimensional UBCSAND model developed at University of British Colombia which utilises isotropic and kinematic hardening mles for primary and secondary yield surfaces to properly take into account accumulation of excess pore water pressure and effect of soil densification during cychc loading. By means of a simplified Rowe's stress-dilatancy theory, the model is capable of modelling liquefaction for different stress paths. It will be presented in this paper how most of the model parameters can be found from the corrected SPT blow count which makes the model easy to use for practical applications. Finally, the model is used for modelling a real boundary value problem and the results are compared with field measurements.

1. I N T R O D U C T I O N

One of the main reasons for the damage of soil stiiictures due to earthquake in saturated condition is liquefaction. Many numerical tools have been developed in the last three decades to assess soil liquefaction. However, prediction of liquefaction is still a challenging task. The simplest way of modelling hquefaction which is still used in practice is done by means of a total stress analysis. I n this type of analysis a dynamic analysis based on total stresses is performed to identify the areas with high potential of liquefaction. By reducing the strength of material to its residual strength and preforming a second total stress analysis, deformation of soil structure will be obtained. In this type of calculation, generation and dissipation of excess pore water pressures are totally ignored in the numerical analysis. The resulted deformation from this type of analysis is always questionable due to oversimplification of the problem. The second type of analysis, which is followed in this paper, is based on an effective stress analysis in which hquefaction is occurred as a result of excess pore water pressure generation. In an effective stress analysis, all input parameters are effective parameters and the total stress is summation of effective stress and excess pore water pressure. The generation of excess pore water pressure is done under undrained conditions by considering volumetric strain and the bulk modulus of water in pores. Although, the effective stress analysis is a natural way of modelling liquefaction, but it needs more information about the mechanical behaviour of soil and the numerical methods used for the analysis. Due to the complex behaviour of soil, a proper constitutive model is needed which should be capable of accumulating volumetric strain and consequently capable of modelling liquefaction. The parameter selection and the use of numerical tools play a significant role in this type of analysis. The main objective of this paper is to show how finite element method as a numerical tool can be used in predicting cyclic liquefaction in soils.

The constitutive model which describes the mechanical behaviour of soil under cycUc loading is UBC3D-PLM model. This model is an attempt of using a simple but powerful plasticity framework to analyse and predict the on set of dynamic induced liquefaction. The UBC3D-PLMis a 3-D extension of the UBCSAND model which was introduced by Puebla et al. The 3D formulation was initially developed and implemented as a user defined soil model in PLAXIS[2] by Tsegaye. The formulation of the model described in this paper is a continuation of the work presented by Petalas et al by adding new features in order to achieve higher accuracy under dynamic and cyclic loading.

After description of the constitutive model, a numerical simulation of a quay wall at Kobe port is presented.

2. F I N I T E E L E M E N T F O R M U L A T I O N

The basic equation for the time-dependent movement of a volume under the influence of a dynamic load is described by:

M - ü - l - C - ü - l - K - u = F (1)

I n which M is the mass matrix; C is the damping matrix; K is the stiffness matrix; u is the displacement vector and F is the dynamic force vector. A widely used method to introduce viscous damping to the system is to utilise Rayleigh damping. According to this method, the damping matrix C is calculated as a linear relationship between the mass matrix M and stiffness matrix K :

C = a^M + j3^K (2)

aji and are the Rayleigh damping coefficients.

Eq. (1) is numerically solved by Newmark integration scheme in PLAXIS, i.e.:

[cgM + c^C + K ] Au =

E^f

where, CQ to c; are the Newmark coefficients which are function of Newmark parameters and time step.

3. D Y N A M I C B O U N D A R Y CONDITIONS

Dynamic boundaries are utilised to reduce the infinite domain (which needs to be analysed) to a finite domain by representing the far-field behaviour of medium. Due to limitations of the finite element method, waves are reflected from both Neumann (prescribed stress) and Dirichlet boundaries (prescribed displacement). Therefore, one of the main tasks in a dynamic analysis is to absorb stresses at artificial boundaries to avoid any reflection of waves (i.e. boundaries of the reduced model). Selection of boundary conditions is usually made based on the problem, accuracy and stability of the boundary condition. I n the following some dynamic boundary conditions are reviewed.

3.1 Tied degrees of freedom

This boundary condition proposed by Zienkiewicz et al is appHed to the lateral sides of a finite element mesh to tie the nodes on the same elevation in order to have the same displacements. It is very stable and can be used with highly advanced soil models.

This type of boundary condition perfectly works i n case of one-dimensional wave propagation but it is unable to absorb the waves propagated from internal sources (loads or waves reflected from stmctures or excavations). Although reflection of waves at boundaries is one of the limitations of this type of boundary condition, it can be used with reasonable accuracy in soils with high damping ratio due to the fact that the reflected waves from internal sources are damped out in the soil material.

3.2 Viscous boundaries

Viscous boundaries can be considered as a Neumann type of boundary conditions where the stresses at boundaries are updated to nullify the reflected stresses. The local absorbing boundary, introduced by Lysmer and Kuhlmeyer, is available in PLAXIS. The idea is to absorb the increment of stress caused by dynamic loads according to the following equations:

(4)

(5)

ll and V are the normal and tangential displacements. (/„ and / are normal and tangential static stresses at the boundary of the main domain respectively.

3.3 Free field boundaries

In a free field boundary, the domain is reduced to the area of interest and the free-field motion is appUed to the boundaries by means of free field elements. A free field element consists of a one dimensional element (in 2D problems) with one way dashpots attached to it (Figure 1). To describe the propagation of wave inside the element, the same mechanical behaviour as the suiTounding soil element in the main domain is used. To absorb the waves reflected from internal stractures (or caused by internal sources), viscous boundaries are considered at the boundary of the main domain (Figure 2 and Figure 3).

Figure 1 Free field elements

(main domain)

viscous boundary oase cluster {external dynamic input)

Figure 2 Free field boundary condition with comphant base (no wave reflection at base)

(extenal dynamic input)

Figure 3 Free field boundary condition with rigid base (wave reflection at base)

The free field motion is transformed from free field elements to the main domain by applying the equivalent forces according to Eq. (6) and Eq. (7). In these equations, effect of viscous boundary condition is also considered at the boundary of the main domain to absorb the waves reflected from the internal structures.

dii'" dt ^llf^^ ' dt + C7' dt ' dt + r (6) (7)

As seen in Figure 2 and Figure 3, free field elements can be attached to the left hand side, right hand side and the bottom of the main domain. I f the base cluster is considered, absorption and application of dynamic input can be done at the same place which is known as compliant base. The stresses in a compliant base are given by: -pVs dii' dii dt dt + cr,: 3 / dt ' dt -\-T (8) (9)

where»" and ii'' are displacements due to upward and downward waves, respectively which can be considered as displacement i n the element and the main domain, respectively.

4. C O N S T I T U T I V E R E L A T I O N

In this section the formulation of the UBC3D-PLM is presented. This model is an extension of the two dimensional UBCSAND model developed at University of British Colombia. The UBC3D-P L M model uses the Mohr-Coulomb yield criterion i n a 3¬ D principal stress space. Moreover, a modified non-associated plastic potential function based on Drucker-Prager criterion is used in order to maintain the assumption of stress-strain coaxiality i n the deviatoric plane for a stress path initiated from the isotropic line. In the proposed model a soil densification nile is added to predict more realistic evolution of excess pore pressures during cyclic loading. A correct procedure of counfing the cycles during dynamic and cyclic loading is implemented i n order to achieve higher accuracy i n the stress paths which do not start from the isotropic line. The densification rale allows the increase of the pore pressures with a decreasing rate during shearing which is observed i n experimental studies.

I n the following the main characterisrics of the model are presented.

4.1 Yield function

The UBC3D-PLM model uses two yield surfaces to model cyclic behaviour of soils, namely primary and secondary yield surfaces (Figure 1). The primary yield surface is based on isotropic hardening and becomes active when the mobiUsed friction angle is equal to the maximum mobilised friction angle that the soil has ever reached. I n this case the current stress ratio is the highest stress ratio in the loading history. A simplified kinematic hardening loile is utilised for the secondary yield surface which becomes active when the mobilised friction angle is less than the maximum mobilised friction angle (Figure 4). This is the case when the current stress ratio is lower than the maximum stress ratio i n the loading history. The distinction between these two yield surfaces is made to be able to have a densification rule (higher rate of hardening) in the secondary yield surface.

The yield surfaces are schematically presented in Figure 4. To explain the effect of both yield surfaces, a stress state is considered on the isotropic axis and both yield surfaces are at the same position. As seen i n the figure, by loading from isotropic stress state both primary and secondary yield surfaces expand according to the same hardening rule. I f the soil is unloaded, the secondary yield surface shrinks and the behaviour of the soil is elastic. By reloading the soil, the secondary yield surface becomes activate and the behaviour will be elasto-plastic again. As soon as the mobilised friction angle reaches the maximum mobilised friction angle, the primary yield surface becomes active again and the behaviour becomes softer.

Failure surface

Q

Primary yield surface ^ - Secondary yieid surface

^ (p'p 1 cp'mob^aï j tp'mob

Effective mean stress, p' Figure 4 Yield surfaces

The well-known Mohr-Coulomb yield function is used to define both yield surfaces. The formulation of the surface is given by:

fm - + c c o t © sm^', (10)

where cf,„ax is the maximum principal effective stress;a',,,/,, is the minimum principal effective stress; c'is the cohesion of the soil; (o'p is the peak effective friction angle of the soil and (p',„„i, is the mobilised friction angle during hardening. The intermediate stress does not influence the yield surface in three dimensional stress space.

4.2 Elastic behaviour

The elastic behaviour which occurs within the secondary yield surface is governed by a stress dependent non-linear rule defined by Puebla et al. The elastic bulk modulus, i i f , and the elastic shear modulus, G^, are stress dependent and are given by the following equations:

G'

(11)

(12)

where Ifg and k'g are the bulk and the shear moduli numbers at a reference stress level, respectively, p' is the mean effective stress.

rik and ng are two parameters to define the rate of stress dependency

of stiffness. The reference stress level, P^, is commonly taken as the atmospheric pressure (/.e.lOOkPa) i n the literature. Pure elastic behaviour is predicted by the model during unloading.

4.3 Plastic potential function

According to the theory of plasticity, the direction of the plastic strain increment (vector) is specified by the gradient of the plastic potential function. In the proposed model, a non-associated flow i-ule based on the Dracker-Prager model is formulated as follows:

6siniy,„„, 3-sin{«^,„

(13)

where g is the plastic potential function; y/,„h is the mobilised dilation angle; cj is the deviatoric stress and /?'is the mean effective stress. It follows from this equation that the direction of the plastic strain increment is perpendicular to the Drucker-Prager surface. This ensures a radial return of stress, i n which Lode angle is constant during the retiun stress point algorithm proposed by Borja and Lee.

The mobilised dilatancy angle, %,„i„ is computed based on the flow rule used by Puebla et al which is derived from the stress-dilatancy theory, linearized and simplified according to energy considerations. The relationship is given as follows:

4.4 Hardening rule

(14)

The hardening rale as introduced by Puebla et al is described as follows:

sin \2

(15)

where Ifg is the plastic shear modulus number; up is a model parameter stands for stress dependency of the plastic shear modulus;

p' is the mean effective stress; is the atmospheric pressure; dA is

the plastic strain increment multiplier and <p'„„ is the ultimate mobilised friction angle to hmit the hyperbolic curve (above equation) and is obtained from the failure ratio, Rf, according to:

- ~

T-smtp,,,,

(16)

The failure ratio, Rf, is always less than 1. 4.5 Densification rule

Based on the mobilised friction angle an unloading-reloading criterion is defined i n the model. A soil densification rule is introduced i n order to have higher accuracy in the predicted evolution of the excess pore pressures following Beaty & Byrne. A secondary yield surface was also introduced i n the model for the secondary loading i n order to ensure a smooth transition into the liquefied state of the soil. The secondary yield surface generates less plastic strains compared to the primary yield surface. Anisotropic hardening rule is used for the primary yield surface, while a simplified kinematic hardening rule is used for the secondary surface.

The plastic shear modulus number hFa during primary loading is identical with the one entered as input parameter. However, during secondary loading it increases as a ftmction of the number of cycles in order to capture the effect of soil densification as follows:

rl> — hP A I "rev

where n is the number of shear stress reversals from loading to unloading or vice versa. Two reversals count for one full cycle.

facdens is a multiplier which is a user input parameter i n order to

calibrate the densification rule. The term koemis a factor between 0.5 and 1.0 to correct the densification lule for loose and non-cohesive soils, i.e. low (A';)6o values.(A'7)(io is the SPT blow counts corrected for overburden pressure and hammer energy ratio, kfcprimaij is the initial value of l^a entered by the user for the primary yield surface.

This rule is the result of calibrating a number of direct simple shear tests. Thus, the calibration factor plays a key role when the user wants to model different stress paths (cychc triaxial tests etc.) and the final value is a matter of judgement according to the most critical stress path for a specific problem. It finally leads to an increase of the excess pore pressure during undrained cyclic loading until the hquefied state is approached. The rate of generation of excess pore pressure decreases by increasing number of cycles which is proven via experiments.

4.6 Post-liquefaction and cyclic mobility

One important issue during the modelling of cyclic liquefaction in sands is the volumetric locking. The evolution of the volumetric strains, after the stress path reaching the yield surface defined by the peak friction angle, becomes constant due to the formulation of the flow nile (in Equation (14)f'„,o(, becomes (»',, and remains constant while (p'cv is also constant).

Due to this issue the stiffness degradation of the soil due to the post-liquefaction behaviour of loose non-cohesive soils or due to the cyclic mobility o f dense non-cohesive sands, which is observed in the experimental studies, cannot be modelled. This limitation is solved i n the formulation of the UBC3D-PLM with the implementation o f an equation which gradually decreases the plastic shear modulus kFa as a function of the generated plastic deviatoric strain during dilation of the soil element.

The stiffness degradation is formulated based on the plastic deviatoric strain related with the dilation of the soil element, due to the deconstioiction of the soil skeleton which occurs during dilative behaviour. This leads to the decreased soil stiffness during contraction which follows after the unloading phase. This behaviour is presented in Figure 5 picturing the process of cyclic mobility of a dense sand. The stiffness degradation is computed as follows:

'^G '^G.priimry ^ with

£ , , , = m i n ( l l O £ , , „ / a c ^ „ J

(18)

(19) Where e^n is accumulation of the plastic deviatoric strain which is generated during dilation of the soil element. With the input parameter/flCpoj, the value of the exponential multiplier term i n the above mentioned equation is limited.

Touching the

Poakrield Surface Piimaiy Loading

« -30 ' ^ ' ^ —

0 10 20 30 40 50 60 70 80 90 100 Effective Vertical S t r e s s , a y y ' (kPa)

Figure 5 Undrained cyclic shear stress path reproduced withUBC3D-PLM for dense sand. Cyclic mobility, stiffness degradation and soil densification are mentioned on the graph

4.7 Undrained behaviour

The increment of the pore water pressure is computed by

n

(20)

where K„ is the bulk modulus of the pore fluid, n is porosity and

de^. is the volumetric strain increment. The bulk modulus of the pore

fluid, when soil is fully saturated is derived from:

= K - K ' = 2G'

l-2v.. l-2v'

Ï + V.. (21)

where K„ and K' denote the undrained and drained bulk moduli of the soil, respectively, v' is the drained Poisson's ratio and K„ is the undrained Poisson's ratio which is assumed to be 0.495. The drained Poisson's ratio can be derived from the elastic parameters of the model as: V =-4.7 3 j i : ' ' - 2 G ' ' 6K'+2G' Model parameters (22)

The input parameters used in the UBC3D-PLM model are summarised here. Table 1 gives the input parameters of the model. Most of the parameters have physical meaning and some are fitting parameters which need to be found by fitting experimental curves. The procedure of determining the parameters is presented i n the section where a real boundary value problem is simulated.

Table 1 Input parameters of the UBC3D-PLM model Parameters Unit Description

V'v (degrees) Peak friction angle

(degrees) Friction angle at constant volume

c' (stress) Effective cohesion

k^B (-) Elastic bulk modulus number

k'a (-) Elastic shear modulus number

I f a (-) Plastic shear modulus number

Ilk (-) Power for stress dependency of

elastic bulk modulus

"g (-) Power for stress dependency of

elastic shear modulus

up (-) Power for stress dependency of

plastic shear modulus

Rf (-) Failure ratio

PA (stress) Reference stress

fa^dens (-) Fitting parameter to adjust densification rale

(-) Fitting parameter to adjust post liquefaction behaviour

(N,)60 (-) Corrected SPT blow counts

5. V A L I D A T I O N O F T H E U B C 3 D - P L M I N E L E M E N T T E S T S

5.1 Monotonic loading

The validation of the UBC3D-PLM i n monotonic loading is presented in this section. The behaviour of loose Syncrude sand is modelled with UBC3D-PLMand the numerical results are compared with experimental data as well as numerical results modelled with the original 2D UBCSAND model pubUshed in Puebla et al, 1997. The input parameters for modelling the monotonic triaxial compression test (TxC) and the monotonic direct simple shear test (DSS) on loose Syncrade sand are given i n Table 2 and Table 3 respectively. The results are presented i n Figure 6 and Figure 7. The results of the UBC3D-PLM are i n a good agreement with the experimental and numerical results of the original 2D UBCSAND.

One of the limitations of the cunent implementation is that the effect of the rotation of the principal stresses and its influence in the stiffness cannot be modelled. Thus, the user has to enter a different input parameter for the plastic shear modulus number i n order to properly model the triaxial and the direct simple shear loading path. This issue has to be taken into account during the calibration of the model according to the most critical stress path expected in a specific case.

Table 2 Input parameters used i n validation of triaxial compression test

Parameters Values Unit

33.7 (degrees) <P'cv 33.0 (degrees) k^B 300 (-) k'o 300 (-) k'c 310 (-) Ilk 0.5 _(-) ng 0.5 _(-) lip 0.5 _(-) Rf 0.95 (-) PA 100 (kPa)

Table 3 Input parameters used in validation of direct shear test Parameters Values Unit

<P'n 33.7 (degrees) <P'cy 33.0 (degrees) k'B 300 (-) k'o 300 (-) k"G 98.3 (-) Ilk 0.5 _(-) ng 0.5 (-) lip 0.5 (-) Rf 0.95 _(-) PA 100 (kPa) 70 -60 i UBC3D : J i 1 ' 1

Puebla et al. 1997 (UBCSAND) | B taboratoiy Results (TxC) j

1 i i

Axial Strain. c_1 (%)

Figure 6 Numerical modelling of loose Syncrude sand under undrained triaxial compression. Comparison of the UBC3D-PLMwith the original 2D UBCSAND and experimental data

Figi

0.0% 0.2% 0.4% 0.6% 0.8% 1.0% 1J% 1.4% 1.6% 1.6% Shear Strain, y j x (%)

ure 7 Numerical modelling of loose Syncrude sand under undrained simple shearing

5.2 Cyclic loading

The behaviour of loose Fraser sand under cyclic direct simple shear is modelled and the numerical results are compared with experimental data as pubUshed by Sriskandakumar. The relative density RD of the tested sand is 40 %. Three different shear stress ratios are used {CSR=0.08, 0.1, 0.12) with the same set of parameters. The vertical applied stress is lOOkPa in all cases. The KQ factor in the numerical calculations is assumed to be 0.46 computed with the well-known Jaky formula. Therefore the initial stresses after consolidation in the two horizontal directions equal 46

kPa.

In Figure 8, Figure 9 and Figure 10 the evolution of excess pore pressure in three stress controlled DSS tests are presented. The input parameters are presented i n Table 4. Even though modelling the onset of liquefaction in the framework of classical plasticity is very complicated, the UBC3D-PLM constihitive model can reproduce the evolution of excess pore pressures towards cyclic liquefaction with adequate accuracy for three different shear stress ratios using the same set of parameters. The updated formulation of the densification rule is helping the model not to predict very steep evolution of the excess pore pressures i n the case of the anisotropic consolidated soil.

One of the main hmitations of the model is presented i n the case of the higher shear stress ratio {CSR=0.\2). The formulation of the UBC3D-PLM cannot take into account the an isotropic consolidation effects during the primary loading which causes higher evolution of excess pore pressures during the first full cycle. Due to this issue the UBC3D-PLMpredicts a slower evolution for this CSR. The slope of the curve in the experimental results is much steeper during the first two half cycles as can be seen in Figure 10. The ability of the UBC3D-PLM of reproducing with high accuracy the cyclic stress paths which are started from the isotropic line was presented in Petalas et al., 2012. The influence of induced an isotropy during the process of anisotropic consolidation is an aspect under research and a new formulation for the framework of classical plasticity is needed.

Table 4 Input parameters for validation of cyclic direct shear test Parameters Values Unit

33.3 (degrees) <P'cv 33.0 (degrees) k'B 848 (-) k'a 594 _(-) kPo 243 _(-) nk 0.5 _(-) ng 0.5 (-) lip 0.4 (-) Rf 0.95 _(-) PA 100 (kPa)

In Figure 11 the influence of the post liquefaction formulation can be seen for the case of the lower shear stress ratio (CSR=O.OS). The predicted total shear strains are in good agreement with the experimental results and prove the liquefled state of the soil.

10 12 14 16 18 20 N u m b e r of C y c l e s

Figure 8 Evolution of excess pore pressures during simple shearing on Fraser sand {RD=40%; CSR=O.Q^; ö;,=100kPa)

0 1 2 3 4 5 6 7 Number of C y c l e s

AkikkK..

experienced permanent deformation and severe damages during the Hyogoken-Nambu seismic event.

The main characteristics of the permanent deformations that observed in the quay wall RC-5 in Rokko Island are presented in Figure 12. The wall displaced towards the sea about 3m to 5m, settled about l m to 2m and titled about 4 to 5 degrees average i n lai et ai., 1998 and Inagaki et al., 1996. There was no collapse or overturning of the caisson wall. Inagaki et al, 1996 reported the in-situ observations after investigation by diving underwater in front of the caisson walls. The reported observations for a quay wall located at the North West corner of Rokko Island are presented in Figure 13 as schematized by Inagaki et al., 1996.

Figure 9 Evolution of excess pore pressures during simple shearing on Fraser sand (RD=40%; CSR=0.l; ö;,=100kPa)

Number of C y c l e s

Figure 10 Evolution of excess pore pressures during simple shearing on Fraser sand (RD=4Q%; CSR=O.U; ö;,=100kPa)

-6% -A% -2% 0%

S h e a r Strain, y {%)

Figure 11: Evolution of shear strain during simple shearing on Fraser sand (RD=40%; CS/?=0.08; ö;,=100kPa)

6. B O U N D A R Y V A L U E P R O B L E M 6.1 Description of the case study

A numerical simulation of a caisson type quay wall under dynamic loading is presented i n this section. The investigated site is located in the region of Kobe i n Japan, on the Rokko Island which was shaken by the Hyogoken-Nambu Earthquake i n 1995. The quay walls both of the Port and Rokko Island of Kobe severely suffered from the earthquake and major damages were recorded after the event.

The ground of the site was created by reclaiming decomposed granite called Masado. A typical cross section of the quay walls constructed at Rokko Island is presented in Figure 12. The soil replacement technique was used i n order to improve the bearing capacity of the foundation layer. Thus, the Alluvial Clay was replaced by the Masado sand and foundation mbble was constructed as well. The soil i n both the foundation rubble and the backfill rubble was also improved. The quay walls i n both the Port and the Rokko Island were constructed using the pseudo static method with low seismic coefficients ranging from 0.10 to 0.15. These walls

Figure 12 Cross section of caisson wall RC-5 in Rokko Island and its residual deformation observed after Kobe 1995 earthquake

(from lai etal, 1998)

After the pfgSi-^"-"""^

g S I i l M ^ r Earthguate

S COLt *Q4II

Afttr the Eofttiauol(e Before Ihe Eorll^

_ I Q Q Drifted arid DepMlt»ri gnli mm

'(Unit In meters)

Figure 13 Deformation of mbble mound beneath the caisson (from lai eM/., 1998)

According to the observations the caisson tilted into the foundation mbble and pushed the mbble out in accordance with the displacements of the caisson wall. The conclusions of the in-situ measurements and site investigations indicate that seismic inertia force, hquefaction of the backfill soil and the deformation of foundation soils beneath the caisson affected the performance of the caisson type quay walls. The backfill soil liquefied at the free field approximately 30 meters further from the walls and the reduction of the bearing capacity of the foundation layer due to excess pore water pressure increase was speculated as a main cause of the damage to the caisson wall in lai etal., 1998 and Inagaki et al, 1996.

6.1 Numerical model and parameter selection

The performance of the quay wall during the Hyogoken-Nambu seismic event is simulated with the use of PLAXIS 2D dynamic module. The target of this sUidy is to vahdate the performance of the implemented finite element formulation in PLAXIS, test the influence of the selected boundary conditions and judge the abiUty

of the UBC3D-PLM constitutive model in simulating this complex behaviour in the framework of classical elasto-plasticity.

Two finite element models are developed with different boundary conditions at the sides. In both cases prescribed acceleration is imposed at base which can be considered as a rigid base. I n the first model the free field elements are used in order to model the wave propagation from free field to main domain and in the second model tied degrees of freedom are used as one of the most widely used method in numerical modelling of dynamic problems. The models are presented in Figure 14 and Figure 15 respectively.

Figure 14 Finite element mesh with free field elements at the sides. (1) Quay wall, (2) backfill nibble, (3) backfiU soil, (4) foundation

mbble, (5) foundation soil, (6) Alluvial clay

Figure 15 Finite element mesh with tied degrees of freedom at the sides. (1) Quay wall, (2) backfill rubble, (3) backfill soil, (4)

foundation mbble, (5) foundation soil, (6) Alluvial clay

One of the advantages of the free field elements is that the size of the finite element mesh can be reduced which improves the calculation time. In the case of tied degrees of freedom the geometry has to be mirrored and the problem has to be fully symmetric in order to properly model the motion due to the earthquake. The number of elements in this case is 5438, whereas, i n the case of the free field elements boundary condition is 3171 and the mesh is slightly finer. 6-noded elements are used in both cases. The results of both models are discussed and the influence of the boundary condition is highlighted.

Five material models are used in this case in order to model the different soil layers both during static and dynamic loading. For the sandy soil layers (2, 3, 4 and 5) the Hardening Soil model is used to generate the initial stresses due to gravity loading and the UBC3D-P L M model is used to calculate the undrained dynamic response. For the Alluvial clay layer (6) the Mohr-Coulomb model is used for both the static and dynamic calculation. The concrete wall is modelled with the Linear Elastic model.

When the free field elements are used, drained soil layers with the Mohr-Coulomb model are considered in order to stabilise the lateral boundaries (Figure 14). These layers have equivalent strength and stiffness parameters as the soil layers next to them. This is important when liquefaction occurs next to the free field element column. This is the case i n the specific problem and with this method stabilisation of the side boundary is succeeded.

The model parameters for the soil layers are extracted based on the reported N-SPT values and the Relative Densities of the sandy layers of the site in Inagaki et al., 1996. The relative density was used to extract the parameters for the Hardening Soil (HS) model for the cases of the backfill and foundation rubble {RD=6Q%) and the backfill and foundation soil according to the coiTelations proposed by Brinkgreve et al. The model parameters are presented in Table 5. The concrete caisson is model as a non-porous linear elastic material with Young's modulus {E) of 30000 kPa and Poisson's ratio (v) of 0.2.

For the undrained dynamic calculation the UBC3D-PLM model is used i n order to properly model the evolution of the excess pore pressures i n the sandy soils and capture the onset of hquefaction. The model parameters are extracted based on N-SPT values reported in Inagaki et al., 1996. Beaty et al proposed a set of equations based on the corrected N-SPT value (Nl)6o for the initial generic calibration of the UBCSAND 904aR model. The input parameters

were calibrated i n order to reproduce the liquefaction triggering behaviour recommended by the 1997 NCEERMSF triggering chart corresponds to earthquakes with magnitudes of about 7.5.

Table 5 Material parameters as used i n the HS model for modelling the generation of the initial stresses i n the gravity loading phase

Parameters R D 60% R D 4 0 % (Backfill (Backfill &foundation &foundation rubble) soil) 36000 24000 i?„,/-^(kPa) 36000 24000 (kPa) 108000 72000 (kPa) m 0.51 (-) 0.58 (-)

Cref 1 (kPa) 1 (kPa)

V 35.5 (degrees) 33 (degrees) V 5.5 (degrees) 3 (degrees)

Makra (2013) revised the proposed equations and highlighted the differences between the original UBCSAND formulation and the UBC3D-PLM as implemented in PLAXIS. The results of the generic calibration of the latter are presented in Figure 16 as compared with semi-empirical correlation for the cyclic strength of non-cohesive soils (CRR). The proposed equations for the generic initial calibration are as following [20]:

Ki = 21,7 X 20 X (23) Kl = Ki X 0,7 (24) 4 = A"! X iNjio X 0.003 + 100 (25) W , = f c . + ^^ -«rf ^ , , = 33 (26) , / . . (^'hhii -15 \ (27) <Pp = f p + niaxj^ö.ö, j •nk = O.5., ng = 0.5, n p = 0.4 ( d s f a u k ) (2^) R,.^l.lxCNOsi'' (29) facd,„, = 0.45 (30)

Figure 16 Values of CRR predicted by UBC3D-PLM and compared to semi-empirical relationships (Maki'a, 2013)

The backfill and foundation nibble is modelled with (Nl)eo equals 20, the backfill soil with (Nl)eo equals 10 and finally the foundation soil with (Wj^o equals 13. The input parameters of the IJBC3D-PLM model are presented i n Table 6. The performance of the model with parameters derived from N-SPT values is discussed in the next section. A parametric study for the influence of the

faCpos,is presented to investigate the behaviour of the foundation

after possible failure.

One of the limitations of the proposed soil densification formulation is the stress path dependency of its formulation. The proposed input parameter fac^g„sis characterising the soil densification during cychc simple shearing. This may lead to inaccuracies i n case of modelling different stress paths.

Table 6 Material parameters as used i n UBC3D-PLM for modelhng the undrained dynamic behaviour of the sandy soils Parameters (Ni)6o=20 (N,)fl,=15 (Ni)6„=10

Layer 2,4 Layer 5 Layer 3

<p'„(dsgvees) 36 35 34 g)'„. (degrees) 27 30 33 k'ni-) 823 748 654 k'ci-) 1176 1070 934 k"G(-) 1512 820 380 "k(-) 0.5 0.5 0.5 ngi-) 0.5 0.5 0.5 np(-) 0.4 0.4 0.4 R,(-) 0.7 0.73 0.78 P,(kPa) 100 100 100 0.45 0.45 0.45 facvmti.-) 0.02 0.2 0.02 (Nl)60{-) 20 15 10

6.2 Results of numerical simulation

In this section results of the numerical simulations are presented in which the free field elements are utilised. Effect of different boundary conditions, namely the tied degrees of freedom option and the free field boundaries will be discussed in a separate section. Figure 17 shows the deformed finite element mesh at the end of the earthquake (i.e. /=30 s). As it can clearly be seen, the quay wall significantly moves and settles due to generation of excess pore water pressures that lead to liquefaction i n soil. By moving the wall towards the sea-side, the nibber deforms which forms a heave i n the seabed. Figure 18 shows the horizontal and vertical displacements for a point selected at the upper sea-side corner of the wall. At the end of earthquake, the sea-side corner of the wall moves 3.2m towards the sea and settles 1.1m which are in agreement with the observation where the waU displacements are reported to be 3m-5m horizontally and lm-2m vertically (e.g. lai et al., 1998).

The horizontal displacement for the cross-section A A , which crosses the foundation and the wall, is plotted in Figure 19. It follows from this figure that the foundation deforms significantly during the earthquake which leads to a significant horizontal displacement in the quay wall. This result is shown i n Figure 20 and it can be seen that the movement of the wall is not due to the slippage i n the interface beneath the quay waU and is mainly due to deformation in the foundation.

10 IS Time (s)

Figure 18 Horizontal and vertical displacements at the top corner of the quay wall during the earthquake

-3 -2.5 -2 -1.5 -1 Horizontal displacement (m)

Figure 19 Horizontal displacement along the cross section A A through the wall during the earthquake

X

Figure 17 Deformed mesh at the end of the earthquake (/=30 s)

Figure 20 Deformed mesh around the quay wall at the end of the earthquake (t=30 s)

The main cause of increasing the displacements i n the quay wall is the development of plastic strains in the foundation which can be shown by simulating the problem with different amount of facp^si, namely 0.02, 0.2 and 0.7. The resulting horizontal displacements are plotted i n Figure 21. The facp„„ parameter controls the plastic behaviour of soil at the post liquefaction regime. The more this parameters is, the less plastic strain is generated.

Figure 21 Effect of post liquefaction parameter,/aCpoj», on the horizontal displacement at the top corner of

the quay wall during the earthquake

Although the computed displacements are in agreement with the field measurements, the predicted rotation of the quay wall (Figure 22) is far from the observed data. As seen i n the figure the angle of tilting reaches almost 0.8 degrees at 1 Is and then decreases to 0.3 at 30s while the observed angle of tilting is about 4 degrees, lai et a!., 1998. Such a poor prediction has also been reported by other researcher (e.g. Dakoulas et al., 2008) which might be due to the lack of enough information about the stiffness of soil in the foundation, the mbber and the layer behind the wall as reported in Dakoulas e / ö / . , 2 0 0 8 .

Figure 22 Rotation of the quay wall during the earthquake

The excess pore water pressures are plotted i n Figure 23. To highlight the liquefied materials, the excess pore water pressure ratio,

PPR, at the end of the earthquake are plotted i n Figure 24. PPR is

defined as:

PPR = Ap,,. _ P'-P'o (31)

Po Po

where, Ap^,. is the excess pore water pressure, p' is the mean effective stress and p'g is the initial mean effective stress (i.e. in the beginning of the dynamic analysis). I n the case in which the total stress is constant, the soil is liquefied when PPR is equal to 1.0.

The maximum values of PPR reached during the dynamic analysis are plotted i n Figure 25 which shows that the soils are mainly hquefied i n the far filed and not beneath the foundation. This is due to the fact that the layers just behind the wall and beneath the wall are denser and consequently they have higher stiffness. These results are i n agreement with the observation.

6.3 Effect of boundary condition

Effect of dynamic boundary conditions is highhghted i n this section. A l l results presented i n the previous section obtained from the simulations i n which two lateral free field boundaries together with a rigid base (prescribed acceleration) are considered. Since the free field boundaries are sensitive to liquefaction, as soon as they get

liquefied, there is not enough resistance at boundary. This may lead to some difficulties i n the numerical convergence. One possibility to prevent the generation of liquefaction at boundary is to consider drained behaviour in the soil elements attached to the free field elements. In a drained element, the bulk modulus of water is assumed to be zero which prevents generation of excess pore water pressure. Another way of modelling this kind of problems is to consider tied degrees of freedom boundaries in which lateral boundaries are tied to each other to produce some confined stresses at boundaries. Figure 26 shows the deformed mesh at t=30 s (at the end of the earthquake).

700 00 600,01 503 OO KO.W 3» to IMOO OW -ICOOO •IWW •TCO-OO -mjy3 •lOCO-CO

Figure 23 Excess pore water pressures at the end of the earthquake (/=30 s)

T

-3

I

Figure 24 Excess pore water pressure ratio (PPR) at the end of the earthquake (/=30 s)

Figure 25 Ivlaximum excess pore water pressure ratios (PPR,„a,^) during the earthquake

f '

Figure 26 Deformed mesh at the end of the earthquake (1=30 s) when tied degree of freedom boundaries are used

The computed horizontal displacements at the upper sea-side corner of the quay wall are plotted i n Figure 27 when the free field or tied degrees of freedom boundaries are used for the lateral boundaries of the domain. As seen, when the tied degrees of freedom boundaries are used, less displacement is predicted (2.7m compared to 3.2m). This is due to short distance between the quay wall and the boundary of the model with tied degrees of freedom and the constraint apphed at that boundary which does not allow the foundation to move freely. This restriction does not exist in the model simulated with the free field elements. Another difference between these two types of boundary conditions is that the tied boundaries do not prevent any wave reflection while i n the case of free field element boundaries are capable of absorbing waves. When the distance between the area of interest and the boundaries is too short, the results might be affected by the restriction considered in the tied boundaries.

E -1

•—^Free field boundaries — Tied degrees of freedom

0 5 10 15 20 25 30 Time (s)

Figure 27 Horizontal displacements at the top corner of the quay when different type of boundary condition is used

7. C O N C L U S I O N

This paper provides an overview of the basic features available in PLAXIS in modelling dynamic problems, especially dynamic liquefaction. Basic finite element formulation of dynamic problems, dynamic boundary conditions and a constiUitive model for modelling liquefaction in soils are presented. The constitutive model is simplified but capable of modelling the onset of liquefaction i n sands which can be used to predict the monotonic and cyclic behaviour of sands with high accuracy. The quay wall i n the region of Kobe i n Japan on the Rokko Island has numerically been simulated to demonstrate all above mentioned features and the capabiUties of the model in predicting the mechanical behaviour of liquefied soils. Since one of the main difficulties in using numerical tools in practice is the determination of the model parameters, by the modelling the quay waU, it is highlighted how the model parameters can be found from the experimental data.

8. R E F E R E N C E

Beaty, M . , Byrne, P. (1998) " A n effective stress model for predicting liquefaction behaviour of sand", Geotechnical Earthquake Engineering and Soil Dynamics I I I ASCE Geotechnical Special Publication No.75, 1:766-777.

Beaty, M . H . , and, P. M . Byrne (2011), UBCSAND constitutive model. Version 904aR, UBCSAND Constitutive model on Itasca U D M Web Site, pp. 69.

Borja, R.I. and S.R. Lee (1990), "Cam clay plasticity, part I : Implicit integration of constitutive relations". Computer Methods in Applied Mechanics and Engineering, V o l 78, pp 49-72.

Brinkgreve, R.B.J., E. Engin, H.K. Engin (2010), Validation of empirical formulas to derive model parameters for sands. In T. Benz, S. Nordal (eds.), 7* European Conference Numerical Methods in Geotechnical Engineering. N U M G E 2010, Trodheim, volume 1, pp. 137-174.

Brinkgreve, R.B.J., E. Engin, and W . M . Swolfs, (2012) PLAXIS, Finite Element code for soil and rock analysis, user's manual. Website: www.plaxis.nl

Dakoulas, P., and G. Gazetas, (2008),"Insight into seismic earth and water pressures against caisson quay walls", Geotechnique, Vol. 58, No. 2, pp. 95-111.

lai, S., K. Ichn, H . L i u , M . Toshikazu(1998), Effective stress analyses of port structures. Special issue of soils and foundations, Japanese Geotechnical Society, pp. 97-114. Inagaki, H., S. Susumu, T. Sugano, H. Yamazaki, T. Inatomi (1996),

Performance of caisson type quay walls at Kobe Port, Special issue of soils and foundations, Japanese Geotechnical Society, pp. 119-136.

Joyner, W. B., and A.T.F. Chen. (1975), "Calculation of nonlinear ground response in earthquake". Bulletin of the Seismological Society of America, V o l . 65, No. 5, pp. 1315¬ 1336.

Kontoe, S. Development of time integration schemes and advanced boundary conditions for dynamic geotechnical analysis.PhD thesis, Imperial College of Science, Technology and Medicine, London, U K (2006).

Lysmer, J.,andR.L.Kuhlemeyer, (1969), "Finite dynamic model for infinite media ", Journal of the Engineering mechanics Division, ASCE, Vol. 95. No. 4, pp. 859-877.

Makia, A. (2013), Evaluation of the UBC3D-PLM constitutive model for prediction of earthquake induced liquefaction on embankment dams, MSc Thesis, T U Delft.

Martin, G., W. Finn, and H . Seed (1975), "Fundamentals of liquefaction under cyclic loading". Journal of the Geotechnical Engineering Division, ASCE 101(GT5), pp. 423-438.

Newmark, N . M . , (1959), " A Method of Computation for Structural Dynamics", Journal of Engineering Mechanics Division, ASCE, Vol. 85, pp. 67-94.

Petalas, A., V . Galavi, and R.B.J. Brinkgreve, (2012), "Validation and verification of a practical constitutive model for predicting liquefaction in sands". Proceedings of the 22"'' European young geotechnical engineers conference, Gothenburg, Sweden, pp. 167-172.

Potts, D.M., and L . Zdravkovic, (1999), Finite element analysis i n geotechnical engineering, theory. Thomas Telford.

Puebla, H., M . Byrne, and, P. PhilHps, (1997),"Analysis of canlexliquefaction embankments prototype and centrifuge models", Canadian Geotechnical Journal 34, pp. 641-657. Rowe, P. W. (1962), "The stress-dilatancy relation for static

equihbrium of an assembly of particles in contact", Proc. Royal Society of London, Mathematical and Physical Sciences, Vol. 269, Series A , 500-557.

Sriskandakumar, S. (2004). Cyclic loading response of frasersand for validation of numerical models simulating centrifugetests. Master's thesis. The University of BritishColumbia,Department of Civil Engineering.

Tsegaye, A., (2010), "Plaxis liquefaction model", external report. PLAXIS knowledge base: www.plaxis.nl

Zienkiewicz, O.C., N . Bicanic, andF.Q. Shen, (1988),"Earthquake input definition and the transmitting boundary conditions". Proceedings Advances i n Computational Nonlinear Mechanics I , Ed. St. Doltnis., Springer-Veriag, pp. 109-138.