Sub-Homogeneous Peridynamic Model for Fracture and Failure Analysis of Roadway Surrounding Rock

Shijun Zhao ,Qing Zhang ,Yusong Miao ,Weizhao Zhang ,Xinbo Zhao and Wei Xu,⋆

1School of Science,Qingdao University of Technology,Qingdao,266520,China

2College of Mechanics and Materials,Hohai University,Nanjing,211100,China

3Xibei Mining Co.,Ltd.,Shandong Energy Group,Xi’an,710018,China

ABSTRACT The surrounding rock of roadways exhibits intricate characteristics of discontinuity and heterogeneity.To address these complexities,this study employs non-local Peridynamics(PD)theory and reconstructs the kernel function to represent accurately the spatial decline of long-range force.Additionally,modifications to the traditional bondbased PD model are made.By considering the micro-structure of coal-rock materials within a uniform discrete model,heterogeneity characterized by bond random pre-breaking is introduced.This approach facilitates the proposal of a novel model capable of handling the random distribution characteristics of material heterogeneity,rendering the PD model suitable for analyzing the deformation and failure of heterogeneous layered coal-rock mass structures.The established numerical model and simulation method,termed the sub-homogeneous PD model,not only incorporates the support effect but also captures accurately the random heterogeneous micro-structure of roadway surrounding rock.The simulation results obtained using this model show good agreement with field measurements from the Fucun coal mine,effectively validating the model’s capability in accurately reproducing the deformation and failure mode of surrounding rock under bolt-supported(anchor cable).The proposed subhomogeneous PD model presents a valuable and effective simulation tool for studying the deformation and failure of roadway surrounding rock in coal mines,offering new insights and potential advancements.

KEYWORDS Roadway surrounding rock;peridynamics;heterogeneous material;fracture analysis;numerical simulation

1 Introduction

The surrounding rock of roadways is a composite,heterogeneous body comprising mineral particles and various defects formed through geological processes.It includes discontinuous structures such as fractures,joints,bedding,and faults [1].These characteristics result in distinct mechanical properties,deformation,and failure behaviors compared to homogeneous materials.Under external loads and other disturbances,pre-existing micro-defects expand and aggregate,forming new macrocracks and energy dissipation.Consequently,the strength and stiffness of the coal-rock mass are significantly reduced,ultimately leading to irreversible fracture failure,extensive roadway deformation,damage to support structures,and even overall roadway instability [2].Over recent decades,advances in computing technology have made numerical methods an essential and valuable tool for investigating rock mechanical damage and failure.However,a key challenge in underground engineering is understanding the mechanical mechanisms behind the deformation and failure of surrounding rock in such scenarios and developing an effective numerical simulation approach capable of accurately capturing these processes.

In order to enhance understanding of roadway surrounding rock instability mechanisms,many researchers have conducted extensive research using numerical methods to investigate damage and failure characteristics.Traditional numerical techniques,such as the finite element method (FEM),require predicting crack location and propagation direction.When cracks propagate and branch,the mesh at the crack tip must be redefined,leading to pronounced mesh dependence in the simulation results.Introducing Goodman joint elements into the FEM framework allows describing discontinuous structures,such as joints and fissures within rock masses.However,numerical instability issues may arise when many joint elements are involved.To address these challenges,Belytschko et al.[3]proposed the extended finite element method (XFEM),introducing discontinuous and crack tip displacement field functions into the conventional finite element displacement mode.However,XFEM requires that cracks remain continuous at the interfaces of adjacent elements,which still makes it challenging to deal with complex failure problems like multiple crack interactions and branching.Another approach,the boundary element method (BEM) [4],treats discontinuous interfaces,such as cracks,as boundary interfaces.However,basic solutions for corresponding differential operators are not guaranteed,limiting their applicability.The discrete element method (DEM) [5] and other discrete medium numerical methods are widely used in geotechnical engineering.Nevertheless,their accuracy is relatively modest,the number of contact simulations is extensive,and numerous artificial assumptions must be introduced,potentially leading to significant discrepancies between simulation results and analytical situations.The numerical manifold method(NMM)[6]employs a finite element mesh as its mathematical coverage.During the process of crack propagation,this mathematical coverage remains unchanged.It is only necessary to place densely distributed nodes at the crack tip without regenerating the mesh.This approach can effectively simulate the progressive cracking process.However,the function corresponding to the control equations of NMM does not always exist,which limits its scope of application.The mesh-free method(MF)[7],utilizing node information and weight functions within compact support domains for local approximation,requires the addition of nodes near the crack tip to facilitate fracture simulation.Although MF methods eliminate the mesh dependency issue,they require the introduction of additional parameters.Moreover,their higher-order continuous approximation functions do not necessarily offer advantages when solving crack propagation problems,and accurately applying essential boundary conditions can be quite challenging.

From the perspective of the mechanism underlying material deformation and failure,coalrock geotechnical materials comprise various micro-defects at the mesoscale,resulting in distinct mechanical behaviors for each mesoscale unit under external loads.The interaction among these mesoscale units within the coal-rock mass leads to an uneven internal stress distribution,producing nonlinear deformation characteristics at the macro-scale.Existing numerical modeling methods for heterogeneous geotechnical materials are complex in modeling and simulation,and traditional numerical methods face limitations in addressing discontinuous problems[8,9].The primary challenge of traditional numerical methods in simulating fracture and related issues stems from discontinuities caused by cracks.Nevertheless,within the framework of continuum mechanics,various difficulties are inevitable.Peridynamics (PD) offers a multi-scale simulation approach capable of overcoming differential operation issues associated with singular or discontinuous nodes in continuum mechanics.PD has surmounted the limitations of spatial scale.Currently,the mechanical constitutive relationships of materials mostly rely on laboratory-scale specimen tests,which do not provide a unified constitutive model between the micro-scale and macro-scale.In contrast,PD utilizes a unified model to describe deformation and failure problems from the atomic scale to the macro-scale.This approach circumvents the complexities of traditional multi-scale mechanics methods for transferring mechanical quantities across different scales.It enables cross-scale numerical computations within a unified framework[10].Other researchers,such as Zhou et al.[11,12],have integrated PD theory into the failure analysis of rock-like materials,examining crack growth under various loading conditions.Ma et al.[13]proposed an improved PD model to capture the stress-strain changes in rock materials and simulated the changes in cracking angles under different pre-crack angles and tensile load conditions.Moreover,Li et al.[14]used the PD method to simulate and study crack propagation and penetration in rocklike materials.Gao et al.[15] implemented a material node dormancy method based on PD to simulate roadway excavation,exploring the deformation and failure characteristics of surrounding rock during this process.Wang et al.[16–18]defined nonlocal stresses to obtain isotropic and deviatoric forces commensurate with deformations and proposed a new failure model to link computational PD with some phenomenological failure criteria,highly relevant for both brittle and quasi-brittle rocks.Furthermore,Fan et al.[19–21] conducted a hybrid PD-SPH (smoothed particle hydrodynamics)approach to simulate geotechnical materials fragmentation under buried explosive loads.However,the studies above on the failure of rock-like materials have not adequately considered material heterogeneity.

This study presents the construction of a quadratic constitutive force kernel function that effectively captures the spatial variation of long-range force.Departing from the micro-structural characteristics of coal-rock materials,this research addresses the heterogeneity of such materials and introduces a novel heterogeneous PD model centered on random pre-breaking “bonds.” The improved model treats the coal seams,roof,and floor layers as multi-layer composite materials.The incorporation of material heterogeneity is accomplished by establishing a pre-breaking “bond”coefficient,enabling the representation of the random heterogeneous mesostructure of surrounding rock in roadways and considering support effects.Subsequently,a corresponding numerical simulation method is developed to effectively analyze the deformation and failure of surrounding rock in extensive discontinuous and heterogeneous geotechnical engineering scenarios.The proposed sub-homogeneous PD model holds great potential in accurately analyzing complex geotechnical systems exhibiting significant discontinuities and material heterogeneity.

2 Peridynamics Theory

2.1 Bond-Based Peridynamic Model

As shown in Fig.1,at a specific time t,there is an interaction force represented by a vector functionfbetween any material nodexoccupying a certain space domainRand all other material nodes within a certain range of the surrounding space δ,then the internal force f caused by the interaction between the nodes within the rangeHxof any node is,

In the space domainR,the equation of motion of PD at any time can be expressed as follows:

Figure 1:Schematic diagram of interaction between nodes in PD theory

The radiusδis called the horizon range,Hxis the collection of all the material nodes within the rangeδ,specifically,

When the distance between the material nodes>δ,the interaction between the material nodes no longer occurs,

whereξ=-xis the relative position vector,η=-uis the relative displacement vector,andf(η,ξ)is the constitutive force function between nodesxand.

For homogeneous materials,the general expression of the constitutive force functionf(η,ξ)between nodes is,

wheref(η,ξ)is a scalar value force vector function,which contains information that can be utilized to describe the mechanical response of materials,and the functionf(η,ξ)must meet,

In addition to the relevant material parameters,the functionf(η,ξ)is only related to|η|,|ξ|and the angle between them.Based on Stokes theorem,there must be differentiable scalar functionswsatisfying,

The scalar functionw(η,ξ)is the potential energy function between nodes,representing the energy densityξstored in the“bond”.

For the prototype micro elastic brittle(PMB)material model proposed by Silling et al.[22],the mechanical properties of the material are related to the micro-moduluscand critical elongations0.PMB model is based on modeling interactions between material nodes with spring-like“bonds,”where these“bonds”can only be stretched and do not consider lateral or torsional deformations.

The potential energyw(η,ξ)stored in the “bond”between nodes,namely,the density of deformation energy,is defined as follows:

The deformation energy density of the material nodexin its horizonδis,

The coefficient indicates that the deformation energy density of each node occupies half of the deformation energy density of related two nodes.

Based on the PMB model,the pairwise force function of micro-elastic material is derivable from the scalar micro potentialw,

whereμ(t,ξ)is a time-related scalar function;sis the elongation of the“bond”between nodes.

In PD,local damage is represented by the ratio of the number of broken bonds to the total number of bonds in the horizon of each node,

2.2 Modification of Constitutive Force Function

In the constitutive function of the PMB model,the interaction between nodes changes in direct proportion to the elongation |ξ| of the “bond” without considering the effect of the length of the“bond” on the interaction force between nodes.The PMB model can be modified by introducing different forms of kernel functions into the constitutive function to eliminate the numerical dispersion of the PD equation.

For the modified model,the constitutive force functionf(η,ξ)can be expressed as follows:

The kernel functionκ(ξ,δ)can represent the characteristics of the spatial decline of the long-range force between nodes.In this study,κ(ξ,δ)was defined as follows:

The corresponding micro-modulus functionc(ξ,δ)can be expressed as follows:

The relationship between micro-modulus functionc(ξ,δ)and mechanical parameters in continuum mechanics can be derived utilizing equivalence between the strain energy density of nodes in PD and the strain energy density in continuum mechanics.For the quadratic polynomial type kernel function used in this paper,its micro-modulus function is as follows:

For 2-D problems,the fracture energyG0in PD is given by,

Substituting Eq.(18)into(17),the approximate expression of the energy release rateG0in PD can be obtained as follows:

An approximate expression of the energy release rateG0can be derived as follows:

Then,the critical elongations0in the modified model of the 2-D problem can be obtained as follows:

2.3 Numerical Solution Method

In the PD simulation system,the model is composed of material nodes.The motion of nodes obeys Newton’s second law.The mechanical response of the model can be obtained by solving the integral equation.Solving the motion equation of the PD simulation model involves time integration(difference scheme),spatial discretization,and integration.The unit volume force generated by the nodexiin the reference configuration under the action of other nodes within its horizonδcan be obtained by summing,and the spatial integral equation can be discretized as follows:

The numerical integration of time can be obtained by the Velocity-Verlet difference scheme.The displacementu(x,t)and velocityof a node are as follows:

The force boundary conditions in traditional theory cannot be directly applied in the PD model,and the external force usually needs to be applied through volume force density.The concentrated force and surface forces are applied on several layers of material nodes at the boundary.If the model is subject to the external forcep,it can be converted into the external force density in PD.

3 Sub-Homogeneous PD Model

3.1 PD Modeling of Heterogeneous Coal-Rock Mass

A substantial number of experiments in rock mechanics indicated that coal-rock geotechnical materials display significant variability in their physical and mechanical parameters,and these properties also demonstrate spatial heterogeneity.The heterogeneity in coal-rock mass materials arises from various factors.Characterizing this heterogeneity during numerical modeling is exceptionally challenging.Many researchers argue that heterogeneity due to micro-structural defects,such as pores and voids,is a crucial factor influencing the deformation and failure of coal-rock mass.Therefore,the geometric configuration and distribution of these micro-defects greatly influence the physical and mechanical properties of the materials.

Based on Eq.(13),ϕ(x,t)∈[0,1]local damage can also be defined as the local damage indexdindexof a material nodexat timet.

whereNbis the number of“bonds”cut off between the node and all other nodes in its horizon range;Nis the total number of“bonds”between the same node and all other nodes in its horizon range.

The material presents micro-defects,like pores and holes,leading to macro-level heterogeneity.In the context of the PD simulation,this heterogeneity appears through the absence of“bonds”between nodes,indicating local damage between them.This study introduces a random pre-breaking“bond”coefficient,symbolized by the initial damage quantity,into the uniform discrete simulation model to effectively characterize micro-defects,including material pores(as shown in Fig.2),and thus capture the material’s heterogeneity.This method is inspired by the analysis of failure in functionally graded materials by Chen et al.[23].

Figure 2:Schematic diagram of bond random breaking in heterogeneity PD model

Following this concept,this study defines the PD model’s heterogeneity as the simulation model’s initial damage amount.The process of setting the random pre-breaking“bond”is depicted in Fig.3.The pre-breaking“bond”index of one node is,

whereM(x)is the degree of heterogeneity of one node.Hence,the rate of random pre-breaking“bond”Mcis the critical heterogeneity of the material.

Figure 3:Flow chart of bond random pre-breaking

For the degree of heterogeneity of material nodesxiandxj,respectively,areM(xi)andM(xj),the probability of the“bond”ξ=xj-xibetween material nodes remains intact is,

For the uniform discrete model,assuming that the degree of heterogeneity in the model is consistent,i.e.,the degree of heterogeneity of each node isM(x),then the probabilitypthat the“bond”between nodes remains complete is as follows:

Then,the damage indexDindexof any node can be expressed as follows:

Compared to existing homogenized PD models,the model in this study can be referred to as a subhomogeneous PD model.In the heterogeneous PD model,homogeneous dispersion is performed to ensure uniform physical and mechanical parameters for materials based on their heterogeneity degree.Establishing the sub-homogeneous PD model involves setting the heterogeneity degree coefficient,symbolized by the pre-breaking “bond”.Importantly,no additional pre-processing is necessary for the simulation model,and the model effectively reflects the random distribution characteristics of material heterogeneity.Moreover,the model’s implementation is straightforward,and the efficiency of heterogeneous modeling is unaffected by the simulation model’s size or the number of nodes,making it especially suitable for large discontinuous and heterogeneous geotechnical engineering applications.

Literature[24]indicated that coal-rock mass in deep strata undergoes compaction during diagenesis,typically resulting in low porosity.Regarding the coal-rock material heterogeneity in this study,Fig.4 illustrates the material heterogeneity degree and the corresponding damage index represented by random pre-breaking “bonds.”It is important to note that the damage in the nephogram does not relate to damage caused by external loads on the model;instead,it represents the proportion of random pre-breaking“bonds,”quantifying the material heterogeneity degree.This serves as the initial condition of the model during the simulation process.

3.2 Numerical Solution Process

Based on the modified PD model,this study develops a simulation program for the subhomogeneous PD model using the FORTRAN language.The fundamental procedure is shown in Fig.5.The pre-processing module involves compiling the model discrete program in MATLAB,enabling the acquisition of crucial information,such as node coordinates.Then,in the post-processing module,the simulation file information is read using Ensight,permitting the output of simulation results,including damage and displacement.

Figure 4:Degree of material heterogeneity and corresponding damage index characterized by bond random pre-breaking

Figure 5:Flow chart of the sub-homogeneous PD model simulation program

4 Analysis of Surrounding Rock Excavation and Support of Roadway

4.1 Computational Model

This research focuses on the No.3 upper 411 haulage roadway at Fucun Coal Mine as the research subject.The simulation model’s dimensions are set to 60 m×40 m×0.1 m(as illustrated in Fig.6),using orthogonal uniform dispersion.The node lattice size,denoted as,is set at 0.1 m,resulting in a total of 232320 nodes.Considering the heterogeneity degree,symbolized by the random pre-breaking“bond” coefficient,and to avoid premature damage at the model boundary,a suitably increased horizon range size was selected asδ=.

Figure 6:PD simulation model of gob-side entry

The Poisson’s ratio for diagenetic minerals is generally 0.2–0.3,as shown in Table 1,cited from the literature [25].In bond-based PD,the Poisson’s ratio for plane stress problems is limited to 1/3,and for plane strain problems,to 1/4.Thus,the influence of Poisson’s ratio on simulation results is not significant.The problem is then simplified as a plane strain problem,with a Poisson’s ratio ofν=0.25 assigned.The physical and mechanical parameters of the material are outlined in Table 2.The haulage roadway,situated along the goaf,is approximately buried at a depth of-500 m.For most crystalline rock masses,the porosity is relatively small,often ranging from 0.1%to 10%.The surrounding rock stratum shows a relatively dense composition with low porosity,resulting in a material heterogeneity of 5%selected for this study.

Table 1:Elastic constants of several common diagenetic minerals

Table 2:Physical and mechanical parameters of coal seam and roof and floor rock layer

In this proposed model,the PD parameters of various interlayer materials are different,and the micro-modulus functionc(ξ,δ)and critical elongations0of composite “bond”can be obtained by corresponding reduction.The micro-modulus functionc(ξ,δ)and critical elongations0of the material are determined by the physical and mechanical parameters of the material.In order to reflect the difference between the two kinds of material,in this paper,Young’s modulusand fracture energy densityof the composite“bond”take the harmonic average of the mechanical parameters of two materials.

The composite “bond”micro-modulus functionand critical elongationcan be obtained as follows:

The macro strength and other mechanical parameters of materials with different degrees of heterogeneity are also different.Chen et al.[26],based on the wave propagation theory under the plane stress state,obtained the relationship between longitudinal wave velocity and elastic modulus of heterogeneous materials through a large number of numerical tests.

whereCLis the longitudinal wave velocity,EMis Young’s modulus,νis Poisson’s ratio,for plane stress problemsν=1/3,for plane strain problemsν=1/4.

The corresponding Young’s modulus values of materials under different degrees of heterogeneity can be obtained through Eq.(34),and the following relationship exists between the Young’s modulusEMof materials with different degrees of heterogeneity and the Young’s modulusEof materials withM=0 through fitting,

whereEis Young’s modulus of material when it is heterogeneousM=0.For coal-rock mass materials,it is the Young’s modulus of the original rock.

Substitute Eq.(35)into Eq.(30)to get,

Similarly,the fracture energyG0of materials with different heterogeneity also follows the following relationship:

whereG0is the fracture energy of material when the degree of heterogeneity isM=0.For coal-rock mass materials,it is the fracture energy of the original rock.

4.2 Consideration of Support Structure

The reinforcement effect of the bolt on the surrounding rock of the roadway is shown in Fig.7.Generally,the tensile strengthof the bolt and anchor cable is far greater than the tensile strengthRtof coal-rock mass,and the reinforcement effect of the bolt(cable)on surrounding rock can be evenly distributed to each node within its scope of action through conversion.Assuming that the support density of the anchor bolt (cable) isρb,that is,the number of anchor bolts (cables) in the unit area,and the pre-stress applied to the anchor bolt(cable)isp,the maximum loadPmaxfor the tensile failure of the solid after anchor bolt(cable)is applied can be obtained as follows:

whereAis the surface area of the roadway;ris the radius of the bolt(cable)body.

Generally,inclined bolts will be applied at the shoulder angle of the roadway,so Eq.(38)can be rewritten as,

Figure 7:Schematic diagram of bolt reinforcement

In PD theory,whether the“bond”cut off is judged by the relationship between the elongation of the“bond”and the critical elongation between nodes.If the critical elongations0of“bond”between coal-rock mass nodes before reinforcement is,the critical elongationof“bond”between coal-rock mass nodes after reinforcement is,

In the simulation process,the local damage of nodes can be obtained from Eq.(13).

4.3 Initial Balance of Model

The numerical simulation process for roadway excavation and support application encompasses three primary steps: (1) establishing self-weight stress balance in the model,(2) executing roadway excavation,and(3)applying support.The study employs the gravity field method to create the initial stress field balance within the simulation model.The model’s bottom vertical displacement is fixed,and the normal horizontal displacement of its sides is also fixed.A balanced state is attained by allowing the model to achieve equilibrium solely through self-weight.Displacement monitoring nodes are strategically positioned symmetrically along the model’s horizontal axis.Once the displacement at these nodes stabilizes,the model has attained its initial balance.Fig.8 shows the vertical displacement at selected monitoring points,revealing that the model reaches equilibrium after 15,000 simulation iterations.After equilibrium,the force density of all nodes in the model is extracted,using these values as the initial force densities for the nodes in the roadway excavation simulation.The node displacements generated during the initial stress field equilibrium are reset to zero,ensuring precise simulations in subsequent excavation and support application phases.

Figure 8:Vertical displacement curve of monitoring node during loading

4.4 Simulation Results and Analysis

The PD parameters,such as the micro-modulus function and critical elongation of nodes in the anchor bolt (cable) area,are determined using Eqs.(40) and (41).Figs.9 and 10 depict the initial roadway support scheme’s cross-section and plan view.The roof support comprisesΦ18/Q500 lefthanded bolts without longitudinal reinforcement,each 2.4 m long,andΦ17.8 × 6000 mm anchor cables.The two side supports consist of similar bolts,each 2.0 m long.However,this scheme severely damages the coal seam on both roadway sides.In order to improve support effectiveness,the spacing between side anchor bolts is reduced based on the original scheme,while the spacing between roof anchor bolts(cables)remains constant,as shown in Fig.11.The row spacing of anchor bolts(cables)remains as in the original scheme.

Figure 9:Cross section of original roadway support(unit:mm)

Figure 10:Scheme of original roadway support(unit:mm)

Figure 11:Cross section of modified roadway support(unit:mm)

Unlike traditional continuum mechanics methods,bond-based PD eschews stress and strain descriptions.The study’s analysis focuses on horizontal displacement changes of the roadway’s sides and the roof’s vertical displacement variations.Figs.12 and 13 show these displacements under the original support scheme: the coal side’s maximum deformation is 35 mm,and the coal pillar side is 59 mm,totaling a maximum deformation of 94 mm.The maximum vertical displacement of the roadway roof is 53 mm.

Figure 12:Horizontal displacement of roadway surrounding rock under original support scheme(unit:m)

Figure 13:Vertical displacement of roadway surrounding rock under original support scheme(unit:m)

Figs.14 and 15 illustrate the displacements under the enhanced support scheme.The maximum deformation of the coal side is 26 mm,and that of the coal pillar side is 42 mm,leading to a total maximum deformation of 68 mm.The maximum vertical displacement of the roadway roof is 51 mm.

Figure 14:Horizontal displacement of roadway surrounding rock under modified support scheme(unit:m)

Figure 15:Vertical displacement of roadway surrounding rock under modified support scheme(unit:m)

Accordingly,reducing the spacing between side bolts and increasing support density significantly mitigates the horizontal movement of the roadway sides.The enhanced support scheme effectively controls roadway deformation during its service period.Fig.16 visually demonstrates the enhanced scheme’s actual support effect.

Figure 16:Actual effect of roadway support

In order to assess the support scheme’s feasibility,monitoring the ground pressure on the support system and surrounding rock is imperative.Observing mine pressure in the stopping roadway involves various instruments to measure macroscopic ground pressure manifestations,such as displacements of the roof,floor,and sides and deformation of the coal seam and support structure.The observation sections’layout appears in Fig.17.

Table 3 presents statistical data on deformations observed in the roof and sides of the roadway in field observations.The results show a maximum deformation of 65 mm on the sides and a roof subsidence of 43 mm.Remarkably,the simulated deformation of the surrounding rock by the subhomogeneous PD model closely aligns with field measurements,meeting the engineering requirements for investigating the deformation and failure characteristics of the surrounding rock and optimizing the support scheme.

Table 3:Displacement statistics of measuring points of gob-side entry surrounding rock(unit:mm)

Figure 17:Layout of observation sections

5 Conclusion

This study simulates the fracture and failure process of roadway surrounding rock subjected to support application using an improved sub-homogeneous PD model.This investigation analyzes the impact of various factors,including kernel function,heterogeneity,and support application,on the fracture pattern and failure mode.

(1) A quadratic constitutive force kernel function was developed to account for the spatial variation of long-range forces.Moreover,considering the heterogeneity of coal-rock mass materials and structures,a random pre-breaking“bond”coefficient was introduced,characterizing the degree of heterogeneity in rock-like materials.Establishing this heterogeneous PD model,based on uniform dispersion and reflecting the random distribution characteristics of material heterogeneity,enabled the simulation of deformation and failure in heterogeneous materials and structures.

(2)Given the layered composition of coal seams,including its roof and floor layers,a PD model was developed for layered coal-rock mass structures.To address underground engineering challenges,this bond-based PD model adopted a method of initial balance.Additionally,a support facility application mode for the PD model was proposed,facilitating the analysis of underground engineering.The proposed PD model effectively captured the random heterogeneous meso-structure of the roadway surrounding rock while considering the influence of support measures.

(3)The PD method was employed to assess coal-rock masses’structural deformation and failure.By employing the sub-homogeneous PD model,the study explored the deformation and failure characteristics of roadway surrounding rock when supported by bolts (cables).The displacement distribution along both sides and the roof of the roadway was meticulously analyzed.A comparative assessment between the simulation outcomes and actual measurements revealed a high degree of consistency,confirming the sub-homogeneous PD model’s efficacy in replicating the discontinuous deformation and failure of diverse coal-rock mass structures.

Acknowledgement:The authors would like to acknowledge the Key Laboratory of Geotechnical Mechanics and Offshore Underground Engineering of Qingdao for promoting this research and providing the laboratory facility and support,the National Natural Science Foundation of China and the Natural Science Foundation of Shandong Province for the financial support.The authors thank Shandong Energy Group for providing observation data for this project.

Funding Statement:The research described in this paper was financially supported by the National Natural Science Foundation of China (Nos.12302264,52104004,12072170,and 12202225),the Natural Science Foundation of Shandong Province (No.ZR2021QA042) and Special Fund for Taishan Scholar Project(No.Tsqn202211180).

Author Contributions:The authors confirm contribution to the paper as follows: study conception and design:Shijun Zhao,Qing Zhang;data collection:Weizhao Zhang;analysis and interpretation of results:Yusong Miao,Xinbo Zhao;draft manuscript preparation:Shijun Zhao,Wei Xu.All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials:The datasets used or analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.