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

Authors: Hongyan Liu; Xiaochen Zhang; Xidong Yan
Journal: Geofluids
 (2021)

---

## Abstract

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 verified with the test result. In all, the proposed model can perfectly reflect the deterioration of the rock mechanical behavior under compression after the freeze-thaw cycles.


---

## Body

## 1. 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 thermo-hydro-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-hydro-salt-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–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.


## 2. 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

Figure1 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.Figure 1
The microcrack propagation model under the frost heaving pressure.2a and 2b are the length and width of the microcrack, respectively. △a and △b are the propagation length of the microcrack along two directions, respectively. ry is the length of the plastic zone near the microcrack tip. △a′ is the equivalent propagation length of the microcrack.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 factorKI is larger than the rock fracture toughness KIC, the microcrack will propagate, which will lead to the release of the elastic strain energy. According to the Griffith energy balance theory, there is
(1)W=E−U,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
(2)W=−U.The workW done by the frost heaving pressure on the inner wall of the microcrack is expressed as
(3)W=2×p×2a×Δb=4paΔb,where p is the frost heaving pressure, △b is the microcrack opening displacement increment, shown in Figure 1.The reduced total potential energyU of the whole system in Equation (2) can be expressed as
(4)U=−2G×Δa,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 △Vi 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:
(5)εv=31−2νiEip,where Ei and vi 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ΔVi′ is
(6)ΔVi′=ΔVi−Viε,where Vi is the volume of water before freezing.According to the relationship of the water volume before and after the phase change, there is(7)πab+ΔVi′=πa+Δab+Δb.The microcrack propagation length can be obtained by combining Equations(3), (4), (6), and (7)
(8)AΔa2+BΔa+C=0,where A=πG, B=π2pab+aG, and C=−2ΔVi′pa.Solving Equation (8) yields
(9)Δa=−2abp+aGπ+2abp+aG2π2+8Gpa×ΔVi′2Gπ.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 lengtha′ of the microcrack is [27]
(10)a′=a+ry,where a is the original microcrack length and ry 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].
(11)ry=12πKIσs2,where σs is the rock yield strength and KI 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
(12)Δa′=Δa+ry.The shape of the microcrack after propagation is shown as the dotted line in Figure1.Afterm freeze-thaw cycles, the microcrack half-length am is
(13)am=am−1+−2am−1bp+am−1Gπ+2am−1bp+am−1G2π2+8Gpam−1×ΔVi′2Gπ.The microcrack propagation corresponds to the rock damage, and accordingly, the rock elastic modulus and compressive strength will also decrease.


### 2.2. The Frost Heaving Pressure

According to Walder and Hallet [28], the frost heaving pressure pi 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
(14)pit=L−TcvsTa1−e−t/τ+p0e−t/τ,where the characteristic time τ is [28]
(15)τ=83π1−νμgavLvs2Rf,where pit is the frost heaving pressure of ice at time t; L−Tc/vsTa=1.1MPa/°C×−Tc. L−Tc is the fusion heat of ice at T=Tc, kJ/mol; vs and vL are the relative volume of ice and water, respectively; when T<−1°C, vL=0.07 and vs=0.93 [29]; Ta is the absolute temperature, 273.15 K; p0 is the initial frost heaving pressure, and according to the experiment result, it is 2 MPa; v is the rock Poisson ratio; μ is the rock shear modulus, MPa; Rf is the flow resistance, Pa·s/m; g is the gravitational acceleration, m/s2; 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 Rf has been introduced by Walder and Hallet [28].Therefore, the variation of the frost heaving pressure with time is shown in Figure2.Figure 2
Variation of the frost heaving pressure with time during one freeze-thaw cycle.The work done by the frost heaving pressure along the microcrack inner wall for one freeze-thaw cycle is(16)W=∫0ti4apitΔbtidt=4aΔbti∫0tipitdt,where ti=43200s and pit can be solved with Equation (14).


## 2.1. The Relationship between the Microcrack Propagation Length and Ice Pressure in One Single Microcrack

Figure1 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.Figure 1
The microcrack propagation model under the frost heaving pressure.2a and 2b are the length and width of the microcrack, respectively. △a and △b are the propagation length of the microcrack along two directions, respectively. ry is the length of the plastic zone near the microcrack tip. △a′ is the equivalent propagation length of the microcrack.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 factorKI is larger than the rock fracture toughness KIC, the microcrack will propagate, which will lead to the release of the elastic strain energy. According to the Griffith energy balance theory, there is
(1)W=E−U,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
(2)W=−U.The workW done by the frost heaving pressure on the inner wall of the microcrack is expressed as
(3)W=2×p×2a×Δb=4paΔb,where p is the frost heaving pressure, △b is the microcrack opening displacement increment, shown in Figure 1.The reduced total potential energyU of the whole system in Equation (2) can be expressed as
(4)U=−2G×Δa,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 △Vi 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:
(5)εv=31−2νiEip,where Ei and vi 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ΔVi′ is
(6)ΔVi′=ΔVi−Viε,where Vi is the volume of water before freezing.According to the relationship of the water volume before and after the phase change, there is(7)πab+ΔVi′=πa+Δab+Δb.The microcrack propagation length can be obtained by combining Equations(3), (4), (6), and (7)
(8)AΔa2+BΔa+C=0,where A=πG, B=π2pab+aG, and C=−2ΔVi′pa.Solving Equation (8) yields
(9)Δa=−2abp+aGπ+2abp+aG2π2+8Gpa×ΔVi′2Gπ.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 lengtha′ of the microcrack is [27]
(10)a′=a+ry,where a is the original microcrack length and ry 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].
(11)ry=12πKIσs2,where σs is the rock yield strength and KI 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
(12)Δa′=Δa+ry.The shape of the microcrack after propagation is shown as the dotted line in Figure1.Afterm freeze-thaw cycles, the microcrack half-length am is
(13)am=am−1+−2am−1bp+am−1Gπ+2am−1bp+am−1G2π2+8Gpam−1×ΔVi′2Gπ.The microcrack propagation corresponds to the rock damage, and accordingly, the rock elastic modulus and compressive strength will also decrease.


## 2.2. The Frost Heaving Pressure

According to Walder and Hallet [28], the frost heaving pressure pi 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
(14)pit=L−TcvsTa1−e−t/τ+p0e−t/τ,where the characteristic time τ is [28]
(15)τ=83π1−νμgavLvs2Rf,where pit is the frost heaving pressure of ice at time t; L−Tc/vsTa=1.1MPa/°C×−Tc. L−Tc is the fusion heat of ice at T=Tc, kJ/mol; vs and vL are the relative volume of ice and water, respectively; when T<−1°C, vL=0.07 and vs=0.93 [29]; Ta is the absolute temperature, 273.15 K; p0 is the initial frost heaving pressure, and according to the experiment result, it is 2 MPa; v is the rock Poisson ratio; μ is the rock shear modulus, MPa; Rf is the flow resistance, Pa·s/m; g is the gravitational acceleration, m/s2; 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 Rf has been introduced by Walder and Hallet [28].Therefore, the variation of the frost heaving pressure with time is shown in Figure2.Figure 2
Variation of the frost heaving pressure with time during one freeze-thaw cycle.The work done by the frost heaving pressure along the microcrack inner wall for one freeze-thaw cycle is(16)W=∫0ti4apitΔbtidt=4aΔbti∫0tipitdt,where ti=43200s and pit can be solved with Equation (14).


## 3. 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 freeze-thaw 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.


### 3.2. 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 damage strain, elastic strain, plastic strain, and additional damage strain induced by the microcrack propagation. It can be expressed as(17)ε=εd+εe+εda+εp,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(18)εd=Cdσ,where Cd is the initial damage compliance matrix.Assume the half-length of theith microcrack becomes am after m freeze-thaw cycles. αi is the orientation of the ith microcrack; then, the initial damage compliance matrix due to one single microcrack is given by [25]
(19)C=C0+Ai−1TΔCiAi−1,where [C0] is the elastic compliance matrix, [Ai] is the transformation matrix, and Ai−1 is its inverse matrix.
(20)Ai=cos2αisin2αi−sin2αisin2αicos2αisin2αi12sin2αi−12sin2αicos2αi,Ai−1=cos2αisin2αisin2αisin2αicos2αi−sin2αi12sin2αi−12sin2αicos2αi,ΔCi=0000CniamiKni2dh000CsiamiKsi2dhwhere Cni, Csi, Kni, and Ksi are the compression transferring coefficient, shear transferring coefficient, normal stiffness, and shear stiffness of the ith crack, respectively. d and h are width and height of the rock sample, respectively.The total compliance matrix includingN microcracks whose half-length am is
(21)C=C0+∑i=1NAi−1TΔCiAi−1.The orientation and size of the microcracks in the rock are assumed to be random, which can be expressed with a probability density functionρa,α=ρaρα. According to the distribution law of the microcrack, its orientation and size satisfy the following normalization condition [30]:
(22)∫aminamaxρada∫0π/2ραsinαdα=1.Assume the total number of the microcracks isNc, and then, the number N of the microcracks with half-length of a can be expressed as
(23)N=Nc∫aminaρada.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
(24)C=C0+Nc∫aminamaxAi−1TΔCiAi−1ρada∫0π/2sinαdα.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.


#### 3.2.2. Elastic Strain

For a plane stress issue, the elastic constitutive relationship of rock is [32]
(25)εe=C0σ,where εe is the elastic strain matrix, εe=ε11ε33ε13T, σ is the stress matrix, σ=σ11σ33σ13, and C0 is the elastic compliance matrix,
(26)C0=1E−νE0−νE1E0001μ,where E, μ, and υ are the rock elastic modulus, shear modulus, and Poisson ratio, respectively.


#### 3.2.3. 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 Figure3.Figure 3
The sliding microcrack model under compression.2a and α are the length and dip angle of the microcrack, respectively. l is the wing crack length. 2h and 2d are the height and width of the model, respectively.Afterm freeze-thaw cycles, the microcrack length becomes am, and for this moment, the first stress intensity factor KI at the microcrack tip is [31, 33]
(27)KI=23τ∗πam,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.
(28)τm=σ1−σ32sin2α,σm=σ1+σ32−σ1−σ32cos2α.WhenKΙ≥KIC (KIC is the rock fracture toughness), the microcrack begins to propagate, and the stress intensity factor KIW at the wing crack tip is [31, 33].
(29)KIW=2aτ∗cosαπl−σ3πl,where
(30)τ∗=12σ1−σ3sin2α−fσ1+σ3+σ1−σ3cos2α.The wing crack will stop propagating whenKΙW≤KIC, 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 macroscopic strain of a rock can be calculated.
(31)ε1∗=8λχcosαE2τ∗cosαπlnla−σ3la−1,ε3∗=χE16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗,where ε1∗ and ε3∗ are the strains along the vertical and horizontal directions, respectively; λ=sinαcosα−fcos2α; γ=−cosαsinα−fsin2α. χ is defined as the initial microcrack density and is expressed as χ=Na2/V. N is the total number of the microcracks in a two-dimensional body of unit thickness whose volume is V=2h×2d.Finally, the total rock strain due to compression considering the normalization condition is(32)ε1=Nc∫aminamaxρada∫0π/2ραsinαdα⋅8λχcosαE×2τ∗cosαπlnla−σ3la−1,ε3=Nc∫aminamaxρada∫0π/2ραsinαdα⋅χE×16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗.


#### 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(33)Fσ,κ=β1I1+J2−κ,Gσ,κ=β2I1+J2−κ,where I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, β1=2sinφ/33−sinφ, β2=2sinψ/33−sinψ, and φ and ψ are the friction angle and dilation angle, respectively. κ is the hardening function, which can be expressed as [35]
(34)κ=σ0+a1λexpa2I1−a3λ,where a1, a2, and a3 are the constants which can be obtained by fitting with the uniaxial compressive stress-strain curve. σ0=6ccosφ/33−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:
(35)cm=c0e−0.072m,where c0 is the original cohesion strength (before freeze-thaw cycle) and cm is the cohesion after the mth freeze-thaw cycle.The plastic strain rate is(36)ε̇p=λ̇∂G∂σ,where λ̇ is a proportion coefficient.


#### 3.2.5. The Number of the Microcracks in a Rock

The probability function [37] of the microcrack length can be expressed as
(37)ρa=−1acexp−aac,amin≤a≤amax,0,otherwise,where ac is the characteristic length of the microcrack, ac=∫aminamaxρaada. Nc is the total number of the microcracks per unit volume, which is determined by the total volume of the microcrack, Vc=2πbNc∫aminamaxρaa2da.Finally, the curve of probability density function of the microcracks is shown in Figure4.Figure 4
Probability density functions of the microcracks.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 Vc of the microcracks can be expressed with the rock void porosity e and the total volume V of the rock sample, namely, Vc=eV.The first stress intensity factorKI of the microcrack with different lengths at the same stress condition can be expressed as
(38)KI=23τ∗πa.Equation (38) can be changed into
(39)3KIC2τ∗21π=acr,where τ∗=1/2σ1−σ3sin2φ−fσ1+σ3+σ1−σ3cos2φ where acr is the critical length of the microcrack which becomes active at the condition of σ1 and σ3 when KΙ=KIC.As shown in Figure5, the microcrack whose length is greater than acr will propagate.Figure 5
Relationship between the microcrack critical length and axial stress.The total number of microcracks that will be actually activated is given byN=Nc∫aminamaxρada, as shown in Figure 6.Figure 6
Relationship between the microcrack initiation number and axial stress.


#### 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, εn+1p, κn+1, rn+1, and hn+1 at tn+1 are determined by integration flow rule and hardening law and Δσtr, σn, κn, rn=∂F/∂σσ=σn, and hn=∂κ/∂λλ=λn at time tn. The main steps are as follows:
(1)
Update strain tensor, and calculate the elastic strain at timetn+1(40)εn+1=εn+Δε,Δσn+1tr=C−1:Δε,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(41)εn+1p=εnp+Δλn+1rn,κn+1=κn∆λn+1hnσn+1=C−1:εn+1−εn+1pFn+1=Fσn+1,κn+1=0where rn=∂G/∂σσ=σn and hn=∂κ/∂κλ=λn.Figure 7
Graphic returns method of semi-implicit integration algorithm.Calculate the plastic internal variable∆λn+1:
(42)an+1k+C−1:∆σn+1k+δλkrn=0,(43)bn+1k−∆κn+1k+δλn+1khn=0,(44)Fn+1k+Fσk:∆σn+1k+Fκk∙∆κn+1k=0,(45)ΔσΔκ=−Aka~k−δλkAkr~n,where
(46)Ak=C00−I,a~k=akbk=00,r~n=rnhn..δλk can be obtained with Equation (45)
(47)δλ=Fk∂F:Ak:r~n,where ∂F=FσFκ.
(48)∆λn+1k+1=Δλnk+δλk.(3)
UpdateNn+1 which is the number of microcracks that begin propagatingAt timetn+1, the updated normal stress is σn+1, under which the critical initiation length acrn+1 of the microcrack is calculated by
(49)acrn+1=1π3KIC2τ∗n+12,where
(50)τ∗n+1=12σ1n+1tr−σ3sin2φ−μσ1n+1tr+σ3+σ1n+1tr−σ3cos2φ,Nn+1=Nc∫acrn+1amaxρada.(4)
Update the additional damage strainεadn+1The additional damage strain matrix ofεad can be expressed as
(51)εad=ε1ε30.Update the length of wing crackln+1(52)KIW=2aτ`n+1cosφπln+1−σ3n+1πln+1,where KIW=KIC.Update the additional strain(53)ε1n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅8λacos2αE×2τ∗n+1cosαπlnln+1a−σ3ln+1a−1,ε3n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅a2E×16τ∗n+1γcos2απlnln+1a+σ3πln+12a2−1−8cosφln+1a−1σ3γ+τ∗n+1.(5)
Update the stress(54)σn+1∗=σn+1−C−1:εadn+1


## 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 freeze-thaw 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.


## 3.2. 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 damage strain, elastic strain, plastic strain, and additional damage strain induced by the microcrack propagation. It can be expressed as(17)ε=εd+εe+εda+εp,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(18)εd=Cdσ,where Cd is the initial damage compliance matrix.Assume the half-length of theith microcrack becomes am after m freeze-thaw cycles. αi is the orientation of the ith microcrack; then, the initial damage compliance matrix due to one single microcrack is given by [25]
(19)C=C0+Ai−1TΔCiAi−1,where [C0] is the elastic compliance matrix, [Ai] is the transformation matrix, and Ai−1 is its inverse matrix.
(20)Ai=cos2αisin2αi−sin2αisin2αicos2αisin2αi12sin2αi−12sin2αicos2αi,Ai−1=cos2αisin2αisin2αisin2αicos2αi−sin2αi12sin2αi−12sin2αicos2αi,ΔCi=0000CniamiKni2dh000CsiamiKsi2dhwhere Cni, Csi, Kni, and Ksi are the compression transferring coefficient, shear transferring coefficient, normal stiffness, and shear stiffness of the ith crack, respectively. d and h are width and height of the rock sample, respectively.The total compliance matrix includingN microcracks whose half-length am is
(21)C=C0+∑i=1NAi−1TΔCiAi−1.The orientation and size of the microcracks in the rock are assumed to be random, which can be expressed with a probability density functionρa,α=ρaρα. According to the distribution law of the microcrack, its orientation and size satisfy the following normalization condition [30]:
(22)∫aminamaxρada∫0π/2ραsinαdα=1.Assume the total number of the microcracks isNc, and then, the number N of the microcracks with half-length of a can be expressed as
(23)N=Nc∫aminaρada.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
(24)C=C0+Nc∫aminamaxAi−1TΔCiAi−1ρada∫0π/2sinαdα.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.


### 3.2.2. Elastic Strain

For a plane stress issue, the elastic constitutive relationship of rock is [32]
(25)εe=C0σ,where εe is the elastic strain matrix, εe=ε11ε33ε13T, σ is the stress matrix, σ=σ11σ33σ13, and C0 is the elastic compliance matrix,
(26)C0=1E−νE0−νE1E0001μ,where E, μ, and υ are the rock elastic modulus, shear modulus, and Poisson ratio, respectively.


### 3.2.3. 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 Figure3.Figure 3
The sliding microcrack model under compression.2a and α are the length and dip angle of the microcrack, respectively. l is the wing crack length. 2h and 2d are the height and width of the model, respectively.Afterm freeze-thaw cycles, the microcrack length becomes am, and for this moment, the first stress intensity factor KI at the microcrack tip is [31, 33]
(27)KI=23τ∗πam,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.
(28)τm=σ1−σ32sin2α,σm=σ1+σ32−σ1−σ32cos2α.WhenKΙ≥KIC (KIC is the rock fracture toughness), the microcrack begins to propagate, and the stress intensity factor KIW at the wing crack tip is [31, 33].
(29)KIW=2aτ∗cosαπl−σ3πl,where
(30)τ∗=12σ1−σ3sin2α−fσ1+σ3+σ1−σ3cos2α.The wing crack will stop propagating whenKΙW≤KIC, 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 macroscopic strain of a rock can be calculated.
(31)ε1∗=8λχcosαE2τ∗cosαπlnla−σ3la−1,ε3∗=χE16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗,where ε1∗ and ε3∗ are the strains along the vertical and horizontal directions, respectively; λ=sinαcosα−fcos2α; γ=−cosαsinα−fsin2α. χ is defined as the initial microcrack density and is expressed as χ=Na2/V. N is the total number of the microcracks in a two-dimensional body of unit thickness whose volume is V=2h×2d.Finally, the total rock strain due to compression considering the normalization condition is(32)ε1=Nc∫aminamaxρada∫0π/2ραsinαdα⋅8λχcosαE×2τ∗cosαπlnla−σ3la−1,ε3=Nc∫aminamaxρada∫0π/2ραsinαdα⋅χE×16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗.


### 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(33)Fσ,κ=β1I1+J2−κ,Gσ,κ=β2I1+J2−κ,where I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, β1=2sinφ/33−sinφ, β2=2sinψ/33−sinψ, and φ and ψ are the friction angle and dilation angle, respectively. κ is the hardening function, which can be expressed as [35]
(34)κ=σ0+a1λexpa2I1−a3λ,where a1, a2, and a3 are the constants which can be obtained by fitting with the uniaxial compressive stress-strain curve. σ0=6ccosφ/33−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:
(35)cm=c0e−0.072m,where c0 is the original cohesion strength (before freeze-thaw cycle) and cm is the cohesion after the mth freeze-thaw cycle.The plastic strain rate is(36)ε̇p=λ̇∂G∂σ,where λ̇ is a proportion coefficient.


### 3.2.5. The Number of the Microcracks in a Rock

The probability function [37] of the microcrack length can be expressed as
(37)ρa=−1acexp−aac,amin≤a≤amax,0,otherwise,where ac is the characteristic length of the microcrack, ac=∫aminamaxρaada. Nc is the total number of the microcracks per unit volume, which is determined by the total volume of the microcrack, Vc=2πbNc∫aminamaxρaa2da.Finally, the curve of probability density function of the microcracks is shown in Figure4.Figure 4
Probability density functions of the microcracks.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 Vc of the microcracks can be expressed with the rock void porosity e and the total volume V of the rock sample, namely, Vc=eV.The first stress intensity factorKI of the microcrack with different lengths at the same stress condition can be expressed as
(38)KI=23τ∗πa.Equation (38) can be changed into
(39)3KIC2τ∗21π=acr,where τ∗=1/2σ1−σ3sin2φ−fσ1+σ3+σ1−σ3cos2φ where acr is the critical length of the microcrack which becomes active at the condition of σ1 and σ3 when KΙ=KIC.As shown in Figure5, the microcrack whose length is greater than acr will propagate.Figure 5
Relationship between the microcrack critical length and axial stress.The total number of microcracks that will be actually activated is given byN=Nc∫aminamaxρada, as shown in Figure 6.Figure 6
Relationship between the microcrack initiation number and axial stress.


### 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, εn+1p, κn+1, rn+1, and hn+1 at tn+1 are determined by integration flow rule and hardening law and Δσtr, σn, κn, rn=∂F/∂σσ=σn, and hn=∂κ/∂λλ=λn at time tn. The main steps are as follows:
(1)
Update strain tensor, and calculate the elastic strain at timetn+1(40)εn+1=εn+Δε,Δσn+1tr=C−1:Δε,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(41)εn+1p=εnp+Δλn+1rn,κn+1=κn∆λn+1hnσn+1=C−1:εn+1−εn+1pFn+1=Fσn+1,κn+1=0where rn=∂G/∂σσ=σn and hn=∂κ/∂κλ=λn.Figure 7
Graphic returns method of semi-implicit integration algorithm.Calculate the plastic internal variable∆λn+1:
(42)an+1k+C−1:∆σn+1k+δλkrn=0,(43)bn+1k−∆κn+1k+δλn+1khn=0,(44)Fn+1k+Fσk:∆σn+1k+Fκk∙∆κn+1k=0,(45)ΔσΔκ=−Aka~k−δλkAkr~n,where
(46)Ak=C00−I,a~k=akbk=00,r~n=rnhn..δλk can be obtained with Equation (45)
(47)δλ=Fk∂F:Ak:r~n,where ∂F=FσFκ.
(48)∆λn+1k+1=Δλnk+δλk.(3)
UpdateNn+1 which is the number of microcracks that begin propagatingAt timetn+1, the updated normal stress is σn+1, under which the critical initiation length acrn+1 of the microcrack is calculated by
(49)acrn+1=1π3KIC2τ∗n+12,where
(50)τ∗n+1=12σ1n+1tr−σ3sin2φ−μσ1n+1tr+σ3+σ1n+1tr−σ3cos2φ,Nn+1=Nc∫acrn+1amaxρada.(4)
Update the additional damage strainεadn+1The additional damage strain matrix ofεad can be expressed as
(51)εad=ε1ε30.Update the length of wing crackln+1(52)KIW=2aτ`n+1cosφπln+1−σ3n+1πln+1,where KIW=KIC.Update the additional strain(53)ε1n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅8λacos2αE×2τ∗n+1cosαπlnln+1a−σ3ln+1a−1,ε3n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅a2E×16τ∗n+1γcos2απlnln+1a+σ3πln+12a2−1−8cosφln+1a−1σ3γ+τ∗n+1.(5)
Update the stress(54)σn+1∗=σn+1−C−1:εadn+1


## 3.2.1. Initial Damage Strain Induced by Freeze-Thaw Cycles

The initial damage strain induced by freeze-thaw cycles is(18)εd=Cdσ,where Cd is the initial damage compliance matrix.Assume the half-length of theith microcrack becomes am after m freeze-thaw cycles. αi is the orientation of the ith microcrack; then, the initial damage compliance matrix due to one single microcrack is given by [25]
(19)C=C0+Ai−1TΔCiAi−1,where [C0] is the elastic compliance matrix, [Ai] is the transformation matrix, and Ai−1 is its inverse matrix.
(20)Ai=cos2αisin2αi−sin2αisin2αicos2αisin2αi12sin2αi−12sin2αicos2αi,Ai−1=cos2αisin2αisin2αisin2αicos2αi−sin2αi12sin2αi−12sin2αicos2αi,ΔCi=0000CniamiKni2dh000CsiamiKsi2dhwhere Cni, Csi, Kni, and Ksi are the compression transferring coefficient, shear transferring coefficient, normal stiffness, and shear stiffness of the ith crack, respectively. d and h are width and height of the rock sample, respectively.The total compliance matrix includingN microcracks whose half-length am is
(21)C=C0+∑i=1NAi−1TΔCiAi−1.The orientation and size of the microcracks in the rock are assumed to be random, which can be expressed with a probability density functionρa,α=ρaρα. According to the distribution law of the microcrack, its orientation and size satisfy the following normalization condition [30]:
(22)∫aminamaxρada∫0π/2ραsinαdα=1.Assume the total number of the microcracks isNc, and then, the number N of the microcracks with half-length of a can be expressed as
(23)N=Nc∫aminaρada.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
(24)C=C0+Nc∫aminamaxAi−1TΔCiAi−1ρada∫0π/2sinαdα.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.


## 3.2.2. Elastic Strain

For a plane stress issue, the elastic constitutive relationship of rock is [32]
(25)εe=C0σ,where εe is the elastic strain matrix, εe=ε11ε33ε13T, σ is the stress matrix, σ=σ11σ33σ13, and C0 is the elastic compliance matrix,
(26)C0=1E−νE0−νE1E0001μ,where E, μ, and υ are the rock elastic modulus, shear modulus, and Poisson ratio, respectively.


## 3.2.3. 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 Figure3.Figure 3
The sliding microcrack model under compression.2a and α are the length and dip angle of the microcrack, respectively. l is the wing crack length. 2h and 2d are the height and width of the model, respectively.Afterm freeze-thaw cycles, the microcrack length becomes am, and for this moment, the first stress intensity factor KI at the microcrack tip is [31, 33]
(27)KI=23τ∗πam,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.
(28)τm=σ1−σ32sin2α,σm=σ1+σ32−σ1−σ32cos2α.WhenKΙ≥KIC (KIC is the rock fracture toughness), the microcrack begins to propagate, and the stress intensity factor KIW at the wing crack tip is [31, 33].
(29)KIW=2aτ∗cosαπl−σ3πl,where
(30)τ∗=12σ1−σ3sin2α−fσ1+σ3+σ1−σ3cos2α.The wing crack will stop propagating whenKΙW≤KIC, 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 macroscopic strain of a rock can be calculated.
(31)ε1∗=8λχcosαE2τ∗cosαπlnla−σ3la−1,ε3∗=χE16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗,where ε1∗ and ε3∗ are the strains along the vertical and horizontal directions, respectively; λ=sinαcosα−fcos2α; γ=−cosαsinα−fsin2α. χ is defined as the initial microcrack density and is expressed as χ=Na2/V. N is the total number of the microcracks in a two-dimensional body of unit thickness whose volume is V=2h×2d.Finally, the total rock strain due to compression considering the normalization condition is(32)ε1=Nc∫aminamaxρada∫0π/2ραsinαdα⋅8λχcosαE×2τ∗cosαπlnla−σ3la−1,ε3=Nc∫aminamaxρada∫0π/2ραsinαdα⋅χE×16τ∗γcos2απlnla+σ3πl2a2−1−8cosαla−1σ3γ+τ∗.


## 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(33)Fσ,κ=β1I1+J2−κ,Gσ,κ=β2I1+J2−κ,where I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, β1=2sinφ/33−sinφ, β2=2sinψ/33−sinψ, and φ and ψ are the friction angle and dilation angle, respectively. κ is the hardening function, which can be expressed as [35]
(34)κ=σ0+a1λexpa2I1−a3λ,where a1, a2, and a3 are the constants which can be obtained by fitting with the uniaxial compressive stress-strain curve. σ0=6ccosφ/33−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:
(35)cm=c0e−0.072m,where c0 is the original cohesion strength (before freeze-thaw cycle) and cm is the cohesion after the mth freeze-thaw cycle.The plastic strain rate is(36)ε̇p=λ̇∂G∂σ,where λ̇ is a proportion coefficient.


## 3.2.5. The Number of the Microcracks in a Rock

The probability function [37] of the microcrack length can be expressed as
(37)ρa=−1acexp−aac,amin≤a≤amax,0,otherwise,where ac is the characteristic length of the microcrack, ac=∫aminamaxρaada. Nc is the total number of the microcracks per unit volume, which is determined by the total volume of the microcrack, Vc=2πbNc∫aminamaxρaa2da.Finally, the curve of probability density function of the microcracks is shown in Figure4.Figure 4
Probability density functions of the microcracks.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 Vc of the microcracks can be expressed with the rock void porosity e and the total volume V of the rock sample, namely, Vc=eV.The first stress intensity factorKI of the microcrack with different lengths at the same stress condition can be expressed as
(38)KI=23τ∗πa.Equation (38) can be changed into
(39)3KIC2τ∗21π=acr,where τ∗=1/2σ1−σ3sin2φ−fσ1+σ3+σ1−σ3cos2φ where acr is the critical length of the microcrack which becomes active at the condition of σ1 and σ3 when KΙ=KIC.As shown in Figure5, the microcrack whose length is greater than acr will propagate.Figure 5
Relationship between the microcrack critical length and axial stress.The total number of microcracks that will be actually activated is given byN=Nc∫aminamaxρada, as shown in Figure 6.Figure 6
Relationship between the microcrack initiation number and axial stress.


## 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, εn+1p, κn+1, rn+1, and hn+1 at tn+1 are determined by integration flow rule and hardening law and Δσtr, σn, κn, rn=∂F/∂σσ=σn, and hn=∂κ/∂λλ=λn at time tn. The main steps are as follows:
(1)
Update strain tensor, and calculate the elastic strain at timetn+1(40)εn+1=εn+Δε,Δσn+1tr=C−1:Δε,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(41)εn+1p=εnp+Δλn+1rn,κn+1=κn∆λn+1hnσn+1=C−1:εn+1−εn+1pFn+1=Fσn+1,κn+1=0where rn=∂G/∂σσ=σn and hn=∂κ/∂κλ=λn.Figure 7
Graphic returns method of semi-implicit integration algorithm.Calculate the plastic internal variable∆λn+1:
(42)an+1k+C−1:∆σn+1k+δλkrn=0,(43)bn+1k−∆κn+1k+δλn+1khn=0,(44)Fn+1k+Fσk:∆σn+1k+Fκk∙∆κn+1k=0,(45)ΔσΔκ=−Aka~k−δλkAkr~n,where
(46)Ak=C00−I,a~k=akbk=00,r~n=rnhn..δλk can be obtained with Equation (45)
(47)δλ=Fk∂F:Ak:r~n,where ∂F=FσFκ.
(48)∆λn+1k+1=Δλnk+δλk.(3)
UpdateNn+1 which is the number of microcracks that begin propagatingAt timetn+1, the updated normal stress is σn+1, under which the critical initiation length acrn+1 of the microcrack is calculated by
(49)acrn+1=1π3KIC2τ∗n+12,where
(50)τ∗n+1=12σ1n+1tr−σ3sin2φ−μσ1n+1tr+σ3+σ1n+1tr−σ3cos2φ,Nn+1=Nc∫acrn+1amaxρada.(4)
Update the additional damage strainεadn+1The additional damage strain matrix ofεad can be expressed as
(51)εad=ε1ε30.Update the length of wing crackln+1(52)KIW=2aτ`n+1cosφπln+1−σ3n+1πln+1,where KIW=KIC.Update the additional strain(53)ε1n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅8λacos2αE×2τ∗n+1cosαπlnln+1a−σ3ln+1a−1,ε3n+1=Nn+1∫acrn+1amaxρada∫0π/2ραsinαdα⋅a2E×16τ∗n+1γcos2απlnln+1a+σ3πln+12a2−1−8cosφln+1a−1σ3γ+τ∗n+1.(5)
Update the stress(54)σn+1∗=σn+1−C−1:εadn+1


## 4. 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.Figure 8
Comparison of the stress-strain relationships for triaxial compression tests with confining pressureσ3=5 MPa by Tan et al. [36] and simulated with the proposed model. The predicted stress-strain relationships agree well with the tested ones.Figure 9
Comparison of the stress-strain relationships for triaxial compression tests with confining pressureσ3=10 MPa by Tan et al. [36] and simulated with the proposed model. The predicted stress-strain relationships agree well with the tested ones.Table 1
The calculation parameters.
E (GPa)vec (MPa)a1 (MPa)a2 (MPa-1)a3ac (m)amin (m)amax (m)36.710.20.006727.435e91e-81001e-41e-70.0012The characteristic lengthac 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 stress-strain 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.


## 5. 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 stress-strain 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



---

*Source: 3177464-2021-07-23.xml*