# An Analytical Solution for the Frost Heaving Force considering the Freeze-Thaw Damage and Transversely Isotropic Characteristics of the Surrounding Rock in Cold-Region Tunnels

Authors: Jiabing Zhang; Xiaohu Zhang; Helin Fu; Yimin Wu; Zhen Huang; Deyue Chen
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Frost damage is a frequent occurrence in cold regions and can threaten the normal use and structural stability of tunnel engineering projects. To accurately determine the frost heaving force and effectively evaluate the frost damage in cold-region tunnels, an analytical solution for the frost heaving force considering the freeze-thaw (F-T) damage and transversely isotropic characteristics of surrounding rock is presented based on complex variable theory and the power series method. The calculation results indicate that the frost heaving force acts on the lining considering that the transversely isotropic characteristics of surrounding rock are significantly greater than those when assuming the surrounding rock is homogeneous isotropic media. This demonstrates that the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered. The frost heaving force continuously increases as the bedding angle increases from 0° to 90°, and the maximum frost heaving force in the Guanjiao tunnel (the rock mass bedding angle is 30°) of the Xining-Geermu Railway in China is approximately 1.04 MPa. In addition, the influence of F-T cycles on the frost heaving force in cold-region tunnels is investigated based on the analytical solution of the frost heaving force proposed in this paper. The frost heaving force acting on the lining decreases with an increasing number of F-T cycles due to the deterioration of the mechanical parameters of the surrounding rock.


---

## Body

## 1. Introduction

The frost heaving of rock is a critical problem that exists in many geotechnical engineering applications, such as tunnel and mining engineering in cold regions [1–3]. With the increasing investment in the construction of infrastructure in the Qinghai-Tibet Plateau in China, tunnel and underground project construction in cold regions has progressively increased [4, 5]. However, most tunnels during the operation period suffer from serious frost damage, as shown in Figure 1, thereby severely affecting the normal use and structural stability of tunnel engineering projects and even causing accidents [6, 7]. When the underground water of tunnel in cold region freezes into ice, the surrounding rock expands in volume, and the deformation pressure is frost heaving force, which is determined by the amount of frost heaving deformation, the stiffness of surrounding rock, and the stiffness of the lining structure. Under this action, frost damage phenomena such as lining cracking will occur [8]. Therefore, the effective calculation of the frost heaving force is significant for engineering applications, helping to provide technical support and guarantee structural safety and scientifically based construction in cold-region tunnels.Figure 1
Frost damage in cold-region tunnels.Many studies have been performed to obtain a reasonable frost heaving force in cold-region tunnels [9–12]. These studies demonstrated that the frost heaving force is primarily induced by the 9% volume expansion from the water-ice phase change in the surrounding rock [13–15]. In addition, a number of theoretical models have been developed in recent decades to predict the frost heaving forces in cold-region tunnels [16–20]. The frost heave of the surrounding rock is assumed to be isotropic in all of the above analytical solutions; however, very few publications are available where the frost heave of the surrounding rock is assumed to be transversely isotropic. In reality, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force [21].In addition, many multifield coupled frost heave models have been proposed to effectively evaluate the frost damage in cold-region tunnels [22, 23]. Lai et al. [24] developed a hydrothermal-mechanical model to describe the freezing process of silty soil. Furthermore, Zhang et al. [25] advanced a hydrothermal-salt-mechanical model in freezing saline soil. However, to the best of the authors’ knowledge, F-T damage to rock masses is also a critical problem during tunnel engineering projects [26, 27], and the frost heaving force acting on the lining changes constantly after the recurrent action of freezing and thawing. Therefore, F-T damage to the surrounding rock should be considered in the frost heaving force calculation.To accurately determine the frost heaving force in cold-region tunnels, an analytical solution for the frost heaving force acting on the lining that considers the F-T damage and the transversely isotropic characteristics of the surrounding rock has been first presented based on complex variable theory and the power series method. The mechanical parameters of surrounding rock with different F-T cycles were then determined by a series of experiments [28]. Finally, the effects of the number of F-T cycles on the frost heaving force in cold-region tunnels were discussed in detail based on the analytical solution of the frost heaving force proposed in this paper.


## 2. Analytical Solution for the Frost Heaving Force Acting on the Lining considering the F-T Damage and Transversely Isotropic Characteristics of the Surrounding Rock

### 2.1. Calculation Model and Basic Assumptions

The calculation model of the frost heaving force considering the F-T damage and transversely isotropic characteristics of the surrounding rock is shown in Figure2(a). The calculation model includes three parts: the lining region (S1), the freeze-thawed surrounding rock (S2), and the unfrozen surrounding rock region (S3). Among them, r1 and r2 are the inner and outer radii of the lining, respectively; r3 is the radius of the freeze-thawed surrounding rock; β is the bedding angle of the surrounding rock; and σx∞, σy∞, and τxy∞ are the stress components acting at infinity. To obtain the analytical solution of the frost heaving force in cold-region tunnels, the following assumptions are introduced: (1) the calculation model for the frost heaving force can be approximated as an infinite plane hole problem, which is not affected by the boundary effect, and the ground stress change gradient can be ignored; (2) the freeze-thawed surrounding rock and unfrozen surrounding rock are heterogeneous with transversely isotropic properties, while the lining is a homogeneous material with isotropic and elastic mechanical properties; (3) the calculation model for the frost heaving force can be regarded as a plane strain problem; and (4) the cross section of the tunnel is equivalent to a circle, and the frost heaving and F-T damage to the lining are not considered.Figure 2
Calculation model of the frost heaving force considering the transversely isotropic characteristics of surrounding rock. (a) Calculation model in the physical plane (z-plane). (b) Calculation model in the image plane (ζ-plane).As shown in Figure2, according to complex variable theory, the mapping function that transforms the outer region of the circular hole in the z-plane into the outer region of the unit circle in the ζ-plane can be given as(1)z=ωζ,=r2ζ,where r2 is the outer radius of the lining.


### 2.2. Mapping Function for the Elastic Mechanic Analysis of Transversely Isotropic Rocks

In anisotropic elastic mechanics, the cross section no longer remains flat after deformation due to the presence of positive shear and shear coupling [29]. To solve this problem, three mapping functions z = ω(ζ), z1 = ω1(ζ1), and z2 = ω2(ζ2) are introduced to transform the outer regions of the circular hole in the z-, z1-, and z2-planes into the outer regions of the unit circle in the ζ-, ζ1-, and ζ2-planes, respectively. The mathematical relationships between the three mapping functions z, z1, and z2 can be given as(2)z1=x+μ1y=γ1z+δ1z¯,z2=x+μ2y=γ2z+δ2z¯,where γ1=1−iμ1/2, γ2=1−iμ2/2, δ1=1+iμ1/2, and δ2=1+iμ2/2.(3)μ1=α1+iβ1,μ2=α2+iβ2,where the constants α1, α2, β1, and β2 can be determined according to the elastic constants of the anisotropic rock mass.For the transversely isotropic plane strain problem, when the coordinatez-axis is perpendicular to one plane of elastic symmetry, the compatibility equation for Airy’s stress function without considering body forces can be given as(4)β11∂4U∂y4−2β16∂4U∂x∂y3+2β12+β66∂4U∂x2∂y2−2β26∂4U∂x3∂y+β22∂4U∂x4=0,(5)βij=aij−ai3aj3a33,i,j=1,2,4,5,6,where aij represents the elastic constants of the transversely isotropic material and a11=1/Ex, a12=−vyx/Ey, a13=−vzx/Ez, a16=0, a22=1/Ey, a23=−vzy/Ez, a26=0, a33=1/Ez, a63=0, and a66=1/Gxy; Ex, Ey, and Ez represent the elastic moduli in the x-, y-, and z-axis directions, respectively; vxy, vzy, and vzy represent Poisson’s ratios; and Gxy represents the shear modulus of the transversely isotropic material. For the transversely isotropic plane strain problem, as shown in Figure 2, the values of the elastic moduli, Poisson’s ratios, and shear moduli should satisfy the following relations: Ey=Ez, vyx=vxy, vzx=vxz, vzy=vyz, Gxy=Gxz, and Gyz=Ey/21+vyz.The solution to equation (4) is related to the roots of the following equation:(6)β11μ4−2β16μ3+2β12+β66μ2−2β26μ+β22=0,where the four roots in equation (6) must be conjugate complex roots (i.e. μ1, μ1¯, μ2, and μ2¯). Here, µ1 and µ2 are exactly the same as those expressed in equation (3). Therefore, the mapping functions z1 and z2 can be calculated by equation (2) if the complex roots µ1 and µ2 are known.Substituting equation (1) into equation (2), we can obtain the following equations:(7)z1=ω1ζ1=γ1ωζ+δ1ωζ¯=γ1r2ζ−1+δ1r2ζ¯−1,z2=ω2ζ2=γ2ωζ+δ2ωζ¯=γ2r2ζ−1+δ2r2ζ¯−1,where ωζ is an analytical function of ζ. However, its conjugate, ωζ¯, is not an analytical function of ζ [29]. Therefore, the mapping functions z1 and z2 in equation (7) are not analytical functions of ζ. To determine the relationships between ζ1, ζ2, and ζ and to determine ω1(ζ1) and ω2(ζ2), we assume that ζ1 = ζ2 = ζ = σ = eiα along the unit circle, which corresponds to the boundary of the hole. Substituting this assumption into equation (7) and replacing σ with ζ1 and ζ2 in the expressions ω1(ζ1) and ω2(ζ2), respectively, we can obtain the following relationships between ζ1, ζ2, and ζ:(8)γ1ζ1+δ1ζ1=γ1ζ−1+δ1ζ¯−1,γ2ζ2+δ2ζ2=γ2ζ−1+δ2ζ¯−1.Therefore,ζ1 and ζ2 can be determined by equation (8) after the value of ζ is determined, and then the three mapping functions ω(ζ), ω1(ζ1), and ω2(ζ2) can be obtained.


### 2.3. Stress and Displacement Components of the Surrounding Rock and Lining

The complex function method is used to determine the stress and displacement components of the surrounding rock and lining, and the two complex potential functions Φ1(z1) and Ψ1(z2) in the unfrozen surrounding rock region (S3) can be expanded by the Laurent series, as follows:(9)Φ1z1=Φ1ω1ζ1=Φ1∗ζ1=B∗ω1ζ1+∑k=0∞akζ1−k,Ψ1z2=Ψ1ω2ζ2=Ψ1∗ζ2=B′∗+iC′∗ω2ζ2+∑k=0∞bkζ2−k,where the complex constant ak=a1k+ia2k, bk=b1k+ib2k, and a1k, a2k, b1k, and b2k are undetermined coefficients. B∗, B′∗, and C′∗ can be determined according to the stress components acting at infinity (i.e., σx∞, σy∞, and τxy∞), as follows:(10)B∗=σx∞+α22+β22σy∞+2α2τxy∞2α2−α12+β22−β12,B′∗=α12−β12−2α1α2σy∞−σx∞−2α2τxy∞2α2−α12+β22−β12,C′∗=α1−α2σx∞+α12−β12−α22−β22τxy∞2β2α2−α12+β22−β12+α2α12−β12−α1α22−β22σy∞2β2α2−α12+β22−β12.The two complex potential functions Φ2(z1) and Ψ2(z2) in the freeze-thawed surrounding rock region (S2) can be expanded by the Laurent series, as follows:(11)Φ2z1=Φ2ω1ζ1=Φ2∗ζ1=∑k=0∞ckζ1k+∑k=1∞dkζ1−k,Ψ2z2=Ψ2ω2ζ2=Ψ2∗ζ2=∑k=0∞ekζ2k+∑k=1∞fkζ2−k,where the complex constants ck=c1k+ic2k, dk=d1k+id2k, ek=e1k+ie2k, and fk=f1k+if2k and c1k, c2k, d1k, d2k, e1k, e2k, f1k, and f2k are undetermined coefficients.Similarly, two complex potential functions Φ3(z) and Ψ3(z) in the lining region (S1) can be expanded by the Laurent series, as follows:(12)Φ3z=Φ3ωζ=Φ3∗ζ=∑k=0∞gkζk+∑k=1∞hkζ−k,Ψ3z=Ψ3ωζ=Ψ3∗ζ=∑k=0∞mkζk+∑k=1∞nkζ−k,where the complex constants gk=g1k+ig2k, hk=h1k+ih2k, mk=m1k+im2k, and nk=n1k+in2k and g1k, g2k, h1k, h2k, m1k, m2k, n1k, and n2k are undetermined coefficients.When the lining is a homogeneous isotropic material, the stress components at any point on the z-plane can be given by the complex potential functions, as follows:(13)σx=2ReΦ3′z−Rez¯Φ′′3z+Ψ3′z,σy=2ReΦ3′z+Rez¯Φ′′3z+Ψ3′z,τxy=Imz¯Φ′′3z+Ψ3′z,where Re is the real part of a complex number and Im is the imaginary part of a complex number.The freeze-thawed surrounding rock and unfrozen surrounding rock are heterogeneous with transversely isotropic materials; hence, the stress components at any point on the z-plane in the freeze-thawed and unfrozen surrounding rock regions can be determined by the complex potential functions, as follows:(14)σx=2Reμ12Φi′z1+μ22Ψi′z2,σy=2ReΦi′z1+Ψi′z2,τxy=−2Reμ1Φi′z1+μ2Ψi′z2,i=1,2.Assuming that the volume of the lining and unfrozen surrounding rock are not affected by the change in the surrounding rock temperature, the displacement components of any point on the z-plane in the lining and unfrozen surrounding rock region can be determined by the complex potential functions, as follows:(15)2Gu+iv=κΦ3z−zΦ3′z¯−Ψ3z¯,u=2Rep1Φ1z1+p2Ψ1z2,v=2Req1Φ1z1+q2Ψ1z2,where G is the shear modulus of the lining. Since the longitudinal length of the tunnel is much larger than the diameter of the section, this problem can be regarded as a plane strain problem, and then κ=3−4ν. p1=β11μ12+β12−β16μ1, p2=β11μ22+β12−β16μ2, q1=β12μ1+β22/μ1−β26, and q2=β12μ2+β22/μ2−β26.Since the volume of the freeze-thawed surrounding rock region is affected by temperature changes, the displacement components at any point on the z-plane in the freeze-thawed surrounding rock region cannot satisfy the expression of a complex variable function. After mathematical deduction, the displacement equation considering the transversely isotropic characteristics of surrounding rock can be obtained as follows:(16)u+iv=2Rep1Φ2z1+p2Ψ2z2+2Req1Φ2z1+q2Ψ2z2i+εz,(17)ε=ν1+2νθ,where ε represents the linear expansion coefficient of surrounding rock in the freeze-thawed surrounding rock region on the isotropic plane y′oz′, ν represents Poisson’s ratio perpendicular to the isotropic plane, and θ represents the volume change rate of surrounding rock in the freeze-thawed surrounding rock region.


### 2.4. Boundary Conditions and Solution of Frost Heaving Force

In terms of the stress and displacement solutions of elastic mechanics, boundary conditions are usually divided into two categories: displacement boundary conditions and stress boundary conditions. The complex potential function must satisfy the continuous boundary conditions of displacement and stress at interfacesρ = ρ2 and ρ = ρ3 and the stress boundary conditions at interface ρ = ρ1. Therefore, when frost heaving of rock occurs, the displacement boundary conditions and the stress boundary conditions at interfaces ρ = ρ1, ρ = ρ2, and ρ = ρ3 are as follows.


#### 2.4.1. Stress Boundary Conditions

(18a)Φ3∗ρ1σ+ωρ1σω′ρ1σ¯Φ3∗′ρ1σ¯+Ψ3∗ρ1σ¯=0,ζ=ζ1=ζ2=ρ1σ,(18b)2ReΦ2∗σ+Ψ2∗σ=ReΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,2Reμ1Φ2∗σ+μ2Ψ2∗σ=ImΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,ζ=ζ1=ζ2=σ,(18c)ReΦ1∗ρ3σ+Ψ1∗ρ3σ=ReΦ2∗ρ3σ+Ψ2∗ρ3σ,Reμ1Φ1∗ρ3σ+μ2Ψ1∗ρ3σ=Reμ1Φ2∗ρ3σ+μ2Ψ2∗ρ3σ,ζ=ζ1=ζ2=ρ3σ.


#### 2.4.2. Displacement Boundary Conditions

(19a)2GB1+iB2+εωσ=κΦ3∗σ−ωσω′σ¯Φ3∗′σ¯−Ψ3∗σ¯,B1=2Rep1Φ2∗σ+p2Ψ2∗σ,B2=2Req1Φ2∗σ+q2Ψ2∗σ,ζ=ζ1=ζ2=σ,(19b)C1+iC2=B3+iB4+εωρ3σ,B3=2Rep1Φ2∗ρ3σ+p2Ψ2∗ρ3σ,B4=2Req1Φ2∗ρ3σ+q2Ψ2∗ρ3σ,C1=2Rep1Φ1∗0ρ3σ+p2Ψ1∗0ρ3σ,C2=2Req1Φ1∗0ρ3σ+q2Ψ1∗0ρ3σ,ζ=ζ1=ζ2=ρ3σ.It should be noted that the expressionsB∗ω1ζ1 and B′∗+iC′∗ω2ζ2 in equation (4) represent the formation displacement component caused by the initial ground stress of the surrounding rock, and thus these expressions should be removed when calculating the formation displacement caused by tunnel excavation in equation (19b). In addition, it is assumed that there is no relative slip between the contact surfaces in equations (17) and (18). From the mapping function (1), we can obtain the following:(20)ωζω′ζ¯=−1ζ3.Based on the power series solution method proposed by Chen [30] and through the combination of the stress boundary conditions and displacement boundary conditions, infinite linear equations can be obtained. The coefficients in the complex potential functions can be obtained by solving the infinite linear equations, and the stress distributions of the surrounding rock and lining can be obtained by equations (13) and (14), respectively. Then, the stress components in the rectangular coordinate system are transformed into the polar coordinate system by equation (21) to obtain the radial stress on the outer boundary of the lining. Finally, the analytical solution of the frost heaving force in cold-region tunnels, which considers the transversely isotropic characteristics of surrounding rock, can be obtained by analysing the radial stress on the outer boundary of the lining before and after frost heaving.(21)σρ=σx+σy2+σx−σy2cos2α+τxysin2α,σα=σx+σy2+σy−σx2cos2α−τxysin2α,τρα=σy−σx2sin2α+τxycos2α.


## 2.1. Calculation Model and Basic Assumptions

The calculation model of the frost heaving force considering the F-T damage and transversely isotropic characteristics of the surrounding rock is shown in Figure2(a). The calculation model includes three parts: the lining region (S1), the freeze-thawed surrounding rock (S2), and the unfrozen surrounding rock region (S3). Among them, r1 and r2 are the inner and outer radii of the lining, respectively; r3 is the radius of the freeze-thawed surrounding rock; β is the bedding angle of the surrounding rock; and σx∞, σy∞, and τxy∞ are the stress components acting at infinity. To obtain the analytical solution of the frost heaving force in cold-region tunnels, the following assumptions are introduced: (1) the calculation model for the frost heaving force can be approximated as an infinite plane hole problem, which is not affected by the boundary effect, and the ground stress change gradient can be ignored; (2) the freeze-thawed surrounding rock and unfrozen surrounding rock are heterogeneous with transversely isotropic properties, while the lining is a homogeneous material with isotropic and elastic mechanical properties; (3) the calculation model for the frost heaving force can be regarded as a plane strain problem; and (4) the cross section of the tunnel is equivalent to a circle, and the frost heaving and F-T damage to the lining are not considered.Figure 2
Calculation model of the frost heaving force considering the transversely isotropic characteristics of surrounding rock. (a) Calculation model in the physical plane (z-plane). (b) Calculation model in the image plane (ζ-plane).As shown in Figure2, according to complex variable theory, the mapping function that transforms the outer region of the circular hole in the z-plane into the outer region of the unit circle in the ζ-plane can be given as(1)z=ωζ,=r2ζ,where r2 is the outer radius of the lining.


## 2.2. Mapping Function for the Elastic Mechanic Analysis of Transversely Isotropic Rocks

In anisotropic elastic mechanics, the cross section no longer remains flat after deformation due to the presence of positive shear and shear coupling [29]. To solve this problem, three mapping functions z = ω(ζ), z1 = ω1(ζ1), and z2 = ω2(ζ2) are introduced to transform the outer regions of the circular hole in the z-, z1-, and z2-planes into the outer regions of the unit circle in the ζ-, ζ1-, and ζ2-planes, respectively. The mathematical relationships between the three mapping functions z, z1, and z2 can be given as(2)z1=x+μ1y=γ1z+δ1z¯,z2=x+μ2y=γ2z+δ2z¯,where γ1=1−iμ1/2, γ2=1−iμ2/2, δ1=1+iμ1/2, and δ2=1+iμ2/2.(3)μ1=α1+iβ1,μ2=α2+iβ2,where the constants α1, α2, β1, and β2 can be determined according to the elastic constants of the anisotropic rock mass.For the transversely isotropic plane strain problem, when the coordinatez-axis is perpendicular to one plane of elastic symmetry, the compatibility equation for Airy’s stress function without considering body forces can be given as(4)β11∂4U∂y4−2β16∂4U∂x∂y3+2β12+β66∂4U∂x2∂y2−2β26∂4U∂x3∂y+β22∂4U∂x4=0,(5)βij=aij−ai3aj3a33,i,j=1,2,4,5,6,where aij represents the elastic constants of the transversely isotropic material and a11=1/Ex, a12=−vyx/Ey, a13=−vzx/Ez, a16=0, a22=1/Ey, a23=−vzy/Ez, a26=0, a33=1/Ez, a63=0, and a66=1/Gxy; Ex, Ey, and Ez represent the elastic moduli in the x-, y-, and z-axis directions, respectively; vxy, vzy, and vzy represent Poisson’s ratios; and Gxy represents the shear modulus of the transversely isotropic material. For the transversely isotropic plane strain problem, as shown in Figure 2, the values of the elastic moduli, Poisson’s ratios, and shear moduli should satisfy the following relations: Ey=Ez, vyx=vxy, vzx=vxz, vzy=vyz, Gxy=Gxz, and Gyz=Ey/21+vyz.The solution to equation (4) is related to the roots of the following equation:(6)β11μ4−2β16μ3+2β12+β66μ2−2β26μ+β22=0,where the four roots in equation (6) must be conjugate complex roots (i.e. μ1, μ1¯, μ2, and μ2¯). Here, µ1 and µ2 are exactly the same as those expressed in equation (3). Therefore, the mapping functions z1 and z2 can be calculated by equation (2) if the complex roots µ1 and µ2 are known.Substituting equation (1) into equation (2), we can obtain the following equations:(7)z1=ω1ζ1=γ1ωζ+δ1ωζ¯=γ1r2ζ−1+δ1r2ζ¯−1,z2=ω2ζ2=γ2ωζ+δ2ωζ¯=γ2r2ζ−1+δ2r2ζ¯−1,where ωζ is an analytical function of ζ. However, its conjugate, ωζ¯, is not an analytical function of ζ [29]. Therefore, the mapping functions z1 and z2 in equation (7) are not analytical functions of ζ. To determine the relationships between ζ1, ζ2, and ζ and to determine ω1(ζ1) and ω2(ζ2), we assume that ζ1 = ζ2 = ζ = σ = eiα along the unit circle, which corresponds to the boundary of the hole. Substituting this assumption into equation (7) and replacing σ with ζ1 and ζ2 in the expressions ω1(ζ1) and ω2(ζ2), respectively, we can obtain the following relationships between ζ1, ζ2, and ζ:(8)γ1ζ1+δ1ζ1=γ1ζ−1+δ1ζ¯−1,γ2ζ2+δ2ζ2=γ2ζ−1+δ2ζ¯−1.Therefore,ζ1 and ζ2 can be determined by equation (8) after the value of ζ is determined, and then the three mapping functions ω(ζ), ω1(ζ1), and ω2(ζ2) can be obtained.


## 2.3. Stress and Displacement Components of the Surrounding Rock and Lining

The complex function method is used to determine the stress and displacement components of the surrounding rock and lining, and the two complex potential functions Φ1(z1) and Ψ1(z2) in the unfrozen surrounding rock region (S3) can be expanded by the Laurent series, as follows:(9)Φ1z1=Φ1ω1ζ1=Φ1∗ζ1=B∗ω1ζ1+∑k=0∞akζ1−k,Ψ1z2=Ψ1ω2ζ2=Ψ1∗ζ2=B′∗+iC′∗ω2ζ2+∑k=0∞bkζ2−k,where the complex constant ak=a1k+ia2k, bk=b1k+ib2k, and a1k, a2k, b1k, and b2k are undetermined coefficients. B∗, B′∗, and C′∗ can be determined according to the stress components acting at infinity (i.e., σx∞, σy∞, and τxy∞), as follows:(10)B∗=σx∞+α22+β22σy∞+2α2τxy∞2α2−α12+β22−β12,B′∗=α12−β12−2α1α2σy∞−σx∞−2α2τxy∞2α2−α12+β22−β12,C′∗=α1−α2σx∞+α12−β12−α22−β22τxy∞2β2α2−α12+β22−β12+α2α12−β12−α1α22−β22σy∞2β2α2−α12+β22−β12.The two complex potential functions Φ2(z1) and Ψ2(z2) in the freeze-thawed surrounding rock region (S2) can be expanded by the Laurent series, as follows:(11)Φ2z1=Φ2ω1ζ1=Φ2∗ζ1=∑k=0∞ckζ1k+∑k=1∞dkζ1−k,Ψ2z2=Ψ2ω2ζ2=Ψ2∗ζ2=∑k=0∞ekζ2k+∑k=1∞fkζ2−k,where the complex constants ck=c1k+ic2k, dk=d1k+id2k, ek=e1k+ie2k, and fk=f1k+if2k and c1k, c2k, d1k, d2k, e1k, e2k, f1k, and f2k are undetermined coefficients.Similarly, two complex potential functions Φ3(z) and Ψ3(z) in the lining region (S1) can be expanded by the Laurent series, as follows:(12)Φ3z=Φ3ωζ=Φ3∗ζ=∑k=0∞gkζk+∑k=1∞hkζ−k,Ψ3z=Ψ3ωζ=Ψ3∗ζ=∑k=0∞mkζk+∑k=1∞nkζ−k,where the complex constants gk=g1k+ig2k, hk=h1k+ih2k, mk=m1k+im2k, and nk=n1k+in2k and g1k, g2k, h1k, h2k, m1k, m2k, n1k, and n2k are undetermined coefficients.When the lining is a homogeneous isotropic material, the stress components at any point on the z-plane can be given by the complex potential functions, as follows:(13)σx=2ReΦ3′z−Rez¯Φ′′3z+Ψ3′z,σy=2ReΦ3′z+Rez¯Φ′′3z+Ψ3′z,τxy=Imz¯Φ′′3z+Ψ3′z,where Re is the real part of a complex number and Im is the imaginary part of a complex number.The freeze-thawed surrounding rock and unfrozen surrounding rock are heterogeneous with transversely isotropic materials; hence, the stress components at any point on the z-plane in the freeze-thawed and unfrozen surrounding rock regions can be determined by the complex potential functions, as follows:(14)σx=2Reμ12Φi′z1+μ22Ψi′z2,σy=2ReΦi′z1+Ψi′z2,τxy=−2Reμ1Φi′z1+μ2Ψi′z2,i=1,2.Assuming that the volume of the lining and unfrozen surrounding rock are not affected by the change in the surrounding rock temperature, the displacement components of any point on the z-plane in the lining and unfrozen surrounding rock region can be determined by the complex potential functions, as follows:(15)2Gu+iv=κΦ3z−zΦ3′z¯−Ψ3z¯,u=2Rep1Φ1z1+p2Ψ1z2,v=2Req1Φ1z1+q2Ψ1z2,where G is the shear modulus of the lining. Since the longitudinal length of the tunnel is much larger than the diameter of the section, this problem can be regarded as a plane strain problem, and then κ=3−4ν. p1=β11μ12+β12−β16μ1, p2=β11μ22+β12−β16μ2, q1=β12μ1+β22/μ1−β26, and q2=β12μ2+β22/μ2−β26.Since the volume of the freeze-thawed surrounding rock region is affected by temperature changes, the displacement components at any point on the z-plane in the freeze-thawed surrounding rock region cannot satisfy the expression of a complex variable function. After mathematical deduction, the displacement equation considering the transversely isotropic characteristics of surrounding rock can be obtained as follows:(16)u+iv=2Rep1Φ2z1+p2Ψ2z2+2Req1Φ2z1+q2Ψ2z2i+εz,(17)ε=ν1+2νθ,where ε represents the linear expansion coefficient of surrounding rock in the freeze-thawed surrounding rock region on the isotropic plane y′oz′, ν represents Poisson’s ratio perpendicular to the isotropic plane, and θ represents the volume change rate of surrounding rock in the freeze-thawed surrounding rock region.


## 2.4. Boundary Conditions and Solution of Frost Heaving Force

In terms of the stress and displacement solutions of elastic mechanics, boundary conditions are usually divided into two categories: displacement boundary conditions and stress boundary conditions. The complex potential function must satisfy the continuous boundary conditions of displacement and stress at interfacesρ = ρ2 and ρ = ρ3 and the stress boundary conditions at interface ρ = ρ1. Therefore, when frost heaving of rock occurs, the displacement boundary conditions and the stress boundary conditions at interfaces ρ = ρ1, ρ = ρ2, and ρ = ρ3 are as follows.


### 2.4.1. Stress Boundary Conditions

(18a)Φ3∗ρ1σ+ωρ1σω′ρ1σ¯Φ3∗′ρ1σ¯+Ψ3∗ρ1σ¯=0,ζ=ζ1=ζ2=ρ1σ,(18b)2ReΦ2∗σ+Ψ2∗σ=ReΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,2Reμ1Φ2∗σ+μ2Ψ2∗σ=ImΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,ζ=ζ1=ζ2=σ,(18c)ReΦ1∗ρ3σ+Ψ1∗ρ3σ=ReΦ2∗ρ3σ+Ψ2∗ρ3σ,Reμ1Φ1∗ρ3σ+μ2Ψ1∗ρ3σ=Reμ1Φ2∗ρ3σ+μ2Ψ2∗ρ3σ,ζ=ζ1=ζ2=ρ3σ.


### 2.4.2. Displacement Boundary Conditions

(19a)2GB1+iB2+εωσ=κΦ3∗σ−ωσω′σ¯Φ3∗′σ¯−Ψ3∗σ¯,B1=2Rep1Φ2∗σ+p2Ψ2∗σ,B2=2Req1Φ2∗σ+q2Ψ2∗σ,ζ=ζ1=ζ2=σ,(19b)C1+iC2=B3+iB4+εωρ3σ,B3=2Rep1Φ2∗ρ3σ+p2Ψ2∗ρ3σ,B4=2Req1Φ2∗ρ3σ+q2Ψ2∗ρ3σ,C1=2Rep1Φ1∗0ρ3σ+p2Ψ1∗0ρ3σ,C2=2Req1Φ1∗0ρ3σ+q2Ψ1∗0ρ3σ,ζ=ζ1=ζ2=ρ3σ.It should be noted that the expressionsB∗ω1ζ1 and B′∗+iC′∗ω2ζ2 in equation (4) represent the formation displacement component caused by the initial ground stress of the surrounding rock, and thus these expressions should be removed when calculating the formation displacement caused by tunnel excavation in equation (19b). In addition, it is assumed that there is no relative slip between the contact surfaces in equations (17) and (18). From the mapping function (1), we can obtain the following:(20)ωζω′ζ¯=−1ζ3.Based on the power series solution method proposed by Chen [30] and through the combination of the stress boundary conditions and displacement boundary conditions, infinite linear equations can be obtained. The coefficients in the complex potential functions can be obtained by solving the infinite linear equations, and the stress distributions of the surrounding rock and lining can be obtained by equations (13) and (14), respectively. Then, the stress components in the rectangular coordinate system are transformed into the polar coordinate system by equation (21) to obtain the radial stress on the outer boundary of the lining. Finally, the analytical solution of the frost heaving force in cold-region tunnels, which considers the transversely isotropic characteristics of surrounding rock, can be obtained by analysing the radial stress on the outer boundary of the lining before and after frost heaving.(21)σρ=σx+σy2+σx−σy2cos2α+τxysin2α,σα=σx+σy2+σy−σx2cos2α−τxysin2α,τρα=σy−σx2sin2α+τxycos2α.


## 2.4.1. Stress Boundary Conditions

(18a)Φ3∗ρ1σ+ωρ1σω′ρ1σ¯Φ3∗′ρ1σ¯+Ψ3∗ρ1σ¯=0,ζ=ζ1=ζ2=ρ1σ,(18b)2ReΦ2∗σ+Ψ2∗σ=ReΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,2Reμ1Φ2∗σ+μ2Ψ2∗σ=ImΦ3∗σ+ωσω′σ¯Φ3∗′σ¯+Ψ3∗σ¯,ζ=ζ1=ζ2=σ,(18c)ReΦ1∗ρ3σ+Ψ1∗ρ3σ=ReΦ2∗ρ3σ+Ψ2∗ρ3σ,Reμ1Φ1∗ρ3σ+μ2Ψ1∗ρ3σ=Reμ1Φ2∗ρ3σ+μ2Ψ2∗ρ3σ,ζ=ζ1=ζ2=ρ3σ.


## 2.4.2. Displacement Boundary Conditions

(19a)2GB1+iB2+εωσ=κΦ3∗σ−ωσω′σ¯Φ3∗′σ¯−Ψ3∗σ¯,B1=2Rep1Φ2∗σ+p2Ψ2∗σ,B2=2Req1Φ2∗σ+q2Ψ2∗σ,ζ=ζ1=ζ2=σ,(19b)C1+iC2=B3+iB4+εωρ3σ,B3=2Rep1Φ2∗ρ3σ+p2Ψ2∗ρ3σ,B4=2Req1Φ2∗ρ3σ+q2Ψ2∗ρ3σ,C1=2Rep1Φ1∗0ρ3σ+p2Ψ1∗0ρ3σ,C2=2Req1Φ1∗0ρ3σ+q2Ψ1∗0ρ3σ,ζ=ζ1=ζ2=ρ3σ.It should be noted that the expressionsB∗ω1ζ1 and B′∗+iC′∗ω2ζ2 in equation (4) represent the formation displacement component caused by the initial ground stress of the surrounding rock, and thus these expressions should be removed when calculating the formation displacement caused by tunnel excavation in equation (19b). In addition, it is assumed that there is no relative slip between the contact surfaces in equations (17) and (18). From the mapping function (1), we can obtain the following:(20)ωζω′ζ¯=−1ζ3.Based on the power series solution method proposed by Chen [30] and through the combination of the stress boundary conditions and displacement boundary conditions, infinite linear equations can be obtained. The coefficients in the complex potential functions can be obtained by solving the infinite linear equations, and the stress distributions of the surrounding rock and lining can be obtained by equations (13) and (14), respectively. Then, the stress components in the rectangular coordinate system are transformed into the polar coordinate system by equation (21) to obtain the radial stress on the outer boundary of the lining. Finally, the analytical solution of the frost heaving force in cold-region tunnels, which considers the transversely isotropic characteristics of surrounding rock, can be obtained by analysing the radial stress on the outer boundary of the lining before and after frost heaving.(21)σρ=σx+σy2+σx−σy2cos2α+τxysin2α,σα=σx+σy2+σy−σx2cos2α−τxysin2α,τρα=σy−σx2sin2α+τxycos2α.


## 3. A Case Study of the Frost Heaving Force considering the F-T Damage and Transversely Isotropic Characteristics of the Surrounding Rock

### 3.1. Determination of the Calculation Model

To study the frost heaving force considering the F-T damage and transversely isotropic characteristics of the surrounding rock, section DyK304 + 370 of the Guanjiao tunnel of the Xining-Geermu Railway in China is taken as a case study. The Guanjiao tunnel is 32.605 km in length and 200 m in depth. The weight of the surrounding rock is 22.0 kN/m3, the bedding angle is 30°, the lateral pressure coefficient is 1.0, and the volume change rate of the surrounding rock under the condition of frost heaving is 9%.In the above theoretical analysis of the frost heaving force, the tunnel section is treated as a circular section, while the actual excavated section of the Guanjiao tunnel is a curved wall section. At present, a noncircular section is usually treated as a circular section by using mathematical transformation. The main treatment methods for circular sections are as follows:(22)r=4h2+b24costan−1b/2h,where r is the radius of the tunnel after equal circle treatment, b is the section span of the original tunnel, and h is the section height of the original tunnel. The radius of the Guanjiao tunnel after equal circular treatment is approximately 6.27 m according to equation (22). Combined with the initial support thickness of 30 cm and the maximum freezing depth of 2.99 m in the Guanjiao tunnel, the radii r1 = 5.97 m, r2 = 6.27 m, and r3 = 9.26 m can be calculated for the z-plane in Section 2.1. Furthermore, the radii ρ1 = 5.97 m, ρ2 = 6.27 m, and ρ3 = 9.26 m can be obtained by equation (1) for the ζ-plane.


### 3.2. Determination of the Mechanical Parameters of Surrounding Rock

Before the calculation of the frost heaving force acting on the lining considering the F-T damage and transversely isotropic characteristics of the surrounding rock, the basic mechanical parameters of the surrounding rock should be determined. In general, the determination of the mechanical parameters of the surrounding rock is relatively complicated, and the methods used are mainly test methods, engineering analogy methods, inverse analysis methods, etc. In this study, a series of uniaxial compression experiments with different numbers of F-T cycles (the temperature control curve for each F-T cycle is shown in Figure3) was conducted to obtain the elastic parameters of surrounding rock [28]. The basic elastic parameter test results for different numbers of F-T cycles are listed in Table 1.(23)μ1=μ1′cosβ−sinβcosβ+μ1′sinβ,μ2=μ2′cosβ−sinβcosβ+μ2′sinβ,where E1 and μ1 are the elastic modulus and Poisson’s ratio in the direction parallel to the bedding surface, respectively, and E2, μ2, and G2 are the elastic modulus, Poisson’s ratio, and shear modulus in the direction perpendicular to the bedding surface, respectively.Figure 3
F-T cycle experimental procedures (adapted from [4]). (a) Characteristics of temperature changes at the location of the project. (b) TDS-300 automatic F-T cycling device. (c) Temperature control curve for each F-T.Table 1
Basic elastic parameters of surrounding rock with different numbers of F-T cycles.
Number of F-T cyclesE1 (GPa)E2 (GPa)ν1ν2G2 (GPa)024.6717.250.2750.2038.4525.3418.370.2670.2119.21521.6515.140.3180.2347.2320.9715.830.3090.2427.461019.4613.440.3380.2606.3118.9213.730.3470.2716.091517.9912.210.3610.2775.8317.2712.070.3750.2845.492016.4011.400.3820.2885.2415.9711.860.3940.2995.083014.8710.250.4540.3084.6315.4910.870.4480.3164.58Combined with the elastic parameter test results in Table1, the elastic constants of the surrounding rock, a′ij, can be obtained for different numbers of F-T cycles in the local coordinate system, as shown in Table 2. Then, the elastic parameters of the surrounding rock are substituted into equations (4) and (5), and the four complex roots of the compatible equation in the local coordinate system under different numbers of F-T cycles can be determined, as shown in Table 3. Finally, the four complex roots of the compatible equation in the global coordinate system under different numbers of F-T cycles can be obtained according to the principle of coordinate transformation (as shown in equation (23)), and the constants α1, α2, β1, and β2 in equation (3) can be determined.Table 2
Elastic constants of surrounding rock,a′ij, for different numbers of F-T cycles.
Number of F-T cyclesa11′a12′a13′a22′a23′a33′a66′00.056−0.0083−0.00830.040−0.01080.0400.113350.065−0.0112−0.01120.047−0.01470.0470.1361100.074−0.0138−0.01380.052−0.01780.0520.1613150.082−0.0159−0.01590.057−0.02090.0570.1767200.086−0.0188−0.01880.062−0.02410.0620.1938300.095−0.0229−0.02290.066−0.02980.0660.2172Table 3
Complex roots of the compatible equation with different numbers of F-T cycles.
Number of F-T cyclesμ1′μ1′¯μ2′μ2′¯00.8087i−0.8087i1.0219i−1.0219i50.7921i−0.7921i1.0413i−1.0413i100.7253i−0.7253i1.1142i−1.1142i150.7349i−0.7349i1.0852i−1.0852i200.7298i−0.7298i1.1093i−1.1093i300.6957i−0.6957i1.1164i−1.1164i


### 3.3. Determination of the Mechanical Parameters of the Lining

The initial support considers the action of only the I-beam support and C25 shotcrete. Among them, the steel arch support parameters are equivalently converted by improving the elastic modulus of the shotcrete. The specific conversion equation is as follows:(24)E=Ec+EsAshd,where E is the elastic modulus of the converted shotcrete, Ec is the elastic modulus of the original shotcrete, Es is the elastic modulus of the steel arch, As is the sectional area of the steel arch, h is the thickness of the original shotcrete, and d is the spacing between steel arches. According to equation (24), the tunnel’s lining elasticity modulus is 34.27 GPa, Poisson’s ratio is 0.25, and the shear modulus is 13.71 GPa.


### 3.4. Calculation and Analysis of the Frost Heaving Force in the Guanjiao Tunnel

With the zero freezing and thawing cycles as an example, the analytical solution of the stress distribution in a transversely isotropic rock mass tunnel can first be obtained according to the boundary conditions in equations (18) and (19). Then, the stress components in the rectangular coordinate system are transformed to a polar coordinate system by equation (21), and the radial stress on the outer boundary of the lining can be obtained. Finally, the frost heaving force considering the transversely isotropic frost heave of surrounding rock can be obtained by comparing and analysing the calculated results for analytical unfrozen heaving and analytical frozen heaving, as shown in Figure 4. The following conclusions can be obtained:(1)
Figure4(a) is the calculation model of the frost heaving force acting on the lining assuming that the surrounding rock is homogeneous isotropic continuous media. The radial stress on the outer boundary of the lining considering the transversely isotropic frost heave of the surrounding rock is significantly greater than that when assuming the surrounding rock is homogeneous isotropic media, and then it gradually increases with increasing bedding orientation of the surrounding rock.(2)
Based on the complex variable theory and power series solution method, the frost heaving force considering the transversely isotropic frost heave of the surrounding rock is calculated to be between 0.21 MPa and 3.59 MPa (0 F-T cycles), and the frost heaving force continuously increases as the bedding angle increases from 0° to 90°. When the bedding orientation is 90°, the frost heaving force reaches a maximum. In addition, the rock mass bedding angle of the Guanjiao tunnel is 30°, and the maximum frost heaving force acting on the lining is approximately 1.04 MPa.(3)
The positions of the maximum frost heaving force acting on the lining differ because the mechanical properties of transversely isotropic rocks exhibit distinct anisotropy and directional dependence. When the surrounding rocks are not layered, the frost heaving force acting on the lining is distributed symmetrically. The frost heaving force is the highest at the arch waist (θ = 0° or 180°) of the tunnel and the lowest at the arch top (θ = 90°) or arch bottom (θ = 270°) of the tunnel (as shown in Figure 4(a)). When the rock mass bedding angle is 0°–15°, the maximum frost heaving force occurs near the arch waist of the tunnel (θ = 0° or 180°) (as shown in Figures 4(b)-4(c)). When the rock mass bedding angle is 45°–60°, the maximum frost heaving force occurs near the right arch shoulder (θ = 40°) or the left arch foot (θ = 220°) (as shown in Figures 4(e)-4(f)). When the rock mass bedding angle is 75°–90°, the maximum frost heaving force occurs near the arch top (θ = 90°) or the arch bottom (θ = 270°) (as shown in Figures 4(g)-4(h)).(4)
Figure4 shows that the frost heaving force acting on the lining considering the transversely isotropic characteristics of surrounding rock is significantly greater than that assuming that the surrounding rock is homogeneous isotropic continuous media. It can be concluded that the anisotropy of the rock mass has a significant influence on the frost heaving force. Additionally, the lining structures of transversely isotropic rock mass tunnels are more prone to frost damage, which makes the maintenance of tunnels very difficult and poses a serious threat to the structure and operational safety of the tunnels. Therefore, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered.Figure 4
Radial stresses on the outer boundary of the lining for different bedding angles of the surrounding rock. (a) Homogeneous isotropic rock mass. (b) Bedding angles of 0°. (c) Bedding angles of 15°. (d) Bedding angles of 30°. (e) Bedding angles of 45°. (f) Bedding angles of 60°. (g) Bedding angles of 75°. (h) Bedding angles of 90°.
(a)(b)(c)(d)(e)(f)(g)(h)Figure5 shows the circumferential stress on the outer boundary of the lining for different bedding angles of the surrounding rock. The circumferential stress on the outer boundary of the lining considering the transversely isotropic characteristics of the surrounding rock is significantly smaller than that when assuming that the surrounding rock is homogeneous isotropic media, and then it gradually decreases with increasing bedding orientation of the surrounding rock. Additionally, the circumferential stress on the outer boundary of the lining significantly increases because of the influence of frost heaving on the rock. When the surrounding rock is not layered, the maximum increment of the outer boundary circumferential stress of the lining is approximately 3.48 MPa (as shown in Figure 5(a)). Moreover, the increment of the outer boundary circumferential stress of the lining decreases with increasing bedding orientation of the surrounding rock (as shown in Figures 5(b)–5(h)).Figure 5
Circumferential stresses on the outer boundary of the lining for different bedding angles of the surrounding rock. (a) Homogeneous isotropic rock mass. (b) Bedding angles of 0°. (c) Bedding angles of 15°. (d) Bedding angles of 30°. (e) Bedding angles of 45°. (f) Bedding angles of 60°. (g) Bedding angles of 75°. (h) Bedding angles of 90°.
(a)(b)(c)(d)(e)(f)(g)(h)


## 3.1. Determination of the Calculation Model

To study the frost heaving force considering the F-T damage and transversely isotropic characteristics of the surrounding rock, section DyK304 + 370 of the Guanjiao tunnel of the Xining-Geermu Railway in China is taken as a case study. The Guanjiao tunnel is 32.605 km in length and 200 m in depth. The weight of the surrounding rock is 22.0 kN/m3, the bedding angle is 30°, the lateral pressure coefficient is 1.0, and the volume change rate of the surrounding rock under the condition of frost heaving is 9%.In the above theoretical analysis of the frost heaving force, the tunnel section is treated as a circular section, while the actual excavated section of the Guanjiao tunnel is a curved wall section. At present, a noncircular section is usually treated as a circular section by using mathematical transformation. The main treatment methods for circular sections are as follows:(22)r=4h2+b24costan−1b/2h,where r is the radius of the tunnel after equal circle treatment, b is the section span of the original tunnel, and h is the section height of the original tunnel. The radius of the Guanjiao tunnel after equal circular treatment is approximately 6.27 m according to equation (22). Combined with the initial support thickness of 30 cm and the maximum freezing depth of 2.99 m in the Guanjiao tunnel, the radii r1 = 5.97 m, r2 = 6.27 m, and r3 = 9.26 m can be calculated for the z-plane in Section 2.1. Furthermore, the radii ρ1 = 5.97 m, ρ2 = 6.27 m, and ρ3 = 9.26 m can be obtained by equation (1) for the ζ-plane.


## 3.2. Determination of the Mechanical Parameters of Surrounding Rock

Before the calculation of the frost heaving force acting on the lining considering the F-T damage and transversely isotropic characteristics of the surrounding rock, the basic mechanical parameters of the surrounding rock should be determined. In general, the determination of the mechanical parameters of the surrounding rock is relatively complicated, and the methods used are mainly test methods, engineering analogy methods, inverse analysis methods, etc. In this study, a series of uniaxial compression experiments with different numbers of F-T cycles (the temperature control curve for each F-T cycle is shown in Figure3) was conducted to obtain the elastic parameters of surrounding rock [28]. The basic elastic parameter test results for different numbers of F-T cycles are listed in Table 1.(23)μ1=μ1′cosβ−sinβcosβ+μ1′sinβ,μ2=μ2′cosβ−sinβcosβ+μ2′sinβ,where E1 and μ1 are the elastic modulus and Poisson’s ratio in the direction parallel to the bedding surface, respectively, and E2, μ2, and G2 are the elastic modulus, Poisson’s ratio, and shear modulus in the direction perpendicular to the bedding surface, respectively.Figure 3
F-T cycle experimental procedures (adapted from [4]). (a) Characteristics of temperature changes at the location of the project. (b) TDS-300 automatic F-T cycling device. (c) Temperature control curve for each F-T.Table 1
Basic elastic parameters of surrounding rock with different numbers of F-T cycles.
Number of F-T cyclesE1 (GPa)E2 (GPa)ν1ν2G2 (GPa)024.6717.250.2750.2038.4525.3418.370.2670.2119.21521.6515.140.3180.2347.2320.9715.830.3090.2427.461019.4613.440.3380.2606.3118.9213.730.3470.2716.091517.9912.210.3610.2775.8317.2712.070.3750.2845.492016.4011.400.3820.2885.2415.9711.860.3940.2995.083014.8710.250.4540.3084.6315.4910.870.4480.3164.58Combined with the elastic parameter test results in Table1, the elastic constants of the surrounding rock, a′ij, can be obtained for different numbers of F-T cycles in the local coordinate system, as shown in Table 2. Then, the elastic parameters of the surrounding rock are substituted into equations (4) and (5), and the four complex roots of the compatible equation in the local coordinate system under different numbers of F-T cycles can be determined, as shown in Table 3. Finally, the four complex roots of the compatible equation in the global coordinate system under different numbers of F-T cycles can be obtained according to the principle of coordinate transformation (as shown in equation (23)), and the constants α1, α2, β1, and β2 in equation (3) can be determined.Table 2
Elastic constants of surrounding rock,a′ij, for different numbers of F-T cycles.
Number of F-T cyclesa11′a12′a13′a22′a23′a33′a66′00.056−0.0083−0.00830.040−0.01080.0400.113350.065−0.0112−0.01120.047−0.01470.0470.1361100.074−0.0138−0.01380.052−0.01780.0520.1613150.082−0.0159−0.01590.057−0.02090.0570.1767200.086−0.0188−0.01880.062−0.02410.0620.1938300.095−0.0229−0.02290.066−0.02980.0660.2172Table 3
Complex roots of the compatible equation with different numbers of F-T cycles.
Number of F-T cyclesμ1′μ1′¯μ2′μ2′¯00.8087i−0.8087i1.0219i−1.0219i50.7921i−0.7921i1.0413i−1.0413i100.7253i−0.7253i1.1142i−1.1142i150.7349i−0.7349i1.0852i−1.0852i200.7298i−0.7298i1.1093i−1.1093i300.6957i−0.6957i1.1164i−1.1164i


## 3.3. Determination of the Mechanical Parameters of the Lining

The initial support considers the action of only the I-beam support and C25 shotcrete. Among them, the steel arch support parameters are equivalently converted by improving the elastic modulus of the shotcrete. The specific conversion equation is as follows:(24)E=Ec+EsAshd,where E is the elastic modulus of the converted shotcrete, Ec is the elastic modulus of the original shotcrete, Es is the elastic modulus of the steel arch, As is the sectional area of the steel arch, h is the thickness of the original shotcrete, and d is the spacing between steel arches. According to equation (24), the tunnel’s lining elasticity modulus is 34.27 GPa, Poisson’s ratio is 0.25, and the shear modulus is 13.71 GPa.


## 3.4. Calculation and Analysis of the Frost Heaving Force in the Guanjiao Tunnel

With the zero freezing and thawing cycles as an example, the analytical solution of the stress distribution in a transversely isotropic rock mass tunnel can first be obtained according to the boundary conditions in equations (18) and (19). Then, the stress components in the rectangular coordinate system are transformed to a polar coordinate system by equation (21), and the radial stress on the outer boundary of the lining can be obtained. Finally, the frost heaving force considering the transversely isotropic frost heave of surrounding rock can be obtained by comparing and analysing the calculated results for analytical unfrozen heaving and analytical frozen heaving, as shown in Figure 4. The following conclusions can be obtained:(1)
Figure4(a) is the calculation model of the frost heaving force acting on the lining assuming that the surrounding rock is homogeneous isotropic continuous media. The radial stress on the outer boundary of the lining considering the transversely isotropic frost heave of the surrounding rock is significantly greater than that when assuming the surrounding rock is homogeneous isotropic media, and then it gradually increases with increasing bedding orientation of the surrounding rock.(2)
Based on the complex variable theory and power series solution method, the frost heaving force considering the transversely isotropic frost heave of the surrounding rock is calculated to be between 0.21 MPa and 3.59 MPa (0 F-T cycles), and the frost heaving force continuously increases as the bedding angle increases from 0° to 90°. When the bedding orientation is 90°, the frost heaving force reaches a maximum. In addition, the rock mass bedding angle of the Guanjiao tunnel is 30°, and the maximum frost heaving force acting on the lining is approximately 1.04 MPa.(3)
The positions of the maximum frost heaving force acting on the lining differ because the mechanical properties of transversely isotropic rocks exhibit distinct anisotropy and directional dependence. When the surrounding rocks are not layered, the frost heaving force acting on the lining is distributed symmetrically. The frost heaving force is the highest at the arch waist (θ = 0° or 180°) of the tunnel and the lowest at the arch top (θ = 90°) or arch bottom (θ = 270°) of the tunnel (as shown in Figure 4(a)). When the rock mass bedding angle is 0°–15°, the maximum frost heaving force occurs near the arch waist of the tunnel (θ = 0° or 180°) (as shown in Figures 4(b)-4(c)). When the rock mass bedding angle is 45°–60°, the maximum frost heaving force occurs near the right arch shoulder (θ = 40°) or the left arch foot (θ = 220°) (as shown in Figures 4(e)-4(f)). When the rock mass bedding angle is 75°–90°, the maximum frost heaving force occurs near the arch top (θ = 90°) or the arch bottom (θ = 270°) (as shown in Figures 4(g)-4(h)).(4)
Figure4 shows that the frost heaving force acting on the lining considering the transversely isotropic characteristics of surrounding rock is significantly greater than that assuming that the surrounding rock is homogeneous isotropic continuous media. It can be concluded that the anisotropy of the rock mass has a significant influence on the frost heaving force. Additionally, the lining structures of transversely isotropic rock mass tunnels are more prone to frost damage, which makes the maintenance of tunnels very difficult and poses a serious threat to the structure and operational safety of the tunnels. Therefore, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered.Figure 4
Radial stresses on the outer boundary of the lining for different bedding angles of the surrounding rock. (a) Homogeneous isotropic rock mass. (b) Bedding angles of 0°. (c) Bedding angles of 15°. (d) Bedding angles of 30°. (e) Bedding angles of 45°. (f) Bedding angles of 60°. (g) Bedding angles of 75°. (h) Bedding angles of 90°.
(a)(b)(c)(d)(e)(f)(g)(h)Figure5 shows the circumferential stress on the outer boundary of the lining for different bedding angles of the surrounding rock. The circumferential stress on the outer boundary of the lining considering the transversely isotropic characteristics of the surrounding rock is significantly smaller than that when assuming that the surrounding rock is homogeneous isotropic media, and then it gradually decreases with increasing bedding orientation of the surrounding rock. Additionally, the circumferential stress on the outer boundary of the lining significantly increases because of the influence of frost heaving on the rock. When the surrounding rock is not layered, the maximum increment of the outer boundary circumferential stress of the lining is approximately 3.48 MPa (as shown in Figure 5(a)). Moreover, the increment of the outer boundary circumferential stress of the lining decreases with increasing bedding orientation of the surrounding rock (as shown in Figures 5(b)–5(h)).Figure 5
Circumferential stresses on the outer boundary of the lining for different bedding angles of the surrounding rock. (a) Homogeneous isotropic rock mass. (b) Bedding angles of 0°. (c) Bedding angles of 15°. (d) Bedding angles of 30°. (e) Bedding angles of 45°. (f) Bedding angles of 60°. (g) Bedding angles of 75°. (h) Bedding angles of 90°.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4. Effect of F-T Cycles on the Frost Heaving Force in Cold-Region Tunnels

To comprehensively understand the effect of F-T cycles on the frost heaving force in cold-region tunnels, different mechanical parameters were used for the surrounding rock after different numbers of F-T cycles were applied. The maximum frost heaving force acting on the lining after F-T cycles can be obtained based on the analytical solution for the frost heaving force proposed in this paper. As revealed in Figure6, the maximum frost heaving force acting on the lining decreases with an increasing number of F-T cycles when the bedding angle of the surrounding rock is constant. For example, when the bedding angle is 30°, the maximum frost heaving forces acting on the lining are 1.04 MPa, 0.98 MPa, 0.94 MPa, 0.91 MPa, 0.89 MPa, and 0.85 MPa for 0, 5, 10, 15, 20, and 30 F-T cycles, respectively. The mechanism for this phenomenon is that the deformation of rock is more likely to occur after the action of external forces due to the rock deterioration gradually increasing as the number of F-T cycles increases [13]. This cushions the lining structure’s stress and results in a gradual reduction in the frost heaving force in cold-region tunnels.Figure 6
Maximum frost heaving force in cold-region tunnels after different numbers of F-T cycles. (a) After 0 F-T cycles. (b) After 5 F-T cycles. (c) After 10 F-T cycles. (d) After 15 F-T cycles. (e) After 20 F-T cycles. (f) After 30 F-T cycles.
(a)(b)(c)(d)(e)(f)


## 5. Discussion

The seasonal frozen and permafrost regions are widely distributed throughout China, accounting for 75.8% of the total land area. Along with the implementation of the strategy of “the Belt and Road,” a large number of tunnels are under construction in cold regions. The frost heaving force is one of the most critical challenges that can cause serious frost damage in cold-region tunnels. Therefore, a reasonable calculation of the frost heaving force that considers the F-T damage and the transversely isotropic characteristics of the surrounding rock is expected to provide references for the design and construction of cold-region tunnels.At present, the commonly used calculation models of the frost heaving force are mainly divided into three categories. The first model assumes that the surrounding rock forms a freezing-thawing circle and the pore water or fissure water in the surrounding rock of the freezing-thawing circle will freeze and expand under a negative temperature condition, thus generating a frost heaving force. The second model suggests that there is a water accumulation space between the initial support and the secondary lining of the tunnel and the pore or fissure water in the surrounding rock will continuously supply to the water accumulation space. Under a negative temperature condition, the water will freeze and expand in volume, and the existence of the primary support and secondary lining restricts such volume expansion, thus forming the frost heaving force. The third model holds that the surrounding rock of any tunnel all has a weathering layer and the water in the weathering layer will freeze and expand under the condition of negative temperature, thus causing the frost heaving force. However, due to the complexity of the tunnel and underground engineering system, the existing calculation models of the frost heaving force are all studied with isotropic surrounding rocks as a precondition. In fact, as the tunnel construction gradually moves to cold regions, the surrounding rock of the tunnel will be in a state of alternating positive and negative temperatures for a long time, leading to a stronger anisotropy of the rock mass. In this case, the anisotropy of the rock mass cannot be ignored [4]. Therefore, an analytical solution for the frost heaving force in cold-region tunnels considering the transversely isotropic characteristics of the surrounding rock, which is presented in this paper, may improve the design techniques and safe operation of a transversely isotropic rock mass tunnel in cold regions.In addition, the tunnel section is treated as an equal circle in the calculation of the frost heaving force in cold-region tunnels, and the model corresponding to the actual tunnel section is not used for the calculation. Due to the difference in the models, these research results are somewhat different from a realistic situation, so it is still necessary to further study the problem of the frost heaving force considering an actual section of tunnel in future studies.


## 6. Summary and Conclusions

To accurately determine the frost heaving force acting on the lining in cold-region tunnels, an elastic analytical solution of the frost heaving force in cold-region tunnels that considers the F-T damage and transversely isotropic characteristics of the surrounding rock has been developed. The following conclusions can be drawn based on the limited research work reported in this paper:(1)
The frost heaving force acting on the lining continuously increases as the bedding angle increases from 0° to 90° when the number of F-T cycles is constant. The rock mass bedding angle of the Guanjiao tunnel is 30°, and the maximum frost heaving force acting on the lining is approximately 1.04 MPa.(2)
The frost heaving force acting on the lining considering the transversely isotropic characteristics of surrounding rock is significantly greater than that when assuming the surrounding rock is homogeneous isotropic continuous media. Therefore, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered.(3)
The positions of the maximum frost heaving force acting on the lining differ because the mechanical properties of transversely isotropic rocks exhibit distinct anisotropy. Additionally, the frost heaving force acting on the lining decreases with an increasing number of F-T cycles due to the deterioration of the mechanical parameters of the surrounding rock.



---

*Source: 6654778-2020-12-03.xml*