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 a p-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.
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 Duncan E-, 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 Figure 1. The collected soil samples were used as a subgrade filler for expressway and showed obvious F-T effects on both strength and deformation behaviors.
The grain composition of the soil sample is shown in Figure 2. 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.
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. Figure 3 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].
(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 Figure 4. 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.
3.2. Deviatoric Stress-Axial Strain Characteristics
Deviatoric stress-axial strain curves under different test conditions are shown in Figure 5. 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.
(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.
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 Table 1. 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.
The volumetric strain-axial strain relationships of the loading stage under different test conditions are shown in Figure 7. 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.
(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 and , where is related to compression and 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
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]:where M1 presents the elliptic shape change, p is the mean normal stress, q is the generalized shear stress, 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 for pr was stated as
Under isostatic triaxial testing conditions, the relationship between p0 and could be expressed in a hyperbolic form. Therefore, the hardening function corresponding to the first yield surface was written as [13, 14]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 equation f2 corresponding to shear expansion was a parabolic function expressed as [13, 14]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, is the hardening parameter of parabolic yield surface, and Ha2 is the hardening function of expansion yield surface.
G was calculated aswhere for the elastic modulus, we had E = 2.0Ei. Generally, the proportion of elastic strain in the total strain is small; therefore, is always set at 0.3. The initial tangent modulus Ei is related to confining pressure, and its empirical equation was derived aswhere 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]
4.2. Constitutive Equations in Increment Form
The three strain components of equation (1) have been determined in this chapter, where elastic strain component included volumetric elastic strain increment and generalized shear elastic strain increment . Calculation equations were as follows:where K is the bulk modulus.
By differentiating both sides of equation (2), we found
The plastic compression component of the first yield surface was consisted of two parts: volume change increment and shear strain increment . According to equations (4) and (10), the following equation was obtained for :
Assuming both yield surfaces to be applicable to associated flow law, the following equation was obtained:
The calculation equation of was obtained by combining equations (11) and (12).
The plastic compression components of the second yield surface included volumetric elastic strain and shear strain increments. Using the same derivation method as the first yield surface, the following equation was obtained for :
The calculation equation of was obtained by combining equations (13) and (14).
According to typical elastic-plastic theory, and dεs were defined aswhere A, B, C, and D were expressed as follows:
4.3. Equations of Triaxial Tests
In conventional consolidated drained triaxial tests, stress state is generally in A3 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:
In strain-controlled triaxial tests, the relationship of axial dεa and volumetric strain increments were
By substituting equation (21) into equation (15), the volumetric strain increment was found to be
By substituting equations (22) and (23) into equation (15), axial strain increment was derived as
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 cycles N 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 of c(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 aswhere α 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 Figure 9, 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:where x0 is the shear strength parameters of non-F-T-experienced specimens and A1∼A3 are the fitting parameters.
5.2. Determination of M1(N)
The empirical equation of M1(N) is derived as [13, 14]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:
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.
5.3. Determination of h(N) and t(N)
h(N) and t(N) were obtained through modifying equation (4) [13, 14]:
When , equation (29) could be transformed as
According to the previous experiences, plastic volumetric strain 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 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.
5.4. Determination of k(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.
5.5. Determination of M2(N) and a(N)
M2 and a were determined by transforming equations (5) and (8) [13, 14]:
In equation (31), was obtained using the following empirical equation:where d is the slope of the curve of ∼ε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/G)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.
Based on the two-yield surface model described in Section 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 Figure 17. 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.
(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 parameters c, φ, 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.
Data Availability
All data, models, and code generated or used during the study appear in the submitted article.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
T. H. conceived the study, carried out data analysis, and wrote the manuscript. D. W. implemented the experiments. H. W. provided the key suggestion and idea.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (Grant no. 41731281) and the Fundamental Research Funds for the Central Universities of China (Grant no. C17JB00520).