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

College of Civil Engineering and Architecture, Guangxi University, Nanning, Guangxi 530004, China Key Laboratory of Disaster Prevention and Structural Safety, Guangxi University, Nanning, Guangxi 530004, China School of Civil Engineering, Central South University, Changsha, Hunan 410075, China National Engineering Laboratory for Construction Technology of High Speed Railway, Central South University, Changsha, Hunan 410075, China


Introduction
e 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][2][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]. erefore, 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.
ese 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][14][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][17][18][19][20]. e 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. erefore, F-T damage to the surrounding rock should be considered in the frost heaving force calculation.
To accurately determine the frost heaving force in coldregion 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. e 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.

Calculation Model and Basic Assumptions.
e calculation model of the frost heaving force considering the F-T damage and transversely isotropic characteristics of the surrounding rock is shown in Figure 2(a). e calculation model includes three parts: the lining region (S 1 ), the freeze-thawed surrounding rock (S 2 ), and the unfrozen surrounding rock region (S 3 ). Among them, r 1 and r 2 are the inner and outer radii of the lining, respectively; r 3 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.
As shown in Figure 2, 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 where r 2 is the outer radius of the lining.
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 coordinate z-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 where a ij represents the elastic constants of the transversely isotropic material and , a 26 � 0, a 33 � 1/E z , a 63 � 0, and a 66 � 1/G xy ; E x , E y , and E z represent the elastic moduli in the x-, y-, and z-axis directions, respectively; v xy , v zy , and v zy represent Poisson's ratios; and G xy 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: . e solution to equation (4) is related to the roots of the following equation: 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). erefore, the mapping functions z 1 and z 2 can be calculated by equation (2) if the complex roots µ 1 and µ 2 are known.

Stress and Displacement Components of the Surrounding Rock and Lining.
e 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 (z 1 ) and Ψ 1 (z 2 ) in the unfrozen surrounding rock region (S 3 ) can be expanded by the Laurent series, as follows: where the complex constant a k � a 1k + ia 2k , b k � b 1k + ib 2k , and a 1k , a 2k , b 1k , and b 2k 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: e two complex potential functions Φ 2 (z 1 ) and Ψ 2 (z 2 ) in the freeze-thawed surrounding rock region (S 2 ) can be expanded by the Laurent series, as follows: where the complex constants e 2k , f 1k , and f 2k are undetermined coefficients. Similarly, two complex potential functions Φ 3 (z) and Ψ 3 (z) in the lining region (S 1 ) can be expanded by the Laurent series, as follows: where the complex constants g k � g 1k + ig 2k , h k � h 1k + ih 2k , m k � m 1k + im 2k , and n k � n 1k + in 2k and g 1k , g 2k , h 1k , h 2k , m 1k , m 2k , n 1k , and n 2k 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: where Re is the real part of a complex number and Im is the imaginary part of a complex number. e 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: 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: 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]. p 1 � β 11 μ 2 1 + 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: 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 freezethawed surrounding rock region.

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. e 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 . erefore, when frost heaving of rock occurs, the displacement boundary conditions and the stress boundary conditions at interfaces ρ � ρ 1 , ρ � ρ 2 , and ρ � ρ 3 are as follows.

Stress Boundary Conditions
Advances in Civil Engineering 5 (18c)

Displacement Boundary Conditions
It should be noted that the expressions B * ω 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: 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. e 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. en, 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. e weight of the surrounding rock is 22.0 kN/m 3 , 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%.

A Case Study of the Frost Heaving Force considering the F-T Damage and
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. e main treatment methods for circular sections are as follows: 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. e 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 r 1 � 5.97 m, r 2 � 6.27 m, and r 3 � 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.

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 Figure 3) was conducted to obtain the elastic parameters of surrounding rock [28]. e basic elastic parameter test results for different numbers of F-T cycles are listed in Table 1.
where E 1 and μ 1 are the elastic modulus and Poisson's ratio in the direction parallel to the bedding surface, respectively, and E 2 , μ 2 , and G 2 are the elastic modulus, Poisson's ratio, and shear modulus in the direction perpendicular to the bedding surface, respectively. Combined with the elastic parameter test results in Table 1, 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.
en, 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.

Determination of the Mechanical Parameters of the Lining.
e 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. e specific conversion equation is as follows: where E is the elastic modulus of the converted shotcrete, E c is the elastic modulus of the original shotcrete, E s is the elastic modulus of the steel arch, A s 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.

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). en, 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. e following conclusions can be obtained: (1) Figure 4(a) is the calculation model of the frost heaving force acting on the lining assuming that the surrounding rock is homogeneous isotropic continuous media. e 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      When the surrounding rocks are not layered, the frost heaving force acting on the lining is distributed symmetrically. e 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) Figure 4 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. erefore, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered. Figure 5 shows the circumferential stress on the outer boundary of the lining for different bedding angles of the surrounding rock. e 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)).

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. e 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 Figure 6, 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. e 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]. is cushions the lining structure's stress and results in a gradual reduction in the frost heaving force in cold-region tunnels.

Discussion
e 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. e frost heaving force is one of the most critical challenges that can cause serious frost damage in cold-region tunnels. erefore, 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. e 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. e 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. e 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]. erefore, 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 coldregion 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.

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. e following conclusions can be drawn based on the limited research work reported in this paper:  0°to 90°when the number of F-T cycles is constant. e 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) e 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. erefore, the transversely isotropic characteristics of surrounding rock have a considerable impact on the frost heaving force and should be considered.
(3) e 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.

Data Availability
All data are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.