Rock Failure Analysis under Dynamic Loading Based on a Micromechanical Damage Model

A micromechanical constitutive damage model accounting for micro-crack interactions was developed for brittle failure of rock materials under compressive dynamic loading. The proposed model incorporates pre-existing flaws and microcracks that have same size with specific orientation. Frictional sliding on micro-cracks leading to inelastic deformation is very influential mechanism resulting in damage occurrence due to nucleation of wing-crack from both sides of pre-existing micro-cracks. Several homogenization schemes including dilute, Mori-Tanaka, self-consistence, Ponte-Castandea & Willis are usually implemented for up-scaling of micro-cracks interactions. In this study the Self-Consistent homogenization Scheme (SCS) was used in the developed damage model in which each micro-crack inside the elliptical inclusion surrounded by homogenized matrix experiences a stress field different from that acts on isolated cracks. Therefore, the difference between global stresses acting on rock material and local stresses experienced by micro-crack inside inclusion yields stress intensity factor (SIF) at the cracks tips which are utilized in the formulation of the dynamic crack growth criterion. Also the damage parameter was defined in term of crack density parameter. The developed model was programmed and used as a separate and new constitutive model in the commercial finite difference software (FLAC). The dynamic uniaxial compressive strength test of a brittle rock was simulated numerically and the simulated stress-strain curves under different imposed strain rates were compared each other. The analysis results show a very good strain rate dependency especially in peak and post-elastic region. The proposed model predicts a macroscopic stress-strain relation and a peak stress (compressive strength) with an associated transition strain rate beyond which the compressive strength of the material becomes highly strain rate sensitive. Also the damage growth process was studied by using the proposed micromechanical damage model and scale law was plotted to distinguish the dynamic and quasi-dynamic loading boundary. Results also show that as the applied strain rate increases, the simulated peak strength increases and the damage evolution becomes slower with strain increment.


Introduction
Brittle materials such as rocks show highly non-linear and complex response to external dynamic loads especially in peak and post-peak region.The comprehension and interpretation of the rock failure mechanism under dynamic loading is important for engineering practices.The measures to promote the failure of rocks in open pit mines and underground spaces rely on the fundamental understanding of the failure mechanisms in rock materials.The degree of issue greatly increases when the rock material is under dynamical loading.The rock materials have intrinsic defects in micro-scale such as in-homogeneities, micro-pores, micro-cracks and grain boundaries mismatches that cause stress concentration and decrease the overall rock strengths [1].These intrinsic defects could be categorized into different families in term of size and orientation.According to the weakest link theory under quasi-static loading condition some of the large flaws in the material are activated that affect significantly overall behavior of material.In other words, small flaws within the material do not activate crack growth, and failure is therefore controlled by large size flaws.In contrast with static loading, under dynamic loading the weakest link theory is not applicable and the whole families of the flaws participate in the mechanical response of brittle solids [2].
To describe damage effects on rock behavior, different kinds of constitutive damage models for the rock materials have been developed for various applications and these models can be categorized into two types.The first type is called phenomenological damage model, and the other type is the micro-mechanical damage model.In the framework of phenomenological models, the material free energy is formulated as a function of a number of internal variables, such as plastic strain and damage variable.
Actual physical mechanisms e.g.unilateral effects due to opening and closure of micro-cracks, coupling between frictional sliding on the micro-crack surfaces and damage evolution, as well as interactions between cracks that have fundamental role in materials overall behavior and damage due to wing-crack nucleation are not considered in this kind of model [3].Therefore, phenomenological damage models might be efficient for computational purposes, but they are generally unable to reproduce the microscale mechanisms.Furthermore, large number of parameters has usually to be determined in phenomenological models, and many of them have no actual physical interpretation.On the other hand, micromechanical damage models have been developed based on homogenization and up-scaling process.This kind of models can describe the main mechanical features of brittle materials to a certain extent and they are generally based on phenomena and mechanisms in microscale rather than on a solid mathematical derivation from macroscale phenomena.Micro-mechanical damage models have ability to consider the physical mechanisms involved in damage process and deterioration of material properties calculated by damage variable representing micro-cracks statistics features [4].
The strength of brittle materials depends largely on the applied loading rate for both tensile and compressive loading.The mechanical behavior of brittle material under dynamic loading can be studied at the experimental scale by the Kolsky bar (or similarly split-Hopkinson bar) test.Under high loading rate, propagation of cracks with large sizes is rate limited, so the stress magnitude in brittle material continues to rise to a high level.When the stress intensity factor of loading reaches the critical value (toughness), wing-cracks also nucleate from smaller flaws.Due to the preferential direction of the wing-crack nucleation from initial flaw tips, anisotropy and volume dilatancy are induced in damaged material [2].As a result, higher absorbed and dissipated energy through the multiple opening wing-crack surfaces leads to the dynamic strength of rock exceeded its static strength which proved by experimental observation [5].
An increase in number and overall length of cracks and a decrease in resulted fragment sizes due to high applied loading rate are widely reported in the experimental studies, which supports this hypothesis [6].Due to lack of plastic deformation related to slipping on closed micro-crack faces, the typical failure mode of brittle rocks is associated with micro-cracks growth and their coalescence.In order to solve an engineering problem under dynamic loading, a micromechanical constitutive model reflecting the actual micro-scale defects behavior and dynamic micro-fracture mechanisms should be proposed.Several researchers studied on materials response to dynamical loading.For instance, Nemat-Nasser and Deng (1994) calculated the mode I stress intensity factor (SIF) for a pre-existing sliding micro-crack under dynamical loading, and they used dilute scheme for micro-cracks homogenization and up-scaling [7].Initial micromechanical damage models are usually limited to dilute distribution of micro-cracks and then are not able to take into account interaction between micro-cracks [8].In order to overcome the shortcomings of dilute distribution, other Eshelby-based homogenization procedures such as Mori-Tanaka (1973), Ponte-Castaneda and Willis (1995) and selfconsistent were proposed later [9].Paliwal and Ramesh (2008) considered micro-cracks interaction by using the self-consistent homogenization scheme to calculate the effective properties for the dynamic behavior of brittle materials under biaxial compression loads [9].Katcoff and Graham-Brady (2014) implemented the Self-consistent micro-mechanical damage model to study the compressive dynamic failure of brittle materials with circular flaws [10].Tonge et al. (2013) applied a similar methodology as part of a three dimensional macro-scale constitutive model of brittle dynamic failure [1].Tonge et al. (2016) focused on the model formulation that captures the strength variability and strain rate sensitivity of brittle materials and presents a statistical approach to assigning the local flaw distribution within a specimen.Their works developed a model for brittle material failure during impact events.They used SCS to homogenization of cracked media and complex method to determine local stress from remote stress acting on specimen [11].
Although Hu et al. (2015) proposed a plastic-damage model based on the self-consistent homogenization scheme, matrix plastic strain due to dislocation in crystalline network is only considered in this model and material degradation and inelastic deformation are more popularly characterized separately [12].The effective homogenized properties of cracked materials are obtained by following an up-scaling method based on the Eshelby inhomogeneous inclusion solution [13].Ayyagari et al.(2018) proposed a three-dimensional generalized anisotropic constitutive model for brittle solids contain of the spatial evolution of planar wing-cracks subjected to dynamic compressive loading [14].
In this paper, according to Figure1 the micro-crack interactions were considered through a crack-matrix-effectivemedium approach based on the self-consistent homogenization scheme to model the brittle failure process of rock.In self-consistent scheme it is assumed the presence of an elliptical inclusion surrounding each individual flaw inside of which the material is undamaged.Material outside of the ellipse has the effective homogenized properties accounting for the damage associated with the entire flaw population.The interactions among different flaw families are accounted through the means of crack-matrix-effective-medium approach, which the self-consistent scheme plays a crucial role in determining the local (effective) stress field around the individual flaws.
Under uniaxial compressive loading, the mode I SIF for a pre-existing sliding micro-crack was calculated, and then dynamic effect on mode I SIF was taken into account.The key factor to study the micro-cracks propagation under dynamic loading is dynamic stress intensity factor in modeling process.In this problem the dynamic stress intensity factor only depends on the current crack tip speed and static stress intensity factor.According to Freund (1972Freund ( , 1990)), the mode-I dynamic stress intensity factor can be determined for propagating flaws by multiplication of a universal function of the propagation rate of micro-cracks tips in static mode I SIF.[15].
Regarding the SCM approach, the constitutive model related to failure of brittle materials can be decomposed into two parts (a) Homogenized mechanical properties (compliance or stiffness tensor) with respect to damage, which is used to identify the local stress state on elliptical inclusion boundary.
(b) Propagation of micro-cracks under subjected loads, including the crack growth after wing-crack nucleation, the interactions among the pre-existing micro-cracks and coalescence to create the macro-scale failure plane.
The aim of this study is to develop the micro-mechanical damage model originally proposed by Paliwal and Ramesh (2008) to take into account the frictional sliding and damage evolution due to micro-cracking under dynamic compressive loading.The proposed model is formulated in the framework of a micromechanical model based on the self-consistent homogenization scheme which was programmed and implemented into a commercial code.For simplicity, the special case of a material whose all flaws have uniform distribution of size and orientation in twodimensional space is assumed throughout this paper.In this paper, the method introduced by Graham-Brady et al. ( 2015) was implemented to calculate the local stress field imposed on the elliptical inclusion containing an isolated micro-crack as a function of the global principal stress tensor.

Theory and Background
Up scaling analyses are usually performed over a representative elementary volume (REV) (area A for the 2-D case) of an elastic brittle material containing distributed pre-existing flaws.In micromechanical damage models, the REV behavior is obtained by homogenization, after calculating local stresses and displacements at crack faces [16,17].Because of the surrounding cracks, each activated flaw experiences a stress field that is inherently different from the stress field of an isolated crack and also different from the external stress field.This effect can be modeled in a selfconsistent scheme, by replacing the material surrounding each flaw (which is enclosed in the ellipsoidal inclusion in Figure 1) with an effective material having the reduced elastic modulus:  ̅ = .[1 − (̅ ).Ω],  ̅ = .[1 − (̅ ).Ω] that is a function of the evolving damage parameter  = ∫ () 2 .This damage parameter accounts for the damage caused by wing cracks of length () associated with each individual flaw length (2) as described by the PDF (()) and the flaw density ().() and (̅ ) functions were described as follows [9]: Figure 1.The schematic illustrations of SCM scheme As regarded, the properties of material inside elliptical inclusion differ from the effective properties of materials around inclusion calculated by SCM scheme, so the local stress field in each micro-crack elliptical inclusion boundary is different from the overall remote stress.This issue leads to sliding on pre-existing micro-crack associated with wingcrack nucleation.
As shown in Figure 2, a micro-crack is assumed with two wing-cracks that propagate from both tips of crack embedded in an elliptical inclusion of pristine isotropic material, which in turn is located in a medium with the effective properties of the micro-cracked solid matrix.The difference in the properties of the elliptical inclusion and the homogenized cracked matrix leads to the crack inside the inclusion responds to a stress field different from the applied macroscopic stress field.Therefore the resulting mismatch of material properties at the elliptical inclusion boundary leads to local stresses fields (  ) within the elliptical inclusion thatdiffers from the far-field global stress ().
One of the most important challenges for employing the self-consistent scheme is determination of local stress field acting on the elliptical inclusion containing micro-cracks.In this paper, the method introduced by Graham-Brady et al. (2015) was implemented in this paper to calculate the local stress field [13].It is noteworthy that another method proposed early by Paliwall and Ramesh (2008), complex method, can be implemented to determine the local stress fields around isolated micro-cracks [6].To calculate the local stress field acting on the elliptical inclusion containing the individual micro-crack, the method proposed by Graham-Brady et al.( 2015) is used as following [13]: Where [] is the 3 × 3 identity tensor, [  ] and [  ] are undamaged stiffnesstensor for material in the ellipsoidal inclusion and stiffnesstensor for the damaged matrix material respectively and they are determined as follows [13]: Where,  is the material Young's modulus and[C * ] is determined as follows [13]: Where [] is the Eshelby tensor for the problem of an ellipsoidal inclusion in a matrix, defined as [13]:  As illustrated in Figure 2, the major axis of the ellipsoidal inclusion is directed parallel to the global applied loading.The other two axes are perpendicular to the major axis.The nine components of the Eshelby tensor can be calculated by the following equation [13]: Where, / is the aspect ratio of the ellipsoid inclusion.The solving of numerical integration proposed by Graham-Brady et al. (2015), by supposing; (/) 2 = , leads to the following results for  1 ,  2 , and  3 [13]: As illustrated in Figure 3, the resolved shear stress on the pre-existing micro-flaw due to the slip-induced gap, b, (mostly called wedging effect at the pre-existing micro-crack) results in a significant increase in the stress intensity factor at the micro-crack tips which is utilized in the formulation of the crack growth dynamics.Equation 4, provides a straightforward relationship between the global principal stresses and the local stresses on each flaw embedded in the elliptical inclusion, which are used to calculate the stress intensity factor at the crack tip associated with each flaw.The general formulation for calculating the mode-I SIFis presented [9,13]: Where,   is the micro-cracks planes cohesion,  is initial orientation of the micro-cracks, 2 is initial size of the microcracks and  is the wing-crack length after micro-crack propagation.The parameter ( * = 0.27) is assumed to avoid singularity of (  ) when the wing-cracks are too small ( = 0).
When the rock materials are subjected to continuous compressive loading, local stress concentration (stress intensity factor) (SIF), at the crack tips increases.At the specific moment, the SIF at crack tips reaches the rock material fracture toughness level, so the tensile wing-cracks with length of  propagate from the tips of a pre-existing micro-cracks [18].This phenomenon is called wing-crack nucleation as shown in Figure 3.After nucleation of tensile wing-cracks, the interaction among the tensile cracks increases as these cracks grow close to each other.This interaction makes local stress concentration at the crack tip and increases the mode-I SIF (  ) that leads to the crack opening displacements and hence affects the overall response of the solid.The crack growth processes leading to damage are profoundly related to such pre-existing flaws at the micro-scale and decrease the material stiffness [1].

Figure 3. Local stress field acting on the isolated micro-crack and wing-cracks nucleation [5]
Here it must be noted that the SIF introduced in Equation 10 is in the static form.However, the dynamic form of SIF is essential to study the micro-cracks propagation under dynamic loading.The Mode-I dynamic SIF for two-dimensional in-plane crack growth has been established by Freund (1972).Under dynamic loading, a micro-crack could propagate with non-uniform speed.According to Freund (1972), the dynamic SIF can be estimated by multiplying the equivalent static SIF (calculated for the identical loads but assuming zero crack tip speed) by a function of the crack tip speed.The function k( ̇) is a universal function of the crack speed, which represents the inertial effect on the crack growth.The function k( ̇) can be written as [15,19]: Where   is the Rayleigh wave speed, and  ̇ is the wing-crack tips propagation velocity.After calculating the function ( ̇), the dynamic SIF is defined in terms of the static SIF as follows [20]: Wing cracks are initiated from the pre-existing flaws when the mode-I SIF at the crack tip reaches the mode-I fracture initiation toughness (FIT) of the material.In this work, the mode-I FIT was assumed to be rate independent and equal to the mode-I fracture toughness   : Combination of the Equations 12 to 14 results in wing-crack propagation velocity expressed as follows [15]: Where   shows the maximum (terminal) speed of a dynamic propagating crack and (, ) are the fitting parameters characterizing the toughness-velocity relation.

Computational Algorithm
Figure 4 shows the computational processes carried out in one step to calculate the wing-crack length, damage parameter and updated stress tensor in developed micromechanical damage model.The strong interaction between crack growth (damage) and micro-flaw nucleation as well as the strong nonlinearity in mechanical response is considered in this proposed model.In this algorithm, a trial elastic stress prediction and then damage correction are performed.Assuming at the end of step , the internal variables and the macroscopic stress and strain tensors have been determined with a prescribed precision after an iterative process.With assuming a new strain increment ∆  (+1) , the calculation done in step ( + 1) is summarized as follows:  Trial elastic prediction: At first it is assumed that the strain increment ∆  (+1) is completely elastic and contributed by the deformation in the matrix phase.Thus,   (+1) =   () + ∆  (+1) ,   (+1) = σ ij () + E  () ∆  (+1) .
 Damage correction: Applying continuous loading over time leads to the condition that SIF is greater than   i.e.   ≥   .Therefore, the wing-cracks propagate over time and the damage evolution can be calculated.Then, the stress tensor is corrected based on the new damage parameter   (+1)  = σ ij () + ∆E  ()   (+1)  Residual behaviour prediction: when the damage parameter in rock material reaches a certain critical value, the stress-strain curve become a roughly horizontal path (residual strength).

Numerical Simulation
To implement the developed micromechanical damage model in numerical simulation of engineering problems, it was programmed within the Fish environment provided in the commercial finite difference software FLAC as a separate and new constitutive model according to the algorithm illustrated in Figure 4.
order to investigate the proposed model performance, the dynamic uniaxial compressive strength test of a brittle rock was simulated numerically.Therefore, it was attempted to simulate numerically the uniaxial compressive strength test condition as closely as possible.The simulated sample geometry and boundary condition were selected similar to the standard test condition reported in ISRM suggestion as shown in Figure 5.The main objective was to reproduce numerically the strain rate sensitivity of the rock stress-strain curve and delve into the sample failure mechanism in the peak and post-peak regions.The dynamic loading was simulated by imposing a velocity field at the top boundary in the y-direction equivalent to the dynamic strain rates including (0.16, 0.24, 0.28  0.32) × 10 5  −1 while a roller Trial stress calculation (with an elastic stress assumption):   ().=   (−1) +  (−1) ∆  () Based on the far-field stress acting on REV homogenized by SCM, calculate the local stress field around inclusion.From the local stress field, the resolved shear stress and SIF could be calculated, from Eqs ( 10) and (11).
Update trial stress: boundary condition is placed at the bottom boundary.There are no constraints on the sides of the simulated sample and the specimen sides are allowed to move in the horizontal and vertical directions.

Figure 5. Geometry and boundary condition of the simulated specimen
The pre-existing flaws characteristic parameters and mechanical properties for a typical rock material used in the numerical simulation are listed in Table 1.For simplicity, the flaws in our study were presumed to have a uniform distribution for flaws size and orientation.Therefore,  = 50.7 ° is selected for orientation because it is the most critical orientation under uniaxial compressive loading.It's noteworthy that the fitting parameters (, ) are selected as which is adopted from Paliwal and Ramesh (2008).
The main objective of the analysis presented in this section was to demonstrate the capability of the proposed model in simulating the rock strength dependency on loading rate.In order to evaluate the general behavior of numerical sample, stress and vertical displacement variables on the upper model boundary were recorded and averaged under compressive loading.
The stress-strain curves are simulated under a range of the imposed strain rates from ̇= 1.6 × 10 4 1  ⁄ to ̇= 3.2 × 10 4 1  ⁄ to investigate the variation of the compressive peak strength by the imposed strain rates.The simulated stress-strain curves under various applied strain rates shown in Figure 6 demonstrates that the applied strain rate affects significantly on the rock material strength and the simulated peak strength is greatly dependent on the applied loading (strain) rate.

Figure 6. The stress-strain curves simulated under various imposed strain rates
Even though the simulated stress-strain curves show linear elastic behavior at low strain values, the flaws begin to be activated by stable crack growth which leads to increase in damage of material resulting in stiffness degradation and decrease in slope of the stress-strain curves.Once the stress state in rock material reaches the peak strength, the slope of the simulated stress-strain curve becomes negative in post-peak region as a result of unstable crack growth.Furthermore, the increment of imposed strain rate with the same pre-existing flaws parameters leads to increase in the simulated peak strength.Therefore, applied strain rate dependency can be clearly observed in the simulated stress-stain curves under dynamic loading condition.The fact that flaws have not enough time to grow and propagate under applied high strain rates explains why the peak strength of material increases under dynamic loading by increment of the imposed strain rate.At high applied strain rate, the crack propagation is overcome by very high rate of loading; therefore, the rock material is able to carry higher stresses before ultimate failure and the strength of the material increases under high dynamic strain rates.As the applied strain rate becomes lower, the damage growth process with strain increment becomes faster.
The area below the post-peak part of the stress-strain diagram represents the energy dissipated by the rock specimen.According to Figure 6, for higher applied strain rates the energy dissipation by material in post-peak region during failure process increases.The residual parts of all the simulated stress-strain curves under various applied strain rates can be seen in Figure 6.While the accumulated damage parameter in rock material reaches the critical value predefined in algorithm, the developed micromechanical damage model reflects the residual behaviour of rock material.
Here, It is worth noting that the dynamic Young's modulus of a rock is often greater than the static one, but to a small extent.For the rocks tested, their dynamic Young's modulus are greater than their static ones by less than 30%, indicating that the difference between static and dynamic Young's moduli is very small, compared with the difference between static and dynamic rock strengths [5].That's why the slopes of curves in Figure 6 are approximately same.
The developed micromechanical model can adequately reproduce many features of the rock behavior such as the linear elastic, hardening prior to the peak strength and eventually softening in post-peak region.The simulated compressive strength is dramatically sensitive to the imposed strain rate especially in dynamic strain rate ranges.This is in agreement with a great number of experiments indicating that compressive rock strength is sensitive to the applied loading rates, particularly in dynamic loading range [18].
The simulated compressive peak strengths are plotted against the imposed strain rates in logarithmic scale to investigate the strain rate-sensitivity of the compressive peak strength predicted based on the micromechanical damage model proposed in this paper as shown in Figure 7.According to Figure 7, variation of the applied strain rate affects the mechanical behavior of rock materials.The simulated compressive peak strength slightly increases with rise in loading rate under quasi-static loading condition, but it dramatically increases with increment of loading rate under dynamic loading condition.While the imposed strain rate exceeds the transitional strain rate representing the boundary between quasi-static and dynamic loading conditions, the compressive peak strength increases significantly.
Since the time variable used in dynamic numerical calculations is real, the damage growth process within the simulated sample can be monitored and recorded over time domain.To do this purpose, under the applied strain rate of (ε̇= 0.32 × 10 5 s −1 ), the damage contours in the simulated rock sample were monitored and plotted at different times and steps.The damage evolution process within the simulated sample is shown in Figure 8 based on the proposed micromechanical damage model.According to Figure 9 in the first phase OA, the rock behavior remains linear elastic, the stress and strain are linearly related to each other and no damage occurs before the activation of sliding along the pre-existing flaws.At this stage, the damage in material is purely from the pre-existing internal flaws and very small.The slope of the simulated axial stress-strain curve is the same as the Young's modulus.In phase AB, the progressive accumulation of the frictional sliding leads to the stress concentration at the pre-existing flaws tips and increase in the stress intensity factor at these flaws.Once the dynamic stress intensity factor reaches the fracture toughness of the material, in the third phase BC the pre-existing micro-cracks nucleate and develop into wing-cracks and the damage evolution law is activated.

Conclusion
In this paper, a micromechanical damage model for brittle rock materials was developed under dynamic compressive loading condition based on frictional sliding along the pre-existing micro-cracks and associated wing cracks sprouting from the tips of the pre-existing flaws.These pre-existing micro-cracks are assumed to be uniformly directed in rock material and all of them have same sizes.Frictional sliding along pre-existing flaws leads to nucleation of wing cracks that causes damage growth in rock material.The proposed micromechanical damage model incorporates the interactions between micro-cracks by means of the self-consistent homogenization scheme that calculate an effective stress field around the micro-cracks which is different from that acts on isolated micro-cracks.The developed micromechanical damage model was programmed and used as a separate and new constitutive model in the commercial finite difference software (FLAC) under different applied strain rates.Variation of the applied strain rate significantly affects the mechanical behavior of brittle materials.The developed micromechanical damage model is able to simulate the complete stress-strain curve under dynamic loading and an associated transitional strain rate beyond which the compressive strength of the rock material becomes highly strain rate sensitive.The simulated compressive peak strength remains nearly constant and insensitive to the imposed strain rates below the transitional strain rate representing the boundary between quasi-static and dynamic loading conditions.While the imposed strain rate exceeds the transitional strain rate, the compressive peak strength increases dramatically by increment of the applied strain rate.Furthermore, the damage evolution in the simulated rock specimen was monitored and recorded under dynamic uniaxial compressive loading condition.The accumulated damage within the simulated rock specimen at low stress level is insignificant, but damage spreads in the simulated rock material over time in peak and post-peak region.

Conflicts of Interest
The authors declare no conflict of interest.

Figure 2 .
Figure 2. The local stress induced in the elliptical inclusion boundary due to global stress acting on homogenized medium [13]

Figure 4 .
Figure 4.The algorithm of a micromechanical damage model for a rock material under a constant dynamic strain rate

Figure 7 .
Figure 7. Variation of the uniaxial compressive strength by the imposed strain rate by scale law

Figure 8 .
Figure 8.The damage growth process in the simulated rock specimen at different times According to Figure8, the damage spreads in the rock sample with increment of time under applied dynamic loading.Figure9shows damage growth in the rock specimen simulated based on the proposed micromechanical damage model under dynamic uniaxial compressive loading.The results in Figures8 and 9show that the developed model is capable of predicting damage induced by frictional sliding along micro-cracks and associated wing-cracks nucleation.

Figure 9 Figure 9 .
Figure 9.The simulated stress-strain curve and damage evolution in the rock specimen