1. Theory and Implementation

1.1. Site-Response Analysis in Geotechnical Earthquake Engineering

Earthquake loads are transferred to the foundations of the structures, in the form of accelerations at the ground-structure interfaces. Therefore, an essential part of earthquake engineering is to determine the surficial ground motion.

In earthquake scenarios, the ruptures of faults generate seismic wave, which is then transmitted and modified by path effects (wave scattering, anelastic attenuation, and geometric spreading) and further site effects when it is transmitted through soil layers before it reaches the ground surface. Site effects are known to have a profound influence on surficial ground motions due to the variabilities in surficial geology and the complexity in its characterization, especially when the site needed to be considered with nonlinearity.

Site response analysis is widely employed to quantify site effects. It simulates the propagation of seismic waves from the bed rock to the ground surface, going through soil deposits.

One-dimensional (1D) site response analyses are simplified methods that are commonly used due to their simplicity and cost efficiency. These simplified methods employ constitutive models based on a variation of the hyperbolic model to represent the initial stress-strain backbone curve.

s³hark has been developed to allow site response analyses to be performed taking into account: (1) the nonlinear anisotropic and hysteretic stress-strain behavior of the soil materials; and (2) the effects of the transient flow of the pore water through the soil strata. The procedures used (field and constitutive equations) are general and applicable to multidimensional situations.

1.2. Basics concepts of dynamics of porous media

Soils consist of an assemblage of particles with different sizes and shapes which form a skeleton whose voids are filled with water and air or gas. The word “soil” therefore implies a mixture of assorted mineral grains with various fluids. Hence, soil in general must be looked at as a one (dry soil) or two (saturated soil) or multiphase (partially saturated soil) material whose state is to be described by the stresses and displacements (velocities) within each phase. There are still great uncertainties on how to deal analytically with partly saturated soils. Attention is therefore restricted in the following to dry and fully saturated soils. The stresses carried by the soil skeleton are conventionally called “effective stresses” in the soil mechanics literature (see e.g., [Ter43]), and those in the fluid phase are called the “pore fluid pressures”.

In a saturated soil, when free drainage conditions prevail, the steady state pore-fluid pressures depend only on the hydraulic conditions and are independent of the soil skeleton response to external loads. Therefore, in that case, a single phase continuum description of soil behavior is certainly adequate. Similarly, a single phase description is also adequate when no drainage (i.e., no flow) conditions prevail. However, in intermediate cases in which some flow can take place, there is an interaction between the skeleton strains and the pore-fluid flow. The solution of these problems requires that soil behavior be analyzed by incorporating the effects of the transient flow of the pore-fluid through the voids, and therefore requires that a two phase continuum formulation be available for porous media. Such a theory was first developed by Biot [Bio55] for an elastic porous medium. However, it is observed experimentally that the stress-strain strength behavior of the soil skeleton is strongly non-linear, anisotropic, hysteretic and path-dependent. An extension to the Biot’s theory into the non-linear anelastic range is therefore necessary in order to analyze the transient response of soil deposits. This extension has acquired considerable importance in recent years due to the increased concern with the dynamic behavior of saturated soil deposits and associated liquefaction of saturated sand deposits under seismic loading conditions. Such an extension of Biot’s formulation [Prevost80] has been developed in recent years. In the extended theory, soil is viewed as a multi-phase medium (Green and Naghdi [GN65], and Eringen and Ingram [EI65][IE67]). General mixture results can be shown through formal linearization of the field and constitutive equations, to reduce to Biot’s linear poroelastic model.

1.3. Balance laws

Darcy’s laws

From experiments conducted in 1857, Darcy found that the specific discharge of a fluid in a porous media is proportional to the loss of water head.

(1.3.1)\[\boldsymbol{q} = -\frac{\kappa}{\mu}(\nabla p - \rho_f\boldsymbol{g})\]

where \(\kappa\) is the intrinsic permeability of the porous material, \(\mu\) is the viscosity of the fluid, and \(g\) is the gravity vector. The value of permeability depends upon the size of the voids. As a first approximation one may consider that the permeability \(\kappa\) is proportional to the square of the particle size.

Effective stress

Stresses applied to a saturated porous media are distributed to the solid skeleton and the pore fluid according to a propotion. The former stresses are responsible for skeletal deformations, this is why they are called effective. Considering that stresses are positive when they are tensile and pressure is positive when it is compressive, the principle of effective stress is written in index notation as

(1.3.2)\[\sigma = \sigma^{\prime} - \alpha p\]

The parameter \(\alpha\) is Biot’s coefficient. In 1941, it is defined as [Bio41]

(1.3.3)\[\alpha = \frac{K}{H}\]

where \(K\) is the drained bulk modulus of the porous media, and \(1/H\) is the poroelastic expansion coefficient. In 1957, this parameter was redefined by Biot and Willis [BW57] in terms of jacketed compressibility and unjackedted compressibility.

(1.3.4)\[\alpha = 1- \frac{\chi_s}{\chi}\]

The jacketed compressibility coefficient, \(\chi\), is the drained compressibility of the porous material. The unjacketed coefficient of compressibility, \(\chi_s\), is the compressibility of the solid phase.

Balance of fluid mass

One of the major principles of the theory of poromechanics is that the mass of the two constituent, solid particles and the fluid, must be conserved. Its formulation is presented as following.

Fig. 1.3.1 Fluid mass conservation of a infinitesimal unit

Consider a porous medium, comprising a solid matrix (an assembly of soil particles, with a continuous pore space). The pore space is filled with a fluid. The average velocity of the fluid is denoted by \(v_f\) and the average velocity of the solid skeleton is denoted by \(v_s\). The densities are denoted by \(\rho_{f}\) and \(\rho_{s}\), respectively, and the porosity by \(n\).

The equations of mass conservation of the infinitesimal unit can be established by considering the flow into and out of the unit (Fig. 1.3.1). The mass of the fluid in the infinitesimal unit with volume \(V\) is \(n\rho_{f}V\). The increment of this mass per unit time can be calculated based on the net flux across the surfaces of the unit. Thus, the fluid mass balance equation is

(1.3.5)\[\frac{\partial (n\rho_{f}V)}{\partial t} + \frac{\partial (n\rho_{f}v_{x})}{\partial x}V + \frac{\partial (n\rho_{f}v_{y})}{\partial y}V + \frac{\partial (n\rho_{f}v_{z})}{\partial z}V + Q_{m}V = 0\]

Rewriting Equation (1.3.5) in vector form leads to

(1.3.6)\[ \frac{\partial (n\rho_f)}{\partial t} + \nabla \cdot (n\rho_{f}\boldsymbol{v}_f) + {Q_m} = 0\]

where \(\boldsymbol{v}_f\) is the vector expression for fluid velocity; \(Q_m\) is the fluid source. In order to decompose \(\frac{\partial (n\rho_f)}{\partial t}\) in Equation (1.3.6), a parameter, fluid compressibility, \(\chi\), is introduced here. Fluid compressibility is related to the change in fluid pressure and the fractional change in fluid volume. Its definition is

(1.3.7)\[\chi = - \frac{1}{V_f} \frac{\Delta_f}{\Delta_p}\]

In a constant fluid, the continuity of mass requires that \(\rho_{fr}V_{fr} = \rho_{fc}V_{fc}\), in which the subscription \(f\) represent ‘fluid’, \(r\) and \(c\) represent the reference and current configuration, respectively. Inserting the differential expressions \(V_f = V_{fr} + \Delta V_f\) and \(\rho_f = \rho_{fr} + \Delta \rho_f\) into the mass continuity equation produces \(\frac{\Delta V_f}{V_{fr}} = - \frac{\Delta \rho_f}{\rho_f}\). In the assumption of small perturbations, i.e., change in \(V\) is neglected, the equation rewrites as

(1.3.8)\[\frac{\partial V_f}{V_f} = - \frac{\partial \rho_f}{\rho_f}\]

Equations (1.3.7) and (1.3.8) facilitate the expression of \(\frac{\partial \rho_f}{\partial t}\) using fluid pressure \(p\)

(1.3.9)\[\frac{\partial \rho_f}{\partial t} = \rho_f \chi \frac{\partial p}{\partial t}\]

Therefore, the fluid mass balance equation (1.3.6) rewrites as

(1.3.10)\[ \frac{\partial n}{\partial t} + n\chi \frac{\partial p}{\partial t} + \nabla \cdot (n\boldsymbol{v}_f) + \frac{Q_m}{\rho_f} = 0\]

Balance of solid mass

In the infinitesimal unit, the mass of the solid skeleton is \(m_s=(1-n)\rho_s V\) and its relative density is \((1-n)\rho_s\). Similar to the mass balance of fluid, the mass balance of solid is

(1.3.11)\[ \frac{\partial [(1-n)\rho_s]}{\partial t} + \nabla \cdot [(1-n)\rho_{f}\boldsymbol{v}_s] = 0\]

Consider the infinitesimal unit is loaded with an isotropic compressive stress \(\Delta P\) under undrained conditions. The mean stresses caused by the loading is

(1.3.12)\[\Delta \sigma = - \Delta P\]

In the following analysis, the load is divided into two stages when applied to the unit. In the first stage, the load only caused an increase of \(\Delta p\) in the pore pressure. In the second, the pore pressure will be kept unchanged, while the pressure on the unit will be increased with the magnitude of \(\Delta P - \Delta p\).

In the first stage, two definitions are useful: The unjacketed bulk compressibility of the solid phase

(1.3.13)\[\chi_s = - \frac{1}{V} \frac{\Delta V}{\Delta p}\]

and the unjacketed fluid compressibility

(1.3.14)\[\chi_{\phi} = - \frac{1}{V_f} \frac{\Delta V_f}{\Delta f}\]

The factional volume change of the solid phase is

(1.3.15)\[\frac{\Delta V_s}{V_s} = \frac{\Delta V}{V_s} - \frac{\Delta V_p}{V_s} = \frac{\Delta V}{(1-n)V} - \frac{n \Delta V_p}{V_p}\]

With Equations (1.3.13) and (1.3.14), the volume change Equation (1.3.15) rewrites

(1.3.16)\[\frac{\Delta V_s}{V_s} = \frac{n \chi_{\phi} \Delta p - \chi_s \Delta p}{1-n}\]

In the second stage, the effective stress applied on the solid skeleton will be \(\Delta \sigma ^{\prime} = (\Delta P - \Delta p)/(1-n)\). The definition of the compressibility of the solid phase \(\chi_{s} = -\frac{1}{V_s} \frac{\Delta V_s}{\Delta \sigma^{\prime}}\) yields

(1.3.17)\[\frac{\Delta V_s}{V_s} = \frac{-\chi_s(\Delta P - \Delta p)}{1-n}\]

Adding Equations (1.3.16) and (1.3.17) and with the help of Equation (1.3.12) yields

(1.3.18)\[\frac{\Delta V_s}{V_s} = \frac{\chi_s \Delta \sigma + n \chi_{\phi} \Delta p }{1-n}\]

Similar to Equation (1.3.8), the continuity of the solid phase requires

(1.3.19)\[\frac{\partial V_s}{V_s} = - \frac{\partial \rho_s}{\rho_s}\]

Combining Equations (1.3.18) and (1.3.19) yields: the constitutive law for the solid phase

(1.3.20)\[\frac{\partial \rho_s}{\partial t} = \frac{\rho_s}{1-n} (-\chi_s \frac{\partial \sigma}{\partial t} -n \chi_{\phi} \frac{\partial p}{\partial t})\]

Inserting Equation (1.3.20) into Equation (1.3.11) gives

(1.3.21)\[-\frac{\partial n}{\partial t} + (-\chi_s \frac{\partial \sigma}{\partial t} -n \chi_{\phi} \frac{\partial p}{\partial t}) +\nabla \cdot \boldsymbol{v}_s = \nabla \cdot (n\boldsymbol{v}_s)\]

Noticing the fluid mass balance Equation (1.3.10) can be expressed with \(\boldsymbol{v}_s\)

(1.3.22)\[ \frac{\partial n}{\partial t} + n\chi \frac{\partial p}{\partial t} + \nabla \cdot [n(\boldsymbol{v}_f - \boldsymbol{v}_s)] + \nabla \cdot (n\boldsymbol{v}_s) + \frac{Q_m}{\rho_f} = 0\]

Substituting Equation (1.3.21) into the fluid mass balance Equation (1.3.22)

(1.3.23)\[\nabla \cdot \boldsymbol{v}_s + n(\chi - \chi_{\phi}) \frac{\partial p}{\partial t} -\chi_s \frac{\partial \sigma}{\partial t} + \frac{Q_m}{\rho_f} = - \nabla \cdot [n(\boldsymbol{v}_f - \boldsymbol{v}_s)]\]

Noticing \(n(\boldsymbol{v}_f - \boldsymbol{v}_s) = \boldsymbol{q}\) and \(\nabla \cdot \boldsymbol{v}_s = \frac{\partial \epsilon}{\partial t}\) is the time derivative of the volumetric strain, Equation (1.3.23) becomes

(1.3.24)\[\frac{\partial \epsilon}{\partial t} + n(\chi - \chi_{\phi}) \frac{\partial p}{\partial t} -\chi_s \frac{\partial \sigma}{\partial t} + \frac{Q_m}{\rho_f} = - \nabla \cdot \boldsymbol{q}\]

Substituting the stress balance Equation (1.3.2) and \(\sigma^{\prime} = \frac{\epsilon}{C}\) into Equation (1.3.24),

(1.3.25)\[\alpha\frac{\partial \epsilon}{\partial t} + S \frac{\partial p}{\partial t} + \frac{Q_m}{\rho_f} = - \nabla \cdot \boldsymbol{q}\]

in which \(S\) is the dimensionless storage coefficient.

Substituting Darcy’s law (1.3.1) into (1.3.25) yields the final form of the fluid mass conversation equation

(1.3.26)\[\alpha\frac{\partial \epsilon}{\partial t} + S \frac{\partial p}{\partial t} + \frac{Q_m}{\rho_f} = \nabla \cdot (\frac{\kappa}{\mu} \nabla p)\]

Note that the divergence of \(\rho_f\boldsymbol{g}\) vanished in (1.3.26).

1.4. Constitutive laws

The constitutive laws that are available in s³hark are listed in Table Table 1.4.1. The derivation of individual constitutive law can be found in the corresponding reference.

Table 1.4.1 Available material models

Material models	Usage	Development Status
ElasticIsotropic	2D/3D	✓
PM4Sand	2D	✓
PM4Silt	2D	✓
PressureIndependMultiYield	2D/3D	✓
PressureDependMultiYield	2D/3D	✓
PressureDependMultiYield02	2D/3D	✓
ManzariDafalias	2D/3D	✓
Borja-Amies	3D	✓

ElasticIsotropic

In the application, choosing Elastic as the material model will bring up the editing tab for this model.

PM4Sand

Choosing PM4Sand as the material model will bring up the editing tab for the PM4Sand model.

The sand plasticity model PM4Sand (version 3) [BZ15] follows the basic framework of the stress-ratio controlled, critical state compatible, bounding surface plasticity model for sand presented by [DM04]. Modifications to the model were developed by [BZ15] to improve its ability to approximate the stress-strain responses important to geotechnical earthquake engineering applications; in essence, the model was calibrated at the equation level to provide for better approximation of the trends observed across a set of experimentally- and case history-based design correlations. The model is shown to provide reasonable approximations of desired behaviors and to be relatively easy to calibrate.

PM4Silt

Choosing PM4Silt as the material model will bring up the editing tab for this model.

PM4Silt [BZ18] is a plasticity model for representing low-plasticity silts and clays in geotechnical earthquake engineering applications. The PM4Silt model builds on the framework of the stress-ratio controlled, critical state compatible, bounding surface plasticity PM4Sand model. Modifications to the model were developed and implemented to improve instability to approximate undrained monotonic and cyclic loading responses of low-plasticity silts and clays, as opposed to those for purely nonplastic silts or sands. Emphasis was given to obtaining reasonable approximations of undrained monotonic shear strengths, undrained cyclic shear strengths, and shear modulus reduction and hysteretic damping responses across a range of initial static shear stress and overburden stress conditions. The model does not include a cap, and therefore is not suited for simulating consolidation settlements or strength evolution with consolidation stress history. The model is cast in terms of the state parameter relative to a linear critical state line in void ratio versus logarithm of mean effective stress. The primary input parameters are the undrained shear strength ratio (or undrained shear strength), the shear modulus coefficient, the contraction rate parameter, and an optional post-strong-shaking shear strength reduction factor.

PressureIndependMultiYield

Choosing PIMY as the material model will bring up the editing tab for this model. PressureIndependMultiYield material [GCEY09] is an elastic-plastic material in which plasticity exhibits only in the deviatoric stress-strain response. The volumetric stress-strain response is linear-elastic and is independent of the deviatoric response. This material is implemented to simulate monotonic or cyclic response of materials whose shear behavior is insensitive to the confinement change. Such materials include, for example, organic soils or clay under fast (undrained) loading conditions.

PressureDependMultiYield and PressureDependMultiYield02

Choosing PDMY or PDMY02 as the material model will bring up the editing tab for these models.

PressureDependMultiYield and PressureDependMultiYield02 materials [YEP03] are elastic-plastic material for simulating the essential response characteristics of pressure sensitive soil mate- rials under general loading conditions. Such characteristics include dilatancy (shear-induced volume contraction or dilation) and non-flow liquefaction (cyclic mobility), typically exhibited in sands or silts during monotonic or cyclic loading.

Manzari-Dafalias

Choosing Manzari-Dafalias as the material model will bring up the editing tab for this model.

Manzari Dafalias [DM04] is a stress-ratio controlled, critical state compatible, sand plasticity model. A fabric-dilatancy related quantity, scalar valued in the triaxial and tensor valued in generalized stress space, which is instrumental in modeling macroscopically the effect of fabric changes during the dilatant phase of deformation on the subsequent contractant response upon load increment reversals, and the ensuing realistic simulation of the sand behavior under undrained cyclic loading. The dependence of the plastic strain rate direction on a modified Lode angle in the multiaxial generalization enables it to produce realistic stress-strain simulations in non-triaxial conditions. A very systematic connection between the simple triaxial and the general multi-axial formulation makes it possible to use correctly the model parameters of the former in the implementation of the latter.

Borja-Amies

Choosing J2Bounding as the material model will bring up the editing tab for this model.

Borja-Amies [BA94] is a total stress-based bounding surface plasticity model for clays developed to accommodate multiaxial stress reversals. The model is constructed based on the idea of a vanishing elastic region undergoing pure translation inside a bounding surface, and an interpolation function for hardening modulus which varies with stress distance of the elastic region from the unloading point. Central to the development of the model are the general criteria for loading and unloading, which are phrased based upon the simple argument that with continued loading the hardening modulus should decrease monotonically with deformation. Combined with numerical integration of the elastoplastic constitutive equations in a form suitable for a robust computer implementation, the model is applied to cohesive soils undergoing undrained stress reversals and cyclic loading. With a suitable choice of the interpolation function for the hardening modulus, it is shown that existing one-dimensional nonlinear laws for soils can be replicated, such as the hyperbolic, exponential, the Davidenkov, and even the Ramberg-Osgood models. Specifically, the appropriateness of the exponential hardening function for cohesive soils is investigated and its parameters determined for some clays and silts for use in dynamic soil-structure interaction modeling.

[Bio41]

Maurice A Biot. General theory of three-dimensional consolidation. Journal of applied physics, 12(2):155–164, 1941.

[Bio55]

Maurice A Biot. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of applied physics, 26(2):182–185, 1955.

[BW57]

Maurice A Biot and DG Willis. The elastic coeff cients of the theory of consolidation. J Appl Mech, 15:594–601, 1957.

[BA94]

Ronaldo I Borja and Alexander P Amies. Multiaxial cyclic plasticity model for clays. Journal of geotechnical engineering, 120(6):1051–1070, 1994.

[BZ15]

(1, 2) R. W Boulanger and K Ziotopoulou. Pm4sand (version 3): a sand plasticity model for earthquake engineering applications. Center for Geotechnical Modeling Report No. UCD/CGM-15/01, Department of Civil and Environmental Engineering, University of California, Davis, California, 2015.

[BZ18]

Ross W Boulanger and Katerina Ziotopoulou. Pm4silt (version 1): a silt plasticity model for earthquake engineering applications. Report No. UCD/CGM-18/01, Center for Geotechnical Modeling, Department of Civil and Environmental Engineering, University of California, Davis, CA, 108 pp., 2018.

[DM04]

(1, 2) Yannis F Dafalias and Majid T Manzari. Simple plasticity sand model accounting for fabric change effects. Journal of Engineering mechanics, 130(6):622–634, 2004.

[EI65]

A Cemal Eringen and John D Ingram. A continuum theory of chemically reacting media—i. International Journal of Engineering Science, 3(2):197–212, 1965.

[GN65]

Albert Edward Green and Paul Mansour Naghdi. A dynamical theory of interacting continua. International journal of engineering Science, 3(2):231–241, 1965.

[GCEY09]

Quan Gu, Joel P Conte, Ahmed Elgamal, and Zhaohui Yang. Finite element response sensitivity analysis of multi-yield-surface j2 plasticity model by direct differentiation method. Computer methods in applied mechanics and engineering, 198(30-32):2272–2285, 2009.

[IE67]

John D Ingram and A Cemal Eringen. A continuum theory of chemically reacting media—ii constitutive equations of reacting fluid mixtures. International Journal of Engineering Science, 5(4):289–322, 1967.

[Prevost80]

Jean H Prévost. Mechanics of continuous porous media. International Journal of Engineering Science, 18(6):787–800, 1980.

[Ter43]

Karl Terzaghi. Theoretical soil mechanics. johnwiley & sons. New York, pages 11–15, 1943.

[YEP03]

Zhaohui Yang, Ahmed Elgamal, and Ender Parra. Computational model for cyclic mobility and associated shear deformation. Journal of Geotechnical and Geoenvironmental Engineering, 129(12):1119–1127, 2003.