Dynamic constitutive model of frozen soil that considers the evolution of volume fraction of ice

A new constitutive model for frozen soils under high strain rate is developed. By taking the frozen soil as a composite material and considering the adiabatic temperature rise and interfacial debonding damage, the nonlinear dynamic response (NDR) of the frozen soil is predicted. At the same time, the relationship between instantaneous temperature and unfrozen water content is given, and an evolution rule of the volume fraction of ice particles is obtained. This relationship shows good agreement with experimental data. Using this new constitutive model, the stress–strain relationship of frozen soil under impact loading at temperatures of − 3 °C, − 8 °C, − 18 °C, and − 28 °C is calculated. There is good agreements between the results based on this new constitutive model and the data of dynamic impact.


Scientific Reports
| (2020) 10:20941 | https://doi.org/10.1038/s41598-020-77955-6 www.nature.com/scientificreports/ consider the nonlinear responses of frozen soil in the initial stage of increasing stress under impact loading, it is inadequate for dealing with the initial nonlinearity of frozen soil. A large number of studies 20,21 have shown that under impact loading, materials are in an adiabatic state and their internal temperature rises sharply. For frozen soil, the adiabatic temperature rise will change the volume fraction of ice and lead to change in the mechanical properties of the frozen soil. This can also explain the NDR of frozen soil materials. Therefore, in order to consider the dynamic mechanical properties of frozen soil with reasonable precision, it is necessary, when establishing a dynamic constitutive model of frozen soil, to consider the adiabatic temperature rise.
In the present study, taking these considerations into account, a dynamic constitutive model is developed to describe the nonlinear responses and the strain softening following the peak stress of frozen soil. To capture the nonlinear responses in the initial deformation stage, frozen soil is regarded as a composite material, wherein ice provides the reinforcing particles and soil is the matrix. Hence, the nonlinear responses of frozen soil before the peak stress are described by superposing the elastic responses of ice and the plastic responses of soil, which can be simulated using an elasto-plastic model with the modified Drucker-Prager (DP) yield criterion 22 . At the same time, in order to describe the phenomenon of softening of frozen soil after hardening, the interfacial debonding damage is taken into account. The adiabatic temperature rise is also considered. The proposed model solves a problem that existed in previous work 19 -inadequate prediction of initial nonlinearity.

Constitutive model
The influences of unfrozen water and air on the mechanical properties of frozen soil are much less than that of ice and soil 23 . Therefore, in this work the influences of unfrozen water and air on the mechanical properties of frozen soil are neglected; frozen soil is regarded as a composite material consisting of the reinforcing particles (ice) and the matrix (soil), with the representative volume element (RVE) of initial frozen soil as shown in Fig. 2 11 .
The stress-strain responses of frozen soil before the peak stress can be described by the following equation: (1) dσ all i = f ice kdσ ice i + f soil dσ soil i Figure 1. The NDR of frozen soil in the initial stages of deformation 11 . www.nature.com/scientificreports/ where dσ all i (I = 1, 2, …, 6) is the differential of the total stress tensor; dσ ice i and dσ soil i are the differentials of the stress tensor acting on ice and soil, respectively. The multipliers f ice and f soil are the volume fractions of ice and soil, respectively, with f ice + f soil = 1 . The parameter k denotes the extent of interfacial debonding damage.
The process of interfacial debonding damage in the frozen soil under impact loading has been described in reference 11 . The uniform stress in the fully bonded (i.e., intact) ice particle is denoted as σ ice i . The effective stress in the particle is simply calculated by [σ ice i ] = kσ ice i . (1) if k = 1 , the particle is intact; (2) if k = 0 , the particle is fully debonded; (3) if 0 < k < 1 , the particle is partially debonded.
The ice particle is isotropic and much stronger than the interface and the soil matrix, so it is assumed to be linearly elastic before being fully debonded. This means that the stress-strain relationship of an ice particle can be written: where C ij is the elastic matrix of ice and dε j is the differential of the total strain tensor (i, j = 1, 2,…, 6).
The strain-rate dependence and nonlinearity of the dynamic stress-strain responses of frozen soil is caused by the plasticity of the soil matrix, which is considered to be homogeneous and isotropic.
According to the associated flow rule, the differential of the plastic strain tensor acting on soil is related to a potential function as follows: In the stress space 24 , where f is the plastic potential function; ∂f ∂ε i controls the direction of plastic deformation, and d is a scalar function of proportionality.
The DP yield criterion is widely used for inhomogeneous materials; it is expressed as where I 1 = σ 1 + σ 2 + σ 3 is the first invariant of the stress tensor; J 2 = 1 2 s ij s ij is the second invariant of the deviatoric stress s ij ; and α and K are material parameters.
Modification of Eq. (5) from the stress space to the strain space gives where I 1e is the first invariant of the strain tensor; J 2e is the second invariant of the deviatoric strain; and ε = 2 9 ε ij ε ij is the effective strain. In Eq. (6), the DP yield criterion consists of a deviatoric deformation term and a linear dilatational term. This indicates that the DP yield criterion considers the dilatation of materials to be linear. However, the dilatation of composites is, in practice, often nonlinear. Therefore, in order to describe the nonlinear behavior of frozen soil more accurately, a quadratic dilatational term is added in the expression of the DP yield criterion: where parameters a , b , and c are associated with the deviatoric strain, linear dilatation, and quadratic dilatation, respectively. If c = 0 , Eq. (7) becomes the conventional DP yield criterion.
Because frozen soil is always regarded as an orthogonal material 25,26 , the differential of plastic work is given by If the relationship between the effective plastic stress and the effective strain is assumed to be where parameters n and A are related to the strain rate and temperature, respectively, the scalar function of proportionality becomes Using Eqs. (1), (2), (4), and (12), the model can be expressed as The above model framework is applicable to multiaxial cases, but because the impact loading is usually uniaxial, it can be simplified. In uniaxial compression tests, the shear strain is zero (i.e., ε 4 = ε 5 = ε 6 = 0 ), and the relationship between normal strains is where ν is the negative Poisson's ratio of frozen soil. Under uniaxial compression, the absolute value of the axial strain ε 1 is equal to the effective strain. Therefore, the potential function can be written as In addition, for an isotropic material, (20) ∂f ∂ε 1 = 1 and dσ 6 dε 6 = 0. www.nature.com/scientificreports/ Combination of Eqs. (19), (21) and (22) gives the value of parameter a as: Equation (15) yields

This means that
The values of parameters a , b and c can be determined by combining Eqs. (17), (23), and (25). Finally, turning to statistics, the probability distribution of the extent of interfacial debonding damage around the ice particles obeys the Weibull's distribution function at macroscopic scale 19,27 .
So the damage parameter k can be written as: where m and α are the parameters of Weibull's distribution function and ε e is the equivalent strain.

The evolution of the volume fraction of ice particles
As previously discussed, the test temperature is an important factor influencing the dynamic mechanical properties of frozen soil. The volume fraction of ice particles f ice in Eq. (13) is a key parameter. It depends, of course, on the test temperature, and can be calculated using the following formula : where S 0 is the water content of frozen soil, s water is the unfrozen water content, ρ water is unfrozen water density, and G and V are the weight and volume of the frozen soil specimen, respectively. It should be noted that this formula does not apply when the temperature is close to 0 °C. In a gas-free soil, the Clapeyron equation is commonly used to relate temperature and liquid (unfrozen water)-ice capillary pressure 28,29 where P water and P ice are unfrozen water and ice pressures, L fusion is the heat of fusion for water ice, ρ water is unfrozen water density, T is the instantaneous temperature, and T 0 is the nominal freezing temperature.
The unfrozen water content s water in this gas-free situation can then be related to the ice-water capillary pressure by where S is the soil moisture retention curve for unfrozen conditions, in this work the data in reference 28 is used; and β is the ratio of ice-water to water-air surface tensions for noncolloidal soils; it equals unity for colloidal soils 30 .
Combining Eqs. (28) and (29), the following equation is obtained: where H is the Heaviside function, which has been used to make Eq. (30) applicable to both frozen and unfrozen conditions. Equation (30) relates the unfrozen water fraction to temperature and to the soil moisture retention curve, which can be obtained by a soil physical-property test. Painter and Karra's research 27 shows that the formula is applicable to frozen soil only if the test temperature is well below 0 °C (< − 2 °C).
The instantaneous temperature T can be calculated using the following formula: and C 66 = 2G. www.nature.com/scientificreports/ where T initial is the initial temperature of frozen soil and T is the adiabatic temperature rise. Adiabatic temperature rise refers to the increase of material temperature caused by the heat generated by plastic deformation of materials, which cannot be rapidly spread to the outside world. The adiabatic temperature rise in the impact process is generally calculated by the following formula 21 where ρ is the density of frozen soil, C v is the specific heat of frozen soil, and η is the work-thermal conversion factor of the material.
According to the principle of meso mechanics, the physical parameters X of granular composites can be calculated using the following formula where f 1 and f 2 are the volume fractions of the enhanced phase and the matrix, respectively; and x 1 and x 2 are the physical parameters of the enhanced phase and the matrix, respectively. Equation (33) can be used to calculate the physical parameters of frozen soil such as density and specific heat.
The calculated adiabatic temperature-rise curve of frozen soil with moisture content of 10% at − 28 °C is shown in Fig. 3.
This figure shows that during the impact process, the adiabatic temperature rise of the frozen soil increases with the increase of strain and strain rate. At 838 s −1 strain rate, the temperature rises about 8 °C during the impact process. This phenomenon obviously has a great influence on the mechanical properties of frozen soil, which has high temperature sensitivity.

Model verification
The above equations have been incorporated into a constitutive model that couples a meso-mechanical equation of composite materials, a damage-evolution equation, an evolution equation of the volume fraction of unfrozen water, and an adiabatic temperature rise equation.
There are many material parameters used in the developed model. The values of those parameters can be obtained in the following ways: www.nature.com/scientificreports/ (1) Damage parameters (m and α) The frozen soil materials in this paper are consistent with those in reference 11 , so the values of damage parameters are also consistent with those in reference 11 ; (2) Physical parameters (ν, G, ρ, ρ water ,β, L fusion , C v ) These parameters can be obtained experimentally.
(3) Working-condition parameters (c, A, η) These parameters are related to the actual working conditions. Their values need to be determined according to details of the actual situation.
By referring to relevant literature as well as actual test and experimental conditions, the model parameters are obtained, as listed at Table 1.
Published experimental results 19 are used here to validate our approach. The dynamic stress-strain curves of frozen soil with moisture contents of 10% and 30% were predicted by the proposed model at − 3 °C, − 8 °C, − 18 °C, and − 28 °C.
The reasons for choosing these temperatures and moisture contents are as follows: (1) The applicable range of the model is temperature < − 2 °C. Therefore, in order to verify the predictive ability of this model, − 3 °C, close to the applicable temperature limit, is chosen for analysis. In actual freezing construction, the lowest temperature of active freezing is − 28 °C, so − 28 °C is chosen as the lowest research temperature. In addition, intermediate temperatures of − 8 and − 18 °C were selected to be included in the analysis.
(2) Freezing construction is generally used in soils with moisture content greater than 10%. Therefore, critical moisture content (10%) and a higher moisture content (30%) are selected for analysis.  Figures 4 and 5 show that for moisture contents of 10% and 30%, most of the predictions are accurate, but the softening rule of the − 18 °C, and − 28 °C with moisture content of 30% curves deviated from the predicted results. It may be caused by excessive ice content in frozen soil, that the evolution of debonding damage is faster than the prediction of the model. But on the whole, the accuracy of prediction results is acceptable. Figures 6  and 7 show that the model also well describes the initial nonlinear characteristics of frozen soil. This means that the proposed model is reasonable and quite accurate. Figures 4 and 6 also show that the slope of the initial stage of "warm" frozen soil increases significantly with increasing strain rate. However, the slope of the initial stage of "cold" frozen soil does not change with increasing strain rate. This is because the mechanical properties of frozen soil depend largely on the ice content. The strength and elastic modulus of ice are very high, which causes it to exhibit a linear elastic response during impact loading. This means that the strain rate has little influence on the mechanical properties of ice. As the temperature decreases, the ice content of frozen soil increases, thereby causing a decrease in the effect of the strain rate on the mechanical properties of frozen soil. This result coincides with the experimental results of Zhao et al. 31 . This indicates that it is very reasonable for the model to take into account the variation of ice content. www.nature.com/scientificreports/

Conclusions
(1) Existing experimental results of frozen soil show that no purely elastic stage occurs in the static and dynamic deformations of frozen soil. Therefore, in constructing the proposed model, it was necessary to consider the plastic behavior of frozen soil in the initial stage of loading. (2) A constitutive model relating temperature, unfrozen water content, adiabatic temperature rise, and interfacial debonding damage is found to provide a reasonable fit to available laboratory data. (3) The results show that the proposed model is capable of simulating the nonlinear and initial plastic characteristics of frozen soil, and can well describe the dynamic deformation of frozen soil at different temperatures. (4) The ranges of temperature, moisture content, and strain rate discussed in this paper are consistent with actual working condition in freezing construction. Therefore, in engineering practice, the research results of this paper can be used to predict the strength of artificial frozen soil produced by freezing construction. Appropriate excavation methods can then be selected according to the known strength of frozen soil. This will help to save construction costs and enhance construction safety. (5) It should be noted that, the proposed constitutive model is validated by uniaxial deformation. If it is used for three-dimensional deformation, further verification is needed.