Analytical modeling of complex contact behavior between rock mass and lining structure

Dingli Zhng, Tong Xu, Hungheng Fng,*, Qin Fng, Liqing Co, Ming Wen

a Key Laboratory for Urban Underground Engineering of the Ministry of Education, Beijing Jiaotong University, Beijing,100044, China

b College of Civil and Transportation Engineering, Shenzhen University, Shenzhen, 518060, China

c China Academy of Transportation Sciences, Beijing,100029, China

Keywords:Loose contact Local non-contact Tunnel lining Complex variable method Analytical mechanical model

A B S T R A C T

1. Introduction

Unreasonable design, inferior materials, and groundwater erosion usually lead to unfavorable contact conditions between rock mass and lining structure (Gao et al., 2014; Meguid and Dang, 2009; Qin et al., 2020). The poor conditions such as cavities and insufficient grouting may result in a discontinuous displacement field and stress concentration in the rock-support system (Ding et al., 2019; Wang et al., 2014). Severe stress concentration can typically cause cracking or damage to the supporting structure, leaving a potential threat to tunnel safety.Therefore, the complex contact behavior between the rock mass and lining structure must be investigated when evaluating the safety of tunnel engineering.

Presently, the mechanical responses induced by poor contact conditions are usually investigated by model tests (Cao et al.,2019; Zhang et al., 2019; Min et al., 2020; Qin et al., 2020), numerical simulation (Jones and Hunt, 2011; Meguid and Kamel,2014; Xiao et al., 2018), and analytical methods (Lu et al., 2015,2019; Sun et al., 2021a, b; Yasuda et al., 2017, 2019; Zhang et al.,2020a). The model test can qualitatively reflect the mechanical behavior of the rock-lining system. However, it requires specialized test equipment and expensive test materials. In addition, it cannot easily measure the contact behavior at the interface. For complex engineering conditions, numerical simulation provides an effective way to calculate the mechanical response of the rock-lining system.However,when dealing with nonlinear contact problems,numerical simulation faces problems such as difficulty in achieving convergence and high computational cost(Hüeber et al.,2008;Fang et al.,2020,2021a,2022).By contrast,the analytical method has a low computational cost and high accuracy, making it more suitable for the preliminary assessment of tunneling projects (Zhang et al., 2020b; c, 2021).

The complex variable method proposed by Muskhelishvili(1963) has been extensively applied in the elastic analysis of the rock-lining interaction (Dong et al., 2019; Exadaktylos and Stavropoulou, 2002; Huo et al., 2006; Kargar et al., 2014, 2015,2020; Lu et al., 2016; Zhang et al., 2018; Fang et al., 2021b). In the initial studies, the concrete lining and surrounding rock are usually assumed to be completely adhered (Kargar et al., 2014), which implies that the displacement and stress at the interface are continuous without any jumps. However, in practice, the adhesion state between rock mass and concrete lining cannot be maintained when the tangential stress is greater than the frictional strength.In recent studies, the interaction between rock mass and lining structure has also been analyzed by assuming a full-sliding contact condition (Lu et al., 2015, 2019). In fact, neither of the above two assumptions is realistic, because the slip and adhesion states may exist simultaneously at the rock-lining interface.

In the case of cavities behind the lining (i.e. local non-contact condition), obtaining analytical solutions for the rock-lining interaction will be further challenging. To solve this problem,Yasuda et al. (2017, 2019) proposed a point matching method to obtain the elastic solution of a circular tunnel with voids behind the lining. In their method, the boundary conditions only need to be satisfied at a set of finite points, which considerably reduces the complexity of calculation.However,their method does not consider the re-closure of cavities behind the lining. On the other hand, in actual projects,engineers usually use backfill grouting to eliminate cavities.If the backfill grouting is insufficient,contact stress transfer between rock mass and lining structure still cannot be effectively achieved (i.e. loose contact condition), resulting in local stress concentration in the rock-lining system (Oh and Ziegler, 2014).Thus far,no analytical models are available to systematically study the insufficient grouting problem.

Although many analytical studies have been conducted to investigate the rock-lining interaction, three key issues remain to be addressed, namely the re-closure of cavities behind the lining, loose contact due to insufficient grouting, and tangential frictional slip at the interface. To this end, we classify the contact modes of the rock-lining system into the dense contact mode,local non-contact mode, and loose contact mode. Then, the corresponding analytical model is established for each contact mode using the complex variable method, where the linear complementarity function and Mohr-Coulomb friction law are utilized to characterize the normal and tangential contact behaviors,respectively. Furthermore, a new method is also proposed to obtain the accurate conformal mapping function of the concrete lining.Finally,the proposed method is used to study the complex contact behavior between rock mass and lining structure.

2. Statistical analysis of the contact state of operating railway tunnels

In the past few years,we have conducted nondestructive testing on more than 100 operating tunnels in China to investigate the contact mode behind the tunnel lining (Zhang et al., 2013). As a nondestructive testing method, we used the ground penetrating radar(GPR),and its basic detection principle is illustrated in Fig.1a.To avoid affecting the traffic in the operating tunnels, the detection vehicles must pass through quickly during the inspection process.For this reason, we only layout the survey line in the longitudinal direction.

According to the detection images shown in Fig. 1, we can classify the contact modes into three categories: dense contact,local non-contact,and loose contact.The first one,dense contact,is the most common mode,where the concrete lining is in full contact with the surrounding rock without any defects. The second one,local non-contact, is mainly caused by over-excavation without backfill grouting,whereas the last one,loose contact,is mainly due to over-excavation with insufficient backfill grouting.Hereafter,we refer to the latter two collectively as poor contact modes.

Fig. 2 shows the statistical results of field test data of contact states behind tunnel linings, where the vertical coordinate is the proportion of the poor contact to the total length of tunnels.We can observe that the height of the poor contact zone is mainly concentrated in the range of 0.1-0.6 m. That is, the probability of poor contact modes occurring on the tunnel vault is greater than that on the tunnel waist and arch foot. The percentage of loose contact is larger than that of local non-contact. This is mainly because most over-excavations can usually be detected and backfilled during the construction process.For more field test data and statistical analysis,please refer to our previous study(Zhang et al.,2013).

Fig.1. Nondestructive testing method: (a) Ground penetrating radar; and (b, c) Detection images (Zhang et al., 2013).

Fig. 2. Statistical results of two poor contact modes: (a) Local non-contact; and (b) Loose contact (Zhang et al., 2013).

3. Analytical expression of complex contact behavior in the rock-lining system

The poor contact state behind the lining can seriously affect the mechanical response of the rock-lining system, even posing a considerable risk to the tunnel safety. Therefore, developing an analytical method for rock-lining contact analysis is crucial. However, solving the following two problems in the analytical aspect is challenging: how to accurately characterize the plastic behavior of geotechnical materials and how to properly model the complex contact interaction between the rock mass and lining structure.Under poor geological conditions, the rock mass usually undergoes plastic deformation after tunnel excavation. Compared with the elastic stage, the plastic deformation stage is a highly nonlinear process. It is almost impossible to use an explicit analytical expression to characterize the stress and displacement distributions in the plastic zone of surrounding rock,except in some simple cases such as the circular tunnel under hydrostatic pressure. Moreover, a combination of contact nonlinearity and material nonlinearity further increases the complexity of the problem. Even in the numerical simulation, the contact behavior of plastic materials often leads to the non-convergence of the iterative calculation. Based on these reasons, this paper mainly addresses the second problem of developing an appropriate model for the complex contact interaction between rock mass and lining structure.

In Section 2, we propose three contact modes according to the field test results. In this section, the corresponding mechanical models will be established with the following assumptions:

(1) The rock-lining system is simplified as a two-dimensional(2D) plane strain model.

(2) The tunnel is located at great depth and subjected to nonhydrostatic stress.

(3) Both the concrete and rock materials are characterized by the linear elastic constitutive model.

3.1. Complex variable method

In the complex variable theory, the stress and displacement at any material point can be determined by two potential functions:?(z)and φ(z).In this way,the stress field and the displacement field induced by tunneling can be written as (Muskhelishvili,1963):

where σρ, σθ, and τρθare the Cauchy stress components;uρ anduθ are the displacement components;α is the dip angle;zdenotes the coordinate value,andz=x+yi,where i is the imaginary unit;?(z)and φ(z)are the two analytical functions;κ=3-4ν,where ν is the Poisson’s ratio;Gis the shear modulus;σ∞ρ， σ∞θand τ∞ρθare the farfield stresses,which can be defined as

where λ is the lateral pressure coefficient, and σ∞is the far-field stress component in the vertical direction.

Using conformal mapping, the arbitrarily shaped boundary can be mapped into a simple circle/ring, as shown in Fig. 3. The conformal transformations can be expressed as follows:

Fig.3. Conformal transformations of the tunnel excavation boundary and the concrete lining.

For the tunnel boundary:

For the concrete lining:

with

wheredk,ek,andfkare the complex coefficients;nis the truncated number of the conformal transformation;ris the inner radius of the mapped lining ring;ζ is the mapping coordinate system; the subscripts R and L represent the rock and lining; and ρ and θ are the radius and angle of the polar coordinate system,respectively.

Since the conformal mapping function above is analytical, the two potential functions (Eqs. (9) and (10)) can be analytically derived. In this case, the potential functions can be expanded into the following Laurent series:

For the rock mass:

For the concrete lining:

where the complex coefficientsak,bk,gk,hk,pk,andqkcan be solved by the given boundary and contact conditions;mis the truncated number of Laurent series; φR(ζ), ψR(ζ), φL(ζ) and ψL(ζ) are the analytical functions in the mapping plane;andDRandDLrepresent the rigid displacements of rock mass and concrete lining, respectively, which can be determined by fixing the displacement of a material point.

3.2. Solving strategy for determining the conformal mapping function

In this section, we only discuss the solving strategy of the conformal mapping function of concrete lining. The mapping function of the tunnel excavation boundary can be obtained in a similar way. The method proposed below can be regarded as an extension of the мелентьев method (Chen,1994).

First, dividing both sides of Eq. (10) by ζ yields

The above expression can be divided into real and imaginary parts as follows:

where the superscript“?”denotes the real part and“??”denotes the imaginary part of each complex coefficient,i.e.ek=.Then,Fourier transform is performed to Eq. (15) on the outer and inner boundaries, respectively:

For the outer boundary:

Through Eqs.(17)and(18),complex coefficients in the mapping function can be determined by

The integrals in the above equations can be obtained by using the fast Fourier transform technique. The mapping relationship between points in the physical plane and the mapping plane is unknown; hence,uandvcannot be given directly. Therefore, the following iterative process should be employed to determine the conformal mapping function:

(1) In the mapping plane,select 2sequally spaced points on the inner and outer boundaries, respectively.

(2) Repeat the first step in the physical plane.

(3) Use the following formula to estimate the mapping inner radius:

wherezouterandzinnerare the physical coordinate values of inner and outer boundaries, respectively.

(4) Assume that the points on the physical plane and those on the mapping plane have a one-to-one correspondence, and substitute them into Eqs. (19)-(23) to calculate all the unknown complex coefficients.

(5) Substitute the mapping coordinates of discrete points into Eq. (10) to update the physical coordinates.Since Eqs. (19)-(22) only use coordinateuto determine the unknown coefficients,the updated physical points may not be located on the lining boundary, as shown in Fig. 4.

(6) Calculate the distances Δzbetween the physical points and the lining boundary. If max(Δz) is less than the predefined tolerance ε, exit the loop; otherwise, correct the physical points using the method given in Fig.4,and then go back to Step 3.

Fig. 4. Correction method of mapping points.

It should be remarked that, as the requirement of fast Fourier transform,the discrete point number(2s)on each boundary needs to be greater than the mapping ordern.

3.3. Mechanical model for different contact modes

Without loss of generality, we simplify the contact analysis model of the rock-lining system as shown in Fig.5.For the sake of simplicity, we define the three independent geometric surfaces as tunnel surface Γ1, lining outer surface Γ2, and lining inner surface Γ3. According to the different contact modes,the contact interface can be further divided into three sub-interfaces: dense contact interface Γp, local non-contact interface Γn, and loose contact interface Γlc.

3.3.1. Analytical expression of normal contact behavior

Here, we first discuss the normal contact behavior of different contact modes. On the dense contact interface, the concrete lining completely adheres to the rock mass. Therefore, the normal displacement and stress are continuous on the interface Γp. The displacement continuity condition can be expressed as

and the stress continuity condition can be expressed as

On the local non-contact interface, the concrete lining is completely separated from the surrounding rock. Therefore, the displacement and stress on the interface Γndo not need to satisfy the continuity condition. In this case, the normal stress on the tunnel boundary and lining outer surface should be set to zero,and it takes the form:On the loose contact interface, the cavity behind the lining is filled with uncompacted grouting material. In general, the elastic modulus of these grouting materials is considerably smaller than that of the rock mass. Here, we consider the grouting materials equivalent to a series of Winkler springs,as shown in Fig.5.In the normal direction, the strain and stress of grouting material can be approximately expressed as

Fig. 5. Contact analysis model of the tunnel-lining system.

with the normal relative displacement Δuρ:

whereEgand ερ are the elastic modulus and normal strain of grouting material,respectively;andgap(z)is the shape function of the cavity.In this way,the governing equation on the loose contact interface can thus be written as

wherekρ is the normal stiffness coefficient of the Winkler spring.It can be given by

3.3.2. Analytical expression of tangential contact behavior

In the tangential direction of contact interfaces, the tangential stress is continuous for all contact modes, which takes the form:

For the dense contact mode, we assume that the rock-lining interface satisfies the following Mohr-Coulomb criteria:

where μ andcare the friction coefficient and cohesion,respectively.Then, the tangential contact behavior on the dense contact interface can be expressed as follows:

For the sticking state:

For the sliding state:

On the local non-contact interface, the tangential stress should be set to zero, which can be written as

On the loose contact interface, the shear stress of grouting material can be approximately expressed as

whereGgand γρθare the shear modulus and shear strain of grouting material, respectively. Then, the tangential contact behavior of the loose contact mode can be written as

wherekθ is the tangential stiffness coefficient; and Δuθ is the tangential relative displacement, which are given as

Furthermore,since the lining inner boundary is not subjected to any external load, the corresponding boundary condition can be written as

3.3.3. Analytical expression of re-contact behavior

It should be noted that the local non-contact mode would locally transform into the re-contact mode when the normal relative displacement is greater than the allowable penetration.In this case,the stress is continuous but there is a jump in the normal displacement. Then, the governing equations (Eqs. (27) and (28))need to be replaced by the following form:

where Γvdenotes the re-contact interface.

The re-contact process of the non-contact interface is a highly nonlinear problem. Therefore, it needs to be solved by a special algorithm, which will be discussed in the next section.

4. Solving process for determining the mechanical responses of the rock-lining system

The present problem is a mixed boundary value problem,where the displacement and stress boundary conditions exist simultaneously on the contact interface. In this case, conventional methods, such as the power series method, cannot be directly utilized to analyze the contact mechanical behavior between rock mass and lining structure.To this end,an efficient solving strategy will be introduced in this section.

The series expansions of the potential functions given in Eqs.(12) and (13) indicate that there are 2(6m+10) unknowns in the current problem.Therefore,at least 2(6m+10)linear equations are required to uniquely determine these unknowns. These equations can be obtained by satisfying the boundary and contact conditions at a set of finite points (Yasuda et al., 2017, 2019). For the circular tunnel, we can select 2m+3 discrete points equally spaced at the lining inner boundary and substitute them into the stress boundary conditions given in Eqs. (46) and (47), obtaining 2(2m+3) basic equations.Similarly,we select 2m+3 discrete points on the contact interface and substitute them into the governing equations of different contact modes given in Eqs.(25)-(45),obtaining 2(4m+6)basic equations. The last two equations can be obtained by fixing the vertical and horizontal displacements of a material point in the rock mass. Compared with circular tunnels, noncircular tunnels have a more complicated mechanical field;therefore,more discrete points are required to solve the unknowns in the potential function.Here, we suggest selecting 2(2m+3) discrete points on each boundary, which results in a larger number of equations than the unknowns. In this case, the QR decomposition or least square method should be used to solve the linear equations.

The basic linear equations of dense and loose contacts can be directly obtained by the above procedure. However, for the local non-contact mode, it involves the partial re-contact between the rock mass and lining structure, thus an additional process is required to determine the normal contact state.Here,we introduce a complementary function(Hintermuller et al.,2003;Hüeber et al.,2008) to detect whether the selected discrete points on the local non-contact interface are re-contacted. This complementary function takes the form:

wherezjis the physical coordinate of thejth discrete point,andEis the elastic modulus.Then,the normal contact state of thejth point can be determined by the following dimensionless criteria:

Finally, the solving strategy for determining the mechanical responses of the rock-lining system can be summarized as follows:

(1) Select 2(2m+3)discrete points on the lining inner boundary and substitute them into the stress boundary conditions(Eqs. (46) and (47)) to obtain the stress equations.

(2) Select 2(2m+3) discrete points on the contact interface and judge the contact mode of each point. Then, establish the basic system equations of normal contact behavior as follows:

(i) If the discrete pointzjis located at the dense contact interface, substitute it into Eqs. (25) and (26);

(ii) If the discrete pointzjis located at the loose contact interface, substitute it into Eqs. (32) and (33);

(iii) If the discrete pointzjis located at the local non-contact interface,assume that this point will not re-contact and substitute it into Eqs. (27) and (28).

(3) Assume that all points are in the sticking state,and substitute them into Eq. (38) to obtain the basic system equations of tangential contact behavior.

(4) Use the QR decomposition or least square method to solve system equations and update the mechanical field of the rock mass and concrete lining. The displacement and stress results should be processed by the Gaussian smoothing method to eliminate numerical errors.

(5) Use the criterion given in Eq. (51) to judge whether the points on the local non-contact interface are re-contacted.If the discrete pointzjis re-contacted, substitute it into Eqs.(48)and(49),else substitute it into Eqs.(27)and(28).Then,update the system equations.

(6) Use the criteria given in Eqs. (36) and (37) to judge the tangential contact state of each discrete point.If the discrete pointzjis in the sliding state,substitute it into Eq. (39), else substitute it into Eq.(38).Then,update the system equations.

(7) If||ru||<εuand||rσ||<εσ,the calculation ends,else go back to Step 4, whereruandrσ represent the residuals of the displacement and stress,and εuand εσ are the two predefined tolerances, respectively.

5. Verification of the proposed method

In this section,we considered a deeply buried noncircular tunnel to verify the effectiveness of our method. The calculation model is plotted in Fig.6,where a local non-contact zone is set in the range of 0.1π ≤β ≤0.3π, a loose contact zone is set in the range of 0.7π ≤β ≤0.9π, and the rest is assumed the dense contact mode. This tunnel is located at great depth and subjected to a horizontal pressure of 1 MPa and vertical pressure of 2 MPa. The elastic modulus and Poisson’s ratio of the rock mass are set to 1.8 GPa and 0.3,whereas those of concrete lining are set to 25 GPa and 0.2, respectively. On the loose and dense contact interfaces, the friction coefficient and cohesion are set to 0.8 and 1.2 MPa,respectively,whereas those of the local non-contact interface are both assumed zero. In addition,the elastic modulus of grouting material is set to 0.18 GPa.

5.1. Accuracy verification

To provide a more convincing presentation, we also used the finite element method (FEM) to solve the current problem. The results of contact stress and relative displacement obtained by the two methods are compared in Fig.7,where the number of Laurent series terms in our analytical method is 250. As shown in the figure, our results agree well with the FEM results. The maximum error of contact stress is located in the transition region of different contact modes,which is mainly attributed to the stress concentration in this region.The maximum error of relative displacement appears in the middle of the poor contact interfaces.Note that since the FEM adopts spatial approximation,the errors discussed here also come from the FEM.On the other hand,compared to the FEM,our analytical method does not require meshing or solving large equation systems.Hence,it has higher computational efficiency and is very suitable for the rapid safety assessment of tunneling projects.

5.2. Sensitivity study

To further investigate the accuracy of our method,the sensitivity of analytical solutions to the numbers of Laurent series terms and discrete points is studied here. The normal contact stress distributions with respect to the number of Laurent series terms,m,are given in Fig. 8a, where the number of discrete points on each boundary is 2(2m+3). It can be found that our analytical solution can converge rapidly to the FEM reference solution with the increase of the number of Laurent series terms.When the number of Laurent series terms is greater than 150,the results obtained by our method are basically consistent with the FEM solution except in the stress concentration region.

Fig. 6. Schematic diagram of the verification example: (a) Calculation model; and (b) Geometric detail.

Fig. 7. Comparison of results from the two methods: (a) Contact stress; and (b) Relative displacement.

Fig. 8. Sensitivity of analytical solutions to the numbers of (a) Laurent series and (b) discrete points.

The normal contact stress distributions with respect to the number of discrete points are given in Fig. 8b. In this example, we consider three cases of discrete points,namely 1.4N,1.7N,andN,on each boundary (whereN= 2m+3), while the number of Laurent series terms is fixed atm=250.As Fig.8b shows,if the number of discrete points was too small, the obtained analytical solution would have a large oscillation error at the tunnel invert.When the number of discrete points is greater than 1.7N,this oscillation error disappears and the result accuracy is satisfactory. Therefore, the suggestion of selecting 2Ndiscrete points on each boundary given in Section 4 is reasonable.

6. Mechanical response of concrete lining under different contact modes

Poor contact conditions have a considerable effect on the safety of the rock-lining system. To analyze this effect, three parametric investigations are performed in this section. The basic material parameters are listed in Table 1.In addition,the calculation models for the three investigations are given in Fig.9,and the geometry of the tunnel is consistent with that in Fig. 6b.

Table 1 Basic material parameters for the investigations under three contact modes.

6.1. Local non-contact

The first example is designed to investigate the effect of local non-contact on the interaction between rock mass and lining structure. For this purpose, a local non-contact zone is predefined at the tunnel vault in the range of 0.3π ≤β ≤0.7π.Moreover,four typical heights of the non-contact zone are considered:1 cm,2 cm,3 cm and 20 cm.

The evolution of the normal contact stress and relative displacement along the rock-lining interface with increasing noncontact zone height is illustrated in Fig.10. As the figure indicates,the local non-contact area only affects the stress distribution in a small range nearby, which is mainly concentrated in the range of 0 ≤β ≤π in this example.When the height Δhincreases from 1 cm to 3 cm, the normal contact stress in the center of the local noncontact zone decreases gradually, while the normal relative displacement increases rapidly. Note that the non-contact zone completely re-contacts when the height Δhis less than 1 cm.In this case,the stress concentration effect in the transition area of contact modes (β = 0.3π and 0.7π) is minimal. Because this special mode cannot be detected by the GPR technology,it may pose a potential threat to tunnel safety.When the height Δhis greater than 3 cm,the non-contact zone is almost no longer closed again. Comparison of the cases of Δh=3 cm and 20 cm indicates almost similar contact stress distributions as well as an extremely small difference in the normal relative displacements. Therefore, if no re-contact occurs,the effect of the non-contact zone height on the rock-lining interaction can be negligible.

Fig.10. Mechanical responses of the rock-lining interface under the effect of local non-contact: (a) Normal contact stress; and (b) Normal relative displacement.

6.2. Loose contact

The second example analyzes the effect of loose contact on the interaction between rock mass and lining structure. Here, we consider two influencing factors: the height of the loose contact zone Δhand the elastic modulus of grouting materialEg.The loose contact zone is also predefined at the tunnel vault in the range of 0.3π ≤β ≤0.7π.

Fig.11a plots the curves of the normal contact stress distribution for Δh= 20 cm, 40 cm, 60 cm, and 80 cm. The results of dense contact are also shown in Fig.11 as reference.Compared with local non-contact(if no re-contact occurs),the effect of loose contact on the normal contact stress is more significant, especially in the transition zone of contact modes (β = 0.3π and 0.7π). Fig. 11b shows the curves of the normal contact stress distribution forEg=0.05 GPa,0.1 GPa,0.5 GPa,and 1 GPa.Clearly,the contact stress curves show no obvious change in shape but a decrease in the stress concentration effect. When the elastic modulus of the grouting material is greater than one-quarter of that of the surrounding rock(Eg≥0.5 GPa andER= 2 GPa in this example), the stress redistribution caused by loose contact is negligible.

6.3. Dense contact with tangential frictional slip

The last example aims to analyze the tangential frictional contact behavior in the rock-lining system. Here, we test four cases with friction coeffciients set as 0,0.4,0.8,and 1.2.Note that no poor contact mode is considered in this example.

Fig.12a illustrates the circumferential stress distribution along the lining inner boundary under different friction coefficients. As can be seen, the maximum circumferential stress appears at β = 1.18π and 1.82π, around the tunnel arch foot. This maximum value is mainly due to the stress concentration caused by the lining geometry.The internal force of the lining near the tunnel arch foot(β = 1.18π and 1.82π) is barely affected by the tangential friction contact behavior. As the friction coeffciient decreases, the circumferential stress at the tunnel vault and invert increases gradually but the stress distribution becomes more uniform. Fig.12b shows the tangential relative displacements along the rock-lining interface under different friction coeffciients. The results show that,with increasing friction coefficient, the sliding area and the maximum sliding displacement gradually decrease. When the friction coefficient is greater than 0.8, the upper part of the lining(0 ≤β ≤π)is in the sticking state,whereas the tunnel invert is still in the sliding state.

Fig. 9. Calculation models for parametric analysis: (a) Local non-contact; (b) Loose contact; and (c) Dense contact with tangential frictional slip.

7. Conclusions

In this study, several analytical mechanical models were developed to model the complex contact behavior between rock mass and lining structure. The complex variable method was used to derive the mechanical field of the rock-lining system and the point matching method was introduced to solve the mixed boundary value problem. Furthermore, several examples are performed to demonstrate the accuracy and practicability of our method. The following conclusions can be drawn from the results of the study:

(1) Based on the field test data, the contact mode can be classified into local non-contact, loose contact, and dense contact. According to the mechanical characteristics of these contact modes, the governing equations of normal and tangential contact behaviors were established:the local noncontact state was characterized by free stress conditions;the loose contact state was modeled with a simplified Winkler model; and the dense contact state was characterized by stress and displacement continuity conditions.

(2) At a small local non-contact zone height,the increase in the surrounding rock deformation caused by the advance of the tunnel face may lead to the closure of the non-contact zone.This special mode cannot be detected by the GPR technology,which may have a potential impact on tunnel safety. At a large height of the local non-contact zone, the surrounding rock and concrete lining are no longer closed again. In this case, the height does not have a significant effect on the contact stress at the rock-lining interface.

Fig.11. Mechanical responses of the rock-lining interface under the effect of loose contact: (a) Height of the loose contact zone; and (b) Elastic modulus of grouting material.

(3) Different from local non-contact, the height of the loose contact zone has a significant effect on the contact stress.With the increase of the height, the stress concentration effect caused by the loose contact mode gradually increases.On the other hand, by increasing the elastic modulus of the grouting material, the stress concentration effect can be reduced effectively.

(4) The tangential frictional contact behavior on the rock-lining interface has a certain effect on the internal force of the lining. By reducing the friction coefficient, a more uniform internal force distribution of the lining can be obtained.

Fig.12. Mechanical responses of the rock-lining system under the effect of dense contact: (a) Circumferential stress; and (b) Tangential relative displacement.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 51738002 and 52108376) and Fundamental Research Funds for the Central Universities(Grant No.2021CZ111).

List of symbols

ak,bk,gk,hk,pk, qkComplex coefficients of analytical functions

dk,ek,fkComplex coefficients of conformal transformations

DR,DLRigid displacements

E,G,ν Elastic modulus, shear modulus, and Poisson’s ratio

Eg,GgElastic and shear modulus of grouting material

gap(z) Shape function of the cavity

kρ,kθ Normal and tangential stiffness coefficients

mTruncated numbers of the Laurent series

nMapping order of the conformal transformation

rInner radius of the mapped lining

uρ,uθDisplacement components

z,ζ Physical and mapping coordinate systems

zjPhysical coordinate values of the discrete points

zouter,zinnerPhysical coordinate of inner and outer boundaries

α Dip angle

Δuρ,ΔuθNormal and tangential relative displacements

ερ,γρθNormal and shear strains

φR(ζ),ψR(ζ),φL(ζ),ψL(ζ) Analytical functions in the mapping plane

κ Material parameter, and κ = 3 - 4ν

λ Lateral pressure coefficient

μ,cFriction coefficient and cohesion

?(z), φ(z)Analytical functions in the physical plane

ρ,θ Radius and angle of the polar coordinate system

σ∞ρ， σ∞θ， τ∞ρθFar-field stresses

σ∞Vertical far-field stress

σр,σθ,τрθCauchy stress components