Soil Dynamic Constitutive Model for Characterizing the Nonlinear-Hysteretic Response

Abstract

Characterization of nonlinear hysteretic responses plays a significant role in predicting soil behaviors. They are mostly described with either simple empirical functions or complex constitutive models. However, the input parameters lack both a physical basis and robustness, and the use of these models is limited to some typical soils. Therefore, there is a need for a simple but robust model that uses only a small number of physically meaningful parameters. This study proposes explicit formulas to capture different nonlinear hysteretic soil responses, including a constitutive model, backbone curve, tangent shear modulus, secant shear modulus, and damping ratio. In particular, the Davidenkov model, with two physically meaningful parameters, is adopted to assess the constitutive relationships of soils under steady-state cyclic loading. The proposed models are validated with resonant column test (RCT) data (shear modulus and damping ratio). This paper finds that the use of the linear characteristic equation to calculate the shear modulus from the resonance frequency in the RCT, which is clearly irrelevant and approximate, produces data interpretation errors.

1. Introduction

After earthquakes, detailed studies of the seismic ground responses are typically conducted, showing that (1) soil deposits exhibit strong dynamic nonlinear responses, and (2) characterization of the soil dynamic properties plays a crucial role in accurately estimating the site amplification characteristics [1,2,3].

Strong ground motion causes the accumulation of plastic deformation during earthquake cycles, and elastic theory cannot be applied to describe the deformational characteristics of the soil deposits (i.e., shear modulus decreases and damping ratio increases with the amplitude of the shear strain). Site response analysis methods with nonlinear responses can be divided into the equivalent linear and nonlinear approaches. The first approach, the equivalent linear approach, approximates the nonlinear cyclic responses of soil samples. The effective shear strain (equal to approximately 65% of the peak strain in the time domain) is iteratively computed by updating the shear modulus and damping, after which the site response is simply estimated in the frequency domain. Given that the appropriate selection of the equivalent soil stiffness and damping for soil layers only represents a particular state in the stress-strain space and not the entire stress-strain evolution, this method fails to predict large plastic deformations of soil columns over the entire duration of a seismic event, and its use is strictly limited to relatively small shear strains or small nonlinearities [4,5,6].

The second approach, in which the soil dynamic characteristics are captured through a nonlinear hysteretic constitutive relationship, can represent the strain-dependent shear modulus and damping ratio. The simplest constitutive relationship uses a model relating the shear stress to the shear strain, whereby the backbone curve is expressed by a hyperbolic function. Several functions with empirical fitting parameters have been proposed to describe the strain-dependent response defined by the backbone curve [7,8,9,10,11]. More sophisticated models for cyclic loading use more fitting parameters so as to precisely reproduce closed hysteresis loops that are influenced by basic soil parameters (plasticity, void ratio, confinement stress) and the imposed shear strain amplitude [12,13,14,15,16,17,18,19,20,21]. While those models are capable of expressing the nonlinear hysteretic response or degree of strain-dependency during steady-state cyclic loading, their input parameters lack both a physical basis and robustness, and the use of such models is quite limited to certain soils. Therefore, a simple but robust model involving a small number of physically meaningful parameters is necessary for the characterization of the hysteretic nonlinear properties (Ockham’s criterion).

This study proposes explicit formulas for capturing certain nonlinear hysteretic soil responses. These include a constitutive model, as well as the backbone curve, tangent shear modulus, secant shear modulus, and damping ratio. The proposed models are examined with resonant column test data (shear modulus and damping ratio). The validity of the linear characteristic equation is addressed during the process of data interpretation.

2. Nonlinear-Hysteretic Constitutive Model

Particulate materials (i.e., rock, sandstone, sand, and sediment) exhibit strong dynamic behaviors, such as end-point memory and a closed loop (hysteresis) in the stress-strain curve, higher harmonic generation in propagating waves, a resonance frequency shift and nonlinear damping in standing waves, and slow dynamics [22,23,24,25,26,27]. These macroscopic nonlinear responses of soils are the consequences of complex frictional mechanisms involving contact between the grains and grain rearrangement under loading-unloading cycles. The variation of the damping ratio is a physical indicator of energy dissipation as quantified by a hysteresis loop at one cycle of deformation.

Intensive efforts have been made in the geophysics community to analytically describe these hysteretic nonlinear behaviors of granular materials. This study adopts the Davidenkov model, which has been experimentally validated, to capture the amplitude-dependent internal friction related to nonlinear hysteretic nonlinear behavior. The phenomenological equation of state was analytically derived, and its parameters were found to have physical meaning. Details can be found in the literature [28]. The equation can be rewritten in terms of the shear component:

$$\tau =G_{\max }\left\{\gamma -\mathrm{sgn}(\dot{\gamma })\frac{\alpha }{n}\left[{\left(\Delta \gamma +\mathrm{sgn}(\dot{\gamma })\cdot \gamma \right)}^n-2^{n-1}{\left(\Delta \gamma \right)}^n\right]\right\}$$

(1)

where G_max is the small-strain shear modulus, γ is the shear strain, ∆γ is the amplitude of the cyclic shear strain, n is the integer exponent related to the model order (n > 1), α is the nonlinear parameter, $\dot{\gamma }\left(=\partial \gamma /\partial t\right)$ is the strain rate, and sgn(x) is the signum function. Numerical solutions that use a tangent formulation involve the tangent shear modulus G_tan, which is defined as the slope of the stress-strain hysteresis loop along the loading path. By taking the derivative of Equation (1), the tangent shear stiffness is obtained:

$$G_{\tan }=\frac{d\tau }{d\gamma }=G_{\max }\left[1-\alpha {\left(\Delta \gamma +\mathrm{sgn}(\dot{\gamma })\cdot \gamma \right)}^{n-1}\right]$$

(2)

It should be noted that the tangent stiffness is not associated with small strain stiffness, as the tangent stiffness G_tan is a mathematical concept that reveals the instantaneous state of fabric changes during a large strain test, while the small strain stiffness G_max is a fabric constant measurement of the stiffness [29,30]. The backbone curve represents the trajectory of the extrema of each hysteresis curve. When the shear strain γ is equal to the strain amplitude ∆γ in Equation (1), the piecewise function of the backbone curve is formulated for the entire strain range as follows:

$${\tau }_{\mathrm{upper}}=G_{\max }\gamma \left[1-\frac{\alpha }{n}{\left(2\gamma \right)}^{n-1}\right] \mathrm{for} \dot{\gamma }>0$$

(3a)

$${\tau }_{\mathrm{lower}}=G_{\max }\gamma \left[1+{\left(-1\right)}^n\frac{\alpha }{n}{\left(2\gamma \right)}^{n-1}\right] \mathrm{for} \dot{\gamma }<0$$

(3b)

The hysteretic nonlinear models involved with two model parameters are simulated under three different strain amplitude levels (Figure 1). The backbone curve consists of two polynomial functions, which intersect at the coordinate origin. Note that, when using the model order of n = 2, the first derivative of each function is equal to the small strain stiffness G_max. In addition, higher imposed strain amplitude produces more energy dissipation and a decreasing slope (the strain-softening effect). Indeed, the hysteretic loops superimposed on the backbone curve follow the Masing rule, which has been experimentally observed in particulate materials under cyclic torsional loading [31].

3. Explicit Formula for the Characterization of the Nonlinear Response

The soil dynamic response is often characterized by the shear modulus and damping ratio of soils using a resonant column and torsional shear (RCTS) testing apparatus. It is a typical practice to express the nonlinear stress-strain behavior of soils in terms of the shear modulus and the damping ratio. The secant type shear modulus (G_sec = τ_upper/γ), which is the slope of the secant drawn from the origin to any specified point on the stress-strain curve, can be easily obtained using the analytical form of the backbone curve in Equation (3a):

$$\frac{G_{\sec }}{G_{\max }}=1-\frac{\alpha }{n}{\left(2\gamma \right)}^{n-1}$$

(4)

Note that the secant modulus is always larger than the tangent modulus, because n is a positive integer. Figure 2 presents the evolution of the shear modulus along the loading path. The individual part in the hysteresis loop is formulated by mechanical features whose maximum shear stress switches between two configurations at the imposed strain amplitude. The phenomenon of discrete memory is illustrated in the A → B → C → D loading path. This loop begins and ends at the same stress-strain point. This process creates bow-tie stiffness in the entire hysteresis loop [24,32]. The results of this simple analysis indicate that the dynamic stiffness has strong dependency on the loading path and imposed strain amplitude. In addition, the secant-type shear stiffness at the given shear strain amplitude is superimposed on its tangent counterpart. It was found that G_sec(∆γ) is equal to the average of the two tangent moduli at $\pm \Delta \gamma$ according to Equations (2) and (4), such that G_sec(∆γ) = G_tan(γ = 0) with the second model order (n = 2).

The inherent viscous characteristics of soils dissipate the energy against the imposed loads. Thus, this relationship allows an estimation of the nonlinear damping ratio with closed hysteresis loop:

$${\zeta }^{NL}(\gamma )=\frac{W_D}{4\pi W_S}$$

(5)

where W_S (= G_sec ∆γ²/2) is the maximum strain energy stored during the cycle and W_D is the closed area related to the hysteresis loop in one cycle of loading. The hysteresis loop is constructed by the upper and lower parts depending on the sign of the strain rate. For example, the upper part is the positive strain rate from +∆γ to −∆γ. The energy dissipation per cycle can be calculated using Equation (1) with two integration sections:

$$\begin{gathered} W_D={\displaystyle \int \tau d\gamma =}{\displaystyle {\int }_{-\Delta \gamma }^{+\Delta \gamma }{G}_{\max }\left\{\gamma -\frac{\alpha }{n}\left[{\left(\Delta \gamma +\gamma \right)}^n-2^{n-1}{\left(\Delta \gamma \right)}^n\right]\right\}d\gamma }+{\displaystyle {\int }_{+\Delta \gamma }^{-\Delta \gamma }{G}_{\max }\left\{\gamma -\frac{\alpha }{n}\left[-{\left(\Delta \gamma -\gamma \right)}^n+2^{n-1}{\left(\Delta \gamma \right)}^n\right]\right\}d\gamma } \\ =\frac{2^{n+1}\left(n-1\right)}{n\left(n+1\right)}{G}_{\max }\alpha {(\Delta \gamma )}^{n+1} \end{gathered}$$

(6)

Thus, the nonlinear damping ratio can be expressed using Equation (5):

$${\zeta }^{NL}(\gamma )=\frac{2^n}{\pi }\frac{\left(n-1\right)\alpha {\gamma }^{n-1}}{n\left(n+1\right)}\frac{G_{\max }}{G_{\sec }}$$

(7)

Theoretically, there should be no energy dissipation below the elastic threshold range for the hysteresis damping function defined by Equation (5). However, even at a very low strain level, some energy dissipation can be measured in laboratory tests. The linear damping ratio ζ^L at very low strain levels has a constant value. At higher strain levels, nonlinearity in the stress-strain relation produces more energy dissipation with an increase in the strain amplitude. Thus, the damping ratio of a granular material is the summation of the linear damping ratio at a small strain level and the nonlinear damping ratio which increases with the strain. The strain-dependent damping function can be rewritten by inserting Equation (4) into Equation (7):

$$\zeta (\gamma )={\zeta }^L+{\zeta }^{NL}(\gamma )={\zeta }^L+\frac{2^n}{\pi }\frac{\left(n-1\right)\alpha {\gamma }^{n-1}}{n\left(n+1\right)}\cdot {\left(1-\frac{\alpha }{n}{\left(2\gamma \right)}^{n-1}\right)}^{-1}$$

(8)

4. Discussion

4.1. Examples

Resonant column test data (shear modulus and damping ratio) sourced from the literature are compared with a model for sand and clay with different confinement stress levels. The parameters of model are selected by least-square fitting of the model to the damping ratio as measured from a resonant column test. The linear damping ratio ζ^L is experimentally constrained with a constant value below elastic threshold regime. It was found that the model order as an indicator of degree of nonlinearity is less sensitive then the nonlinear parameter. Moreover, the values only range from 2 to 3 (41 cases are tested). Previous studies revealed that the exponent of the hysteresis is n = 2 for many polycrystalline metals (copper, zinc, and lead) and particulate materials (rock, sandstone, and sand), and that these materials require quadratic correction by a constitutive stress-strain equation in order to describe the symmetrical hysteresis response [24,33]. Therefore, this study uses the second model order. The data points and fitted models are plotted on the semi-log scale (Figure 3). Higher confinement stress extends the elastic threshold region with lower α values. Clay is expected to be less flexible or less nonlinear than sand and the application of confinement stress would reduce the nonlinearity. These outcomes are reflected in the trends of the α values. These examples demonstrate that the α parameter contains physical, not ad-hoc, fitting parameters.

4.2. Is the Linear Characteristic Equation Valid for Data Interpretation?

The secant shear modulus curve, when computed with the optimal α parameter, deviates somewhat from the experimental results obtained through the conventional procedure, which relies on the linear characteristic equation. The conventional RC test provides both the resonance frequency and damping ratio as a function of the shear strain, with the shear modulus then calculated from the resonance frequency using the characteristic equation for the linear vibration of a column-mass system. Indeed, the hysteretic nonlinear property of the vibration problem, as manifested by the strong strain-dependency, is ignored during the data interpretation process. Therefore, these deviations indicate how much error the use of a linear characteristic equation can introduce into the shear modulus. As an alternative means of describing the amplitude-dependent frequency shift, a nonlinear characteristic equation can be derived by incorporating Equation (1) into the moment balance equation together with the relevant boundary conditions, yet a nonlinear vibration analysis is quite challenging. Given that the resonance frequency shift is inherently related to the increase in the damping ratio, it is recommended that the strain-dependent shear stiffness be estimated using Equation (4), with the optimal nonlinear parameter defined by the measured damping curve.

5. Conclusions

The characterization of the dynamic properties of soils plays a crucial role in accurately estimating the degree of site amplification. While previously suggested models have been proposed to capture the nonlinear-hysteretic during seismic loading, they had physically meaningless or redundant input parameters. This paper employs a two-parameter Davidenkov model to characterize the strain-dependent properties of soil under steady-state cyclic loading. Explicit formulas are derived to describe the backbone curve, tangent shear modulus, secant shear modulus, and damping ratio. The models are tested with the resonant column test data (the shear modulus and damping ratio). The results show that the proposed models can capture the strain-dependent soil dynamic responses and help to infer the nonlinear hysteretic stress-strain relation. The results demonstrate that the shear modulus degradation curve unfortunately includes data interpretation errors because a linear characteristic equation is used to estimate the shear modulus from the resonant frequency. Thus, it is recommended that the strain-dependent shear stiffness be characterized using the models proposed here, which utilize optimal nonlinear parameters.