The numerical simulations are performed using the finite element (FE) software ABAQUS. A plastic damage constitution model is adopted with a criterion of Drucker-Prager yield [11], which was verified to be reasonable for study the problems of strain softening and progressive failure of soft rock mass [12] and recited as following:

If the stress estimation was beyond yield criterion, it appeared that the material damage initiated with plastic strain. According to strain equivalent law proposed by Lemaitre [13], the Cauchy stress σ _ij will be replaced by effective stress , while material stiffness is degraded progressively according to the specified damage evolution response. The relationships can be defined as following form:

(1)

where, is the renewed elastic stiffness, D is the value of damage variable, which generally ranges from 0 to 1. Considering the residual strength of the material, the maximum value of damage variable is set as 0.8 according to the uniaxial compression test of similar rock materials [12].

Thus, the initialization and evolution of damage variable are critical specifications of the model. The initiative damage is determined by the equivalent plastic strain ε₀pl, which is obtained by Equations (3). The damage variable will evolve with a linear form defined by equivalent plastic strain, which is given as:

(2)

where, ε _f^pl is equivalent plastic strain when the damage variable reaches the maximum. In the failure process of the model material, there are linear elastic, plastic hardening, and plastic softening phases. The relations between strain values at different phases are proposed as follows [14]:

(3)

where, the uniaxial compression strength σ_c can be obtained from Table 1 and the maximum elastic strain ε_e can be automatically passed by ABAQUS during calculation iteration.

A kind of 8-nodes linear brick element is introduced for rock mass surrounding tunnel, whose parameters are adopted according to the Table 1. An elastic structure element Truss3D adhere to ABAQUS is selected for simulating the anchor bolts with the elastic modulus of 300 GPa, Poisson’s ratio of 0.3 and geometric parameters presented in Fig. (1).

According to the size of tunnel, the study dimension simulated by numerical method is 80m×75m×20m. The model of finite element is constructed for studying the failure mechanism of tunnel collapse by mapping mesh as shown in Fig. (3), which includes 8160 elements and 9972 nodes. The lateral boundaries are only constrained along normal directions and the bottom one is restricted by three directions.

Fig. (3). FE model for numerical calculations.

For studying the support effects of anchor bolts of tunnel, some other numerical simulations are carried out by changing the geometric parameters. The schemes are presented in Table 2.

The process of tunnel collapse is reproduced by the numerical simulation as shown in Fig. (4), which is completed by 1000 sub steps in all after excavation. The red zones mean the failure zone with the complete damage of rock mass (D=0.8), the blue zones mean no damage of the rock mass (D=0), and the other zones indicate the partial damage of rock mass (0.8>D>0). At initial stage as shown in Fig. (4a), the failure zones distribute with some gaps at the positions of anchor bolts and the zones will develop as the calculation sub step increases as shown in Fig. (4b). Finally as shown in Fig. (4c), the terminal failure zones are almost continuous concentrated inner about 4m to tunnel surface which will bring on the tunnel collapse, while the damage zones are analogously distributed but greatly more than it.

(a)sub step=100 (b)sub step =500 (c)sub step =1000

The final failure zone of rock mass surrounding tunnel is compared between numerical simulation and field investigation, which is shown in Fig. (5). They are very similar to each other, especially for the part on the top of tunnel crown, which is always made up of the loosening load to tunnel structure. The comparison shows that the method of numerical simulation is reasonable and believable.

Fig. (4). Isograms distribution of damage variable in different sub step.

Fig. (5). Comparison of failure zone between field investigation and numerical simulation.

Generally speaking, tunnel excavation will break the initial stress balance of ground and cause stress redistribution [15]. According to the results of ground stress redistribution, the inner surrounding rock mass of tunnel can be divided into: stress loosening zone (Δσ_θ < 0, Δσ_r < 0), stress bearing zone (Δσ_θ > 0, Δσ_r < 0) and in-situ stress zone (Δσ_θ → 0, Δσ_r → 0), as presented in Fig. (6).

Fig. (6). Schematic diagram of stress redistribution of surrounding rock mass.

The tangential stress increment Δσ_θ and radial one Δσ_r of tunnel hance are presented in Fig. (7). Their distributions of variation with the depth greatly agree with the rules of Huang’s results as shown in Fig. (6) [15]. The stress loosening zone derived from ground is about 4 m according to the rule of stress increment Δσ_θ < 0.

Fig. (7). Variations of tangential stress increment along depth under different overload.

The damage zone, failure zone and loosening zone with the different length of bolts are compared as shown in Fig. (8). It is not difficult to find out that the 3 kinds of zones have the analogous rules: the more is the length, the less is the zone area. However, the zones variation with length of bolt become less obvious when the length is more than 5 m. Therefore, the length of more than 5 m is usually not expected, which also meets the requirements of Chinese design codes of road tunnel [10].

Fig. (8). Comparisons of zones between the cases of different length of bolts.

The tangential stress increments Δσ_θ of tunnel hance with the different length of bolts are presented in Fig. (9). Since the rock mass are strengthened by anchor bolts, the peak stress increment move inner tunnel, which means there will be less loosening zone with more length of bolts. Since the results of length of 5 m are almost same as the case of 8 m, the support effect of bolt is very little when the length is over 5m. Consequently, the length of bolts is not expected to be over 5 m in practice for not wasting money.

Fig. (9). Variations of tangential stress increment along depth in the case of different length of bolts.

The damage zone, failure zone and loosening zone with different spacing of bolts are presented in Table 3. The 3 kinds of zone also become less with less spacing of bolts. In fact, the results of spacing of 1m are almost the same as the one of 0.5 m. It indicates the support effect of bolt is very little when the spacing of bolts is less than 1m, which is unnecessary and uneconomical in practice. It agrees with the Zou and Wang’s results from model tests in laboratory, which suggested that the ratio between the length and spacing of bolts should be less than 5 [16].

Fig. (10). Variations of tangential stress increment along depth in the case of different anchor distance.

The tangential stress increments Δσ_θ of tunnel hance with the different spacing of bolts are presented in Fig. (10). With the decrease in spacing of bolts, the stress increments move inner tunnel because the strength of rock mass is enhanced by bolts. By comparison, the effects of bolts spacing are more notable than the effects of bolt length under other same conditions.