#### Abstract

A standard stress path triaxial test system was applied to carry out conventional triaxial shearing tests for gravelly sands under confining pressures ranging from 50 kPa to 400 kPa at the initial relative densities of 0.15, 0.35, 0.55, and 0.75, respectively. The test results show that all the samples of gravelly sand present strain hardening and shear contraction during the process of shearing test. Additionally, gravelly sands are significantly affected by the initial relative density. The hardening degree of gravelly sand samples rises in line with increasing initial relative densities during shearing tests. When initial relative densities *D*_{r} are at 0.15 and 0.35, the volume shrinkage of samples decreases with the increasing confining pressures. Instead, when initial relative densities *D*_{r} are at 0.55 and 0.75, the volume shrinkage of samples increases with the growth of confining pressures. To describe these triaxial shearing mechanical properties of gravelly sands, a higher-order dilatancy equation was proposed based on the concept of a super yield surface. A constitutive model which can describe the mechanical properties of gravelly sand was established when the associated flow laws were applied to compare with the results of the triaxial shearing test under the consolidated drained condition. The comparison results showed that the proposed model can reflect the strain hardening and shear contraction characteristics of gravelly sands from low to high confining pressures under different initial relative densities.

#### 1. Introduction

The Qinghai-Tibet Highway (QTH) covers 550 km across the Plateau Permafrost Region (PPR) with harsh ecological environment and complex geological conditions. Frequent highway damage such as thaw subsidence, wave terrain, and longitudinal cracks develops along the highway. Even with perennial maintenance, the area of subgrade subsidence on the highway accounts for approximately 19% of the total area, and the depths of subsidence for several sections of the highway exceed 20 cm, with 60 cm in maximum. The subgrade subsidence significantly affects the safety of QTH traffic [1]. The Qinghai-Tibet Expressway (QTE), the only national expressway network in the Tibet Autonomous Region, is the last section of the Beijing-Tibet Expressway subordinated to the “71118 Expressway Network” of China. However, owing to the insufficient support system for expressway construction technology in PPR, the construction of the expressway in the Tibetan Plateau does not start. Therefore, an in-depth study on countermeasure for subgrade subsidence is urgently needed. Long-term engineering practices suggested that the key for the successful construction of the QTE lies in understanding and controlling subgrade damage. Hence, it is necessary to study the mechanical properties of gravelly sand, a subgrade material in the Qinghai-Tibet Plateau.

An understanding of the strength and dilatancy of sand is critical for the study of strength and deformation in soil mechanics [2]. With experimental studies, researchers have examined the impacts of multiple factors on the mechanical properties of sandy soil, such as confining pressure, temperature, loading rate, stress path, and initial relative density. Based on the results of drained triaxial tests with confining pressures changing from 0.2 to 6.4 MPa, Lu et al. concluded that stress-strain characteristics of sand changed from strain softening to strain hardening, and the dilatancy characteristics transformed from dilatancy to contraction [3]. On the basis of confining pressures ranging from 25 to 150 MPa, Martin et al. studied the mechanical properties of sandy soil. The results showed that the strength and strain hardening degree of the sandy soil increased with confining pressures [4]. A series of triaxial compression experiments on sandy soil from –0.5 to –6°C were conducted by Lai et al. The test results showed that the strength of frozen sandy soil had significant temperature dependence: it increased as temperature decreased [5]. With exploring the mechanical behaviors of sandy soil via triaxial compression tests with different stress paths, Guo and Han found the peak strength of sandy soil successively reduced under the stress paths of CTC, TC, and RTC, and the dilatancy decreased successively in stress paths of RTC, TC, and CTC [6]. Some scholars have also explored the influence of stress path on the compression and particle fragmentation of sand soil [7]. Chai et al. found that the shear strength of calcareous sand decreased first and then increased in a shear rate from 0.1 to 2.4 mm/min in direct shear tests [8]. A shear test on frozen sandy soil by Gao et al. at different loading rates (0.005%, 0.010%, 0.015%, and 0.020% min^{−1}) indicated that as loading rates increased, the peak strength of the frozen sandy soil increased; instead, the failure strain decreased [9]. With drained triaxial tests on saturated sandy soils in the laboratory, Zhu et al. found that as the initial relative densities *D*_{r} increased, the brittleness of the sandy soil increased, the failure strain decreased, and the strength increased in different degrees [10]. These achievements contributed to our understanding of the mechanical properties of sandy soil.

Based on the nonassociated flow law, Lade developed the Lade–Duncan model suitable for sandy soil [11]. On the basis of the Modified Cambridge Model [12], Yao et al. revised the hardening parameter and proposed a unified hardening model [13, 14]. They then extended the application of this unified hardening model to sandy soil [15]. In the framework of critical state soil mechanics, Been et al. proposed a state-dependent dilatancy theory, where a state parameter is defined as the difference between the void ratio of the current effective average normal stress and the critical void ratio. And the state parameter can reflect the current state of sandy soil [16]. Gajo and Wood introduced the state parameter *ψ* into the dilatancy equation and established a seven-trend model for sandy soil [17]. Through introducing the state parameter into the dilatancy equation of the boundary surface model, Li et al. established an elastoplastic constitutive model for sandy soil [18]. Dafalias and Manzari subsequently introduced the state parameter into the constitutive model to develop an elastoplastic constitutive model with double yield surfaces, which was applied to simulate monotonic loading characteristics of sandy soil under different drainage conditions [19]. Given the critical state theory and material state-dependent dilatancy theory, Huang et al. established a double yield surface constitutive model for sandy soil with different initial relative densities [20]. In the introduction of the state parameter into the constitutive model, Cai and Li proposed a critical state constitutive model for sandy soil based on state-dependent dilatancy theory [21, 22]. With introducing void ratio as a state parameter into the modified Rowe dilatancy equation, Wan and Guo developed a critical state constitutive model for sandy soil [23]. The foregoing constitutive models are widely used in simulating the stress-strain characteristics and dilatancy of sandy soil. They also provide ideas for other scholars to establish constitutive model of sand soil considering various influencing factors [24].

It is known that the mechanical properties of geotechnical materials are significantly affected by the physical constituents of the materials. Characteristics of particle breakage and granular constitution of local sandy soil are not identical to those of sandy soils in other regions due to the strong solar radiation and the long-term freeze-thaw cycles. Additionally, ice in the frozen soil melts in the warmer seasons, thereby forming thawed soil with high water content, which influences the mechanical properties of sandy soil. In view of such circumstances and considering the theory for prevention and control of subgrade damage in permafrost region, this paper takes the subgrade material of QTH as research object. The effects of the initial relative density on the mechanical properties of gravelly sand were investigated based on conventional triaxial compression tests under drained condition. Material parameters were applied to modify the dilatancy equation of the Modified Cambridge Model to construct a constitutive model for gravelly sand in the Qinghai-Tibet Plateau, considering the effects of its initial relative density. Finally, the validity of the model was verified according to triaxial test results. The results are of great significance: they reveal the subsidence mechanisms of the frozen soil subgrade and provide a reference for the construction and control of subgrade failure of the QTE.

#### 2. Test Materials and Methods

##### 2.1. Test Materials and Soil Sample Preparation

The test soil was from the Na Qu section. During field sampling, it was found that the ice particles in the pores of frozen sandy soil melted under the influence of rising temperatures in the warm season, thereby forming melted soil with high moisture content, as shown in Figure 1. The fresh soil in the photograph is pale yellow in color. During the sieving process, it was found that the soil sample contained gravel particles with size greater than 2 mm. Therefore, these soil samples were distinguished to be gravelly sand [25]. The particle grading curve of gravelly sand is shown in Figure 2, and its physical parameters are listed in Table 1.

Test samples are prepared according to a standard procedure followed by the Specification of Soil Tests (GB/T50123-1999) issued by the Ministry of Water Resources, People’s Republic of China. Before a test, the sample preparation process is described below [26]. First, the air-dried gravelly sand was sifted soil particles with diameter ranging from 0.075 to 10 mm using sieve meshes. Then, the different qualities of gravelly sands were prepared according to their initial relative densities *D*_{r} (0.15, 0.35, 0.55, and 0.75). Next, the soil samples were filled in copper molds of size 50 × 100 mm. The samples were compacted in three layers, and the height of each layer was strictly controlled at 33.4 mm. The initial moisture content *ω*, saturated water content *ω*_{s}, total mass *m*, and filling quality of each layer are listed in Table 2.

##### 2.2. Test Apparatus and Methods

The tests were performed using a standard stress path triaxial test system produced by the company GDS, UK, as shown in Figure 3. The testing apparatus comprised a triaxial pressure chamber, cell pressure controller, back pressure controller, axial pressure controller, and sensor. The ranges of confining pressure and back pressure were 0–1.3 MPa and 0–3 MPa, respectively. The values of maximum axial load and maximum axial displacements were 7 kN and 25 mm, respectively. The test process is described as follows [27].

(1) The soil sample was first placed in the base of the pressure chamber. Water was then injected into the triaxial pressure chamber until the top of the sand sample was submerged. (2) Soil sample saturation process (carbon dioxide saturation, hydraulic saturation, and back pressure saturation) was as follows. First, a confining pressure of 20 kPa was applied to the soil sample using the cell pressure controller. Next, the carbon dioxide cylinder was opened and CO_{2} passed slowly through the sand sample from the bottom to the top at a pressure of 8 kPa to emerge from the top. This process lasted for approximately 60 min. After the carbon dioxide saturation process, distilled water was injected into the sample from the bottom to the top, and carbon dioxide was squeezed out of the pores by the distilled water. Hydraulic saturation was considered completed when the water inflow and water output were identical. Next, the confining pressure in the triaxial pressure chamber was increased by 30 kPa to conduct B detection on the soil sample (B value is the ratio of pore pressure increment to confining pressure increment, which can reflect the current saturation degree of the soil sample). If the B value did not reach 0.95, back pressure saturation was performed by applying a back pressure, which was 20 kPa lower than the current confining pressure. After the back pressure volume was stable, the foregoing B detection and back pressure saturation process were repeated until the B value reached 0.95. Then, the soil sample saturation process was completed. (3) Consolidation of soil sample: the soil samples were consolidated after the B value reached 0.95. The consolidation pressure values were set to 50, 100, 200, and 400 kPa. During the consolidation process, the drainage valve was opened, and consolidation was considered completed when the back pressure volume changed by <5 mm^{3}in 5 min. (4) Loading: the conventional triaxial stress path was used to load the samples under the condition of consolidated drainage. The load path is illustrated in Figure 4. The loading rate was set to 0.024 mm/min according to the Specification of Soil Tests (GB/T 50123-1999). The test was terminated when the axial strain reached 15%. (5) Unloading: after the test was completed, the pressure was unloaded according to the principle “from the inside to the outside.”

**(a)**

**(b)**

#### 3. Results and Discussion

Figure 5 illustrates triaxial shear test results of gravelly sand with initial relative densities of 0.15, 0.35, 0.55, and 0.75 and confining pressures of 50, 100, 200, and 400 kPa, respectively. As shown in Figure 5(a), all the samples of gravelly sand exhibited strain hardening during the process of shearing test. Additionally, the degree of hardening and peak stress increased with an increase in initial relative density under a certain confining pressure. Figure 5(b) shows that all the samples of gravelly sand exhibited shear contraction during the shearing process. Under a certain confining pressure, all the volume shrinkage of gravelly sand samples decreases with the increase of initial relative densities. When initial relative densities *D*_{r} were 0.15 and 0.35, the volume deformation of samples decreased with the increasing confining pressures. Instead, when initial relative densities *D*_{r} were 0.55 and 0.75, the volume deformation of samples rose with the growth of confining pressures. Through comparison, it was found that the test results of gravelly sand were consistent with the results of other scholars [28]. All the samples of sand present strain hardening and shear contraction in the range of 30 kPa to 100 kPa.

**(a)**

**(b)**

#### 4. Constitutive Model of Gravelly Sand

All the stresses considered in this study are effective stresses. The effective average normal stress *p*′, shear stress *q*, and shear stress ratio *η* can be defined as follows:where *σ*_{1}, *σ*_{2}, *σ*_{3} are principal stresses; are effective principal stresses.

##### 4.1. Higher-Order Dilatancy Equation

Figure 6 shows a schematic for the dilatancy of gravelly sand. The generalized shear stress is equal to the deviatoric stress under triaxial compression conditions. As illustrated in Figure 6(a), the deviatoric stress *q* and the volumetric strain in gravelly sand increased with the axial strain in the initial loading stage. After reaching the point , the deviatoric stress *q* and volumetric strain become stable and approach horizontal curves as the axial strain *ε*_{1} increases [29]. Figure 6(b) presents the relationship between the plastic dilatancy factor *d* and the shear stress ratio *η* (*q*/*p*′) of gravelly sand. As evident, the dilatancy equation of gravelly sand is nonlinear. Starting from the current stress point A, the shear stress ratio *η* increased gradually, and the plastic dilatancy factor *d* decreased gradually with loading. Subsequently, stress point moved from point A to point B gradually and being into the critical state where the soil shows irregular plastic shear deformation, its plastic volumetric strain increment is equal to 0, and shear stress ratio *η* is equal to *M* [30]. During the loading process, the plastic dilatancy factor *d* is greater than 0, and the volume deformation in gravelly sand underwent shrinkage.

**(a)**

**(b)**

With the recommendation of K. H. Roscoe, the plastic work increment of normal consolidated remolding soil during shear failure satisfies the following equation [31]:

In the Cam-Clay model, and when shear failure occurs in the soil sample. Substituting these into equation (4) yieldswhere *M* is the critical state ratio.

The above equation can be used to give the dilatancy equation of the Cam-Clay model:

The Cam-Clay model cannot reflect the idea that the plastic shear strain increment is 0 in isotropic consolidation. To overcome this situation, Roscoe and Burland modified the energy equation of the Cam-Clay model [12]:

Equation (7) was used to obtain the dilatancy equation of the modified Cam-Clay model as follows:

Based on the above studies, we modified equation (8) to obtain a new high-order dilatancy equation:where *s* and *t* are the material parameters. To consider the influence of initial relative density *D*_{r}, we assumed that *s* and *t* are functions of initial relative density *D*_{r}.

From equation (9), we can assert the following:(1)When *s* = 2 and *t* = 1, the proposed higher-order dilatancy equation can be reduced to equation (8) of the modified Cam-Clay model(2)When *s* = 1 and *t* = 0.5, the proposed higher-order dilatancy equation can be reduced to equation (6) of the Cam-Clay model

The effects of the material parameters *s* and *t* on the proposed high-order dilatancy equation (equation (9)) are shown in Figure 7. Figure 7(a) illustrates the relationship between parameter *s* and the plastic dilatancy factor *d*. Figure 7(a) shows that, at the same shear stress ratio *η*, the plastic dilatancy factor *d* increases as parameter *s* increases. The shear shrinkage tendency of gravelly sand decreases as parameter *s* increases. Figure 7(b) shows the influence of parameter *t* on the plastic dilatancy factor *d*. It can be observed that, at the same shear stress ratio *η*, the plastic dilatancy factor *d* of the soil decreases gradually as parameter *t* increases. The shear shrinkage tendency of gravelly sand increases with parameter *t*.

**(a)**

**(b)**

##### 4.2. Higher-Order Yield Function

Different types of soils have different yield surface shapes. For example, Kaolin soil has an oval yield surface, Weald clay has a round yield surface, and the yield surface of London clay is hyperbolic [32]. In this study, the hypothesis of the modified Cam-Clay model was applied to determine the yield surface of gravelly sand within the framework of critical state soil mechanics. In the *p*-*q* coordinate system, the principal axis of stress was in the same direction as that of the plastic strain increment, which was orthogonal to the plastic potential surface of the material satisfying the normality rule. Based on this, a new plastic potential function was obtained by integrating the higher-order dilatancy equation of gravelly sand; the higher-order yield function of gravelly sand can be obtained while combining the associated flow rules with the plastic potential function. The derivation of the higher-order yield function is as follows.

Roscoe and Burland reported that the integral form of the plastic potential function of the modified Cambridge model can be expressed as [12]where *p*_{0} is the reference effective mean principal stress.

Substituting (9) into (10) yields

When the shear stress ratio *η* = 0, the effective average normal stress *p*′ = *p*_{x}, where it is the effective mean principal stress at the intersection of the plastic potential function and *p*-axis. By transforming equation (11), we can obtain

Equation (12) can be used to express the plastic potential function as follows:

Combining with the associated flow rules, the higher-order yield function of gravelly sand can be obtained from equation (13):

The higher-order yield function equation (14) has a similar form as the modified Cam-Clay model and contains two material parameters *s* and *t* related to the initial relative density. When *s* = 2 and *t* = 1, equation (14) can be directly reduced to the yield function of the modified Cam-Clay model When *s* = 1 and *t* = 0.5, equation (14) can be reduced to the yield function of the Cam-Clay model

The effects of parameters *s* and *t* on the shape of the high-order yield surface are demonstrated in Figure 8. Figure 8(a) shows the relationship between parameter *s* and the high-order yield surface. With an increase in parameter *s*, the higher-order yield surface expanded, and the shape of the yield surface gradually develops from a water drop to an ellipse. Figure 8(b) shows the relationship between parameter *t* and the high-order yield surface. With an increase in parameter *t*, the yield surface contracted to the left, and the effective average normal stress *p*′ associated with the intersection of the yield surface and critical state line decreased gradually.

**(a)**

**(b)**

When the yield function was established with different dilatancy equations, the trajectories of the yield surfaces were different. Figure 9 illustrates the *p*-*q* coordinate system for the higher-order yield surface and other yield surfaces. As shown in the figure, the trajectory of the higher-order yield surface is between the yield surface of the Cam-Clay model and the modified Cam-Clay model.

##### 4.3. Hardening Parameter

In the constitutive model, hardening parameters are used to determine the magnitude of plastic strain increment due to a given increment in the stress [33]. Conventional hardening parameters include the plastic work [11], plastic volumetric strain [12], and the uniform hardening parameter *H* [34]. The results of triaxial shear tests for gravelly sand state that the plastic volumetric strain is taken as the hardening parameter of gravelly sand. The hardening parameters were determined as follows.

Figure 10 shows the test results for gravelly sand loading and unloading in the *e*-ln *p* coordinate system. As shown, the slopes of the loading and unloading curves are *λ* and *k*, respectively. When the load increased from *p*_{0} to *p*_{x}, the pore ratio changed bywhere *e*_{0} is the initial pore ratio of the sample.

The volumetric strain comprises the elastic volumetric strain and the plastic volumetric strain :

From Figure 10, the volumetric strain in gravelly sand can be expressed as follows:

The elastic volumetric strain can be determined via an unloading test:

By substituting equations (17) and (18) into equation (16), we can determine the plastic volume strain as follows:

From equation (19), the following expression can be derived:

By substituting equation (20) into equation (14), we can obtain the higher-order yield function:

From equation (21), we have

##### 4.4. Constitutive Relationship

###### 4.4.1. Elastic Strain Increment

The elastic strain increment of the proposed model is identical to that of modified Cam-Clay model. In the modified Cam-Clay model, the elastic strain increment comprises the elastic volumetric strain increment and elastic shear strain increment , as follows:

The elastic strain increment is calculated using generalized three-dimensional Hooke’s law [35, 36] and can be expressed as follows:

By assuming and substituting equation (26) into equations (24) and (25), we obtainwhere *E* represents the elastic modulus and represents Poisson’s ratio [33].

###### 4.4.2. Plastic Strain Increment

The plastic strain increment of the proposed model consists of the plastic volumetric strain increment and plastic shear strain increment :

By assuming and substituting *C*_{p} into equation (21), the yield function *f* of the gravelly sand can be obtained as follows:

By integrating the yield function in equation (29), we obtain

The derivative of the plastic volumetric strain is calculated using equation (29):

Substituting equation (31) into equation (30) yields

The derivative of the effective mean principal stress *p* is calculated using equation (29):

The derivative of the shear stress *q* is calculated using equation (29):

Substituting (33) and (34) into (32) yields

By substituting equation (35) into equation (9), we can determine the increment in the plastic shear strain , as follows:

###### 4.4.3. Stress-Strain Relationship

In compliance with the classical elastic-plastic theory, we can obtain the expression for the elastic-plastic matrix. The constitutive stress-strain equation can be written in an incremental form as follows:

By substituting equations (27) and (35) and (36) into equation (37), the matrix elements *D*_{pp}, *D*_{pq}, *D*_{qp}, and *D*_{qq} in (37) are expressed as follows:

The modified Cam-Clay model for the matrix elements includes the critical state stress ratio *M*, Poisson’s ratio , compression coefficient *λ*, and rebound coefficient *k*. In comparison to the modified Cambridge model, the constitutive model of gravelly sand adds the parameters *s* and *t*, which are related to the initial relative density *D*_{r}. Thus, the constitutive model of gravelly sand can better reflect its stress-strain and the volumetric deformation characteristics.

#### 5. Parameter Determination and Preliminary Verification of Model

##### 5.1. Model Parameter Determination

A triaxial shear test was conducted on gravelly sand at initial relative densities of 0.15, 0.35, 0.55, and 0.75 and confining pressures of 50, 100, 200, and 400 kPa, respectively. Then, according to the results of the triaxial shear test, the parameters of the model were determined, and the model was preliminarily verified.

There were six parameters in the model: the critical state ratio *M* was related to the dilatancy equation; the compression coefficient *λ* and rebound coefficient *k* were related to the elastic volume strain and the hardening law; Poisson’s ratio was related to the elastic strain increment [37]; and the parameters *s* and *t* were related to the initial relative density. Among them, the parameter *M* denotes the slope of the critical state line, and the method for determining this parameter is described in Section 5.1.1. The parameters *λ* and *k* denote the slopes of the compression and rebound curves in the loading and unloading tests, respectively, and they are determined via the method described in Section 5.1.2. The parameters *s* and *t* were obtained by the test results obtained from the conventional triaxial test, using the method described in Section 5.1.3. The sensitivity analyses of these parameters are detailed in Section 5.3, where we examined how the changes in these variables influenced the stress-strain and volumetric curves of gravelly sand. The values of the above parameters under different test conditions based on the triaxial shear test results for gravelly sand are listed in Table 3.

###### 5.1.1. Critical State Stress Ratio *M*

In the axial shear tests, the values of *p* and *q* for different confining pressure failure points (maximum principal stress difference) were plotted in the *p*-*q* coordinate system. The slope of the line is the stress ratio of the critical state *M* and can be obtained via linear regression [38]. Gravelly sand with initial relative densities *D*_{r} of 0.15, 0.35, 0.55, and 0.75 was tested, and the *q* and *p* values corresponding to confining pressures of 50, 100, 200, and 400 kPa were plotted in the *q*-*p* plane. The slopes of *q*-*p* curves from the coordinate systems are shown in Figure 11. The *M* values of gravelly sand with initial relative densities of 0.15, 0.35, 0.55, and 0.75 are 1.34, 1.37, 1.4, and 1.44, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

###### 5.1.2. Compression Coefficient *λ* and Rebound Coefficient *k*

The compression coefficient *λ* and rebound coefficient *k* of gravelly sand can be determined via loading and unloading tests, respectively [39]. The *p* of the loading curve and the corresponding void ratio *e* were plotted on the *e*-ln *p* plane, and the slope of the straight line, which was obtained through linear regression, denoted the compression coefficient *λ*. Similarly, the *p* of the unloading curve and the corresponding void ratio *e* were plotted on the e-ln *p* plane, and the slope of this straight line, which was also obtained via linear regression, denoted the rebound coefficient *k*. In this study, loading and unloading tests were performed on gravelly sand with initial relative densities *D*_{r} of 0.15, 0.35, 0.55, and 0.75, and the test results were plotted in the *e*-ln *p* plane. As shown in Figure 12, when the initial relative densities *D*_{r} of gravelly sand were 0.15, 0.35, 0.55, and 0.75, the values of the compression coefficient *λ* were 0.058, 0.05, 0.04, and 0.025, and those of the rebound coefficient *k* were 0.0075, 0.004, 0.0035, and 0.003, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

###### 5.1.3. Material Parameters *s* and *t*

The material parameters *s* and *t* were related to the initial relative density, which could be obtained by fitting the results of the triaxial shear test for gravelly sand. The following empirical formula was devised and explains how the parameters *s* and *t* depend on the initial relative density *D*_{r.}:

It can be concluded from equation (42) that when the initial relative densities of gravelly sand were 0.15, 0.35, 0.55, and 0.75, the values of parameter *s* were 1.01, 1.4, 2, and 3 and those of *t* were 0.4, 0.75, 1.5, and 2.3, respectively.

##### 5.2. Preliminary Verification

To compare the fitting effect of the proposed model, the simulation results of the modified Cambridge model and the proposed constitutive model of gravelly sand were compared with the test results for initial relative densities of 0.15, 0.35, 0.55, and 0.75 and confining pressures of 50, 100, 200, and 400 kPa, respectively. The comparison results are shown in Figure 13. To facilitate observation and comparison, we defined *M*1 in Figure 13 as the calculation result of the modified Cambridge model and *M*2 as the calculation result of the proposed constitutive model for gravelly sand. As shown in Figure 13, the proposed constitutive model of gravelly sand could simulate the stress-strain and the volume deformation characteristics of gravelly sand under the initial relative densities *D*_{r} of 0.15, 0.35, 0.55, and 0.75; the accuracy was found to be better than that of the modified Cambridge model.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

**(m)**

**(n)**

**(o)**

**(p)**

##### 5.3. Sensitivity Analyses of Model Parameters

###### 5.3.1. Effects of Parameter *M*

The critical state stress ratio *M* was determined by the mechanical properties of the material and test results of the triaxial shear experiment. Figure 14 illustrates the relationship between the initial relative densities and the critical state stress ratio *M*. As shown in the figure, *M* increases with the initial relative densities *D*_{r} of gravelly sand. Through further sorting, the relationship between the critical state stress ratio *M* and the initial relative density *D*_{r} was found to be

Taking the test results for the situation where the initial relative density and confining pressure were 0.15 and 100 kPa, respectively, as an example, the stress-strain relationship and volumetric strain-axial strain relationship of gravelly sand were simulated using different model parameters *M*. The results are shown in Figure 15. With an increase in *M*, the hardening degree of the stress-strain curve and shrinkage degree of the volumetric strain-axial strain curves increased.

**(a)**

**(b)**

###### 5.3.2. Effects of Parameters *λ* and k

The compression coefficient *λ* and rebound coefficient *k* have definite physical meaning and can be determined via loading and unloading tests, respectively. Figure 16 shows how the initial pore ratio *e*_{0} (*e*_{0} = pore/ sand), compression coefficient *λ*, and rebound coefficient *k* are related to the initial relative density *D*_{r}.

**(a)**

**(b)**

**(c)**

As shown in Figure 16(a), the initial pore ratio *e*_{0} decreased with an increase in the initial relative density *D*_{r}. The relationship between the initial pore ratio *e*_{0} and initial relative density *D*_{r} can be fitted as follows:

Figure 16(b) shows that the compression coefficient *λ* decreases with an increase in the initial relative density *D*_{r}. The fitting relationship between the compression coefficient *λ* and initial relative density *D*_{r} is

As shown in Figure 16(c), the rebound coefficient *k* decreased as the initial relative density increased. The fitting relationship between the initial relative density and *k* is

Taking the test results for the situation where the initial relative density and 100 confining pressure were 0.15 and 100 kPa, respectively, as an example, the stress-strain relationship and volumetric strain-axial strain relationship of gravelly sand were simulated using different model parameters *λ* and *k*. The results are illustrated in Figure 17. Figures 17(a) and 17(b) show that as *λ* increases for gravelly sand, the hardening degree of the stress-strain curve and the shrinkage degree of the volumetric strain-axial strain curve decreased. As shown in Figures 17(c) and 17(d), *k* had minor influence on the stress-strain curve. With an increase in *k*, the shrinkage degree of the volumetric strain-axial strain curves decreased for gravelly sand.

**(a)**

**(b)**

**(c)**

**(d)**

###### 5.3.3. Effects of Parameters *s* and *t*

The parameters *s* and *t* were found to be dependent on the initial relative density. The relationships between the initial relative density and parameters *s* and *t* are demonstrated in Figure 18. As shown, both *s* and *t* increased with the initial relative density.

**(a)**

**(b)**

Taking the test results for the situation where the initial relative density and confining pressure were 0.15 and 100 kPa, respectively, as an example, the stress-strain relationship and volumetric strain-axial strain relationship of gravelly sand were simulated using different model parameters *s* and *t*. The results are illustrated in Figure 19. Figures 19(a) and 19(b) show that *s* significantly influenced the stress-strain and volumetric strain-axial strain curves of gravelly sand. With an increase in parameter *s*, the hardening degree of stress-strain curves of gravelly sand decreased and the initial tangent modulus increases. With an increase in parameter *s*, the shrinkage degree of the volumetric strain-axial strain curves decreased and gradually tended to dilatancy. As shown in Figures 19(c) and 19(d), with an increase in parameter *t*, the hardening degree of the stress-strain curves and the shrinkage degree of the volumetric strain-axial strain curves for gravelly sand increased.

**(a)**

**(b)**

**(c)**

**(d)**

#### 6. Conclusions

(1)The effects of initial relative density on the mechanical properties of gravelly sand are not negligible. When initial relative densities *D*_{r} were 0.15, 0.35, 0.55, and 0.75, all the samples of gravelly sand corresponded to strain hardening and shear contraction. Additionally, under a certain confining pressure, the hardening degree and peak stress of gravelly sand samples increased with an increase in initial relative density, and the volume shrinkage in soil sample decreased with the increase of initial relative density. In a certain initial relative density, the hardening degree and peak stress of samples increased with an increase in confining pressures. When initial relative densities *D*_{r} were 0.15 and 0.35, the volume deformation of soil sample decreased with the increasing confining pressure. Instead, when initial relative densities *D*_{r} were 0.55 and 0.75, the volume deformation of samples rose with the growth of confining pressure.(2)The shear contraction of gravelly sand is significantly affected by parameters *s* and *t.* At the same shear stress ratio *η*, the plastic dilatancy factor *d* increases and the volume shrinkage tendency of gravelly sand samples decreases with increasing parameter *s.* Similarly, the plastic dilatancy factor *d* of gravelly sand samples decreases and the volume shrinkage tendency of gravelly sand samples increases with the incremental parameter *t*(3)A law exists between initial relative density and parameter. The critical state stress ratio *M* rose with the growth of initial relative density, while the compressibility coefficient *λ* and rebound coefficient *k* decreased with initial relative density increasing(4)A higher-order dilatancy equation was proposed based on the concept of a super yield surface. Built on the associated flow laws, a constitutive model describing the mechanical properties for gravelly sand was also established. The proposed constitutive model of gravelly sand could simulate the stress-strain and volume deformation characteristics of gravelly sand under the initial relative densities *D*_{r} of 0.15, 0.35, 0.55, and 0.75; the accuracy was found to be better than that of the modified Cambridge model

#### Data Availability

Since the experiment was completed with the support of Sichuan Agricultural University, the data used to support the results of this study can be obtained from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Authors’ Contributions

Dongjie Zhang and Fei Luo contributed equally.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (41701063 and 41672304).