# Experimental Testing and Numerical Modelling of Mechanical Behaviors of Silty Clay under Freezing-Thawing Cycles

Authors: Tianfei Hu; Dawei Liu; Hui Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Seasonal freezing-thawing cycle is one of the most common physical weathering processes in cold regions, which can significantly affect the mechanical behaviors of soil. In this paper, a series of freezing-thawing (F-T) cycle and consolidated drained triaxial tests have been carried out on silty clay samples collected from Tibetan Plateau. To do so, a modified numerical model was developed taking into F-T effect. Test results showed that the stress-strain curves of original soil specimens presented strain hardening characteristics, accompanied with shear shrinkage. In F-T experienced specimens, volumetric strain in triaxial loading stage was gradually increased, while failure strength was decreased. Elliptic and parabolic functions were selected in numerical modelling to describe volume and shear yield surfaces on ap-q plane, respectively. Moreover, a double-yield surface constitutive model was developed to describe relationships among deviatoric stress, axial strain, and volumetric strain. Furthermore, equations for model parameters with the number of F-T cycles as variable were derived based on the triaxial test results which were then substituted into the established model to take into account the effects of F-T cycles. Finally, numerical results were validated with experimental findings.


---

## Body

## 1. Introduction

In cold regions, ground bases, and soil foundations such as subgrades, generally suffer seasonal and diurnal freezing-thawing (F-T) cycles due to atmospheric temperature changes. F-T process causes soil structure degradation, strength reduction, settlement, etc. Numerical constitutive models can quantitatively predict the mechanical behaviors of soil [1]. Therefore, numerical modelling of soil mechanical behaviors is one of the main geotechnical engineering research issues in cold regions.Currently, soil constitutive models mainly include nonlinear elastic and elastoplastic types. Nonlinear elastic models include DuncanE-v, E-B, and K-G models, none of which can describe soil dilatancy and the influences of both compress and shear loadings [2]. Elastoplastic models mainly refer to the Cam-clay model which was first proposed by Roscoe and Schofield [3] and a series of modified single-yield surface models developed on the basis of the original one. Among them, Asaoka et al. [4] established a modified model capable of reflecting historical stress path and soil structure characteristics, by embedding super-consolidation and structural parameters into the Cam-clay model yield function to form upper and lower loading yield surfaces in the p-q plain. Based on state-related dilatancy theory, Li and Dafalias [5] suggested several parameters for the modification of the Cam-clay model to make it capable of describing shear contraction and dilatancy characteristics of soil. Moreover, Yin and Graham [6] developed a viscoelastic-plastic constitutive model considering time effect based on equivalent time concept. Furthermore, Huang et al. [7] established the Tsinghua model, and Li [8] improved it. They also experimentally determined the direction of plastic strain increment under each stress state and found yield surface according to the associated flow law which is only appropriate for soil strain hardening process. Yao [9] proposed UH models suitable for both super- and normal consolidated soils, based on the Cam-clay model and lower loading surface concept. To sum up, the single-yield surface model is a “hat” model. Therefore, it can only reflect soil shear shrinkage. However, there may still be yielding phenomena within yield surfaces in stress space.To address the shortcomings of single-yield surface models, two-yield and multiple-yield surface models were developed. Lade [10] first elaborated two-yield surface concept and further defined volumetric and shear yield surfaces as failure and plastic expansion yield surfaces, respectively; the latter one being a bullet-shaped cone. Shen [11] combined the advantages of Duncan–Chang and Cam-clay models and proposed a two-yield surface model called Nanshui model, which could simulate the stress path characteristics of decreasing confining or average pressures. Zheng and Chen [12] developed multiple-yield surfaces based on the vector principle and introduced the concepts of volume yield surface, shear yield surface, and shear yield surface of θσ direction. Yin et al. [13, 14] proposed yield and hardening functions relating to compress and shear loadings, respectively, and then established another two-yield surface model using correlated flow law. Huang et al. [15] introduced a shape parameter into the Cam-clay model and developed a two-yield surface model which was able to reflect soil strain softening. Several other yield surface models and boundary theory-based models have also been developed [16].At present, one of the popular research fields is to improve models to make them applicable for different soil properties, environments, and boundary conditions based on the above models [17]. Particularly, for soil constructions in cold regions, F-T process of soil causes drastic changes in its physical and mechanical characteristics. Many forms of visco-elastoplastic and elastoplastic models [18, 19] have been developed based on the mechanical behaviors of frozen soil. Meanwhile, the thawing of frozen soil is a major source of engineering problems. Regarding soil thawing, Liu et al. [20] proposed two new factors based on the Duncan model, namely, deviatoric stress reduction ratio and the number of F-T cycles, to predict the stress-strain relationship of soil under F-T cycles. Cui et al. [21] fitted volume and shear hardening functions with logarithmic and exponential functions, respectively, and evaluated F-T effect by describing the variation law of hardening function coefficient. Chang et al. [22, 23] established a two-yield surface model considering F-T effect by fitting the equations of hardening function coefficients as a cubic polynomial with the number of F-T cycles. The above studies showed that soil F-T effect could be simply determined by selecting reasonable yield surfaces and equations of functional coefficients based on typical two-yield surface constitutive models. However, the applicability of this model considering F-T effect remains to be further analyzed.In this paper, a series of F-T cycle and consolidated drained triaxial tests were carried out on the silty clay samples collected from Tibetan Plateau, where soil ground suffers continuous F-T cycles. A two-yield surface constitutive model was first established based on the conventional Nashui Model. Then, model parameter equations by considering the number of F-T cycles as a variable were derived according to triaxial test results. Finally, a modified model was proposed and validated to predict the deviatoric stress, axial strain, and volumetric strain of soil under F-T cycles.


## 2. Soil Specimens and Test Scheme

### 2.1. Soil Materials

The soil samples employed in the experiments were collected from the newly constructed Gonghe-Yushu Expressway, which locates in Tibetan Plateau, China, and is the first expressway in the world to be constructed in a permafrost region, as shown in Figure1. The collected soil samples were used as a subgrade filler for expressway and showed obvious F-T effects on both strength and deformation behaviors.Figure 1
Location of soil collection site along Gonghe-Yushu Expressway.The grain composition of the soil sample is shown in Figure2. Soil was turned over, dried, and passed through a 2 mm sieve before testing, according to the preparation requirements of soil samples as advised in reference [24]. The soil sample had maximum dry density of 1.828 g/cm3, optimal moisture content of 14.8%, relative density of 2.64, liquid limit of 28.0%, and plasticity index of 10.3. Based on the above values, the soil sample was specified as silty clay.Figure 2
Grain size distribution of the tested soil.


### 2.2. Preparation of Specimens and Test Scheme

According to subgrade compaction requirements, the specimen preparation standards for moisture content and compaction degree were considered to be 14.8% and 95%, respectively. Figure3 shows the details of the simulation of subgrade fillers in real operational conditions. Cylinder specimens with diameter 39.1 mm and height 80 mm were prepared through the layered compaction method. First, F-T cycle tests were performed in an incubator with alternate constant temperatures of −5°C and 20°C with a lasting time of 12 hours at each temperature to ensure the complete freezing and thawing of specimens, as shown in Figure 3(a). F-T-experienced specimens are shown in Figure 3(b), which were covered with plastic layers to prevent water evaporation or absorption. Then, in order to test the pore water pressure accurately to further study the mechanical and deformation characteristics of tested soils and to obtain the potential most unfavorable results, the thawed specimens were taken out of the incubator for vacuum saturation in a sealed vessel, as shown in Figure 3(c). Finally, F-T-experienced saturated specimens underwent consolidated drained triaxial tests on a triaxial loading apparatus made by GDS Instrument Company, UK, as shown in Figure 3(d). Specifically, the F-T-experienced sample was further saturated through a back-pressure saturation process, which was conducted by a stepwise (20 kPa) method and continued until the ratio of pore water pressure increment to confining pressure increment was larger than 98%. The drain valve of the GDS apparatus was then opened, and the consolidation process was carried out under a predetermined confining pressure, until the pore water pressure dissipated by more than 95%. At last, the loading test started with a constant axial shearing rate of 0.01%/min, and it was continued until the ultimate axial strain of each sample reached 20%. All the above testing methods are referred to the provisions of reference [24].Figure 3
Details of the experimental process: (a) the incubator used in F-T cycle tests; (b) the specimens subjected to F-T cycles; (c) the sealed vessel for vacuum saturation; (d) GDS triaxial apparatus.
(a)
(b)
(c)
(d)In this study, the number of F-T cycles and confining pressure were selected as test variables in F-T cycle and consolidated drained triaxial tests, respectively. F-T cycle numbers were selected as 0, 1, 3, 6, 9, and 12 times, and confining pressures for both consolidating and triaxial loading process were set at 100, 200, 300, and 400 kPa. The axial strain rate was set at 0.01%/min, and the test was terminated when the axial strain reached 20%. A total of 24 specimens were tested.


## 2.1. Soil Materials

The soil samples employed in the experiments were collected from the newly constructed Gonghe-Yushu Expressway, which locates in Tibetan Plateau, China, and is the first expressway in the world to be constructed in a permafrost region, as shown in Figure1. The collected soil samples were used as a subgrade filler for expressway and showed obvious F-T effects on both strength and deformation behaviors.Figure 1
Location of soil collection site along Gonghe-Yushu Expressway.The grain composition of the soil sample is shown in Figure2. Soil was turned over, dried, and passed through a 2 mm sieve before testing, according to the preparation requirements of soil samples as advised in reference [24]. The soil sample had maximum dry density of 1.828 g/cm3, optimal moisture content of 14.8%, relative density of 2.64, liquid limit of 28.0%, and plasticity index of 10.3. Based on the above values, the soil sample was specified as silty clay.Figure 2
Grain size distribution of the tested soil.


## 2.2. Preparation of Specimens and Test Scheme

According to subgrade compaction requirements, the specimen preparation standards for moisture content and compaction degree were considered to be 14.8% and 95%, respectively. Figure3 shows the details of the simulation of subgrade fillers in real operational conditions. Cylinder specimens with diameter 39.1 mm and height 80 mm were prepared through the layered compaction method. First, F-T cycle tests were performed in an incubator with alternate constant temperatures of −5°C and 20°C with a lasting time of 12 hours at each temperature to ensure the complete freezing and thawing of specimens, as shown in Figure 3(a). F-T-experienced specimens are shown in Figure 3(b), which were covered with plastic layers to prevent water evaporation or absorption. Then, in order to test the pore water pressure accurately to further study the mechanical and deformation characteristics of tested soils and to obtain the potential most unfavorable results, the thawed specimens were taken out of the incubator for vacuum saturation in a sealed vessel, as shown in Figure 3(c). Finally, F-T-experienced saturated specimens underwent consolidated drained triaxial tests on a triaxial loading apparatus made by GDS Instrument Company, UK, as shown in Figure 3(d). Specifically, the F-T-experienced sample was further saturated through a back-pressure saturation process, which was conducted by a stepwise (20 kPa) method and continued until the ratio of pore water pressure increment to confining pressure increment was larger than 98%. The drain valve of the GDS apparatus was then opened, and the consolidation process was carried out under a predetermined confining pressure, until the pore water pressure dissipated by more than 95%. At last, the loading test started with a constant axial shearing rate of 0.01%/min, and it was continued until the ultimate axial strain of each sample reached 20%. All the above testing methods are referred to the provisions of reference [24].Figure 3
Details of the experimental process: (a) the incubator used in F-T cycle tests; (b) the specimens subjected to F-T cycles; (c) the sealed vessel for vacuum saturation; (d) GDS triaxial apparatus.
(a)
(b)
(c)
(d)In this study, the number of F-T cycles and confining pressure were selected as test variables in F-T cycle and consolidated drained triaxial tests, respectively. F-T cycle numbers were selected as 0, 1, 3, 6, 9, and 12 times, and confining pressures for both consolidating and triaxial loading process were set at 100, 200, 300, and 400 kPa. The axial strain rate was set at 0.01%/min, and the test was terminated when the axial strain reached 20%. A total of 24 specimens were tested.


## 3. Test Results and Analysis

### 3.1. Volume Change of Specimens Subjected to F-T Cycles

The volume change of specimens with the number of F-T cycles is illustrated in Figure4. It can be seen that the specimen volume was increased by increasing the number of F-T cycles. The ultimate volume increasing rate can be 3.63% at the F-T cycle numbers of 12. It is because the volume of liquid water in soil pores increases by 9% after freezing [25]. Pore water freezing process destroys the original soil structure and the connection of soil particles which cannot be restored during subsequent thawing. As a result, continuous degradation on soil mechanical behaviors was observed in specimens.Figure 4
Volume change of specimens with different F-T cycle numbers.


### 3.2. Deviatoric Stress-Axial Strain Characteristics

Deviatoric stress-axial strain curves under different test conditions are shown in Figure5. It can be seen that all stress-strain curves showed strain hardening and the deviatoric stress was gradually increased with the increase of axial strain. In addition, at low confining pressures, stress-strain curves showed certain weak strain softening characteristics. Meanwhile, the spatial position of the stress-strain curve was gradually decreased with the increase of F-T cycle numbers.Figure 5
Deviatoric stress-axial strain curves of specimens for different F-T cycle numbers at confining pressure of (a) 100 kPa; (b) 200 kPa; (c) 300 kPa; (d) 400 kPa.
(a)
(b)
(c)
(d)From the perspective of practice, it is feasible to use the peak value of the deviatoric stress as the failure criterion. For simplicity, the failure strength of each specimen was obtained based on its deviatoric stress value corresponding to axial strain of 15% as advised in reference [24]. The change rate of failure strength is displayed in Figure 6. It can be seen in the figure that failure strength was decreased with the increase of the number of F-T cycles. Meanwhile, the change rate of failure strength was decreased by increasing confining pressure, which indicated that confining pressure could weaken F-T effect.Figure 6
Change rate of failure strength of specimens at different F-T cycle numbers and confining pressures.


### 3.3. Volumetric Strain-Axial Strain Characteristics

Consolidated drained triaxial test results for pore water drainage volume after different F-T cycle numbers are shown in Table1. It can be seen that drainage volume was increased by increasing both F-T cycle number and confining pressure. It can be seen from Figure 6 and Table 1 that confining pressure at both triaxial loading stage and consolidation stage is beneficial to reduce the structural deterioration of the specimen in the F-T cycling test. Meanwhile, drainage volume increase due to F-T cycles revealed that the unfavorable degradation of soil structure was increased with the number of F-T cycles. Nonetheless, soil structure was recovered to some extent during consolidating process, and recovery degree was increased by confining pressure.Table 1
Drainage volume of specimens during consolidating process under different F-T cycle numbers and confining pressures (cm3).
Confining pressure (kPa)
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
100
1.29
1.29
1.32
1.31
1.36
1.45
200
1.99
2.29
2.37
2.32
2.31
2.40
300
2.49
2.50
2.71
2.71
2.95
3.03
400
3.03
3.03
3.30
3.32
3.41
3.60The volumetric strain-axial strain relationships of the loading stage under different test conditions are shown in Figure7. It can be seen that specimens basically showed shear shrinkage property, and volumetric strain was gradually increased by increasing axial strain under similar conditions. Except for non-F-T-experienced specimens, they showed a dilatation behavior in the early loading stage under the confining pressure of 100 kPa. Meanwhile, with the increase of F-T cycle numbers and confining pressures, volumetric strain corresponding to similar axial strain was also increased. However, at a confining pressure of 400 kPa, the volumetric strain corresponding to the same axial strain showed an opposite decreasing trend compared with that at 300 kPa, which was because of more significant consolidation drainage under the former test condition.Figure 7
The volumetric strain-axial strain curves of specimens with different F-T cycle numbers at confining pressure of (a) 100 kPa; (b) 200 kPa; (c) 300 kPa; (d) 400 kPa.
(a)
(b)
(c)
(d)


## 3.1. Volume Change of Specimens Subjected to F-T Cycles

The volume change of specimens with the number of F-T cycles is illustrated in Figure4. It can be seen that the specimen volume was increased by increasing the number of F-T cycles. The ultimate volume increasing rate can be 3.63% at the F-T cycle numbers of 12. It is because the volume of liquid water in soil pores increases by 9% after freezing [25]. Pore water freezing process destroys the original soil structure and the connection of soil particles which cannot be restored during subsequent thawing. As a result, continuous degradation on soil mechanical behaviors was observed in specimens.Figure 4
Volume change of specimens with different F-T cycle numbers.


## 3.2. Deviatoric Stress-Axial Strain Characteristics

Deviatoric stress-axial strain curves under different test conditions are shown in Figure5. It can be seen that all stress-strain curves showed strain hardening and the deviatoric stress was gradually increased with the increase of axial strain. In addition, at low confining pressures, stress-strain curves showed certain weak strain softening characteristics. Meanwhile, the spatial position of the stress-strain curve was gradually decreased with the increase of F-T cycle numbers.Figure 5
Deviatoric stress-axial strain curves of specimens for different F-T cycle numbers at confining pressure of (a) 100 kPa; (b) 200 kPa; (c) 300 kPa; (d) 400 kPa.
(a)
(b)
(c)
(d)From the perspective of practice, it is feasible to use the peak value of the deviatoric stress as the failure criterion. For simplicity, the failure strength of each specimen was obtained based on its deviatoric stress value corresponding to axial strain of 15% as advised in reference [24]. The change rate of failure strength is displayed in Figure 6. It can be seen in the figure that failure strength was decreased with the increase of the number of F-T cycles. Meanwhile, the change rate of failure strength was decreased by increasing confining pressure, which indicated that confining pressure could weaken F-T effect.Figure 6
Change rate of failure strength of specimens at different F-T cycle numbers and confining pressures.


## 3.3. Volumetric Strain-Axial Strain Characteristics

Consolidated drained triaxial test results for pore water drainage volume after different F-T cycle numbers are shown in Table1. It can be seen that drainage volume was increased by increasing both F-T cycle number and confining pressure. It can be seen from Figure 6 and Table 1 that confining pressure at both triaxial loading stage and consolidation stage is beneficial to reduce the structural deterioration of the specimen in the F-T cycling test. Meanwhile, drainage volume increase due to F-T cycles revealed that the unfavorable degradation of soil structure was increased with the number of F-T cycles. Nonetheless, soil structure was recovered to some extent during consolidating process, and recovery degree was increased by confining pressure.Table 1
Drainage volume of specimens during consolidating process under different F-T cycle numbers and confining pressures (cm3).
Confining pressure (kPa)
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
100
1.29
1.29
1.32
1.31
1.36
1.45
200
1.99
2.29
2.37
2.32
2.31
2.40
300
2.49
2.50
2.71
2.71
2.95
3.03
400
3.03
3.03
3.30
3.32
3.41
3.60The volumetric strain-axial strain relationships of the loading stage under different test conditions are shown in Figure7. It can be seen that specimens basically showed shear shrinkage property, and volumetric strain was gradually increased by increasing axial strain under similar conditions. Except for non-F-T-experienced specimens, they showed a dilatation behavior in the early loading stage under the confining pressure of 100 kPa. Meanwhile, with the increase of F-T cycle numbers and confining pressures, volumetric strain corresponding to similar axial strain was also increased. However, at a confining pressure of 400 kPa, the volumetric strain corresponding to the same axial strain showed an opposite decreasing trend compared with that at 300 kPa, which was because of more significant consolidation drainage under the former test condition.Figure 7
The volumetric strain-axial strain curves of specimens with different F-T cycle numbers at confining pressure of (a) 100 kPa; (b) 200 kPa; (c) 300 kPa; (d) 400 kPa.
(a)
(b)
(c)
(d)


## 4. Two-Yield Surface Constitutive Model

In elliptic-parabolic two-yield surface constitutive model for soil [13, 14], yield surface is regarded as the boundary of elastic deformation regions which are composed of elliptic and power functions. Soil deformation consists of elastic dεe and plastic dεp deformations. dεp is composed of dεvp and dεsp, where dεvp is related to compression and dεsp is related to expansion. Both types of plastic deformation are assumed to be associated with flow law. In Figure 8, f1 and f2 represent the yield locus and N is the current point. Hence, the p-q plane is divided into four regions of A0, A1, A2, and A3, where A0 is only characterized by elastic strain, A1 is related to plastic compression yield, A2 is related to plastic expansion yield, and A3 is related to both compression and expansion yields simultaneously. At point N, plastic increment direction was the vector sum of different plastic strain increments. As a result, the composition of total strain dε was(1)dε=dεe+dεvp+dεsp.Figure 8
Elliptic-parabolic two-yield surface model [13, 14].


### 4.1. Yield Associated with Compression and Expansion

The yield surface corresponding to compression adopts elliptic shape in the modified Cam-clay model [2] and improved yield equation f1 is still in the form of the elliptic function. The main equations and parameters were expressed as follows [13, 14]:(2)f1σij,Ha=p+q2M12p+pr−Ha1εvp=0,where M1 presents the elliptic shape change, p is the mean normal stress, q is the generalized shear stress, εvp is the hardening parameter of elliptic yield surface, pr is the absolute value of the intercept of failure line on the p-axis in the p-q plane, and Ha1 is the hardening function of compression yield surface.The calculation method forpr was stated as(3)pr=ccotφ.Under isostatic triaxial testing conditions, the relationship betweenp0 and εvp could be expressed in a hyperbolic form. Therefore, the hardening function corresponding to the first yield surface was written as [13, 14](4)Ha1εvp=p0=hεvp1−tεvppa,where p0 is the horizontal coordinate of the intersection point of yield locus f1 and axis p and h and t are the dimensionless parameters to describe the relationship between confining pressure and plastic volumetric strain.Yield equationf2 corresponding to shear expansion was a parabolic function expressed as [13, 14](5)f2σij,Ha=aqGqM2p+pr−q−Ha2εsp=0,where a represents the strain ratio of yield surface f2 to total plastic strain; therefore, high values of a correspond to dilation. G is the elastic shear modulus, εsp is the hardening parameter of parabolic yield surface, and Ha2 is the hardening function of expansion yield surface.G was calculated as(6)G=E21+υ,where for the elastic modulus, we had E = 2.0Ei. Generally, the proportion of elastic strain in the total strain is small; therefore, vis always set at 0.3. The initial tangent modulus Ei is related to confining pressure, and its empirical equation was derived as(7)Ei=kpaσ3pan,where pa is the standard atmospheric pressure, which is 101.3 kPa, and k and n are the intercept and slope of the line lg(Ei/pa)∼lg(σ3/pa).The hardening function corresponding to the second yield surface was obtained to be [13, 14](8)Ha2εsp=εsp=∫dεsp=∫23⋅dε1p−dε2p2+dε2p−dε3p2+dε3p−dε1p2.


### 4.2. Constitutive Equations in Increment Form

The three strain components of equation (1) have been determined in this chapter, where elastic strain component dεe included volumetric elastic strain increment dεve and generalized shear elastic strain increment dεse. Calculation equations were as follows:(9)dεve=dpK=31−2υdp2G1+υ,dεse=dq3G,where K is the bulk modulus.By differentiating both sides of equation (2), we found(10)∂f1∂pdp+∂f1∂qdq=∂Ha1εvp1∂εvp1dεvp1.The plastic compression component of the first yield surface was consisted of two parts: volume change incrementdεvp1 and shear strain increment dεsp1. According to equations (4) and (10), the following equation was obtained for dεvp1:(11)dεvp1=hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2dp+2qM12p+prdq.Assuming both yield surfaces to be applicable to associated flow law, the following equation was obtained:(12)dεsp1dεvp1=∂f1/∂q∂f1/∂p=2qp+prM12p+pr2−q2.The calculation equation ofdεsp1 was obtained by combining equations (11) and (12).The plastic compression components of the second yield surface included volumetric elastic straindεvp2 and shear strain dεsp2 increments. Using the same derivation method as the first yield surface, the following equation was obtained for dεsp2:(13)dεsp2=aG⋅qM2p+pr−qdq+a2G⋅q1/2p+prdq−q3/2dpM2p+pr−q3/2,(14)dεvp2dεsp2=−M2q3M2p+pr−2q.The calculation equation ofdεvp2 was obtained by combining equations (13) and (14).According to typical elastic-plastic theory,dεv and dεs were defined as(15)dεvdεs=ABCDdpdq,where A, B, C, and D were expressed as follows:(16)A=31−2υ2G1+υ+hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2+−M2q3M2p+pr−2q×a2G⋅−q3/2M2p+pr−q3/2,(17)B=hpahpa+tp+q2/M12p+pr2×2qM12p+pr+−M2q3M2p+pr−2q×aG⋅qM2p+pr−q+a2G⋅q1/2p+prM2p+pr−q3/2,(18)C=2qp+prM12p+pr2−q2×hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2+a2G⋅−q3/2M2p+pr−q3/2,(19)D=13G+2qp+prM12p+pr2−q2×hpahpa+tp+q2/M12p+pr2×2qM12p+pr+aG⋅qM2p+pr−q+a2G⋅q1/2p+prM2p+pr−q3/2.


### 4.3. Equations of Triaxial Tests

In conventional consolidated drained triaxial tests, stress state is generally inA3 region, as illustrated in Figure 8. The relationships between p and q and increments dp and dq are given by equations (20) and (21), respectively:(20)p=q3+σ3,(21)dp=dq3.In strain-controlled triaxial tests, the relationship of axialdεa and volumetric dεv strain increments were(22)dεa=dεs+13dεv.By substituting equation (21) into equation (15), the volumetric strain increment was found to be(23)dεv=13A+Bdq.By substituting equations (22) and (23) into equation (15), axial strain increment was derived as(24)dεa=19A+13B+13C+Ddq.Finally, the incremental constitutive equation of the relationships among deviatoric stress, axial strain, and volumetric strain of soil were obtained by combining equations (16)∼(19), (23), and (24). There are 9 parameters in the above numerical model including c, φ, M1, h, t, k, n, M2, and a.


## 4.1. Yield Associated with Compression and Expansion

The yield surface corresponding to compression adopts elliptic shape in the modified Cam-clay model [2] and improved yield equation f1 is still in the form of the elliptic function. The main equations and parameters were expressed as follows [13, 14]:(2)f1σij,Ha=p+q2M12p+pr−Ha1εvp=0,where M1 presents the elliptic shape change, p is the mean normal stress, q is the generalized shear stress, εvp is the hardening parameter of elliptic yield surface, pr is the absolute value of the intercept of failure line on the p-axis in the p-q plane, and Ha1 is the hardening function of compression yield surface.The calculation method forpr was stated as(3)pr=ccotφ.Under isostatic triaxial testing conditions, the relationship betweenp0 and εvp could be expressed in a hyperbolic form. Therefore, the hardening function corresponding to the first yield surface was written as [13, 14](4)Ha1εvp=p0=hεvp1−tεvppa,where p0 is the horizontal coordinate of the intersection point of yield locus f1 and axis p and h and t are the dimensionless parameters to describe the relationship between confining pressure and plastic volumetric strain.Yield equationf2 corresponding to shear expansion was a parabolic function expressed as [13, 14](5)f2σij,Ha=aqGqM2p+pr−q−Ha2εsp=0,where a represents the strain ratio of yield surface f2 to total plastic strain; therefore, high values of a correspond to dilation. G is the elastic shear modulus, εsp is the hardening parameter of parabolic yield surface, and Ha2 is the hardening function of expansion yield surface.G was calculated as(6)G=E21+υ,where for the elastic modulus, we had E = 2.0Ei. Generally, the proportion of elastic strain in the total strain is small; therefore, vis always set at 0.3. The initial tangent modulus Ei is related to confining pressure, and its empirical equation was derived as(7)Ei=kpaσ3pan,where pa is the standard atmospheric pressure, which is 101.3 kPa, and k and n are the intercept and slope of the line lg(Ei/pa)∼lg(σ3/pa).The hardening function corresponding to the second yield surface was obtained to be [13, 14](8)Ha2εsp=εsp=∫dεsp=∫23⋅dε1p−dε2p2+dε2p−dε3p2+dε3p−dε1p2.


## 4.2. Constitutive Equations in Increment Form

The three strain components of equation (1) have been determined in this chapter, where elastic strain component dεe included volumetric elastic strain increment dεve and generalized shear elastic strain increment dεse. Calculation equations were as follows:(9)dεve=dpK=31−2υdp2G1+υ,dεse=dq3G,where K is the bulk modulus.By differentiating both sides of equation (2), we found(10)∂f1∂pdp+∂f1∂qdq=∂Ha1εvp1∂εvp1dεvp1.The plastic compression component of the first yield surface was consisted of two parts: volume change incrementdεvp1 and shear strain increment dεsp1. According to equations (4) and (10), the following equation was obtained for dεvp1:(11)dεvp1=hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2dp+2qM12p+prdq.Assuming both yield surfaces to be applicable to associated flow law, the following equation was obtained:(12)dεsp1dεvp1=∂f1/∂q∂f1/∂p=2qp+prM12p+pr2−q2.The calculation equation ofdεsp1 was obtained by combining equations (11) and (12).The plastic compression components of the second yield surface included volumetric elastic straindεvp2 and shear strain dεsp2 increments. Using the same derivation method as the first yield surface, the following equation was obtained for dεsp2:(13)dεsp2=aG⋅qM2p+pr−qdq+a2G⋅q1/2p+prdq−q3/2dpM2p+pr−q3/2,(14)dεvp2dεsp2=−M2q3M2p+pr−2q.The calculation equation ofdεvp2 was obtained by combining equations (13) and (14).According to typical elastic-plastic theory,dεv and dεs were defined as(15)dεvdεs=ABCDdpdq,where A, B, C, and D were expressed as follows:(16)A=31−2υ2G1+υ+hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2+−M2q3M2p+pr−2q×a2G⋅−q3/2M2p+pr−q3/2,(17)B=hpahpa+tp+q2/M12p+pr2×2qM12p+pr+−M2q3M2p+pr−2q×aG⋅qM2p+pr−q+a2G⋅q1/2p+prM2p+pr−q3/2,(18)C=2qp+prM12p+pr2−q2×hpahpa+tp+q2/M12p+pr2×1−q2M12p+pr2+a2G⋅−q3/2M2p+pr−q3/2,(19)D=13G+2qp+prM12p+pr2−q2×hpahpa+tp+q2/M12p+pr2×2qM12p+pr+aG⋅qM2p+pr−q+a2G⋅q1/2p+prM2p+pr−q3/2.


## 4.3. Equations of Triaxial Tests

In conventional consolidated drained triaxial tests, stress state is generally inA3 region, as illustrated in Figure 8. The relationships between p and q and increments dp and dq are given by equations (20) and (21), respectively:(20)p=q3+σ3,(21)dp=dq3.In strain-controlled triaxial tests, the relationship of axialdεa and volumetric dεv strain increments were(22)dεa=dεs+13dεv.By substituting equation (21) into equation (15), the volumetric strain increment was found to be(23)dεv=13A+Bdq.By substituting equations (22) and (23) into equation (15), axial strain increment was derived as(24)dεa=19A+13B+13C+Ddq.Finally, the incremental constitutive equation of the relationships among deviatoric stress, axial strain, and volumetric strain of soil were obtained by combining equations (16)∼(19), (23), and (24). There are 9 parameters in the above numerical model including c, φ, M1, h, t, k, n, M2, and a.


## 5. Functional Expressions for Model Parameters

According to the test results of consolidated drained triaxial tests, the model parameter variation rules for the number of F-T cyclesN and corresponding function expressions were determined, namely, c(N), φ(N), M1(N), h(N), t(N), k(N), n(N), M2(N), and a(N).


### 5.1. Determination ofc(N) and φ(N)

The internal friction angleφ and cohesion c were calculated according to the slope and intercept of the tangent line of Mohr’s stress circle as(25)φ=sin−1tanα,c=dcosφ,where α is the slope angle and d is the intercept of the envelope.The tangent lines of Mohr’s stress circles for different F-T cycle numbers are shown in Figure9, and shear strength parameters are summarized in Table 2. It was seen that both c and φ were decreased with the number of F-T cycles. The logistic function, as given by equation (26), was used to fit the regression relationships among φ, c, and number of F-T cycles, as shown in Figure 10:(26)y=x0−A11+N/A2A3+A1,where x0 is the shear strength parameters of non-F-T-experienced specimens and A1∼A3 are the fitting parameters.Figure 9
Variation of tangent lines of Mohr’s stress circles for different F-T cycle numbers.Table 2
Shear strength parameters for different F-T cycle numbers.
Shear strength
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
φ (°)
31.97
31.50
30.95
30.56
30.17
29.87
c (kPa)
21.05
19.85
19.45
19.39
19.30
19.24Figure 10
Variation of shear strength parameters for different F-T cycle numbers.


### 5.2. Determination ofM1(N)

The empirical equation ofM1(N) is derived as [13, 14](27)M1=1+0.25β2M,where β is the ratio of volumetric strain to axial strain at stress level 75% and β values vary with confining pressure; therefore, its average value was adopted. M is the slope of failure line as shown in Figure 8, which can be obtained by equation (28) combined with triaxial compression test results:(28)M=6sinφ3−sinφ.The values ofβ, M, and M1 are defined in Table 3. It can be seen that β was gradually increased with the number of F-T cycles, while M and M1 were decreased. The regression relationship between M1 and the number of F-T cycles was described by the logistic function, as shown in Figure 11.Table 3
The values ofβ, M, and M1 for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
β
0.460
0.503
0.534
0.546
0.574
0.616
M
1.286
1.265
1.241
1.224
1.207
1.194
M
1
1.340
1.328
1.312
1.296
1.286
1.284Figure 11
The variation ofM1 for different F-T cycle numbers.


### 5.3. Determination ofh(N) and t(N)

h(N) and t(N) were obtained through modifying equation (4) [13, 14]:(29)dεvp=hpadp0hpa+tp02.WhenBp=dp0/dεvp, equation (29) could be transformed as(30)Bppa=h+thp0pa.According to the previous experiences, plastic volumetric strainεvp accounts for about 65% of volumetric strain by increasing the stress level from 0 to 50%. Based on equation (2), p0 was obtained through p and q at the stress level of 50%. Then, Bp was calculated according to εvp and p0, specifically. The intercept and slope of the curve of (Bp/pa)0.5∼p0/pa were h0.5 and t/h0.5, respectively; therefore, h and t were also obtained.Curves of (Bp/pa)0.5∼p0/pa for different F-T cycle numbers are shown in Figure 12, and h and t values are given in Table 4. It was seen that h was gradually decreased with the number of F-T cycles while t was increased. The relationships among h, t, and the number of F-T cycles were fitted by the logistic function, as shown in Figure 13.Figure 12
Variation ofp0/pa∼(Bp/p0)0.5 with F-T cycle numbers.Table 4
Values ofh and t for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
h
884.23
603.59
486.60
385.73
302.80
197.91
t
6.538
6.756
7.548
7.755
7.970
7.778Figure 13
Variation ofh and t with F-T cycle numbers.


### 5.4. Determination ofk(N) and n(N)

The curve of lg (σ3/pa)∼lg (Ei/pa) based on triaxial test results is shown in Figure 14. It was seen that the spatial positions of fitting lines moved downwards by increasing the number of F-T cycles. The values of intercept k and slope n are given in Table 5. As can be seen, k was decreased with the number of F-T cycles while n was increased. The relationships among k, n, and the number of F-T cycles were described by the logistic function, as illustrated in Figure 15.Figure 14
Variation of lg (σ3/pa)∼lg (Ei/pa) with F-T cycle numbers.Table 5
The values ofk and n under different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
lgk
1.566
1.492
1.450
1.404
1.388
1.374
n
0.478
0.636
0.703
0.723
0.767
0.758
k
36.81
31.05
28.18
25.35
24.43
23.66
R
2
0.998
0.975
0.944
0.977
0.967
0.963Figure 15
Variations ofk and n with F-T cycle numbers.


### 5.5. Determination ofM2(N) and a(N)

M
2 and a were determined by transforming equations (5) and (8) [13, 14]:(31)p+prq=a2M2qεspG2+1M2.In equation (31), εsp was obtained using the following empirical equation:(32)εsp=0.3−0.1dεa,where d is the slope of the curve of εv∼εa at a stress levels of 75%∼95%. The above valuing method of d is actually an empirical one, and its reliability has been verified [13, 14]. Negative values of d meant that there was dilation. The values of d in each test curve were different, and therefore, its mean values were adopted.The intercept of the curve of (p + pr)/q∼(q/εspG)2 was 1/M2, and its slope was a2/M2; therefore, M2 and a were also obtained. The values of M2 and a for different F-T cycle numbers are summarized in Table 6. It was seen that M2 and a were decreased with the number of F-T cycles and the relationships among M2, a, and the number of F-T cycles were described by the logistic function, as shown in Figure 16.Table 6
Values ofM2 and a for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
M
2
1.486
1.466
1.448
1.438
1.430
1.437
a
0.034
0.020
0.012
0.008
0.006
0.003Figure 16
Variation ofM2 and a with F-T cycle numbers.Based on the two-yield surface model described inSection 4, model parameters were fitted as functional expressions with the number of F-T cycles N as variable. A modified model was then developed considering the influences of confining pressure and F-T cycles. Equations (23) and (24) are incremental equations of deviatoric stress, axial strain, and volumetric strain. Among them, coefficients A, B, C, and D were calculated using equations (16)∼(19). The calculation methods of parameters c(N), φ(N), M1(N), h(N), t(N), k(N), n(N), M2(N), and a(N) are shown in Figures 10, 11, 13, 15, and 16, respectively.


## 5.1. Determination ofc(N) and φ(N)

The internal friction angleφ and cohesion c were calculated according to the slope and intercept of the tangent line of Mohr’s stress circle as(25)φ=sin−1tanα,c=dcosφ,where α is the slope angle and d is the intercept of the envelope.The tangent lines of Mohr’s stress circles for different F-T cycle numbers are shown in Figure9, and shear strength parameters are summarized in Table 2. It was seen that both c and φ were decreased with the number of F-T cycles. The logistic function, as given by equation (26), was used to fit the regression relationships among φ, c, and number of F-T cycles, as shown in Figure 10:(26)y=x0−A11+N/A2A3+A1,where x0 is the shear strength parameters of non-F-T-experienced specimens and A1∼A3 are the fitting parameters.Figure 9
Variation of tangent lines of Mohr’s stress circles for different F-T cycle numbers.Table 2
Shear strength parameters for different F-T cycle numbers.
Shear strength
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
φ (°)
31.97
31.50
30.95
30.56
30.17
29.87
c (kPa)
21.05
19.85
19.45
19.39
19.30
19.24Figure 10
Variation of shear strength parameters for different F-T cycle numbers.


## 5.2. Determination ofM1(N)

The empirical equation ofM1(N) is derived as [13, 14](27)M1=1+0.25β2M,where β is the ratio of volumetric strain to axial strain at stress level 75% and β values vary with confining pressure; therefore, its average value was adopted. M is the slope of failure line as shown in Figure 8, which can be obtained by equation (28) combined with triaxial compression test results:(28)M=6sinφ3−sinφ.The values ofβ, M, and M1 are defined in Table 3. It can be seen that β was gradually increased with the number of F-T cycles, while M and M1 were decreased. The regression relationship between M1 and the number of F-T cycles was described by the logistic function, as shown in Figure 11.Table 3
The values ofβ, M, and M1 for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
β
0.460
0.503
0.534
0.546
0.574
0.616
M
1.286
1.265
1.241
1.224
1.207
1.194
M
1
1.340
1.328
1.312
1.296
1.286
1.284Figure 11
The variation ofM1 for different F-T cycle numbers.


## 5.3. Determination ofh(N) and t(N)

h(N) and t(N) were obtained through modifying equation (4) [13, 14]:(29)dεvp=hpadp0hpa+tp02.WhenBp=dp0/dεvp, equation (29) could be transformed as(30)Bppa=h+thp0pa.According to the previous experiences, plastic volumetric strainεvp accounts for about 65% of volumetric strain by increasing the stress level from 0 to 50%. Based on equation (2), p0 was obtained through p and q at the stress level of 50%. Then, Bp was calculated according to εvp and p0, specifically. The intercept and slope of the curve of (Bp/pa)0.5∼p0/pa were h0.5 and t/h0.5, respectively; therefore, h and t were also obtained.Curves of (Bp/pa)0.5∼p0/pa for different F-T cycle numbers are shown in Figure 12, and h and t values are given in Table 4. It was seen that h was gradually decreased with the number of F-T cycles while t was increased. The relationships among h, t, and the number of F-T cycles were fitted by the logistic function, as shown in Figure 13.Figure 12
Variation ofp0/pa∼(Bp/p0)0.5 with F-T cycle numbers.Table 4
Values ofh and t for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
h
884.23
603.59
486.60
385.73
302.80
197.91
t
6.538
6.756
7.548
7.755
7.970
7.778Figure 13
Variation ofh and t with F-T cycle numbers.


## 5.4. Determination ofk(N) and n(N)

The curve of lg (σ3/pa)∼lg (Ei/pa) based on triaxial test results is shown in Figure 14. It was seen that the spatial positions of fitting lines moved downwards by increasing the number of F-T cycles. The values of intercept k and slope n are given in Table 5. As can be seen, k was decreased with the number of F-T cycles while n was increased. The relationships among k, n, and the number of F-T cycles were described by the logistic function, as illustrated in Figure 15.Figure 14
Variation of lg (σ3/pa)∼lg (Ei/pa) with F-T cycle numbers.Table 5
The values ofk and n under different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
lgk
1.566
1.492
1.450
1.404
1.388
1.374
n
0.478
0.636
0.703
0.723
0.767
0.758
k
36.81
31.05
28.18
25.35
24.43
23.66
R
2
0.998
0.975
0.944
0.977
0.967
0.963Figure 15
Variations ofk and n with F-T cycle numbers.


## 5.5. Determination ofM2(N) and a(N)

M
2 and a were determined by transforming equations (5) and (8) [13, 14]:(31)p+prq=a2M2qεspG2+1M2.In equation (31), εsp was obtained using the following empirical equation:(32)εsp=0.3−0.1dεa,where d is the slope of the curve of εv∼εa at a stress levels of 75%∼95%. The above valuing method of d is actually an empirical one, and its reliability has been verified [13, 14]. Negative values of d meant that there was dilation. The values of d in each test curve were different, and therefore, its mean values were adopted.The intercept of the curve of (p + pr)/q∼(q/εspG)2 was 1/M2, and its slope was a2/M2; therefore, M2 and a were also obtained. The values of M2 and a for different F-T cycle numbers are summarized in Table 6. It was seen that M2 and a were decreased with the number of F-T cycles and the relationships among M2, a, and the number of F-T cycles were described by the logistic function, as shown in Figure 16.Table 6
Values ofM2 and a for different F-T cycle numbers.
Parameters
0 F-T
1 F-T
3 F-T
6 F-T
9 F-T
12 F-T
M
2
1.486
1.466
1.448
1.438
1.430
1.437
a
0.034
0.020
0.012
0.008
0.006
0.003Figure 16
Variation ofM2 and a with F-T cycle numbers.Based on the two-yield surface model described inSection 4, model parameters were fitted as functional expressions with the number of F-T cycles N as variable. A modified model was then developed considering the influences of confining pressure and F-T cycles. Equations (23) and (24) are incremental equations of deviatoric stress, axial strain, and volumetric strain. Among them, coefficients A, B, C, and D were calculated using equations (16)∼(19). The calculation methods of parameters c(N), φ(N), M1(N), h(N), t(N), k(N), n(N), M2(N), and a(N) are shown in Figures 10, 11, 13, 15, and 16, respectively.


## 6. Validation of the Modified Two-Yield Surface Model

The above incremental constitutive equations were calculated using a spreadsheet. Firstly, according to strain-controlled triaxial test results, equation (24) was used to calculate deviator stress increment dq corresponding to a constant axial strain increment value of 0.125%. Then, according to the obtained dq and equation (23), the volumetric strain corresponding to each level of axial strain was obtained.The relationships among deviatoric stress, axial strain, and volumetric strain under different conditions were obtained through the calculation program shown in Figure17. Among them, the results obtained after 4 F-T cycles were a supplementary verification test. It was seen that numerical results were basically consistent with experimental findings. For the same axial strain, deviatoric stress was decreased with the number of F-T cycles, while the volumetric strain was increased. It was demonstrated that the modified numerical model could predict the variations of deformation and strength characteristics of silty clays exposed to F-T cycles.Figure 17
Comparison of numerical and experimental values: (a) deviatoric stress and axial strain; (b) volumetric strain and axial strain.
(a)
(b)


## 7. Conclusions

(1)
The volume of silty clay was increased with the number of F-T cycles, which degraded soil mechanical behaviors to some extent. In the consolidated drained triaxial test, the discharge amount of pore water in the consolidating stage was increased with the number of F-T cycles. In the triaxial loading stage, specimens basically showed a shear contraction property. Meanwhile, volumetric strain was increased with the number of F-T cycles. Confining pressure in both consolidating and triaxial loading stages weakened the unfavorable effects of F-T cycles.
(2)
The incremental constitutive equations of deviatoric stress, axial strain, and volumetric strain increments have been derived in this paper by introducing the associated flow law to a typical two-yield surface model. Under the effect of F-T cycles, model parametersc, φ, h, k, M1, M2, and a were decreased while t and n were increased. Moreover, the logistic function was used to fit the regression relationship between the above model parameters and the number of F-T cycles. Furthermore, a modified numerical model considering F-T effects was developed by substituting model parameter expressions to incremental constitutive equations.
(3)
A program for the calculation of modified constitutive equations was developed. The comparison between numerical and test results showed that the proposed modified two-yield surface model well predicted the effect of F-T cycles on soil mechanical behaviors.



---

*Source: 9693156-2020-02-20.xml*