Analysis of the damage mechanism of strainbursts by a global-local modeling approach

Jun Wng,Derek B.Apel,*,Artur Dyzko,Andrzej Wlentek,Stnisłw Prusek,Huwei Xu,Chong Wei

a School of Mining and Petroleum Engineering,University of Alberta,Edmonton,Alberta,T6G 2R3,Canada

b Jastrze˛bska Spółka We˛glowa S.A.Group,Aleja Jana Pawła II 4,Jastrze˛bie-Zdrój,44-330,Poland

c Central Mining Institute,Plac Gwarków,Katowice,140-166,Poland

Keywords:Strainburst Numerical modeling Damage mechanism Finite difference method (FDM)Discrete element method (DEM)Underground mining

ABSTRACT Strainburst is the most common type of rockbursts.The research of strainburst damage mechanisms is helpful to improve and optimize the rock support design in the burst-prone ground.In this study,an improved global-local modeling approach was first adopted to study strainburst damage mechanisms.The extracted stresses induced by multiple excavations from a three-dimensional (3D) global model established by fast Lagrangian analysis of continua in 3 dimensions (FLAC3D) are used as boundary conditions for a two-dimensional(2D)local model of a deep roadway built by universal distinct element code (UDEC) to simulate realistic stress loading paths and conduct a detailed analysis of rockburst damage from both micro and macro perspectives.The results suggest that the deformation and damage level of the roadway gradually increase with the growth of surrounding rock stress caused by the superposition of mining-or excavation-induced stresses of the panel and nearby roadways.The significant increase of surrounding rock stresses will result in more accumulated strain energy in two sidewalls,providing a necessary condition for the strainburst occurrence in the dynamic stage.The strainburst damage mechanism for the study site combines three types of damage: rock ejection,rock bulking,and rockfall.During the strainburst,initiation,propagation,and development of tensile cracks play a crucial role in controlling macroscopic failure of surrounding rock masses,although the shear crack always accounts for the main proportion of damage levels.The deformation and damage level of the roadway during a strainburst positively correlate with the increasing peak particle velocities (PPVs).The yielding steel arch might not dissipate kinetic energy and mitigate strainburst damage effectively due to the limited energy absorption capacity.The principles to control and mitigate strainburst damage are proposed in this paper.This study presents a systematic framework to investigate strainburst damage mechanisms using the global-local modeling approach.

1.Introduction

Rockburst is a frequently encountered dynamic rock failure phenomenon in deep underground excavations.It is associated with a rapid release of energy stored in surrounding rock masses of underground openings and is usually characterized by the sudden and violent ejection of rock materials.It can damage equipment and even result in injuries and fatalities (Pechmann et al.,2008;Durrheim,2010; Naji et al.,2018; Yang and Xiang,2018; Pu et al.,2019).To date,rockburst events have been reported in all mining countries (e.g.South Africa,Australia,China,India,Canada,USA,Poland,Russia,and Chile).Some civil engineering projects such as deep tunnels in Switzerland,Norway,Iran,Peru,and China also have suffered rockburst problems (Farrokh and Rostami,2009;Zhang et al.,2012; Dammyr,2016).Unfortunately,as excavation activities continue progressing to greater depths,the frequency and severity of rockburst events increase due to the higher in situ stress and more complex geological environments(Kaiser and Cai,2012;Manouchehrian and Cai,2018; Gao et al.,2019a),which are undisputed evidence that rockburst is a severe and global problem that urges much work to be done to prevent and mitigate rockburst damage.

Ortlepp and Stacey (1994) suggested that it is essential to distinguish the“source mechanism”and the“damage mechanism”of rockbursts,because they are often not the exact mechanism and might have different locations.The source mechanism is the trigger factor that induces rockbursts.In general,rockbursts can be grouped into three types: strainburst,pillar burst,and fault-slip burst based on different source mechanisms (Hedley,1992; Kaiser and Cai,2012; Cai,2013).Strainburst is the most common rockburst type(Zhang et al.,2012;Cai,2013).It happens because of the excavation-induced stress concentration and a “soft” loading environment around fractured rock masses (Kaiser and Cai,2012).The rock can fail locally at excavation boundaries in an unstable and violent manner.The damage mechanism refers to failure modes(e.g.rock ejection,rock bulking,rockfall,rock buckling,and shear displacement) induced by rockbursts.The research of rockburst damage can provide insight into the understanding of initiation,development,extent,and types of failure within surrounding rock masses during rockbursts.The rock failure type is an essential criterion to select rational support elements (e.g.cable bolts or rockbolts),and the damage severity can affect the capacity,extent,and intensity of the support system (e.g.strength and length of rockbolts).Thus,the understanding of rockburst damage mechanisms is helpful to improve and optimize the rock support design in the burst-prone ground to control and mitigate rockburst damage.The focus of this study is on strainburst damage mechanisms.

At present,significant efforts have been made to study strainburst damage mechanisms by experimental,field surveying,and numerical methods.Various engineering disciplines have employed numerical modeling,ranging from bioengineering(Apel et al.,2021)to mining engineering.Compared with field surveying and experimental tests,numerical methods have the advantages of low cost,safety,time-saving,and flexibility.Numerical modeling can also obtain more information(Salamon,1993).In the past few decades,many numerical methods and programs have been developed and used to deal with rock mechanics problems (Wang et al.,2021a).These methods have become a valuable and frequently used way to study the damage mechanism of strainbursts.Zhang et al.(2011,2013) proposed a failure approaching index (FAI) and then programed it in FLAC3D (fast Lagrangian analysis of continua in 3 dimensions) to estimate the potential rockburst damage zones in deep powerhouses and tunnels.Hosseini (2016) used the CA3 program to study the strainburst damage where discrete particles modeled the rock,and finite elements simulated the frame structure.Feng et al.(2017,2019)employed ELFEN to analyze the effects of structural planes and stress regimes on the process of rockbursts occurring in a circular opening.Manouchehrian and Cai (2018) built a heterogeneous model in ABAQUS via Python programming and then simulated rockburst damage near faults in deep excavations.Gao et al.(2019a,b) proposed a novel bonded block method in UDEC (universal distinct element code)to model strainbursts and studied rockburst damage by monitoring failed cracks.Vazaios et al.(2019)employed IRAZU (finite element method (FEM)-discrete element method(DEM) program) to study the effect of a different number of preexisting joints (discrete fracture networks (DFNs)) on strainburst phenomena in a hard rock excavation under high magnitude stresses.Guo et al.(2019)adopted FLAC3D numerical simulation to investigate the failure laws of the roadway surrounding rocks during rockbursts considering different ratios between principal stresses.The research work of those scholars has made many significant achievements and provided good references for understanding the damage mechanisms of strainbursts under different conditions.

However,most current research belongs to the parametric study(e.g.different factors related to excavation geometry,stress scenario,discontinuity,and material property) without considering the gradual stress concentration or energy accumulation resulting from nearby mining or excavation activities,indicating that the influences of the realistic stress loading path on strainburst occurrence were ignored.Furthermore,most studies’ analyses of rock mass fracturing or damage are qualitative,although it is meaningful to use a quantity to describe different damage levels meticulously.Therefore,it is worth developing or employing a new modeling approach to capture the realistic stress loading path and investigate rock mass fracturing or damage quantitively during strainbursts.For example,Edelbro (2008) developed the globallocal modeling approach as an efficient and effective method to study the failure and deformation of roadways.In this approach,the various stresses from excavating multiple openings can be extracted accurately from an overall two-dimensional (2D) or three-dimensional (3D) model and then used as boundary conditions for 2D local models of roadways.

Besides,the computational efficiency is significantly enhanced because of the usage of elastic properties in the global model,and the detailed analysis of 2D local models can be accomplished with a high mesh resolution and elastoplastic relationships.Hence,this approach might have great potential to investigate the damage mechanism of rockbursts,especially the strainbursts,since they usually occur in tunnels or roadways.However,to the authors’knowledge,no studies in the literature have been reported to use this approach to analyze strainburst damage.

The presented research used back analysis to simulate the state of mining-or excavation-induced stresses before the rockburst event that happened at the Zofiowka Mine in Poland on May 5,2018 (Adam,2018).The rockburst event had a seismic energy of 2 × 109J,causing tremendous damage to roadways and even fatalities.The location where the rockburst occurred had a complex geological structure with many faults.The complex geological and mining environments resulted in high concentrated stresses in roadway surrounding rock masses,which provides a good case example for studying strainburst damage mechanisms.

In this study,an improved global-local modeling approach is first employed to study strainburst damage mechanisms.First,a 3D finite difference numerical model is built via FLAC3D software.The elastoplastic constitutive relationship of rock masses is used in 3D modeling to capture the actual and accurate mining-or excavation-induced stress distribution.Then,the extracted stresses induced by multiple excavations from the 3D global model are used as boundary conditions for a 2D local model of a deep roadway built by a DEM code UDEC.The damage mechanism of strainbursts is investigated from both micro and macro perspectives.Finally,the effects of the current support system are initially evaluated,and the control principles of rockburst supporting are proposed according to the analysis of strainburst damage mechanisms.This study presents a systematic framework to investigate strainburst damage mechanisms using the improved global-local modeling approach.

2.Engineering overview

2.1.Geology and geotechnical overview

The Zofiówka Mine is located around 307 km southwest of Warsaw in Poland.The mine produces around 3.7 million tonnes of coal per year.The analyzed area is in Section H,southeast of the mine,and it is surrounded by large faults,e.g.Central fault(throw:10-11 m),Eastern fault (throw: 15-20 m),and Jastrze˛bie fault(throw:15-60 m).There are two main coal seams in Section H:coal seams 409/3 and 409/4.The buried depth of coal seams is 815-990 m and the dip angle is 5°-10°.The coal seams 409/3 and 409/4 have the thicknesses of 1.8-2.7 m and 4-5.3 m,respectively.The roof and floor are shales and sandstones.Fig.1 shows the rock strata histogram (G-42 borehole).The mine uses the retreating longwall mining method to extract coal.Stoping of Panels H-4 and H-6 in coal seam 409/3 had been finished during 2016-2017.Headgate and tailgate were advancing to prepare the Panel H-4 in coal seam 409/4 in 2018(Fig.2).The roadways were supported by the yielding steel arch (ŁP12/4/V32) and welded wire mesh.

Fig.1.Geological column in the study site.

Fig.2.Location of Section H and the rockburst event.

A violent rockburst event happened in the intersection area of H-4 tailgate and H-10 main drift(Fig.2).The rockburst event has a seismic energy of 2 × 109J and its Richter magnitude is 3(Bieniawski,1987).The seismic energy is calculated using squared and integrated velocity records of radiated body waves at a short distance from a point source,and assuming a focal sphere at 500 m from the source (Mutke et al.,2016).Significant rockfall and floor heaving are observed in the H-10 main drift with a total length of around 300 m (Wang et al.,2021b).Unfortunately,this accident resulted in four injuries and five fatalities.

2.2.Rock mass properties

The uniaxial compressive strength (UCS) test and Brazilian tension test were conducted to determine the mechanical parameters of intact rocks following ISRM standards (Fairhurst and Hudson,1999).Then,the rock mass properties were obtained by scaling the mechanical parameters of intact rocks using the generalized Hoek-Brown criterion (Marinos and Hoek,2000) and geological strength index (GSI) system (Fig.1).The values of GSI were determined from the quantitative GSI chart (Marinos and Hoek,2000) according to the inspection of rock mass structures and the surface condition of discontinuities.The rock mass properties used in numerical modeling are listed in Table 1,which can also be referred to literature (Małkowski et al.,2017; Szott et al.,2018; Małkowski and OstrowskI,2019).The UCS and deformation modulus (Em) of rock masses were estimated from the following equations (Hoek et al.,2002; Hoek and Diederichs,2006):

where mb,s,and a are the constants for rock masses;σciis the UCS of intact rocks;Eiis the Young’s modulus of intact rocks;σcmis the UCS of rock masses;Emstands for the deformation modulus of rock masses;and D is a factor representing the disturbance level of rock masses caused by different tunneling methods.For the study site,D is assumed to be zero considering that the mechanical tunneling results in minimal disturbance to confined rock masses(Hoek et al.,2002).The calculated results of UCS and Emof rock masses are also summarized in Table 1.

Table 1 Physical and mechanical parameters of rock masses.

3.Numerical modeling

3.1.Global-local modeling approach

The global-local modeling approach was initially proposed to study the brittle failure of footwall drifts in an underground hard rock mine (Edelbro,2008; Edelbro et al.,2012).In this approach,stresses induced by mining multiple stopes were extracted from an overall 3D model and then were input into 2D local models of drifts as boundary conditions.Due to elastic properties’ employment in the global 2D model,the computational efficiency is significantly enhanced,and the detailed analysis of 2D local models can be accomplished with a high mesh resolution and elastoplastic relationships.Edelbro et al.(2012)reported that the simulated failure of drifts agreed well with the in situ observations,indicating that the used approach is an efficient and robust way to analyze the damage and stability of tunnels.After that,the approach was used to analyze the stress and deformation of gateways in longwall panels.Basarir et al.(2019) established a 3D global model by FLAC3D to extract the mining-induced stresses implemented into 2D local models built by RS2 (Rocscience,2015) for the stability analysis of gateways.The modeled deformation of gateways in 2D local models matched well with the monitored roadway convergences,which further verified the effectiveness of this modeling approach.

The global-local modeling approach might have great potential to investigate the damage mechanism of rockbursts,especially the strainbursts,which usually occur in tunnels or roadways.First,the 3D global model can capture the realistic stress loading path before a strainburst occurs.Second,rock mass fracturing or damage during strainbursts can be investigated quantitively and in detail in 2D local models with the employment of elastoplastic relationships and high mesh resolution.Hence,this approach was employed in this study.

However,there are two significant improvements in the used modeling approach.First,the elastoplastic constitutive relationships of rock masses rather than elastic models are used in the 3D modeling to capture the redistribution of mining-or excavationinduced stress more realistically and accurately.It seems that we go back again,because the elastoplastic constitutive relationships are employed,which might increase the calculation cost significantly.Nevertheless,over the last two decades,with rapid developments in computer equipment,the calculation speed of numerical programs (e.g.FEM and FDM codes) has been significantly increased.For instance,it takes around 0.15 s for the FLAC3D(Version 7.0,Itasca,2019) to process 624,000 zones in 100 timesteps on a computer with an Intel i9-9820 CPU at 3.30 GHz (10 cores) due to better multi-threading and improved algorithms(Hazzard,2020).Thus,the usage of elastoplastic properties for global 3D modeling is no longer an inefficient and luxury thing at present.

Second,a 2D DEM code UDEC (Version 7.0,Itasca,2020) was employed in this research,while the 2D FEM code RS2 was used in the previous studies.The rockburst failure is a process that rock masses suffer a continuous small deformation first and then a large discontinuous deformation in a very short time.However,continuum methods,such as boundary element mthod (BEM),FEM,and finite difference method (FDM),are unbale to model the discontinuous deformation behavior of rock masses.Hence,researchers preferred to adopt DEM and DEM-related hybrid methods to model the discontinuous deformation of rock masses during rockbursts(Grisi et al.,2016; Feng et al.,2017,2019; Li et al.,2018; Gao et al.,2019a,b; Vazaios et al.,2019).The rock mass consists of distinct blocks and contacts in the block DEM.The blocks will slide or separate with each other when the stress exceeds contact’s strengths.Besides,the blocks can move,rotate,and deform according to Newton’s second law and related constitutive relationships (Jiang,2017).Thus,the DEM is suitable for modeling the discontinuous deformation of rock masses to reproduce rockburst damage processes(Wang et al.,2021a).It might be argued that only a single 3D DEM code can solve this problem.However,the main concern is that the calculation speed of current 3D DEM codes(e.g.3DEC)is much slower than those of 3D FEM and FDM,and 2D DEM codes(Yang et al.,2018).Thus,the combination of 3D FDM and 2D DEM (global-continuous and local-discontinuous) can be accepted as a more efficient and rational approach to the current investigation of strainburst damage mechanisms.

3.2.Global model

3.2.1.Model setup

A 3D finite difference numerical model was established using FLAC3D software,as shown in Fig.3.Instead of all coal basins,only the most critical parts,including H-4 headgate,H-4 tailgate,H-10 main drift,and a part of Panel H-6 were built in the model,as this study aims to investigate the strainburst damage mechanism for a roadway.It could be argued that the major faults,e.g.Central fault,Eastern fault,and Jastrze˛bie fault,may contribute to the redistribution of stress surrounding the main drift.Tsesarsky et al.(2013)stated that the discontinuities could be assumed continuous over areas much greater than any excavation.In the study site,the scale of excavations (tunneling faces of the gateway) is much smaller than that of main faults,with a length of a few hundred meters to more than 1 km.Therefore,it is assumed that faults are not included in the 3D global model,which significantly reduces the size and zone number of the model and saves the computational cost.The model’s top boundary was free,and vertical stress of 21.87 MPa was applied to it to model the overburden weight.This value is determined by multiplying the unit weight of overburden(assuming 0.027 MN/m3) and the distance from model’s top boundary to the ground surface (810 m).The horizontal displacement on model’s side boundaries is not allowed.Both horizontal and vertical displacements are constrained on the bottom boundary.The initial stresses were set based on in situ stress measurements from the study site.To ensure the efficiency and accuracy of computation,the roadways and surrounding rock masses were discretized into small zones with a side length of 1 m and the remaining domain was set with larger sizes.The model has 1,066,921 zones in total,which is determined according to a mesh convergence study (Sepehri et al.,2020).

Fig.3.The layout of the numerical model.

3.2.2.Constitutive model and material properties

(1) Constitutive model of rock masses

The Hoek-Brown failure criterion is suitable for modeling the mechanical behavior of jointed rock masses (Hoek and Brown,1997),and it has been widely used in numerical modeling for the stress and deformation analysis in many Polish coal mines (e.g.Małkowski et al.,2017; Małkowski and OstrowskI,2019).Besides,the brittle failure behavior of rock masses (e.g.spalling and rock ejection)during strainbursts was studied in UDEC modeling rather than FLAC3D in this study.Therefore,the mechanical behavior of rock masses in FLAC3D is governed by the Hoek-Brown failure criterion,which is expressed as follows (Hoek and Brown,1997):

where σ1and σ3are the major and minor principal stresses,respectively.It should be noted that the pore pressure of rocks was not considered in this study.

(2) Constitutive model and material properties of the gob

In coal mining,the overburden stress carried by the mined coal is then supported by surrounding rock masses and gob materials,resulting in the redistribution of mining-induced stresses.Hence,it is essential to choose a suitable constitutive model to simulate the mechanical behavior of gob materials and reproduce the realistic state of mining-induced stresses.

The null and elastic models were often used to simulate gobs(Kose and Cebi,1988; Cheng et al.,2010; Jiang et al.,2012; Basarir et al.,2015,2019).However,the gob materials are not even simulated when using the null model,and the elastic model cannot capture the strain-hardening behavior of compacted gob materials.Some research has recently verified that the double yield (DY)model can simulate the strain-stiffening behavior of gob materials(Yavuz,2004; Qiu et al.,2019; Wang et al.,2020).This study also adopted the DY model to simulate the gradual compaction of caved rock materials in the gob.

The DY model requires the cap pressure and material parameters of the gob.The Salamon model can calculate the cap pressure(Salamon,1990).The gob material parameters can be obtained according to a detailed calibration procedure that can be referred to Wang et al.(2021b).Following this procedure,the cap pressure and gob material parameters are finally determined.Fig.4 shows that simulated gob stress agrees well with the cap pressure calculated by the Salamon model.The used gob material parameters are listed in Table 2.

Table 2 Gob material parameters in numerical modeling.

3.2.3.Validation of the global model

(1) Validation of gob modeling

As shown in Fig.5,a stress monitoring line is arranged in the gob and coal seam 409/3 to plot the vertical stress distribution.The vertical stress in the gob gradually rises to 92.59% of the in situ vertical stress at 196 m(0.22 times the average depth of coal seam 409/3)from the edge to the center.The vertical stress in the gob will be recovered to the in situ stress level at a distance of 0.2-0.3 times the mining depth according to field measurements (Wilson and Carr,1982; Campoli et al.,1993).Hence,the simulated results are in good agreement with field measurements,verifying the effectiveness of the DY model and gob material parameters.

Fig.4.Stress-strain curves of the numerical model and Salamon model.

(2) Validation of the full model

The full model was verified in this research by comparing the simulated deformation and failure depth of surrounding rock masses of the headgate with field measurement data.As shown in Fig.6a and b,the simulated convergences of roof-to-floor(350 mm)and rib-to-rib (490 mm) are within the measured roadway deformation range(Wang et al.,2021b).Additionally,Fig.6c and d shows that the simulated failure depths of the roof(7.02-7.93 m)and floor(2.64 m)match well with the monitored rock damage depths(8.2-8.5 m and 2.7 m).Thus,the used material properties,constitutive relationships,and boundary conditions in the global model and the model geometry are validated.

Fig.5.Vertical stress distribution in the gob and coal seam 409/3.

3.2.4.Mining-and excavation-induced stresses

The global model has calculation steps as follows: (1) initial stress equilibrium; (2) mine out Panel H-6 and conduct gob modeling; (3) tunnel H-4 headgate and a part of H-4 tailgate; and(4)tunnel H-4 tailgate at 5 m/step with 11 excavation steps in total when it is near to H-10 main drift.Fig.3 shows the excavation directions of the panel and roadways.The normal and shear stresses at the monitoring point (see Fig.3) of the main drift in the global model were recorded at different excavation steps.The extracted stresses are shown in Fig.7,whose orientations have been transferred in the 2D local model for further usage.The excavation of Panel H-6 and a significant part of gateroads significantly influences the redistribution of the stresses in the main drift.After that,the stresses increase slowly with the driving of the tailgate.This is because the distance from the tunneling face to the monitoring point is great,and the influence of the excavation-induced effect on stress concentration is smaller than that of the mininginduced effect.Finally,there is a sudden surge of stresses,and the direction of shear stress is also deflected when the last excavation step (Ex11) is finished.This is due to the tunneling of the tailgate from the main drift.It should be noted that the plane strain assumption for the main drift in the 2D local model is no longer valid during the Ex11 excavation step.Therefore,the extracted stresses in this step were not considered for subsequent simulation in the local model.

3.3.Local model

3.3.1.Model setup

The 2D local model was extracted from the global model’s Y-Z plane along the tailgate axis,as shown in Fig.3.A widely used 2D DEM code UDEC was adopted to construct the local model for conducting a detailed analysis of the damage mechanism of strainbursts.Fig.8 shows the local model,which is built according to the lithology(see Fig.1)and the designed size of the main drift.A Trigon approach (Gao and Stead,2014) was applied to generate blocks,as this approach can be capable of reproducing the natural fracturing processes and dynamic mechanical behavior of rock masses without adopting complicated constitutive models(Stavrou et al.,2019).In the Trigon approach,a rock mass is composed of triangular blocks and contacts (Gao et al.,2015).The fracturing process can be exhibited either by the sliding or opening of contacts.In the 2D local model,the average blocks’edge lengths in two coal seams and nearby clay shale are set to 0.3 m.The block size with a range of 0.2-0.5 m was sufficiently fine to model the failure behavior of surrounding rock masses for a 2D model (Gao,2013;Gao and Stead,2014; Zang et al.,2020).The average blocks’ edge lengths in the upper clay shale,sandy shale,and fine-grained sandstone were set to 0.5 m,0.5 m,and 1 m,respectively.The average blocks’ edge lengths in the floor are 0.3 m and 1 m.A graded increasing edge length of blocks can avoid the resulting loss of simulation accuracy and enhance the calculation’s reliability.The horizontal displacement on the model’s side boundaries is not allowed.Both horizontal and vertical displacements are constrained on the bottom boundary.The normal and shear stresses shown in Fig.7 were applied to model boundaries.

3.3.2.Constitutive model and material parameters

The elastic constitutive model was chosen for blocks composed of finite difference zones.The Coulomb slip model was used for contacts.A spring-rider simulates the behavior of contacts,and the model deformation occurs when the contact stress is smaller than the contact strength,which is governed by the elastic modulus of blocks and contact stiffness;contact failure occurs when the stress exceeds its shear or tensile strength,and then blocks will slide or separate with each other (Chen et al.,2016).The constitutive behavior of contacts is shown in Fig.9.

Fig.6.(a-c) Simulated vertical and horizontal displacement contours and failure zone of the headgate,respectively; and (d) Monitored damage depth.

Fig.8.2D local model: (a) Static stage; and (b) Dynamic stage.

I n the Trigon approach,the macro properties of rock masses(e.g.UCS and Em)depend on the micro parameters(e.g.contact strength and stiffness)of blocks and contacts(Gao,2013;Gao et al.,2015).In this research,these micro parameters were calibrated with rock mass properties (see Table 1) via simulated UCS tests (Gao et al.,2015).To eliminate block size’s effects on simulation accuracy,the UCS test model had a large scale(4 m×8 m)(Yang et al.,2017)and identical block size with the main drift model.However,there is a problem that different block sizes were employed for the rock strata with the same lithology(e.g.block size of 0.3 m and 0.5 m for clay shale,and 0.3 m,0.5 m,and 1 m for sandy shale),which means that different material parameters might be used even for the same lithology.A sensitivity study was conducted and has shown that the effect of block size on simulated rock mass properties can be negligible.The loading in tests was simulated by applying a velocity of 0.1 m/s to the top platen’s surface (Fig.10).This loading rate is slow enough to avoid the dynamic responses of models after conducting a sensitivity study.The initial micro parameters were first assumed according to rock mass parameters.Then,the modeling of UCS tests was conducted iteratively with the adjustment of micro parameters until simulated results were identical with targeted material properties.Gao(2013)and Tan and Konietzky(2014)gave a detailed calibration process of micro parameters.The simulated failure modes and stress-strain curves of rock mass samples are illustrated in Fig.10.The main failure modes of rock mass samples are tensile (axial splitting) and tensile-shear failure,which are consistent with typical failure modes of rock masses under no or low confining pressures (Diederichs,2007).Table 3 shows the calibrated micro parameters of blocks and contacts.The targeted and simulated Emand UCS errors are all less than 3% (Table 4),suggesting that the targeted values match well with the calibrated rock mass parameters.Thus,the calibrated micro parameters in Table 3 could be used to further numerical analysis of the damage mechanism of strainbursts.

Fig.9.Constitutive behavior of contacts(after Yang et al.,2017).cj,φj,and σjt are the cohesion,internal friction angle,and tensile strength,respectively;Kn and Ks are the normal and shear stiffness,respectively;Δσn and Δun are the effective normal stress increment and normal displacement increment,respectively;σn and τs are the normal and shear stresses,respectively; and Δues is the elastic shear displacement increment.

Fig.10.Simulated failure modes and stress-strain curves of rock mass samples.

3.3.3.Simulation scheme

The modeling of the damage mechanism of strainbursts was performed using three stages:

(1) Stage I (static stage): The in situ stress field was applied to the model,and the geostatic equilibrium was achieved.Then,the main drift was excavated by deleting the blocks.Finally,adequate calculation steps were run to ensure gradual and slow release of surrounding rock stresses(Gao et al.,2015).

(2) Stage II (static stage): The varied normal and shear stresses extracted from the 3D global model (Fig.7) were applied to local model boundaries.Adequate calculation steps were run to stabilize the model.

(3) Stage III (dynamic stage): The dynamic mode was activated.Studies have shown that the sudden release of elastic energy stored in rock masses is often caused by dynamic loads,resulting in severe rockburst accidents (Ortlepp and Stacey,1994; Lu et al.,2015; Masny et al.,2017; Kong et al.,2019;Ji et al.,2021).Zhu et al.(2010)reported that the disturbance of the external dynamic load is one of the key factors to trigger rockbursts.Kaiser and Cai (2012) proposed that strainbursts can be mining-induced due to static stress changes caused by nearby mining,or they can be dynamically induced due to a dynamic stress increase caused by a remote seismic event.Diederichs (2018) stated that many rockburst damages in mining environments were based on the primary mechanism of a remote seismic event triggered by large-scale mining activities,while the primary rockburst source is the roadway surrounding rock mass in geotechnical engineering projects.The source mechanism of rockbursts associated with the seismic waves and high static stress is currently the most common in Polish underground coal mines (Mutke et al.,2015).Therefore,it is necessary and meaningful to study the effects of dynamic loads on strainburst damage mechanisms.The dynamic load is exerted on mining works in a fashion of vibration wave or stress wave resulting from rock fracturing,fault-slip,blasting,mechanical vibration,etc.A seismic wave caused by fault-slip was assumed to be the dynamic load source because post-event observations show the activation of the Jastrze˛bie and Eastern faults.Mutke et al.(2015) and Kong et al.(2019)reported that the rockburst potential of roadways is positively correlated to the peak particle velocity (PPV) of vibration waves and buried depths,as shown in Table 5.The statistical analysis of rockbursts in the Upper Silesian Coal Basin(USCB)suggests that PPVs are mainly within 0.05-1 m/s(Fig.11a and b).Additionally,rockbursts are usually related to seismic waves characterized by low frequencies from 10 Hz to 30 Hz (Fig.11c).Therefore,seven PPVs of 0.2 m/s,0.4 m/s,0.6 m/s,0.8 m/s,1 m/s,1.2 m/s and 1.4 m/s are adopted to simulate different dynamic loads (σn=2(ρCp)V)(Itasca,2020).

Table 3 Calibrated micro parameters of blocks and contacts.

Table 4 Comparison between the targeted and simulated rock mass parameters.

Table 5 Evaluation of the influence of seismicity on rockburst potential(Mutke et al.,2015;Kong et al.,2019).

Fig.11.(a)One hundred and twenty seismic events associated with rockbursts in USCB(1988-2006);(b)Thirty two seismic events associated with rockbursts in Polish coal mines(2003-2012);and(c)The relationship between frequencies and seismic energies of rockbursts(seismic events which caused significant failure in roadways)in USCB during 1998-2006 (R is the distance from damaged working to seismic sources,and Es is the seismic energy) (Mutke,2008; Mutke et al.,2016).

The frequency is assumed to be 20 Hz,and the time is 120 ms(a vibration period plus a quiet time of 70 ms).The seismic waveform is simplified to be a sine wave since any complex stress wave can be obtained by the Fourier transform of simple sine waves(Liu,2017).The waveforms with different PPVs are shown in Fig.12.The dynamic calculation used the viscous boundary to avoid propagating waves’ reflection and allow the necessary energy radiation.A recommended Rayleigh damping of 0.5%was used(Itasca,2020).This value is suitable for dynamic analyses that involve large block deformation or great joint displacement.Then,a series of seismic waves was applied to the model’s right boundary to investigate the dynamic responses of the roadway.Ideally,the dynamic mechanical properties of rock masses and joints and related constitutive relationships should be employed (Gao et al.,2019b).Considering that those data are not available,and the research’s objective is to show the effectiveness and advance of the global-local modeling approach in studying strainburst damage mechanisms,the properties of rock masses and joints and constitutive relationships were therefore not changed.The simulation procedures are shown in Fig.13.

Fig.12.Waveforms of seismic waves with different PPVs.

Fig.13.Flowchart of the simulation process.

4.Results and discussion

4.1.Static stage

4.1.1.Displacement analysis

The deformation of surrounding rock masses of the main drift during different excavation steps is shown in Figs.14 and 15.The severe roadway deformation is observed when the main drift is not supported.The primary deformation forms of the main drift are roof subsidence and the shrinkage of two sidewalls.After excavating the main drift,the roof subsidence and deformation of two sidewalls are 510.5 mm,361 mm,and 347.3 mm,respectively.The minimum deformation(102.9 mm)appears in the floor due to the relatively high strength of the floor stratum.The extraction of Panel H-6 creates a superposition of the mining-induced stress and initial surrounding rock stress induced by excavation.Thus,the main drift suffers severer deformation during this stage.The displacement of the roof,two sidewalls,and floor are 610 mm,571.2 mm,350.7 mm,and 104.6 mm,and the growth rates are 19.49%,58.22%,1.37%,and 1.65%,respectively.A considerable increment of roof subsidence and deformation of the left sidewall is observed.The deformation difference between two sidewalls is great (220.5 mm),indicating an asymmetric convergence phenomenon.When a part of the gateways was excavated,the deformation of the roof,two sidewalls,and floor were 674.7 mm,709.9 mm,359.7 mm,and 104.6 mm,and the growth rates were 10.61%,24.28%,2.57%,and 0%,indicating a moderate increase of deformations.This is because the main drift is very far from the tailgate’s tunneling face and the excavationinduced effects are much smaller than mining.After the tunneling of the tailgate was finished,the final deformation of the roof,two sidewalls,and floor are 683.9 mm,709.9 mm,382.2 mm,and 104.6 mm,respectively.

Fig.15 shows that the roadway deformation,especially the displacements in roof and sidewalls,gradually increases with different excavation steps.This is due to the growth of surrounding rock stress resulting from the superposition of mining-or excavation-induced stress of the panel and roadways.Extraction of Panel H-6 has the most significant influence on the roadway deformation,suggesting that the deformation values are related to the magnitude of mining-or excavation-induced stress.Besides,the deformation becomes more uneven (see displacement vector maps in Fig.14).The roof and left sidewall suffer more deformation than the right sidewall and floor.The nearby excavations result in the least influence on the floor deformation.The surrounding rock masses of the main drift will lose their stability if no control measures are adopted in time.

Fig.14.Distribution of displacement(in m)and velocity(in m/s)of the main drift during four excavation steps:(i)Excavation of the main drift;(ii)Excavation of the Panel H-6;(iii)Excavation of a part of roadways; and (iv) The 10th excavation step (Ex10) of the tailgate.

Fig.15.Comparison of the deformation at different positions of the main drift during four excavation steps.

Fig.14c shows that some rock blocks at the excavation periphery suffer a large displacement (up to 1.2 m).Thus,it is interesting to discriminate whether a strainburst occurs or not.Fig.14d shows the velocity distribution of rock blocks.The velocity of all rock blocks is 0-3 m/s.The physical equation V2=2gu in Newton’s laws can analyze the velocity inversely based on the gravitational acceleration g and displacement u.If g and u are taken as 9.81 m/s2and 0-1.2 m,respectively,then the velocity will be within 0-4.85 m/s,identical to the velocity of rock blocks.Therefore,the movement style of rock blocks is the natural fall or caving from failed rock masses.This indicates that a strainburst does not occur during the static stage.

4.1.2.Fracture evolution and failure process

Studies have shown that many rock engineering accidents,including rockbursts,are due to the weakening of rock mass strengths resulting from the initiation and development of internal fractures (Esterhuizen et al.,2011; Chen et al.,2016; Gao et al.,2019a,b; Wu et al.,2019a; Zang et al.,2020).Therefore,the fracture evolution and failure process were analyzed to study the damage mechanism of roadway surrounding rock masses during static and dynamic stages.A function was developed using the FISH language in UDEC to record the length and number of failed contacts (including tensile and shear failure) in the local model.A damage variable was then defined in the self-developed FISH function according to the ratio of failed contacts’length to the total contact length in the local model (Gao et al.,2015):

where Dc,Dt,and Dsare the total,tensile,and shear damage levels,respectively; and Lc,Lt,and Lsare the total contact length and the lengths of failed contacts in tensile and shear failure,respectively.

Fig.16 illustrates the crack development in surrounding rock masses during four excavation stages.Fig.17 illustrates the variation of damage levels of the main drift.The comparison of different damage levels is shown in Fig.18.As shown in Fig.16a,tensile cracks are mainly observed in the opening periphery due to radial stress release and the concentration of tangential stress after excavation(Yu et al.,2015).The depth of tensile cracks is 0.3-1.3 m.The shear cracks are well distributed and extend deeper into rock masses,indicating that rock blocks’ slipping gradually develops from the roadway surface to depth.The depth of shear cracks is 0.8-1.8 m.Both tensile and shear damages undergo a rapid increase.Fig.18 shows that the shear damage(1.7%)is the dominant failure mode,accounting for 80.95%of the total damage,while the tensile damage (0.4%) accounts for 19.05%.The ratio of shear to tensile damage is 4.25.After the extraction of Panel H-6,the extent of shear and tensile cracks increased significantly while the damage depth grew slightly.The shear and tensile damage levels are 1.96%and 0.55%,with growth rates of 15.29%and 37.5%,respectively.The ratio of shear to tensile damage is 3.56,suggesting that the tensile failure gradually plays a key role in the stability of surrounding rock masses.The depth of shear cracks is 1-1.9 m,increased by 0.1-0.2 m; and the depth of tensile cracks is 0.4-1.5 m,increased by 0.1-0.2 m.The extent and depth of tensile and shear cracks grow little in the next stage.The evolution of shear and tensile damage gradually becomes steady.The ratio of shear (2%) to tensile (0.6%)damage is decreased to 3.33.When the tunneling of the tailgate is finished,the extent of shear and tensile cracks increases,and the failure depth reaches the maximum values (1.4-2.8 m for shear cracks and 0.7-2 m for tensile cracks).The ratio of shear(2.12%)to tensile(0.62%)damage is 3.41,indicating a slightly decreased ratio of tensile damage.

Fig.16.Evolution of cracks and macroscopic failure patterns of the main drift during four excavation steps: (i) Excavation of the main drift; (ii) Excavation of the Panel H-6; (iii)Excavation of a part of roadways;and(iv)The 10th excavation step of the tailgate.slip_n-Shear failure now;slip_p-Shear failure in the past;tens_n-Tensile failure now;and tens_p-Tensile failure in the past.

Fig.17.Simulated damage evolution of the main drift during four excavation steps.

Fig.18.Comparison of the damage levels of the main drift during four excavation steps.

As shown in Fig.16d,macroscopic fractures are observed around the periphery of the main drift after the initiation,growth,and expansion of microcracks.The main failure mechanisms are the collapse of the roof and spalling of sidewalls.When the extraction of Panel H-6 influenced the main drift,more rock blocks were detached from failed rock masses,producing a huge fractured zone.During the last two stages,macro fractures continue to grow and expand due to the further development of microcracks in depth.However,the failure mechanism changes from the collapse of the roof and spalling of sidewalls to the continuous fall of detached rock blocks.Therefore,the rock and coal masses beyond the fractured zone are still intact,having the capability of continuous stress concentration and energy accumulation for the subsequent occurrence of strainbursts.

4.1.3.Stress analysis

The concentrated stress significantly influences rockburst occurrence and damage in underground openings (Li et al.,2018).Thus,the maximum principal stress surrounding the main drift induced by excavation and mining was analyzed in this study.Fig.19 shows the contour of maximum principal stresses during different excavation steps.After excavation,there is an apparent stress relaxation zone(red and orange areas in Fig.19a)around the main drift.This is because surrounding rock masses fail due to the release of radial stress and the concentration of tangential stress.The stress concentration (green areas) mainly occurs around the stress relaxation zone,and the significant stress concentration(blue areas)is observed in two bottom corners and the roof of the main drift.When the Panel H-6 was extracted,the range of the stress relaxation zone was enlarged as more rock blocks were detached from failed rock masses owing to the influence of mininginduced stress.Compared with the previous stage,the range of stress concentration zone(green areas)is also enlarged because of mining-and excavation-induced stress superposition.After excavating a part of roadways,stress relaxation and concentration zones are further expanded.When the tunneling of the tailgate was finished,the distribution of maximum principal stresses was similar to that in the former stage,and the range of the stress relaxation zone increased a little due to minor excavation-induced effects.

To further study the evolution law of mining-and excavationinduced stress concentration around the main drift,a monitoring line was arranged on the main drift’s right sidewall to record maximum principal stresses’ variation (Fig.8).The simulation results are illustrated in Fig.20.Stress relaxation and increased areas(Fig.20a)are produced due to rock masses’failure and excavationand mining-induced stress superposition,respectively.The average values of maximum principal stresses in the stress increased area during four stages are 30.48 MPa,30.64 MPa,31.68 MPa,and 31.88 MPa(Fig.20b),respectively,showing the gradual increase of surrounding rock stresses in deep zones.Additionally,the peak stresses during the four stages are 37.46 MPa,37.95 MPa,39.22 MPa,and 39.56 MPa,respectively,and the stress concentration coefficients (in situ vertical stress over peak stress) are 1.54,1.56,1.61,and 1.63,respectively.The growth rates are 1.3%,4.55%,and 5.84%,respectively,compared with that in the first stage,suggesting an overall increasing peak stress trend (Fig.20b).The significant increase of surrounding rock stresses will result in more accumulated strain energy in two sidewalls,providing a necessary condition for the strainburst occurrence in the dynamic stage.

4.2.Dynamic stage

4.2.1.Occurrence process of a strainburst

The configuration with a PPV of 1.4 m/s was used as an example to demonstrate the occurrence procedure of a strainburst caused by a seismic wave.The simulated velocity vectors of rock blocks during the strainburst are shown in Fig.21.When the time was 20 ms,the rock masses in the main drift’s right sidewall were initially influenced by the incident seismic wave,while the main drift has not been affected by dynamic stresses.When the time was 40 ms,the surrounding rock masses in the right sidewall and roof were influenced by dynamic stresses and ejected suddenly.The maximum velocity of ejected rock blocks was around 20-28 m/s.When the time was 60 ms,the seismic wave affected the rock masses in all directions,and some rock blocks were ejected from the floor with a velocity of 20-30 m/s.From 60 ms to 80 ms,more rock blocks were ejected from the roof,and a few were initially ejected from the left sidewall.When the time was 80 ms,the velocity range of newly ejected rock blocks was 30-39 m/s,which is greater than before.It can also be observed that the rocks in the floor having a velocity of 15-20 m/s were moving into the opening.When the time was 100 ms and 120 ms,more rock blocks were ejected from the roof and sidewalls,and the maximum velocity was 35-45 m/s.The red arrows in Fig.21f present the movement direction of most rock blocks: from the right top sidewall,roof,left sidewall,and floor to the roadway center.

It might be argued that the simulated velocity of some ejected rocks is higher than that observed in the field(typically lower than 10 m/s).This can be partly attributed to the fact that no supporting systems were considered in the current simulation.Thus,the deformation energy of rock masses can be sufficiently released and transferred to the kinetic energy of ejected rocks.The ejection velocity is also positively correlated to PPVs.The higher the PPV is,the greater the ejection velocity is (Qin and Mao,2008).Besides,considering that the seismic data are not available,and the research’s objective is to show the effectiveness and advance of the global-local modeling approach in studying the damage mechanism of strainbursts,the seismic loading in this study is thus a simplified situation.The seismic waves were directly applied to the model’s right boundary.The distance from the seismic source to the excavation boundary is very small compared to the sites’ cases.Hence,the attenuation of seismic waves in numerical modeling is weaker than that in the field.Another cause is that there is no energy dissipation when two contact faces are divided.

Further studies (e.g.setting residual values of contacts or selecting more representative constitutive models)will be done by authors to address this problem.However,it should be noted that the average velocity of ejected rocks is 16.6 m/s after conducting a statistical analysis.In addition,some researchers have confirmed that the ejection velocity with an order of 10 m/s or greater is not unusual (Ortlepp and Stacey,1994; McGarr,1997).

In summary,the rock blocks with a high velocity during a strainburst might be ejected from all directions oriented to the roadway profile even the seismic wave originates from a specific direction.The rapid bulking or heaving of the floor can occur,which needs to be paid much attention to since the treatment of the floor is usually the most overlooked part in rockburst support.

Fig.19.The maximum principal stress (in MPa) surrounding the main drift during four excavation steps: (i) Excavation of the main drift; (ii) Excavation of the Panel H-6; (iii)Excavation of a part of roadways; and (iv) The 10th excavation step of the tailgate.

Fig.20.(a)Variation of surrounding rock stress during four excavation steps;and(b)Variation of stress concentration coefficient and the average stress in the stress increased area during four excavation steps.

Fig.21.Simulated velocity vectors of rock blocks during the occurrence process of a strainburst.

Fig.22 shows the evolution of cracks and macroscopic failure patterns during the strainburst.Fig.23 shows the variation of damage levels with the dynamic calculation time.When the time was 20 ms,there was no noticeable increase in tensile and shear damages than in the static stage.From 20 ms to 40 ms,the tensile damage suffered a rapid growth while the shear damage decreased slightly because some shear cracks disappeared due to the detachment of rock blocks.When the time was 40 ms,tensile cracks in the floor continued to initiate,propagate,and develop,and gradually formed more macroscopic fractures.Besides,tensile and shear damages also developed into deep areas within the roof and floor.Some ejected rock blocks from the right sidewall were observed.After 40 ms,both tensile and shear damage increased rapidly.From 40 ms to 60 ms,the tensile damage was increased from 1.17%to 1.49%,with a growth rate of 27.35%.The shear damage was increased from 1.71% to 3.14%,with a growth rate of 83.62%.When the time was 60 ms,more rock blocks were ejected from the right sidewall,macroscopic fractures began to occur in deep areas,and the floor was more fractured than before.After that,the tensile damage continued to increase almost linearly,while the shear damage also grew significantly,although it might suffer a few fluctuations.When the time was 80 ms,new macroscopic fractures were observed in the roof with the propagation and development of microcracks,and some rock blocks were ejected from the left sidewall.When the time was 120 ms,tensile and shear damage levels were 2.57% and 4.68%,increased by 119.66% and 173.68%,respectively,compared with those at 40 ms.As shown in Fig.22d,the damage mechanism for this case is the combination of three types of damage:rock ejection,rock bulking,and rockfall,which is consistent with the complexity of damage mechanisms of many violent rockburst events(Cai,2013;Manouchehrian and Cai,2018;Lu et al.,2019;Vazaios et al.,2019).Hence,rock supporting systems must possess multiple functions,such as dissipating kinetic energy and retaining ejected rock fragments,to resist and mitigate different types of rockburst damage.

Fig.22.Evolution of cracks and macroscopic failure patterns of the main drift during the occurrence process of a strainburst.

Fig.23.Simulated damage evolution of the main drift during the occurrence process of a strainburst.

Additionally,it is interesting to note that the range of the macroscopic failure zone is the same as that of the tensile damage(see Fig.22a and d).Thus,the initiation,propagation,and development of tensile cracks play a key role in controlling macroscopic failure of surrounding rock masses,although the shear crack always accounts for the main proportion of damage levels during the entire occurrence process of the strainburst.This finding agrees that the wave impedance of rock masses and air differ greatly that most stress waves are reflected at the surface,resulting in the tensile failure of surrounding rock masses(Wu et al.,2019b).

It is also interesting that many shear cracks(see purple cracks in Fig.22)accumulate within clay shale.This can be attributed to the great difference in wave impedance between coal seams(409/3 and 409/4)and clay shale.The greater the contrast in wave impedance,the smaller the seismic energy transmitted through the rock layer interface.Besides,because of the constraint of two coal seams,the“trapped”seismic energy in clay shale is difficult to transfer to the kinetic energy of ejected rocks as observed near the excavation boundary.Thus,the shear sliding of contacts occurs to consume the seismic energy in clay shale.

The comparison between the simulated failure patterns and in situ observations is shown in Fig.24.It can be seen that the rockfall and floor heaving induced by a strainburst can be realistically captured by the local model,suggesting the rationality and capability of the used method in modeling strainburst damage.However,it should also be noted that much more rocks are ejected from roof and sidewalls in simulation than in the field.This is because no support(steel arch and mesh)was modeled at this stage.Thus,the stress can be fully released from the roof and ribs and ejected rock materials are not held by any support elements.

Fig.24.Simulated failure patterns of a strainburst and field observations of rockburst damage.

Fig.25 shows the contour of maximum principal stresses during the strainburst.The significant stress concentration (blue areas)mainly appeared in two bottom corners and the roof of the main drift when the time was 20 ms.There were no significant changes in stress distribution compared to that in the static stage.This is because most surrounding rock masses were not affected by dynamic stresses during the early stage of the strainburst process.When the time was 40 ms,the significant stress concentration zone (blue areas) in the roof and bottom corners was enlarged because clay shale and sandy shale have relatively high strength,and a part of strata did not fail in this stage prone to stress concentration.After that,significant stress concentration was no longer observed in the roof with the further development of cracks and fractures,while it gradually developed in the deep floor surrounding fractured zones.When the time was 120 ms,a large stress relaxation zone(orange and red areas)was produced due to rock masses’ ejection,bulking,and fracturing.The significant stress concentration mainly occurred in the floor and has a“u” shape.Therefore,it can be concluded that the stress concentration gradually transfers from the roof and bottom corners to the floor.As shown in Fig.25,the macroscopic fracturing mainly happened in the floor’s shallow part.Thus,the rock masses in deep floor areas were competent,which provides an environment for further stress concentration.Hence,it can be anticipated that the rapid rock bulking in the floor will continue to develop if the dynamic time is sufficient or the floor is not hard.This anticipation agrees with many facts that significant floor heaving(up to 2 m)is often observed in rockburst events induced by seismic waves(Mutke et al.,2009; Stacey and Rojas,2013; Prusek and Masny,2015).

Fig.25.Distribution of the maximum principal stress surrounding the main drift during a strainburst.

The microscopic mechanism of crack initiation and development can be further explained based on the stress analysis.Rock masses fail predominantly by extensile fracturing (e.g.spalling)under no or low confining pressure (Gao et al.,2019b).Hence,during the early stage of the strainburst,the tensile failure of cracks mainly occurs surrounding the opening boundary because of the stress relaxation after excavation (see orange and red areas in Fig.25).The initiation of tensile cracks is constrained at depth due to the increasing confining pressure(stress concentration zone,see blue and green areas in Fig.25),but the shear fracturing of rock masses happens under high confining pressures.Thus,the depth of shear damage is greater than that of tensile damage.With the expansion of the stress relaxation zone,more cracks fail in tension in this area.The range and depth of shear cracks also grow because the stress concentration zone transfers to deep areas around the tensile failure zone.

4.2.2.Influence of PPVs

(1) Displacement analysis

Figs.26 and 27 show the displacement patterns of the main drift affected by different PPVs.The roadway deformation aggravated due to dynamic stresses.The magnitude of floor heaving and sidewall convergence rise significantly with the growth of PPVs,while the roof subsidence keeps relatively constant.However,the primary deformation forms are roof subsidence and the shrinkage of two sidewalls when the PPV is less than 0.8 m/s.After that,the floor heaving accounts for an essential part of the deformation.The floor heaving value is increased from 528 mm (PPV=0.8 m/s) to 920 mm(PPV=1.4 m/s),and the ratio of floor heaving to the total deformation is increased from 12.6%to 18.5%,with a growth rate of 46.8%.Hence,the rock support should consider the data of seismic activities in this area.Treatment might not be adopted in the floor when the PPV is low.However,additional support should be conducted in the floor when the roadway is located in areas where seismic events are active and often have a great magnitude.Two alternative measures are rockbolting the floor and using inverted arches (Yang et al.,2017).In this regard,numerical modeling is suggested to evaluate the performance of support measures during the design stage as it has benefits of low cost,safety,time-saving,and flexibility.

Fig.26.Distribution of displacement and velocity of the main drift affected by different PPVs.

Fig.27.Comparison of the deformation in different positions of the main drift affected by different PPVs.

Besides,it can also be observed that the difference between the displacement in two sidewalls is more significant with the increasing PPVs.For instance,the displacement in the right sidewall is 905 mm more than that in the left sidewall(PPV=1.4 m/s),while the displacement difference between the two sidewalls is 690 mm (PPV=0.2 m/s),suggesting that a more non-uniform displacement pattern occurs when the seismic event has a large magnitude.This phenomenon also indicates that the sidewall deformation is strongly correlated with seismic waves’ incident direction if the PPV is great.It can be noted that there are several fluctuations in the variation of roadway deformations.This can be attributed to the fact that some monitoring points are located at rock blocks that have been detached and might be collided by others.Nevertheless,the roadway profile continuously shrinks with increasing seismic magnitudes (Fig.26).Therefore,the support strength should be high enough to resist the large deformation induced by dynamic stresses in areas where seismic activities occur frequently.

Fig.26 also shows the velocity distribution of rock blocks during the 10th excavation step of the tailgate (step iv).The velocity of some rock blocks can reach 20-40 m/s even when the PPV is 0.2 m/s.These rock blocks are mainly located at the right sidewall.With the increase of PPVs,the range of the region with high velocities grows.Additionally,this region is gradually enlarged from the right sidewall to the roof and the left sidewall,suggesting that the areas with high rockburst risks are more significant with the increasing PPVs.Thus,it can be concluded that the asymmetric support scheme could be used when the seismic source is known (e.g.a nearby minor fault) and the monitored PPVs are usually low.However,if the monitored PPV is often great or the seismic source is unknown,equal support strength should be adopted in both sides.

(2) Fracture evolution and damage analysis

Fig.28 shows cracks in roadway surrounding rock masses concerning different PPVs.The comparison of damage levels is shown in Fig.29.Both the extent and level of damage grow with the increasing PPVs.The tensile damage mainly occurs in roof and two sidewalls when the PPV is small (less than 0.6 m/s).Tensile cracks are concentrated in floor’s shallow part,and the failure depth is less than 1.2 m.When the PPV is significantly more than 0.6 m/s,tensile cracks begin to appear in deep areas of the floor,and the damage depth can reach 2.4 m (PPV=1.4 m/s).Compared with tensile damage,the shear damage occurs in the roof,floor and two sidewalls,and its extent and depth are large.For instance,when the PPV is 1.4 m/s,the tensile cracks’depth in the roof is about 4.3 m,while the shear cracks’ depth is 5.5 m.The tensile cracks’ depth in the floor is 2.4 m,while the shear cracks’ depth is 3.7 m.As shown in Fig.29,the shear damage is always the dominant failure mode and increases significantly with the growth of PPVs.The tensile damage decreases first and then grows.The shear to tensile damage ratio increases from 1.25 to 1.82,suggesting that more slipping of rock masses occurs with the increasing dynamic stresses.As a result,the roadway deformation becomes more ununiform(see displacement vector maps in Fig.26).Therefore,the cable bolt is recommended in burst-prone areas since it can provide high pre-tension stress to strengthen fractured rocks,resist the shear deformation,and have greater control extent than regular rockbolts.

Fig.28 also shows the macroscopic failure patterns of the main drift.When the PPV is less than 0.6 m/s,macroscopic fractures are mainly observed in the roof and two sidewalls.With the growth of PPVs,more macroscopic fractures appear in the floor and gradually extend into deep areas.The floor heaving is more obvious when the PPV is higher.In conclusion,the damage mechanism of the strainburst is the combination of three types of damage: rock ejection,rock bulking,and rockfall when the PPV is high,while there is no obvious rock bulking when the PPV is low.Therefore,the objects of rock support can be adjusted according to the analysis and summarization of a large amount of seismic data in a specific area.However,it is always conservative and safe for engineering practices to consider all the damage mechanisms for support designs.Additionally,the macroscopic fashion of surrounding rock masses is still the tensile failure,and it mainly occurs in shallower areas compared with the shear failure.Hence,the rockbolts and cable bolts with enough pre-stress should restrain the initiation and development of tensile cracks,thereby reducing damage during a strainburst.

Fig.28.Distribution of cracks and macroscopic failure patterns of the main drift affected by different PPVs.

Fig.29.Comparison of the damage levels of the main drift affected by different PPVs.

(3) Energy evolution

The severity of strainbursts is related to the magnitude of the kinetic energy of ejected rock materials.The kinetic energy is one part of the total released energy that the whole supporting system(e.g.rockbolt,cable bolt,liner and steel mesh) must absorb to reduce rockburst risks (Raffaldi et al.,2017).Therefore,the influence of PPVs on the change and distribution of kinetic energy was investigated in this study.The kinetic energy of rock blocks was captured by the FISH language programing in UDEC using the following formula:

where m and V are the mass and velocity of rock blocks at the current time step,respectively.

Fig.30 shows the variation of kinetic energy with time influenced by different PPVs.The peak kinetic energy positively correlates with PPVs,indicating that seismic waves with higher PPVs cause severer rockbursts.However,the variation of kinetic energy with time depends on PPV values.The kinetic energy almost linearly increases with time from 0 ms to 87 ms and is gradually stable(PPV=0.2 m/s).When the PPV is higher than 0.2 m/s,the kinetic energy first grows to a peak value and then gradually decreases to a certain extent with time.The drop in kinetic energy is because some initially ejected rock blocks have been settled on the floor.It can also be observed that the higher the PPV is,the less the time needed to reach the peak kinetic energy,which is reasonable because of the equal distance from the seismic source to the main drift.This phenomenon suggests that rock ejection that is influenced by higher PPVs is sudden and violent,while that of lower PPVs is relatively slow.

Fig.30.Comparison of the kinetic energy of the whole rock system influenced by different PPVs.

The distribution of kinetic energy influenced by different PPVs is shown in Fig.31.The kinetic energy pattern is similar to velocity(see Fig.26).The kinetic energy of some rock blocks can exceed 14 kJ even when the PPV is 0.2 m/s.These rock blocks are mainly ejected from the right sidewall.This phenomenon matches some local strainburst damage observed in the field (e.g.ejected rocks and rockbolts from a limited area (Zhang et al.,2012,2013)).With the increase of PPVs,the range of the region with great velocities grows,and thus more rock blocks possess high kinetic energy.Besides,this region is also gradually enlarged from the right sidewall to the roof and the left sidewall,which indicates that the areas with high rockburst risks are expanded with the increasing PPVs.Assume that the tensile yield force of resin-grouted rebar is 160 kN,and the maximum allowed deformation is 25 mm(Stillborg,1994).The energy absorption capacity for this rebar is only 4 kJ,which is 17% of the kinetic energy (23.4 kJ) of a 0.3 m × 0.3 m × 1 m coal block with an ejection velocity of 20 m/s.Thus,the ejected coal block has a velocity of 16.6 m/s,which is still very harmful to mine workers and equipment.Therefore,yielding support elements(e.g.He-bolt,Roofex,and D-bolt) are essential and should be used to absorb the ejected rocks’ kinetic energy effectively.

Fig.31.The kinetic energy (in J) distribution of the whole rock system influenced by different PPVs.

It should be noted that the detached blocks might penetrate too far into another in numerical modeling,especially in the dynamic stage,which can cause too great overlap of contacts.Some alternative approaches to handle this problem can be referred to Itasca(2020).

5.Suggestions for rockburst supporting

5.1.Evaluation of the current support system

Rockburst support systems must resist dynamic loads and great deformation caused by rock fracturing,dilation,and ejection,which are different from conventional or standard supports aiming to control gravity-induced rockfalls and manage shallow areas of loose rock (Kaiser and Cai,2012).Kaiser and Cai (2012) and Cai(2013) summarized three key and indispensable functions of rockburst supporting: reinforcement,retaining,and holding.The aims of reinforcement(e.g.employing fully grouted rebar and cable bolts) are to prevent the propagation and expansion of fractures and to strengthen rock masses by increasing their cohesion and friction angle,thereby producing a bearing structure or anchorage body to support itself and outside rock masses (Gao et al.,2009).Retaining elements (e.g.wire mesh,steel arch,and shotcrete) are used to avoid unraveling fractured rocks.The holding function can be fulfilled by tying retaining elements and loose rocks back to stable areas in depth to dissipate the kinetic energy induced by rock ejection and to prevent gravity-induced rockfall.Yielding support elements should be used since they can accommodate the great deformation and absorb kinetic energy.

The main drift was supported by yielding steel arches and welded wire mesh.These support elements can accomplish the retaining function,and the shrink characteristic of the yielding steel arch (fulfilled by the frictional sliding between connecting pieces and steel arches) can absorb an amount of kinetic energy(Horyl and Snuparek,2012).However,the combination of steel arch and wire mesh is a passive support fashion.The high stiffness of steel arches can limit roadway convergence to a certain extent,but the excavation-induced fracturing is not effectively restrained due to the lack of active compressive force and reinforcement acting on fractured and deteriorated surrounding rock masses (see Fig.16).Although steel arches and wire meshes can retain the loose rocks,the bearing load of these elements is significantly increased because they are subject to the whole loads of loose rocks that are not held with deeper stable rock masses.Therefore,the capacity of energy absorption of the yielding steel arch might not be sufficient enough to dissipate kinetic energy during a rockburst when a part of capacity has been consumed owing to the large deformation and loads of fractured rocks because of high in situ stresses,mininginduced stresses,and poor rock conditions (Gao et al.,2009).The kinetic energy from ejected rocks will then be transferred to the high strain energy accumulated in yielding steel arches that might be overloaded and lose their functions.This is why severe rockburst damage is still widespread in many deep roadways that use yielding steel arches and wire meshes as the main support elements (Fig.32).

Fig.32.Pictures showing severe failure of steel arches,damaged belt conveyer,rockfall,and floor heaving in deep roadways after rockbursts:(a)From Horyl and Snuparek(2012);and (b,c) From Mutke et al.(2009).

5.2.Principles of rockburst supporting

It is also imperative to figure out the possible rockburst damage mechanisms to select reasonable support elements for accomplishing three critical functions of rockburst supporting(Kaiser and Cai,2012).The used modeling approach in this study can anticipate the damage mechanisms of strainbursts (the type,location or range,level,and energy)from both micro and macro perspectives,making it possible to improve and optimize the support design.Taking the case study site as an example,the following principles of rockburst supporting are proposed according to the analysis of simulation results of strainburst damage mechanisms (see also Fig.33):

Fig.33.Principles of rockburst supporting based on the analysis of strainburst damage mechanisms.

(1) Support area.The seismic activities in planned construction areas should be recorded and analyzed in advance.If the seismic source is known (e.g.a nearby minor fault) and monitored PPVs are usually low,treatment might not be adopted in the floor,and the asymmetric support scheme could be used.Although the asymmetric support scheme is a common practice in entries or drifts in longwall mining,it should be used with caution considering the greater risk and more uncertainties in the burst-prone ground.The feasibility and effectiveness of this strategy needs to be verified with more cases.In contrast,when the roadway is located in areas where seismic events are active and often have a great magnitude or the seismic source is unknown,the equal support strength should be used in both two sidewalls,and the floor needs to be treated with additional support measures to avoid rapid bulking or floor heaving.

(2) Using yielding rockbolts.Kaiser at al.(1996)reported that in rockburst-prone mines,ejection velocities below 1.5 m/s could be handled by standard ground support,but additional support is required for higher velocities.The yielding rockbolts with high strengths (e.g.D-bolt) are recommended.This type of yielding rockbolts can bear high loads and accommodate large deformation during rockbursts.Also,it can absorb more kinetic energy than those with low strength and conventional rockbolts (Fig.34).Enough pre-stress should be applied to surrounding rocks when installing rockbolts to control separation and fracturing induced by the development of tensile cracks,and to increase the confining pressure,which has been proven to be effective in improving the physical and mechanical performances of fractured rocks(Yang et al.,2017).The use of yielding rockbolts can fulfill reinforcement and a part of the holding function.

Fig.34.Comparison between the conventional rockbolt and yielding rockbolts with low and high strengths (after Li,2020).

(3) Using cable bolts.The cable bolts with high strengths and pre-stress should be used to restrain the development of cracks and resist large deformations,thereby strengthening surrounding rock masses and maintaining their integrity(reinforcement function).Cable bolts usually have a great length,which can hold retaining elements and shallow rocks back to stable areas in depth.Besides,the cable bolts also have a relatively high elongation rate (4%-7%,Kang et al.,2009) to absorb the deformation energy induced by rock bulking.

(4) High-strength retaining elements.The wire mesh with high tensile strength is recommended since this mesh type can deaccelerate ejected rocks and absorb kinetic energy effectively because of its strength and flexibility (Roth et al.,2007).

(5) Improving the connection between rockbolts and meshes.Cai (2013) proposed that the support system often loses its effectiveness due to the weakness of the bolt-mesh linkage rather than the insufficient capacity of rockbolts.Thus,he recommended that the relatively large plates (minimum 150 mm × 150 mm) should be used to connect rockbolts to wire meshes.

6.Conclusions

This paper presents an improved global-local modeling approach to study strainburst damage mechanisms.The extracted stresses induced by multiple nearby excavations from a calibrated 3D FDM global model are used as boundary conditions for a 2D DEM local model of a roadway.The damage mechanism of strainbursts is investigated from both micro and macro perspectives.The main conclusions are as follows:

(1) The used approach provides a more insightful understanding of the influences of the realistic stress loading path on strainburst occurrence.The roadway surrounding rock stress increases remarkably because of the superposition of excavation-induced stresses of nearby panels and roadways.This results in more accumulated strain energy in two sidewalls,providing a necessary condition for the strainburst occurrence in the dynamic stage.Having a better understanding of the realistic stress loading path can also help to identify stress concentration zones and adopt distress measures in time to mitigate rockburst risks.

(2) The strainburst damage mechanism of a violent rockburst event in the Zofiówka Mine,including rock ejection,bulking,and rockfall,is successfully captured in the local model,which confirms the rationality and capability of the modeling approach.The rock mass fracturing or damage during the strainburst is investigated quantitively using a self-defined damage variable in UDEC.During the strainburst,the initiation,propagation,and development of tensile cracks play a crucial role in controlling macroscopic failure of surrounding rock masses,although the shear crack always accounts for the main proportion of damage levels.The rockbolts and cable bolts with enough pre-stress should be used to restrain the initiation and development of tensile cracks,thereby reducing the level of damage.

(3) The deformation and damage level of the roadway during the strainburst have a positive correlation with the increasing PPVs.A more non-uniform displacement pattern occurs when the seismic event has a large magnitude.The asymmetric support scheme is an alternative strategy when the seismic source is known and the monitored PPVs are usually low.

(4) The yielding steel arch might not dissipate kinetic energy and mitigate strainburst damage effectively owing to the limited capacity of energy absorption.Based on the findings in this paper,the principles (Fig.33) to control and mitigate strainburst damage are proposed.

The improved modeling approach in this study provides an effective and potential tool for understanding strainburst damage mechanisms and improving rockburst supporting.Further studies(e.g.using available seismic data,modeling support elements,setting residual values of contacts,and selecting more representative constitutive models) will be done by authors to improve the accuracy and robustness of the approach.

Declaration of competing interest

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

Acknowledgments

The authors are grateful to Itasca Consulting Group,Inc.and Itasca Educational Partnership for permitting the software usage.Support from China Scholarship Council is also acknowledged by the first author.