Constitutive model for rock damage
The inherent intricate composition and structure of rock materials make them inherently susceptible entities35. Therefore, the utilization of damage mechanics and the development of appropriate damage constitutive models have emerged as effective strategies for investigating rock deformation. In terms of micro-damage mechanics, the representation of the number of internal damage microelements within rocks during the loading process is denoted as Nd. To quantitatively measure this parameter, a statistical damage variable D is introduced, which is defined as the ratio of the number of damaged microelements to the total number of microelements N.
$$D=\frac{{N}_{d}}{N}.$$
(24)
The failure of a microelement in rock occurs once its stress value reaches the fracture strength. Assuming a normal distribution for the strength of internal microelements in the rock, including the maximum stress values, the failure probability, denoted as p(σ), can be defined within a stress interval (σ, σ + dσ).
$$p\left(\sigma \right)=\frac{1}{\sqrt{2\pi }m}\text{exp}\left(-\frac{{\left(\sigma -n\right)}^{2}}{2{m}^{2}}\right),$$
(25)
where the function p(σ) represents the probability density function that describes the strength of the microelement, m denotes the mean value and n represents the standard deviation of the distribution.
When the external load reaches a stress level σ1, microelements within the rock that have a strength lower than this level will fail. The cumulative number of failed microelements, denoted as Nd, represents the sum of previously failed microelements in the preceding intervals.
$${N}_{d}={\int }_{-\infty }^{{\sigma }_{1}}Np\left(\sigma \right)d\sigma =NP\left({\sigma }_{1}\right),$$
(26)
with
$$P\left({\sigma }_{1}\right)=\frac{1}{\sqrt{2\pi }\sigma }{\int }_{-\infty }^{{\sigma }_{1}}\text{exp}\left(-\frac{{\left({\sigma }_{1}-n\right)}^{2}}{2{m}^{2}}\right)dt=\Phi \left(\frac{{\sigma }_{1}-n}{m}\right),$$
(27)
where P(σ1) represents the statistical probability function that describes the strength of microelements, and Φ represents the standard normal distribution function.
The statistical damage variable D can be given by:
$$D=P\left(\sigma \right).$$
(28)
The aforementioned equations, derived from statistical strength theory, illustrates the progression of rock damage. With increasing stress level of σ, the corresponding function value P(σ) also increases gradually. Additionally, the statistical probability function P(σ) ranges from 0 to 1, reflecting the variation of the damage variable (D).
The stress-strain relationship of brittle rocks follows a two-stage division, governed by the evolution law of rock damage36. Assuming ε0 represents the strain at the initiation of damage, the constitutive equation for rock damage, derived through the principle of strain equivalence, can be formulated as:
$$\upsigma =\left\{\begin{array}{c}E\varepsilon , 0\le \varepsilon \le {\varepsilon }_{0}\\ E\left(1-D\right)\varepsilon , \varepsilon >{\varepsilon }_{0}\end{array}\right.,$$
(29)
where the damage evolution model is given by:
$$\text{D}=\left\{\begin{array}{c}0, 0\le \varepsilon \le {\varepsilon }_{0}\\ \frac{1}{\sqrt{2\pi }\sigma }{\int }_{-\infty }^{\sigma }\text{exp}\left(-\frac{{\left(\sigma -n\right)}^{2}}{2{m}^{2}}\right)dt=\Phi \left(\frac{\sigma -n}{m}\right), \varepsilon >{\varepsilon }_{0}\end{array}\right..$$
(30)
Damage intensity index and PCF correction
The Griffith fracture theory is a fundamental explanation of material fracture behavior. It assumes that materials are homogeneous and isotropic, providing insights into crack initiation and suggesting that the critical stress for fracture is related to the initial crack length. However, this theory does not fully account for the direction of crack propagation. In reality, materials often exhibit microstructural variations, anisotropy, and multiple defects, which have a significant impact on fracture behavior. Moreover, the mode of crack propagation is influenced by various factors, including material properties, loading conditions, and stress concentration points37.
Crack propagation is an integral part of the rock failure process, and the evolution of rock damage can be divided into five distinct stages4,38. The rock damage evolution (blue curve) and the stress-strain relationship during the damage process (red curve) are presented in Fig. 5.
Figure 5
Rock damage evolution process and stress-strain curve.
(1) Compaction stage (oa): Under compressive loads, microvoids and cracks within rocks undergo closure as a result of the applied pressure.
(2) Linear elastic stage (ab): During this stage, rocks display elastic deformation when subjected to external forces. Elastic deformation is reversible, causing the rocks to revert to their original shape and volume once the external force is removed. Rocks exhibit high rigidity and strength in this stage.
(3) Plastic stage (bc): When the external force surpasses the strength limit of the rocks, plastic deformation takes place. In the plastic stage, rocks undergo volume contraction, yielding, and deformation without immediate fracture. The rocks lose the ability for complete elastic recovery in this stage.
(4) Fracture stage (cd): As the external force continues to increase past the fracture strength of the rocks, cracks and damage begin to manifest. In this stage, cracks inside the rocks gradually expand and connect, forming distinguishable fracture surfaces.
(5) Collapse stage (de): Once the fractures and cracks in the rock propagate to a certain extent, the overall strength of the rock mass decreases, leading to rock collapse and failure. In the collapse stage, the stability of the rocks is entirely compromised, resulting in the loss of cohesion among rock blocks and substantial destruction and deformation.
The evolution law of rock damage encompasses various stages, starting from the initiation of cracks to their propagation, leading to rock instability and failure. Sustaining continuous crack expansion requires increased energy input due to surface energy requirements. Nevertheless, the unpredictable nature of crack propagation may result in a CF exceeding the crack initiation stage.
To provide a more precise understanding of rock failure behavior, based on the research on the rock damage constitutive model, this paper incorporates the concept of a damage intensity index (Ke). This index allows for the adjustment of the cutting force (CF) at the onset of crack formation in rocks, considering the energy input during the rock failure process and the prevailing conditions during unstable failure.
The area under the stress-strain curve represents the total work performed by external forces on the rock, encompassing both elastic and plastic deformations (See Fig. 5). Considering the influence of rock brittleness on the initial crack size, the elastic energy normalization factor Kf related to rock brittleness is introduced in this paper based on Griffith’s theory, which is used to correct the fracture stress. The corrected fracture stress (σL’) is incorporated into the damage constitutive equation of the rock in order to compute the damage intensity index Ke.
The damage intensity index Ke is defined as the ratio between the energy needed for crack initiation and the total energy necessary for rock failure, σ is characterized by a damage evolution model that satisfies a normal distribution:
$$K_{e} = \frac{{\mathop \smallint \nolimits_{0}^{{\varepsilon_{e} }} e}}{{\mathop \smallint \nolimits_{e}^{{\varepsilon_{e} }} \sigma d\varepsilon }} = \frac{{\mathop \smallint \nolimits_{0}^{{\varepsilon_{C} }} E_{0} \left( {1 – D} \right)\varepsilon d\varepsilon }}{{\mathop \smallint \nolimits_{0}^{{\varepsilon_{L} }} E_{0} \left( {1 – D} \right)\varepsilon d\varepsilon }}.$$
(31)
In summary, the modified PCF can be represented by the following equation:
$$PCF = K_{e} \cdot \frac{{2K_{Ic} }}{{\sqrt {\pi \sqrt(3){{\frac{2}{{\pi^{2} }}\mathop \smallint \nolimits_{0}^{d} p_{l}^{2} \mathop \smallint \nolimits_{0}^{a} \cos^{2} \theta \cdot \sqrt {1 + \tan^{2} \theta_{2} } \cdot K_{f} dxdl}}} }}\mathop \smallint \limits_{0}^{d} \sin \varphi p_{l} \mathop \smallint \limits_{0}^{a} \cos \theta \sqrt {1 + \tan^{2} \theta_{2} } dxdl.$$
(32)
Based on the Griffith theory and rock damage evolution theory, there is a strong correlation between the rock brittleness index (C) and the elastic energy normalization factor (Kf). Taking the elastic modulus (E) and uniaxial compressive strength (UCS) of rocks as references, the relationship between C and Kf is fitted to 20 sets of experimental data of full-size rock cuts selected on the basis of the research results of Belgin. This relationship is then extended and applied to the correction of PCF.
A comprehensive overview of the crucial parameters pertaining to the rock samples, test conditions, and corresponding results observed across the selected 20 test cases is presented in Table 1. The tests are conducted using Sandvik S-35/80H conical pick with a tip angle of 80 degrees and an alloy head diameter of 20 mm. The cutting angle is set at 55 degrees, and the cutting speed is 12.7 cm/s14. In Table 1, E is the modulus of elasticity, UCS is the uniaxial compressive strength, BTS is the uniaxial tensile strength, KIc is the type I fracture toughness, εC is the critical strain, γ is the mounting angle, d is the cutting depth, PCF0 is the pre-correction peak cutting force, PCF is the corrected peak cutting force, PCFexp is the experimental value, Ke is the damage intensity index, Kf is the elastic energy normalization factor, C is the rock brittleness index.
Table 1 Rock characteristic parameters in cutting tests and calculation results of PCF.
A nonlinear equation system is constructed based on the UCS, εC, and the damage evolution model (See Eq. 30) of the rock samples’ parameters from the 20 cutting tests. By solving for the distribution coefficients (m and n) of the normal distribution, both the stress–strain curve during the rock damage process and the evolution curve of the damage variable (Fig. 5) are established. Through the integration of Eqs. (31, 32), and the obtained damage evolution model, the corrected PCF can be determined.
For the Griffith theory and rock damage evolution theory, Kf is only related to the rock type and independent of the cutting depth. Based on the strong correlation between Kf and C, there exists an optimal Kf and Ke that accurately reflect the rock’s failure characteristics under different rock brittleness conditions, leading to the precise correction of PCF. By solving Eqs. (14, 15, 22), an expression for fracture stress incorporating Kf can be derived.
$${{\upsigma }_{L}}{\prime}=\frac{2{K}_{Ic}}{\sqrt{\pi \sqrt(3){\frac{2}{{\pi }^{2}}{\int }_{0}^{d}{p}_{l}^{2}{\int }_{0}^{a}{cos}^{2}\theta \cdot \sqrt{1+{tan}^{2}{\theta }_{2}}\cdot {K}_{f}dxdl}}}.$$
(33)
The fracture stress (σL’) corrected in relation to Kf should not surpass the UCS of the rock. To ensure this, the Ke is determined within a suitable range of Kf for each rock sample, resulting in the correction of PCF. The theoretical PCF values derived from the 20 experimental sets showed the highest accuracy compared to the experimental PCF values. The corresponding data are meticulously recorded in Table 1, facilitating a comprehensive evaluation of the fitting accuracy.
The brittleness index considered in this study incorporates the elastic modulus, uniaxial compressive strength, and tensile strength as contributing factors in Eq. (34). A lower value of the brittleness index indicates a higher degree of brittleness in the rock. The relationship between C and Kf can be determined using Allometric function fitting (See Fig. 6).
$$C=\frac{{e}^\frac{E}{20000}\cdot {\sigma }_{c}-{\sigma }_{t}}{{e}^\frac{E}{20000}\cdot {\sigma }_{c}+{\sigma }_{t}},$$
(34)
where σt represents the tensile strength of the rock, σc refers to the uniaxial compressive strength of the rock.
Figure 6
Energy correction function fitted by C and Kf.
The elastic energy normalization factor can be expressed in terms of the brittleness index C as follows:
$${K}_{f}=\text{a}\cdot {C}^{b}=0.00747\cdot {\left(\frac{{e}^\frac{E}{20000}\cdot {\sigma }_{c}-{\sigma }_{t}}{{e}^\frac{E}{20000}\cdot {\sigma }_{c}+{\sigma }_{t}}\right)}^{-22.75684}.$$
(35)
Through the correlation analysis, the determination coefficient is 0.89, indicating a strong predictive capability of the independent variable C in the regression model for the dependent variable Kf. Furthermore, a two-tailed t-test is performed on the two parameters with a confidence level of 0.05, as presented in Table 2. The degree of freedom of the sample is 8 and the level of significance is 0.05. The critical value tL is determined to be 2.306 based on the t-test critical value table. Significantly, both the absolute values of the t-statistics for the two parameters exceed tL, while the corresponding probability (p-value) for the larger absolute t-value is considerably smaller than 0.025. These findings suggest that both parameters are crucial in the model, playing significant roles in interpreting and fitting the data. In Table 2, SE is the standard error.
Table 2 t-value test of energy correction function.
The energy correction function, obtained from the relationship between C and Kf, along with the determination coefficient and t-test results, clearly demonstrate the strong correlation exhibited by the C and the Kf. For rocks with known lithological parameters, the brittleness index can be solved to characterize the rock’s brittleness, and the elastic energy normalization factor Kf can be calculated through this function. In further analysis, in conjunction with the rock’s constitutive equation, the damage intensity index Ke can be used to modify the theoretical PCF. The solutions for C, Kf, and Ke for each rock sample are listed in Table 1.
To verify the reliability of the damage intensity index Ke, the discrete element method is used to perform microparameter matching with each lithological parameter under the same cutting conditions. Subsequently, discrete element simulations (PFC2D) are performed using conical picks with linear cutting (Fig. 7). During each cutting cycle, the cutting force acting on the pick at the point of first crack initiation is considered as the fracture CF. Furthermore, the highest CF observed within one cycle of rock failure is regarded as the PCF. Over 5 cutting cycles, the fracture CF, PCF, and their averages are recorded. Finally, the value of Ke is determined based on the simulations results.
Figure 7
Simulation of conical pick linear cutting test using PFC.
For each cutting stage of the simulation process, as the conical pick compresses the rock, the CF increases linearly while the number of cracks remains constant, corresponding to the elastic stage. By detecting the change in the number of cracks, the CF corresponding to the initial stage of continuous crack growth is identified as the fracture CF, which coincides with the computed fracture stress in this work. The PCF during the process of crack propagation and rock block detachment inside the rock is recorded as illustrated in Fig. 8. The ratio of the average PCF to the average fracture CF over 5 cycles is then calculated as the corresponding Ke for the rock sample, consistent with the calculation process of the theoretical method. Based on the PFC simulation results for 20 different cutting conditions, it can be observed that for cutting depths of 5 and 9 mm, the simulated results of the damage intensity index exhibit the same trend as the theoretical solutions, with a maximum error of less than 15%, as depicted in Fig. 9. This finding demonstrates the agreement between simulation and theoretical calculations, providing initial validation for the reliability of employing the energy correction function (C- Kf) to modify the PCF method.
Figure 8
Variation of cutting force with cutting distance.
Figure 9
Comparison between theoretical and simulated calculation of Ke.
In order to investigate the effect of C, Kf and Ke on the prediction accuracy of PCF, the sensitivity of PCF to the three parameters is calculated numerically by combining Eqs. (31,32,34,35). C, Kf and Ke are varied up and down by 10%, and the corresponding PCF values are calculated, and the sensitivity results are recorded in Table 3 and Fig. 10.
Table 3 Sensitivity analysis of PCF to C, Kf and Ke.Figure 10
Comparison between theoreticsal and simulated calculation of C, Kf and Ke.
Drawing upon the aforementioned findings, a significant influence of increased rock brittleness (C) on the PCF is observed, leading to a considerable surge in its values. Furthermore, the impact of Ke on PCF calculations appears to be comparatively minor. Additionally, Eq. (32) and Fig. 10 suggest that Ke serves as a proportional correction factor for PCF. In summary, rock brittleness (C) has the utmost significance on PCF prediction outcomes, followed by Ke, whereas Kf exerts the least influence on PCF prediction results.
The correlation between the theoretical uncorrected and corrected PCF values and the experimental PCF values is illustrated in Fig. 11. To evaluate the predictive capabilities of the method proposed in this study against existing models, this paper introduce the Evans model Eq. (36), Roxborough model Eq. (37), and Goktan semi-empirical model Eq. (38) and the corresponding PCF values for the samples based on these models are calculated. The correlations between the predicted values of each model and the experimental values are presented in Table 4.
$$PCF=\frac{16\pi {{\sigma }_{t}}^{2}\cdot {d}^{2}}{{\text{cos}}^{2}\alpha \cdot {\sigma }_{c}},$$
(36)
where σt represents the tensile strength of the rock, d denotes the cutting depth, and σc refers to the uniaxial compressive strength of the rock.
$$PCF=\frac{16\pi {\sigma }_{\text{c}}{d}^{2}{\sigma }_{\text{t}}^{2}}{(2{\sigma }_{\text{t}}+({\sigma }_{\text{c}}\text{cos}\alpha /((1+\text{tan}f)/\text{tan}\alpha )){)}^{2}},$$
(37)
where f denotes the friction angle between the pick and the rock.
Figure 11
Results of correlative analysis between theoretical model calculation results and experimental results. (a) the uncorrected method of this paper, (b) the corrected method of this paper, (c) Evans model, (d) Roxborough model, (e) Goktan Semi-empirical model.
Table 4 Comparison of prediction accuracy among different models.
$$PCF=\frac{4\pi {\sigma }_{\text{t}}{d}^{2}{\text{sin}}^{2}(\alpha +f)}{\text{cos}(\alpha +f)}.$$
(38)
Based on the data presented in Table 4, the uncorrected method proposed in this paper shows superior performance in terms of R2, RMSE, and RSS when compared to other models. Furthermore, the corrected method demonstrates a significant improvement over the uncorrected method. The specific analysis is as follows:
(1). Coefficient of determination (R2): The uncorrected method in this study has a coefficient of determination of 0.88776, which increased to 0.91104 after correction. The corrected method demonstrates superior performance compared to the uncorrected one, exhibiting a 2.6% increase in the coefficient of determination. This implies that the corrected method is capable of explaining a higher proportion of the variance in the dependent variable.
(2). Root mean square error (RMSE): The initial method had an RMSE of 3.67451, which decreased to 3.27134 after correction. The corrected method has a smaller RMSE compared to the uncorrected method, indicating higher predictive accuracy and reduced differences between the predicted and actual values.
(3). Residual sum of squares (RSS): The initial method had an RSS of 243.03676, which decreased to 192.63038 after correction. This reduction of 20.7% in the RSS implies an improved model fit and decreased prediction errors in the corrected method.
The corrected method exhibits an increase in R2 and a decrease in both RMSE and RSS, indicating enhanced accuracy and effectiveness in predicting and explaining the data. Overall, the corrected method demonstrates superior performance compared to other models.
To further validate the validity of the linear relationship between the theoretical calculation results obtained from the regression analysis and the experimental values, a variance analysis was conducted at a 95% confidence level. In the analysis, the experimental results and theoretical calculation values are treated as the dependent and independent variables, respectively. The analysis results, presented in Table 5, demonstrate that the P-value in the F-test for the linear relationship between the calculation results of the five models and the experimental results is significantly less than 0.05. This indicates that the linear correlations obtained from the regression analysis are all valid. In Table 5, MSD is the mean squared deviation.
Table 5 Analysis of variance for models.
According to the aforementioned results, the corrected method demonstrates superior predictive performance compared to other models, as it exhibits closer alignment with actual data. Moreover, the corrected method shows significant improvement over the uncorrected one.
The improvement is mainly caused by the following factors:
(1). Refinement of fracture conditions: Existing theories assume specific bedding structures or follow specific rock fragmentation patterns, and only consider one or two rock parameters to differentiate rock properties. This limited approach results in a lack of universality in predictive theories of PCF. In contrast, the Griffith theory used in this study determines the crack initiation length for rocks of different types from an energy perspective. It incorporates multiple rock parameters such as elastic modulus, KIc, and UCS to consider the influence of rock properties on elastic energy and surface energy. As a result, it exhibits enhanced adaptability to rocks with diverse bedding structures and fragmentation forms.
(2). Accurate modeling of conical picks: The conical pick consists of a pick tip, a weld, and a pick body, with the pick tip having a curved profile and the pick body having a straight profile. Existing theories simplify the conical pick as a two-dimensional pick tip model within a plane, primarily considering the cutting stresses between the tip and the rock. In contrast, the model in this study is based on the projection contour method combined with a simplified elliptical profile. It comprehensively analyzes the stress distribution in both the axial and radial directions of the pick, thereby accurately reproducing the cutting stresses on the conical surface of the pick to the greatest extent.
(3). Realistic stress distribution function: Based on the research results of D. L. Sikarskie stress-optical experiments on the pick surface, the D. L. Sikarskie stress distribution function is the closest match to the experimental results. The method in this study relies on this theory when analyzing the axial forces on the pick teeth, while the theories proposed by Evans, et al., assume a linear decrease in teeth tip forces, which is an idealized hypothesis lacking experimental evidence.
(4). Energy-based crack propagation law: The rock fragmentation process varies significantly for different rock types and bedding structures, and the crack propagation process is highly unstable. Existing theories founded on specific fragmentation laws, resulting in the instability of predictive outcomes. In contrast, the method in this study evaluates the energy relationship between rock crack initiation and propagation stages based on the Griffith theory and the destruction intensity index, grounded in physical principles. Therefore, it can accurately describe the rock crack initiation and propagation processes.
Based on the above analysis, rocks with varying lithologies and bedding structures demonstrate substantial distinctions in the fracturing process, characterized by highly unstable crack propagation. Existing theories proposed by Evans et al., with their reliance on specific fragmentation patterns, exhibit significant inaccuracies. In contrast, the energy relationship governing the onset of rock fracture and subsequent crack propagation stages can be more accurately described by utilizing the Griffith theory and damage intensity index.
Reliability verification of energy correction function
Based on the research in Sect. “Damage intensity index and PCF correction”, the elastic energy normalization factor can be mathematically expressed by utilizing the brittleness index (See Eq. 35).
This section presents the verification of the reliability of the correction process based on 20 data sets. Table 6 lists the rock sample parameters, test conditions, and experimental results used, where rock numbers 1–12 are quoted from the Belgin study and rock numbers 13–20 are obtained by rock mechanical properties test and single-tooth cutting tests, see Figs. 12, 13.
Table 6 Rock characteristic parameters in the validation set and calculation results of PCF.Figure 12
Rock mechanical properties test. (a) the TENSON concrete pressure testing machine, (b) the uniaxial compression test of gritstone, (c) the brazilian splitting test of gritstone.
Figure 12 shows the material performance test of gritstone, which is carried out on the TENSON concrete pressure test machine. Figure 12b shows the uniaxial compression test of red gritstone with loading speed of 0.2 mm/min and specimen size of Φ50 × 100 mm. Figure 12c conducts the Brazilian splitting test under axial loading of disc-shaped specimen with specimen size of Φ50 × 25 mm. The single-tooth cutting test is shown in Fig. 13. Single-tooth cutting experimental bench with hydraulic drive unit, drive the rock body to do straight-line movement, while the tool is fixed mounted on the gantry, to achieve a fixed depth of linear cutting. The tests are conducted using AID-U95-25 conical pick with a tip angle of 90 degrees and an alloy head diameter of 15.5 mm. The cutting angle is set at 40 degrees, and the cutting speed is 12.7 cm/s.
Figure 13
Single-tooth cutting test. (a) the Single-tooth cutting experimental bench, (b) the gritstone cutting test.
For each rock sample, the damage evolution model, which follows a normal distribution, is derived using UCS and critical strain (εC). The elastic energy normalization factor (Kf) is calculated using Eq. (35), while the damage intensity index (Ke) is determined by Eq. (31). The results of the calculations are summarized in Table 6.
To assess the accuracy and performance of the regression model, correlation analysis is conducted on five models. The mean absolute error (MAE) and root mean square error (RMSE) are computed to measure the disparity between the calculated values and experimental values. Results are presented in Table 7. Among the validation set, the current method demonstrates superior predictive ability compared to other models. The uncorrected method in this study achieves a coefficient of determination (R2) of 0.81591. After applying the correction function derived in Sect. “Damage intensity index and PCF correction”, the R2 increases to 0.90404, indicating a 10.8% improvement over the initial value. Besides, the mean absolute error (MAE) decreases by 34.8%, and the root mean square error (RMSE) decreases by 20.8%, indicating a better model fit to the experimental data.
Table 7 Comparative analysis of prediction accuracy in the validation set across different models.
In conclusion, based on statistical analysis and evaluation metrics, it is evident that the current method surpasses other models in terms of its predictive capability for the validation set. Moreover, the corrected method exhibits enhanced performance, characterized by smaller prediction errors and better fitting, consequently improving the predictive accuracy. This result further validates the reliability of the PCF correction method using the energy correction function with the damage intensity index described in Sect. “Damage intensity index and PCF correction”.
