Semi-analytical solution of frost heave force for circular tunnels in strain-softening rock mass

The freezing method is often adopted to reinforce the surrounding formation during deep foundation pit excavations. According to the traditional theory of freezing wall design in coal mine construction engineering, it is necessary to design a continuous, thicker freezing wall. The multi-turn combination cylinder freezing wall structure can optimize the freezing design, but the strength and stability of the wall structure are theoretical and technical issues that urgently need to be researched and solved. In this work, the mechanical characteristics of frozen walls of multi-coil composite cylinders are analyzed using theoretical methods. Based on the elastoplastic theory, the freezing wall of a multi-turn combined cylinder is simplified to an axisymmetric plane problem, and the theoretical calculation model of the freezing expansion force, generated by the freezing process, around the perimeter of a circular hole under an axisymmetric plane strain model is developed. The corresponding variable model for the displacement immobility point of the freezing circle is proposed, and the progressive formation of the freezing wall is achieved through iterative calculations. The theoretical calculation formulas for the freezing force generated by the freezing process around a circular hole under the plane strain model are derived considering the gradual formation process of freezing expansion. The results of the research provide the theoretical basis for the construction of permafrost deep foundation pit engineering.

1 Introduction

With the country’s implementation of the Western development strategy, more and more tunnels have been built in complex mountainous areas (Fu et al., 2024; Ji et al., 2024; Zhang T. et al., 2024), especially in cold regions. Unfortunately, low temperatures, which may be lower than −40°C, bring a series of problems for geotechnical engineering, such as frost heaving. This significantly affects the physical and mechanical behaviors of the surrounding rock and lining structure. To ensure safety, the frost heaving acting on the lining structure should be determined first.

Numerous studies on tunnel excavation in cold regions have been carried out in the last few years. The results show that the main cause for frost heave is water migration and volume expansion of water-rich geotechnical materials. Some scholars have proposed analytical solutions for circular openings using the Laplace transform method. Ji et al. (2017) investigated the frost heaving-induced pressure of frost-susceptible silty clay under various restraints and thermal gradients. The results indicate that the frost heaving-induced pressure decreases rapidly with the relaxation of restraints, while this downward trend slows down as the frost heaving ratio increases.

The constitutive equation directly determines the displacement direction induced by frost heaving behavior, which is the main factor of frost heave force. In a previous study, the surrounding rock was considered to undergo isotropic expansion during the freezing process, where the elastic modulus of frozen rock was assumed to deteriorate. However, this is obviously different from the actual conditions. Researchers assumed the average radius of the frozen region as zero displacement position for circular tunnels. Hence, the frost heaving displacement point inside the tunnel and that of the outer part move toward infinity. Based on the finite element method, Feng et al. (2017) proposed a supplemental equation relating the frost heaving force, frost heaving displacement, and the inner and outer radii. Considering the predominant frost heaving along the heat conduction direction, Zhang et al. (2023) derived analytical solutions for the frost heaving force for circular tunnels, in which isotropic expansion is only a particular case.

Moreover, rock-like materials generally manifest progressive failure characteristics after the peak strength is reached during the loading process, known as strain-softening behavior (Shi et al., 2023a; Shi et al., 2023b; Shi et al., 2023c; Zhang et al., 2023; Zhang W. et al., 2024; Gu et al., 2024). Considering the strain-softening behavior, Lee and Pietruszczal (2008), Wang et al. (2010), and Zhang et al. (2012) studied the stress and displacement of circular tunnels using numerical and analytical methods. These studies show that the volume plastic strain gradually increases due to the formation of macro cracks (Jing et al., 2020; Wu et al., 2020; Wu et al., 2022; Wu et al., 2024; Yuan et al., 2024), so the porosity of the surrounding rock should depend on both its initial porosity and the volume plastic strain. Moreover, the expansion of frost-affected rock mass is induced by the volume expansion of water during its transition from liquid to solid. In other words, the linear coefficient of linear expansion should be considered a function of volume plastic strain. On the other hand, the freezing of the surrounding rock is a progressive process. When the constitutive model is related to the geometry of the frozen region, the frost heaving force and displacement would show significant differences with the assumption of the freezing region forming simultaneously. However, existing research studies have demonstrated a relative lack of focus on this particular area. In this paper, based on the elastoplastic theory, the progressive formation process of frost heave is considered. The theoretical formula for the freezing force generated by the freezing process around a circular hole under the plane strain model is derived to provide guidance for the design and construction of freezing walls in ultra-large diameter deep foundation pit projects.

2 Problem description

Figure 1 shows a circular tunnel of radius R₀ excavated in an infinite isotropic rock mass subject to a hydrostatic stress p₀. The linear structure, with a thickness of t, is immediately applied after excavation with a deformation allowance of u₀. Once the supporting pressure provided by the linear structure is lower than the critical value p_c, a plastic region with a radius R_p is developed around the tunnel, and then, the water fills the cracks of the plastic rock mass. The frost rock mass is progressively formed from the excavation surface to the inner during heat conduction. The radius of the frost region is R_f, which can be determined using the field temperature test of the surrounding rock. In this study, the frost behavior is assumed after the excavation effect, and the frost heave does not change the failure state of the surrounding rock, so the total stress is composed of excavation force and frost heave stress.

Figure 1

Figure 1. Calculation mechanical model.

2.1 Governing equations

For the axisymmetric plane strain problem, the stress should satisfy the following equilibrium equation without considering the body force:

$\frac{\mathrm{d}{\sigma }_r}{\mathrm{d}r}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=0.(1)$

The failure state of the surrounding rock is governed by the yield criterion. The nonlinear Hoek–Brown (HB) criterion, which takes the influence of the geological strength index (GSI) into account, is widely used for geotechnical engineering. In the polar coordinate, the generalized HB criterion in the ith annulus can be expressed as

${\sigma }_{\theta }={\sigma }_r+H\left({\sigma }_r,{\gamma }^p\right),(2)$

where $H\left({\sigma }_r,{\gamma }^p\right)={\sigma }_c{\left(m{\sigma }_r/{\sigma }_c+s\right)}^a$, σ_c is the uniaxial compression strength of the intact rock, and m and s are strength parameters.

The strain–displacement compatibility equation can be expressed as follows:

${\epsilon }_r=\frac{\mathrm{d}u}{\mathrm{d}r},{\epsilon }_{\theta }=\frac{u}{r}.(3)$

For the elastoplastic problem, the total strain of plastic rock mass is composed of the elastic strain and plastic strain, i.e.,

$\left\{\begin{array}{l}{\epsilon }_r={\epsilon }_r^e+{\epsilon }_r^p\\ {\epsilon }_{\theta }={\epsilon }_{\theta }^e+{\epsilon }_{\theta }^p\end{array}\right..(4)$

The elastic strain should satisfy Hooke’s law, and it can be expressed as Equation 5:

$\left[\begin{array}{l}{\epsilon }_r^e\\ {\epsilon }_{\theta }^e\end{array}\right]=\frac{1}{2G}\left[\begin{array}{cc}1-\nu & -\nu \\ -\nu & 1-\nu \end{array}\right]\left[\begin{array}{l}{\sigma }_r\\ {\sigma }_{\theta }\end{array}\right].(5)$

For the geotechnical material, the M-C plastic potential function is commonly used. The relation between the radial and circumferential plastic strain can be obtained as Equation 6:

${\epsilon }_r^p=-\beta {\epsilon }_{\theta }^p,(6)$

where β = (1 + sinψ)/(1-sinψ) and ψ is the dilation angle.

Theoretically speaking, the frost heave is mainly induced by the decrease in the inner radius of the frost region. When a free frost heave happens for a circular annulus, the axial length increases and the inner radius decreases, so there exists a zero displacement radius $\overline{R}$, which is located between the inner and outer radii. Assuming the coefficient of linear expansion ξ, the displacement induced by free frost of the ith annulus can be expressed as $\xi \left(\overline{R}-r\right)$, and the displacement due to the axial constraint effect can be expressed as $-\mu \xi r$, so the total radial displacement induced by free frost heave under plane strain conditions can be expressed as $\xi \overline{R}-\left(1+\mu \right)\xi r$, and the radial strain induced by frost heave can be obtained using the geometric equation from Equation 3 as $-\left(1+\mu \right)\xi$. Additionally, additional radial and circumferential forces also occur due to the expansion of frost-affected rock mass and the constraint of the linear structure and unfrozen rock, so the total radial and circumferential strains induced by free heave can be expressed as Equation 7:

$\left\{\begin{array}{l}{\epsilon }_r^f=\frac{1}{2G^f}\left[\left(1-{\nu }^f\right){\sigma }_r^f-{\nu }^f{\sigma }_{\theta }^f\right]-\left(1+{\nu }^f\right)\xi \\ {\epsilon }_{\theta }^f=\frac{1}{2G^f}\left[\left(1-{\nu }^f\right){\sigma }_{\theta }^f-{\nu }^f{\sigma }_r^f\right]+\xi \overline{R}/r-\left(1+{\nu }^f\right)\xi \end{array}\right.,(7)$

where ${\sigma }_r^f$ and ${\sigma }_{\theta }^f$ are the radial stress and circumferential stress induced by frost heave, respectively.

2.2 Evolutional law of parameters

With respect to the strain-softening rock mass, the plastic shear strain ${\gamma }^p$ is widely employed as the softening index to describe the evolution law of parameters. The real evolutional law of strength and deformation parameters can be determined using a series of triaxial compression tests. The parameter value can be obtained in a piecewise manner. For the sake of simplicity, a linear evolutional law of strength and deformation parameters is assumed, and it can be expressed as Equation 8:

$\omega =\left\{\begin{array}{ll}{\omega }_p-\left({\omega }_p-{\omega }_r\right)\frac{{\gamma }^p}{{\gamma }^{p*}}& 0<{\gamma }^p<{\gamma }^{p*}\\ {\omega }_r& {\gamma }^p\ge {\gamma }^{p*}\end{array}\right.,(8)$

where ω denotes any one of m, s, σ_c, a, ψ, E, and v in the ith annulus; ω_p and ω_r correspond to the elastic and residual state values, respectively; and γ^p* is the critical plastic shear strain, beyond which the surrounding rock transitions into the residual state.

It is important to note that the frost heave mainly occurs due to the expansion of water when its temperature decreases below 0. Therefore, the porosity ratio depends on the frost heave coefficient for the water-rich surrounding rock. For the strain-softening surrounding rock, the macro cracks gradually increase once the crack initiation stress is reached. These cracks are much larger than the potential pores of intact rock, so the frost heave ratio also depends on the failure degree, which can be described by the volume plastic strain for strain-softening rock mass. For the plane strain problem, the volume plastic strain equals the plastic shear strain due to the zero plastic strain for the intermediate strain when the M-C plastic potential function is employed. In this study, the expansion coefficient is assumed to be linear and related to the plastic shear strain (Equation 9):

$\xi =b{\gamma }^p+{\xi }_0,(9)$

where ξ is the frost heave coefficient of the surrounding rock; ξ₀ is the frost heave ratio of the elastic rock; and b is a fitting coefficient.

3 Solutions of stress and displacement

Since the frost heave ratio depends on the failure degree, the stress and displacement induced by excavation should be solved first. Subsequently, a methodical approach involving a progressive frost heave process, layer by layer, is employed to derive theoretical solutions for the stress and displacement within the structure.

3.1 Stress and displacement induced by excavation

On the elastic–plastic boundary, the stress and strain under plane strain conditions can be easily obtained by use of the peak strength parameters and the initial in situ stress p₀. The explicit analytical expression of σ_k can only be obtained when a = 0.5 (Park and Kim, 2006):

${\sigma }_k=\frac{m_p{\sigma }_{cp}}{8}-\frac{1}{8}\sqrt{m_p^2{\sigma }_{cp}^2+16m_p{\sigma }_{cp}{p}_0+16s_p{\sigma }_{cp}^2}+p_0.(10)$

For the common condition of a > 0.5, only the implicit calculation equation of Equation 10 can be obtained. However, it can be easily calculated using the Newton method with an initial value of Equation 11.

$2\left(p_0-{\sigma }_k\right)={\sigma }_{cp}{\left(m_p{\sigma }_k/{\sigma }_{cp}+s_p\right)}^{a_p}.(11)$

Then, the stress and strain at the elastic–plastic boundary can be obtained as Equations 12, 13:

$\left[\begin{array}{l}{\sigma }_{r\left(k\right)}\\ {\sigma }_{\theta \left(k\right)}\end{array}\right]=\left[\begin{array}{c}{\sigma }_k\\ 2p_0-{\sigma }_k\end{array}\right],(12)$

$\left[\begin{array}{l}{\epsilon }_{r\left(k\right)}\\ {\epsilon }_{\theta \left(k\right)}\end{array}\right]=\frac{1}{2G_p}\left[\begin{array}{c}{\sigma }_k-p_0\\ p_0-{\sigma }_k\end{array}\right].(13)$

For the elastoplastic problem of plane strain circular openings, the radial stress gradually increases from the supporting pressure σ₀ to the initial stress p₀. By dividing the plastic region into n annuli with equal increments, the radial stress increment between the inner and outer regions can be obtained as Equation 14:

$\Delta \sigma =\frac{{\sigma }_0-{\sigma }_k}{k}.(14)$

For the strain-softening surrounding rock, the strength and deformation should obey the evolutional laws with the plastic shear strain. Meanwhile, the plastic shear strain also depends on the strength and deformation parameters. The difficulty is that they should be compatible at any point, which leads to a non-analytical solution. However, if the annulus number is large enough, the thickness of each annulus would be significantly small. The increment in shear plastic strain at the inner and outer radii for each annulus would also be significantly less. In this condition, the plastic shear strain at the outer radius can be employed to represent the annulus’s softening index. Then, the strength and deformation parameters of the ith annulus can be determined once the shear plastic strain at r = R_{i + 1} for the (i + 1)th annulus is obtained. Using the yield criterion of Equation 2, the circumferential stress at r = R_i+1 can be expressed as follows:

${\sigma }_{\theta \left(i\right)}={\sigma }_i+H\left({\sigma }_{i+1},{\gamma }_{i+1}^p\right).(15)$

Then, the circumferential stress increment can also be approximated as

$\Delta {\sigma }_{\theta \left(i\right)}={\sigma }_{\theta \left(i\right)}-{\sigma }_{\theta \left(i+1\right)}.(16)$

The equilibrium equation of Equation 1 can be expressed differentially by representing the variables with average values:

$\frac{{\sigma }_{r\left(i-1\right)}-{\sigma }_{\left(i\right)}}{{\rho }_{i-1}-{\rho }_i}=\frac{H\left({\sigma }_{r\left(i-1\right)},{\gamma }_i^p\right)+H\left({\sigma }_{r\left(i\right)},{\gamma }_i^p\right)}{{\rho }_{i-1}+{\rho }_i},(17)$

where ρ_i = R_i/R_p. Then, the dimensionless radius ρ_i can be explicitly expressed as

${\rho }_{i-1}=\frac{H\left({\sigma }_{r\left(i-1\right)},{\gamma }_i^p\right)+H\left({\sigma }_{r\left(i\right)},{\gamma }_i^p\right)+\Delta \sigma }{H\left({\sigma }_{r\left(i-1\right)},{\gamma }_i^p\right)+H\left({\sigma }_{r\left(i\right)},{\gamma }_i^p\right)-\Delta \sigma }{\rho }_i.(18)$

By using Hooke’s law, the compatibility equation can be reformulated as Equation 19:

$\frac{d{\epsilon }_{\theta \left(i\right)}^p}{d\rho }+\frac{{\epsilon }_{\theta \left(i\right)}^p-{\epsilon }_{r\left(i\right)}^p}{\rho }=-\frac{d{\epsilon }_{\theta \left(i\right)}^e}{d\rho }-\frac{H\left({\sigma }_{\left(i\right)},{\gamma }_{i-1}^p\right)}{2G_{i}r}.(19)$

When the average value of strain and H is employed, Equation 4 can also be reformulated in a differential way. Then, the circumferential plastic strain increment can be obtained as Equation 20:

$\Delta {\epsilon }_{\theta \left(i-1\right)}^p=\frac{-\frac{\Delta {\epsilon }_{\theta \left(i-1\right)}^e}{\Delta {\rho }_{i-1}}-\frac{H\left({\sigma }_{i-1},{\gamma }_i^p\right)+H\left({\sigma }_i,{\gamma }_i^p\right)}{2G_{\left(i\right)}\left({\rho }_{i-1}+{\rho }_i\right)}-\frac{2\left({\epsilon }_{\theta \left(i\right)}^p-{\epsilon }_{r\left(i\right)}^p\right)}{{\rho }_{i-1}+{\rho }_i}}{\frac{1}{\Delta {\rho }_{i-1}}+\frac{1+k_{\left(i\right)}}{{\rho }_{i-1}+{\rho }_i}},(20)$

$\frac{\Delta {\epsilon }_{\theta i}^p}{\Delta {\rho }_i}+\frac{2\left({\epsilon }_{\theta i-1}^p-{\epsilon }_{ri-1}^p\right)+\Delta {\epsilon }_{\theta i}^p-\Delta {\epsilon }_{ri}^p}{{\rho }_{i-1}+{\rho }_i}=-\frac{\Delta {\epsilon }_{\theta i}^e}{\Delta {\rho }_i}-\frac{H\left({\sigma }_{i-1},{\gamma }_i^p\right)+H\left({\sigma }_i,{\gamma }_i^p\right)}{2G_i\left({\rho }_{i-1}+{\rho }_i\right)},(21)$

where k_i = (1 + sinψ_i)/(1-sinψ_i).

Then, the radial plastic strain, shear plastic strain, volume plastic strain, and total strain at R_i-1 are updated as

$\left\{\begin{array}{l}{\epsilon }_{\theta i-1}^p={\epsilon }_{\theta i}^p+\Delta {\epsilon }_{\theta i-1}^p\\ {\epsilon }_{ri-1}^p={\epsilon }_{ri}^p-k_i\Delta {\epsilon }_{\theta i-1}^p\end{array}\right.,(22)$

$\left\{\begin{array}{l}{\epsilon }_{ri-1}={\epsilon }_{ri}+\Delta {\epsilon }_{ri-1}^e+\Delta {\epsilon }_{ri-1}^p\\ {\epsilon }_{\theta i-1}={\epsilon }_{\theta i}+\Delta {\epsilon }_{\theta i-1}^e+\Delta {\epsilon }_{\theta i-1}^p\end{array}\right.,(23)$

${\gamma }_{ri-1}^p={\gamma }_{ri}^p+\Delta {\epsilon }_{\theta i-1}^p-\Delta {\epsilon }_{ri-1}^p.(24)$

All the stresses and strains can be calculated recursively until the radial stress decreases to the supporting pressure σ₀. The radius and displacement of the plastic region can be obtained as follows:

$R_i=\frac{R_0{\rho }_i}{{\rho }_1},(25)$

$u_i=\frac{R_0{\rho }_i{\epsilon }_{\theta i}}{{\rho }_1}.(26)$

The real supporting pressure provided by the lining structure due to the displacement of the excavation surface can be easily obtained

${\sigma }_L=\frac{E_L\left({\overline{R}}_0^2-{\overline{R}}_L^2\right)\left(u_0-{\overline{u}}_0\right)}{{\overline{R}}_0^2\left(1+v_L\right)\left[1-2v_L+{\left({\overline{R}}_L/{\overline{R}}_0\right)}^2\right]},(27)$

where ${\overline{R}}_0=R_0-{\overline{u}}_0$ and ${\overline{R}}_L=R_0-{\overline{u}}_0-t$.

If the real supporting pressure σ_L equals the assumed supporting pressure σ₀, the calculation results would be the actual solution. Otherwise, the assumed supporting pressure can be easily calculated using the secant method or dichotomy method.

3.2 Stress and displacement induced by frost heaving

Without loss of generality, the frost front radius R_f is assumed to be larger than the plastic radius R_p. The plastic region uses the same partition method as strain-softening behavior, and the region R_p < r < R_f is divided into m annuli. The elastic stress and displacement of the outer annulus induced by frost heaving can be easily obtained as follows:

$\left\{\begin{array}{l}{\sigma }_{r\left(k+m\right)}={\sigma }_{k+m}^f{\left(R_f/r\right)}^2\\ {\sigma }_{\theta \left(k+m\right)}=-{\sigma }_{k+m}^f{\left(R_f/r\right)}^2\\ u_{\left(k+m\right)}=-{\sigma }_{k+m}^f{R}_f^2/\left(2r{G}_{k+m}\right)\end{array}\right..(28)$

For the annulus in the frost region, the closed-form analytical solution for the stress and displacement can be uniformly expressed as follows:

$\left\{\begin{array}{l}{\sigma }_{r\left(i\right)}^f=-A_i{R}_i^2/r^2+C_i\\ {\sigma }_{\theta \left(i\right)}^f=A_i{R}_i^2/r^2+C_i\\ u_{\left(i\right)}^f=r\left[A_i{R}_i^2/r^2+\left(1-2v_i\right)C_i\right]/\left(2G_i\right)+{\xi }_i\left[{\overline{R}}_i-r\left(1+v_i\right)\right]\end{array}\right.,(29)$

where $A_i=\frac{{\sigma }_i^f-{\sigma }_{i-1}^f}{{\rho }_i^{f2}-1}$, $C_i=\frac{{\sigma }_i^f{\rho }_i^{f2}-{\sigma }_{i-1}^f}{{\rho }_i^{f2}-1}$, and ${\rho }_i^f=R_i/R_{i-1}$.

Using the displacement continuum condition of the adjacent annulus, the radial stress induced by frozen heaving can be recursively expressed as Equations 30, 31:

${\sigma }_{k+m-1}^f=\frac{{\sigma }_{k+m}^f\left[\left({\rho }_{k+m}^{f2}-1\right)n_{k+m}+\left(1-2v_{k+m}\right){\rho }_{k+m}^{f2}+1\right]+2G_{k+m}{\xi }_{k+m}\left(1+v_{k+m}\right)\left({\rho }_{k+m}^{f2}-1\right)\left(a/rou+b-1\right)}{2\left(1-v_{k+m}\right)},(30)$

${\sigma }_{i-1}^f=\frac{{\sigma }_i^f\left[1+\left(1-2v_i\right){\rho }_i^{f2}\right]}{K_i\left({\rho }_i^{f2}-1\right)}+\frac{2G_i\left(1+v_i\right){\xi }_i\left(a/rou+b-1\right)}{K_i}-\frac{n_i}{n_{i+1}{K}_i}\left[A_{i+1}{\rho }_{i+1}^{f2}+\left(1-2v_{i+1}\right)C_{i+1}^2\right]-\frac{2G_i\left(1+v_{i+1}\right){\xi }_{i+1}\left(a-1+brou\right)}{K_i}.(31)$

For the sake of safety, the lining structure is only allowed to work in the elastic state. The stress and displacement can also be obtained as Equation 32:

$\left\{\begin{array}{l}{\sigma }_{rL}^f={\sigma }_0^f\frac{1-{\left(R_L/r\right)}^2}{1-{\left(R_L/R_0\right)}^2}\\ {\sigma }_{\theta L}^f={\sigma }_0^f\frac{1+{\left(R_L/r\right)}^2}{1-{\left(R_L/R_0\right)}^2}\\ u_L^f=\frac{{\sigma }_0^f{R}_0^2\left(1+v_L\right)\left[1-2v_L+{\left(R_L/R_0\right)}^2\right]}{E_L\left(R_0^2-R_L^2\right)}\end{array}\right..(32)$

3.3 Calculation process

To calculate the frost heave force, the failure degree of the surrounding rock should be determined first. The detailed calculation flow is listed in the following steps:

Step 1: Input data of material properties, geometry, and the section number k and then space the radial stresses into k sections using Equation 1.

Step 2: Calculate the dimensionless radius ρ_i using Equations 22, 23 of the ith annulus, and then determine the material parameters of the adjacent inner annulus by Equations 10, 11 using the radial stress at ρ = ρ_i.

Step 3: Judge whether the shear plastic strain satisfies the critical state using Equation 21. If not, i = i-1, and go to step 2, or go to step 4.

Step 4: Calculate the dimensionless radial stresses and displacement of plastic rock mass using Equations 15, 16, 24, 25 for M-C rock mass and Equations 17, 18, 28, 29 for H-B rock mass;

Step 5: Calculate the shear plastic strain using Equations 26, 30 for M-C and H-B rock mass, respectively, and then determine the material parameters of the adjacent inner annulus using Equations 10, 11;

Step 6: Repeat step 5 until the innermost annulus; calculate the residual radius using Equation 27; and output the results.

4 Examples

4.1 Verifications for strain-softening rock mass

To verify the reliability of the analytical solution, the normally strained softened rock body is considered the research object, where the coefficients and friction angles are plastic strain functions. The typical input data used by Alonso et al. (2003), Lee and Pietruszczak (2008), and Wang et al. (2010) are employed, and the input data are visualized in Figure 2 and provided in Table 1. The outer radius of the elastic region is 1000R₀. The GRC curves are shown in Figure 2. The relative error between the two analytic solutions is relatively small, at only 1.42%. Moreover, it is approximately consistent with the analytical results of Zhang et al. (2012). So, the solutions are in good agreement with both the numerical and analytical results.

Figure 2

Figure 2. Ground response curve of surrounding soils induced by excavation.

Table 1

Table 1. Values of the calculated parameters.

4.2 Analysis of freezing and expansion forces in excavated tunnels

Once the excavation of the cavern is completed, the stresses and deformations in the surrounding geotechnical medium produce secondary adjustments. However, artificial or natural freezing behavior occurs after the adjustment of surrounding rock stresses and deformations, and this section exclusively examines the freezing forces generated by freezing around the perimeter of the media. According to the mechanical properties of commonly utilized concrete, the elastic modulus E_L = 28 GPa, Poisson’s ratio v_L = 0.2, and the frost heaving coefficient of the medium are set as 0.055°C⁻¹. In the analysis of the freezing force of the lining structure due to frost heave, it is generally assumed that there is no displacement in the frozen zone due to frost heave behavior at the radius of the midpoint. This assumption serves as a supplementary equation for further calculating the freezing force exerted on the lining structure. However, the freezing behavior of geotechnical media develops gradually in stages, with each expansion of the frozen zone being constrained. It is more reasonable to analyze the assembly form caused by the expansion of the newly generated frozen zone, which means that it is considered that the displacement immobility point after the freezing and expansion of the newborn frozen zone is generated within the scope of the region. The effect of internal radial stress after media freezing under different ground freezing coefficients is illustrated in Figure 3. In this case, the maximum freezing expansion within the medium occurs within the frozen zone, and the freezing expansion force imposed on the lining structure is not at the point of maximum radial stress. When the coefficients of linear expansion for media freezing expansion were 0.0055, 0.0035, and 0.0015, the maximum radial stress induced by freezing expansion within the media was recorded as 18.55 MPa, 11.80 MPa, and 5.06 MPa, respectively. Conversely, the lining structure was exposed to stresses of only 10.62 MPa, 6.76 MPa, and 2.90 MPa, which are 57.26% of the maximum radial frost heave stress.

Figure 3

Figure 3. Distribution of the freezing force inside the medium for different freezing coefficients.

Figure 4 illustrates the influence of different freezing linear expansion coefficients of the medium on the freezing and expansion forces applied to the lining structure. It can be observed that the freezing force increases linearly with an increase in the freezing linear expansion coefficient. Specifically, each unit increase in the linear expansion coefficient results in a 1940 MPa surge in the freezing force sustained by the masonry structure. The freezing force on the masonry structure increases by 1935.83 MPa per unit coefficient of linear expansion. Therefore, the media freezing coefficient has a significant effect on the freezing force of the lining structure. Moreover, lining structures with different resistance characteristics generate distinct freezing and expansion forces, and the typical frost heave force characteristics of lining structures with different lining thicknesses are shown in Figure 5. As the thickness of the lining increases, the freezing force acting on the structure exhibits a nonlinear increase trend, and the growth rate of the freezing force decreases gradually with an increase in thickness. When the thickness of the lining structure increased from 10 cm to 50 cm, there was a corresponding increase in the freezing force exerted on it, from 2.90 MPa to 7.55 MPa, with a growth rate of 160.34%. Therefore, for weak strata, a thicker lining structure is often required to ensure the stability of the surrounding geotechnical body. If this type of formation structure develops frost heave behavior, the frost heave hazard tends to be greater than that in hard rock formations.

Figure 4

Figure 4. Frost heave of lining under different linear strain coefficients.

Figure 5

Figure 5. Frost heave of lining with different lining thicknesses.

4.3 Discussion

This paper conducts a theoretical analysis and research on the mechanical properties of multi-turn composite cylinder freezing walls. In future research, we will apply the theoretical structures obtained in this paper to practical engineering. We plan to use the temperature–stress field sequential coupling computational functionality of ANSYS finite element analysis software to establish a plane strain model of frozen soil interspersed with unfrozen soil and an axisymmetric vertical excavation model. This paper investigates the mechanical characteristics of single-, two-, and three-ring combined cylinder freezing walls, respectively. The influence of freezing deformation, horizontal displacement of excavation, and equivalent stress and strain under different unfrozen soil conditions is analyzed. The research results can provide a theoretical basis for the design and construction of freezing walls in super-large diameter deep foundation pit projects and support the numerical simulation of multi-turn combined cylinder freezing walls. The ultimate aim is to achieve a safe, reliable, economical, rational, and environmentally friendly design and construction technology for multi-ring composite pipe-freezing walls.

5 Conclusion

In this paper, a theoretical calculation model of the freezing expansion force generated by the freezing process around a circular hole under the axisymmetric plane strain model was developed. A corresponding variable model for the immobile point of the freezing circle displacement was proposed, and the progressive process of frost heave was realized through calculations. The conclusions drawn from this work are listed as follows:

(1) A variable model for the displacement of the frozen ring’s fixed point is proposed, taking into account the progressive formation process of freezing and expansion around the cavern. The theoretical formula for calculating the freezing force generated by the freezing process around a circular hole under the plane strain model is derived.

(2) The maximum frost heave force point of the medium around the circular hole usually occurs at a certain point inside the medium. The freezing force on the lining structure tends to increase linearly with the medium’s expansion coefficient.

(3) The thicker the permanent building structure, the greater the freezing force it experiences when the ground freezes. This is highly detrimental to the practical construction using the freezing method.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

LS: data curation and writing–original draft. MX: formal analysis and writing–review and editing. PW: validation and writing–review and editing. KW: project administration and writing–review and editing. XC: software and writing–review and editing. HL: resources and writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Natural Science Research Project of Colleges and Universities in the Anhui Province (grant number 2022AH050963).

Conflict of interest

Authors MX, PW, and KW were employed by Gu-Cheng Coal Mine of Shanxi Lu’an Mining Group.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alonso, E., Alejano, L. R., Varas, F., Fdez-Manin, G., and Carranza-Torres, C. (2003). Ground response curves for rock masses exhibiting strain-softening behaviour. Int. J.Numer. Anal. Methods.Geomech. 27 (13), 1153–1185. doi:10.1002/nag.315

CrossRef Full Text | Google Scholar

Feng, Q., Liu, W. W., and Jiang, B. S. (2017). Analytical solution for the stress and deformation of rock surrounding a cold-regional tunnel under unequal compression. Cold Reg. Sci. Tech. 139, 1–10. doi:10.1016/j.coldregions.2017.04.003

CrossRef Full Text | Google Scholar

Fu, Q., Yang, J., Gao, Y., Song, H., Liu, Y., et al. (2024). Combined blasting for protection of gob-side roadway with thick and hard roof. J. Rock Mech. Geotech. Eng. 16, 3165–3180. doi:10.1016/j.jrmge.2023.11.027

CrossRef Full Text | Google Scholar

Gu, Q. X., Zhang, Q., Ye, S. Z., Dai, W., and Li, T. (2024). Shear fracture behavior and damage constitutive model of rock joints considering the effect of pre-peak cyclic loading. Theor. Appl. Fract. Mech. 130, 104289. doi:10.1016/j.tafmec.2024.104289

CrossRef Full Text | Google Scholar

Ji, S., Lai, X., Cui, F., Liu, Y., Pan, R., and Karlovšek, J. (2024). The failure of edge-cracked hard roof in underground mining: an analytical study. Int. J. Rock. Mech. Min. Sci. 183, 105934. doi:10.1016/j.ijrmms.2024.105934

CrossRef Full Text | Google Scholar

Ji, Y. K., Zhou, G. Q., Zhao, X. D., Wang, J., Wang, T., Lai, Z., et al. (2017). On the frost heaving-induced pressure response and its dropping power-law behaviors of freezing soils under various restraints Cold Reg. Sci. Tech. 142, 25–33. doi:10.1016/j.coldregions.2017.07.005

CrossRef Full Text | Google Scholar

Jing, H. W., Wu, J. Y., Yin, Q., and Wang, K. (2020). Deformation and failure characteristics of anchorage structure of surrounding rock in deep roadway. Int. J. Rock. Mech. Min. Sci. 30, 593–604. doi:10.1016/j.ijmst.2020.06.003

CrossRef Full Text | Google Scholar

Lee, Y. K., and Pietruszczak, S. (2008). A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn.Underger. Sp. Tech. 23 (5), 588–599. doi:10.1016/j.tust.2007.11.002

CrossRef Full Text | Google Scholar

Park, K. H., and Kim, Y. J. (2006). Analytical solution for a circular opening in an elastic-brittle-plastic rock. Int. J. Rock. Mech. Min. Sci. 43 (4), 616–622. doi:10.1016/j.ijrmms.2005.11.004

CrossRef Full Text | Google Scholar

Shi, H., Chen, W., Zhang, H., and Lei, S. (2023b). A novel obtaining method and mesoscopic mechanism of pseudo-shear strength parameter evolution of sandstone. Environ. Earth Sci. 82, 60. doi:10.1007/s12665-023-10748-y

CrossRef Full Text | Google Scholar

Shi, H., Chen, W., Zhang, H., Song, L., Ming, L., Wang, M., et al. (2023a). Dynamic strength characteristics of fractured rock mass. Eng. Fract. Mech. 292, 109678. doi:10.1016/j.engfracmech.2023.109678

CrossRef Full Text | Google Scholar

Shi, H., Zhang, H., Chen, W., Song, L., and Li, M. (2023c). Pull-out debonding characteristics of rockbolt with prefabricated cracks in rock: a numerical study based on particle flow code. Comput. Part. Mech. 11, 29–53. doi:10.1007/s40571-023-00607-9

CrossRef Full Text | Google Scholar

Wang, S. L., Yin, X. T., Tang, H., and Ge, X. R. (2010). A new approach for analyzing circular tunnel in strain-softening rock masses. Int. J. Rock. Mech. Min. Sci. 47 (1), 170–178. doi:10.1016/j.ijrmms.2009.02.011

CrossRef Full Text | Google Scholar

Wu, J. Y., Jing, H. W., Gao, Y., Meng, Q. B., Yin, Q., and Du, Y. (2022). Effects of carbon nanotube dosage and aggregate size distribution on mechanical property and microstructure of cemented rockfill. Cem. Concr. Compos. 127, 104408. doi:10.1016/j.cemconcomp.2022.104408

CrossRef Full Text | Google Scholar

Wu, J. Y., Jing, H. W., Yin, Q., Yu, L., Meng, B., and Li, S. (2020). Strength prediction model considering material, ultrasonic and stress of cemented waste rock backfill for recycling gangue. J. Clean. Prod. 276, 123189. doi:10.1016/j.jclepro.2020.123189

CrossRef Full Text | Google Scholar

Wu, J. Y., Wong, H. S., Zhang, H., Yin, Q., Jing, H. W., and Ma, D. (2024). Improvement of cemented rockfill by premixing low-alkalinity activator and fly ash for recycling gangue and partially replacing cement. Cem. Concr. Compos. 145, 105345. doi:10.1016/j.cemconcomp.2023.105345

CrossRef Full Text | Google Scholar

Yuan, W., Lai, X. P., Pei, Z. H., Sifan, L., Haibin, L., Wei, W., et al. (2024). Experimental and numerical study on the anchoring mechanism of an anchor bolt considering its lateral restraint effect. Rock Mech. Rock Eng. 57, 9617–9633. doi:10.1007/s00603-024-03981-9

CrossRef Full Text | Google Scholar

Zhang, Q., Jiang, B. S., Wang, S. L., Ge, X. R., and Zhang, H. Q. (2012). Elasto-plastic analysis of a circular opening in strain-softening rock mass. Int. J. Rock Mech. Min. Sci. 50 (1), 38–46. doi:10.1016/j.ijrmms.2011.11.011

CrossRef Full Text | Google Scholar

Zhang, T., Gao, Y., Liu, S. Y., Yu, J., and Zhang, C. (2024). Experimental and numerical studies of breakage and fractal characteristics of silica sands in high-pressure triaxial tests. Adv. Powder Technol. 35 (7), 104548. doi:10.1016/j.apt.2024.104548

CrossRef Full Text | Google Scholar

Zhang, T., Zhang, C., Song, F. N., Zou, J. Q., Gao, Y., and Yang, W. H. (2023). Breakage behavior of silica sands during high-pressure triaxial loading using X-ray microtomography. Acta Geotech. 18 (10), 5195–5211. doi:10.1007/s11440-023-01866-9

CrossRef Full Text | Google Scholar

Zhang, W., Hu, L., Yao, Z.-B., Xiong, Y. R., Zhao, J., Ma, T., et al. (2024). In-situ and experimental investigations of the failure characteristics of surrounding rock through granites with biotite interlayers in a tunnel. Eng. Geol. 343, 107816. doi:10.1016/j.enggeo.2024.107816

CrossRef Full Text | Google Scholar