Comprehensive Liquefaction Constitutive Models

The state of practice for seismic analysis involving liquefiable materials is currently experiencing a shift from using empirical schemes (first developed in the 1970s) to simulate liquefaction, to time-marching numerical methods incorporating liquefaction constitutive models currently at various stages of development. There are several ways liquefaction behavior is included in numerical methods, ranging from total-stress empirical schemes to estimate liquefaction conditions (e.g., the UBCTOT model—see Beaty and Byrne 2000), to simple effective-stress shear-volume coupling schemes (e.g., the Finn model, and the “Roth” model—see Dawson et al. 2001), to more comprehensive constitutive models (e.g., the UBCSAND model—see Byrne et al. 1995; and bounding surface models such as the model described in Wang 1990, and the model described in Papadimitriou et al. 2001) that address cyclic shearing via kinematic hardening.

To help practicing engineers choose a procedure best suited to their needs, a selection of approaches is outlined below, ranging from simple to elaborate in terms of complexity and model parameter determination. This is not a comprehensive list, but is a selection that illustrates the different types of liquefaction models that have been developed and used in two-dimensional FLAC, and could be used in FLAC3D. [1] One important factor to keep in mind is that, in engineering practice, the use of a very complicated model for liquefaction analyses is often hardly justified, considering the many uncertainties with respect to soil properties and earthquake motions, and the numerous approximations that must be made (see Dawson et al. 2001). Before describing the individual models, it will be helpful to review the current state of practice for liquefaction analysis. (See Byrne et al. 2006 for further discussion on state-of-practice analysis.)

State of Practice — The standard practice approach for liquefaction analysis of earthquake loading is based on a total-stress analysis in which it is assumed that the liquefiable soil remains undrained at the in-situ void ratio (Byrne and Wijewickreme 2006). Typically, this analysis approach is divided into three steps:

  1. Triggering Evaluation: Typically, an equivalent-linear elastic, dynamic analysis (such as SHAKE) using strain-compatible moduli and damping is conducted for the design earthquake. The cyclic stress ratio (CSR) [2] is evaluated from the numerical simulation and compared to the value of the cyclic stress resistance that the soil has because of its density (CRR) [3], derived from empirical curves. A factor of safety against triggering liquefaction is evaluated using the ratio of CRR and CSR (e.g., see Byrne and Anderson 1991, and Youd et al. 2001).

  2. Flow Slide Assessment: Post-liquefaction (undrained) strengths are assigned in zones predicted to liquefy from the triggering evaluation analysis, and a standard limit-equilibrium analysis is carried out to evaluate the factor of safety against a flow slide. Post liquefaction strengths may be derived from penetration resistance (corrected blow count \((N_1)_{60}\)) using empirical charts (e.g., see Seed and Harder 1990, and Olson and Stark 2002).

  3. Seismic Displacements: Displacements are evaluated using the Newmark approach (see Newmark 1965). In this step, the potential sliding block of soil is modeled as a rigid mass resting on an inclined plane. The design time history of acceleration is applied at the base, and the equation of motion is solved to obtain the displacement of the mass caused by the shaking.

The main shortcomings of the standard practice approach are that the three aspects of liquefaction (triggering, flow slide, and deformation) are treated sequentially, when in reality they may interact locally in various zones of the soil structure and affect the overall behavior of the soil mass. Also, no direct account is made of excess pore-pressure redistribution and dissipation.

Total-Stress Synthesized Procedure

The synthesized procedure of Beaty and Byrne (2000) uses FLAC and the UBCTOT constitutive model to combine the three steps of the standard practice approach (triggering, flow slide, and estimate of liquefaction-induced displacements) into one single analysis. The procedure, which assumes undrained behavior, uses a (two-dimensional) total-stress approach to liquefaction analysis and relies on adjustment of liquefied element properties (stiffness and strength) at the instant of triggering of liquefaction. The main features of the UBCTOT model are summarized below.

A seismic analysis using the synthesized procedure starts from a static state of equilibrium for the FLAC model. The seismic analysis is conducted in total-stress space. UBCTOT uses Mohr-Coulomb elasto-plastic logic with zero friction and a value of cohesion equal to the undrained shear strength, in combination with Rayleigh damping. The elastic shear modulus is assigned a value of \(G_{max}\) multiplied by a modulus reduction factor (MRF). Unlike equivalent-linear methods, this approach is not iterative, and appropriate values of MRF and damping are selected at the start of the seismic analysis. Triggering of liquefaction is based on changes of shear stress on the horizontal plane, \(\tau_{xy}\). The irregular shear-stress history caused by the earthquake is interpreted in each FLAC zone as a succession of half-cycles with the contribution to triggering determined by the maximum value of \(\tau_{cyc}\), defined as the difference between \(\tau_{xy}\) and the initial horizontal shear-stress prior to earthquake loading (i.e., the static bias). A cumulative damage approach is used to combine the effects of each half-cycle. The approach converts the nonuniform history into an equivalent series of uniform stress cycles with amplitude equal to \(\tau_{15}\) (i.e., the value of \(\tau_{cyc}\) required to cause liquefaction in 15 cycles, which is approximately the number of cycles in a magnitude 7.5 earthquake). This is done using an empirical chart giving the cyclic stress ratio (CSR) versus cycles to liquefaction. Several property changes are imposed when liquefaction is detected in a FLAC zone: a residual shear strength is assigned; a reduced loading stiffness is used; and unloading uses a stiffer modulus than loading, according to a bilinear stiffness model. Also, a hydrostatic stress state is imposed when a zone experiences a shear-stress reversal. Finally, reduced viscous damping is assigned in a liquefied zone. The constitutive model also has logic to account for anisotropy in stiffness and strength. See Beaty (2001) for additional information.

The UBCTOT model removes some of the limitations associated with the sequential approach to problem solving used in the state-of-practice procedure, while relying on similar empirical charts for triggering of liquefaction and residual strength.

The following are some of the drawbacks of the model:

  • the use of equivalent modulus ratio that may not capture the pre-liquefaction phase well;

  • the cyclic shear stresses are accounted for on the horizontal plane only;

  • the simplified manner in which the undrained shear strength is specified;

  • pore pressure is not taken into account explicitly; and

  • liquefaction due to monotonic loading is not considered.

Loosely Coupled Effective-Stress Procedure

The Roth model is a loosely coupled effective-stress constitutive model to generate pore pressure from shear stress cycles using the Seed cyclic stress approach (Seed and Idriss 1971). This model is built around the standard FLAC Mohr-Coulomb model. The (two-dimensional) model counts shear stress cycles by tracking the shear stress acting on horizontal planes (\(\tau_{xy}\)) and looking for stress reversals. The cyclic stress ratio (CSR) of each cycle is measured, and this is used to compute the incremental “damage” that is then translated into an increment of excess pore pressure. The procedure is “loosely coupled” because pore pressures are only computed after each 1/2 cycle of strain or stress as the analysis proceeds. The model incorporates residual strength (a critical parameter for seismic stability analyses) by using a two-segment failure envelope consisting of a residual cohesion value and zero friction angle that is extended to meet with the traditional Mohr-Coulomb failure envelope.

The model is simple, robust, and practice-oriented; it is based on the widely accepted cyclic-stress approach with input parameters readily obtainable from routine field investigations. (Note that liquefaction due to monotonic loading is not considered.) A disadvantage of the model is that liquefaction-induced consolidation settlements are not captured, because the actual physical mechanism of liquefaction, whereby pore pressure is generated through contraction of the soil skeleton, is bypassed. The model is applicable to problems where slope movements due to reduced shear strength are the main concern (such as seismic stability of dams and waterfront retaining structures), while shaking-induced consolidation settlements are of secondary importance. See Roth et al. (1991), Inel et al. (1993), Roth et al. (1993), and Perlea et al. (2008) for some field applications.

The Roth model is similar to the built-in Finn model, which is also considered a loosely coupled effective-stress model. The primary difference is that in the Finn model, the volumetric strains induced by cyclic loading are evaluated based on an experimental curve of irrecoverable volumetric strain versus number of constant amplitude cycles. Pore pressures are then generated from these volumetric strains, as well as from contraction of the soil skeleton. Also, the Finn model in FLAC (and FLAC3D ), at present, does not include a post-liquefaction residual strength.

Fully Coupled Effective-Stress Procedure

The UBCSAND model is a fully coupled (kinematic hardening) effective-stress constitutive model to predict seismic response and liquefaction of cohesionless soils in plane strain problems. The model uses an elasto-plastic formulation based on an assumed hyperbolic relation between stress ratio and plastic shear strain, similar to the Duncan and Chang (1970) formulation. It is applicable for monotonic as well as cyclic loading (e.g., see Byrne et al. 2003, 2006).

The model implementation is a modified form of the built-in Mohr-Coulomb model in FLAC that accounts for a strain-hardening frictional behavior, neglects cohesion, and applies to plane strain conditions. The hardening law is a hyperbolic function of plastic shear strain. Unloading is assumed to be nonlinear elastic, with bulk and shear modulus as functions of mean (in-plane) effective stress. Stress reversal is detected by a change of sign in horizontal shear stress, \(\tau_{xy}\). Reloading is elasto-plastic, with the yield locus reset to the value at the reversal point. Plastic flow is nonassociated; the logic is based on a variation of Rowe’s stress-dilatancy theory. According to this theory, there is a constant-volume stress ratio, \(\phi_{cv}\), below which the material contracts (i.e., for mobilized friction, \(\phi_m\), smaller than \(\phi_{cv}\)), while for higher stress ratios (i.e., for \(\phi_m > \phi_{cv}\)), the material dilates. The effect of relative density is addressed through the choice of material properties. Most properties are calibrated to field experience as well as centrifuge tests and are conveniently related to blow count, \((N_1)_{60}\).

The model is able to capture the stiff pre-liquefaction stage, the onset of liquefaction at the appropriate number of cycles, and the much softer post-liquefaction response observed in cyclic, simple shear-constant volume tests.

The coupled effective-stress approach corrects many drawbacks of the previous approaches. Although most parameters are related to blow count and rely on a growing body of data and experience, it is always good practice to check on model parameters for each layer (using numerical simulation of a simple shear test) to verify that if it has to liquefy in \(N\) cycles according to the field data during dynamic loading, it will. Also, because comparison with standard procedures may not be straightforward, it is recommended that the model be used under the supervision of an experienced practitioner. With time, and with the increase of its usage, the model should become more prominent and be used for problems ranging from simple to complex with little effort.

The primary disadvantages are: 1) the logic for detection of stress reversal is based on horizontal shear stress only; and 2) the present formulation applies only to two-dimensional analysis.

Fully Coupled Effective-Stress Bounding-Surface Procedure

Bounding surface plasticity provides a framework to account for cyclic stress reversal in two and three dimensions (see Dafalias 1986, and Wang 1990). The models developed by Wang (1990) (herein named the WANG model) and Papadimitriou et al. (2001) (herein named the PAPADIMITRIOU model) are two (kinematic-hardening) constitutive models that have been implemented in FLAC based on that logic.

WANG Model — The WANG model is an effective stress, bounding-surface hypoplasticity model for (cohesionless) soil that is capable of reproducing, in detail, typical monotonic and cyclic, drained and undrained, hardening and softening behavior observed in classical laboratory tests on initially dense and loose soils (Wang 1990). (The term hypoplasticity characterizes the dependence of loading and plastic strain-rate directions on stress-rate direction.)

The model formulation includes a noncircular pyramidal failure (bounding) surface, a loading surface, a surface of phase transformation (at which contractive behavior changes to dilative during shearing), and a critical state surface (defining an ultimate state in which the sand deforms at constant volume under constant stress). The three-dimensional effective-stress model requires the specification of 15 constants for a given sand (eight parameters are required for two-dimensional analysis). One disadvantage is that model calibration is an arduous task because most model constants are not related to properties with which the practitioner is familiar, and the body of available parameter data is not yet sufficiently well-developed. Also, comparison to state-of-practice analysis is not straightforward. The relation between cyclic stress ratio, number of cycles to liquefaction and normalized blow count can be used to calibrate the model (Wang et al. 2001), but there is no direct relation between empirical rules used in standard practice and theoretical laws used in the model theory. The WANG model is a sophisticated research tool for laboratory-scale experiments. Application of the model to study boundary-value problems at the field scale should probably not be attempted without the advice of the model developer, whose assistance may be required for model calibration, interpretation of results, and support on possible issues with numerical implementation.

PAPADIMITRIOU Model — The PAPADIMITRIOU model is an effective stress, bounding-surface model for loose and dense sand that is based on critical state elasto-plasticity (Papadimitriou et al. 2001, Papadimitriou and Bouckovalas 2002). The model applies to monotonic as well as cyclic loading (in two and three dimensions) of non-cohesive soils under small and large strains. The model uses a kinematic hardening noncircular cone as the (loading) yield surface. In addition to bounding and dilatancy (marking the transition between contractive and dilatants behavior) surfaces, the model also contains a critical-state surface (defining an ultimate state in which the sand deforms at constant volume under a constant shear and confining stress). The model contains a total of 14 parameters: 11 of the parameters can be derived from in-situ and laboratory tests, while the remaining three must be derived indirectly via trial-and-error simulations of drained and undrained laboratory tests. (Note that each parameter set is independent of initial and drainage conditions, as well as cyclic shear-strain amplitude.)

The model has the capability to reproduce, qualitatively, the characteristic behavior observed in cyclic experiments, including the degradation of shear modulus and increase of hysteretic damping with cyclic shear strain amplitude, the shear and volumetric strain accumulation at a decreasing rate with increasing number of cycles, and the increase in liquefaction resistance with density. Comparisons with centrifuge experiments have been made (Andrianopoulos et al. 2006a), and the ability of the model to study a practical problem of geotechnical earthquake engineering has been demonstrated (Andrianopoulos et al. 2006b). The disadvantages are the model calibration, which is a rather tedious procedure and requires a test database not readily available in most cases, and the long computational time required for the solution of practical problems. Also, comparison to standard practice is not straightforward. Application of the model to study boundary-value problems at the field scale is not recommended at this time without the assistance of the model developer for model calibration, interpretation of results, and eventual support for issues related to numerical implementation.