Evolution Prediction of Hysteresis Behavior of Sand under Cyclic Loading

: Soil cyclic degradation is a serious issue that should be considered in engineering design and maintenance. The hysteretic response causes strength degradation and excessive settlement of soil structure in engineered and natural geosystems. Hysteresis is essentially the coupling deformation of elastic and plastic components during reloading and unloading processes. Conventional hysteretic models for sand in the elastoplastic framework rely highly on yield surface or potential surface evolution and fall short on complexity and inaccuracy. This study proposes a decoupling method to describe the hysteretic response of sand. In contrast to the conventional elastoplastic approach, this decoupling method can directly decouple the elastic and plastic components by determining the boundary between the elastic strain extension domain and the plastic strain extension domain for each stress cycle. In this way, the elastic and plastic stress–strain branches during cyclic loading can be separately obtained, and the corresponding elastic and plastic parameters are employed to characterize mechanical behavior. With the respective evolution of elastic and plastic strains, the hysteretic behavior of sand is reproduced by combining all the branches. Finally, this decoupling method is validated by three conventional cyclic loading tests. The predictions are consistent with the experiments, indicating that the decoupling method is generally effective in describing the hysteretic behavior under cyclic loading. This decoupling method provides new insight to obtain elastic and plastic deformation behaviors separately, without recourse to complicated plastic surface and hardening law.


Introduction
The cyclic degradation of soil structure caused by cyclic loads is a serious issue that should be considered in engineering design and maintenance, which has a potentially high risk in terms of safety [1,2]. Excessive settlement and excess pore pressure accumulation are the ultimate direct results of soil cyclic degradation, hence deteriorating the long-term functionality of soil structure in engineered and natural geosystems [3,4]. For example, pavement cracking from fatigue is caused by the deterioration of the subgrade under long-term cyclic loading with coupling-resilient deformation and accumulated plastic deformation [5,6]. Soil cyclic degradation can also endanger the stability of a foundation in undrained conditions, in which the effective stress decreases with the accumulation of pore water pressures due to cyclic shearing [7]. The soil cyclic degradation process is presented well by the hysteresis response. The hysteretic response is a significant mechanical behavior performance during the cyclic degradation process, which is often regarded as an extensive research object. Therefore, it is of great significance to describe the hysteresis behavior accurately in cyclic loading. regarded as an extensive research object. Therefore, it is of great significance to describe the hysteresis behavior accurately in cyclic loading.
Hysteresis constitutes coupling elastic and plastic behaviors in the unload-reload process. The conventional elastoplastic models from static analysis assume purely elastic deformation within the interior of the yield surface in response to external stress change. They work well in the static analysis but are not applicable in the cyclic analysis. One typical feature in hysteresis is the phenomenon that the unload-reload cycle with constant stress amplitude can result in a gradual accumulation of permanent strain. The different behavior in cyclic loading inevitably limits the application of the conventional elastoplastic constitutive models on the hysteretic response. For example, in classical criticalstate models, as well as in modified-state surface approach models, large plastic strains occur only on primary loading (the first loading), and elastic strains are assumed on subsequent unload-reload cycles within the yield surface [8][9][10][11]. To better describe the hysteresis response, extensive endeavors have been made by introducing the hardening modulus theories, such as multiple yield surface plasticity [12,13], kinematic hardening [14,15], sub-loading surface [16], and finally boundary surface plasticity [17][18][19] and generalized plasticity [20,21]. In those continuum-based theories, hardening modulus is associated with either specified yield surface, potential surface, or directly defined, which obscures the deformation mechanism of particular soils in hysteretic response. Particles soils are found to yield in the unloading process, with a relatively large stress increment exceeding the transient elastic limit [22][23][24]. Hence, the evaluation of the elastic parameters directly in the unloading process seriously affects the accuracy of the subsequent plastic parameter calibration. Moreover, the pseudo-elasticity obtained during the unloading process depends on the stress history, resulting in the arbitrariness of the parameters of the hardening modulus calibration. In fact, the constitutive theory focuses on the yield surface or potential surface in the elastoplastic framework, and cannot directly apply to soils compared with continuum materials [25]. Thus, the hysteretic behavior needs to be understood directly from the evolution process of elastic and plastic deformation during reloadunload cycles.
The reload-unload cycle constitutes the loading and unloading branches in each cycle, and each branch has a coupling of elastic and plastic deformations. It is noted that plastic deformation (strain) develops during unloading in a contractive nature [21]. Each axial strain increment in loading or unloading is the sum of elastic and plastic strain components, and the strain increments accumulate with increasing stress levels. Therefore, the accumulated strain can be separated into the elastic strain extension domain and the plastic strain extension domain. These two domains are demarcated by loading or unloading elastic curves. If this elastic curve is obtained, the elastic strain and plastic strain at any stress level can be readily calculated, as schematically shown in Figure 1. Combining the elastic and plastic parts of each cycle, a continuous hysteretic response is reproduced immediately. Motivated by the need to characterize the hysteresis process accurately, the objective of this study was to propose a decoupling method to accurately obtain the strain components on both sides of the elastic curve at any stress level during each stress cycle. As opposed to the aforementioned methods, the advantage of our method lies in directly describing the elastic and plastic behaviors for all branches in cyclic loading without assumptions about yielding surface expansion or hardening laws. The elastic and plastic strains are obtained directly. First, the elastic dividing curve is determined to decouple the elastic and plastic components during each loading or unloading branch. Then, empirical formulas on elastic and plastic parameters are employed to describe the corresponding variations of strains; therewith, the hysteretic behavior of sand is reproduced by combining all the branches. Finally, this decoupling method is validated by three conventional cyclic loading tests.

Determining the Elastic Curve
The elastic curve delineates the boundary between the elastic strain extension domain and the plastic strain extension domain at all stress levels (see Figure 1). The region between the elastic curve and the total stress-strain curve bounds the plastic stress-strain response. Therefore, it is significantly important to capture the elastic curve.
Determining the elastic curve is essential for obtaining the elastic response of the soil. Tatsuoka and Shibuya [26] showed that the deformation of geomaterials (soils) measured at strains smaller than about 0.001% is essentially strain-rate-independent and recoverable, i.e., elastic. Generally, it is assumed that the elastic domain exists at the strain below 0.001% [26][27][28][29], which is acknowledged in a small strain regime. However, in large strain regime, the elastic domain continues to extend beyond the 0.001% limit during shearing, and it is also coupled with plastic strain. Elastic behavior in the extended elastic domain is usually characterized by the triaxial compression test, in which a single or multiple unload-reload cycles having small amplitudes of cyclic strain with the order of 10 −5 or even 10 −6 are applied at different stress or strain points [22,[30][31][32][33]. However, this measurement cannot fully characterize the continuous nonlinear elastic behavior of the soil due to the limited stress points. Another commonly used method for testing elastic properties is dynamic testing. In the dynamic testing, the elastic properties can be evaluated by resonant column tests [34][35][36][37] or wave propagation tests [38][39][40], which often treat the test material as elastic and isotropic [41,42]. The dynamic testing generally measures the small strain elastic properties of low strain levels under the isotropic stress state, but direct application to the extended elastic domain is unreasonable. In this regard, within a large strain regime, an alternative quasi-elastic curve is obtained from the quasi-elastic domain to approximate the real elastic curve. The quasi-elasticity, rather than pure elasticity mentioned here, takes into account the energy dissipation in the nonzero area of the closed hysteresis loop. It is evidenced by experiments and numerical simulations [43][44][45][46] that the soil gradually exhibits elastic behavior under cyclic loading; therewith, the quasi-elastic domain can be obtained when the specimen reached a deformation steady state. In the deformation steady state, the stress-strain curve is a closed hysteretic loop showing the elastic behavior, and the shape of the hysteresis loop does not change dramatically in response to the subsequent stress cycles [46]. In this study, if the strain increment of two adjacent cycles is lower than 10 −5 (the normal elastic strain limit), the specimen is considered to reach the deformation steady state, and the quasi-elastic curve is obtained. It is noted that the issue of whether this quasi-elastic curve can be directly regarded as the elastic curve mentioned above still needs further demonstration regarding the influence of cyclic stress history on the evolution of elastic behavior of sand.
Extensive studies have shown that soils show a gradual hardening behavior under cyclic loading. However, few studies have focused on the evolution of elastic behavior in cyclic loading. Karg and Haegeman [3] stated that the elastic properties of the soil remain almost constant during cyclic loading within the small strain regime. Xia, et al. [47] also draw a similar conclusion based on the periodic static triaxial test that the axial elastic strain increment within small strain regime is a constant in the subsequent stress cycles. Furthermore, Xia, et al. [48] experimentally compared the differences between the quasielastic curves obtained at different deformation steady states, showing that the quasi-elastic curves are not sensitive to the influence of cyclic stress history. Pradhan, et al. [49] directly assume that elasticity remains constant during cyclic loading based on simplification considerations. Similarly, the same assumption is made here, that the elastic behavior is independent of the loading cycles. Therefore, the quasi-elastic curve in the deformation steady state can be shifted to any stress cycle, as shown in Figure 2.
almost constant during cyclic loading within the small strain re draw a similar conclusion based on the periodic static triaxial strain increment within small strain regime is a constant in the Furthermore, Xia, et al. [48] experimentally compared the differ elastic curves obtained at different deformation steady states, sho tic curves are not sensitive to the influence of cyclic stress his directly assume that elasticity remains constant during cyclic lo cation considerations. Similarly, the same assumption is made he ior is independent of the loading cycles. Therefore, the quasi-e mation steady state can be shifted to any stress cycle, as shown i

Decoupled Plastic Stress-strain Curve
The total stress-strain relation is decoupled into elastic and matically shown in Figure 3. The elastic branch is obtained direc steady state and remains unique. Thereafter, the plastic stressobtained by subtracting from the total stress-strain relation. It plastic branch in different cycles is different depending on the elastic and plastic parameters are employed to describe the elasti ments, respectively. Elastic parameters are the (axial) elastic mod (different from Poisson's ratio); plastic parameters are (axial) pla tic strain ratio vp. They are defined as the tangent slope of the c shown in Figure 3. Besides, elastic bulk modulus K and plastic b in isotropic cyclic loading conditions.

Decoupled Plastic Stress-Strain Curve
The total stress-strain relation is decoupled into elastic and plastic branches, as schematically shown in Figure 3. The elastic branch is obtained directly from the deformation steady state and remains unique. Thereafter, the plastic stress-strain relation is readily obtained by subtracting from the total stress-strain relation. It should be noted that the plastic branch in different cycles is different depending on the hardening behavior. The elastic and plastic parameters are employed to describe the elastic and plastic strain increments, respectively. Elastic parameters are the (axial) elastic modulus E and strain ratio v (different from Poisson's ratio); plastic parameters are (axial) plastic modulus E p and plastic strain ratio v p . They are defined as the tangent slope of the corresponding curves, as shown in Figure 3. Besides, elastic bulk modulus K and plastic bulk modulus K p are used in isotropic cyclic loading conditions.  Figure 4 shows the procedures of the decoupling method. Axial incremental s ∆ε is the sum of the elastic strain increment ∆εe and plastic strain increment ∆εp. The step is to obtain the alternative quasi-elastic curve from the deformation steady state alternative elastic strain increment at the deformation steady state is denoted as (∆εe each stress increment. Then, shift the quasi-elastic curve to the origin of each cycle. ∆ approximated at each stress level by (∆εe)′. The corresponding plastic strain can be c lated as ∆εp = ∆ε − (∆εe)′. The separated elastic and plastic strain increments are then pressed as the function of elastic and plastic parameters. Immediately, total strains ca achieved by integrating the elastic and plastic parameters over the input stress his Finally, evolution prediction of hysteresis behavior of sand under cyclic loading is re duced by combining all the branches. This data iteration process can be implemented ily through Excel. The following sections demonstrate the effectiveness of the decoup method considering different cyclic loading conditions.  Figure 4 shows the procedures of the decoupling method. Axial incremental strain ∆ε is the sum of the elastic strain increment ∆ε e and plastic strain increment ∆ε p . The first step is to obtain the alternative quasi-elastic curve from the deformation steady state. The alternative elastic strain increment at the deformation steady state is denoted as (∆ε e ) for each stress increment. Then, shift the quasi-elastic curve to the origin of each cycle. ∆ε e is approximated at each stress level by (∆ε e ) . The corresponding plastic strain can be calculated as ∆ε p = ∆ε − (∆ε e ) . The separated elastic and plastic strain increments are then expressed as the function of elastic and plastic parameters. Immediately, total strains can be achieved by integrating the elastic and plastic parameters over the input stress history. Finally, evolution prediction of hysteresis behavior of sand under cyclic loading is reproduced by combining all the branches. This data iteration process can be implemented easily through Excel. The following sections demonstrate the effectiveness of the decoupling method considering different cyclic loading conditions.

Elastic Parameters E and v
Due to the difficulty in reaching a complete elastic state, many alternative el rameters are defined to approximate the elastic properties [46]. To quantify so mation, an accurate description of E and its dependence on influencing factors sary. An empirical formulation for this dependent relation is an effective wa widely used. Some researchers [50][51][52] have studied E and proposed different e formulations for its dependency on influencing factors. Generally, it is believed defined for major elastic principal strain increments in a certain direction, is a function of the normal stress in that direction [22,32,53,54].
In this study, a new representation of E is used to formulate the elastic stres relation in standard triaxial (axisymmetric) stress conditions, which is formulat the deformation steady state [47,55], as follows: where E0 is a constant which is related to the initial state, σ1 is the axial stress, material parameter, and pa is the atmospheric pressure. In comparison with E, the empirical formula on v is rare. It is believed that v is not a constant but increases

Elastic Parameters E and v
Due to the difficulty in reaching a complete elastic state, many alternative elastic parameters are defined to approximate the elastic properties [46]. To quantify soil deformation, an accurate description of E and its dependence on influencing factors is necessary. An empirical formulation for this dependent relation is an effective way and is widely used. Some researchers [50][51][52] have studied E and proposed different empirical formulations for its dependency on influencing factors. Generally, it is believed that E, defined for major elastic principal strain increments in a certain direction, is a unique function of the normal stress in that direction [22,32,53,54].
In this study, a new representation of E is used to formulate the elastic stress-strain relation in standard triaxial (axisymmetric) stress conditions, which is formulated from the deformation steady state [47,55], as follows: where E 0 is a constant which is related to the initial state, σ 1 is the axial stress, α is the material parameter, and p a is the atmospheric pressure. In comparison with E, the existing empirical formula on v is rare. It is believed that v is not a constant but increases with the stress level in the shear direction [32,[56][57][58], sometimes even exceeding 0.5 [56,58]. It is observed that the strain ratio increases linearly with the deviatoric stress in the deformation steady state [59], which is composed of two parts, as follows: where v 0 is the initial strain ratio, β is the coefficient, and q = σ 1 − σ 3 is deviatoric stress. The initial strain ratio was determined by confining pressure and initial void ratio, while the incremental part is related to deviatoric stress. Thus, four elastic parameters can be obtained from the deformation steady state.

Plastic Parameters E p and v p
As important as elastic parameters, both E p and v p play key roles in accurately describing plastic responses. According to the decoupled plastic branches, in standard triaxial stress condition, E p in reloading can be expressed as [55]: where q u is reloading ultimate deviatoric stress, b (<0) is reloading coefficient, h(N) is the reloading plastic modulus hardening function in which δ is a constant and N is the number of cycles. On the other hand, although the plastic behavior in the reloading is different from that in the unloading, the plastic modulus changes linearly with the deviator stress, that is, it decreases linearly in the reloading and increases linearly in unloading, but corresponds to different parameters, as follows: where q u is unloading ultimate deviatoric stress, b (>0) is the unloading coefficient, and h (N) is the unloading plastic modulus hardening function in which δ is a constant. The existing study on v p is limited. By definition, the plastic strain ratio is calculated by plastic strains. Generally, in continuum mechanics, it is assumed that the plastic volume is incompressible, and the corresponding v p is a constant of 0.5 [60]. However, this is not accurate for porous materials, because the plastic deformation is mainly due to the changes in the relative position between the particles, which will inevitably cause changes in the void volume; consequently, v p of porous materials is not constant. Few studies have reported v p of microporous metallic foams; however, v p of granular materials (soils) has not been reported. In this paper, in standard triaxial stress condition, the empirically formulated v p in reloading is expressed as follows: where v 0,p is the initial plastic strain ratio in reloading and p ref is the reference pressure (in this paper p ref = 1 kPa), µ, ω, and λ are constants, and g(N) is the plastic strain ratio hardening function in reloading. Moreover, the plastic strain ratio in unloading is expressed by the same formulation, as follows: where v 0,p is the initial plastic strain ratio in unloading, µ , ω , and λ are constants, and g(N) is the plastic strain ratio hardening function in unloading.

Simulation
By combining Equations (1), (3) and (4), the axial stress-strain relations are reproduced, and the radial-axial strain relation is reproduced by combining Equations (2), (5) and (6). The empirical model parameters are calibrated through the Solver module in Excel, where the solution iteration method is a nonlinear GRG. The details of the experimental results selected can be found in [46]. The simulations and experimental results are shown in Figure 5. The predictions are consistent with the experiments, indicating that the decoupling method is effective in describing the hysteretic behavior under cyclic triaxial stress conditions. The constants of the empirical model are shown in Table 1.

Application in the Cyclic Isotropic Compression Test (CICT)
Isotropic compression tests are also very common, providing a continuous relati between deformation modulus (bulk modulus) and pressure (confining pressure). In t cyclic isotropic compression test (CICT), the hysteresis loops indicate that, in both reloa ing and unloading, the deformation is a coupling of elastic and plastic components,

Application in the Cyclic Isotropic Compression Test (CICT)
Isotropic compression tests are also very common, providing a continuous relation between deformation modulus (bulk modulus) and pressure (confining pressure). In the cyclic isotropic compression test (CICT), the hysteresis loops indicate that, in both reloading and unloading, the deformation is a coupling of elastic and plastic components, as shown in Figure 6. The elastic and plastic components are calculated using the elastic bulk modulus K and plastic bulk modulus K p , respectively. Additionally, the K is obtained from the elastic curve corresponding to the tangent slope of the elastic branch at the deformation steady state, while the K p is obtained from the decoupled plastic stress-strain curve corresponding to the tangent slope of the plastic branch. To present the performance of this decoupling method in describing the hysteresis response during the CICT, the empirical formulas for K and K p need to be formulated first. Researchers have proposed some empirical formulas for volumetric stiffness. For example, Qubain, et al. [61] proposed a quasi-linear elastic constitutive model to describe the behavior of sand well below failure based on isotropic compression tests. García and Medina [62] performed detailed simulations of stress-strain relations in unconsolidated granular packs under unload-reload cycles, in which the power function was used to describe the bulk elastic modulus which scales with pressure with a 1/2 power law exponent in the limited cycle (i.e., deformation steady state). Figure 7a shows the variation of K with pressure, where (K/pa) 2 changes linearly with mean pressure. Therefore, K is expressed as: where p is mean pressure (p = 1/3(σ1 + 2σ3)), and σ1, σ1, and σ3 are principal stresses. In the isotropic compression test, p = σ1 = σ2 = σ3 holds. ϕ and φ are constants. The details of the experiment for the cyclic isotropic compression test are presented in [63]. Researchers have proposed some empirical formulas for volumetric stiffness. For example, Qubain, et al. [61] proposed a quasi-linear elastic constitutive model to describe the behavior of sand well below failure based on isotropic compression tests. García and Medina [62] performed detailed simulations of stress-strain relations in unconsolidated granular packs under unload-reload cycles, in which the power function was used to describe the bulk elastic modulus which scales with pressure with a 1/2 power law exponent in the limited cycle (i.e., deformation steady state). Figure 7a shows the variation of K with pressure, where (K/p a ) 2 changes linearly with mean pressure. Therefore, K is expressed as: where p is mean pressure (p = 1/3(σ 1 + 2σ 3 )), and σ 1 , σ 1 , and σ 3 are principal stresses. In the isotropic compression test, p = σ 1 = σ 2 = σ 3 holds. φ and ϕ are constants. The details of the experiment for the cyclic isotropic compression test are presented in [63].
On the other hand, the decoupled bulk plastic modulus is expressed as: where ψ and ξ are reloading bulk constants. Similarly, the bulk plastic modulus in unloading is expressed as: where ψ and ξ are unloading bulk constants. The variation of bulk plastic modulus with pressure is plotted in Figure 7b. By combining Equations (8) and (9), the prediction is presented in Figure 8, which is consistent with the experimental results. The constants of the empirical model are shown in Table 2. Note that the large deviation has occurred during the final phase of the first unloading cycle, where the prediction highly underestimates the experiment. The possible reason is that the unloading empirical model has large errors at low stress levels in the primary unloading, resulting in smaller unloading strain. However, this deviation gradually weakens with increasing cycle number, indicating that the unloading empirical model also couples with the effect of cyclic stress history.
Researchers have proposed some empirical formulas for volumetric stiffness. For example, Qubain, et al. [61] proposed a quasi-linear elastic constitutive model to describe the behavior of sand well below failure based on isotropic compression tests. García and Medina [62] performed detailed simulations of stress-strain relations in unconsolidated granular packs under unload-reload cycles, in which the power function was used to describe the bulk elastic modulus which scales with pressure with a 1/2 power law exponent in the limited cycle (i.e., deformation steady state). Figure 7a shows the variation of K with pressure, where (K/pa) 2 changes linearly with mean pressure. Therefore, K is expressed as: where p is mean pressure (p = 1/3(σ1 + 2σ3)), and σ1, σ1, and σ3 are principal stresses. In the isotropic compression test, p = σ1 = σ2 = σ3 holds. ϕ and φ are constants. The details of the experiment for the cyclic isotropic compression test are presented in [63]. On the other hand, the decoupled bulk plastic modulus is expressed as: where ψ and ξ are reloading bulk constants. Similarly, the bulk plastic modulus in unloading is expressed as: where ψ′ and ξ′ are unloading bulk constants. The variation of bulk plastic modu pressure is plotted in Figure 7b. By combining Equations (8) and (9), the predictio sented in Figure 8, which is consistent with the experimental results. The constan empirical model are shown in Table 2. Note that the large deviation has occurred the final phase of the first unloading cycle, where the prediction highly undere the experiment. The possible reason is that the unloading empirical model has larg at low stress levels in the primary unloading, resulting in smaller unloading strai ever, this deviation gradually weakens with increasing cycle number, indicating unloading empirical model also couples with the effect of cyclic stress history.

Application in Cyclic Oedometric Compression Tests (COCT)
The oedometric compression test can be regarded as a special triaxial test i the confining pressure is variable in proportion to the axial stress. The elastic and moduli presented in triaxial stress conditions are employed herein. Note that t plastic modulus includes the influence of the confining pressure but is expressed axial stress, as follows:

Application in Cyclic Oedometric Compression Tests (COCT)
The oedometric compression test can be regarded as a special triaxial test in which the confining pressure is variable in proportion to the axial stress. The elastic and plastic moduli presented in triaxial stress conditions are employed herein. Note that the axial plastic modulus includes the influence of the confining pressure but is expressed by the axial stress, as follows: where b 1 , b 2 , and n are the reloading constants, and b 1 , b 2 , and n are the unloading constants. The performance of the empirical formulas in the oedometric compression test is shown in Figure 9, indicating that the decoupling method is also effective in describing the hysteretic behavior under cyclic oedometric compression stress condition. The details of the experiment for the cyclic oedometric compression test are presented in [47]. The constants of the empirical model are shown in Table 3.

Limitations of the Decoupling Method
Although this decoupling method abandons the yield surface as the study object and can accurately describe the strains in each cycle, there are still some limitations. It can be seen from the above verification cases that the decoupling method is suitable for a single changing stress condition (here, the axial direction), such as in CTT, but the prediction results do not match well with the experimental results in the cases with changing multidirectional stress. The possible reason is that the empirical formulas established for the elastic and plastic parameters are not the most effective. Therefore, the deviation of prediction is closely related to the determination of the quasi-elastic curve, as well as the empirical formulas of elastic and plastic parameters for each branch. Considering more complex cyclic loading patterns is necessary to further verify the effectiveness of the decoupling method in the future.

Implications of This Work
By the decoupling method, the evolution of plastic strain with stress level and stress history is directly presented, which has implications for the development of elastic-plastic theory in assuming the hardening laws for soils. The descriptions of the hysteresis process can also facilitate the understanding of the soil shakedown limit analysis. Generally, the most likely cyclic stress is selected as the plastic shakedown limit by observing the trend of permanent strain accumulation, during which the reference and judgment criteria for selection are extremely important. At present, the reference or judgment criteria are em-

Limitations of the Decoupling Method
Although this decoupling method abandons the yield surface as the study object and can accurately describe the strains in each cycle, there are still some limitations. It can be seen from the above verification cases that the decoupling method is suitable for a single changing stress condition (here, the axial direction), such as in CTT, but the prediction results do not match well with the experimental results in the cases with changing multidirectional stress. The possible reason is that the empirical formulas established for the elastic and plastic parameters are not the most effective. Therefore, the deviation of prediction is closely related to the determination of the quasi-elastic curve, as well as the empirical formulas of elastic and plastic parameters for each branch. Considering more complex cyclic loading patterns is necessary to further verify the effectiveness of the decoupling method in the future.

Implications of This Work
By the decoupling method, the evolution of plastic strain with stress level and stress history is directly presented, which has implications for the development of elastic-plastic theory in assuming the hardening laws for soils. The descriptions of the hysteresis process can also facilitate the understanding of the soil shakedown limit analysis. Generally, the most likely cyclic stress is selected as the plastic shakedown limit by observing the trend of permanent strain accumulation, during which the reference and judgment criteria for selection are extremely important. At present, the reference or judgment criteria are empirical and variable, which is inseparable from the fact that researchers only focus on the final result of the permanent strain accumulation and ignore the mechanical response during the shakedown process in cyclic loading. Due to the lack of understanding of the evolution of the plastic mechanical behaviors with stress cycles, the mechanical properties of the hysteresis response or shakedown process cannot be revealed, in essence, which in turn leads to the lack of a unified and effective criterion when evaluating the plastic shakedown limit. With the decoupling method presented in this study, research on this topic has more possibilities. Nevertheless, further verification work is needed.

Conclusions
In this study, a decoupling method for describing the hysteresis response of sand is proposed, in which the elastic and plastic strain components are directly decoupled by the boundary between the elastic strain extension domain and the plastic strain extension domain. This boundary can be obtained directly in a deformation steady state. From the performance in three application cases, the decoupling method is generally effective in describing the hysteretic behavior under cyclic loading. However, the deviation of prediction is closely related to the determination of the quasi-elastic curve, as well as the empirical formulas of elastic and plastic parameters for each branch. Considering more complex cyclic loading patterns is necessary to further verify the effectiveness of the decoupling method in the future.

Data Availability Statement:
The cyclic loading test data used to support the findings of this study are available from the corresponding author upon reasonable request.