Theoretical Analysis of Rock Blasting Damage in Construction of Tunnels Closely Under-Passing Sewage Box Culverts

Abstract

With the large-scale construction of urban traffic tunnels in China, it has become common to underpass existing buildings and structures such as sewage box culverts and pipelines using the drilling-blasting method. How to analyze accurately the blasting damage of surrounding rock and reasonably determine the safe distance between tunnel and box culvert or pipelines is an urgent issue to be solved. In this paper, the Cowper-Symonds plastic kinetic hardening model was improved using both rock initial damage degree and damage modification coefficient considering rock residual strength. The proposed model was implemented into LS-DYNA. The proposed damage model was used to evaluate the blasting construction of rock tunnels closely under-passing sewage box culverts. The results of numerical simulation using the proposed damage model shows that the blasting damage range of rock with a damage degree of more than 0.5 very significantly reduces from 1.0 m to 0.3 m as the spacing between the box culvert and the tunnel increases from 1.0 m to 4.0 m, and the evolution process of rock blasting damage can be well-presented. Moreover, the safe distance between tunnel and box culvert in blasting construction can be reasonably determined to be no less than 4.0 m. The findings in this paper could be significant for guiding the blasting construction of rock tunnels closely under-passing sewage box culverts.

1. Introduction

Tunnels often are constructed using the drilling-blasting method in many eastern Chinese cities. The blasting construction of rock tunnels must inevitably underpass existing buildings and structures such as sewage box culverts and pipelines. When the tunnel constructed using the drilling-blasting method underpasses through the structure of box culverts, the impact caused by rock blasting will damage the box culvert structure adjacent to the tunnel [1]. Blasting damage in the construction of tunnels is usually regulated in infrastructure contracts because it can influence the lifecycle cost and quality for tunneling projects [2]. The most critical issue is how to reasonably determine the safe distance between tunnel and existing box culvert in tunnel blasting construction. The control of explosion-induced crushed and cracked zones is of the utmost importance in the rock explosion design [3].

In numerical simulations of tunnel excavation, rock damage-plasticity models are often used to study the influence of both effects on the evolution of deformation and stress distribution in rock masses [4,5]. The research on the constitutive models of rock blasting damage is a developing process. The blasting damage model of rock (the GK model) [6] was used to simulate the process of rock blasting. The GK model was improved upon using the definition of crack density distribution with the effect of loading time and damage variables with the fracture probability [7,8]. Using the damage mechanics theory of rock can better simulate the formation and development process from microscopic cracks to macroscopic cracks of rock with native defects [9,10]. The rock tension-compression damage model [11] and tension-compression-shear damage model [12] were developed. The damage process of the surrounding rock under multiple full-face blasts is simulated using rock blast-induced damage model integrating the tensile damage model with the Drucker-Prager yield condition [13]. The rock blasting damage constitutive model coupling the tension-compression weighted damage variable and the classical PLASTIC_ KINEMATIC model was used to investigate the rock damage evolution law in slotting-blasting construction [14].

With the development of computer technology, numerical simulation technology has become one of effective methods to investigate the damage evolution of rock in the process of blasting excavation construction [15,16,17,18,19,20]. Surrounding rock damage is controlled by accurately describing rock characteristics such as Young’s modulus and Poisson’s ratio in tunnel blasting construction [17], and by optimizing charge structures and the number of blast-hole rows [18]. The damage degree of the blasting crushing zone was determined using rock fracture theory [19]. The Johnson-Holmquist-Rock (JHR) damage model is proposed to simulate the failure process of rock under blasting load [20]. Elastoplastic cellular automaton was employed to simulate blasting-induced elastoplastic dynamic responses [21].

Rock explosion is a complex problem because neither the rock characteristics nor the explosion waves could be accurately estimated. A first-order reliability analysis can predict rock blasting damage [22].

Through the use of ANSYS/LS-DYNA to establish a three-dimensional finite element model for slope and tunnel, the minimum safe distance from tunnel face to slope surface was determined according to the first strength theory and a linear statistical relationship between peak particle velocity and rock effective tensile stress [23]. Using the 3D numerical model established by LS-DYNA, the safety vibration velocity threshold of each vulnerable area was calculated based on the maximum tensile stress criterion during the whole process of tunnel blasting excavation [24]. The displacement of surrounding rock in tunnel construction can be predicted or estimated using the artificial neural network [25,26], numerical analysis, fuzzy logic, and statistical analysis models [26], and the gene expression programming [27]. Moreover, the plastic zone surrounding rock in tunnel excavation can be determined using numerical analysis, the fuzzy inference system, and the multivariate regression model [28].

At present, the minimum safe distance between sewage box culverts in operation and the tunnel constructed using the drilling-blasting method is usually determined by the peak blasting vibration velocity. However, in actual blasting construction, it is often difficult to monitor the blasting vibration velocity of the operating sewage box culvert structure embedded underground. Therefore, it is a very effective and feasible method to use numerical simulation calculation to determine the minimum safe distance based on the threshold value of the rock blasting damage degree of the surrounding rock. In addition, current blasting damage models rarely consider both the initial damage and the residual strength of rock.

In this paper, the Cowper-Symonds kinetic hardening model was improved using the rock initial damage degree and damage modification coefficient considering rock residual strength. Using the proposed damage model, blasting damage characteristics are investigated when the initial damage degree of the surrounding rock is different, and the residual strength of the damaged rock can be considered simultaneously. Then, dangerous parts of the box culvert structure are theoretically determined. Finally, the safe distance between tunnel and the existing box culvert in blasting construction can be determined using the blasting damage degree of the rock.

The novelty of this paper is that the Cowper-Symonds kinetic hardening damage model was improved upon by considering the rock initial damage degree and its residual strength after blasting. In addition, the safe distance between tunnel and box culvert in tunnel blasting construction can be determined reasonably using the rock blasting damage degree of the rock.

2. Methodology

2.1. Calculation of Total Tensional Strain

According to the continuous damage theory, Yang (1996) established a blasting damage model of rock [7]. In the Yang model, the equivalent tensile strain θ of the microelement is defined as the sum of the tensile strains (${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$) along the three principal directions, described by the logarithmic strain. When the equivalent tensile strain of microbody exceeds the critical value, the tensile failure occurs.

The logarithmic principal strain is denoted by Formula (1) [7]:

$$e_i=\ln {\lambda }_i, i = 1, 2, 3$$

(1)

where, ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are engineering strains along the first principal strain direction, the second principal strain direction, and the third principal strain direction, respectively.

The total tensional strain θ is denoted by Formula (2):

$$\theta ={\displaystyle \sum e_i^t, i = 1, 2, 3}$$

(2)

where $e_i^t=0$ if the principal strain is compressional; $e_i^t=e_i$ if the principal strain is tensile.

2.2. Strain Rate-Dependence and Cumulative Damage

The Yang’s model assumed that flaws and cracks exist within rock mass [7]. Under tension, these flaws and cracks may grow, and new flaws and cracks may be generated, leading to the reduction of the stiffness and strength of rock. It was assumed that the increase of crack density is controlled by the tensile strain of rock. Formula (3) was proposed to calculate the increase of the crack density [7,9,29]:

$$\mathrm{d}{C}_{\mathrm{d}}/\mathrm{d}t=\alpha {(\theta -{\theta }_{\mathrm{c}})}^{\beta }$$

(3)

where $C_{\mathrm{d}}$ is crack density; t is time; α and β are material constants; θ is equivalent tensile strain; ${\theta }_{\mathrm{c}}$ is critical tensile strain, and can be determined by the uniaxial dynamic tensile test, ${\theta }_{\mathrm{c}}=\ln (1+{\epsilon }_t)$; ${\epsilon }_t$ is ultimate tensile strain. When $\theta <{\theta }_{\mathrm{c}}$, crack density does not increase.

Damage variable D can be denoted by Formula (4):

$$D=1-{\mathrm{e}}^{-C_{\mathrm{d}}^2}$$

(4)

The relationship between the relative damage variable and the acoustic wave velocity for rock mass was denoted by Formula (5) when $0\le D\le 1$ and $c\le c_0$ [30]:

$$D=1-{(c/c_0)}^2$$

(5)

where c is the acoustic wave velocity of damaged rock after being loaded; $c_0$ is the initial acoustic wave velocity of rock.

When rock is the undamaged material before being loaded, D can be denoted by Formula (6):

$$D=1-{(c/\widehat{c})}^2$$

(6)

where $\widehat{c}$ is the acoustic wave velocity of undamaged rock. In this paper, rock mass with an integrity index of more than 0.8 is assumed as the undamaged material.

The initial damage value of rock $D_0$ can be calculated by Formula (7):

$$D_0=1-{(c_0/\widehat{c})}^2$$

(7)

The initial integrity index of rock mass ${\eta }_0$ is denoted by Formula (8):

$${\eta }_0={(c_0/\widehat{c})}^2$$

(8)

Then, the approximate relationship between initial damage value and initial integrity coefficient of rock mass is expressed by Formula (9):

$$D_0=1-{\eta }_0$$

(9)

Considering the effect of initial crack density, $D_0$ can be denoted by Formula (10) [9,29].

$$D_0=1-{\mathrm{e}}^{-C_{\mathrm{d}0}^2}$$

(10)

where $C_{\mathrm{d}0}$ is the initial crack density of the rock mass.

Considering the initial crack density of the rock, Formula (3) can be expressed by Formula (11) [9,29]:

$$C_{\mathrm{d}}={\displaystyle {\int }_{{\theta }_{\mathrm{c}}}^{\theta }\frac{\alpha }{\dot{\theta }}}{(\theta -{\theta }_{\mathrm{c}})}^{\beta }\mathrm{d}\theta +C_{\mathrm{d}0}$$

(11)

where $\dot{\theta }$ is the equivalent tensile strain rate; $C_{\mathrm{d}0}=\sqrt{-\ln {\eta }_0}$.

2.3. Equivalent Bulk Modulus and Shear Modulus of Damaged Rock

It is assumed that the Poisson’s ratio of damaged rock is the same as that of undamaged rock. In the blasting damaged process of rock, Young’s modulus and the shear modulus can be obtained by Formulae (12) and (13), respectively [7,31].

$$E_{\mathrm{d}}=E(1-D)/\left(3\left(1-2\nu \right)\right)$$

(12)

$$G_{\mathrm{d}}=G(1-D)$$

(13)

where $E_{\mathrm{d}}$ is Young’s modulus of damaged rock; $G_{\mathrm{d}}$ is shear modulus of damaged rock; E is Young’s modulus of undamaged rock; $\nu$ is Poisson’s ratio of rock; G is shear modulus of undamaged rock.

Young’s modulus and initial shear modulus with initial damage can be obtained by Formulae (14) and (15), respectively [7,31]:

$$E_0=\widehat{E}(1-D_0)/\left(3\left(1-2\nu \right)\right)$$

(14)

$$G_0=\widehat{G}(1-D_0)$$

(15)

where $E_0$ is initial Young’s modulus of rock with initial damage; $G_0$ is initial shear modulus of rock with initial damage; $\widehat{E}$ is equivalent Young’s modulus of undamaged rock; $\widehat{G}$ is equivalent shear modulus of undamaged rock.

When Formula (14) is substituted into Formula (12), Young’s modulus of damaged rock can be obtained by Formula (16) using initial Young’s modulus, initial damage value, and damage variable of rock [29]:

$$E_{\mathrm{d}}=E_0(1-D)/\left(3\left(1-2\nu \right)(1-D_0)\right)$$

(16)

When Formula (15) is substituted into Formula (13), the shear modulus of damaged rock can be obtained by Formula (17) using initial shear modulus, initial damage value, and damage variable of rock [29]:

$$G_{\mathrm{d}}=G_0(1-D)/(1-D_0)$$

(17)

2.4. Modification Coefficient of Damage Variable

After the rock is completely broken, there is still some residual strength [29,30]. The existing blasting damage models of rock rarely consider the characteristics that compressive stress and shear stress are still transmitting after the rock has been broken completely (i.e., the damage degree of rock is close to 1). The statistical damage-softening constitutive model considering the post-peak residual strength of rock was established based on the Hoek-Brown yield criterion [32]. The damage variable modification coefficient was introduced into the damage model proposed by Yumlu (1995) [29]. In this paper, the modification coefficient of damage variable $\delta$ is denoted by Formula (18) [32]:

$$\delta =1-\sqrt{\frac{{\sigma }_{\mathrm{r}}}{{\sigma }_{\mathrm{lc}}}}$$

(18)

where $\delta$ is the modification coefficient of damage variable; ${\sigma }_{\mathrm{r}}$ is the post-peak residual strength of rock; ${\sigma }_{\mathrm{lc}}$ is peak strength of rock.

Through Formula (18), Formulas (12) and (13) can be transformed to Formulae (19) and (20), respectively:

$$E_{\mathrm{d}}=E_0\left(1-\delta D\right)/\left(3\left(1-2\nu \right)\left(1-\delta D_0\right)\right)$$

(19)

$$G_{\mathrm{d}}=G_0(1-\delta D)/(1-\delta D_0)$$

(20)

2.5. Constitutive Model Considering Initial Damage and Residual Strength of Rock

When total stress is decomposed into average stress or hydrostatic pressure and deviatoric stress, the elastic damage constitutive relation of rock can be expressed by Formulae (21) and (22):

$${\sigma }_{\mathrm{v}}=E_{\mathrm{d}}{\epsilon }_{\mathrm{v}}/\left(3\left(1-2\nu \right)\right), {\sigma }_{\mathrm{m}}=E_{\mathrm{d}}{\epsilon }_{\mathrm{v}}/\left(9\left(1-2\nu \right)\right)$$

(21)

$$S_{ij}=2G_{\mathrm{d}}{e}_{jj}$$

(22)

where ${\sigma }_{\mathrm{v}}$ is bulk stress, ${\sigma }_{\mathrm{v}}={\sigma }_1+{\sigma }_2+{\sigma }_3$; ${\sigma }_{\mathrm{m}}$ is average stress or hydrostatic pressure, ${\sigma }_{\mathrm{m}}=\frac{1}{3}({\sigma }_x+{\sigma }_y+{\sigma }_z)$; ${\epsilon }_{\mathrm{v}}$ is bulk strain, ${\epsilon }_{\mathrm{v}}={\epsilon }_1+{\epsilon }_2+{\epsilon }_3$; ${\epsilon }_{\mathrm{m}}$ is average strain, ${\epsilon }_{\mathrm{m}}=\frac{1}{3}({\epsilon }_x+{\epsilon }_y+{\epsilon }_z)$; $S_{ij}$ is the component of deviatoric stress tensor; and $e_{jj}$ is the component of deviatoric strain tensor.

When Formula (19) is substituted into Formula (21), the model considering initial damage and residual strength of rock can be expressed by Formula (23). When Formula (20) is substituted into Formula (22), the constitutive model considering initial damage and residual strength of rock can be expressed by Formula (24):

$${\sigma }_{\mathrm{m}}=E_0{\epsilon }_{\mathrm{m}}(1-\delta D)/\left(3\left(1-2\nu \right)(1-\delta D_0)\right)$$

(23)

$$S_{ij}=2G_0{e}_{jj}(1-\delta D)/(1-\delta D_0)$$

(24)

where $D_0\le D\le 1$.

When denoted in increments, Formulae (23) and (24) can be transformed to Formulae (25) and (26), respectively:

$${\dot{\epsilon }}_{\mathrm{m}}=3\left(1-2\nu \right)(1-\delta D_0)\left({\dot{\sigma }}_{\mathrm{m}}(1-\delta D)+\dot{D}\delta {\sigma }_{\mathrm{m}}\right)/\left(E_0{(1-\delta D)}^2\right) (\theta >{\theta }_{\mathrm{c}})$$

(25)

$${\dot{e}}_{jj}=\frac{{\dot{S}}_{ij}(1-\delta D)+\dot{D}\delta S_{ij}}{2G_0{(1-\delta D)}^2/(1-\delta D_0)} (\theta >{\theta }_{\mathrm{c}})$$

(26)

where ${\dot{\epsilon }}_{\mathrm{v}}$ is average strain rate; ${\dot{\sigma }}_{\mathrm{m}}$ is average stress rate; ${\dot{e}}_{jj}$ is the component of deviatoric strain rate tensor; and ${\dot{S}}_{ij}$ is deviatoric stress rate.

2.6. Cowper-Symonds Plastic Kinetic Hardening Model

The stress-strain curves of rock under dynamic load show obvious strain rate-related properties. The limit yield strength of rock obviously increases with the increase of strain rate. In this paper, the Cowper-Symonds plastic kinetic hardening model was used to denote the stress-strain relation of rock.

The yield function $\varphi$ is denoted by Formula (27):

$$\varphi ={\sigma }_i-{\sigma }_{\mathrm{y}}$$

(27)

where ${\sigma }_i$ is stress strength, ${\sigma }_i=\sqrt{\frac{3}{2}(S_{ij}-{\alpha }_{ij})(S_{ij}-{\alpha }_{ij})}$; ${\sigma }_{\mathrm{y}}$ is yield stress; $S_{ij}$ is the component of deviatoric stress tensor, $S_{ij}={\sigma }_{ij}-\frac{1}{3}{\sigma }_{\mathrm{kk}}$; ${\alpha }_{ij}$ is movement tensor at the center of yield surface, $\Delta {\alpha }_{ij}=\frac{2}{3}(1-\beta )E_{\mathrm{p}}{\dot{\epsilon }}_{ij}^{\mathrm{p}}\Delta t$; ${\sigma }_{\mathrm{kk}}$ is the first invariant of stress tensor; $\Delta {\alpha }_{ij}$ is the movement tensor increment at the center of yield surface; β is hardening parameter of rock, and β = 1 when rock belongs to isotropic hardening; $E_{\mathrm{p}}$ is plastic hardening modulus, $E_{\mathrm{p}}=\frac{E_0{E}_{\tan }}{E_0-E_{\tan }}$; $E_0$ is Young’s modulus; $E_{\tan }$ is tangent modulus; ${\dot{\epsilon }}_{ij}^{\mathrm{p}}$ is plastic strain rate; and $\Delta t$ is time increment.

When rock is isotropic material, the effect of movement tensor is not considered. Yield stress is denoted by Formula (28):

$${\sigma }_{\mathrm{y}}=\left[1+{(\frac{\dot{\epsilon }}{\mathrm{C}})}^{1/\mathrm{P}}\right]({\sigma }_{\mathrm{y}0}+\beta E_{\mathrm{p}}{\epsilon }_{\mathrm{p}}^{\mathrm{eff}})$$

(28)

where $\dot{\epsilon }$ is strain rate; C and P are the constants of the Cowper-Symonds plastic kinetic hardening model; ${\sigma }_{\mathrm{y}0}$ is initial yield stress; ${\epsilon }_{\mathrm{p}}^{\mathrm{eff}}$ is effective plastic deformation, ${\epsilon }_{\mathrm{p}}^{\mathrm{eff}}={\displaystyle {\int }_0^t\mathrm{d}}{\epsilon }_{\mathrm{p}}^{\mathrm{eff}}$; $\mathrm{d}{\epsilon }_{\mathrm{p}}^{\mathrm{eff}}=\sqrt{\frac{2}{3}\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}}$; ${\epsilon }_{ij}^{\mathrm{p}}$ is plastic strain.

3. Numerical Algorithm of the Proposed Damage Model

The proposed blasting damage model of rock was implemented into the commercial software LS-DYNA as a user-defined material model [33], shown in Figure 1. The coupling process between the Cowper-Symonds plastic kinetic hardening model and the proposed damage model can be summarized as follows:

(1) According to ${\sigma }_{ij}^n$ and $\Delta \epsilon$ obtained in the former incremental step, stress in the current incremental step is determined using elastic prediction, where $C_{cjkl}$ is elastic tensor and $\Delta {\epsilon }_{kl}$ is strain increment.

(2) The value of damage variable is updated according to the given algorithm.

(3) Judge whether to yield according to the given yield condition. If rock yields, plastic correction is performed according to the given plastic potential function.

(4) Adjust stress value according to damage value.

In Figure 1, ${\overline{\sigma }}_{ij}^{n+1}$ is effective stress tensor; ${\sigma }_{ij}^{\ast n+1}$ is trial stress tensor; ${\sigma }_{ij}^n$ and ${\sigma }_{ij}^{n+1}$ are stress tensor with damage, respectively; $C_{cjkl}$ is elastic tensor; $\Delta \epsilon$ and $\Delta {\epsilon }_{kl}$ are strain increment; $\Delta C_{\mathrm{d}}$ is crack density increment.

The coupling process between the Cowper-Symonds kinetic hardening model and the proposed blasting damage model (shown in Figure 1) can be summarized as follows.

3.1. Elastic Prediction

According to the plastic increment theory, total strain can be decomposed into the increment of elastic strain and that of plastic strain for any time step (see Formula (29)):

$$\Delta {\epsilon }_{ij}=\Delta {\epsilon }_{ij}^{\mathrm{e}}+\Delta {\epsilon }_{ij}^{\mathrm{p}} (i = 1, 2, 3)$$

(29)

where $\Delta {\epsilon }_{ij}$ is strain increment; $\Delta {\epsilon }_{ij}^{\mathrm{e}}$ is increment of elastic strain; and $\Delta {\epsilon }_{ij}^{\mathrm{p}}$ is increment of plastic strain.

According to Hooke’s law, the relationships among the deviations of elastic stress, hydrostatic stress and strain are denoted by Formula (30):

$$\{\begin{array}{l}{\sigma }_x-{\sigma }_{\mathrm{m}}=2G({\epsilon }_x-{\epsilon }_{\mathrm{m}}), {\tau }_{xy}=G{\gamma }_{xy}\\ {\sigma }_y-{\sigma }_{\mathrm{m}}=2G({\epsilon }_y-{\epsilon }_{\mathrm{m}}), {\tau }_{yz}=G{\gamma }_{yz}\\ {\sigma }_z-{\sigma }_{\mathrm{m}}=2G({\epsilon }_z-{\epsilon }_{\mathrm{m}}), {\tau }_{zx}=G{\gamma }_{zx}\end{array}$$

(30)

where ${\tau }_{xy}$, ${\tau }_{yz}$ and ${\tau }_{zx}$ are shear stress; and ${\gamma }_{xy}$, ${\gamma }_{yz}$ and ${\gamma }_{zx}$ are shear strain.

Assuming that all strains are elastic, Formula (30) is expressed as Formula (31) in incremental form, and Formula (21) is expressed as Formula (32) in incremental form:

$$\{\begin{array}{l}{\sigma }_x^{\ast n+1}={\overline{\sigma }}_x^n+\Delta {\sigma }_{\mathrm{m}}^n+2G(\Delta {\epsilon }_x-\Delta {\epsilon }_{\mathrm{m}}), {\tau }_{xy}^{\ast n+1}={\overline{\tau }}_{xy}+G\Delta {\gamma }_{xy}\\ {\sigma }_y^{\ast n+1}={\overline{\sigma }}_y^n+\Delta {\sigma }_{\mathrm{m}}^n+2G(\Delta {\epsilon }_y-\Delta {\epsilon }_{\mathrm{m}}), {\tau }_{yz}^{\ast n+1}={\overline{\tau }}_{yz}+G\Delta {\gamma }_{yz}\\ {\sigma }_z^{\ast n+1}={\overline{\sigma }}_z^n+\Delta {\sigma }_{\mathrm{m}}^n+2G(\Delta {\epsilon }_z-\Delta {\epsilon }_{\mathrm{m}}), {\tau }_{zx}^{\ast n+1}={\overline{\tau }}_{zx}+G\Delta {\gamma }_{zx}\end{array}$$

(31)

$$\Delta {\sigma }_{\mathrm{m}}=K_{\mathrm{d}}\Delta {\epsilon }_{\mathrm{m}}$$

(32)

where ${\overline{\sigma }}_x^n$, ${\overline{\sigma }}_y^n$, ${\overline{\sigma }}_z^n$, ${\overline{\tau }}_{xy}$, ${\overline{\tau }}_{yz}$ and ${\overline{\tau }}_{zx}$ are stresses at the previous calculation step after being converted to effective stress space; ${\sigma }_x^{\ast n+1}$, ${\sigma }_y^{\ast n+1}$, ${\sigma }_z^{\ast n+1}$, ${\tau }_{xy}^{\ast n+1}$, ${\tau }_{yz}^{\ast n+1}$, ${\tau }_{zx}^{\ast n+1}$ are stresses of elastic test; $\Delta {\sigma }_{\mathrm{m}}^n$ is increment of average stress, $\Delta {\sigma }_{\mathrm{m}}^n=\frac{1}{3}\left(\Delta {\sigma }_x^n+\Delta {\sigma }_y^n+\Delta {\sigma }_z^n\right)$; $\Delta {\sigma }_x^n$, $\Delta {\sigma }_y^n$ and $\Delta {\sigma }_z^n$ are increments of stress; $\Delta {\epsilon }_{\mathrm{m}}$ is increment of average strain, $\Delta {\epsilon }_{\mathrm{m}}=\frac{1}{3}\left(\Delta {\epsilon }_x+\Delta {\epsilon }_y+\Delta {\epsilon }_z\right)$; $\Delta {\epsilon }_x$, $\Delta {\epsilon }_y$ and $\Delta {\epsilon }_z$ are increments of strain; and $\Delta {\gamma }_{xy}$, $\Delta {\gamma }_{yz}$ and $\Delta {\gamma }_{zx}$ are increments of shear strain.

3.2. Calculation of the Damage Variable

According to Formulae (1) and (2), equivalent tensile strain can be obtained by Formula (33):

$$\theta ={\displaystyle \sum _{i=1}^3\ln \left[1+\frac{{\epsilon }_i+\left|{\epsilon }_i\right|}{2}\right]}$$

(33)

Rock damage occurs only when equivalent tensile strain is greater than critical tensile strain. To improve the calculation accuracy of equivalent tensile strain, the equivalent tensile strain is recorded for each load step and transmitted to the next cycle. The calculation algorithm of the damage variable is described as follows.

(1) When ${\theta }_n>{\theta }_{\mathrm{c}}$, tensile damage of element occurs, and is obtained by Formulae (4) and (11):

(2) Calculate the current equivalent tensile strain by Formula (34):

$$\theta =\{\begin{array}{l}\left({\theta }_n+{\theta }_{n+1}\right)/2, {\theta }_{n+1}>{\theta }_{\mathrm{c}}\\ \left({\theta }_n+{\theta }_{\mathrm{c}}\right)/2, {\theta }_{n+1}\le {\theta }_{\mathrm{c}}\end{array}$$

(34)

where ${\theta }_n$ is equivalent tensile strain in the nth load step; and ${\theta }_{\mathrm{c}}$ is critical tensile strain.

(3) Calculate crack density increment $\Delta C_{\mathrm{d}}^{n+1}$ according to Formula (11) and update the history variable of crack density $C_{\mathrm{d}}^{n+1}=C_{\mathrm{d}}^n+\Delta C_{\mathrm{d}}^{n+1}$.

(4) Update damage variable D using Formula (4).

(5) When ${\theta }_n\le {\theta }_{\mathrm{c}}$, do not calculate tensile damage and directly calculate the Cowper-Symonds kinetic hardening model.

3.3. Plastic Correction

If some element is at elastic state in the current calculation step, stress tensor of elastic test ${\overline{\sigma }}_{ij}^{n+1}={\sigma }_{ij}^{\ast n+1}$.

If some element is at plastic state, stress of elastic prediction and plastic hardening function should be echoed to the yield surface by means of radial regression.

The yield function $\varphi$ is calculated using Formula (27). The Cowper-Symonds kinetic hardening model considers the increase of rock yield stress under the condition of high strain rate. Therefore, it is necessary to calculate yield stress ${\sigma }_y^n$ in the current load step using Formula (28).

When rock enters yield state ($\varphi >0$), effective plastic strain is modified using Formula (35):

$${\epsilon }_{\mathrm{p}}^{\mathrm{eff}n+1}={\epsilon }_{\mathrm{p}}^{\mathrm{eff}n}+\Delta {\epsilon }_{\mathrm{p}}^{\mathrm{eff}}={\epsilon }_{\mathrm{p}}^{\mathrm{eff}n}+\frac{{\sigma }_i-{\sigma }_y^n}{3G+E_{\mathrm{p}}}$$

(35)

Finally, deviatoric stress is decreased by Formula (36) and return to new yield surface:

$$S_{ij}^{n+1}=S_{ij}^{\ast }-\frac{3G\Delta {\epsilon }_{\mathrm{p}}^{\mathrm{eff}}}{{\sigma }_i}{S}_{ij}^{\ast }$$

(36)

where $S_{ij}^{\ast }$ is the component of deviatoric stress tensor before decreasing.

3.4. Stress Update

Update stress using Formula (37):

$${\sigma }_{ij}^{n+1}={\overline{\sigma }}_{ij}^{n+1}(1-\delta D)/(1-\delta D_0)$$

(37)

4. Validation of the Availability of the Proposed Model

4.1. Numerical Model of the Uniaxial Compression Test

Numerical simulation of a uniaxial compression test is carried out using the proposed damage model. To verify the availability of the proposed model, numerical simulation results are compared with those calculated by Yang (1996) [7].

The model with length, width, and height all 0.1 m is established to simulate a uniaxial compression test. SOLID 164 is used to mesh the model. The relationship between stress and strain of rock and the development of rock damage variable under the specific constant strain rate are investigated by applying the constant velocity to the node along the Z (vertical) direction at the top of the model.

In the case of uniaxial compression test, compressive stress is ${\sigma }_1$, and equivalent tensile strain is θ denoted by Formula (38):

$$\theta ={\epsilon }_2+{\epsilon }_3$$

(38)

where ${\epsilon }_2={\epsilon }_3=\nu {\sigma }_1/E$.

The main calculation parameters of rock in the uniaxial compression test are listed in

Table 1. Main calculation parameters of rock in uniaxial compression test [7].

E (GPa)	$\mathit{\nu }$	${\mathit{\sigma }}_{\mathbf{0}}$ (MPa)	α	β	${\mathit{\theta }}_{\mathit{c}}$
51.8	0.33	215	3.15 × 106	1.0	0.0002

, which was used in the paper of Yang (1996) [7].

4.2. Rock Blasting Damage Evolution and Softening under Compression Conditions

Figure 2 is the stress-strain curve and damage curve obtained using Yang’s model [7] and it was obtained using the proposed model at a strain rate equal to 100 s⁻¹. The curves in Figure 2 are based only on calculations of Yang’s model and the proposed model.

From Figure 2, it can be noted that the evolution curve of the damage variable with strain is monotonically increasing. With the increase of axial stress, ${\epsilon }_2$ and ${\epsilon }_3$ become tensile strains, and rock damage accumulates gradually to D = 1, when rock is softened significantly until it cannot withstand the applied load. D is about 0.22 when the rock started to soften significantly under stress, seen in Figure 2. This can be taken as a criterion of rock blasting damage and failure.

Compared with the calculation results using Yang’s model [7] and those of the proposed model, it can be concluded that the proposed model can be used to simulate rock damage evolution process under compression, as seen in Figure 2.

4.3. Influence of Strain Rate on Rock Blasting Damage

Since rock deformation characteristics are different for different strain rates, the blasting damage model should have obvious rate correlation. From Figure 3, it can be noted that rock peak stress increases obviously with an increase of strain rate. Rock peak stress increases about 650% as the strain rate varies from 10 to 1000 s⁻¹, that is, there is a significant rate modification phenomenon.

Since the time-dependent damage variable D is related to strain rate, the proposed model can reflect that the higher the strain rate, the more obvious the lag phenomenon of rock softening.

4.4. Influence of Damage Modification Coefficient on Stress-Strain Relationship

From Figure 4, it can be noted that when $\delta$ is reduced from 1.0 to 0.5, rock residual strength gradually increases. Rock residual strength is 66.8% of peak strength when $\delta$ = 0.5, but only 9.59% of peak strength when $\delta$ = 1.0. That is, rock residual strength is very sensitive to the damage modification coefficient. Damage modification coefficient can successfully reflect the residual strength characteristics of rock.

4.5. Influence of Initial Damage on Rock Stress-Strain Relationship

From Figure 5, it can be noted that the stress-strain curve of rock changes obviously with the increase of its initial damage. When strain is less than 0.001, stress-strain curves present an approximately linear change. When strain is more than 0.001, rocks with greater initial damage will soften earlier. When the initial damage degree of rock varies from 0.05 to 0.2, its stress peak value is from 26% to 39% and significantly lower than that of initial undamaged rock.

5. Case Study

Taking the tunnel under-passing sewage box culvert of Shenzhen’s Metro, China, as an example, blasting damage range, blasting damage degree of tunnel surrounding rock, and dangerous parts of box culvert structure are analyzed using the proposed blasting damage model. In addition, the safe distance between the top of the tunnel and the existing box culvert in the tunnel blasting construction is investigated. The tunnel is excavated using the drill-blasting method.

5.1. Calculation Model

Surrounding rocks of the tunnel and sewage box culvert are medium weathered granite in the calculation model. Cross-section dimensions of the calculation model are shown in Figure 6.

Main calculation parameters of rock are listed in

Table 2. Main calculation parameters of rock [7,9].

$\mathit{\rho }$/kg/m3	E/GPa	$\mathit{\nu }$	${\mathit{\sigma }}_{\mathbf{0}}$/MPa	${\mathit{E}}_{\mathbf{tan}}$/GPa	C /s−1	P	α	β	$\mathit{\delta }$
2600	30	0.27	100	50	2.5	4.0	3.15 × 106	1.0	0.9

. In Table 2, $\rho$ is the density of surrounding rock; E is Young’s modulus; $\nu$ is Poisson’s ratio; ${\sigma }_0$ is yield stress; $E_{\tan }$ is tangent modulus; C and P are the constants of Cowper-Symonds plastic kinetic hardening model; α and β are hardening parameters of rock; and $\delta$ is damage modification coefficient of rock. $\rho$, E, $\nu$, ${\sigma }_0$, $E_{\tan }$, and $\delta$ were determined by laboratory tests; C and P are obtained from Xie et al., 2016 [12]; α and β are obtained from Yang et al., 1996 [9].

The length of the model along the tunnel longitudinally is 45.0 m; the length of the model is 55.0 m along the horizontal direction; the height of the model is 35.2~40.0 m. The thickness of the soil overlying the box culvert is 5.0 m. The spacing between the box culvert and the tunnel is 1.0~2.0 m.

The model consists of 4360 elements and 5301 nodes. According to the site conditions, two groups of vertical planes of the model are fixed at the horizontal direction, and the base plane (xy plane) is fixed at the z direction. The gravity force is applied to take the effect of the self-weight stress field into account.

5.2. Blasting Load and Its Application Method

The purpose of tunnel blasting numerical simulation in this paper is to obtain the damage range of surrounding rock caused by blasting construction, and it does not discuss the crushing mechanism and fracture effect of rock near blast-holes. Therefore, it is reasonable and feasible to simulate the blasting effect by applying simplified equivalent loads.

In tunnel blasting numerical simulation, the shapes of blast-holes are not reflected, and blasting pressure load is equivalently applied to the axis plane of blast-holes through conversion [30]. And explosion pressure load is approximated by triangular explosion load, as shown in Figure 7. In Figure 7, $P_{\max }$ is peak pressure; $t_{\mathsf{\tau }}$ is the time corresponding to the rise of load; $t_{\mathrm{s}}$ is total load time.

$P_{\max }$ can be calculated by Formula (39) [34]:

$$P_{\max }=\frac{1}{8}{\rho }_{\mathrm{e}}{V}^2{k}_{\mathrm{d}}^6\eta$$

(39)

where ${\rho }_{\mathrm{e}}$ is charge density; V is explosive detonation velocity; $k_{\mathrm{d}}$ is borehole charge decoupling coefficient, $k_{\mathrm{d}}=\frac{d_{\mathrm{b}}}{d_{\mathrm{c}}}$, $d_{\mathrm{b}}$ and $d_{\mathrm{c}}$ are diameter of blast-hole and diameter of equivalent explosive coil, respectively; $\eta$ is the multiple of pressure increase value when detonating gas hits blast-hole wall, $\eta =8\sim 11$.

Based on the Saint-Venant principle in mechanics, the equivalent blasting load time-history curve is applied to the plane determined by the central line of blast-hole and the axis of blast-holes in the same row (see Figure 8). Equivalent pressure is determined according to the above principle, and is calculated by Formula (40) [34].

$$P_{\mathrm{e}}=\frac{2r_0}{a}{P}_0$$

(40)

where $P_{\mathrm{e}}$ is equivalent pressure; $P_0$ is blasting force acting on single blast-hole wall, $P_0=P_{\max }$; $r_0$ is blast-hole radius; and a is blast-hole spacing.

$t_{\mathsf{\tau }}$ can be calculated by Formula (41):

$$t_{\mathsf{\tau }}=\frac{L_1}{V}$$

(41)

where $L_1$ is charge length.

In the paper, $t_{\mathsf{\tau }}$ = 170 μs, $t_{\mathrm{s}}$ = 1000 μs [35]. Explosive parameters are listed in

Table 3. Explosive parameters.

Type of Explosive	${\mathit{\rho }}_{\mathbf{e}}$ (kg/m3)	V (m/s)
2# rock emulsion explosive	950	3500

. The diameter of blast-holes is 42 mm; the diameter of the charges is 32 mm.

According to the acoustic wave velocity of rock, $D_0$ = 0.375 is obtained by Formula (7) in this study. Combined with our experience and considering the safety of the box culvert, the threshold value of rock damage corresponding to the edge of the blasting damage impact range can be determined by Formula (42):

$$D_{\lim }=0.20+0.80D_0$$

(42)

where $D_{\lim }$ is the threshold value of rock blasting damage; $D_0$ is initial damage of rock.

Formula (42) can be used to determine the damage accumulation range of rock caused by blasting. According to Formula (42), the blasting damage threshold of surrounding rock around the box culvert structure is 0.50. In the process of tunnel blasting construction, the bottom of the box culvert is obviously the most dangerous. Therefore, the surrounding rock elements close to the bottom of the box culvert structure were selected for the blasting damage analysis, as shown in Figure 9.

5.3. Results and Discussion

Figure 10a is a three-dimensional blasting damage diagram of surrounding rock when d = 1 m. Figure 10b is cross-section diagram of surrounding rock blasting damage when d = 1 m. d is the spacing between the box culvert and the tunnel.

From Figure 10a,b, it can be noted that the thickness of blasting-damaged rock beneath the box culvert with $D>0.9$ is approximately 0.5 m; the thickness of blasting-damaged rock overlying the tunnel with $D>0.8$ is approximately 20~22 cm.

From Figure 10a, it can be noted that the thickness of the blasting-damaged rock between the box culvert and the tunnel with $D<0.7$ is less than 30 cm. From Figure 10b, it can be noted that the thickness of the surrounding rock between the box culvert and the tunnel with $D<0.7$ is approximately 20 cm; the thickness of the surrounding rock between the box culvert and the tunnel with $D\ge 0.5$ is 1.0 m.

Figure 11a is three-dimensional blasting damage diagram of surrounding rock when d = 2 m. Figure 11b is cross-section diagram of surrounding rock blasting damage when d = 2 m.

From Figure 11a,b, it can be noted that the thickness of blasting-damaged surrounding rock beneath the box culvert with $D>0.9$, is less than 15 cm, located from 1.2 m to 1.35 m beneath the box culvert. The thickness of blasting-damaged surrounding rock beneath the box culvert with $D>0.7$, is less than 45 cm, located from 0.9 m to 1.35 m beneath the box culvert (see in Figure 11a). The thickness of blasting-damaged surrounding rock overlying the tunnel with $D<0.6$, is approximately 60 cm (see in Figure 11a). The thickness of blasting-damaged surrounding rock beneath the box culvert with $D\le 0.5$, is more than 55 cm (see in Figure 11a).

From Figure 11b, it can be noted that the thickness of blasting-damaged surrounding rock beneath the box culvert with $D>0.8$, is less than 20 cm, located from 1.2 m to 1.4 m beneath the box culvert. The thickness of blasting-damaged surrounding rock overlying the tunnel with $D\le 0.5$, is more than 55 cm (see in Figure 11b).

Comparing Figure 10a with Figure 11a, it can be found that the blasting damage range of surrounding rock and blasting damage degree of surrounding rock both very significantly reduce as the spacing between the box culvert and the tunnel increases from 1.0 m to 2.0 m. In addition, blasting damage of surrounding rock exists the concentrated phenomenon beneath the box culvert when the spacing between the box culvert and the tunnel is 1.0 m.

Figure 12a is the three-dimensional blasting damage diagram of surrounding rock when d = 3.0 m. Figure 12b is a cross-section diagram of surrounding rock blasting damage when d = 3.0 m.

From Figure 12a,b, it can be noted that the thickness of blasting-damaged surrounding rock beneath the box culvert with $D<0.5$, is more than 2.35 m. The thickness of blasting-damaged surrounding rock beneath the box culvert with $D=0.5\sim 0.7$ is 0.65 m, located from 0.65 m to 1.3 m overlying the tunnel (see in Figure 12a,b).

Figure 13a is three-dimensional blasting damage diagram of surrounding rock when d = 4.0 m. Figure 13b is cross-section diagram of surrounding rock blasting damage when d = 4.0 m.

From Figure 13a,b, it can be noted that the thickness of blasting-damaged surrounding rock beneath the box culvert with $D<0.5$, is more than 3.7 m. The thickness of blasting-damaged surrounding rock beneath the box culvert with $D=0.5\sim 0.7$ is 0.3 m, located from 0.45 m to 0.75 m overlying the tunnel (see in Figure 13a,b).

Comparing Figure 10 with Figure 11, Figure 12 and Figure 13, it can be found that when the tunnel underpasses the box culvert at a short distance (d ≤ 2.0 m), the blasting damage of surrounding rock between the box culvert and the tunnel is highly related to d; especially when d ≤ 2.0 m, the box culvert structure has a significant “amplification effect” on the surrounding rock blasting damage, and the amplification effect attenuates with the increase of d; when d ≥ 4.0 m, it can be considered that the box culvert structure has very little influence on the surrounding rock blasting damage of the tunnel.

The “amplification effect” phenomenon on surrounding rock blasting damage may be caused by the facts as follows: the inner space of the box culvert may create an additional free face for detonation holes if the spacing between the box culvert and the tunnel is too small (d ≤ 2.0 m). Furthermore, the compressive wave generated by blasting reaches the box culvert and is reflected to form the tensile wave, which acts together with the subsequent stress wave to accelerate the damage development speed of surrounding rock in this zone.

From Figure 14, it can be noted that when the distance between the box culvert and the tunnel is 3~4 m, the blasting damage of surrounding rock beneath the box culvert does not exceed $D_{\lim }$. Thence, when d > 3.0 m, the surrounding rock around the box culvert are essentially not damaged by blasting. From Figure 15, it can be noted that when the blasting damage of surrounding rock beneath the box culvert shows an upwardly concave parabolic change with d.

From Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15, it can be suggested that the safe distance between the tunnel and the existing box culvert in blasting construction should be not less than 4.0 m. In addition, it can be suggested that the pre-reinforcement of surrounding rock overlying the tunnel should be implemented before tunnel blasting when d < 2.0 m. The sewage box culvert was successfully under-passed using the pre-reinforcement of pipe sheds before tunnel blasting, together with laying out one empty hole without explosives between every two blasting holes with explosives for the peripheral holes of tunnel.

6. Conclusions

• The Cowper-Symonds kinetic hardening model improved using initial damage degree and damage modification coefficient considering rock residual strength, and can be used to simulate the damage evolution process of rock and softening of rock caused by damage growth under compression.
• Blasting damage range of surrounding rock with D ≥ 0.5 very significantly reduces from 1.0 m to 0.3 m as the spacing between the box culvert and the tunnel increases from 1.0 m to 4.0 m.
• Blasting damage of surrounding rock exists some concentrated phenomena beneath the box culvert when the spacing between the box culvert and the tunnel is less than 2.0 m.
• Using the proposed model can reasonably determine the safe distance between the tunnel and the existing box culvert in blasting construction, which should be not less than 4.0 m.
• This paper lacks comparative analysis and validation with other methods such as the model using rock tension-compression coupling damage algorithm. Further studies are required for these deficiencies as well as the validation of numerical simulation by field damage detection.
• Notwithstanding its limitation, the simulation using the proposed model conducted in the paper is still very significant as a reference to guiding tunnel drill-blasting construction closely under-passing sewage box culverts.