Water Infusion on the Stability of Coal Specimen under Di ﬀ erent Static Stress Conditions

: Underground coal mines are frequently subjected to water infusion, resulting in many mining hazards. This study investigated the e ﬀ ect of water infusion on the stress and energy evolution characteristics of coal specimens representing isolated pillars under di ﬀ erent initial axial stress conditions using the discrete element method. A water infusion distribution model was developed, in which random functions were employed to describe water distribution for the purpose of realizing the dispersion of results for a better reliability. Based on the results, a stress-level classiﬁcation was presented to evaluate the water e ﬀ ect on pillars’ instability. For the investigated coal specimens, the water weakening e ﬀ ect on stress and energy remains stable when the axial geo-stress on pillars is less than 65% of uniaxial compressive strength (UCS). In contrast, when the axial stress coe ﬃ cient is greater than 65%, pillars become unstable eventually. A higher axial stress coe ﬃ cient is more likely to introduce a lower critical instability point of the water saturation coe ﬃ cient for pillars in the process of water infusion. However, the instability point remains random to some extent for specimens following the same water distribution rule under the identical test condition. Two instability types, which also happened randomly, were observed in the numerical results for damaged coal specimens under di ﬀ erent water saturation coe ﬃ cients and axial geo-stresses, namely free-falling and step-falling. rate v s increased exponentially with k s increasing. Therefore, the stress decreasing rate is believed to be a strong indicator for predicting pillar instability induced by water infusion. In practical engineering, the critical axial stress decreasing rate, which is about 8.3 × 10 −2 MPa per 1% of S c in this study, can be determined according to laboratory experiments and numerical simulations. For a pillar with stress sensors mounted inside it, its stability can be assessed and predicted by comparing its current stress decreasing rate and the critical value, which is of great significance for engineering.


Introduction
In underground coal mines, the stability of coal pillars and the surrounding rock is a key factor controlling the stability of the cavity group system. Groundwater, as an important environment factor, significantly weakens the coal strength and stiffness, leading further to a series of mining disasters, such as gob collapse, pillar instability, roof fall and surface subsidence [1][2][3][4][5][6][7][8][9][10]. When mines are abandoned, rock pillars are commonly in moisture conditions because of the rising water table. In addition, water injection is found to be effective for the mitigation and prevention of coal bump and rockburst [11][12][13][14][15][16]. Therefore, it is worthwhile to study the mechanical behavior of coal rocks under watery conditions.
Massive efforts have been made to explore the effect of water on the mechanical properties of coal [17][18][19][20][21][22][23][24][25][26]. For example, White and Mazurkiewicz [17] conducted uniaxial compressive tests on Nemo coal specimens with different water contents. They found that with the increase of water content, the compressive strength and Young's modulus of the Nemo coal specimen decrease. Yao et al. [22] carried out similar tests on Xishahe coal specimens. They discovered that the peak stress and modulus during elastic, strain softening and post-peak regimes decrease as water content increases while the failure strain increases. Based on the equivalent strain hypothesis and statistical damage theory, they derived a constitutive model for the coal specimens with different water contents [19]. Yao et al. [20] obtained the damage evolution characteristics of coal specimens by means of acoustic emission technique.
They found that the stress thresholds of coal are not affected by water, and they also stated that the presence of water promotes the shear failure in coal transformed from the split failure. Gu et al. [26] determined the dynamic compressive strength of coal specimens under moisture conditions by using a split Hopkinson pressure bar. They suggested that both of water content and porosity strongly affect the dynamic strength of coal.
In the majority of the above-mentioned studies, coal specimens were loaded under an increasing stress. However, the geo-stress acting on a pillar remains constant for a long time. The instability of a pillar in a stable stress field is essentially the result of the development and coalescence of inside micro-cracks, which is more likely to be creep deformation when the roof has a lower stiffness and to be stress relaxation when the pillar stands under a hard roof. By using some novel algorithms like the neural network and genetic algorithm, researchers are able to provide optimal estimation of the rheological parameters, which contributes to the prediction of the long-term rheological behavior of rock engineering [27]. As the degradation of the pillar is strongly sensitive to groundwater saturation, studying water's effect on rock and coal specimens is of great significance for a better understanding of the stability of pillars suffering from moisture or water immersion under geo-stress. However, this is impossible to determine by obtaining the whole instability process of pillars in underground mines. On the one hand, the instability of underground pillars is caused by combined factors including not only water and geo-stress but other complicated factors like defects, dynamic disturbance, weathering and so on, which makes it impracticable to quantitatively determine the water effect on pillar stability. On the other hand, the instability of pillars is unpredictable to some extent, which usually takes a long time, from months to years. Moreover, it is difficult to determine water distribution and its effects on the stability of rock engineering using the current monitoring techniques like microseismic monitoring. Coupled with the dangerous underground environment and limitation of monitoring techniques, it has made it difficult to reveal the realistic instability mechanism of water-effected pillars underground. In addition, laboratory experiments for such rheological mechanical behavior can be also time-consuming and expensive and fail to obtain some valuable parameters due to technical limitations. Taking the above-mentioned factors into consideration, numerical simulation seems to be a preferable approach to revealing the instability mechanism of pillar under water infusion at different geo-stress levels.
The purpose of this study is to investigate the effect of water infusion on the mechanical behavior of coal subjected to different initial stress. The mechanical conditions of an isolated pillar with hard roof and floor were simplified as a specimen sandwiched between two rigid platforms in the simulation section. The characteristics of the stress and energy release of coal pillars under different water saturation coefficients and initial geo-stresses were obtained and compared. This study is expected to give a better understanding of pillar stability under the influence of water and present an assessment and prediction method for the instability of hydrous pillar based on stress condition and water saturation coefficient.

Numerical Model in PFC
In this study, the two-dimensional particle flow code (PFC2D) based on the discrete element method (DEM) method was used to realize the simulation of mechanic behavior of rock specimens under water infusion. In PFC2D, the numerical specimen is composed of an assembly of rigid particles with different sizes. These particles are linked with contacts to represent materials like rock mass by means of an internal force and moment. The movements of particles are computed based on Newton's laws of motion. Here the parallel bond model (PBM) developed by Potyondy and Cundall [28] was employed for its excellent performance of reproducing the realistic mechanical behavior of rock-like materials, which has been proved by many studies [29][30][31][32][33].

Definition of Saturation Coefficient in Simulation
In numerous studies, saturation coefficient (S c ) is used to quantify the amount of water in a rock specimen. It is defined as the ratio of the current water content (w s ) of a rock specimen to the maximum water content (w m ): The saturation coefficient evidently ranges from 0 for absolutely dry samples to 1 for completely saturated samples.
Previous research suggested that the water-induced reduction in the strength of rock materials is primarily attributed to the physical and chemical deterioration of cement material between grains [34][35][36][37]. Therefore, we realized the water-weakening effect on coal behavior under the framework of the particle flow code according to the approach proposed by Jiang et al. [38,39], in which the bond strength and stiffness of the contact between two particles are reduced. In the numerical simulation, we set potential water-weakening contacts (PWCs) randomly distributed in the model. When the sample is in dry state, the PWCs have the same properties as the normal contacts. However, with the increase of saturation coefficient, the PWCs will be activated and converted into weakening contacts increasingly. Herein, the current water content (w s ) and the maximum water content (w m ) are defined as: where N w is the number of the activated weakening contacts; N p is the number of the potential water-weakening contacts, and N a is the total number of contacts in the model. Thus, the saturation coefficient in simulation model can be expressed as: Zhou et al. [40] studied the water distribution in cylindrical rock sample during soaking process by nuclear magnetic resonance technique. They reported that water permeates from the end surface into the core and the uneven water distribution significantly affects rock behavior. Therefore, the water distribution should be taken into consideration. For rock pillars exposed to underground water or moist atmosphere with low-permeability roof and floor, water mainly infuse into pillars from their side surface. To study water effect in such cases, the top and bottom surfaces of specimens were sealed to make sure water only infuses into the core from the side surfaces in both numerical and laboratory tests.
As shown in Equation (4), in the PFC2D model, the saturation coefficient of the model can be determined as: where s i = n wi/n pi is the saturation coefficient in the ith section line, n wi and n pi is the number of activated weakening contacts and potential water-weakening contacts in the ith section line, According to water distribution characteristics described by Zhou et al. [40], we assume that (1) for the dry state (Sc = 0), the saturation coefficient in each vertical line is 0. For the fully saturated state (Sc = 1), the saturation coefficient in each vertical line is 1; (2) for the partially saturated state, decreases monotonically from the surface to the core; (3) water distribution becomes more and more uniform with the increase of overall water saturation coefficient. Hence, the relation between and can be expressed as: where m is a coefficient changing with the overall water saturation coefficient (Sc) of the model. For a given Sc, m is a constant. After substituting Equation (5) into Equation (4), we obtain the expression of the overall water saturation coefficient for a rectangular model in PFC 2D as: The variation in si versus Di in different Sc is shown in Figure 2.  According to water distribution characteristics described by Zhou et al. [40], we assume that (1) for the dry state (S c = 0), the saturation coefficient in each vertical line is 0. For the fully saturated state (S c = 1), the saturation coefficient in each vertical line is 1; (2) for the partially saturated state, s i decreases monotonically from the surface to the core; (3) water distribution becomes more and more uniform with the increase of overall water saturation coefficient. Hence, the relation between S i and D i can be expressed as: where m is a coefficient changing with the overall water saturation coefficient (S c ) of the model. For a given S c , m is a constant. After substituting Equation (5) into Equation (4), we obtain the expression of the overall water saturation coefficient for a rectangular model in PFC 2D as: The variation in s i versus D i in different S c is shown in Figure 2. According to water distribution characteristics described by Zhou et al. [40], we assume that (1) for the dry state (Sc = 0), the saturation coefficient in each vertical line is 0. For the fully saturated state (Sc = 1), the saturation coefficient in each vertical line is 1; (2) for the partially saturated state, decreases monotonically from the surface to the core; (3) water distribution becomes more and more uniform with the increase of overall water saturation coefficient. Hence, the relation between and can be expressed as: where m is a coefficient changing with the overall water saturation coefficient (Sc) of the model. For a given Sc, m is a constant. After substituting Equation (5) into Equation (4), we obtain the expression of the overall water saturation coefficient for a rectangular model in PFC 2D as: The variation in si versus Di in different Sc is shown in Figure 2.  In order to validate the accuracy of the water distribution model, several cubic models were established with different saturation coefficients. The randomness function was employed, each modelling was repeated three times (cases 1 to 3). The simulated saturation coefficients and theoretical values are given in Table 1. It can be seen that the simulated saturation coefficients for all models are almost identical to the theoretical values with quite low average deviations from 0.10% to 3.00%. The water distribution also strictly follows the theoretical rule. For example, as shown in Figure 3, when the overall water saturation coefficient is 0.3, the s i obtained from the PFC model strongly agrees well with the theoretical curve calculated by Equation (4). Table 1. Saturation coefficients obtained in the particle flow code (PFC) model. In order to validate the accuracy of the water distribution model, several cubic models were established with different saturation coefficients. The randomness function was employed, each modelling was repeated three times (cases 1 to 3). The simulated saturation coefficients and theoretical values are given in Table 1. It can be seen that the simulated saturation coefficients for all models are almost identical to the theoretical values with quite low average deviations from 0.10% to 3.00%. The water distribution also strictly follows the theoretical rule. For example, as shown in Figure 3, when the overall water saturation coefficient is 0.3, the si obtained from the PFC model strongly agrees well with the theoretical curve calculated by Equation (4). Table 1. Saturation coefficients obtained in the particle flow code (PFC) model.

Model Setup
The numerical model was established as shown in Figure 4. The axial stress coefficient ks is defined as ks = σ/UCS, where σ is the current axial stress. Specimens were assumed to be weakened by water infusing from surface to core to simulate the mechanical weakening behavior of pre-stressed pillar under water infusion.

Model Setup
The numerical model was established as shown in Figure 4. The axial stress coefficient k s is defined as k s = σ/UCS, where σ is the current axial stress. Specimens were assumed to be weakened by water infusing from surface to core to simulate the mechanical weakening behavior of pre-stressed pillar under water infusion.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 15 In order to validate the accuracy of the water distribution model, several cubic models were established with different saturation coefficients. The randomness function was employed, each modelling was repeated three times (cases 1 to 3). The simulated saturation coefficients and theoretical values are given in Table 1. It can be seen that the simulated saturation coefficients for all models are almost identical to the theoretical values with quite low average deviations from 0.10% to 3.00%. The water distribution also strictly follows the theoretical rule. For example, as shown in Figure 3, when the overall water saturation coefficient is 0.3, the si obtained from the PFC model strongly agrees well with the theoretical curve calculated by Equation (4). Table 1. Saturation coefficients obtained in the particle flow code (PFC) model.

Model Setup
The numerical model was established as shown in Figure 4. The axial stress coefficient ks is defined as ks = σ/UCS, where σ is the current axial stress. Specimens were assumed to be weakened by water infusing from surface to core to simulate the mechanical weakening behavior of pre-stressed pillar under water infusion.  From Figure 5, the simulation procedure involves the following procedures: Step 1: Build the model with the dimension of 50 × 100 mm. We set PWC contacts randomly distributed in the model to simulate the cement in coal. The maximum water content of the coal specimen is 0.17, and the total number of contacts in the model is 54,945. Accordingly, the number of PWC is 9341; Step 2: The numerical model is compressed until its axial stress σ reaches the desired value for a given k s ; Step 3: Fix the loading plate. Increase the water saturation coefficient of the model with a constant rate of 0.01/step by activating the PWC according to the correlation between the S c and D i ; Step 4: The calculation is performed until the stress in the model begins to stay constant, which is regarded as a sign that the specimen has reached mechanical stability. However, if the specimen experiences overall instability during calculation, the calculation will be terminated and the simulation is completed; Step 5: Repeat Step 3 and 4 continually until the model either experiences instability or reaches its maximum saturation coefficient.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 15 From Figure 5, the simulation procedure involves the following procedures: Step 1: Build the model with the dimension of 50 × 100 mm. We set PWC contacts randomly distributed in the model to simulate the cement in coal. The maximum water content of the coal specimen is 0.17, and the total number of contacts in the model is 54,945. Accordingly, the number of PWC is 9341; Step 2: The numerical model is compressed until its axial stress σ reaches the desired value for a given ks; Step 3: Fix the loading plate. Increase the water saturation coefficient of the model with a constant rate of 0.01/step by activating the PWC according to the correlation between the Sc and Di; Step 4: The calculation is performed until the stress in the model begins to stay constant, which is regarded as a sign that the specimen has reached mechanical stability. However, if the specimen experiences overall instability during calculation, the calculation will be terminated and the simulation is completed; Step 5: Repeat Step 3 and 4 continually until the model either experiences instability or reaches its maximum saturation coefficient.

Calibration of Particle Parameters
Standard uniaxial compressive laboratory experiments and numerical simulations were carried out in order to calibrate the optimal numerical parameters of the model. Figure 6 shows the failure patterns of dry and saturated coal specimens in laboratory tests and their corresponding stress-strain curves. Table 2 lists the critical mechanical properties of the coal specimens. It can be seen that, the

Calibration of Particle Parameters
Standard uniaxial compressive laboratory experiments and numerical simulations were carried out in order to calibrate the optimal numerical parameters of the model. Figure 6 shows the failure patterns of dry and saturated coal specimens in laboratory tests and their corresponding stress-strain curves. Table 2 lists the critical mechanical properties of the coal specimens. It can be seen that, the UCS and young's modulus of the fully saturated coal specimen reduced by 69.4% and 51.2%, respectively. UCS and young's modulus of the fully saturated coal specimen reduced by 69.4% and 51.2%, respectively.  Sufficient trial and error tests were performed through numerical static uniaxial compressive tests in PFC model to calibrate the microscopic parameters, such that the macroscopic mechanical behavior of the model matches that of the coal specimens well. We found that for the parallel bond contact model, there is a near-linear relation between macro-mechanical properties and relevant microscopic parameters in terms of strength and deformation within given limits. The microscopic strength and deformation parameters are also nearly independent to each other. Therefore, the rough range of the microscopic parameters for strength and deformation can be determined separately using linear interpolation method and then a much accurate range is determined for the next calibration. The optimal parameters will be determined after several times of iteration. The final calibration parameters of the normal contacts are listed in Table 3. Similarly, the "trial and error" method was also employed to calculate the microscopic parameters for the weakening contacts based on the laboratory result of the fully saturated specimen. Finally, it was determined that the deformation effective modulus and bond deformation effective modulus are only 6% of those of the normal contacts. Additionally, the cohesion and tensile strength are 9% of those of the normal contacts. Other parameters are not changed. Note: μ: particle friction coefficient; Rmax: maximum particle radius; Rmin: maximum particle radius; ρ: particle density; Ec: deformability effective modulus; k * : normal-to-shear stiffness ratio; ̅ : bond deformability effective modulus; ̅ * : bond normal-to-shear stiffness ratio; ̅ : cohesion; ̅t: tensile strength.  Sufficient trial and error tests were performed through numerical static uniaxial compressive tests in PFC model to calibrate the microscopic parameters, such that the macroscopic mechanical behavior of the model matches that of the coal specimens well. We found that for the parallel bond contact model, there is a near-linear relation between macro-mechanical properties and relevant microscopic parameters in terms of strength and deformation within given limits. The microscopic strength and deformation parameters are also nearly independent to each other. Therefore, the rough range of the microscopic parameters for strength and deformation can be determined separately using linear interpolation method and then a much accurate range is determined for the next calibration. The optimal parameters will be determined after several times of iteration. The final calibration parameters of the normal contacts are listed in Table 3. Similarly, the "trial and error" method was also employed to calculate the microscopic parameters for the weakening contacts based on the laboratory result of the fully saturated specimen. Finally, it was determined that the deformation effective modulus and bond deformation effective modulus are only 6% of those of the normal contacts. Additionally, the cohesion and tensile strength are 9% of those of the normal contacts. Other parameters are not changed. Note: µ: particle friction coefficient; R max : maximum particle radius; R min : maximum particle radius; ρ: particle density; E c : deformability effective modulus; k * : normal-to-shear stiffness ratio; E c : bond deformability effective modulus; k * : bond normal-to-shear stiffness ratio; c: cohesion; σ t : tensile strength.

Results
We repeated each case for five times to eliminate the data discreteness induced by the random-distributed PWC in the model. The simulated results are presented in three sub-sections: (1) stress variation, (2) energy evolution and (3) failure pattern. Figure 7 shows the stress variation versus saturation coefficient under different initial stress conditions. For the specimen with lower initial stress (ks < 0.65), its bearing stress decreases linearly until the specimen reached saturation in all cases, the specimen remains stable. When the instability occurs, the stress in the specimen suddenly drops. Generally, two stress falling modes, namely free-falling and step-falling, can be concluded for those damaged specimens. For specimens suffered free-falling, the stress decreases linearly at first then dropped rapidly to a very low level close to 0 when the saturation coefficient reached a certain value, indicating that the instability mode is violent and drastic. For specimens experienced step-fall instability, the stress drops several times in a relatively mild manner.

Stress Variation
We repeated each case for five times to eliminate the data discreteness induced by the randomdistributed PWC in the model. The simulated results are presented in three sub-sections: (1) stress variation, (2) energy evolution and (3) failure pattern. Figure 7 shows the stress variation versus saturation coefficient under different initial stress conditions. For the specimen with lower initial stress (ks < 0.65), its bearing stress decreases linearly until the specimen reached saturation in all cases, the specimen remains stable. When the instability occurs, the stress in the specimen suddenly drops. Generally, two stress falling modes, namely freefalling and step-falling, can be concluded for those damaged specimens. For specimens suffered freefalling, the stress decreases linearly at first then dropped rapidly to a very low level close to 0 when the saturation coefficient reached a certain value, indicating that the instability mode is violent and drastic. For specimens experienced step-fall instability, the stress drops several times in a relatively mild manner.

Stress Variation
It should be noted that there is a certain randomness for the instability mode. The failure evolution can be dominated by either free-fall instability or step-fall instability for different specimens under the same condition. For example, when ks is 0.75, specimens in Cases 1, 3 and 4 were considered to be characterized by step-fall instability where stress dropped and then remained stable for some further saturation increase until the specimens reached total instability. However, the stress dropped dramatically to zero once the water saturation coefficient reached a certain point in the other cases where specimens were characterized by free-fall instability. We can also find that the ks of 0.65 is the instability critical point in all cases. For specimens with ks less than 0.65, they all released the stress stably until they reached saturation. However, when ks was equal or greater than 0.65, the specimens collapsed with dramatic stress drop.  It should be noted that there is a certain randomness for the instability mode. The failure evolution can be dominated by either free-fall instability or step-fall instability for different specimens under the same condition. For example, when k s is 0.75, specimens in Cases 1, 3 and 4 were considered to be characterized by step-fall instability where stress dropped and then remained stable for some further saturation increase until the specimens reached total instability. However, the stress dropped dramatically to zero once the water saturation coefficient reached a certain point in the other cases where specimens were characterized by free-fall instability. We can also find that the k s of 0.65 is the instability critical point in all cases. For specimens with k s less than 0.65, they all released the stress stably until they reached saturation. However, when k s was equal or greater than 0.65, the specimens collapsed with dramatic stress drop. To understand the coupled effects of pre-axial stress and water saturation coefficient, we introduced the critical instability water saturation coefficient (S cr ), i.e., the water saturation coefficient at the stress drop point. Figure 8 shows that S cr decreases linearly with the increase of pre-stress when k s is not less than 0.65. For example, under the k s of 0.65, S cr is 0.86. When k s increases to 0.80, S cr significantly drops to 0.41. This suggests that under higher the axial stress, little water saturation coefficient will result in the instability of rock pillar.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 15 To understand the coupled effects of pre-axial stress and water saturation coefficient, we introduced the critical instability water saturation coefficient (Scr), i.e., the water saturation coefficient at the stress drop point. Figure 8 shows that Scr decreases linearly with the increase of pre-stress when ks is not less than 0.65. For example, under the ks of 0.65, Scr is 0.86. When ks increases to 0.80, Scr significantly drops to 0.41. This suggests that under higher the axial stress, little water saturation coefficient will result in the instability of rock pillar. To better describe the weakening effect of water on axial stress, here we define the stress decreasing rate vs as the decrement of axial stress when the water saturation coefficient increased 1% during the linear stable stress-reduction period. For example, the linear stable stress-reduction period for the specimen with ks = 0.80 in Case 1 was considered to be where the saturation coefficient increased from 0.00 to 0.35.
As mentioned above, when ks was greater than 0.65, the specimens finally reached instability with the increase of saturation coefficient. However, when ks was less than or equal to 0.65, the weakening rate increased gently and stably until the specimen was fully saturated. It can be seen from Figure 9 that the stress decreasing rate vs increased exponentially with ks increasing. Therefore, the stress decreasing rate is believed to be a strong indicator for predicting pillar instability induced by water infusion. In practical engineering, the critical axial stress decreasing rate, which is about 8.3 × 10 −2 MPa per 1% of Sc in this study, can be determined according to laboratory experiments and numerical simulations. For a pillar with stress sensors mounted inside it, its stability can be assessed and predicted by comparing its current stress decreasing rate and the critical value, which is of great significance for engineering.  To better describe the weakening effect of water on axial stress, here we define the stress decreasing rate v s as the decrement of axial stress when the water saturation coefficient increased 1% during the linear stable stress-reduction period. For example, the linear stable stress-reduction period for the specimen with k s = 0.80 in Case 1 was considered to be where the saturation coefficient increased from 0.00 to 0.35.
As mentioned above, when k s was greater than 0.65, the specimens finally reached instability with the increase of saturation coefficient. However, when k s was less than or equal to 0.65, the weakening rate increased gently and stably until the specimen was fully saturated. It can be seen from Figure 9 that the stress decreasing rate v s increased exponentially with k s increasing. Therefore, the stress decreasing rate is believed to be a strong indicator for predicting pillar instability induced by water infusion. In practical engineering, the critical axial stress decreasing rate, which is about 8.3 × 10 −2 MPa per 1% of S c in this study, can be determined according to laboratory experiments and numerical simulations. For a pillar with stress sensors mounted inside it, its stability can be assessed and predicted by comparing its current stress decreasing rate and the critical value, which is of great significance for engineering.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 15 To understand the coupled effects of pre-axial stress and water saturation coefficient, we introduced the critical instability water saturation coefficient (Scr), i.e., the water saturation coefficient at the stress drop point. Figure 8 shows that Scr decreases linearly with the increase of pre-stress when ks is not less than 0.65. For example, under the ks of 0.65, Scr is 0.86. When ks increases to 0.80, Scr significantly drops to 0.41. This suggests that under higher the axial stress, little water saturation coefficient will result in the instability of rock pillar. To better describe the weakening effect of water on axial stress, here we define the stress decreasing rate vs as the decrement of axial stress when the water saturation coefficient increased 1% during the linear stable stress-reduction period. For example, the linear stable stress-reduction period for the specimen with ks = 0.80 in Case 1 was considered to be where the saturation coefficient increased from 0.00 to 0.35.
As mentioned above, when ks was greater than 0.65, the specimens finally reached instability with the increase of saturation coefficient. However, when ks was less than or equal to 0.65, the weakening rate increased gently and stably until the specimen was fully saturated. It can be seen from Figure 9 that the stress decreasing rate vs increased exponentially with ks increasing. Therefore, the stress decreasing rate is believed to be a strong indicator for predicting pillar instability induced by water infusion. In practical engineering, the critical axial stress decreasing rate, which is about 8.3 × 10 −2 MPa per 1% of Sc in this study, can be determined according to laboratory experiments and numerical simulations. For a pillar with stress sensors mounted inside it, its stability can be assessed and predicted by comparing its current stress decreasing rate and the critical value, which is of great significance for engineering.

Energy Evolution
The evolution curves of strain energy per unit volume (e) versus water saturation coefficient under different initial stress conditions are presented in Figure 10. It shows that the evolution of strain energy per unit volume agrees quite well with that of axial stress, which indicates that strain energy stored in a specimen is associated with its stress condition.

Energy Evolution
The evolution curves of strain energy per unit volume (e) versus water saturation coefficient under different initial stress conditions are presented in Figure 10. It shows that the evolution of strain energy per unit volume agrees quite well with that of axial stress, which indicates that strain energy stored in a specimen is associated with its stress condition. As similar to vs, here the strain energy releasing rate ve is defined as the rate of released strain energy per unit volume when the water saturation increased 1% during the linear stable strain energy-reduction period. As shown in Figure 11, there is a perfect exponential relationship between ve and ks. As the instability of a specimen is caused by the release of its stored strain energy, the strain energy releasing rate is a direct reflection of its instability propensity. The strain energy releasing rate remained at a low level when the axial stress coefficient was less or equal to 0.4, and then increased obviously. When the axial stress coefficient reached the critical value of 0.65, the strain energy releasing rate began to accelerate rapidly with the increase of water saturation coefficient. As similar to v s , here the strain energy releasing rate v e is defined as the rate of released strain energy per unit volume when the water saturation increased 1% during the linear stable strain energy-reduction period. As shown in Figure 11, there is a perfect exponential relationship between v e and k s . As the instability of a specimen is caused by the release of its stored strain energy, the strain energy releasing rate is a direct reflection of its instability propensity. The strain energy releasing rate remained at a low level when the axial stress coefficient was less or equal to 0.4, and then increased obviously. When the axial stress coefficient reached the critical value of 0.65, the strain energy releasing rate began to accelerate rapidly with the increase of water saturation coefficient.
From Figures 10 and 11 it can be seen that in all cases, the results of both v s and v e are almost identical when axial stress coefficient was under the critical value of 0.65 and diverging when axial stress coefficient was over the critical value. This demonstrates that the randomness of water distribution at the micro-level has little effect on the specimens' stress and energy characteristics under lower initial axial stress condition while this effect became obvious under higher initial axial stress condition. The critical stress decreasing rate mentioned-above can fluctuate within a wider range for an engineering case and therefore a higher safety factor is required to ensure the stability of the pillar.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 15 From Figures 10 and 11 it can be seen that in all cases, the results of both vs and ve are almost identical when axial stress coefficient was under the critical value of 0.65 and diverging when axial stress coefficient was over the critical value. This demonstrates that the randomness of water distribution at the micro-level has little effect on the specimens' stress and energy characteristics under lower initial axial stress condition while this effect became obvious under higher initial axial stress condition. The critical stress decreasing rate mentioned-above can fluctuate within a wider range for an engineering case and therefore a higher safety factor is required to ensure the stability of the pillar. Figure 11. ve evolution curves with ks increasing.

Failure Evolution
The final failure patterns for all specimens are shown in Figure 12. Similar failure patterns were observed in all specimens. In the majority of cases, the sides of specimens, which were water-rich areas, were damaged by splitting failure with fragments. Though cracks and fragments distribution varied from case to case, all the failure patterns were dominated by shear failure through the specimens. In general, specimens under higher initial axial stress got damaged more severely with more cracks and fragments. Comparing them with strain energy releasing rate (Figure 12), it can be indicated that failure intensity highly depends on the release of strain energy. A higher strain energy releasing rate means that more strain energy is released with the increase of water saturation coefficient, which may contribute to increasing the propensity of sudden instability or coal bump.
As shown in Figure 13, the failure evolution of the specimen with ks = 0.80 in Case 5 was chosen for detailed analysis of failure evolution. In the beginning, initial cracks appeared and concentrated on both the sides. Later on, cracks propagated into the core and several fragments were split from the specimen. When the water saturation coefficient reached 0.41, a shear macro-crack strip finally ran through the specimen and resulted in the overall instability. Figure 11. v e evolution curves with k s increasing.

Failure Evolution
The final failure patterns for all specimens are shown in Figure 12. Similar failure patterns were observed in all specimens. In the majority of cases, the sides of specimens, which were water-rich areas, were damaged by splitting failure with fragments. Though cracks and fragments distribution varied from case to case, all the failure patterns were dominated by shear failure through the specimens. In general, specimens under higher initial axial stress got damaged more severely with more cracks and fragments. Comparing them with strain energy releasing rate (Figure 12), it can be indicated that failure intensity highly depends on the release of strain energy. A higher strain energy releasing rate means that more strain energy is released with the increase of water saturation coefficient, which may contribute to increasing the propensity of sudden instability or coal bump.

Discussion
The axial stress and strain energy within the specimens are more sensitive to water saturation coefficient under a higher initial axial stress condition. With water infusion, three initial stress levels can be roughly classified for the investigated coal specimens as follows: (1) When ks is less than 0.45, i.e., the axial geo-stress is less than 40% of the UCS of rock pillar, water has limited effect on the stress and energy for coal at such a low stress level. The axial stress seems to slightly decrease with the increase of water saturation coefficient, indicating the coal As shown in Figure 13, the failure evolution of the specimen with k s = 0.80 in Case 5 was chosen for detailed analysis of failure evolution. In the beginning, initial cracks appeared and concentrated on both the sides. Later on, cracks propagated into the core and several fragments were split from the specimen. When the water saturation coefficient reached 0.41, a shear macro-crack strip finally ran through the specimen and resulted in the overall instability.

Discussion
The axial stress and strain energy within the specimens are more sensitive to water saturation coefficient under a higher initial axial stress condition. With water infusion, three initial stress levels can be roughly classified for the investigated coal specimens as follows:

Discussion
The axial stress and strain energy within the specimens are more sensitive to water saturation coefficient under a higher initial axial stress condition. With water infusion, three initial stress levels can be roughly classified for the investigated coal specimens as follows: (1) When k s is less than 0.45, i.e., the axial geo-stress is less than 40% of the UCS of rock pillar, water has limited effect on the stress and energy for coal at such a low stress level. The axial stress seems to slightly decrease with the increase of water saturation coefficient, indicating the coal pillar gradually degrades. However, the coal pillar keeps stable regardless of the water saturation coefficient. (2) When k s ranges from 0.40 to 0.65, though the specimens still remain stable despite of increasing water saturation coefficient, the strain energy releasing rate began increasing rapidly. In deep coal mining, rockburst accompanying fast energy release commonly happens in such a geo-stress level [41][42][43][44]. From Figure 10, water can release the stored energy within the specimen without triggering instability. Hence, water injection is often used for rockburst prevention [5,11,12]. (3) When k s exceeds 0.65, the stress drops abruptly. With the increase of k s , the critical instability water saturation coefficient decreases. This means that under such geo-stresses, the coal pillar is unstable and easier to be failed by water infusion.
It should be noted that rock and coal with different mechanical properties may have different sensitivities to axial stress coefficient under water infusion, which leads to different threshold values for the classification.

Conclusions
In the present study, the coupled effect of water infusion and axial stress on rock stability was numerically investigated. The main conclusions can be drawn as follows: (1) For all damaged specimens in the simulations, the decrease of stress and energy can be classified into two manners, namely free-falling and step-falling. For the former, the stress and energy of coal specimen decreases linearly and stably and then drops to zero suddenly. The catastrophic and violent instability is very dangerous without obvious early warning, such as coal bump and rockburst. For the latter, specimens experience several steps of stress and energy drop. (2) The stress decreasing rate (v s ) is suggested to be an effective index to assess the stability of pillar. There is a clear critical point of v s over which the pillar will enter into an unstable stage with a high instable risk. In practical engineering, when the v s of a pillar is remarkably greater than the other stable pillars, the objective one can be determined as unstable. In addition, exact assessment and prediction of the stability of pillars can be expected through field monitoring and the laboratory experiments. (3) It is known that reliable laboratory results should be repeated. For numerical simulations where randomness is taken into consideration, reiterations are also expected to be required. As shown in this study, the distribution of PWC is random at a micro-level though it follows a definite law. Consequently, the average results show a strong regularity of the effects of initial axial stress coefficient and water saturation coefficient on the stress and energy characteristics of the specimen. However, the results of the same simulation are different every time. As the water distribution on a minute scale is exactly different for different underground pillars even with the same mechanical properties and water saturation coefficient, the randomness adopted in this study is reasonable and agrees with the practical engineering. The regular average results of those repeated cases should be correct. Therefore, for the numerical simulation studying such rock engineering problems, randomness function and repeated cases are suggested for better outcomes.
Author Contributions: Z.Z. and L.T. conceived and designed the work; L.T. and X.C. proposed and participated in the numerical simulation; L.T. and X.C. wrote the paper. All authors have read and agreed to the published version of the manuscript.