Modeling and Simulation of Dynamic Unloading of Prestressed Rockmass

Liang Wu,Xiaorui Xiang,Yang Chen,Karrech Ali ,Junru Zhou, and Ming Chen

1 College of Science,Wuhan University of Science and Technology,Wuhan,430065,China.

2 School of Civil and Resource Engineering,the University of Western Australia,35 Stirling Hwy,Crawley WA 6009,Australia.

3 Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering Ministry of Education,Wuhan University,Wuhan,430072,China.

Abstract:During the excavation of deep rock,a sudden change in boundary conditions will cause the in-situ stress on the excavation surface to release instantaneously.This disturbance propagates in the form of an unloading stress wave,which will enlarge the damage field of surrounding rock.In this paper,the dynamic unloading problem of the insitu stress in deep rock excavation is studied using theoretical,numerical,and experimental methods.First,the dynamic unloading process of rock is analyzed through adopting the wave equation,and the equivalent viscous damping coefficient of the material is taken into consideration.Calculations show that there is significant tensile strain in the rock bar when the strain rate is above 10-1 s-1.With an increase in the length or damping coefficient,the wave state will change from an underdamped to an overdamped state.Second,implicit and explicit solvers of the finite element method are employed to simulate rock unloading processes,which can be used to verify the theoretical results from one-dimensional to three-dimensional stress states.Finally,the dynamic unloading experiment of a onedimensional bar is used to further verify the validity and accuracy of the theoretical analysis.

Keywords:Deep rock mass,dynamic unloading,damping,strain rate,wave equation.

1 Introduction

At present,drilling and blasting methods are widely used in deep rock excavation projects.Under the conditions of drilling and blasting,the stress on the excavation surface is released within a few milliseconds,accompanied by explosive detonation,cracking and fragmentation of rock mass and the formation of excavation free surface.Under high insitu stress conditions,excavation will induce strong stress adjustment in the rock mass.The disturbance propagates rapidly from the excavation boundary to the interior of rock mass in the form of unloading wave,and causes the rock mass to loosen or deteriorate in the excavation disturbance area.[Chen,Lu,Zhou et al.(2008);Lu,Peng,Xu et al.(2007);Lu,Yang,Yan et al.(2012)].In addition,this process may cause a variety of geological engineering disasters (such as rock bursts and sudden deformations),which will affect the progress of a project and even threaten the safety of workers and equipment [Galybin and Odintsev (1992);Jiang,Liu,Zhao et al.(2005)].The study on the mechanical mechanism and effect of dynamic unloading caused by excavation is of great significance to the stability and safety evaluation of underground excavation engineering.

The systematic study on the mechanism and effect of dynamic unloading of rock began in the 1970s.In 1966,Cook et al.[Clay,Cook,Cook et al.(1966)] discovered that the sudden release of in-situ stress during the excavation of rock resulted in overrelaxation of the rock and generated tensile stress in the rock mass.Recently,many researchers have engaged in research in this field.Bauch et al.[Bauch and Lempp (2004)] carried out an indoor test on German red sandstone and found out that when the confining pressure is constant and the axial pressure is suddenly unloaded,the rock member has a rapid cracking failure,and the stored strain energy is released quickly.Lu et al.[Lu,Yang,Yan et al.(2012)] studied the process of release of in-situ stress induced by rock mass excavation with the method of drill and blast and determining the duration of the release of in-situ stress,it is found that the release of in-situ stress induced strain rate can reach a magnitude of 10-1-101/s or higher if the initial in-situ stress has a level of 20-50 MPa.In the review of dynamic testing techniques and dynamic mechanical properties of rock materials,Zhang et al.[Zhang and Zhao (2014)] explored various dynamic mechanical properties and corresponding fracture behaviors of rock.Li et al.[Cao,Li,Tao et al.(2016)] adopted a numerical model to simulate the dynamic responses around an existing tunnel under unloading disturbance forces.Fan et al.[Fan,Lu,Yan et al.(2015)] analyzed the changes induced by quasi-static unloading of in-situ stress and the transient release of in-situ stress for the case of a circular excavation under the conditions of in-situ hydrostatic stress.Tao et al.[Li,Li,Tao et al.(2013)] adopted a finite element program (LS-DYNA) to investigate the unloading failure mechanisms of hard rock in a confined state from strain energy density rates.Li et al.[Cao,Li,Zhou et al.(2014)] used the particle flow code PFC2D to analyze the unloading mechanisms of brittle rock under different stress paths.However,many uncertainties related to the stress unloading process still need to be investigated.

In this paper,the combination of theory,numerical and experimental methods is used to study the dynamic unloading problem of in-situ stress in the excavation of deep rock mass with the direction perpendicular to the excavation face.Firstly,taking the blasting excavation of floor of hydropower station as the engineering background of onedimensional rock dynamic unloading mechanics model,the dynamic unloading process of rock is analyzed by wave equation,and the equivalent viscous damping coefficient of material is considered.Secondly,the separation of variables method and the Duhamel principle are adopted to solve the stress wave equation,and the dynamic strain equation for a rock bar in the unloading process is obtained.Additionally,the influence of the rock bar length,damping coefficient,and stress unloading path on the dynamic unloading process is discussed.The effects of unloading strain rate on failure of the rock mass under different unloading time and initial stress values is also analyzed.Then,considering the confining pressure environment,the one-dimensional wave equation is extended to threedimensional wave equation and the numerical solution is calculated by the difference method.The dynamic strain data under different confining pressure are analyzed,and the relationship between confining pressure and one-dimensional unloading model is established.Finally,the finite element method is used to simulate the penetrating excavation process of the floor of the underground powerhouse of the hydropower station,which verifies the applicability of the theoretical model in relevant projects.The correctness and accuracy of the theoretical solution are further verified by the dynamic unloading test of rock bar.

2 Theoretical model of dynamic unloading of prestressed rock mass

2.1 Establishment of mechanical model

In the excavation of high side wall and floor of underground powerhouse of hydropower station,the direction of initial in-situ stress release is mainly perpendicular to the free surface of excavation.Therefore,the initial in-situ stress unloading problem for threedimensional rock excavation can be simplified in theory to a one-dimensional rock bar.As shown in Fig.1,the side length of the floor blasting area of underground powerhouse of Baihetan Hydropower Station is more than 25 m.The floor of underground powerhouse adopts cylindrical borehole charging structure and one-time through blasting in excavation area.The one-dimensional theoretical model shown in Fig.2 is a rock bar of length L with one fixed end and one free end.Initially,under the effect of the load σ (0),the free end of the bar generates a displacement b and remains static.Starting from t=0,a new free surface is formed due to excavation,and the initial stress is gradually removed in accordance with the function σ (t).Influenced by many factors,the dynamic damping behavior is quite complex.To investigate the influences of damping forces on the dynamic stress and strain of the rock bar,the rock damping is simplified to a linear viscous damped system,which is mathematically easier to address.

Figure 1:The floor blasting area of underground powerhouse of Baihetan Hydropower Station

Figure 2:One-dimensional mechanical model of initial stress unloading

The governing equation for the motion of particles in a rock bar is the one-dimensional wave equation [Arfken,Ruby and Weber (1999)].

where b=σ(0)L/E.

The displacement of the fixed end of the bar is always 0,and this boundary condition can be described by the Dirichlet condition:

The free end of the bar is subjected to the stress σ(t) and can be described by the Neumann condition:

2.2 Equation solution

This is a problem with nonhomogeneous boundary conditions that can be solved.First,the equations to resolve this problem are:

Therefore,differential equations are obtained by substituting (6) into (1)～(5):

Using the superposition principle,Eq.(7) can be simplified to solve the following two problems:

This leads to the solution for the original problem:

The separation of variables method can be used to solve (I).This is shown by the equation below:

Substituting (11) into (8):

Therefore:

And:

Eq.(13) is often referred to as the intrinsic value problem or the Sturm-Liouville problem [Arfken,Ruby and Weber (1999)].In combination with the boundary conditions in (8),only whenλ＞0,X(x) has a nonzero solution,whose intrinsic value is:

The corresponding intrinsic function is:

Substituting (15) for (14):

Eq.(17) is a second-order constant coefficient homogeneous linear differential equation.The form of the solution is related to the roots of the characteristic equation (t2+ζt+a2λk=0).We define ?k=ζ2-4a2λkas the form of the corresponding solution,which has the following three cases according to the positive and negative values of ?k[Department of Applied Mathematics (2014)],as in Tab.1.

Table 1:General solutions corresponding to different characteristic roots

The three general solutions are similar to the damped free vibration solutions of linear systems [Pahud and Pedro (1991)];?k＜0 ,?k＞0,and ?k=0 correspond to the underdamped,over-damped,and critical-damped states of a linear system,respectively.

The case ?k＜0 is solved here (other cases are similar,and it is only necessary to substitute the characteristic root into the general solution corresponding to Tab.1;the steps of the solution are consistent).The corresponding characteristic root is:

Then:

Combining equations (11),(16) and (19) results in:

Using the two initial conditions in (8),the following solutions are obtained:

Thus:

We then solve (II).Duhamel’s principle (the homogeneous impulse principle) is used [Tan and Zhong (2007)].If the functionφ(x,t,τ) is the solution of the following equations,a new system,(II),can be written:

This results in:

16. Unhappy girl that I am: The heroine s repentance58 comes because of her poor status and, perhaps, actually seeing the wealth of King Thrushbeard makes her realize what she refused. Von Franz states that such regret is typical of the animus driven woman, and is pseudo guilt59 and sterile60 (173).Return to place in story.

This is the solution of equations (II),and thus,only problem (III) needs to be solved.Set t′=t-τ;the following can be obtained from Eq.(23):

The method of solving Eq.(25) is the same as that for solving (I):

By combining the two initial conditions in (25),the following solutions are obtained:

Substituting t′=t-τ into (26) and (27),and combining (24),we obtain:

In summary,the solution of the original Eq.(10),can be obtained by combining (22) and (28).Additionally,the dynamic strain at different locations of the rock bar can be defined by

3 Dynamic unloading mechanism of rock mass

3.1 Propagation process of the unloading wave

During the excavation of a deep rock mass,a sudden change in boundary conditions causes the initial stress on the excavation surface to be removed instantaneously.This disturbance propagates in the form of a stress wave (unloading wave) in the rock mass.Its wavelength isl=a?ts,wheretsis the length of the unloading time.According to the one-dimensional wave model,the time required for unloading waves to propagate from one end of the bar to the other ista=L/a=0.27 ms if the length of rock bar is 1.2 m.The unloading wave is reflected at the boundary,and the reflected wave propagates in the direction opposite to that of the incident wave;a superposition zone of stress waves is formed,and the length of the superimposed zone isl/2.According to the reflection theory of elastic waves at the fixed end [Wang (2005)],since the displacement of the fixed end is always 0,the properties and strength of the reflected wave are exactly the same as those of the incident wave,and the intensity of the wave in the superimposed zone doubles.However,the properties of the reflected wave at the free surface are opposite to those of the incident wave,and the strength in the superimposed zone weakens,so the free end stress after unloading is always zero.

Figure 3:Propagation process of the unloading wave (the tensile wave is the bar and the compression wave is below)

Whents＜ta,the unloading wavelength is shorter than the bar;if the unloading wave is a linear unloading rectangular wave,the stress state and its change at each moment of the dynamic unloading process are shown in Fig.3.The dynamic strain curve at different distances from the fixed end of the rock bar is shown in Fig.4.In this figure,A corresponds to Process 1,and the unloading wave (tensile wave) travels from the free end to the fixed end;when the unloading wave passes,the initial compressive strain attenuates to 0 in accordance with the unloading curve.B corresponds to Processes 2-4,where the unloading wave is reflected at the fixed end,the reflected and incident waves are superimposed,and the wave intensity increases.The 0 m and 0.3 m positions are both in the left superimposed zone,so the slope (absolute value) of the dynamic strain curve increases.Additionally,the wave at the end is immediately superimposed when the reflected wave is generated until the incident wave is fully reflected.The wave at a certain distance from the end of the bar is superimposed only between the time that reflection occurs and the time that the incident wave is completely reflected.C corresponds to Process 5,where the reflected wave propagates to the right,as the front unloading wave has reduced the compressive strain to 0 everywhere,and the reflected tensile wave causes rod strain when passing.D corresponds to Processes 6～8;the unloading wave is reflected at the free end,and the reflected wave,which is a compressional wave,is opposite to the incident wave,and the resultant intensity is 0 after superimposition.Similar to B,the wave at the 1.2 m position (free end) is superimposed immediately when the reflection is generated until the completion of the reflection.Therefore,after unloading,the free end surface is in a state of no stress.The wave at a certain distance from the end of the bar is superimposed only between the time that reflection occurs and the time that the incident wave is completely reflected,so the tensile strain slightly decreases compared with the non-superimposed area.E corresponds to Process 9.The reflected compressional wave travels to the fixed end and reduces the tensile strain around the bar to 0.Fcorresponds to Processes 10～12,where the compressional wave is reflected and superimposed at the fixed end,increasing the slope (absolute value) of the curve in the superimposed area.G corresponds to Process 13;the reflected compressional wave propagating to the free end causes compressive strain in the bar.H corresponds to Processes 14～16,where the compressional wave is reflected and changes into a tensile wave at the free end.I corresponds to Process 1.……In this manner,the unloading wave circulates along the bar.Because of the existence of damping,the propagation process of the unloading wave is complete when the initial strain energy is completely dissipated.In the process of the transition between compressive strain and tensile strain,the fixed end plays an important role in the reflection of the unloading wave.

Figure 4:Dynamic strain curves at different locations on rock bar (ts=0.2 ms)

In the above analysis,the unloading wave propagation process is based on the conditionts＜ta;i.e.,the wave lengthlis shorter than the bar lengthL.With an increase ints,lis gradually longer thanL.Then,there are two differences from previous analysis.① Since the superimposed zone is longer thanL/2,the middle section of the bar will be in both the left and the right superimposed zones.② Becauselis longer thanL,reflection at the left and the right ends may occur simultaneously.Whents≥ta,the superimposed area length is longer thanL,and each position on the bar is in the superimposed area.At this time,the propagation of the unloading wave and its effects on the bar become very complicated.

3.2 .The effect of unloading strain rate

According to the above theoretical analysis,if the rock lengthL,the basic mechanical parameters of the rock material,the material damping,the initial stress,and the unloading function are all known,the dynamic strain of the rock bar during the unloading process can be obtained.An example calculation is given here:the length of granite rock bar is 1.2 m,and the basic mechanical parameters of the rock sample are measured as shown in Tab.2.

Table 2:Basic mechanical parameters of rock sample

Accurate damping coefficients of rock materials need to be measured by specific experimental methods.Considering the universality of the theory,this paper,based on the experimental data given in the literature [Nie,Xu and Ren (2011)],uses the damping factorc=5 kN?s?mm-1for analysis.

There are three typical forms of the initial stress unloading functionσ(t),namely,linear attenuation,trigonometric function attenuation,and exponential attenuation:

The initial stress is attenuated from the zero moment and,at thetsmoment,is attenuated to 0.As is shown in Fig.5,curves 1,2 and 3 correspond to linear function attenuation,trigonometric function attenuation and exponential function attenuation,respectively:

Figure 5:Unloading curve of initial stress

In the unloading process,the average strain rate is:

3.2.1 Calculation results for different unloading time durations

If the initial stress isσ0=10 MPa and the unloading time is 0.5 ms,1.0 ms,2.0 ms,4.0 ms and 8 ms,the corresponding average strain rates are 0.4 s-1,0.2 s-1,0.1 s-1,0.05 s-1and 0.025 s-1,respectively.Figs.6(a),6(b) and 6(c) are the dynamic strain curves at a distance of 0.4 m from the fixed end and correspond to the cases of exponential attenuation,linear attenuation and trigonometric function attenuation,respectively.

Figure 6:Dynamic strain curve of bar under different strain rate (initial stress fixed,unloading time changed)

3.2.2 Calculation results for different initial stress levels

If the unloading time is 2 ms and the initial stresses are 40 MPa,20 MPa,10 MPa,5 MPa,and 2.5 MPa,the corresponding average strain rates are 0.4 s-1,0.2 s-1,0.1 s-1,0.05 s-1and 0.025 s-1,respectively.Figs.7(a),7(b) and 7(c) show the dynamic strain curves at a distance of 0.4 m from the fixed end and correspond to the cases of exponential attenuation,linear attenuation and trigonometric function attenuation,respectively.

Figure 7:Dynamic strain curve of bar under different strain rate (unloading time fixed,initial stress changed)

From the above calculations,it can be seen that the initial strain is converted into tensile strain during the initial stress unloading process,and higher strain rates produce higher tensile strain peaks.It is generally believed that,when the strain rate is in the range of 10-5s-1～10-1s-1,a quasi-static process exists,and when the strain rate is greater than 10-1s-1,a dynamic process exists [Li (2014)].When the strain rate is greater than 10-1s-1,the peak value of the tensile strain is clear;when the strain rate is less than 10-1s-1,i.e.,when the unloading belongs to a quasi-static process,the peak value of the tensile strain is smaller and no tensile strain occurs.As the tensile strength of rock is much lower than its compressive strength,the tensile stress generated during unloading may exceed the dynamic tensile strength of the rock and thus result in damage.

3.2.3 Comparative analysis of unloading curves

From the above analysis,it can be seen that,under the same initial stress and unloading time,the dynamic strain fluctuation and attenuation law of the three unloading modes are consistent.However,the peak values of tensile strain under different unloading modes are quite different.In Fig.6,the peak value of tensile strain under exponential attenuation is the highest,and that under linear attenuation is the lowest.This is consistent with the conclusion of a previous study [Lu,Yan,Zhao et al.(2016)].Higher unloading rates produce stronger dynamic unloading effects,while the unloading rates of the first half of the exponential unloading function and the latter half of the trigonometric unloading function are higher than that of the linear unloading function.

In Fig.7,the minimum strain peak value is also generated by linear attenuation,and the peak values of the trigonometric function and the exponential attenuation function are similar,while that of the trigonometric attenuation function is slightly higher than that of the exponential attenuation function.By further analyzing Fig.7,the reasons leading to such results are as follows:

For the first half of the exponential attenuation curve,the unloading rate is high,and the unloading rate in the second half is very low (even close to quasi-static unloading),so the dynamic effect caused by exponential unloading is mainly concentrated in the first half.For example,for the curve with a strain rate of 0.4 s-1in Fig.7(a),the first peak of tensile strain appears at 0.7 ms;this is mainly caused by the first half of the exponential unloading curve,as the initial stress has not attenuated to 0 at this time and is still greater than 6.26 MPa (the free end is at 6.26 MPa,and the fixed end stress value is slightly greater).The superimposition of the residual initial compressive stress and the tensile stress generated by unloading leads to a decrease in the peak value of tensile strain.

The unloading rate of the trigonometric function is higher in the second half,and the dynamic effect is mainly concentrated in the second half.Figs.7(b) and 7(c) show that the first peak value of the tensile stress caused by both unloading of the trigonometric function and linear unloading occurs after the initial stress of the unloading is reduced to 0,and the peak value of the tensile strain does not decrease because of the superimposition.This can be deduced from an analysis of the propagation process of the unloading wave in the bar.This can be deduced from an analysis of the propagation process of the unloading wave in the bar,as discussed in the previous chapter.

3.3 Other factors affecting the tensile strain of rock masses

In addition to the strain rate and unloading function,rock lengths and damping parameters also affect the tensile strain of a rock mass.

3.3.1 Effect of bar length

The elastic wave theory can clearly describe the propagation process of an unloading wave in a rock bar and the mechanism of formation of tensile strain.The analysis shows that tensile strain is generated only when the unloading wave is reflected from the fixed end.If the bar is infinite,i.e.,there is no fixed end,there will be no tensile strain in the bar.In addition,because of the existence of damping,when the unloading wave propagates a certain distance along the infinite bar,the wave intensity decreases gradually and eventually degenerates into quasi-static unloading.A theoretical explanation can be obtained from the solution of the wave equation.When solving Eq.(17),according to the different values of the characteristic equation roots,three different solutions can be obtained,which are similar to the three kinds of damping states of the damped free vibration of the linear system.For these three conditions,when the value of ?k=ζ2-4a2λkis less than,greater than,and equal to 0,three states can be described:underdamped fluctuation,overdamped fluctuation and critically damped fluctuation,respectively.?kis related not only to the damping coefficient and the material (whereais the wave velocity,which depends on the modulus and density of the material) but also to λk.By formula (15),λkincreases with an increase ink.Thus,regardless of changes in other parameters,as long askis sufficiently large,?k＜0 will always be the case.In addition,λkdecreases with an increase in the length of the rock barL.When the damping coefficient and material parameters are constant (i.e.,under uniform rock conditions),there are a limited number of ?k＞0 cases when the rod length grows,and at this point,the solution of Eq.(17) represents an overdamped state.The final solution to the transient unloading problem of a rock mass takes the form of an infinite series (superimposition of various modes),and every item (each order mode) of the series has the composition of the solution to the Eq.(17).Thus,the longer the rock bar is,the larger the number of overdamped states included in the solution and the closer the structure fluctuation is to the overdamped fluctuation state.It is well known that the underdamped vibration of a linear system is a reciprocating attenuation vibration,the overdamped vibration is nonreciprocating attenuation vibration,and the onedimensional fluctuation in the transient unloading problem is similar.The strain in the underdamped fluctuating rock bar alternates between compressive strain and tensile strain,while an overdamped wave will lose this characteristic.

Figure 8:Dynamic strain curves at fixed end of bars with different lengths

Fig.8 shows the dynamic strain curves near the fixed end of bars of differing lengths (under linear unloading conditions,the strain rate is 0.2 s-1and the damping coefficient isc=5 kN?s?mm-1).The results show that,with an increase in the bar length,the tensile strain peak decreases gradually,and the number of wave cycles decreases correspondingly.When calculating ?kwith different bar lengths,we find that,when the length ofLranges from 1.2 m to 6.0 m,the value ofkis arbitrary and ?kis less than 0;i.e.,all the solutions are underdamped fluctuations,so tensile strain is generated in the curve.When the bar length is 7.2 m,?1＞0 and starting fromk=2,?kis less than 0,and the first solution (i.e.,the first order modal damping) is an over-damped fluctuation state.Therefore,when the length of the bar is 7.2 m,the dynamic strain curve is directly attenuated to 0,no tensile strain occurs,and the dynamic strain curve shows no fluctuation.

3.3.2 Effect of material damping

According to the literature [Nie,Xu and Ren (2011)],the viscous damping coefficient is given by c = 5 kN?s?mm-1,7.5 kN?s?mm-1,10 kN?s?mm-1,12.5 kN?s?mm-1and 15.0 kN?s?mm-1.When the bar length is 2.4 m,the dynamic strain curves at the fixed end are shown in Fig.9.With an increase in damping,the strain rate at the fixed end decreases slightly,and the tensile strain peak decreases gradually.In addition,the damping coefficient corresponds to different conditions as well.For the case ofc=15 kN?s?mm-1,?1＜0,the dynamic strain curve of the rock bar is directly attenuated to 0,no tensional strain occurs,and the curve is free of fluctuations.The result is similar to those for different bar lengths.

Figure 9:Dynamic strain curve at fixed end of bars under different damping factors

3.4 Dynamic unloading effect under three-dimensional stress state

One-dimensional wave theory analysis clearly illustrates the propagation process of unloading wave in rock rod,and fully explains that the initial compressive strain (force) will be transformed into tensile strain (force) during dynamic unloading.However,the above theoretical analysis is based on one-dimensional wave equation,without considering the three-dimensional stress state of rock mass.In order to further reveal the unloading effect of rock mass under complex stress conditions,the governing equation (lame equation) of particle motion in rock mass under three-dimensional Cartesian coordinates is established based on wave equation and considering the confining pressure environment of rock mass:

In the formula,ρis density,u、vandware particle displacement components in x,y and Z directions,tis time,λ and μ are Lame constants,?2are Laplace operators,and X,Y and Z are force density components in x,y and z directions respectively.It is difficult to get the theoretical solution of the above equation by analytic method.Therefore,it is advisable to adopt numerical method.When the boundary and initial conditions are known,the finite difference method can be used to solve (31) numerical solutions by combining the geometric and physical equations in elasticity.

While rock masses encountered in engineering are generally under triaxial stress conditions,the one-dimensional wave function theory should be improved to solve practical problems.As is shown in Fig.10,the shapes of different curves remain unchanged,meansing that during the transient unloading process,with an increase in lateral stress,the dynamic strain at the fixed end increases and the lateral stress has no influence on the changing trend of the curve.

Figure 10:Dynamic strain curve at fixed end of rock bar under three-dimensional stress state (σx0=10 MPa,ts=0.5 ms)

According to the physical equation of space problem in elasticity,the influence of lateral pressure (y,z direction pressure) on the axial (x direction) strain is as follows:

According to Eq.(32),the strain generated by direct stress in Fig.10 is calculated from Eq.(31).Compared with Fig.10,the calculated results agree well with the results from theoretical analysis.It can be concluded that the theoretical analysis is suitable for the dynamic unloading problem of rock masses under triaxial stress conditions.First,during the transient unloading process,rock masses will be considered as having unidirectional stress states,and the strain induced by unidirectional stress is calculated from the theoretical analysis;then,the lateral stress will be considered,and Eq.(32) will be applied to calculate the total strain in practical projects.Even though the dynamic adjustment of lateral stress is not considered,the proposed method does make progress for the application of the dynamic unloading problem of the initial geostress,which is dominated by the direction perpendicular to the excavation face.

4 Numerical simulation of dynamic unloading

The well-known software for finite element analysis ANSYS/LS-DYNA is adopted as the preferred numerical simulation method.A rich material model library of ANSYS/LSDYNA has been configured to simulate the nonlinear mechanical properties of rock masses.The implicit-explicit solution mode of ANSYS/LS-DYNA is applied to simulate the process of dynamic unloading to investigate the fracture mechanisms.Among these,the explicit solver is based on the dynamic equation to differentiate the time;the solution of the implicit solver is independent of time and is based on the static equilibrium equation for an iterative solution.The explicit LS-DYNA method is ideal for dealing with transient dynamics problems over short periods,but it is not as efficient as the implicit ANSYS solution method when dealing with static problems.Therefore,in simulating transient offload events after such static loading,these problems can be solved in an implicit-explicit sequence by using these two methods.Moreover,the theoretical analysis should be verified by numerical simulation.In the following section,the implicit calculation mode of ANSYS is employed to first analyze the stress state of the rock bar under initial pressure.The displacement data of model nodes obtained from the implicit calculation are used to initialize the stress state for the explicit calculation.Then,the explicit solver is initiated to calculate transient excavation unloading.In the calculation,the excavation area is excavated mainly by adding a failure unit in LS-DYNA,that is,setting the unit failure time in the *MAT_ADD_EROSION keyword,and directly deleting the unit of the excavation area when the excavation time arrives.

The numerical simulation of the one-time through-excavation of the underground powerhouse floor of Baihetan Hydropower Station shown in Fig.1 is carried out to analyze the dynamic strain on the retained rock mass.In order to improve the calculation efficiency,a scaled scale model is established,as shown in Fig.11.The rock mass in the model is an ideal linear elastic material,and its basic mechanical parameters are shown in Tab.2.The model boundary conditions are full constrained at the bottom and non-reflective boundaries at the sides.

Figure 11:Unloading calculation model for excavation of floor of underground powerhouse of hydropower station

Initial in-situ stress isσ0=17.5 MPa，and the rock mass in the through area is excavated fromt=0to simulate dynamic unloading.The dynamic strain curves of measuring points A and B are shown as follows:

Figure 12:Unloading dynamic strain curve of underground powerhouse floor excavation of hydropower station

As shown in Fig.12,the dynamic strain curve of dynamic unloading of the floor excavation of underground powerhouse of hydropower station has a high similarity with the dynamic strain curve calculated by theoretical one-dimensional theory,which shows that it is feasible to simplify this kind of problem according to the one-dimensional bar unloading model.

5 Dynamic unloading test verification

A test platform for dynamic loading and unloading is constructed based on a theoretical mechanical model.The test platform consists of several systems:the integrated framework,the static loading system,the dynamic loading system,the structural support,and the structure to fix the end of a rock bar and the test system.The layout of these systems is shown as Fig.13.The static loading system is located at one end of the integrated framework.The dynamic loading system is located between the rock bar and the static loading system,which is used to provide dynamic unloading by explosion,lateral sliding,or crushing.The structural support is located in the middle of the integrated framework to control the lateral displacement of the rock bar.The structure for fixing the end is located at the other end of the integrated framework to seal and fix the rock bar with concrete.A dynamic strain indicator (uT3408FRS-DY) is used in the test system to measure the dynamic strain induced by the axial unloading from the preloaded rock bar.

Figure 13:Physical drawing of axial loading and unloading device of rock bar

Figure 14:Schematic diagram of unloading experiment

As shown in Fig.13 and Fig.14,the experimental specimen consists of granite with a cross section of 5 cm×5 cm and a length of 1.5 m.The basic mechanical parameters are shown in Table 1.The rock bar is fixed in the experimental platform,and one end is a 0.3 m length of the rock bar sealed in concrete as the fixed end.The actual length of the rock bar to be measured is 1.2 m,and two strain gauges are arranged along the axis of the rock bar.The #1 strain gauge is located at a distance of 0.65 m from the fixed end of the rock rod.The #2 strain gauge is located at a distance of 1.05 m from the fixed end of the rock.A brittle object (such as a glass bottle) is placed between the static loading system and the rock bar.When the axial preloading stress reaches the compressive strength limit of the brittle object applied to the rock bar by a static loading system,the brittle object is destroyed,and the rock bar boundary abruptly changes,and the rock bar then experiences a sudden unloading effect.When the preload of the rock bar is released,the dynamic strain signals in the rock during the experiment are measured by the strain gauge on the rock.Finally,the time-strain curves of the stress and strain at different positions of the rock are obtained by the dynamic strain signal and can be compared with the theoretical and simulation results.

According to the experimental results,the initial unloading pressure and the unloading duration during the unloading process are obtained.Then,the data from the experiment are defined as the unloading condition for theoretical analysis to calculate the dynamic strain.An experiment is conducted to apply a load to the rock bar until the glass bottle fractures.The experimental data indicate that,when the bottle is crushed,the strain of rock bar reaches 350 με.The initial pressure is 17.5 MPa,the unloading time is 1 ms and the unloading strain rate is approximately 0.35 s-1.The change in the dynamic strain versus time is shown in Fig.15(a) in detail.With increasing time the strain increases more rapidly.In other words,the unloading strain rate increases gradually.The dynamic strains from the theoretical and experimental analyses are compared in Fig.15.Regardless of whether the theoretical or experimental results are considered,the unloading strain at a point 1.05 m from the free end of the bar is very small,as shown in Fig.15(a);as shown in Fig.15(b),the unloading strain at a point 0.65 m from the fixed end is greater.Overall,the entire trend and the attenuation law of strain versus time are verified by the experiment.Moreover,the experimental results match well with the theoretical results at some positions.However,there are some discrepancies between the theoretical and experimental results.The Reasons for these discrepancies are listed as follows:

(1) One end of the rock bar is fixed rigidly in the theoretical analysis,while in the experiment,the end is fixed by sealing 0.3 m of the rock bar in concrete.Therefore,the boundary condition of the fixed end in the experimental case is much more complicated than that of theoretical analysis.

(2) The rock bar is assumed to be made of homogeneous and isotropic elastic material in the theoretical analysis.

(3) The vibration of the integrated framework in the experimental case has an influence on the dynamic strain.

(4) The trigonometric unloading function is applied in the theoretical analysis,whereas the unloading function in the experimental case is only similar to the trigonometric function.

(5) There are inevitable operational errors in the experiment.For example,the initial load is displaced slightly from the bar axis.

There are discrepancies in the rock properties,boundary conditions,and loading functions between the experimental and theoretical analyses.In spite of these inevitable discrepancies,the experimental signal wave and the attenuation of strain agree well with those of the theoretical analysis.The dynamic strain induced by dynamic unloading can be calculated by the wave equation,and the results are satisfactory.

Figure 15:Dynamic strain curves of test results and theoretical results

6 Conclusions

In this paper,a theoretical mechanical model for the unloading of a deep rock mass is established.The rationality of the theoretical model is verified by numerical calculations and experiments,and the following conclusions are obtained:

1) The model clearly reveals the propagation process of an unloading wave in a rock bar and shows that the generation of tensile strain in the rock bar is due to the reflection of the unloading wave by the fixed boundary.

2) Through analyzing the rock strain rate,we find that,when the magnitude of unloading strain rate is greater than 10-1s-1,unloading is a dynamic process,which causes significant tensile strain,resulting in fracturing of the rock mass.

3) Under the condition that the initial stress and the unloading time are the same,the tensile strain generated by the exponential unloading function is the most significant;the trigonometric function is second in importance;and linear strain is the least important.However,when the unloading time increases to a certain level,the strain peak caused by the exponential function may decrease slightly due to superimposition with the initial stress.In the experiment,the unloading curve shape for the initial stress is closer to that of the trigonometric function.

4) By studying the influences of the bar length and the damping coefficient on tensile strain,the fluctuation of the bar in the transient unloading process gradually changes from underdamped to overdamped as the length of rock bar increases or the damping coefficient increases;the dynamic strain loses reciprocity,and no tensile strain is generated.

5) Both the finite element simulation and the experimental results show that the onedimensional wave equation can be used to describe the transient unloading problem of a deep rock mass and that accurate results can be obtained.The theoretical model can effectively solve the dynamic unloading problem of initial ground stress,which occurs dominantly in the direction perpendicular to the excavation face.

In addition,it should be noted that the application scope of one-dimensional unloading model is limited,such as vertical slope,high side wall and excavation of the floor of underground powerhouse of hydropower station.In these problems,excavation will form a plane,and unloading wave can be approximated as plane wave and propagate to the interior of rock mass.At this time,the three-dimensional unloading problem can be simplified into one-dimensional unloading model,and the stress wave equation can be solved directly in Cartesian coordinate system,which can achieve better results.

Acknowledgements:The research work has received funding from the National Natural Science Foundation of China (Grant Nos.51479147,51779193).This work was supported by the Major Program of Technological Innovation of Hubei Province (Grant No.2017ACA102).