A Damage Constitutive Model for a Rock under Compression after Freeze-Thaw Cycles Based on the Micromechanics

The freeze-thaw cycles will cause continuous damage to the rock, which is much related to the microcrack length, rock permeability, and frost heaving pressure. However, the failure mechanism of the rock under compression after freeze-thaw cycles is not very clear; therefore, it is studied with the damage theory here. First of all, according to the hydraulic pressure theory, the relationship between the frost heaving pressure and the microcrack propagation length in one single microcrack is established based on the elastoplastic mechanics and fracture theory. Second, by assuming the total strain of the rock under compression is comprised of the initial damage strain, elastic strain, additional damage strain, and plastic damage strain, a constitutive model for a rock based on the deformation and propagation of the microcrack under compression after freeze-thaw cycles is established. Finally, the proposed model is veri ﬁ ed with the test result. In all, the proposed model can perfectly re ﬂ ect the deterioration of the rock mechanical behavior under compression after the freeze-thaw cycles.


Introduction
Rock deterioration under freeze-thaw cycles is a hot issue in many engineering projects in cold regions. This deterioration proceeds with the freeze-thaw cycles, and the rock will gradually lose its stiffness and strength.
Water in the microcrack is assumed to be the key factor leading to the rock deterioration under freeze-thaw cycles. The deterioration mechanism of the porous media under freeze-thaw cycles was firstly described by Powers [1]. According to his theory, the frost heaving pressure will be generated by 9% volume expansion when water in the closed microcrack freezes into ice. It will make the microcrack propagate and cause damage to the rock. Thereafter, lots of laboratory tests have been done by many researchers, for instance, Altindag et al. [2], Momeni et al. [3], Nicholson and Nicholson [4], Zhang et al. [5], Wang et al. [6], Tounsi et al. [7], and Tang et al. [8]. All these results indicated that with the increasing number of freeze-thaw cycles, rock will deteriorate and degrade to some degree in terms of the compressive strength, elastic modulus, mass density, and so on.
No, many research works have been done in this field. Neaupane and Yamabeb [9] proposed a nonlinear elastoplastic constitutive relationship and a two-dimensional (plane stress) numerical modeling on the basis of the finite element method. With the continuum mechanics, Exadaktylos [10] established a coupled model for the saturated porous rocks under freeze-thaw cycles which can be used for analyzing the preliminary thawing experiments on the porous sandstone. According to the freeze-thaw cyclic fatigue tests on sandstone and shale, Tang et al. [8] assumed that the rock damage process included two coupled parts, e.g., the damage induced by the freeze-thaw cycles and the damage evolution caused by stress erosion. And then, a frost damage constitutive model is accordingly established. Based on the mass conservation law, energy conservation law, and the principle of static equilibrium, Kang et al. [11] studied the thermohydro-mechanical coupling mechanism and then proposed a new THM (thermo-hydro-mechanical) coupling model by considering phase change effect, in which the water migration caused by segregation potential and temperature gradient is described. By considering the coupling effects among fluid flow, heat transfer, crystallization, and deformation in porous media, Wu et al. [12] established a thermo-hydrosalt-mechanical coupled model for fully saturated porous media with phase change. From the viewpoint of the energy conservation law, mass conservation law, and the principle of static equilibrium considering water/ice phase change, Huang et al. [13] set up a fully coupled THM model and verified its validity with the laboratory test. Based on a microstructure-based random finite element model for the frozen soil, Dong and Yu [14] established a holistic model to simulate the temperature, stress, and deformation in frozen soil and implement a model to simulate frost heave and stress on water pipelines. In order to consider the migration of unfrozen water during freezing, Wang et al. [15] proposed a general thermo-mechanical-water migration coupled constitutive model to model mechanical degradation of rocks subjected to freeze-thaw cycles. Fan et al. [16] established a universal damage constitutive model under freeze-thaw and loading conditions based on the statistical damage constitutive model. Meanwhile, it can be seen that plasticity theory [17,18] is often adopted to describe the nonlinear mechanical behavior of rock-like materials under freeze-thaw cycles.
Rock properties including the microcrack size distribution and permeability also have significant influence on rock mass mechanical behavior [19][20][21][22]. Hori and Moriniro [23] treated the shape of the microcrack as an ellipse in rock and set up a micromechanical model for the microscopic process. Although they assumed that the rock damage was induced by the microcrack propagation due to the water freezing and movement, they did not discuss the influence of microcrack distribution and permeability on rock deterioration. Therefore, the main objective of this research is to present an elastoplastic damage model for the rock with random distribution of the microcrack and the plastic yield criterion of the homogeneous medium combined with the micromechanical damage model to simulate the rock plastic deformation. The total strain of the rock is assumed to be comprised of initial damage strain, elastic strain, additional damage strain, and plastic damage strain [24], where the initial damage strain caused by the freeze-thaw cycles is calculated by the initial compliance matrix which is the function of freeze-thaw cycles. Meanwhile, the Drucker-Prager criterion is adopted to describe the plastic behavior of rock under compression, in which the microcrack radius is assumed to obey an exponential law [25]. Finally, the validity of the proposed model is verified with the experiment results.

Propagation of One Single Microcrack under
Frost Heaving Pressure 2.1. The Relationship between the Microcrack Propagation Length and Ice Pressure in One Single Microcrack. Figure 1 illustrates a two-dimensional propagation model for one single microcrack under the frost heaving pressure. The microcrack is an ellipse, and its propagation under the frost heaving pressure will lead to the rock damage. When freezing, the frost heaving pressure p acts on the inner wall of the microcrack normally and evenly.
The following assumptions are made in this study: (1) the microcrack is always elliptical during the whole freeze-thaw cycles. That is to say, the shape of the microcrack is the same; only its size changes; (2) the rock particle is assumed to be unchanged; (3) the microcrack is always saturated, and the microcrack propagation obeys the linear elastic fracture mechanics; and (4) the propagation process of the microcrack is stable.
Water will change into ice with the volume expansion when temperature decreases to a certain degree. But because of the constraint of the microcrack wall, the stress induced by the ice volume expansion will act on the microcrack inner wall which will produce the elastic strain energy in the rock. When the stress intensity factor K I is larger than the rock fracture toughness K IC , the microcrack will propagate, which will lead to the release of the elastic strain energy. According to the Griffith energy balance theory, there is where W, E, are U are work done by the frost heaving pressure, elastic strain energy stored in the rock, and the reduced total potential energy of the whole system, respectively. Assume the elastic strain energy releases completely during the microcrack propagation, Equation (1) can be written into The work W done by the frost heaving pressure on the inner wall of the microcrack is expressed as where p is the frost heaving pressure, △b is the microcrack opening displacement increment, shown in Figure 1. The reduced total potential energy U of the whole system in Equation (2) can be expressed as

Geofluids
where △a is the microcrack propagation length and G is the microcrack Griffith energy release rate. There is a relationship between the volume of water and that of the ice [26]. Without considering the constraint of the microcrack wall, we assume the ice expansion volume is △V i under the free condition. However, in practice, the ice is loaded by the stress p, and accordingly, the corresponding volumetric strain ε v can be calculated: where E i and v i are the ice elastic modulus and Poisson's ratio, and here, they are assumed to be 600 MPa and 0.3, respectively. Then, the actual volume increment ΔV i ′ is where V i is the volume of water before freezing. According to the relationship of the water volume before and after the phase change, there is The microcrack propagation length can be obtained by combining Equations(3), (4), (6), and (7) where A = πG, B = πð2pab + aGÞ, and C = −2ΔV i ′pa. Solving Equation (8) yields For the rock especially the soft rock, the plastic zone will be formed near the microcrack tip when the microcrack wall is loaded by the frost heaving pressure. So, in order to satisfy the requirement of the linear elastic fracture mechanics, it is necessary to deal with the plastic zone near the microcrack tip with the equivalent method. The plastic region reduces the stiffness of the rock which is equivalent to a longer microcrack. The equivalent propagation length a′ of the microcrack is [27] where a is the original microcrack length and r y is the length of the plastic zone near the microcrack tip shown in Figure 1, which can be expressed as for a plane stress issue [27].
where σ s is the rock yield strength and K I is the first stress intensity factor at the microcrack tip. So after amendment, the actual propagation length Δa ′ of the microcrack shown in Figure 1 is The shape of the microcrack after propagation is shown as the dotted line in Figure 1.
After m freeze-thaw cycles, the microcrack half-length a m is The microcrack propagation corresponds to the rock damage, and accordingly, the rock elastic modulus and compressive strength will also decrease.

The Frost Heaving Pressure. According to Walder and
Hallet [28], the frost heaving pressure p i is related to the duration time of the low temperature, the value of temperature, the volume of ice and water, and flow resistance, which can be calculated by where the characteristic time τ is [28] τ = 8 3π where p i ðtÞ is the frost heaving pressure of ice at time t; Lð −T c Þ/v s T a = 1:1 MPa/°C × ð−T c Þ. Lð−T c Þ is the fusion heat of ice at T = T c , kJ/mol; v s and v L are the relative volume of ice and water, respectively; when T < −1°C, v L = 0:07 and v s = 0:93 [29]; T a is the absolute temperature, 273.15 K; p 0 is the initial frost heaving pressure, and according to the experiment result, it is 2 MPa; v is the rock Poisson ratio; μ is the 3 Geofluids rock shear modulus, MPa; R f is the flow resistance, Pa·s/m; g is the gravitational acceleration, m/s 2 ; a is the microcrack half-length. Here, the freezing and thawing temperatures are adopted to be -20°C and 20°C, respectively. The time of one freeze-thaw cycle is 12 h. The detailed calculation method of R f has been introduced by Walder and Hallet [28].
Therefore, the variation of the frost heaving pressure with time is shown in Figure 2.
The work done by the frost heaving pressure along the microcrack inner wall for one freeze-thaw cycle is where t i = 43200s and p i ðtÞ can be solved with Equation (14).

Establishment of the Rock Constitutive
Model under Compression after Freeze-Thaw Cycles 3.1. Establishment of the Rock Constitutive Model. It is assumed that the water/ice phase change in freeze-thaw condition is the main reason leading to rock deterioration. The frost heaving pressure is generated by 9% volume expansion of freezing water in closed microcrack. The pressure makes the microcrack propagate and when the temperature rises, the melt water will go into the newly formed microcracks. The repeated freeze-thaw cycles cause continuous damage to the rock. Based on this viewpoint, the propagation of one single microcrack under the frost heaving pressure is studied, and the relationship between the propagation length of the microcrack and the frost heaving pressure is obtained.
Because new damage continuously occurs under freezethaw cycles, the elastoplastic theory is adopted to study the rock damage mechanical behavior, and finally, a new constitutive model for a rock based on the deformation and propagation of microcracks under compression after freeze-thaw cycles is proposed.

Strain Decomposition.
In this proposed model, the total strain of the rock under compression can be decomposed into the following four components such as the initial dam-age strain, elastic strain, plastic strain, and additional damage strain induced by the microcrack propagation. It can be expressed as where ε is the total strain, ε d is the initial damage strain, ε e is the elastic strain, ε da is the additional damage strain, and ε p is the plastic strain. Their calculation methods are discussed below.
3.2.1. Initial Damage Strain Induced by Freeze-Thaw Cycles.
The initial damage strain induced by freeze-thaw cycles is where ½C d is the initial damage compliance matrix. Assume the half-length of the i th microcrack becomes a m after m freeze-thaw cycles. α i is the orientation of the i th microcrack; then, the initial damage compliance matrix due to one single microcrack is given by [25] C where [C 0 ] is the elastic compliance matrix, [A i ] is the transformation matrix, and ½A i −1 is its inverse matrix.
where C i n , C i s , K i n , and K i s are the compression transferring coefficient, shear transferring coefficient, normal stiffness, and shear stiffness of the i th crack, respectively. d and h are width and height of the rock sample, respectively.
The total compliance matrix including N microcracks whose half-length a m is The orientation and size of the microcracks in the rock are assumed to be random, which can be expressed with a 4 Geofluids probability density function ρða, αÞ = ρðaÞρðαÞ. According to the distribution law of the microcrack, its orientation and size satisfy the following normalization condition [30]: Assume the total number of the microcracks is N c , and then, the number N of the microcracks with half-length of a can be expressed as Assume the orientation of the microcrack evenly distributes in all directions [31], and according to the definition of the density function, ρðαÞ = 1 is obtained. Therefore, the total compliance matrix including all microcracks is With Equation (24), the initial damage compliance matrix for different freeze-thaw cycles can be calculated and then, the initial damage strain can be finally obtained.

Elastic Strain.
For a plane stress issue, the elastic constitutive relationship of rock is [32] where ½ε e is the elastic strain matrix, ½ε e = ½ε 11 ε 33 ε 13 T , ½σ is the stress matrix, ½σ = σ 11 σ 33 σ 13 ½ , and ½C 0 is the elastic compliance matrix, where E, μ, and υ are the rock elastic modulus, shear modulus, and Poisson ratio, respectively.

Additional Damage Strain due to Compression.
Under compression, the microcrack will firstly close, and then, the friction occurs along the microcrack face. When the shear stress along the microcrack face is larger than the friction, the wing crack will initiate and propagate from the microcrack tip. The propagation of the microcrack will lead to the decrease in the rock elastic modulus, strength, and increase in the rock permeability. The sliding microcrack model under compression is shown in Figure 3.
After m freeze-thaw cycles, the microcrack length becomes a m , and for this moment, the first stress intensity factor K I at the microcrack tip is [31,33] where τ * is the effective shear stress on the microcrack face, τ * = τ m − f σ m , τ m and σ m are the shear and normal stresses on the microcrack face, respectively, and f is the friction coefficient of the microcrack face.
When K Ι ≥ K IC (K IC is the rock fracture toughness), the microcrack begins to propagate, and the stress intensity factor K W I at the wing crack tip is [31,33]. where Þcos 2αÞ: The wing crack will stop propagating when K W Ι ≤ K IC , so l can be calculated with Equation (29).
The additional damage of a rock is due to the wing crack propagation. According to Li and Lajtai [34], the where ε * 1 and ε * 3 are the strains along the vertical and horizontal directions, respectively; λ = sin α cos α − f cos 2 α; γ = −cos α sin α − f sin 2 α. χ is defined as the initial microcrack density and is expressed as χ = Na 2 /V. N is the total number of the microcracks in a two-dimensional body of unit thickness whose volume is V = 2 h × 2d.
Finally, the total rock strain due to compression considering the normalization condition is 3.2.4. The Plastic Strain due to Compression. The Drucker-Prager model is adopted to describe the rock plastic behavior, whose yield function and plastic potential function are where ð3 − sin φÞ, β 2 = 2 sin ψ/ ffiffi ffi 3 p ð 3 − sin ψÞ, and φ and ψ are the friction angle and dilation angle, respectively. κ is the hardening function, which can be expressed as [35] where a 1 , a 2 , and a 3 are the constants which can be obtained by fitting with the uniaxial compressive stressstrain curve. σ 0 = 6c cos φ/ ffiffi ffi 3 p ð3 − sin φÞ; c is the rock cohesion.
According to Tan et al. [36], the rock internal friction angle is basically the same, but the rock cohesion decreases with the freeze-thaw cycles and obeys the following exponential function: where c 0 is the original cohesion strength (before freeze-thaw cycle) and cðmÞ is the cohesion after the m th freeze-thaw cycle.
The plastic strain rate is where _ λ is a proportion coefficient.

The Number of the Microcracks in a Rock.
The probability function [37] of the microcrack length can be expressed as where a c is the characteristic length of the microcrack, a c = Ð a max a min ρðaÞada. N c is the total number of the microcracks per unit volume, which is determined by the total volume of the microcrack, V c = 2πbN c Ð a max a min ρðaÞa 2 da.
Finally, the curve of probability density function of the microcracks is shown in Figure 4.
According to the experiment result of Rostásy et al. [38], the total number of the microcracks basically remains the same after freeze-thaw cycles. The longer microcracks will propagate, while the shorter ones will close under the extrusion of other microcracks. Therefore, it is assumed that the total volume of the microcrack is the same during the freeze-thaw cycles and compression. So the total volume V c of the microcracks can be expressed with the rock void porosity e and the total volume V of the rock sample, namely, V c = eV.

Geofluids
The first stress intensity factor K I of the microcrack with different lengths at the same stress condition can be expressed as Equation (38) can be changed into where τ * = ð1/2Þ½ðσ 1 − σ 3 Þ sin 2φ − f ðσ 1 + σ 3 + ðσ 1 − σ 3 Þ cos 2φÞ where a cr is the critical length of the microcrack which becomes active at the condition of σ 1 and σ 3 when As shown in Figure 5, the microcrack whose length is greater than a cr will propagate.
The total number of microcracks that will be actually activated is given by N = N c Ð a max a min ρðaÞda, as shown in Figure 6.
3.2.6. The Numerical Algorithm of the Proposed Model. The calculation of the rock plastic strain at different times can be calculated with the semi-implicit return graphical algorithm [17], shown in Figure 7. The plastic variables σ n+1 , ε p n+1 , κ n+1 , r n+1 , and h n+1 at t n+1 are determined by integration flow rule and hardening law and Δσ tr , σ n , κ n , r n = ∂F/∂σj σ=σ n , and h n = ∂κ/∂λj λ=λ n at time t n . The main steps are as follows: (1) Update strain tensor, and calculate the elastic strain at time t n+1 ε n+1 = ε n + Δε, where C is the initial elastic compliance matrix including the initial damage compliance tensor caused by freeze-thaw and elastic compliance tensor.
(2) Update stresses and plastic strain and hardening function κ n+1 ε p n+1 = ε p n + Δλ n+1 r n , κ n+1 = κ n Δλ n+1 h n where r n = ∂G/∂σj σ=σ n and h n = ∂κ/∂κj λ=λ n . Calculate the plastic internal variable Δλ n+1 : where (3) Update N n+1 which is the number of microcracks that begin propagating At time t n+1 , the updated normal stress is σ n+1 , under which the critical initiation length ða cr Þ n+1 of the microcrack is calculated by where τ * ð Þ n+1 = 1 2 (4) Update the additional damage strain ðε ad Þ n+1 The additional damage strain matrix of ε ad can be expressed as Update the length of wing crack l n+1 where K W I = K IC . Update the additional strain

Verification of the Proposed Model
In order to verify the proposed model, the experiment by Tan et al. [36] is taken for an example. The rock type is granite, which is obtained from Galongla mountain in Tibet of China, where a highway tunnel passes through the mountain and it is very cold in winter. The tested samples are prepared as cylinders with 50 mm in diameter and 100 mm high. The compression tests with confining pressure 5 MPa and 10 MPa on the rock sample are done with a multifunction rock mechanics test machine, and the corresponding stress-strain curves are shown in Figures 8 and 9. The calculation parameters are shown in Table 1.
The characteristic length a c of the microcrack under 0, 50, and 100 freeze-thaw cycles is solved to be 0.1 mm, 0.41 mm and 1.64 mm, respectively, with Equation (13). According to the proposed model, the complete stressstrain curve of the rock under compression with confining pressure σ 3 = 5 MPa/10 MPa is shown in Figures 8 and 9. It can be seen that the simulated stress-strain curves agree well with the tested ones especially when σ 3 = 10 MPa. Meanwhile, with the increasing freeze-thaw cycles, both the climax strength and slope of the stress-strain curve decrease; that is to say, the rock compressive strength and elastic modulus both decrease. It indicates that the freeze-thaw cycles have much effect on the rock mechanical behavior.

Conclusions
(1) Based on the fracture mechanics, the calculation method of the microcrack propagation length induced by the freeze-thaw cycles is proposed. Meanwhile, the variation of the frost heaving pressure with time during one freeze-thaw cycle is also calculated (2) In the framework of the fracture and damage mechanics, the total strain of the rock under compression after freeze-thaw cycles can be decomposed into the initial damage strain, elastic strain, plastic strain, and additional damage strain. And their calculation methods are discussed in detail. Finally, a constitutive model for a rock based on the deformation and propagation of microcracks under compression after freeze-thaw cycles is established (3) By utilizing the semi-implicit algorithm, the stressstrain relationship of the proposed model is calculated. The comparison of the theoretical results of the proposed method and test ones shows that they agree well with each other. Overall, the proposed method provides a new way to simulate the mechanical behavior of a rock under compression after freeze-thaw cycles

Data Availability
All data generated or analyzed during this study are included in this manuscript.

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