Soil Dynamic Constitutive Considering Post-Liquefaction Deformation and Reversible Pore-Water Pressure

Abstract

In the seismic response analysis of liquefiable sites, the existing soil dynamic constitutive model is challenging to simulate saturated sand’s post-liquefaction deformation, and the current pore-water pressure buildup model cannot reflect the decrease in the actual pore-water pressure under unloading stress. We aim at these problems to propose a feasible and straightforward time-domain post-liquefaction deformation constitutive model through experimental analysis and theoretical research, consisting of reversible pore-water pressure. According to the dynamic triaxial test data, the regularities of large deformation stress and strain behavior of the saturated sand after liquefaction are obtained, and the corresponding loading and unloading criteria are summarized. Combined with the effective stress constitutive model proposed by the author, a soil dynamic constitutive that can describe saturated sand’s post-liquefaction deformation path is obtained. According to the test results, the model can simulate the deformation of saturated sand during the whole liquefaction process. The self-developed program Soilresp1D realized the dynamic response analysis of the liquefiable site, and the results were compared with the experimental results. It shows that the model based on the effective stress-modified logarithmic dynamic skeleton and post-liquefaction deformation constitutive can be directly applied to the dynamic response analysis of the liquefiable site.

1. Introduction

Some important marine and coastal infrastructures, such as undersea tunnels, cross-sea bridges, offshore oil platforms, artificial reclamation islands, and infrastructure on coral reefs, are located on or within the liquefiable soil layer [1,2,3,4]. Many liquefaction damage phenomena of deep sand (10–26 m) have occurred in recent large earthquakes [4,5,6,7]. Several factors affect the occurrence of liquefaction, such as the load amplitude, soil type, initial shear stress, shear strain amplitude, age, and hydraulic conditions [8,9]. When determining the seismic fortification parameters of such sites, a soil dynamic constitutive model considering sand liquefaction should be established [10,11,12]. Large deformation caused by liquefaction mainly causes the disasters induced by the liquefaction of saturated sand layers [12,13,14,15]. Therefore, selecting a reasonable dynamic constitutive model suitable for saturated sand is essential when studying liquefaction [16,17,18,19]. Some nonlinear or elastic-plastic constitutive models used for liquefaction analysis have been proposed successively [7,20,21,22,23,24,25]. However, these models can only simulate the dynamic response of a small deformation before the initial liquefaction and are not suitable for simulating the dynamic response of a large deformation after the initial liquefaction. The soil’s stress–strain characteristics and boundary value conditions will change when a severe liquefaction and large deformation occur under earthquake action [7,26,27,28]. Elgamal and Yang developed a plasticity-based constitutive model with emphasis on simulating the cyclic mobility response mechanism and associated pattern of shear strain accumulation; this constitutive model is incorporated into a two-phase (solid-fluid), fully coupled finite element code (PDMY02) [22,29,30]. The model incorporates shear-induced contractive, perfectly plastic, and dilative response phases implemented through an appropriate non-associative flow rule. At the same time, Dafalias and Manzari proposed a stress-ratio controlled, critical state compatible sand plasticity model (SANISAND) in the triaxial and generalized stress space [31]. They illustrated the model’s simulative ability by a comparison with the data over an extensive range of pressures and densities.

Following the basic framework of SANISAND, Boulanger and Ziotopoulou presented the sand plasticity model PM4Sand (the model follows the basic framework of the stress-ratio controlled, critical state compatible, and bounding surface plasticity model for sand presented by Dafalias and Manzari (2004) [32]) for geotechnical earthquake engineering applications, which provides reasonable approximations of the desired behaviors and is relatively easy to calibrate. Asgaria applied the model to the numerical simulation of the liquefaction site [33]. Shamoto and Zhang constructed an elastoplastic constitutive model suitable for the dynamic reaction calculation of the sand liquefaction with a large deformation. Still, it relied on sizeable finite element calculation software and had a poor general applicability [7,26,34]. Shong proposes a completely explicit finite element method for solving dynamic u-p equations of fluid-saturated porous media [35,36]. These studies provide a theoretical basis for establishing a unified constitutive model for a post-liquefaction deformation. However, these models are almost realized by empirically reducing the modulus or introducing mathematical techniques lacking the objective physical mechanisms. The description of a time-domain constitutive relationship changes before and after the liquefaction is not perfect.

There are few studies on the time-domain post-liquefaction deformation constitutive for the liquefiable site seismic response analysis. It is necessary to construct a model that can be directly applied to the liquefiable site and describe the development process of pore-water pressure and the post-liquefaction stress–strain relationship. Additionally, in light of that, the main objectives of this study are as follows: (1) summarizing the fundamental law of a large deformation and post-liquefaction stress–strain relationship basic formula; (2) deciphering the loading and unloading criteria for the large liquefaction deformation of saturated sand; (3) combining with the primary stress–strain relationship formula and the loading and unloading criterion, establishing a large deformation constitutive model of saturated sand liquefaction with a simple expression and accessible parameters, and calibrating the parameter; and (4) embedding the large deformation constitutive model into the site seismic response analysis self-compiled program Soilresp1D (based on the boundary conditions, initial conditions, motion equation, and LDSCM, a 1D sites time-domain nonlinear seismic response analysis program Soilresp1D by Microsoft Visual c++ 6.0 platform was developed. The program can easily and quickly calculate the seismic response at any site depth [37].) and developing a set of numerical analysis methods and calculation programs for the large deformation analysis of the liquefiable site.

2. Model Formulations

2.1. Demonstration of Test Results

New methods and protocols should be described in detail, while well-established methods can be briefly described and appropriately cited. In the cyclic shear liquefaction test, when the effective stress of saturated sand is zero for the first time, the soil is liquefied, and the state of the soil is an initial liquefaction, which is called before the liquefaction and after the liquefaction [19,38]. More studies have been conducted on the soil liquefaction tests and post-liquefaction character [39,40,41,42,43]. Based on the pore-water pressure model, which can simulate the reversible pore water pressure proposed by the author, combined with a 1D time-domain dynamic nonlinear constitutive model based on the logarithmic dynamic skeleton, the effective stress constitutive model of saturated sand under a strong ground motion before liquefaction is established [37,44]. However, this model can only simulate the dynamic reaction of a small deformation before sand liquefaction but cannot affect the large deformation after the sand liquefaction.

In this paper, we studied the stress–strain curves of saturated sand before and after liquefaction under the condition of a thecyclic triaxial test. The dynamic triaxial control system of DYNTTS-60KN was used to apply a vibration stress load to saturated coral sand. The changes in the pore water pressure were observed and recorded, or when the shear deformation of the sample reaches a critical value, it is judged that liquefaction occurs. Some cyclic triaxial test was conducted on saturated coral with a relative density of 45%. We performed a stress-controlled test (constant CSR) until reaching initial liquefaction and then switched to a strain-controlled test. Before reaching initial liquefaction, the effective confining pressure was 100 kPa, CSR was 0.325, 0.3, 0.25, 0.22, and the loading frequency was 0.5 Hz, respectively. The results are shown in Figure 1, Figure 2, Figure 3 and Figure 4, where (a) is the overall stress–strain curve; (b) is the time history of the pore pressure change; (c) is the stress–strain curve before liquefaction; (d) is the stress–strain curve after liquefaction; and (c) and (d) are obtained by the decomposition of (a). According to the test results, the cyclic loads required for sand to reach the initial liquefaction state under various CSRs are 11, 35, 37, and 85 cycles, respectively.

2.2. Analysis of Test Results

As can be seen from the overall stress–strain curve in Figure 1a, Figure 2a, Figure 3a and Figure 4a, with the increase in the number of cycles, each cycle’s maximum shear strain amplitude gradually increases [24,40]. In contrast, the shear modulus gradually decreases, and the shear stress amplitude of the soil in the initial liquefaction state gradually decreases. The time history of the pore-water pressure shown in Figure 1b, Figure 2b, Figure 3b and Figure 4b shows that the pore-water pressure in saturated coral sand increases monotonously during a continuous accumulation and fluctuates periodically with the transformation of shear stress. When the pore-water pressure reaches the confining pressure, the soil presents zero effective stress and an initial liquefaction.

Figure 1c,d, Figure 2c,d, Figure 3c,d and Figure 4c,d, respectively, give the stress–strain curves of the soils before and after liquefaction. Because the mechanical properties of saturated sand are inconsistent before and after liquefaction, the shape of the stress–strain curves is also different. By comparing Figure 1c,d, Figure 2c,d, Figure 3c,d and Figure 4c,d, it shows that the skeleton curve of the upper half before liquefaction is convex upward, while that of the upper half after liquefaction is concave downward. Figure 1c, Figure 2c, Figure 3c and Figure 4c shows that the shear modulus and shear strength are related to effective stress due to compressive hardening characteristics. Before the initial liquefaction, the shear modulus decreases gradually as the effective stress decreases.

Figure 1d, Figure 2d, Figure 3d and Figure 4d shows that the soil properties changed after liquefaction, and the soil behaves like a fluid. After liquefaction, the pore-water pressure still fluctuates, and a zero effective stress state appears twice in one cycle. According to the cyclic stress–strain curves after liquefaction, it can be divided into unloading and loading states. After liquefaction, the initial shear modulus decreases significantly with the number of cycles, which is more evident than before liquefaction, and the cyclic shear strain peak of the soil after liquefaction increases significantly. Unlike before the liquefaction, the peak shear stress of each cycle of soil decreases gradually after liquefaction, and the peak shear stress of the soil cannot reach the value of the stress applied. The experimental phenomenon above is typical according to the undrained cyclic triaxial test results of saturated coral sand under different CSR conditions in Figure 1, Figure 2, Figure 3 and Figure 4.

Based on the above analysis, the variation characteristics of the shear modulus, shear stress, and shear strain of the soil before and after the liquefaction are summarized as follows: (1) before liquefaction, the initial shear modulus of each cycle decreased with the decrease in the effective stress. After liquefaction, it decreased significantly with the increase in the cycles. (2) Before liquefaction, the soil’s cyclic shear stress amplitude was constant and did not increase with the increase in the cyclic times. After liquefaction, it decreases significantly with the increase in the cycles. (3) Before liquefaction, the amplitude of the cyclic shear strain increased with the decrease in the effective stress. After liquefaction, it increases significantly with the increase in the cyclic number. (4) The skeleton curve in the upper half of the stress–strain curve is convex before liquefaction and concave after liquefaction, which shows an antisymmetric relationship.

Therefore, the abrupt change in the stress–strain relationship occurs when the effective stress reaches zero states for the first time. The key to simulating the stress–strain response is to judge the soil’s liquefaction state and the stress–strain curve’s development law. The above is only a limited description of the experimental results, which need to be expressed by a certain functional expression or mathematical model to establish a large deformation constitutive model.

2.3. Determine Loading and Unloading Criteria

The stress–strain curves of saturated sand during the whole liquefaction process can be divided into two stages: small deformation before liquefaction and large deformation after liquefaction. The stress–strain relationship of small deformation before the liquefaction has been described in the literature [37,44]. The effective stress constitutive model can simulate the reversible pore-water pressure, combined with a 1D time-domain dynamic nonlinear constitutive model (LDSCM) based on the logarithmic dynamic skeleton. It shows as in Equations (1) and (2), Equation (1) is the effective stress model; Equation (2) is the pore-water pressure model. The LDSCM LSC-based introduces the concepts of a modified dynamic skeleton curve and the damping ratio degradation coefficient, determined by the test data G/Gmax-γ and ζ-γ of soil ($a,b$, $a_0, b_0$), which can fully consider the damping effect in the soil’s dynamics problems. The pore-water pressure model shown in Equation (2), containing reversible pore-water pressure, can reflect the increase in the pore-water pressure under stress loading and its decrease under stress unloading, and ($c_{1,0}, c_{1,a}, c_{1,b}, A_{4,0})$ are main material constants of the model, which are related to the density of the soil. The derivation process of Equations (1) and (2) is rather complicated. Equation (1) can be seen from Equations (1)–(12) in [37], and Equation (2) can be seen from Formulas (1)–(4) in [44]. Since this part is not the focus of this article, it will not be repeated here.

$$\left\{\begin{array}{c}\tau \left(\gamma \right)=\{\begin{array}{c}K\left({\gamma }_0\right)·[\pm 2ln\left(1+\frac{b^{\prime }}{a^{\prime }}\left|\frac{\gamma -{\gamma }_c}{2}\right|\right)/b^{\prime }-\frac{\pm {\tau }_m-{\tau }_c}{\pm {\gamma }_m-{\gamma }_c}\hfill \\ ·\left(\gamma -{\gamma }_c\right)]+\frac{\pm {\tau }_m-{\tau }_c}{\pm {\gamma }_m-{\gamma }_c}·\left(\gamma -{\gamma }_c\right)+{\tau }_c \left|\gamma \right|\le {\gamma }_m\hfill \\ \frac{\pm \mathit{ln}\left(1\pm b\gamma /a\right)}{b} \left|\gamma \right|\ge {\gamma }_m\hfill \end{array}\hfill \\ K\left({\gamma }_0\right)=\frac{\pi b{\gamma }_0^2}{2\left(a_0+b_0{\gamma }_0\right)\left[\left(2a+b{\gamma }_0\right)-2b{\gamma }_0/\mathit{ln}\left(1+b{\gamma }_0/a\right)\right]}\hfill \\ G_{\mathit{max},N}=\frac{1}{a}=G_{\mathit{max},0}\left(1-U_N\right)=G_{max,0}\left(1-\frac{u_{re,N}+u_{ir,N}}{{\overline{\sigma }}_0}\right)\hfill \\ {\tau }_{ult,N}=\frac{1}{b}={\tau }_{ult,0}\left(1-U_N\right)={\tau }_{ult,0}\left(1-\frac{u_{re,N}+u_{ir,N}}{{\overline{\sigma }}_0}\right)\hfill \end{array}\right.$$

(1)

$$\{\begin{array}{c}\begin{array}{cc}{u}_{ir,0}\hfill & =0\hfill \\ u_{N }\hfill & =u_{re,N}+u_{ir,N}\hfill \\ & =\Delta u_{ir,N}\left(\frac{\tau }{{\tau }_{N,max}}-1\right)+u_{ir,N-1}+\Delta u_{ir,N}\hfill \\ & =\Delta u_{ir,N}\frac{\tau }{{\tau }_{N,max}}+u_{ir,N-1},N=1,2,3,\dots \hfill \end{array}\hfill \\ \Delta u_{ir,N}=\frac{c_{1,0}}{\sqrt{N_{ep}}}{\left(\frac{{\tau }_N}{{\stackrel{-}{\sigma }}_{N-1}}\right)}^{A_{4,0}}·\left[1-c_{1,a}{\left(K_c-1\right)}^{c_{1,b}}\right]\left({\overline{\sigma }}_0-u_{ir,N-1}\right)\hfill \\ N_{ep} ={\displaystyle {\displaystyle \sum }_{i=1}^N}{\left[\frac{{\tau }_i}{{\tau }_N}\right]}^{\alpha }\hfill \end{array}$$

(2)

where $u_{re,N}$ is the reversible pore-water pressure in the Nth stress cycle; $\Delta u_{ir,N}$ is the irreversible pore-water pressure increment in the Nth stress cycle; $u_{ir,N-1}$ is the irreversible pore-water pressure in the N−1st stress cycle; ${\overline{\sigma }}_0$ is the initial effective confining pressure; $c_{1,0}, c_{1,a}, c_{1,b}, A_{4,0}$ are the material constants, which are related to the density of the soil; $K_c$ is the consolidation ratio; $N_{ep}$ is the number of equivalent actions; ${\tau }_i$ is the shear stress amplitude in the Ith stress cycle $\left(1\le i\le N\right)$;$\alpha$ is the material constant; $G_{max,N}$,${\tau }_{ult,N}$ are the maximum shear modulus and ultimate shear stress after the N−1th stress cycle; $\left|\gamma \right|$ is the absolute value of the shear strain; and ${\tau }_c$ and ${\gamma }_c$ are the values of the shear stress and the shear strain, respectively, at the last reversal. ${\tau }_m,{\gamma }_m$ are the values of the shear stress and the shear strain at the biggest reversal, which are positive; $K\left({\gamma }_0\right)$ is the damping degradation coefficient; $a,b$ are related to the maximum shear modulus $G_{max}=1/a$ and the ultimate shear stress${\tau }_{ult}=1/b$ can be obtained by fitting the $G/G_0-\gamma$ experimental data; $a^{\prime }, b^{\prime }$ are the parameters derived from $a,b$ and ${\tau }_m, {\gamma }_m$; and $a_0, b_0$ can be obtained by fitting the $\lambda -\gamma$ experimental data. See the literature for the values of all the parameters [37,44].

After liquefaction, the shear modulus of the saturated sand decreases sharply, and the shear stress cannot reach the applied stress amplitude, resulting in a large shear deformation. In this summary, the regularities of the large deformation of saturated sand after the liquefaction are summarized. Their basic function expressions are determined through decomposition analysis of the large deformation stress–strain curves obtained from the liquefaction tests.

2.3.1. The Basic Function of Stress–Strain Relationship after Liquefaction

The soil constitutive model determination advance determines the skeleton curve and the basic functional relationship. Figure 5 shows the large deformation stress–strain of saturated sand under different cycles (N) of action measured by the test. As shown in Figure 6, the stress–strain curve (A-B-C-D-A) of each cycle can be decomposed into three processes: forward loading, unloading, and reloading. It can be disassembled into four stages: forward loading AB, first unloading stage BC, second unloading stage CD, and then loading DA. The following is the determination method of the basic functional relationship of the four stages.

• Initial loading curve AB

First, to determine the basic function of the skeleton curve with a large deformation, the relation of the AB segment should be determined.

It is not difficult to see from Figure 7a that section AB after the liquefaction is symmetric with the curve before the liquefaction (τ = γ/(a + bγ)) about the line τ = γ. Then, γ and τ can be interchanged in the stress–strain relation τ = γ/(a + bγ) of the hyperbolic model before the liquefaction, and the basic function of the liquefaction is shown in Equation (3).

$$\gamma =\frac{\tau }{a+b\tau }$$

(3)

It can be written as:

$$\tau =\frac{\gamma }{1/a-\gamma b/a}=\frac{\gamma }{a_1-b_1\gamma }$$

(4)

An inverse hyperbolic function (IHF) can express the initial skeleton curve. The parameter $a_1=1/a$ is the initial tangential modulus, as shown in Figure 7c. $a_1$ gradually decreases with the increase in the cycles. The difference between the skeleton curves of large deformation and small deformation also lies in the asymptotes. Before liquefaction, there is a transverse asymptotic line $\tau =1/b$, which means the ultimate shear strength of the soil. According to Equation (4), after liquefaction, there is a longitudinal asymptote $\gamma =a_1/b_1$, whose physical meaning is the ultimate shear strain of saturated sand after liquefaction. With the increase in the number of cycles and the complexity of the soil deformation, the initial skeleton curves of each cycle generally cannot reach the origin. The general function of the loading curve is:

$$\tau =\frac{\gamma -{\gamma }_d}{a-b\left(\gamma -{\gamma }_d\right)}-{\tau }_d$$

(5)

where the coordinate of A is $\left({\gamma }_d,{\tau }_d\right)$. As shown in Figure 7c, ${\gamma }_d$ is approximately zero, and ${\tau }_d$ decreases linearly with the accumulation of the cycles. Figure 7b has an excellent fitting result through Equation (5), which can be used as the initial skeleton curve of the large deformation.

• Unloading curve BCD

As shown in Figure 7d, the unloading curve BCD is divided into two sections: in unloading BC, the pore-water pressure decreases, the soil stress decreases sharply, and the stress–strain relationship presents a nonlinear state; when in unloading CD, the pore-water pressure continues to fall, and the stress–strain relationship decreases linearly. Since the shape of BC is familiar to AB, it can also be expressed by the IHF of the AB above:

$$\tau =\frac{\gamma -{\gamma }_{d1}}{a_1-b_2\left(\gamma -{\gamma }_{d1}\right)}-{\tau }_{d1}$$

(6)

where the coordinate of B is $\left({\gamma }_c, {\tau }_c\right)$ and the coordinate of C is $\left({\gamma }_{d1},{\tau }_{d1}\right)$. Point B has a large slope, while point C is the turning point of the two unloading stages and has a small slope.

Because CD is tangent to BC at C, CD can be expressed as follows:

$$\tau =(\gamma -{\gamma }_{d1})/a_1-{\tau }_{d1}$$

(7)

• Reloading curve DA

As shown in Figure 7e, the stress at the initial loading point D $\left({\gamma }_c^{\prime },{\tau }_c^{\prime }\right)$ gradually decreases with the accumulation of the number of cycles, and the strain value increases substantially. The stress–strain curve of the DA section can be expressed as:

$$\tau =\frac{\gamma -{\gamma }_d}{a_2-b_3\left(\gamma -{\gamma }_d\right)}-{\tau }_d$$

(8)

where $\left({\gamma }_d,{\tau }_d\right)$ is the stress–strain of the initial point A of the next cycle. As the number of cycles increases, A’s strain can be regarded as a constant value and the stress decreases linearly. When the curve reaches the initial point of the next cycle, the reloading stage ends, the stress cycle ends, and the next cycle begins.

To sum up, BC and AB are similar in shape, DA and AB are symmetric about zero, CD is a straight line, and all the curves are connected by the control points (A, B, C, D). The two loading stages are similar in shape and can be expressed in a functional form. The unloading stage presents different shapes, divided into two unloading stages: (1) the unloading stress–strain curve has the maximum tangent modulus in the initial unloading curve, and the tangent modulus of the curve decreases gradually with the decrease in the strain. (2) In the second stage of the unloading curve, the stress–strain curve decreases linearly.

2.3.2. Loading and Unloading Criteria

The liquefaction of pore-water pressure in the stress cycle has the corresponding fluctuation. The effective stress in the zero value fluctuates up and down. In this paper, when the effective stress reaches zero, the irreversible pore-water pressure does not rise, and the reversible pore-water pressure fluctuates up and down with the change in the cyclic load. The relationship between the reversible pore-water pressure and cyclic stress is similar to before liquefaction. That is, the modulus and strength of the soil after liquefaction are not modified by the pore-water pressure. Based on the above analysis, the following loading and unloading criteria are determined:

The stress–strain relationship satisfies IHF in the first cyclic forward loading process AB after the liquefaction. The curve of AB and the soil before the liquefaction is tangent at point A.

(1)

The BC segment satisfies IHF. The slope of point C is equal to the slope of point A. The stress–strain value of point C is related to point B. The relation can be obtained by fitting the test data.

(2)

CD section stress–strain relationship is linear. The line is tangent to curve BC at point C.

(3)

DA satisfies IHF and goes through the initial point A of the next cycle. This curve is tangent to the initial loading curve of the next cycle.

(4)

At the beginning of the new cycle, the tangent modulus at point A is related to soil properties and the number of cycles. The curve is symmetric with the DA curve at about point A.

In the above cyclic loading and unloading process, the tangent modulus $a_1$ of point A and point C should be obtained by fitting the test data after liquefaction, and the other parameters ($b_1,b_2,b_3$) can be determined by the inflection point in the stress–strain history.

The large post-liquefaction deformation stress–strain relationship of saturated sand is related to the soil properties and the number of stress cycles after liquefaction. The lower the stress amplitude of the soil, the lower the reduction rate of the tangent modulus of the stress–strain curve is, and the more intensive the peak points B and D of each cycle are.

2.4. The Constitutive Equation

2.4.1. Calibration Parameter

To get the relationship between $a_1$ and the cyclic action times N, according to the results of the cyclic triaxial test of the saturated coral with a relative density of 45%, the four groups of CSRs were 0.325, 0.3, 0.25, and 0.22. As shown in Figure 8, the variation in $a_1$ with N after the liquefaction is drawn. The four points in each group of cycles, respectively, represent the test results of four groups of different CSRs. The eight points in the same color are the fitting points of the eight cycles after the sand liquefaction under the action of the same CSR. The abscissa is the number of cycles, and the ordinate is the ratio $a/a_{1,n}$of the soil pre-liquefaction parameter $a$ to the tangent modulus of each cycle $a_{1,n}$.

The above figure shows that the tangential modulus $a_{1,n}$ of the soil after liquefaction has a consistent relationship with the number of cycles, which does not change with the change in the CSR and is a linear reduction relationship with a reduction rate of about 0.21. The modulus tends to zero at each CSR at the eighth cycle after the liquefaction.

Figure 9a shows the strain ratio at point C and point B in each cycle under different CSRs. Four points in each cycle number represent the test results of the four groups of different CSRs. The points in the same color are the fitting points of eight cycles after liquefaction. The abscissa is the number of cycles, and the ordinate $\frac{{\gamma }_{d1}}{{\gamma }_c}$ is the strain ratio at point C and point B. Figure 9b shows the stress ratio ${\tau }_{d1}/{\tau }_c$ at point C and point B in different CSRs.

The stress–strain ratio at points C and B of each cycle can be approximately regarded as a constant value, and the relationship does not change with the change in the CSR. The tangent modulus of the curve at point C is consistent with the initial tangent modulus of the cycle.

At this point, we have completed the calibration of the parameters of the large deformation constitutive model of saturated sand. According to the results, the concrete calculation formula of the large deformation constitutive model can be determined by combining the loading and unloading criterion and the primary function of the stress–strain relationship.

2.4.2. Constitutive Relation

After saturated sand is subjected to cyclic loading, the pore-water pressure increment can be obtained by introducing the model. Then, the pore-water pressure increment can be accumulated to obtain the effective stress. The stress–strain relationship at this stage is shown in Equation (1). When the effective stress is zero, the soil layer is considered liquefied, and the stress–strain relationship after liquefaction follows the five large deformation loading and unloading criteria. The stress–strain constitutive relationship of saturated sand after liquefaction is obtained by combining the basic function formula of the stress–strain relationship, as shown in Equation (9):

$$\begin{array}{c}\begin{array}{c}\tau =\frac{\gamma -{\gamma }_d}{a-b\left(\gamma -{\gamma }_d\right)}-{\tau }_d,n=1,\Delta \gamma \ge 0,\gamma \ge {\gamma }_d\\ \tau =\frac{\gamma -{\gamma }_{d1}}{a_1-b_2\left(\gamma -{\gamma }_{d1}\right)}-{\tau }_{d1},\Delta \gamma <0,\gamma \ge {\gamma }_{d1}\\ \tau =\frac{\left(\gamma -{\gamma }_{d1}\right)}{a_1}-{\tau }_{d1}, \Delta \gamma <0,\gamma <{\gamma }_{d1}\\ \tau =\frac{\gamma -{\gamma }_d}{a_1-b_3\left(\gamma -{\gamma }_d\right)}-{\tau }_d,\Delta \gamma \ge 0,\gamma <{\gamma }_d\\ \tau =\frac{\gamma -{\gamma }_d}{a_1-b_1\left(\gamma -{\gamma }_d\right)}-{\tau }_d,n>1,\Delta \gamma \ge 0,\gamma \ge {\gamma }_d\\ \end{array}\end{array}\}$$

(9)

where:

$$\begin{array}{c}{b}_2=\frac{a_1}{{\gamma }_c-{\gamma }_{d1}}-\frac{1}{{\tau }_c-{\tau }_{d1}}\\ b_3=\frac{a_1}{{\gamma }_d-{\gamma }_c^{\prime }}+\frac{1}{\left({\tau }_c^{\prime }-{\tau }_d\right)}\\ b_1=\frac{a_1}{{\gamma }_d-{\gamma }_c^{\prime }}-\frac{1}{2{\tau }_d-{\tau }_c^{\prime } }\end{array}\}$$

(10)

where n is the number of cycles after liquefaction, the tangential modulus $a_{1,n}$ has a consistent relationship with n, and there is a linear reduction relationship with a reduction rate of about 0.21. $a_{1,n}$ tends to zero at each CSR at the eighth cycle after liquefaction. $A\left({\gamma }_d,{\tau }_d\right), B\left({\gamma }_c, {\tau }_c\right), C\left({\gamma }_{d1}, {\tau }_{d1}\right),D\left({\gamma }_c^{\prime }, {\tau }_c^{\prime }\right)$ and the stress–strain ratio at C and B of each cycle can be approximately regarded as a constant value, and the relationship does not change with the change in the CSR. The tangent modulus of the curve at point C is consistent with$a_{1,n}$. The parameters $(b_1,b_2,b_3)$ can be determined by the inflection point in the stress–strain history.

Equations (1) and (9) constitute the saturated sand’s nonlinear dynamic constitutive equations before and after the liquefaction. Their constitutive parameters are few, which can be obtained from undrained cyclic shear test data.

3. Calculation Method

3.1. Time-Domain Seismic Response Analysis

Computational methods to analyze the liquefaction problems solve the equations of motion. The stresses and displacements of the solid particle framework satisfy the equations of the motion and conservation of the mass. The two classic approaches used to discretize a domain in space for a numerical solution are finite differences and finite elements [45,46,47]. The time-domain finite differences method is used to analyze site nonlinear seismic response [48], but its use is limited due to the soil’s lack of post-liquefaction stress–strain relationship. Given this, we develop a program, Soilresp1D, and the large deformation constitutive is embedded into the program. It can calculate the 1D liquefiable site seismic response and analyze the dynamic reaction of saturated sand after liquefaction. In Soilresp1D, finite element and difference methods are applied, and the boundary condition is a multi-transmitting formula [49,50]. The initial conditions are:

$$\begin{array}{c}\begin{array}{c}\tau (0, z)=0\\ \gamma (0,z)=0\\ \upsilon (0,z)=0\end{array}\\ \tau (\mathrm{t},0)=0\end{array}\}$$

(11)

In the formula, $\upsilon (t,z)$ is the velocity; $\tau (t,z)$ and $\gamma (t,z)$ are the shear stress and shear strain; and $\rho \left(z\right)$ is the density. This paper uses the staggered grid method to solve the initial-boundary value problem (Liao, 2002). Based on the free boundary condition Equation (11), the central difference discrete form of the motion balance equation is as follows:

$$\begin{array}{c}{\upsilon }_1^{p+1}={\upsilon }_1^p+\frac{\mathsf{\Delta }t}{m_1}{\tau }_1^p\\ {\upsilon }_n^{p+1}={\upsilon }_n^p+\frac{\mathsf{\Delta }t}{m_n}\left({\tau }_n^p-{\tau }_{n-1}^p\right),n=2,3,\dots ,N\\ m_n=\frac{1}{2}\left({\rho }_n{h}_n+{\rho }_{n-1}{h}_{n-1}\right),n=2,3,\dots ,N\\ m_1=\frac{1}{2}{\rho }_1{h}_1\end{array}\}$$

(12)

The central difference discrete form of the constraint conditions is as follows:

$${\gamma }_n^{p+1}=\frac{\Delta t}{h_n}\left({\upsilon }_{n+1}^{p+1}-{\upsilon }_n^{p+1}\right)+{\gamma }_n^p,n=1,\dots ,N$$

(13)

The relationship between the shear stress and shear strain proposed in this paper can be expressed as:

$${\tau }_n^{p+1}=\tau \left({\gamma }_n^{p+1}\right)$$

(14)

Equations (12)~(14) give an explicit recursive calculation method for the node velocity of each layer and the stress–strain relationship between the layers.

The recursive formula of a boundary node is as follows:

$${\upsilon }_{N+1}^{p+1}={\upsilon }_N^p-{\upsilon }_{I,N}^p+{\upsilon }_{I,N+1}^{p+1}$$

(15)

where ${\upsilon }_{N+1}^{p+1}$ is the velocity of the artificial boundary node at time p + 1, ${\upsilon }_N^p$ is the velocity of the node N at time p, ${\upsilon }_{I,N}^p$ is the velocity of the incident wave at time p at the node N, and ${\upsilon }_{I,N+1}^{p+1}$ is the velocity of the incident wave at time p + 1 at the artificial boundary.

The above time-domain analysis method is a second-order calculation precision, which is conditionally stable. The stability condition is:

$$\Delta \mathrm{t}\le \min \left(\frac{h_n}{c_n}\right),n=1,\dots ,N$$

(16)

where $c_n$is the shear wave velocity of the medium in the NTH layer, and $h_n$ is the thickness of the layer, determined by the following equation:

$$h_n\le \left(\frac{1}{6}~\frac{1}{10}\right)T_{min}{c}_n$$

(17)

where $T_{min}$ is the shortest period of input fluctuation with an engineering signification.

According to the given input ground motion, the seismic response of the soil layer at any depth can be calculated by Soilresp1D.

3.2. Algorithm Implementation Process

In the case of the known soil parameters, combined with the above large deformation constitutive relation of saturated sand, the time-domain nonlinear dynamic response analysis method can be used to gradually calculate the saturated sand layer’s dynamic response at any time at any section. The detailed operation process of the program can be summarized as follows:

• Read site information. For non-liquefiable clay and other soils, the constitutive parameters were $a,b$. So, the saturated sand and the liquefaction of the former constitutive parameters$a,b$, are the liquefied constitutive parameters after $a/a_{1,n}$, ${\gamma }_{d1}/{\gamma }_c$, and ${\tau }_{d1}/{\tau }_c$.
• According to the site model, the stratification of the site, the calculated thickness of each layer, and the calculated time step are determined. Record the number of layers and time nodes.
• According to the staggered grid method, the interlayer node velocity and strain are calculated by displaying the recursive formula.
• According to the stress–strain change trend, judge whether it is the stress–strain inflection point or peak point. If so, record the stress–strain value and calculate the state parameters required in the transportation process.
• They are determining the category of the soil layer. If it is non-liquefiable soil, according to Equation (1), the stress is calculated according to Equation (1). If the saturated sand is not liquefied, the pore-water pressure is calculated according to Equation (2), and the relevant parameters are modified. The stress is calculated according to the effective stress constitutive Equation (1). In liquefied saturated sand, the stress is calculated through the large deformation constitutive Equation (9).

Determine whether this calculation step is the last time node of the topsoil layer. If not, continue to calculate the amount of motion state of the upper soil layer at the next moment. If it is the last time node of the topsoil layer, the calculation ends.

4. Verification

To illustrate our method’s reasonability, Soilresp1D is used to analyze the nonlinear seismic response of the experimental model site. By comparing the numerical simulation results and experimental results, the verification of ours constitutive in this paper is demonstrated.

The overlayer of the site consists of clay, saturated sand, and underlying bedrock 32 m thick; the site parameters are given in

Table 1. Site calculation model.

Soil	Thickness (m)	${\mathbf{V}}_{\mathbf{s}} (\mathbf{m}/\mathbf{s})$	Density (t/m3)
Clay	6	120–142	1.95
Pine fine sand	3	142–153	1.49
Clay	21	153–231	1.95
Bedrock	2	511	2.65

, and the nonlinear characteristic parameters are listed in

Table 2. Nonlinear characteristic parameters of the site.

Soil	${\mathit{a}}_1\mathit{a}\mathit{n}$	${\mathit{\gamma }}_{\mathit{d}}/{\mathit{\gamma }}_{\mathit{c}}$	${\mathit{\tau }}_{\mathit{d}}/{\mathit{\tau }}_{\mathit{c}}$	$\mathit{b}/\mathit{a}$	${\mathit{a}}^1\left({10}^{-3}\right)$	${\mathit{b}}^1$	${\mathit{c}}_{1,0}$	${\mathit{c}}_{1,\mathit{a}}$	${\mathit{c}}_{1,\mathit{b}}$	${\mathit{A}}_{4,0}$	${\mathit{C}}_3$	${\mathit{B}}_3$
pine fine sand	0.34	0.75	0.11	1934	2.24	5.20	33.71	0.38	0.56	2.61	0.49	1.0
clay				1160	0.82	5.90

. In the numerical simulation, the distance $\Delta z$ is 1 m, and the time step is 0.0025 s. 7–9 m from the top to the bottom is the saturated pine fine sand. The input sine wavelength is 48 s and the period is 1 s, with a peak acceleration of 70 cm/s². The undrained cyclic triaxial test of the saturated Nanjingpine fine sand is carried out to compare with the test results. In the test, the effective confining pressure is 100 kPa, the CSR is 0.13, and the loading frequency is 1 Hz.

The stress of the sand at 7 m is consistent with that of the test sand. Now, the simulation results of this soil layer are compared with the test results, as shown in Figure 10. The calculated pore-water pressure is consistent with the measured, which can reflect how the irreversible pore-water pressure gradually increases with the accumulation of cycles and reflect the fluctuation of the reversible pore-water pressure under the action of each cycle load. The sand layer at 7 m began to liquefy at 16.5 s, while the sand at 8 m and 9 m did not liquefy, and the pore-water pressure rose slowly after the liquefaction at 7 m. Due to the small confining pressure of the upper saturated sand layer, it is also liquefied. After the liquefaction of the upper layer, the pore-water pressure of the lower saturated sand layer tends to be gentle, which is consistent with the actual observation phenomenon.

Figure 10b: the time-series peak value of ground acceleration is 124 cm/s², and the amplification coefficient is 1.77. The ground acceleration attenuates significantly after the liquefaction. Figure 10c: the stress–strain results before the liquefaction show that the numerical simulation results are similar to the test results. The shear modulus decreases with the increase in the number of cycles. Figure 10d shows the large-deformation stress–strain hysteresis loop after the liquefaction, and its variation trend is similar to the test results.

The numerical simulation results show that the amplitude of the ground acceleration reaction (PGA) before the liquefaction is significantly larger than that after the liquefaction. The liquefaction of saturated sand has a noticeable damping effect on the dynamic response of the site, which is consistent with some in situ observations.

5. Conclusions

This paper studied the stress–strain variation, loading, and unloading criterion, and the primary function formula of the saturated sand after the liquefaction, according to the dynamic triaxial test of saturated sand, and proposed a systematic post-liquefaction deformation constitutive model. This paper also set up a time-domain analysis method for simulating the dynamic response of saturated sand before and after liquefaction, combined with the staggered grid method, the multi-transmitting boundary condition, and the effective stress constitutive model before the liquefaction. A 1D time-domain site dynamic response analysis program, Soilresp1D, was developed using the method. According to the experimental results, the main findings are as follows:

(1)

The stress–strain curve after the liquefaction can be divided into four loading and unloading stages. IHF can describe the two loading and the first unloading stages. The second unloading stage is linear.

(2)

The tangent modulus of the unloading curve connection point is equal to that of the initial loading point. The two unloading curves are tangent at the connection point. The reloading curve is symmetric with the initial loading curve of the next cycle about the initial loading point.

(3)

The tangent modulus of the initial loading point decreases linearly with the increase in the number of cycles. The test results show that the initial loading curve of the soil is close to the level after eight times of cyclic loading.

The numerical simulation results show that the liquefaction has a damping effect on the dynamic response of the site, which verifies the model’s reliability. This model can be used to analyze the post-liquefaction deformation seismic response in liquefiable sites. It is worth noting that the post-liquefaction deformation constitutive has been verified by the dynamic response under symmetric cyclic loads, which can provide a basis for the study of large deformation constitutive suitable for asymmetric ground motion loads.