Abstract
This study employed a modified elastoplastic constitutive model that can systematically describe the monotonic and cyclic mechanical behaviors of typical marine soils combining the subloading, normal, and superloading yield surfaces, in the seismic response analysis of three-dimensional (3D) marine site. New evolution equations for stress-induced anisotropy development and the change in the overconsolidation of soils were proposed. This model can describe the unified behaviour of unsaturated soil and saturated soil using independent state variables and can uniquely describe the multiple mechanical properties of soils under general stress states, without changing the parameter values using the transform stress method. An effective stress-based, fully coupled, explicit finite element–finite difference method was established based on this model and three-phase field theory. A finite deformation analysis was presented by introducing the Green-Naghdi rate tensor. The simulation and analysis indicated that the proposed method was sufficient for simulating the seismic disaster process of 3D marine sites. The results suggested that the ground motion intensity would increase due to the local uneven complex topography and site effect and also provided the temporal and spatial distribution of landslide and collapse at the specific location of the marine site.
1. Introduction
The seismic safety of marine sites in strong earthquake is receiving increased attention because of the potential serious hazard of earthquakes to offshore engineering. Three-dimensional (3D) elastoplastic seismic response analysis of marine sites is a focus of and frontier topic in geotechnical earthquake engineering. The investigation and analysis of the typical earthquake damage phenomena and failure characteristics of 3D marine sites in large earthquake are of great significance.
In the seismic response analysis of 3D marine sites, one of the key points is to build a constitutive model that can suitably describe the mechanical properties of marine soils; in particular, a reasonable description of the dynamic characteristics of marine soils is the most important task. Considering the coexistence of unsaturated and saturated soils on an island, developing a unified constitutive model that can describe unsaturated and saturated soil continuously is difficult.
The mechanical behaviours of soil are dependent on the state conditions of the soil (overconsolidation (influence of density), structure (whether the soil is disturbed or not, Asaoka et al. [1]), anisotropy (mainly stress induced), and degree of saturation) as well as the loading conditions (drained or undrained condition, monotonic or cyclic loading, and other loading paths aside from the conventional triaxial condition). Hashiguchi and Ueno [2] proposed the concept of “subloading,” which made the description of soil overconsolidation possible. The concept of “superloading,” proposed by Asaoka et al. [1], together with the concept of subloading, not only allows the description of overconsolidation but also explains the effect of the soil structure. Zhang et al. [3] introduced a new approach to describe the stress-induced anisotropy. They noted that the change in density is influenced by not only plastic stretching and elastic unloading but also stress-induced anisotropy. Based on their model, the mechanical behaviour of soil subjected to different loadings under different drainage conditions was simulated. The above constitutive models are all based on the Cam-Clay model (Roscoe et al. [4]) or the modified Cam-Clay model (Schofield and Wroth [5]). Therefore, the physical meaning of the model is easy to understand, and relatively few parameters need to be employed. Based on these models, this study establishes a more widely applicable constitutive model.
Considerable research has been performed experimentally, empirically, and theoretically on the mechanical behaviours of unsaturated soils, which are widely found on islands. Several constitutive models for unsaturated soils have been developed since the work of Alonso et al. [6]. These researchers proposed the Barcelona basic model (BBM) using LC conception, and this model is typically regarded as one of the basic models for unsaturated soils. Constitutive models based on the framework of the BBM have been presented by Kohgo et al. [7, 8], Cui and Delage [9], and Sun et al. [10], among others. In these models, the stress-suction-strain relations of unsaturated soil were considered explicitly, whereas the degree of saturation was not considered directly. Constitutive models considering the influence of the degree of saturation were also proposed in the same period by Kato et al. [11], Gallipoli et al. [12], Muraleetharan et al. [13], and Sheng et al. [14].
The majority of the models in the literature use stress and suction as independent state variables. Ohno et al. [15] made an assumption that the void ratio-logarithmic effective stress relation e-lnp is dependent on the degree of saturation based on various experimental results and modelled the dependency with a simple relation. In the model proposed by Ohno et al. [15], only the effective stress, which was defined by a combination of net stress, suction, and “effective degree of saturation,” was used as a state variable, and a Cam-Clay-type model for unsaturated soil was proposed based on the effective stress. The advantage of using effective stress or skeleton stress instead of net stress and suction as independent state variables in the modelling of unsaturated soils is that it is considerably easier and smoother to describe the behaviours of soil from an unsaturated state to a saturated state and vice versa if the moisture characteristics of the soil are properly described and incorporated into the constitutive model. Based on this type of model, a numerical analysis related to boundary value problems of saturated and unsaturated grounds can be performed in a unique manner. Pioneering research on this topic was conducted by Kanazawa et al. [16].
Zhang and Ikariya [17] proposed a constitutive model for unsaturated soil that used skeleton stress and degree of saturation as independent state variables and can accurately describe the influence of the degree of saturation. The constitutive model described the behaviour of unsaturated soil as well as saturated soil, as the skeleton stress can smoothly shift to effective stress if the degree of saturation changes from an unsaturated condition to a saturated condition. Overconsolidation, one of the main features of soils that were discussed in the models for saturated soils (Zhang et al. [3]), was also considered together with the degree of saturation. Research on unsaturated soil considering the overconsolidation effect for saturated soils includes the work by Kohgo et al. [8] which employs the subloading concept and the work by Russell and Khalili [18] regarding the concept of the bounding surface.
Other mechanical features, such as the soil structure formed during the sedimentary process of the soil (Asaoka et al. [1]) and stress-induced anisotropy (Zhang et al. [3]), are often discussed for saturated soil but are rarely incorporated into the existing models. These factors can be incorporated into the proposed model within the framework of Zhang and Ikariya [17]. The primary objectives of this study are to (1) modify the constitutive model considering the structure and stress-induced anisotropy of marine soils and propose new evolution equations for the development of stress-induced anisotropy and the change in overconsolidation; (2) make the above model uniquely describe the multiple mechanical properties of soils under general 3D stress states without changing the parameter values; (3) introduce the improved model into soil-water-air three-phase coupled field theory, conduct an effective stress-based, fully coupled, explicit finite element–finite difference method, (4) perform an elastoplastic seismic response analysis of a 3D typical marine site, and discuss the local uneven complex topography and site effect (Borcherdt [19] and Anderson et al. [20]). These results will help bridge information gaps regarding the seismic safety of marine sites.
2. Unified Description of Soil Behaviour
2.1. State Line
Based on the experimental results, a quantitative relation for the void ratio-logarithmic mean skeleton stress e-lnp is established using the degree of saturation as a state variable. It is assumed that a normally consolidated line in the unsaturated state (N.C.L.S.) is parallel to the normally consolidated line in the saturated state (N.C.L.) but in a higher position than the N.C.L. [17], which means that unsaturated soil can maintain higher void ratio than saturated soil under the same mean skeleton stress. It is also assumed that the critical state line in the unsaturated state (C.S.L.S.) is parallel to critical state line in the saturated state (C.S.L.) (Roscoe et al. [4]) but in a higher position than the C.S.L.
The N.C.L.S. and C.S.L.S. are expressed as follows:where and are the void ratios of the N.C.L.S. and C.S.L.S., with a reference mean skeleton stress and a certain degree of saturation . is the compression index and has the same value for the saturated and unsaturated states. and are the mean skeleton stress and second invariant of the deviatory skeleton stress tensor, respectively. is the stress ratio, and is the stress ratio at the critical state and has the same value for the saturated and unsaturated states (Cui and Delage [9]).
2.2. Yield Surface
Based on the subloading (Hashiguchi and Ueno [2]), normal, and superloading (Asaoka et al. [1]) yield surfaces in the plane (shown in Figure 1) and the cyclic mobility model (Zhang et al. [3], Zhang et al. [21], and Ye et al. [22]), a new yielding function can be written as shown in (2) (Li [23]). Using the transform stress method (Yao et al. [24]), this model can uniquely describe the multiple mechanical properties of soils under general stress states, without changing the parameter values. The model is suitable for 3D static-dynamic analysis for complex 3D marine sites.where , is the reference mean skeleton stress; is the state variable of the degree of saturation (Zhang and Ikariya [17]); is the void ratio at the reference mean skeleton stress; is the plastic volume strain. The similarity ratio of the normal yield surface to the superloading surface and the similarity ratio of the subloading surface to the superloading surface in the transform stress space are given aswhere (), (), and () represent the present stress state, the corresponding normally consolidated stress state, and the structured stress state in the transform stress plane, respectively. The variables involved in (2) and (3) are defined aswhere is the deviatory stress tensor, is the anisotropic stress tensor, and is the swelling index.
2.3. Consistency Equation for the Yield Surface
Considering the advantages of the associated flow rule (Hashiguchi and Chen [25] and Zhang et al. [3]), the consistency equation for the new yield surface can be given as
In (5), evolution equations must be provided for the development of the state variables, namely, the anisotropic stress tensor, degree of structure, degree of overconsolidation, and degree of saturation.
2.3.1. Evolution Rule for the Anisotropic Stress Tensor
Different from the work by Hashiguchi and Chen [25] and Asaoka et al. [26], the following evolution rule is used for the anisotropic stress tensor here:where is the evolution parameter of anisotropy and is a limitation parameter, the value of which varies from 0.90 to 0.99 according to the definition of “completely liquefaction” in the engineering sense. A value of that is closer to 1.0 indicates that the mean effective stress will be closer to zero in the cyclic mobility region. Without losing physical meaning, we fixed the value of at 0.95 throughout this paper. In (6), the artificial limitation on the development of anisotropy originally proposed by Hashiguchi and Chen [25] is no longer necessary. The increment in anisotropy, which plays a critical role in the evolution rule for the overconsolidation, will be discussed in detail later.
The increment of the evolution rule for the anisotropic stress tensor can be written aswhere provides a natural physical limitation on the development of anisotropy automatically; according to (6), the development of anisotropy will stop at the state when .
2.3.2. Evolution Rule for the Degree of the Structure
The following evolution rule for the degree of the structure , which is the same as that proposed by Asaoka et al. [1, 26, 27], is adopted:where is a parameter that controls the rate of the collapse of the structure during shearing.
2.3.3. Evolution Rule for the Degree of Overconsolidation
In the present model, the changing rate of overconsolidation is assumed to be controlled by two factors, namely, the plastic component of stretching, which was employed as the only factor in the work by Asaoka et al. [26], and the increment in anisotropy. For degree of overconsolidation , the new definition is as follows:where is the degradation parameter of overconsolidation and and are two fitting parameters. In general, and . is proportional to the norm of the plastic component of stretching .
The term is used because, in the region close to the zero point, the effective mean stress is extremely small and the change in the effective stress is slow compared to the rapid development of plastic strain. Therefore, the changing rate of overconsolidation must be reduced. When the effective mean stress is extremely large, that is, , the change of the effective stress is also slow.
2.3.4. Evolution Rule for the Degree of Saturation
The mechanical behaviours of unsaturated soil are strongly influenced by the hydraulic characteristics, which are expressed as the relationship between the degree of saturation and suction. A proper model of the moisture characteristic curve (MCC) should be proposed in advance to describe the mechanical behaviour of unsaturated soil, particularly the transition from the unsaturated state to the saturated state. Zhang and Ikariya [17] proposed a new MCC model that can properly consider hydraulic hysteresis using the boundary surface concept. For convenience, we will employ this MCC model.
The evolution rule for degree of saturation iswhere and are the degrees of saturation of the residual and saturated conditions, respectively. Equation (10) indicates that changes linearly with the degree of saturation. is the void ratios at the N.C.L.S. under the reference mean skeleton stress when the degrees of saturation is in the residual state; that is, .
Due to the associated flow rulethe plastic volumetric strain increment and plastic shear strain increment can be expressed as
The following new relation can be obtained by substituting (7), (8), (9), (10), and (12) into (5):which results in the positive valuable that can be written aswhere
The strain increment can be divided into elastic and plastic parts as follows:
Using Hooke’s theory with stiffness tensor , the incremental stress tensor can be expressed as
Then, we obtain
We define the loading criteria as
The denominator of (18) is typically positive (Asaoka et al. [28] and Asaoka et al. [26]); therefore, is equivalent to the following relation:
Substituting (18) into (11) yields
Thus, we obtainwhere
Equation (24) will be used in the next soil-water-air coupled three-phase field theory.
The above improved model can describe the loose to dense soils subjected to monotonic loading under undrained and drained condition. It can also describe the loose to dense soils subjected to cyclic loading under undrained and drained conditions.
3. Modified Soil-Water-Air Coupled Three-Phase Field Theory
Based on Zienkiewicz’s theory, Oka et al. [29] proposed a coupled field theory:(I)Displacements of the solid and pore water pressure are chosen as unknown variables.(II)The finite element method (FEM) is used to discretize the motion equation for the mixture.(III)The finite difference method (FDM) is used to discretize the continuity equation for the pore fluid.
We modified the above theory and established the equations of three-phase field theory based on a finite deformation algorithm.
Based on mixture theory (Darcy [30], Truesdell and Toupin [31], Borja [32], and Fredlund [33]), the balance equation for soil-water-air coupled three-phase is expressed as (Li [23])where is the appearance density of the mixture, is the total stress tensor, and is body force. can be expressed asin which the appearance densities of the solid phase , water phase , and air phase are defined aswhere , and are the densities of the solid phase, water phase, and air phase, respectively; is the porosity; is the degree of saturation of soil.
Based on Bishop theory, total stress tensor can be written aswhere is the effective stress and , and are the mean pore pressure, pore water pressure, and pore air pressure, respectively.
The continuum equation for water is expressed as
The continuum equation for air is expressed as
In (29) and (30), , is displacement vector of solid phase; and are coefficients of permeability for water and air, respectively; is the volume elastic coefficient of the water phase; is the volume elastic coefficient of the air phase; and are the excessive pore water pressure and excessive pore air pressure, respectively (Li [23]).
After the spatial discretization, (25), (29), and (30) can be represented as (31), (32), and (33), respectively:where different shape functions are used for displacement and pore pressure (water and air) in the FEM discretization. The variables of displacement are given at the nodes, and the variables of pore pressure (water and air) are given at the gravitational centre and are denoted by and . is the nodal displacement vector of all nodes in one element. , and are the calculated parameters (Li [23]).
Based on the Newmark-β method, the new weak form formulations for the FEM are as follows:where , and are the calculated parameters (Li [23]).
Based on (34), the program is capable of solving repeated soil-water-air coupled static and dynamic boundary value problem using an infinitesimal deformation algorithm or finite deformation algorithm. The constitutive model and computational methods have been verified strictly (Zhang et al. [21], Zhang and Ikariya [17], Ye et al. [22], and Li [23]).
4. Application
4.1. Numerical Model
A 3D numerical model of a typical submarine site is shown in Figure 2. The horizontal dimension of the model is m and is divided equally into elements in the and directions. The sea level elevation is ±0.000 m, the elevation of the maximum height of the soil at the monitoring point (MP) 15 is −5.000 m, and the elevation of the model bottom is −59.226 m. The vertical dimension of the soil and bedrock (54.226 m) is divided into 16 elements. Therefore, the total number of elements is 160,000. The numbers in circles shown in Figure 2 are the monitoring points for acceleration, velocity, and displacement.
As shown in Figure 2, three monitoring surfaces (MS) and three monitoring lines (ML) are selected to present the results. MS1, including MPs 0, 1, 2, 3, and 4, is located at the bottom of the model; MS2, including MPs 5, 6, 7, 8, and 9, is located at the top of the bedrock; MS3, including MPs 10, 11, 12, 13, and 14, is located at the top of the model. ML1 is composed of MPs 0, 5, and 10 in a vertical line; ML2 is composed of MPs 15–19, 10, and 20–24; ML3, which is vertical to ML2, is composed of MPs 25–27, 10, and 28–30. MPs 5 and 0 are the projection points of point 10 on the top of the bedrock and the input surface. MPs 1, 6, and 11 are on the same vertical line, and the same situation applies to points (2, 7, and 12), (3, 8, and 13), and (4, 9, and 14).
4.2. Material Parameters
Experimental data obtained from static and cyclic triaxial compression tests were used to determine the parameters of the marine soil [1–5, 17]. The material parameters of marine soil involved in our model are listed in Tables 1, 2, and 3. The comparison between the simulated stress-strain-dilatancy relations of soil sample under monotonic/cyclic loading in conventional triaxial tests is shown in Figures 3 and 4, respectively. The bedrock is assumed to be elastic, and its parameters are shown in Table 4.
4.3. Input Ground Motion
The amplitude, spectral characteristics, and duration of the input ground motion have considerable influences on the seismic response of a marine site. The real three-component ground motions, shown in in Figure 5(a), are selected as the input ground motions in the analysis. The original records were filtered and baseline-corrected; the velocity time history before and that after adjustment are shown in Figure 5. As shown in Figure 5, the portion in the blue dashed rectangular box was the selected input time histories. Before inputting the model, the velocity time history should be preliminarily reduced by half and transformed into stress time history (Jing et al. [34]). EW is in the east-west direction (positive: east), NS is in the north-south direction (positive: north), and UD is in the up-down direction (positive: up).
(a) Original records of acceleration
(b) Acceleration after adjustment
(c) Velocity
(d) Displacement
4.4. Boundary Condition
For static analysis, the four vertical boundaries were simply supported, and the bottom boundary was fixed. However, in the dynamic analysis, the four vertical boundary conditions were free-field boundaries (shown in Figure 6), and the bottom boundary condition was a viscous damping boundary (Lysmer and Kuhlemeyer [35] and Jing et al. [34]).
4.5. Results of the Static Analysis
The effects of water pressure and the pore water pressure cannot be neglected in the seismic response analysis of marine site. As shown in Figure 7(a), the pore water pressure varies with the depth. The initial stress, which was generated by the self-weight stress, is shown in Figure 7(b).
(a) Pore water pressure (Pa), WLE: +5.000 m, IGMSR: 1
(b) -direction initial stress (Pa), WLE: +5.000 m, IGMSR: 1
4.6. Results of the Seismic Response Analysis
4.6.1. Spatial and Temporal Distribution of Displacement
Figure 8 shows the contours of displacement at the following dynamic times: 5, 10, 20, 30, 40, 50, 60, and 67 s. The deformation increases from the bottom to top of the model. The residual deformations of the MPs on the top are considerably larger than the corresponding residual deformations at the foot of the hill. The numerical simulation results correspond well with the field survey.
(a) Time: 5 s, WLE: +5.000 m, IGMSR: 1
(b) Time: 10 s, WLE: +5.000 m, IGMSR: 1
(c) Time: 20 s, WLE: +5.000 m, IGMSR: 1
(d) Time: 30 s, WLE: +5.000 m, IGMSR: 1
(e) Time: 40 s, WLE: +5.000 m, IGMSR: 1
(f) Time: 50 s, WLE: +5.000 m, IGMSR: 1
(g) Time: 60 s, WLE: +5.000 m, IGMSR: 1
(h) Time: 67 s, WLE: +5.000 m, IGMSR: 1
The spatial and temporal distribution of the displacements shown in Figure 8 illustrate that the slope has a large displacement. In the 3D simulation, we can identify the range of sliding, which is critical to our seismic design and reinforcement of marine slope sites.
4.6.2. Dynamic Time History
Typical acceleration, velocity, and displacement time histories at MPs 0, 5, 10, and 23 are shown in Figures 9–11.
The PGD, PGV, and PGA at MPs 0–30 are presented in Table 5.
Figure 12 shows the PGA of the MPs on MSs 1, 2, and 3. The PGA amplifies with increases in the altitude of the MPs.
Figure 13 shows the calculated PGA of the MPs on MLs 1, 2, and 3. The PGA increases gradually as the seismic wave propagates upward, as shown in Figure 13(a). As shown in Table 5 and Figure 13(b), the PGA increases gradually from the centre to both ends of ML 2. Table 5 and Figure 13(c) illustrate that the change of PGA on ML 3 is small.
(a)
(b)
(c)
The PGA varied with the elevation, and earthquake ground motions at the top of a valley are strengthened. This illustrates the notable topographic effect. As shown in Table 5, the PGA, PGV, and PGD on the slope are larger than those on ML 3, and these values decrease to a normal level at even a small distance away from the slope body (Li [23]). The seismic damage is more severe on the slope.
As shown in Figure 13, the PGA magnification at the top of the slope is approximately a factor of 2, thus leading to serious collapse at the slope. The real 3D complex numerical model can simulate the phenomenon of landslide and collapse at specific locations (Li [23]). The uneven distributions of ground motions on a marine site caused by topographic effects are useful for the design and construction of underwater facilities.
5. Conclusion
A new elastoplastic constitutive model with rotation hardening, which can systematically describe the monotonic and cyclic mechanical behaviours of typical marine soils combining the subloading, normal, and superloading yield surfaces, is employed in the seismic response analysis of a 3D marine site. New evolution equations for the stress-induced anisotropy and degree of overconsolidation are proposed. This model can describe the unified behaviour of unsaturated soil and saturated soil by using degree of saturation as independent state variable. By using the transform stress method, our model uniquely describes the overall mechanical properties of marine soils under general stress states, without changing the parameter values. Based on this model and three-phase field theory, an effective stress-based, fully coupled, explicit finite element–finite difference method is developed. The FEM and explicit integration method are applied in spatial and temporal discretization. An improved artificial boundary is adopted in the analysis. A finite deformation analysis is performed by introducing the Green-Naghdi rate tensor. This model is suitable for the static-dynamic analysis of complex 3D sites.
A 3D numerical model for a typical marine sites is developed. The results indicate that the proposed method can accurately simulate a seismic disaster at a marine site. The operating status before the earthquake, the damage state after the earthquake, and the damage generation process are investigated. The displacement, velocity, and acceleration time histories suggest that the ground motion intensity will increase due to local uneven complex topography and site effects. The results also demonstrate the temporal and spatial distributions of landslides and collapse at the specific marine site. Local site and topographic effects control the damage level and pattern on the marine site. The ground motion is most easily focused at the top of the slope, where the wave finally reached. At that location, the ground motion is magnified approximately two times compared to normal values and will lead to serious collapse at the slope.
Conflicts of Interest
All authors have disclosed any actual or potential conflicts of interest including any financial, personal, or other relationships with other people or organizations within three years of beginning the submitted work, which could inappropriately influence, or be perceived to influence, their work.
Acknowledgments
This work was supported by the Scientific Research Fund of Institute of Engineering Mechanics, China Earthquake Administration (Grants nos. 2016A02, 2014B03, and 2017B10); National Key Research and Development Program of China (Grant no. 2016YFC0800205); National Natural Science Foundation of China (Grants nos. 51408566 and 51438004).