3-D Numerical Study on Progressive Failure Characteristics of Marbles under Unloading Conditions

: 3-D particle-based discrete element method (PB-DEM) was employed to numerically study the mechanical and progressive failure characteristics of pre-ﬁssured marble specimens under conventional triaxial unloading conditions. The microscopic parameters of PB-DEM for marble materials were calibrated using comparison with the previous experimental data. To systematically investigate the mechanical properties and the progressive failure characteristics of pre-ﬁssured marble specimens under the unloading conditions, a series of numerical specimens were simulated. The e ﬀ ects of ﬁssure geometric conditions, initial conﬁning pressures, and unloading rates on the mechanical and failure behaviors were explored via simulations. The present numerical results indicate that peak strength increased as the initial conﬁning pressures increased or the unloading rate decreased. Crack coalescence types and the ultimate failure modes in the pre-ﬁssured marble specimens were signiﬁcantly a ﬀ ected by the unloading stress paths. The present numerical results provide a better understanding of unloading mechanical and failure characteristics to scientists and engineers in rock mechanics and rock engineering.


Introduction
In nature, geologic tectonic activities result in the formation of different scale of discontinuous structures, such as pores, fissures, joints, and faults, in rock masses, which have significant influence on the mechanical properties, stability, and safety of rock masses in geotechnical engineering, such as mineral engineering, foundation pit engineering, tunnel engineering, and deep energy engineering [1][2][3].On the other hand, rock masses with discontinuities are located in a three-dimensional state in situ before being excavation.Since there exist different high in-situ geo-stresses, i.e., σ 1 > σ 2 > σ 3 0, excavations in rock masses result in stress redistribution around the cavities, which can be considered as an unloading process.The unloading conditions significantly affect the mechanical behaviors of rock masses, which may lead to slabbing/spalling and rock-bursts [1,2].Therefore, it is important and necessary to investigate the mechanical and failure behaviors of rock masses under unloading conditions, which is not only beneficial to understanding the unloading failure mechanism of rock masses, but also contributes to geotechnical engineering projects.
In the past decades, numerous laboratory investigations have been conducted to explore the deformation and failure characteristics of rocks, especially granite and marble, under unloading conditions in uniaxial and conventional triaxial compression tests.There is a great randomness in the discontinuous structure of natural rock.For the reproducibility of the experiment, pre-fissured rock materials are often adopted in experiments.Lajtai [4] defined wing and secondary cracks based on the fracture initiation and propagation characteristics in laboratory experimental observations.Some scholars also investigated crack propagation and coalescence behaviors and mechanism of rock specimens with two or multiple pre-existing fissures under uniaxial and conventional triaxial compressive loading conditions [5][6][7][8][9][10][11][12][13][14][15][16][17][18].However, laboratory experiments for rock and pre-fissured rock specimens under compressive loading conditions are expensive, time consuming, and not environmentally-friendly.Furthermore, the internal failure characteristics and deep physically-meaningful phenomena are difficult to observe due to lack of monitoring techniques [19,20].
In the past decades, as computer and storage technologies develop, different kinds of numerical methods, which can be classified into continuum and discontinuous approaches, were proposed to investigate the fracture problems in solids.Compared with the laboratory tests, numerical methodologies are cheap, repeatable, and easier to reveal the failure mechanism.For the continuum numerical methods, finite element with remeshing techniques [21,22], smeared elements [23][24][25], extended finite element method [26,27], peridynamics [28][29][30][31][32][33][34][35], and phase field method [36,37] have been successfully developed to investigate the failure evolutions of rock and rock masses under compressive loading.For example, Wang et al. [30] developed stress-based failure criteria in the non-ordinary state-based peridynamics to investigate the crack initiation, propagation, and coalescence in the pre-fissured rock specimens under uniaxial and biaxial compressive loading.Zhou et al. [38] proposed a 3D numerical mesoscale model to describe the time-independent deformation and fracturing of brittle rocks that accounts for material heterogeneity and local material degradation.Zhou et al. [36,37] implemented the phase field model into COMSOL to study the fracture behavior of rock specimens with pre-existing fissures subjected to the uniaxial compressive loads.
For discontinuous numerical methods, discrete element method [39][40][41][42][43][44][45][46][47], displacement discontinuity method [48][49][50], numerical manifold method [51][52][53], etc., were developed to simulate the failure processes of brittle solids.For instance, Zhang and Wang [39] applied discrete element method to investigate the crack initiation, propagation, and coalescence in rock-like specimens with two flaws under uniaxial compression.Wu and Wong [52] developed the enriched numerical manifold method to model cracking behaviors of rock masses.Li et al. [54,55] proposed the grain-based discrete element method to explore dynamic behaviors of fracturing and damage evolution of rock materials at the grain scale.Shang et al. [41] adopted three-dimensional discrete element method to conduct investigation of the fracture behavior of thermally degraded rocks.Song et al. [45,46] adopted a typical discrete element method code to simulate the progressive fatigue failure of concrete subjected to the multi-level cyclic loading, the micro-scale cracks can be characterized and the results show consistency with lab testing.Discrete element method, which discretizes the computational domain into a collection of particles with different shapes, are commonly adopted to simulate the failure of geomaterials.Compared with other numerical methods, the random collections of particles in discrete element methods are much closer to rocks and soils.
Although there exist many literatures on discrete element simulations of rock failure behaviors [39][40][41][42][43][44][45][46][47][56][57][58][59], the existing investigations on failure characteristics of rock and rock masses are mainly focused on the single compressive loading conditions.Failure characteristics and the mechanisms of rock and rock masses under unloading conditions are rarely explored using discrete element simulations, which are necessary to improve understanding the disaster mechanism of rock masses during tunnel excavations and mineral mining, as shown in Figure 1.Li et al. [60] used the PFC2D to carry out a series of studies on the unloading process of pre-cracked marble samples based on the physics experiment.The results reveal the effects of unloading rate and crack inclination on the strength and cracking performance of pre-cracked marble samples.However, the rock mass unloading process in engineering is more complicated, and requires further study.Therefore, numerical explorations on failure characteristics of rock masses under loading and unloading conditions are necessary, which is not only meaningful to revealing the unloading failure mechanism for scientist, but also beneficial to geotechnical engineers.
Appl.Sci.2020, 10, x FOR PEER REVIEW 3 of 33 conditions are necessary, which is not only meaningful to revealing the unloading failure mechanism for scientist, but also beneficial to geotechnical engineers.Most current researches focus on mechanical and failure behaviors of intact rock specimens under unloading conditions, but there are few systematic studies on the progressive failure processes of pre-fissured rock specimens under unloading conditions [61].Although the failure characteristics with different unloading rates are summarized, the influence of initial confining pressures on mechanical and failure behaviors is not carried out.The existing literature provides that failure mechanism of rock specimens under unloading conditions is very different from that under loading conditions, which indicates the necessary to research the unloading mechanical and failure characteristics of pre-fissured rock specimens.
In this work, particle-based discrete element method (PB-DEM) was applied to investigate the mechanical failure characteristics of pre-fissured rock specimens under the triaxial compression unloading conditions.Microscopic parameters in PB-DEM for marble were calibrated using trial-anderror approach based on the previous experiments.Then, different fissured marble specimens with various fissure length , fissure inclination angle , ligament angle , and ligament length  under the conventional triaxial loading and unloading conditions were simulated using PB-DEM.Effects of fissure distribution, initial confining pressures, and unloading rates on the mechanical and failure characteristics of fissured marble specimens were numerically investigated.The crack types and crack coalescence modes in fissured marble specimens under the conventional triaxial compressive unloading conditions were summarized.Numerical results can be used to distinguish the spatial distribution of tensile and shear cracks, which is difficult to obtain directly in the experiment.
The remains of the article are organized as follows.In Section 2, PB-DEM is briefly described and its microscopic parameters for marble specimen are calibrated.In Section 3, a series of fissured marble specimens under the triaxial compression loading and unloading conditions is simulated using PB-DEM.Effects of fissure distribution, confining pressures and unloading rates on mechanical and failure behaviors is numerically explored in Section 4. Conclusions of the whole manuscript are drawn in Section 5. Most current researches focus on mechanical and failure behaviors of intact rock specimens under unloading conditions, but there are few systematic studies on the progressive failure processes of pre-fissured rock specimens under unloading conditions [61].Although the failure characteristics with different unloading rates are summarized, the influence of initial confining pressures on mechanical and failure behaviors is not carried out.The existing literature provides that failure mechanism of rock specimens under unloading conditions is very different from that under loading conditions, which indicates the necessary to research the unloading mechanical and failure characteristics of pre-fissured rock specimens.
In this work, particle-based discrete element method (PB-DEM) was applied to investigate the mechanical failure characteristics of pre-fissured rock specimens under the triaxial compression unloading conditions.Microscopic parameters in PB-DEM for marble were calibrated using trial-and-error approach based on the previous experiments.Then, different fissured marble specimens with various fissure length a, fissure inclination angle α, ligament angle β, and ligament length b under the conventional triaxial loading and unloading conditions were simulated using PB-DEM.Effects of fissure distribution, initial confining pressures, and unloading rates on the mechanical and failure characteristics of fissured marble specimens were numerically investigated.The crack types and crack coalescence modes in fissured marble specimens under the conventional triaxial compressive unloading conditions were summarized.Numerical results can be used to distinguish the spatial distribution of tensile and shear cracks, which is difficult to obtain directly in the experiment.
The remains of the article are organized as follows.In Section 2, PB-DEM is briefly described and its microscopic parameters for marble specimen are calibrated.In Section 3, a series of fissured marble specimens under the triaxial compression loading and unloading conditions is simulated using PB-DEM.Effects of fissure distribution, confining pressures and unloading rates on mechanical and failure behaviors is numerically explored in Section 4. Conclusions of the whole manuscript are drawn in Section 5.

Particle-Based Discrete Element Method (PB-DEM)
Particle-based discrete element method (PB-DEM) is one typical method of DEM, which reflects mechanical and failure behaviors using the motion and interaction among the clusters of rigid particles [61][62][63].The interaction system is treated as a dynamic equilibrium process generated by the motion of wall and particles, propagated through the contact between particles, as shown in Figure 2. It can be observed from Figure 2 that PB-DEM consists of linear contact bond model (CBM) and linear parallel bond model (LPBM), where LPBM is used for the contact relationship between rock particles.

Particle-Based Discrete Element Method (PB-DEM)
Particle-based discrete element method (PB-DEM) is one typical method of DEM, which reflects mechanical and failure behaviors using the motion and interaction among the clusters of rigid particles [61][62][63].The interaction system is treated as a dynamic equilibrium process generated by the motion of wall and particles, propagated through the contact between particles, as shown in Figure 2. It can be observed from Figure 2  When breakage of bonds between particles occurs, bond stiffness between particles is instantaneously removed.Meanwhile, the contact stiffness remains active as long as the particles remain in contact.Therefore, parallel bond model is considered to be a more realistic contact model for rock-like materials because of the reduction in material stiffness during the process of bond breaking [63].

Parameter Calibration
In PB-DEM, microscopic parameters commonly cannot directly correspond to the macroscopic parameters of geo-materials.Thus, microscopic parameters in PB-DEM require calibration through comparison of failure patterns and mechanical responses between numerical simulations the laboratory physical tests.
According to the previous studies [39,61,62,64], a series of uniaxial and triaxial compression experiments with different parameters was employed.The calibration of the marble model parameters was performed using a trial-and-error method that compares the numerical simulation results with the laboratory experiments results [65].The schematic diagram of the compression test is presented in Figure 3a,b.When breakage of bonds between particles occurs, bond stiffness between particles is instantaneously removed.Meanwhile, the contact stiffness remains active as long as the particles remain in contact.Therefore, parallel bond model is considered to be a more realistic contact model for rock-like materials because of the reduction in material stiffness during the process of bond breaking [63].

Parameter Calibration
In PB-DEM, microscopic parameters commonly cannot directly correspond to the macroscopic parameters of geo-materials.Thus, microscopic parameters in PB-DEM require calibration through comparison of failure patterns and mechanical responses between numerical simulations the laboratory physical tests.
According to the previous studies [39,61,62,64], a series of uniaxial and triaxial compression experiments with different parameters was employed.The calibration of the marble model parameters was performed using a trial-and-error method that compares the numerical simulation results with the laboratory experiments results [65].The schematic diagram of the compression test is presented in Figure 3a,b.To avoid influence of model interface size on macroscopic parameters, cylindrical and rectangular models were used for calibration experiments.The ultimate failure patterns of intact marble specimens obtained from the present numerical predictions and the previous experimental observations are presented in Figure 4.It can be observed from Figure 4c that the main crack in the previous laboratory experiments results from the diagonal direction of the sample and a secondary  To avoid influence of model interface size on macroscopic parameters, cylindrical and rectangular models were used for calibration experiments.The ultimate failure patterns of intact marble specimens obtained from the present numerical predictions and the previous experimental observations are presented in Figure 4.It can be observed from Figure 4c that the main crack in the previous laboratory experiments results from the diagonal direction of the sample and a secondary crack in the opposite direction, which is consistent with the present numerical experimental predictions, as shown in Figure 4d,e.The flaking damage at the top and bottom of the sample in the numerical predictions, as shown in Figure 4a,b, is in agreement with the bottom rock spalling presented in Figure 4c.The uniaxial compressive axial stress-strain curves of marble specimens obtained from the present numerical simulations and the previous laboratory experiments are shown in Figure 5.It can be observed from Figure 5 that the marble specimen performs significantly brittle failure after peak strength.The stress-strain curve of a physical experiment has a concave shape due to rock compaction at the beginning of the curve, which is different from the numerical predictions.The numerical model is dense inside due to the expansion during particle generation.Therefore, the axial stress-strain curve directly enters the elastic phase.Under the condition of confining pressure of 30 MPa, the results of experiments and numerical predictions are consistent.Moreover, it can be observed from Figure 5 that peak strength of 85.9 (160.5)MPa and deformation modulus (E) of 35.9 (53.2) GPa obtained from the present numerical predictions are similar to values of 83.5 (154.4)MPa and 36.0 (53.6)GPa (E) in the laboratory tests.Since the unloading speed was a large value in the experiment, the contact damper coefficient was set to 0.5 in order to solve the dynamic effect of particles.The mode of damping is full normal and full shear, so the interaction between particles is viscous.The Poisson ratio of the intact sample obtained by parameter calibration was 0.21, but no clear corresponding parameter was found in previous laboratory experiments [65].Thus, the micro-parameters after calibration of particle and contact are listed in Table 1, which can be used for a series of loading and unloading experiments on marble specimens in the present study.

Numerical Model Setup
The complete marble model was generated using the microscopic parameters above, with a model size of 100 × 50 × 50 mm 3 , containing 164,890 particles with a radius between 0.5 and 0.8 mm, the particles were connected by 820,396 parallel bond models.The geometry of the pre-cracked rock sample in the simulation was simulated by deleting the particle elements on the designated surface.The loading and geometrical conditions of the pre-existing fissures in the marble specimens are plotted in Figure 3b,c.Moreover, it can be observed from Figure 3c that a is crack length, α is crack inclination angle, β is ligament inclination angle, and b is ligament length.
The details of geometries of pre-fissured marble specimens are shown in Figure 6 and Table 2.It can be found from Figure 6 and Table 2 that crack inclination angles were equal to 45 • , 30 • , 60 • , and 75 • in the pre-fissured marble specimens of Type-A, Type-B, Type-C, and Type-D, respectively.In the pre-fissured marble specimens of Type-A ~Type-D, crack lengths were fixed as 24 mm, and ligament lengths were fixed as 33 mm, which resulted in various ligament inclination angles in these marble specimens, as shown in Figure 6 and Table 2.As shown in Figure 6 and Table 2, for Type-A, Type-E, Type-F, the initial fissure inclined angles and initial fissure lengths are fixed as 45 • and 24 mm.The ligament lengths between two pre-existing fissures are 33 mm, 48 mm, and 24 mm in Type-A, Type-E, and Type-F marble specimens, respectively.Similarly, it can be also observed from Figure 6 and Table 2 that the initial fissure inclined angle and ligament length between the two pre-existing fissures were fixed as 45 • and 33 mm in Type-A, Type-G, and Type-H marble specimens.Meanwhile, the initial lengths of pre-existing fissures in Type-A, Type-G and Type-H marble specimens equated to 24 mm, 16 mm, and 20 mm, respectively.
mm.The ligament lengths between two pre-existing fissures are 33 mm, 48 mm, and 24 mm in Type-A, Type-E, and Type-F marble specimens, respectively.Similarly, it can be also observed from Figure 6 and Table 2 that the initial fissure inclined angle and ligament length between the two pre-existing fissures were fixed as 45° and 33 mm in Type-A, Type-G, and Type-H marble specimens.Meanwhile, the initial lengths of pre-existing fissures in Type-A, Type-G and Type-H marble specimens equated to 24 mm, 16 mm, and 20 mm, respectively.Table 2. Geometric and loading parameters of pre-fissured marble specimens.The stress redistribution caused by rock unloading usually manifests as the stress increase in one direction and a release in the orthogonal direction.Therefore, the stress phase of rock unloading is divided as follows: (i) Initial stress phase: At this stage, vertical stress σ z and lateral stress σ x , σ y were simultaneously increased to take the stress to initial stress σ 0 , and then the velocity of the particle is reset.(ii) Stress release phase: At this stage, the lateral stress σ x , σ y were unloaded and the vertical stress σ z increases until the sample reaches the peak stress state.

Geometries of
The vertical loading velocity was applied at a rate of 0.2 mm/s during initial stress phase.The confining pressure maintained the same value as the velocity pressure through the electro-hydraulic servo control system.In the subsequent stress phase, the vertical loading velocity was continuously loaded, and the stress path of the confining pressure is expressed as: where V L is the speed of wall loading; E is the elasticity modulus; t is the time of the experimental; t 0 is the moment when the loading stops; σ 0 is the initial confining pressure of the design; V U is the speed of wall unloading; L 0 is the height of the sample.
The unloading stress path is linear during the stress release phase.Different unloading rates and initial stresses were chosen to study the response characteristics of the rock sample during unloading.Based on the previous literatures [21], the unloading period was controlled in the vicinity of 2 m/s in the experiment.Since the unloading process was speed-controlled and the initial stress level was different, the unloading period was different.In each group of unloading experiments, the unloading rates (V U ) were 0.01 mm/s, 0.1 mm/s, 0.2 mm/s, 0.3 mm/s, and 0.4 mm/s.The experiment in which the unloading rate was 0.01 mm/s was used to study the confined compression simulation for comparison.In this paper, due to the influence of pre-cracking, the strength difference between different groups of samples was large, which is not conducive to uniform confining pressure control.Therefore, the initial confining pressure (σ 0 ) was adjusted according to the uniaxial strength of the specimen.The 60%, 80%, 100%, and 120% of the uniaxial strength were taken as the initial confining pressure.In summary, the stress path contained four initial confining pressures and four unloading rates as shown in Table 2 and Figure 1b.

Uniaxial Compressive Simulations Results
Based on the calibrated microscopic material parameters in Section 3 and geometrical parameters in Table 2, eight different kinds of pre-fissured marble specimens were constructed using discrete element method.In order to study the mechanical properties and failure behaviors of pre-fissured marble specimens subjected to uniaxial compressive loads, numerical results were firstly compared with the previous experimental data [66], as shown in Table 3.It can be found from Table 3 that numerical results are in good agreement with the previous experimental data in laboratory [66].The difference of mechanical properties of Type-A and Type-B pre-fissured marble specimens between numerical results between experimental data [65] was mainly caused by many uncertain factors in laboratory tests.The 3-D final failure patterns of the uniaxial compression of the pre-fissured marble specimens are plotted in Figure 7.In Figure 7, the red portion indicates shear micro-cracks, and the green portion indicates tensile micro-cracks.It can be found from Figure 7 that the final failure modes in different types of pre-fissured marble specimens under uniaxial compressive loading conditions were mainly dominated by tensile cracks.Meanwhile, the shear cracks were mainly distributed around the pre-existing fissure tips.Furthermore, 3-D final failure patterns also demonstrates that the ultimate failure modes of different marble specimens belong to splitting tensile failure mode.Since, the present researches are mainly focused on the mechanical and failure characteristics of pre-fissured marble specimens under unloading conditions, numerical results of failure behaviors of marble specimens can be considered as the comparing references.The 3-D final failure patterns of the uniaxial compression of the pre-fissured marble specimens are plotted in Figure 7.In Figure 7, the red portion indicates shear micro-cracks, and the green portion indicates tensile micro-cracks.It can be found from Figure 7 that the final failure modes in different types of pre-fissured marble specimens under uniaxial compressive loading conditions were mainly dominated by tensile cracks.Meanwhile, the shear cracks were mainly distributed around the preexisting fissure tips.Furthermore, 3-D final failure patterns also demonstrates that the ultimate failure modes of different marble specimens belong to splitting tensile failure mode.Since, the present researches are mainly focused on the mechanical and failure characteristics of pre-fissured marble specimens under unloading conditions, numerical results of failure behaviors of marble specimens can be considered as the comparing references.

Conventional Triaxial Compression Simulations
A series of conventional triaxial compression simulations were then carried out on eight types of marble specimens based on the loading condition listed in Table 2.The 3-D final failure patterns of eight pre-fissured marble specimens in the conventional triaxial compressive simulations were performed in Figure 8.It can be observed from Figure 8 that the micro-cracks in pre-fissured marble specimens under conventional triaxial compression conditions are mainly shear cracks, which occur in rock-bridge area in rock specimens.When the initial fissure length  and ligament length  are fixed, it can be observed from Figure 8a-d that crack coalescence type varied from shear crack coalescence to mixed tensile-shear crack coalescence as the initial fissure inclination angle increased.Furthermore, the final failure mode of pre-fissured marble specimens under conventional triaxial compressive loading conditions transformed from shear failure mode to the mixed tensile-shear failure mode with the increase of the initial fissure inclination angle, as shown in Figure 8a-d.When the initial fissure inclination angle  and initial fissure length  were fixed, it can be found from Figure 8a,e,f that the mixed tensile-shear crack coalescence appeared in the rock bridge region, which resulted in the ultimate mixed tensile-shear failure of pre-fissured marble specimens.When the initial fissure inclination angle  and ligament length  were fixed, it can be observed from Figure

Conventional Triaxial Compression Simulations
A series of conventional triaxial compression simulations were then carried out on eight types of marble specimens based on the loading condition listed in Table 2.The 3-D final failure patterns of eight pre-fissured marble specimens in the conventional triaxial compressive simulations were performed in Figure 8.It can be observed from Figure 8 that the micro-cracks in pre-fissured marble specimens under conventional triaxial compression conditions are mainly shear cracks, which occur in rock-bridge area in rock specimens.When the initial fissure length a and ligament length b are fixed, it can be observed from Figure 8a-d that crack coalescence type varied from shear crack coalescence to mixed tensile-shear crack coalescence as the initial fissure inclination angle increased.Furthermore, the final failure mode of pre-fissured marble specimens under conventional triaxial compressive loading conditions transformed from shear failure mode to the mixed tensile-shear failure mode with the increase of the initial fissure inclination angle, as shown in Figure 8a-d.When the initial fissure inclination angle α and initial fissure length a were fixed, it can be found from Figure 8a,e,f that the mixed tensile-shear crack coalescence appeared in the rock bridge region, which resulted in the ultimate mixed tensile-shear failure of pre-fissured marble specimens.When the initial fissure inclination angle α and ligament length b were fixed, it can be observed from Figure 8a,g,h that crack coalescence type in the rock bridge region belonged to the mixed tensile-shear mode.However, the ultimate shear failure mode of pre-fissured marble specimens transformed from shear failure mode to mixed tensile-shear failure mode with the increase of the initial fissure length a.This phenomenon was caused by the inclined secondary shear cracks initiate from pre-fissure tips and propagated towards the upper and lower boundaries when the initial fissure length was equal to 16 mm and 20 mm, as shown in Figure 8a,g,h.
Appl.Sci.2020, 10, x FOR PEER REVIEW 12 of 33 However, the ultimate shear failure mode of pre-fissured marble specimens transformed from shear failure mode to mixed tensile-shear failure mode with the increase of the initial fissure length .This phenomenon was caused by the inclined secondary shear cracks initiate from pre-fissure tips and propagated towards the upper and lower boundaries when the initial fissure length was equal to 16 mm and 20 mm, as shown in Figure 8a,g,h.

Confining Pressure Unloading Simulations
To investigate the mechanical and failure behaviors of pre-fissured marble specimens under loading and unloading conditions, eight types of pre-fissured marble specimens were simulated using PB-DEM.The geometrical conditions and material properties in these numerical simulations were the same as the previous ones.In these simulations, the unloading rate (  ) was adopted as 0.3 mm/s and the initial confining pressure ( 0 ) was taken as .0 .
The 3-D ultimate failure patterns of pre-fissured marble specimens under unloading conditions are presented in Figure 9.It can be observed from Figure 9 that the number of shear cracks in the

Confining Pressure Unloading Simulations
To investigate the mechanical and failure behaviors of pre-fissured marble specimens under loading and unloading conditions, eight types of pre-fissured marble specimens were simulated using PB-DEM.The geometrical conditions and material properties in these numerical simulations were the same as the previous ones.In these simulations, the unloading rate (V U ) was adopted as 0.3 mm/s and the initial confining pressure (σ 0 ) was taken as 1.0σ ucs .
The 3-D ultimate failure patterns of pre-fissured marble specimens under unloading conditions are presented in Figure 9.It can be observed from Figure 9 that the number of shear cracks in the unloading simulation results was significantly larger than that in the uniaxial compression simulations.The master cracks developed at the tip of the pre-crack, and the cracks that originated at the tip of the pre-crack had a composite feature of tension and shear.For the model with larger ligament inclination angle β, the development of cracks was mainly concentrated in the ligament, which was agreement with the previous experimental observations [67-71].As ligament inclination angle β decreased, the number of shear cracks gradually increased, as shown in Figure 9a,e,f.Similar to the 3-D failure patterns in the uniaxial compression simulations, anti-wing cracks appeared in the numerical specimens of Type-B, Type-E, and Type-G.The same as the uniaxial compression condition, the anti-wing cracks appeared in the models of Type B, E, G, and H, with the difference that the dip angle at the end of the anti-wing crack and the number of shear cracks were larger.The failure characteristics of the samples change from splitting failure to composite failure compared with the results of uniaxial compression simulation, which was mainly manifested in the increase of crack deflection angle of each group of wing crack and anti-wing crack, and the increase of ligament cracking such as Type B, D, and E. The most significant effect of unloading experiment was that the confining pressure restricted the development of tension cracks represented by wing cracks.Similar phenomena can be observed in the previous experimental As ligament inclination angle β decreased, the number of shear cracks gradually increased, as shown in Figure 9a,e,f.Similar to the 3-D failure patterns in the uniaxial compression simulations, anti-wing cracks appeared in the numerical specimens of Type-B, Type-E, and Type-G.The same as the uniaxial compression condition, the anti-wing cracks appeared in the models of Type B, E, G, and H, with the difference that the dip angle at the end of the anti-wing crack and the number of shear cracks were larger.The failure characteristics of the samples change from splitting failure to composite failure compared with the results of uniaxial compression simulation, which was mainly manifested in the increase of crack deflection angle of each group of wing crack and anti-wing crack, and the increase of ligament cracking such as Type B, D, and E. The most significant effect of unloading experiment was that the confining pressure restricted the development of tension cracks represented by wing cracks.Similar phenomena can be observed in the previous experimental observations [70,71].
The deviatoric stress-displacement curve and the microcrack growth curve of each sample are plotted in Figure 10.It can be seen that the peak strength of the sample was significantly improved compared to uniaxial compression and unloading conditions.Among these, the model with a small β value such as Type B, F, and G did not show a drop in peak strength.Although the models under different loading conditions were quite different in strength, the total number of microcracks in the sample destruction process was the same.
The deviator stress-axial strain curves and micro-crack evolution curves of the Type A to Type D pre-fissured marble specimens are recorded in Figure 10.It can be seen that the pack stress of each of the specimens under unloading conditions had a certain degree of improvement compared to uniaxial compression condition.Unlike the linear increase and the sudden drop of the deviatoric stress curve under uniaxial compression, the curve under unloading conditions were more complicated.For a better explanation, the unloading experiment process was divided into the confining pressure loading stage, the elastic stage, the yield stage, and the destruction stage according to the characteristics of the deviatoric stress growth, corresponding to the OA, AB, BC, and CD phases of the curve of Figure 10, respectively.During the loading phase, the triaxial pressures were approximately equal and the deviatoric stress was approximately equal to zero.In the elastic phase, the vertical stress increased, the confining pressure unloaded, and the deviatoric stress increased linearly.Different types of specimens exhibited different deviatoric stress characteristics during the yield phase.One of the manifestations was that the partial stress continued to increase linearly, but the growth rate decreased, and fell rapidly after reaching the peak stress, such as Type C, D, and E. The reason was that the specimens were not broken after the unloading was completed, and entered the uniaxial compression state.Another manifestation was that the partial stress curve fluctuated and did not fall after reaching the peak stress in the yielding phase, such as Type A, B, F, G, and H.These types of specimens reached the theoretical peak stress during the unloading process, but the deviatoric stress did not fall due to the confining pressure limit.In the later stages, the deviatoric stress naturally decreased with loading and the specimens lost their capacity.
It can be found from Figure 10 that cracks in the simulation specimens were not generated synchronously with the deformation, but were a gradual evolution process.The result of the Type A will be described below as an example.The number of microcracks had different growth characteristics depending on the stage in which the deviatoric stress curve was located.The loading and unloading phases were distinguished based on the previous experimental studies [70,71].In the loading phase, the OA phase, the number of cracks generated was small with a maximum of 16 times per calculation step, and the crack generation time was dispersed.After entering the elastic phase, the AB phase, the number of microcracks generated per unit time began to increase, with a maximum of 39 times per calculation step.The generation of cracks was continuous in this phase, but the total amount of microcracks generated was still small.During the yield phase (BC phase), the number of microcracks began to increase significantly, with a maximum of 347 micro-cracks per calculation step.Corresponding to this: Each specimen generated a series of macroscopic cracks that expanded rapidly along the end of the pre-fissures and underwent irreversible plastic damage.In the failure phase (CD phase), the growth rate of micro-cracks decreased significantly, and the newly generated micro-cracks were mainly distributed around the pre-existing fissure tips.The deviator stress-axial strain curves and micro-crack evolution curves of the Type A to Type D pre-fissured marble specimens are recorded in Figure 10.It can be seen that the pack stress of each

Discussions
In this section, influences of pre-existing fissures' distributions, initial confining pressures, and unloading rates on the mechanical and failure characteristics, including (i) 3-D final failure patterns, (ii) number of microscopic cracks, (iii) peak strength, and (iv) effective modulus, of pre-fissured marble specimens under unloading conditions were numerically investigated.

Influences of Preexisting Cracks' Distributions
Based on the previous numerical results and experimental data [32, 65,66], there were five different cracks, i.e., wing crack, horsetail crack, anti-wing crack, inclined secondary crack, and coplanar secondary crack, observed from the 3-D numerical simulations by discrete element method in this research.The schematic diagram of crack type is shown in Figure 11.Wing crack, horsetail crack, and anti-wing crack belong to the tensile cracks.On the other hand, the inclined secondary crack belongs to the mixed tensile-shear crack and coplanar secondary crack is a type of shear crack, as shown in Figure 11.At the same time, three different modes of crack coalescence, including tensile crack coalescence, mixed tensile-shear crack coalescence, and shear crack coalescence, in the rock bridge areas are also drawn in Figure 12.To better understand the influence of initial fissure length , initial fissure inclination angle , ligament length , and ligament inclination angle  on the mechanical behaviors of pre-fissured marble specimens under unloading conditions, relationships between numerically predicted axial deviatoric stress and other geometric factors are summarized in Figure 13.It can be observed from Figure 13 that the numerically predicted axial deviatoric stress decreased as values of the geometrical parameters, i.e., , , , , increased.As shown in Figure 13, influences of pre-existing fissure length  and ligament length  on the peak strength was smaller than that of initial fissure To better understand the influence of initial fissure length , initial fissure inclination angle , ligament length , and ligament inclination angle  on the mechanical behaviors of pre-fissured marble specimens under unloading conditions, relationships between numerically predicted axial deviatoric stress and other geometric factors are summarized in Figure 13.It can be observed from Figure 13 that the numerically predicted axial deviatoric stress decreased as values of the geometrical parameters, i.e., , , , , increased.As shown in Figure 13, influences of pre-existing fissure To better understand the influence of initial fissure length a, initial fissure inclination angle α, ligament length b, and ligament inclination angle β on the mechanical behaviors of pre-fissured marble specimens under unloading conditions, relationships between numerically predicted axial deviatoric stress and other geometric factors are summarized in Figure 13.It can be observed from Figure 13 that the numerically predicted axial deviatoric stress decreased as values of the geometrical parameters, i.e., a, b, α, β, increased.As shown in Figure 13, influences of pre-existing fissure length a and ligament length b on the peak strength was smaller than that of initial fissure inclination angle α and ligament inclination angle β.Moreover, it can be also found from Figure 13 that the unloading rate V U had an important influence on the deviatoric strength effect.With the increase of the unloading rate, the peak value of axial deviatoric stress decreased, as shown in Figure 13.This shows the same pattern as the results of previous studies [8,10,20].However, due to the different materials and loading conditions, the geometrical parameters of the sample had different effects on the partial stress strength.
Appl.Sci.2020, 10, x FOR PEER REVIEW 17 of 33 13.This shows the same pattern as the results of previous studies [8,10,20].However, due to the different materials and loading conditions, the geometrical parameters of the sample had different effects on the partial stress strength.

Influence of Confining Pressures
To investigate the effect of confining pressures on the mechanical and failure behaviors of prefissured marble specimens, eight different type of pre-fissured marble specimens subjected to different confining pressures in the unloading simulations.The different confining pressures are equal to 0%  , 0%  , 00%  , and 20%  , as shown in Table 2.The 3-D ultimate failure patterns of Type-A ~ Type-D pre-fissured marble specimens under different confining pressures in the unloading simulations are presented in Figure 14.

Influence of Confining Pressures
To investigate the effect of confining pressures on the mechanical and failure behaviors of pre-fissured marble specimens, eight different type of pre-fissured marble specimens subjected to different confining pressures in the unloading simulations.The different confining pressures are equal to 60%σ ucs , 80%σ ucs , 100%σ ucs , and 120%σ ucs , as shown in Table 2.The 3-D ultimate failure patterns of Type-A ~Type-D pre-fissured marble specimens under different confining pressures in the unloading simulations are presented in Figure 14.The magnitude of initial confining pressures in the unloading simulations had significant influence on the generation of wing cracks.It can be observed from Figure 14 that crack initiation angles of wing cracks gradually decreased as the initial confining pressures increase from 0%  to 0%  .When the initial confining pressure continued to increase, wing cracks initiated from the tips of the pre-existing fissures gradually changed to horsetail cracks, as shown in Figure 14.For the pre-fissured marble specimens with low inclination angles, i.e., Type-A and Type-B, it can be found from Figure 14 that crack coalescence mode in the rock bridge region transformed from the tensile crack coalescence to the mixed tensile-shear crack coalescence as the initial confining pressure increase.However, for the pre-fissured marble specimens with large inclination angles, i.e., Type-C and Type-D, crack coalescence mode in the rock bridge region belonged to the mixed crack coalescence, which was not affected by the variation of initial confining pressures, as shown in Figure 14.Furthermore, Figure 14 also indicates that the number of shear cracks increased with the increase of the initial confining pressures.These phenomena demonstrate that shear cracks may be induced by the initial confining pressures in the unloading simulations.
Figure 15 presents numerically-predicted axial deviatoric stress-axial strain curves of different marble specimens under unloading conditions with various initial confining pressures.It can be observed from Figure 15 that initial confining pressure had important influence on the mechanical and deformation behaviors of pre-fissured marble specimens, which will be analyzed in detail as follows.The magnitude of initial confining pressures in the unloading simulations had significant influence on the generation of wing cracks.It can be observed from Figure 14 that crack initiation angles of wing cracks gradually decreased as the initial confining pressures increase from 60%σ ucs to 80%σ ucs .When the initial confining pressure continued to increase, wing cracks initiated from the tips of the pre-existing fissures gradually changed to horsetail cracks, as shown in Figure 14.For the pre-fissured marble specimens with low inclination angles, i.e., Type-A and Type-B, it can be found from Figure 14 that crack coalescence mode in the rock bridge region transformed from the tensile crack coalescence to the mixed tensile-shear crack coalescence as the initial confining pressure increase.However, for the pre-fissured marble specimens with large inclination angles, i.e., Type-C and Type-D, crack coalescence mode in the rock bridge region belonged to the mixed crack coalescence, which was not affected by the variation of initial confining pressures, as shown in Figure 14.Furthermore, Figure 14 also indicates that the number of shear cracks increased with the increase of the initial confining pressures.These phenomena demonstrate that shear cracks may be induced by the initial confining pressures in the unloading simulations.
Figure 15 presents numerically-predicted axial deviatoric stress-axial strain curves of different marble specimens under unloading conditions with various initial confining pressures.It can be observed from Figure 15 that initial confining pressure had important influence on the mechanical and deformation behaviors of pre-fissured marble specimens, which will be analyzed in detail as follows.The magnitude of initial confining pressures in the unloading simulations had significant influence on the generation of wing cracks.It can be observed from Figure 14 that crack initiation angles of wing cracks gradually decreased as the initial confining pressures increase from 60% to 80% .When the initial confining pressure continued to increase, wing cracks initiated from the tips of the pre-existing fissures gradually changed to horsetail cracks, as shown in Figure 14.For the pre-fissured marble specimens with low inclination angles, i.e., Type-A and Type-B, it can be found from Figure 14 that crack coalescence mode in the rock bridge region transformed from the tensile crack coalescence to the mixed tensile-shear crack coalescence as the initial confining pressure increase.However, for the pre-fissured marble specimens with large inclination angles, i.e., Type-C and Type-D, crack coalescence mode in the rock bridge region belonged to the mixed crack coalescence, which was not affected by the variation of initial confining pressures, as shown in Figure 14.Furthermore, Figure 14 also indicates that the number of shear cracks increased with the increase of the initial confining pressures.These phenomena demonstrate that shear cracks may be induced by the initial confining pressures in the unloading simulations.
Figure 15 presents numerically-predicted axial deviatoric stress-axial strain curves of different marble specimens under unloading conditions with various initial confining pressures.It can be observed from Figure 15  Each type of marble had different axial deviatoric stress growth rates during the confining pressure loading phase.At the elastic stage, confining pressure decreased while the vertical stress increased, and the elastic deformation dominated the axial deviatoric stress-axial strain curve.Although the various marble specimens exhibited different equivalent elastic moduli at this stage, the equivalent elastic moduli of the same set under different initial confining pressure conditions were identical to each other.Some irreversible processes also occurred in this phase, such as the initiation of a certain amount of micro-cracks and the extension of the tip crack, but the linear stress behavior of the sample did not change.Due to the stress concentration at the crack tip, the macro Each type of marble had different axial deviatoric stress growth rates during the confining pressure loading phase.At the elastic stage, confining pressure decreased while the vertical stress increased, and the elastic deformation dominated the axial deviatoric stress-axial strain curve.Although the various marble specimens exhibited different equivalent elastic moduli at this stage, the equivalent elastic moduli of the same set under different initial confining pressure conditions were identical to each other.Some irreversible processes also occurred in this phase, such as the initiation of a certain amount of micro-cracks and the extension of the tip crack, but the linear stress behavior of the sample did not change.Due to the stress concentration at the crack tip, the macro crack rapidly expanded under increasing deviatoric stress.As the crack expanded further, a portion of the lower stress sample entered the yield phase.The most obvious feature of the yield stage was the sudden change in the slope of the partial stress curve.At this stage, the confining pressure had been unloaded to 0, and the drop rate of different curves was almost the same.It can be inferred that the form of the curve drop was affected by the confining pressure in the current state instead of the stress path, which is the same as the results of the previous unloading studies [21,60,61].Although the curve was still linear at this time, the sample produced a large number of micro-cracks and plastic damage occurred.This stage is common in models with low strength such Type-C ~Type-F or models with low initial confining pressure.It can be seen that confining pressure had an impact on rock failure mode.For a model with a large stress level, there was no obvious period of yielding, or it can be considered that the yielding stage overlaps with the stage of failure.After entering the failure stage, the marble sample formed a master crack that supported the axial stress by the friction between the particles on the surface.The deviatoric stress of each sample decreased rapidly, but due to the lateral force limitation, the characteristic of the curve reduction was different from the phenomenon of the intensity drop that occurs in the uniaxial compression simulation.
Figure 16 presents the number of tensile and shear cracks during the failure evolution processes in the unloading simulations with various initial confining pressures.It can be seen from Figure 16 that the large initial confining pressure can delay the generation time of micro-cracks, which indicates that large confining pressure can increase the peak strength and improve the integrity of marble specimens at the early elastic stage.In addition, number of tensile cracks was much larger than that of shear cracks.Meanwhile, as the initial confining pressure increased, the number of shear cracks significantly increased, which indicates that confining pressures promote the occurrence of shear cracks.Two different types of micro-cracks mainly appeared at the yielding stage.In contrast to tensile cracks, shear cracks did not increase after pre-fissured marble specimens at the failure stage.The influence of initial confining pressures on the peak strength and equivalent elastic modulus of different marble specimens is presented in Figure 17.The influencing degree of the initial confining pressure on the peak strength varied depending on the pre-crack of the sample.The peak strength of Type-C ~ Type-E was hardly affected by initial confining pressure, while the other samples increased with the increase of the initial confining pressure.The increase in strength was mainly due to the confining pressure that limited the development of the sample crack, but did not affect the elastic modulus of the sample.Therefore, the equivalent elastic modulus under different initial confining pressures were substantially the same of the respective sets of models as shown in Figure 17.The influence of initial confining pressures on the peak strength and equivalent elastic modulus of different marble specimens is presented in Figure 17.The influencing degree of the initial confining pressure on the peak strength varied depending on the pre-crack of the sample.The peak strength of Type-C ~Type-E was hardly affected by initial confining pressure, while the other samples increased with the increase of the initial confining pressure.The increase in strength was mainly due to the confining pressure that limited the development of the sample crack, but did not affect the elastic modulus of the sample.Therefore, the equivalent elastic modulus under different initial confining pressures were substantially the same of the respective sets of models as shown in Figure 17.For wing cracks, the unloading rate had a significant effect on the inclination angle.The behavior was as follows: As the unloading rate increased, the inclination of the wing crack increased, which was more obvious in Type A and Type B. Unlike the initial confining pressure, the unloading rate also had a significant effect on the length of the wing crack.The length of the wing crack in each model increased as the unloading rate increased.For the cracks in the rock ligament, the unloading rate and the initial confining pressure showed similar effects: As the unloading rate decreased, the number of the shear microcracks in the ligament was significantly increased.Because for the specimens at the same strain state, the reduction in the unloading rate was somewhat equivalent to the increase in confining pressure.For the Type B model, in the process of reducing the unloading rate, not only did the number of shear cracks in the ligament increase, but also the cracks in the ligament became the main cracks of the sample.Different from the influence of initial confining pressure, in terms of anti-wing crack, the unloading rate had little influence on the angle of the crack, but the effect was mainly as follows: As the unloading rate decreased, the shear microcrack number rose.Macroscopically, the reduction of the unloading rate limited the development of tensile cracks, such as wing cracks and anti-wing cracks, which made the shear cracks develop space, and made the macro tensile cracks change into shear crack.Similar phenomena can be observed in the previous experimental observations [72].

Influence of Unloading Rates
Figure 19 shows the tensile crack and shear crack number curves for the Type A model under unloading conditions with various unloading rates.It can be seen from the figure that as the unloading rate decreased, the generation process of micro-cracks was delayed.The effect of the unloading rate on the crack initiation and propagation was also studied through the previous laboratory experiments using acoustic emission (AE) monitoring [71,72].However, the change in the unloading rate did not affect the total number of tensile cracks and the shape of the curve.The shear cracks increased linearly when the unloading rate decreaed.Since the volume change rate of the sample was negative when the unloading rate was small enough, although the sample was completely destroyed, the number of micro-cracks continued to increase with the axis loading.Microscopically, as the unloading rate decreased, the number of the shear micro-cracks in the model increased, while the number of tensile micro-cracks did not change significantly.For wing cracks, the unloading rate had a significant effect on the inclination angle.The behavior was as follows: As the unloading rate increased, the inclination of the wing crack increased, which was more obvious in Type A and Type B. Unlike the initial confining pressure, the unloading rate also had a significant effect on the length of the wing crack.The length of the wing crack in each model increased as the unloading rate increased.For the cracks in the rock ligament, the unloading rate and the initial confining pressure showed similar effects: As the unloading rate decreased, the number of the shear microcracks in the ligament was significantly increased.Because for the specimens at the same strain state, the reduction in the unloading rate was somewhat equivalent to the increase in confining pressure.For the Type B model, in the process of reducing the unloading rate, not only did the number of shear cracks in the ligament increase, but also the cracks in the ligament became the main cracks of the sample.Different from the influence of initial confining pressure, in terms of anti-wing crack, the unloading rate had little influence on the angle of the crack, but the effect was mainly as follows: As the unloading rate decreased, the shear microcrack number rose.Macroscopically, the reduction of the unloading rate limited the development of tensile cracks, such as wing cracks and anti-wing cracks, which made the shear cracks develop space, and made the macro tensile cracks change into shear crack.Similar phenomena can be observed in the previous experimental observations [72].
Figure 19 shows the tensile crack and shear crack number curves for the Type A model under unloading conditions with various unloading rates.It can be seen from the figure that as the unloading rate decreased, the generation process of micro-cracks was delayed.The effect of the unloading rate on the crack initiation and propagation was also studied through the previous laboratory experiments using acoustic emission (AE) monitoring [71,72].However, the change in the unloading rate did not affect the total number of tensile cracks and the shape of the curve.The shear cracks increased linearly when the unloading rate decreaed.Since the volume change rate of the sample was negative when the unloading rate was small enough, although the sample was completely destroyed, the number of micro-cracks continued to increase with the axis loading.Microscopically, as the unloading rate decreased, the number of the shear micro-cracks in the model increased, while the number of tensile micro-cracks did not change significantly.Figure 20 presents deviatoric stress-axial strain curves of the pre-fissured marble specimens under unloading conditions with various unloading rates.From the curve characteristics, the model's deviatoric stress growth rate was significantly different with various unloading rates after entering the elastic phase: The faster the unloading rate, the greater the deviatoric stress growth rate.It is worth noting that the growth rate of the deviatoric stress curve at different unloading rates became uniform when the model entered the yielding phase, such as with Type C and D. Similar to the results of experiments under different initial confining pressures, only those samples with lower stress levels had obvious yielding stages with various unloading rates.For those models with low unloading rate and high stress level, the yielding stage and the failure stage roughly overlapped.The stress-strain curve of the defective specimen exhibited a sudden change in slope at a low confining pressure, which is consistent with the phenomena in previous studies [14,39].In terms of strength, as the unloading rate decreased, the ultimate strength of the sample also increased.The lower the unloading rate, the more significant the effect on the strength.Due to the mismatch between the confining pressure unloading rate and the axial loading rate, the volumetric compression ratio of the sample tends to be greater or less than one.The result is that the model with large unloading rate exhibited brittle failure similar to that under uniaxial compression in the failure stage due to the loss of confining pressure to deformation.The model with a small unloading rate did not fall even if it is completely destroyed.This can be seen from the graph in which the unloading rate is zero.Figure 20 presents deviatoric stress-axial strain curves of the pre-fissured marble specimens under unloading conditions with various unloading rates.From the curve characteristics, the model's deviatoric stress growth rate was significantly different with various unloading rates after entering the elastic phase: The faster the unloading rate, the greater the deviatoric stress growth rate.It is worth noting that the growth rate of the deviatoric stress curve at different unloading rates became uniform when the model entered the yielding phase, such as with Type C and D. Similar to the results of experiments under different initial confining pressures, only those samples with lower stress levels had obvious yielding stages with various unloading rates.For those models with low unloading rate and high stress level, the yielding stage and the failure stage roughly overlapped.The stress-strain curve of the defective specimen exhibited a sudden change in slope at a low confining pressure, which is consistent with the phenomena in previous studies [14,39].In terms of strength, as the unloading rate decreased, the ultimate strength of the sample also increased.The lower the unloading rate, the more significant the effect on the strength.Due to the mismatch between the confining pressure unloading rate and the axial loading rate, the volumetric compression ratio of the sample tends to be greater or less than one.The result is that the model with large unloading rate exhibited brittle failure similar to that under uniaxial compression in the failure stage due to the loss of confining pressure to deformation.The model with a small unloading rate did not fall even if it is completely destroyed.This can be seen from the graph in which the unloading rate is zero.Figure 20 presents deviatoric stress-axial strain curves of the pre-fissured marble specimens under unloading conditions with various unloading rates.From the curve characteristics, the model's deviatoric stress growth rate was significantly different with various unloading rates after entering the elastic phase: The faster the unloading rate, the greater the deviatoric stress growth rate.It is worth noting that the growth rate of the deviatoric stress curve at different unloading rates became uniform when the model entered the yielding phase, such as with Type C and D. Similar to the results of experiments under different initial confining pressures, only those samples with lower stress levels had obvious yielding stages with various unloading rates.For those models with low unloading rate and high stress level, the yielding stage and the failure stage roughly overlapped.The stress-strain curve of the defective specimen exhibited a sudden change in slope at a low confining pressure, which is consistent with the phenomena in previous studies [14,39].In terms of strength, as the unloading rate decreased, the ultimate strength of the sample also increased.The lower the unloading rate, the more significant the effect on the strength.Due to the mismatch between the confining pressure unloading rate and the axial loading rate, the volumetric compression ratio of the sample tends to be greater or less than one.The result is that the model with large unloading rate exhibited brittle failure similar to that under uniaxial compression in the failure stage due to the loss of confining pressure to deformation.The model with a small unloading rate did not fall even if it is completely destroyed.This can be seen from the graph in which the unloading rate is zero.Figure 21 shows the peak strength and equivalent elastic modulus of each specimens under unloading conditions with various unloading rates.It can be seen from the figure that under the same initial confining pressure, the peak strength of the pre-fissured specimen increased as the unloading rate decreased.Moreover, the higher the strength of the sample and the lower the unloading rate, the greater the increase in peak intensity.In terms of equivalent elastic modulus, the unloading rate is shown on some models: The greater the unloading rate, the smaller the equivalent modulus, especially in Type-B and Type-G.In the same period of time, the lower the unloading rate, the higher the confining pressure level, and the lateral deformation is limited.At this condition, the compression Figure 21 shows the peak strength and equivalent elastic modulus of each specimens under unloading conditions with various unloading rates.It can be seen from the figure that under the same initial confining pressure, the peak strength of the pre-fissured specimen increased as the unloading rate decreased.Moreover, the higher the strength of the sample and the lower the unloading rate, the greater the increase in peak intensity.In terms of equivalent elastic modulus, the unloading rate is shown on some models: The greater the unloading rate, the smaller the equivalent modulus, especially in Type-B and Type-G.In the same period of time, the lower the unloading rate, the higher the confining pressure level, and the lateral deformation is limited.At this condition, the compression behavior of the sample changed from axial compression to volume compression, and the equivalent modulus changed from elastic modulus to bulk modulus.Generally, the bulk modulus of rock is between 60%-80% of the elastic modulus, so the corresponding phenomenon of unloading rate and equivalent modulus can be reasonably explained.behavior of the sample changed from axial compression to volume compression, and the equivalent modulus changed from elastic modulus to bulk modulus.Generally, the bulk modulus of rock is between 60%-80% of the elastic modulus, so the corresponding phenomenon of unloading rate and equivalent modulus can be reasonably explained.

Conclusions
Three-dimensional (3-D) progressive mechanical and failure characteristics of pre-fissured marble specimens under the conventional triaxial unloading conditions were explored using 3-D DEM.The microscopic parameters were calibrated through comparison with the existing experimental data.Eight different types of pre-fissured rock specimens with different fissure distributions in the conventional unloading numerical tests were simulated.Effects of initial confining pressures and unloading rates on the mechanical and characteristics of pre-fissured rock specimens under the unloading conditions are numerically investigated.Some rules are drawn as follows: 1 Effects of fissure length 2a and ligament length b on the conventional peak strength is less than that of fissure inclination angle α and ligament inclination angle β Influences of pre-existing fissure length a, which is agreement with the analogous experimental results [20,69].The increase of initial confining pressures leads to the increasing number of shear cracks, which results in mixed tensile-shear and shear crack coalescence in the rock bridge region.
3 Initial confining pressure in the conventional triaxial unloading simulations has little effect on the equivalent elastic modulus.Meanwhile, the equivalent elastic modulus is significantly affected by the unloading rates.

Figure 1 .
Figure 1.Illustration of failure mechanism of: (a) Rock masses under unloading conditions; (b) stress path of unloading experiment.

Figure 1 .
Figure 1.Illustration of failure mechanism of: (a) Rock masses under unloading conditions; (b) stress path of unloading experiment.

Figure 3 .
Figure 3. Schematics of intact rock specimen under (a) uniaxial compression and (b) triaxial compression.(c) Geometry of pre-fissured marble specimens.

Figure 3 .
Figure 3. Schematics of intact rock specimen under (a) uniaxial compression and (b) triaxial compression.(c) Geometry of pre-fissured marble specimens.

Figure 5 .
Figure 5.Comparison of axial stress-strain curves obtained from the numerical simulations and the previous experimental data [65].

Figure 7 .
Figure 7. 3-D final failure patterns of the pre-fissured marble specimens in the uniaxial compressive simulations.

Figure 7 .
Figure 7. 3-D final failure patterns of the pre-fissured marble specimens in the uniaxial compressive simulations.

Figure 8 .
Figure 8. 3-D final failure patterns of the pre-fissured marble specimens in the conventional triaxial compressive simulations.

Figure 8 .
Figure 8. 3-D final failure patterns of the pre-fissured marble specimens in the conventional triaxial compressive simulations.

Figure 9 . 3 -
Figure 9. 3-D final failure patterns of the pre-fissured marble specimens in the confining pressure unloading simulations.

Figure 9 .
Figure 9. 3-D final failure patterns of the pre-fissured marble specimens in the confining pressure unloading simulations.

Figure 10 .
Figure10.Numerically-predicted deviatoric stress-axial strain curves and micro-crack evolution curves of the pre-fissured marble specimens under loading and unloading conditions.

Figure 10 .
Figure10.Numerically-predicted deviatoric stress-axial strain curves and micro-crack evolution curves of the pre-fissured marble specimens under loading and unloading conditions.

Figure 11 .Figure 12 .
Figure 11.Five different kinds of crack appear in the present numerical simulations: (a) Tensile cracks and (b) shear cracks.

Figure 11 .Figure 11 .Figure 12 .
Figure 11.Five different kinds of crack appear in the present numerical simulations: (a) Tensile cracks and (b) shear cracks.

Figure 12 .
Figure 12.Three modes of crack coalescence in the rock bridge region.(a) Tensile crack coalescence, (b) mixed tensile-shear crack coalescence, and (c) shear crack coalescence.

Figure 13 .
Figure 13.Influences of (a) crack length, (b) crack inclination angle, (c) ligament length, and (d) ligament inclination angle on the numerically predicted peak values of axial deviatoric stresses.

Figure 13 .
Figure 13.Influences of (a) crack length, (b) crack inclination angle, (c) ligament length, and (d) ligament inclination angle on the numerically predicted peak values of axial deviatoric stresses.

Figure 15 .
Figure15.The numerically-predicted axial deviatoric stress-axial strain curves of different marble specimens under unloading conditions with various initial confining pressures.

Figure 16 .
Figure 16.Influence of initial confining pressures on the number of microscopic cracks.

Figure 16 .
Figure 16.Influence of initial confining pressures on the number of microscopic cracks.

Figure 17 .
Figure 17.Effect of initial confining pressures on (a) peak strength and (b) elastic modulus.Figure 17.Effect of initial confining pressures on (a) peak strength and (b) elastic modulus.

Figure 17 .
Figure 17.Effect of initial confining pressures on (a) peak strength and (b) elastic modulus.Figure 17.Effect of initial confining pressures on (a) peak strength and (b) elastic modulus.

Figure 18 Figure 18 .
Figure 18 illustrates the 3-D final failure patterns of Type-A ~Type-D marble specimens under unloading conditions with various unloading rates.Compared with the initial confining pressure,

Figure 19 .
Figure 19.Influence of unloading rates of confining pressures on number of microscopic cracks.

Figure 19 .
Figure 19.Influence of unloading rates of confining pressures on number of microscopic cracks.

Figure 19 .
Figure 19.Influence of unloading rates of confining pressures on number of microscopic cracks.

Figure 20 .
Figure 20.The numerically-predicted axial deviatoric stress-axial strain curves of pre-fissured marble specimens under unloading conditions with different rates.

Figure 21 .
Figure 21.Effect of unloading rate of confining pressures on the peak strength and elastic modulus of marble.

Figure 21 .
Figure 21.Effect of unloading rate of confining pressures on the peak strength and elastic modulus of marble.

4
As the unloading rates increases in the conventional triaxial unloading simulations, the associated equivalent elastic modulus decreases, which indicates that bulk modulus of marble is smaller than the elastic modulus Author Contributions: Y.Z.: Conceptualization, validation, investigation, resources, data curation, project administration, funding acquisition; S.L.: software, validation, writing-original draft preparation; M.K.: methodology, validation, formal analysis, writing-review and editing; Z.W.: Conceptualization; validation; visualization; supervision.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by National Natural Science Foundation of China, grant number of 51674151.

Table 1 .
Microscopic mechanical parameters of marbles in the discrete element simulations.

Table 1 .
Microscopic mechanical parameters of marbles in the discrete element simulations.

Table 2 .
Geometric and loading parameters of pre-fissured marble specimens.

Table 3 .
Comparison of macro-mechanical properties of marbles.