A Novel Damage Model for Strata Layers and Coal Mass

Coal burst occurrences are affected by a range of mining and geological factors. Excessive slipping between the strata layers may release a considerable amount of strain energy, which can be destructive. A competent strata is also more vulnerable to riveting a large amount of strain energy. If the stored energy in the rigid roof reaches a certain level, it will be released suddenly which can create a serious dynamic reaction leading to coal burst incidents. In this paper, a new damage model based on the modified thermomechanical continuum constitutive model in coal mass and the contact layers between the rock and coal mass is proposed. The original continuum constitutive model was initially developed for the cemented granular materials. The application of the modified continuum constitutive model is the key aspect to understand the momentum energy between the coal–rock interactions. The transformed energy between the coal mass and different strata layers will be analytically demonstrated as a function of the rock/joint quality interaction conditions. The failure and post failure in the coal mass and coal–rock joint interaction will be classified by the coal mass crushing, coal–rock interaction damage and fragment reorganisation. The outcomes of this paper will help to forecast the possibility of the coal burst occurrence based on the interaction between the coal mass and the strata layers in a coal mine.


Introduction
There are critical locations in coal mines where coal burst can take place. A number of researchers have classified the main sources of coal burst based on the particular underground geometry of the mine as well as the geological conditions of the coal mine [1][2][3][4][5]. According to [1][2][3][4], the stiff roof, working face and roadways are these locations where there may be stress concentrations. Huang et al. [5] indicated that strata instabilites, including large roof deformation and coal burst, usually happens in the goaf-side entry of the longwall panel under a massive stiff roof. They discussed the energy storage capacity of overhanding competent roof strata for potential coal burst occurrences due to the sudden main roof breakage. The authors also introduced directional hydraulic fracturing to cut hanging hard roof to release the stored strain energy. The presence of a massive hard roof can cause coal fracturing when it reaches the weighting span, resulting in high dynamic stresses [6,7]. Bräuner [8] also reported that the majority of coal bursts occurred in the roof and ribs of coal mines. There is an extensive site investigation in the literature to determine critical regions for coal burst and found that the main cause of coal burst is related to the instability of the interaction between the coal and rock mass, individually, around the working face in the deep strata. Marcak [9] reported that a roof which deforms due to localised activities in an underground mine is a location prone to stress concentration as potential sources of coal burst. Further studiy [10] was carried out to evaluate the effect hard roof rupture or slip instability with the sudden release of large amount of elastic energy. Rock burst was divided into coal-roof type (impact pressure type) and roof type (impact type) that will be caused by violent shock easily. The excessive stiff roof displacement and the effect of the local stress changes on the interaction between the coal and stiff roof due to seismic activity were discussed by [11]. Periodic weighting during longwall mining can cause a cantilever beam reaction in a coal roof and may result in the release of a significant amount of the strain energy due to the hard suspended roof failure [12].
Apart from the competent roof strata, the presence of discontinunites has also been studied in the literature. Ryder [13] discussed critical regions in coal mines where there are significant geotechnical discontinuities such as joints, dykes, bedding planes and shear zones. Rock masses with discontinuities can create stress concentrations where there is a potential area for coal burst. Shear failure in the fault-slip due to the extraordinarily low friction is a main reason for the release of a significant amount of the kinetic energy. Also, when the rock mass is overstressed due to the effect of the massive overburden strata and the coal mass stresses reach the allocated uniaxial compressive strength (UCS), this condition may cause fracturing of the coal as a brittle material [14][15][16][17][18][19][20]. A number of studies [1,2,21,22] also indicated a coal burst may occur when the total stress which can be induced by a combination of the static and dynamic loading reaches a particular crucial stress level. There are a number of coal burst classifications proposed based on various approaches. Pan [23] indicated that coal burst can be classified based on the characteristic compression of the rock, the fracture of the roof and the fault potential activation type. Another approach is to classify coal burst based on the converting and accumulation of the strain energy. He [23] classified coal burst based on the single and compound conversion of the elastic energy to kinetic energy. There are two major stages in this coal burst occurrence. The first stage is static stresses reaching their highest point. Then, it is followed by the dynamic load distribution. Stress changes are also one of the main reasons for coal burst occurrences, induced by mining activities [24]. According to [17,[25][26][27] coal bursts can be classified as follows: strain burst, buckling, face pillar burst, shear rupture and fault-slip burst. In this classification, buckling is a particular type of strain burst, and shear rupturing is part of fault-slip burst. Coal burst can be induced by a remote seismic event which can cause damage in rock mass [27][28][29]. It is also indicated that rock mass failure usually takes place around excavation regions. The excavation region is a region where induced stresses, due to the effect of the great depth of the strata layers as well as mining activities, might exceed the ultimate strength of the rock mass. Fault-slip can also cause a significant coal burst near the excavation area. The most critical condition of a fault-slip is where there is a combination of increasing shear stresses and decreasing normal stresses. This combination of stresses can cause a strong coal burst, releasing a considerable amount of seismic energy [16,30]. Ground motion and vibration can cause significant damage in the rock mass near the excavation area. A high level of ground motion can be considered a dynamic trigger for coal burst [30][31][32].
Based on past studies in the available literature, the mechanisms of coal burst are still not yet fully understood. Although researchers have proposed several theories based on energy balance, system stiffness, or coal strength etc., the actual mechanisms are still being debated. So far, no universally acceptable method or criterion exists for forecasting a coal burst. Also, coal is a size dependent material, which means that by changing of the geometrical properties of a coal mass, the mechanical behaviour of a coal would be significantly affected. Therefore, a range of contributing parameters should be examined to have a sound understanding of the energy balance in a coal mass. This would be a key parameter to categorise the coal burst prone area. This study was carried out to address this knowledge gap. It should be noted that there are number of significant limitations to precisely forecast coal burst incident based on the energy methods due to the crucial sensitivity of the energy calculations on the value of uncertain materials input parameters.

Review of the Developed Damage Models in Brittle Materials
Edelen and Green [33], Edelen and Laws [34] classified generalized continuum mechanics theories into three categories. One category might be recognized by the hypothesis that more than three degree of freedom can be assigned to a particle; secondly, the higher gradients of the calculated displacement as well velocity can be concluded from the fundamental postulates while the third category includes those theories that indicate the likelihood of nonlocal interactions. Edelen and Laws [34] focused on a theory of the third category, focusing on the possibility of nonlocal interactions. Lemaitre [35] provided the background information on continuum damage mechanics in the framework of thermodynamics by providing some individual examples of the key equations for predicting the ductile, creep and fatigue damages. After discussing the key aspects to simulate and calculate the macro-crack initiation, by respecting to the non-coupled or coupled strain damage constitutive models, a number of numerical simulations were presented according to the fracture limits, surface damage and creep fatigue interaction [36]. The fundamental equations to simulate the coupling interaction between strain and damage behaviours, with respect to the framework of the thermodynamics of irreversible processes was developed by Lemaitre [37]. Sammis and King [38] developed a viscous regularisation of strain based, rate-independent damage model using a structure analogous to visco-plasticity, which produces retardation of micro-crack growth at higher strain rates. The suggested algorithms resulted in obtaining a simple and efficient stress update procedures. An anisotropic damage model for concrete by considering the general framework and the influence of the internal variable theory of thermodynamics was developed by Yazdani and Schreyer [39]. Bažant and Pijaudier-Cabot [40] indicated that distributed softening damage might be experimentally investigated on concrete specimens of a rectangular cross section. Ju [41] presented a novel coupled elastoplastic damage theory based on the energy approach concept. The proposed formulation was based on the irreversible thermodynamics and internal state variable theory for ductile and brittle materials. A coupled theory of elasticity and continuum damage mechanics was formulated by Kattan and Voyiadjis [42]. The suggested damage variable relied on the average material degradation which may reflect the different types of damage at the microscale level. Bažant [43] mentioned after imposing two homogenising conditions, it is shown that damage is a nonlocal variable that is a function of the averaged (global) strain from a certain neighbourhood of the given point. Ladeveze and LeDantec [44] proposed a model of the mechanical behaviour of a simple ply of a fibrous-composite laminate. Pijaudier-Cabot and Benallal [45] suggested that the local continuum acts as a lower bound of the global continuum while localisation, defined as a bifurcation into homogeneous modes of deformation, is concerned. The direction of the localisation band is the same. Hansen and Schreyer [46] developed a unified framework for coupled elastoplastic and damage theories. A rigorous thermodynamic process is pursued which is adequately universal to consist of anisotropic plasticity and anisotropic damage classifications. The approach of effective stress is the significant method for engaging these concepts. Yield and damage functions, built from homogeneous functions of degree one, are demonstrated to fulfill thermodynamic conditions [47][48][49].
Luccioni and Oller [50] presented a constitutive model that couples plasticity and damage. The suggested theory is thermodynamically dependable and appears from an overview of generic plasticity theory and isotropic damage assumption of Kachanov. The constitutive approach presented simultaneously solves the issue of progress of plastic strains and stiffness deprivation. softening and damage steady situations are both complied in each incremental computational times [51]. Collins and Houlsby [52] suggested the plasticity simulation of the brittle materials based on the thermodynamic outline. The framework had capacity to rigorously simulate both friction and non-associated flow, and the link between the suggested phenomena was illustrated. The suggested structure was focused on the evolution of the formulation models from theories about the form of the potential energy and a dissipation formulation [53,54]. Jirasek [55] computed nonlocal constitutive models used in simulations of damage and fracture processes of quasi-brittle materials. Jirasek [55] also presented a number of global analyses available in the literature based on the type of unpredictable parameters subjected to global averaging. Polizzotto and Borino [56] elaborated nonlocal and ascent plasticity in a unique procedure within an appropriate thermodynamic structure according to the global feature introduced. The procedure of the failure in the simulated material was governed by isotropic hardening/softening. Polizzotto and Borino [57] discussed a gradient-dependent associative plasticity theory which presents adequate thermodynamic consistency. Furthermore, the internal energy density was considered dependent on a nonlocal scalar internal variable apart from the external variables. An anisotropic damage evolution law based on the assumption of maximum dissipation was used to develop a constitutive model for plain concrete in the framework of multi-surface damage elasto-plasticity by Meschke, Lackner [58]. Lee and Fenves [59] developed a new plastic-damage model for concrete subjected to cyclic loading using the concepts of fracture-energy-based damage and stiffness degradation in continuum damage mechanics. Borino and Fuschi [60] attempted to launch elastic-plastic rate-independent material models with isotropic hardening/softening of global condition. Houlsby and Puzrin [61] discussed a procedure of an elastic-plastic formulation for rate independent materials which was developed according to the use of thermodynamic potentials. A number of energy functions usually used in thermodynamics such as inner energy, Helmholtz free energy, enthalpy and Gibbs free energy were defined to formulate the suggested model. Zhao [62] described that brittle material strength subjected to dynamic loads can be approximately simulated by the Mohr-Coulomb criterion, at low confining pressure range. Li [63] proposed two energy equivalence principles based on the geometric definitions of the damage restriction and thermodynamic standards, which could provide relationships between the macroscopic and valuable descriptions of stress and strain. Comi [64] proposed a global isotropic damage approach which can be used to forecast the behaviour of brittle materials up to failure. Picandet and Khelidj [65] established a relationship between the amplify in permeability and the applied-strain/yield-strain ratio. A thermodynamically reliable global procedure for damaging materials was also suggested by Benvenuti and Borino [66]. Jirásek and Patzák [67] introduced the simulation method of strain localisation simulations based on the global continuum damage models of the integral type. Borino and Failla [68] presented a thermodynamically inclusive formulation for global damage models. Global models were recognized as a theoretically efficient and computationally competent approach to rise above the shortcomings arising in continuum media with softening. Jirásek and Rolshoven [69] devoted special attention to the thermodynamic aspects of global plasticity, individually to the consistent development of the concept of general standard materials to global continua. Simone and Wells [70] described the combined continuous-discontinuous failure in a regularized strain-softening continuum. The continuum was classified through the introduction of gradient terms into the constitutive equations. Combined models allow for the analysis of the entire failure process, from diffuse micro-cracking to localised macro-cracks. Rodrıguez-Ferran, Morata [71] presented a well-organised and consistent view for the numerical modelling of failure with nonlocal damage models. Salari and Saeb [72] suggested a triaxial formulation for elastoplastic behaviour of the rocks, which includes for tensile damage. They formulated the constitutive in the framework based on the continuum thermodynamics using internal variables. Peerlings and Massart [73] derived a new gradient damage model extracted from thermodynamic statements. Nonstandard terms were added to the free energy potential for this purpose. A trend to suggest an outline for global damage-plastic models, combining local plasticity with integral-type nonlocal damage, was proposed by Grassl and Jirásek [74]. The main derived damage evolution was governed by plastic strain. Makowski, Stumpf [75] suggested the equilibrium conditions of material forces together with the classical equilibrium conditions of physical forces. Wang and Li [76] described a method from damage to fracture with a continuous procedure. The description of the continuous damage procedure, defined as a damage-softening probabilistic constitutive model for the brittle material, was suggested based on the Weibull distribution of mesoscopic element strength. Comi, Mariani [77] proposed an integrated strategy to numerically model damage evolution and crack transmission in the brittle rocks. At the first stage of the crack propagation, the simulated brittle material was taken into account to be macroscopically integer and was modelled by a non-symmetric isotropic global damage model. Choinska and Khelidj [78] investigated damage-temperature-stress level-permeability interactions in the brittle concrete structures. It was indicated that the same trend indicated by Picandet et al. [65] where expanding the interaction between permeability and damage permeability amplifies for a normal and a high-performance concrete for a normal compressive damage, in the before failure stage. Nguyen and Houlsby [79] suggested an evolution of a coupled damage-plasticity constitutive formulation for brittle material. Voyiadjis and Kattan [80] discussed and compared various definitions of the damage variables. Individually, special attention was devoted to a new damage variable which was determined in terms of the elastic stiffness of the material. Both the scalar and tensorial cases were researched. Zhou and Zhu [81] derived an elasto-plastic damage formulation model with double yield surfaces based on permanent thermodynamics speculation and damage mechanics approach. Feng and Yu [82] described the methods for characterising progressive damage of micro-cracks and introduced the concept of direction field of micro-crack expansion. Ferrero et al. [83] described mathematical assessment of the rock-joint damage evaluated by a photogrammetric review of the discontinuity while the cyclic rock test was undertaken. Duriez et al. [84] classified a continuously nonlinear formulation interaction to express the mechanical behaviour of in-filled rock discontinuities. The suggested formulation by Duriez et al. [84] was calibrated by taking advantage of a discrete element model, which was compared with experimental data. Li et al. [85] introduced a fundamental theory of damage mechanics to be concerned with the deformation of strain softening for rocks. A new damage variable D was determined by Liu et al. [86] which could describe the degradation process of sandstone samples upon dynamic cyclic loading. Song et al. [87] studied the damage evolution process and crack development of rocks under cyclic uniaxial compression. Song et al. [88] reported that digital image correlation is a useful technique to examine the advancement of displacement and strain fields in the different brittle materials. It can explain the entire procedure of damage accretion, crack opening and propagation of the brittle rock subjected to the various applied loadings. Subsequently, the damage factor was computed according to the deviation of apparent strain, as well as the localisation factor which was determined to explain the level of deformation localisation [89]. Pourhosseini and Shabanimashcool [90] explained a constitutive model to describe the nonlinear behaviour of intact rocks under static loading. Liu et al. [91] developed a new damage constitutive model was based on dissipated energy to illustrate the reaction of the brittle materials under cyclic loading. Vu, Mir [92] proposed a thermodynamically inclusive novel formulation for establishing constitutive models for engineering materials.

Analytical Approach
This section proposes a new damage model based on the modified thermomechanical continuum constitutive model in coal mass and the contact layers between the rock and coal mass. The original continuum constitutive model developed by [93] for cemented granular materials discussed the micro-mechanical interaction between the cemented and non-cemented grains. If coal mass can be assumed to perform same as the grains in a numerical model, and joint/cleat are same as the cemented materials between the grains, a new damage model can be implemented inside the modified thermomechanical continuum constitutive model for brittle granular materials, which is a function of the level of coal joint/cleat density. In this model, the internal variable Breakage, B, can be introduced as the area ratio. Bp is the breakage potential, that is the total area included between the initial and the final distribution, and Bt is the area within the initial and current distribution, meaning the amount of crushing that the grains have undergone [94,95]. By integrating the probability density function pdf f g (x g , B) along the coal mass size x g from the smallest dimension x g m to the largest x g M , the cumulative probability density function cdpf is F g (x g , B). Thus, according to [94,95]: Figure 1 demonstrates the proposed damage function range in a coal mass (B), which is a function of the joint/cleat density per unit of metre, D. D1 to D20 are presenting different damage function based on the different joint/cleat densities. The same concept which was explained in the B can be applied for the D.  Figure 2 illustrates the damage model in the cumulative probability damage model in the coal-rock interface, which is the function of the coal mass damage model as given in Equation (5), and P1 to P20 are presenting the different cumulative probability of the coal-rock interface due to different possible damage intervals.
where the statistically homogenised variable can be evaluated through Equation (8). In this individual case, it was assumed that: According to the Helmholtz free energy in brittle granular material, the internal stored energy of a coal mass as a brittle material in the granular phase can be presented by: where Ψ g and Ψ c are the respective specific energies, which can be calculated as under the curve of the stress strain of the uniaxial test under crushing conditions. Also, ϕ g and ϕ c are composite volume fractions, which are dependent on the coal particle size distributions in the rheology of coal-water slurries. The percentage composed of coal particles with different sizes can be used to determine the effect of volume fraction on the viscosity of the slurry. The density of the Helmholtz free energy of a coal mass size x g , Ψ g is: where, the notation i for initial and u for ultimate is presented. By considering an energy is expected to scale with the coal mass size x g in order of the coal size can be demonstrated by: where ξ g Ψ is the energy split function which distributes the stored energy among coal mass grains according to coal grain sizes x g , Ψ g r is the energy stored in a reference coal mass grain size, x g r is a function of the elastic strain ε e ij while φ denotes the porosity of the coal. The total elastic strain of the coal mass phase can be expressed by: a total−strain = a 2 0 + a 2 1 + a 2 2 (15) According to Einav [94,95], the energy split function is a function of the non-dimensional power in the form: Individually, for a coal mass the energy split would be a function of the joint and cleat density as well as the joint dip direction. According to the proportion of the coal grain sizes x g and reference coal mass grain size x g r is a non-dimensional function. According to Einav [94] and Einav [95]: where ϑ g is the grading index in a coal mass: When a reference coal mass size is expediently considered as: where V g i (X g ) and V g u (X g ) are the initial and ultimate volume of the coal mass grain sizes. To calculate the transmitted energy inside the rock/roof strata, the Helmholtz free energy stored in the coal-rock interface can be expressed by Einav [94,95] and Tengattini, Das [93]: Based on the common elements in the coal-rock interaction sizes, the energy rating can be demonstrated by: The momentum energy between the coal-rock interactions is a function of the rock/joint quality interaction conditions. Thus, where RQD is Rock quality designation according to Bieniawski [96]: J n is Joint set number (number of discontinuity sets), J r is Joint roughness number for critically oriented joint set (roughness of discontinuity surfaces) J a is Joint alteration number for critically oriented joint set (degree of alteration or weathering and filling of discontinuity surfaces), J w is Joint water reduction number (pressure and inflow rates of water within discontinuities), SRF is stress reduction factor (presence of shear zones, stress concentrations, squeezing or swelling rocks). According to Tengattini, Das [93], it can be assumed that the energy split function can be described by a non-dimensional power law: it can be assumed that the energy split function can be described by a non-dimensional power law: The strength of a joint or block surface is quantified and represented by a factor called Joint Condition Factor (J C ) following Barton and Bandis [97] and Palmstrøm [98], and is defined as: where J w is the joint waviness, J s is the joint smoothness, and J a represents joint alteration. Thus, where, for simplicity and accordingly, the reference coal-rock interaction area X c r is For persistent joints, block volume V 0 is given by the equation: where s i = spacing between joints in each set; γ i = angles between the joint sets [27]. Therefore, Analogous to the breakage variable, a grading index is introduced for the damage: Thus, evaluating the Helmholtz free energy for the joint can be presented by Equation (31) [93] and Equation (32): A convenient form to represent the two laws of thermodynamics of rate independent materials in isothermal conditions is [94]: where • Ψ is the rate of change of the Helmholtz free energy, Ψ, Φ is the increment of energy dissipation and W is the increment of mechanical work done on coal or the coal-rock interface.
where σ ij and ε ij are the stress and strain rate tensors applied to the coal mass as well as the coal-rock interface.
The inelastic conditions in the coal and coal-rock joint interaction are coal mass crushing, coal-rock interaction damage and fragment reorganisation. The rate of dissipation can therefore be considered as the sum of these three components: As it was indicated, the relationship between the damage in the coal mass and different strata layers as well as different key energy factors are presented by Equation (53). As it is demonstrated, there is a non-linear relationship between the damage parameters and the key energy parameters. Further explanation will be provided in the following section.

Discussion of the Results
The proportion of damage D/B between the coal mass and the damage in the strata layers in the proposed equations can play a significant role, which is evident in Equation (53). In order to illustrate the relevant interaction between the strain, kinetic and dissipated energy, different values were assigned to the damage ratio D/B between the coal mass and the strata layers as shown in Figures 3-5. Previous literature demonstrate the interaction between the different types of energy due to the effect of the damage in a coal mass and strata layers, where the D/B ratio varies between 2 n and 10 n . The D/B can also determine whether the strain energy would be mostly converted to kinetic or dissipated energy. This is a notable factor to determine in forecasting coal burst prone regions in a coal mine. The proposed parameter can also determine the interaction between the stiffness of the roof and the stiffness of the coal mass to determine sudden energy absorption or releasing between the coal mass and the roof. Thus, D/B ratio can be a significant factor when it comes to forecasting a local brittle failure in the rock/coal interaction and the coal masses. Since the D/B ratio can be directly relevant to the stiffness ratio between the coal mass and the strata layer, it could be taken into account as a criterion to evaluate the effect of the high stiff roof on the converted kinetic energy due to the sudden brittle failure. In general, the interaction between the stiffness of the coal layers and the strata can play a key role to determine the failure mode. In the current study a novel solution has been proposed to determine a new damage model in both coal mass and the strata layers based on the energy approach concept.

Conclusions
Based on the proposed energy based equations, the proportion of damage D/B between the coal mass and the damage in the strata layers (roof layer) has been developed in this paper. It is defined as the function of the interaction between the strain, dissipated and kinetic energies. It was assumed that the D/B ratio varies between 2 n and 10 n . Since the D/B ratio can also present the stiffness ratio between the strata layer (roof layer) and the coal mass, it can also determine the interaction between the high/flexible stiffness of the roof and the stiffness of the coal mass to determine the possibility of the local brittle failure in the coal/rock interface, which is a critical consideration in occurrence of coal bursts.. The energy conversion from the strain energy to the kinetic and dissipated energy was smoothly taken place when the D/B ratio was decreasing from 10 n to 2 n This clearly demonstrated that the stiffer strata will have a more brittle failure compare to a weaker and less stiff unit. In a competent strata, most of the strain energy is suddenly converted to the kinetic energy due to the strata breakage; on the other hand, in a weaker strata, the strain energy tends to be converted to the dissipated energy.