The presence of joints and other types of discontinuities has a significant effect on the mechanical properties of rock, especially for tensile properties to fundamentally influence the stability of rock excavations. The main challenge associated with the experimental research on jointed rock lies in the difficulty to carry out amount of direct tensile tests for analysis of the effect of joint geometric parameters on mechanical properties. In this study, a particle flow model was established by utilizing the flat-joint contact model (FJM) to represent the rock materials. After microscopic parameter calibration, 53 sets of the numerical model were used for investigating the relationship between jointed geometric parameters and tensile mechanical properties. The results show that the crack initiation is related to trace length l and joint angle β, and the tensile-shear crack will appear as β increase. The uniaxial tension strength σt and β had first a weak negative correlation and then a positive correlation as the β increases, which was consistent with mathematical calculations. Furthermore, the relative importance (RI) analysis showed that the β plays a decisive role among the joint geometric parameters for affecting σt, and the effect factors of σt were joint angle β, length l, density n, and aperture d in that order. The present research can be utilized for multiple purposes in the field of jointed rock engineering, such as prediction of surrounding rock instability analysis and estimating the variable values in the inversion analysis in practical engineering projects.

In most of rock engineering cases, there are a large number of discontinuities (such as fissures, joints, and weak surfaces) in rocks because of geological movement and the tectonic stress field [1]. In particular for surrounding rock mass of underground coal mines that exhibit a geologically stratified structure due to the sedimentary effect of coal-forming process, a large number of joints are spread throughout the coal mining roadway [2]. The rock mechanical properties will degrade seriously due to these joints; otherwise, a lot of new cracks will form near the initial joints under loading, and the new cracks may propagate or coalesce with other cracks, which in turn leads to a further degradation in the rock mechanical properties [1, 3]. Therefore, understanding the relationship between joint geometric parameters and rock mechanic behavior is of importance in coal engineering.

Tensile stress and tensile stress gradients are often encountered in underground opening (such as roadway, tunnel, and deep basement) although the general in situ stress state for underground rocks is normally compressive [4]. Due to the much smaller tensile strength of rocks as compared to their compressive strength, tensile failure is the main and significant failure mode of rocks [5]. Consequently, tensile strength is an important material parameter for rocks. In addition, more comprehensive explorations of rock compression tests have been performed by scholars [6]. Jiefan et al. [7] conducted uniaxial compression loading tests on marble specimens containing single joints and observed the crack development stages and determine the crack nature (tensile/shear) in the process. Cao et al. [1] according to conducted servo-controlled uniaxial compression loading on rock-like specimens containing multifissures indicate that the strength of a multifissure specimen is affected by the fissure angle and number of fissures. Xu et al. [8] investigated the effect of printing layer thickness on the mechanical behaviors of rock-like specimens with sand powder printing by uniaxial compression testing, demonstrating that the increased layer thickness of the 3D print results in the uniaxial compression strength (UCS) decrease. At present, the tensile strength can be achieved in a couple ways: (1) Brazilian test, the most popular indirect tension test for rocks due to ease in alignment and low cost in experimentation [911]; (2) three-point bending (TPB), a standard for determining the flexural tensile strength of materials such as natural and artificial building stones, rocks, cement, and concretes [4, 12, 13]; and (3) semicircular bending (SCB), which is modified by TPB and is one of the famous indirect methods for determining the flexural tensile strength of rocks [14, 15]. But the tensile strength is defined as the failure stress under pure uniaxial tensile loading, and the indirect test showed substantial limitations regarding mechanical properties and failure pattern for joint rock; the direct tension method (e.g., direct pull test) is thus the most suitable method for the tensile research of joint rock [4, 15]. However, in practice, the following challenges remain in laboratories for conducting direct uniaxial tensile tests: (1) unfavorable stress concentration over the grip and (2) bending moments due to noncoaxial gripping and curvature of the specimen. Various attempts have been made in this regard [4, 16]. The numerical simulation is a good candidate for the study of tensile behavior and failure pattern of rock.

Saksala [17] developed a numerical approach, studied the detrimental effect of thermally induced cracking on the tensile strength of granitic rock, and demonstrated the validity of the proposed approach as the experimental deterioration of the tensile. Bahaaddini et al. [18] and Li et al. [19] believed that particle flow code (PFC) presents advantageous features for tracing the asperity degradation and development of cracks inside the intact material and investigated the differential effect of joint geometric parameters on the failure pattern under direct shear and uniaxial tension, respectively. It is well known that the parallel bonded model in PFC yields unrealistically low ratios of compressive to tensile strength and unrealistically low friction angles for rock [20]; then, this problem is tackled by the flat-joint model (FJM) analysis proposed by Potyondy [21] in 2012, and the FJM is widely used in the analysis of brittle rock (such as coal, sandstone, and marble) [2224].

The relationship between the main microparameters and overall mechanical properties was investigated in the present study, and the relationship is used for rapid and effective calibration of the numerical model. In addition, a large number of jointed models are created and utilized to analyze the joint weakening effect of rock tensile mechanical properties. Finally, a quantitative analysis of the relationship between joint geometric parameters and mechanical properties was carried out using relative importance (RI) analysis, which provides a basis for estimating the variable values in the inversion analysis in practical engineering projects. The present research can be utilized for multiple purposes in the field of jointed rock engineering, such as prediction of surrounding rock instability and slope stability analysis.

2.1. Calibration of Model Parameters

On the basis of the PFC2D software package, a particle flow model was established by utilizing the flat-joint contact model (FJM) to represent the intact rock materials, and a cylindrical specimen with 100 mm in height and 50 mm in diameter is built in PFC2D. Then, the simulation of discontinuous joints was achieved by deleting the disks to simulate, and the joints are arranged at uniform intervals (the joint interval, that is, the distance between the edges of two parallel joints, is defined as D, and D=5mm); the distance between the edges of two parallel joints is from the specimen body to both ends, as shown in Figure 1, where the aperture (d), trace length (l), density (n), and dip angle (β) will be used for subsequent research of joint rock mass. During the calibration process, the loading rate applied in the model simulations was, respectively, chosen as 0.05 m/s and 0.0025 m/s for the compression test and direct tension test [25, 26], the numerical model consisted of uniformly distributed particles with the radius ranging from 0.26 mm (Rmin) to 0.39 mm (Rmin), and the random seed was set to 10001. Finally, each model contains about 14,698 particles and 39,549 contacts, which vary slightly among different joint sets.

The calibration of microscopic parameters is the most important basic work during the establishment of the numerical model. The trial-and-error process is the most commonly used method in the calibration, that is, by continuously adjusting the microscopic parameters in the numerical model to match the results of the physical test. The microparameters of the model consist of two parts: the particle microparameters (Table 1) and the contact model microparameters (Table 2). Expanding on previous studies [2628], the effect of the main microparameters on overall mechanical properties was investigated in the present study, and the descriptions and effects of the main microparameters are listed in Tables 1 and 2.

According to the trial-and-error process, the microscopic parameters were continuously adjusted in the uniaxial compression and uniaxial tensile test to match the uniaxial compressive stress-strain curve (focusing on comparing uniaxial compressive strength σc and uniaxial compressive elastic modulus Ec), tensile strength σt, tensile elastic modulus Et, etc., in which Ec and Et are the slope of the line connecting the origin and the point corresponding to 50% of the tension and compress the stress-strain curve, respectively. Eventually, a satisfactory set of microparameters was obtained and is displayed in Table 3.

The main microscopic parameters shown in Table 3 were used in the uniaxial compress test of numerical simulation, and the stress-strain curve is obtained and shown in Figure 2, which showed a good agreement with coal (taken from the roof of the Liangbaosi mine). The σc obtained by numerical simulation is 10.50 MPa, which is only 0.11 MPa different from the experimental results (10.61 MPa). The Ec obtained by numerical simulation is 1.40 GPa, which is only 0.23 GPa different from the experimental results (1.17 GPa).

2.2. Validation Study and Tensile Mechanical Property Analysis of Intact Specimens

The present paper is mainly to study the tensile mechanical properties of coal specimens through numerical simulation. Therefore, in addition to comparing the stress-strain curves under uniaxial compression, the calibration of parameters should also focus on comparing the ratio of compressing strength to tensile strength (σc/σt) between the numerical model and coal, and study of the mechanical behavior under direct tension should be taken seriously. The stress-strain curve and the statistics of acoustic emission (AE) counts during the tensile failure process are shown in Figure 3. It can be seen that the stress increases linearly with the loading, but the AE counts do not appear during the stress growth. Until the stress increases to about 95% of the peak intensity, the specimen begins to have AE counts and continues to load, the stress increases to the peak intensity, the specimen appears brittle failure, energy is released, and a large number of AE counts are generated; the stress-strain curve also falls. The uniaxial tensile strength is 0.49 MPa, and σc/σt of the numerical simulation is calculated as 21.43, which is close to the coal [2931]. The tensile elastic modulus is 1.37 GPa, and the ratio of compressing elastic modulus to tensile elastic modulus (Ec/Et) is calculated as 1.02, which complied with the regular for high σc/σt and low Ec/Et.

The numerical test adopts displacement loading, and the displacement distribution with respect to the failure process of the specimen under the test process represents the transfer process of the internal load of the specimen to some extent as shown in Figure 4. And the tensile f compression process can be divided into stage I: linear elastic stage, stage II: peak failure stage, and stage III: postpeak stage of load. As shown in Figure 4, it can be seen that before reaching the peak stress, the specimen is basically in stage I. At this time, the internal stress of the specimen increases linearly, and the load is uniformly transmitted from both ends to the middle inside the specimen. And due to the characteristics of being symmetric and deformable, the displacement of the particles inside the specimen tends to gradually decrease from the two ends to the middle. As the specimen reached stage II, cracks appear inside the specimen, and the particle displacement no longer changes uniformly along the axis of the specimen from the vicinity of the crack, but the displacement of particles farther from the crack still changes uniformly. Subsequently, the specimen reached stage III. The cracks rapidly expanded and persisted throughout the specimen resulting in complete failure of the specimen. The entire failure stage of the specimen is very short and shows obvious brittle failure. Different from the uniaxial compression failure process of the specimen, the specimen has almost no residual strength after the peak stress, and the crack propagation process after the peak is also more rapid, which is similar to the real rock [32, 33].

Previous studies have indicated that the geometric parameters, especially joint angle (β), have a significant effect on the mechanical properties of the specimen [3436]. Therefore, the research in this part focuses on the β, and the research of aperture (d) trace length (l), and density (n) is combined with the β for analysis, in which a total of 49 sets of specimens with different joint geometric parameters were built, and the joint geometric parameter settings are shown in Table 4.

3.1. The Effect of β on the Tensile Mechanical Properties of Rock

For the study of the effect of the joint angle on the tensile properties of rock, seven models with β of 0°, 15°, 30°, 45°, 60°, 75°, and 90° were established in PFC2D, and then, direct tensile tests were conducted on these models to obtain a set of curves of peak strength with β and the failure pattern of the specimens when β was taken at different values as shown in Figure 5. It is found that the tensile strength σt of the specimen is obviously positively correlated with the joint angle β, and there is a steeper increase in σt when the β is large. For example, the σt only increases by 16.7% from 0.18 MPa to 0.21 MPa as β increases from 0° to 45° (three gradients), but the σt increases by 39% from 0.31 MPa to 0.43 MPa as β increases from 75° to 90° (only one gradient). At the same time, the failure pattern of the specimen also changes with the various β. When β60°, the failure pattern is mainly tensile failure, while when β75°, the specimen starts to show tensile-shear failure; especially whenβ=90°, a large number of shear cracks appear after the specimen is completely failure.

Figure 6 shows the displacement distribution with respect to various β, peak failure and final failure are the failure pattern at the stress reaching σt, and the residual strength decreases to 30% σt, respectively. It can be seen that the displacement of the particles at both ends of the specimen is still uniform until the specimen is failure, but in the middle of the specimen, the displacement of the particles is affected by the prefabricated joints, which shows an obvious circular change, indicating that the displacement of the materials near the joints under tensile stress is more different. And with the increase in β, the influence of the joints on the displacement of the particles inside the specimen becomes weaker. In particular when β increases to 90°, the circular change of displacement disappears, and the displacement of the particles gradually transitions from the two ends of the specimen to the middle and becomes smaller.

3.2. The Effect of d on the Tensile Mechanical Properties of Rock

Figure 7 shows the relationship between σt and β under different d. It can be seen that β has a significant effect on σt. Taking d=1mm as an example, as β increases from 0° to 30°, the σt decreases in a negatively correlated manner to 0.17 MPa, with a drop of 22.7%. However, the correlation between σt and β becomes positive as β increases from 30° to 90°; the σt of β=90° specimen reaches 0.47 MPa, which is dramatically increased with 176.5% from β=30°. Similar affecting patterns can also be found with different d from d=1.5mm to d=2.5mm.

Figure 8 shows the displacement distribution with respect to various d, and it can be seen that the joint aperture d has little effect on the failure pattern of the specimen; in particular when the angle is small, such as when β=0°, the specimen deformation is nearly unaffected by β no matter in the stage of peak failure or final failure. But when β increases, d will also have a certain effect on specimen deformation. Taking β=90° as an example, the increase in d will create deformation on both sides in the stage of peak failure, but the final failure will become more uniform as d increases.

3.3. The Effect of l on the Tensile Mechanical Properties of Rock

Figure 9 shows the relationship between σt and β under different l. It can be seen that both β and l have a significant effect on σt. The σc showed a decreasing tendency first and then an increase with the increase in β under lower l. As d increases, σt turns to be significantly positively correlated with β, and the growth rate of σt increases dramatically with the growth of β. Taking l=20mm as an example, the σt only increases by 18.2% from 0.11 MPa to 0.13 MPa as β increases from 0° to 45°, but the σt of β=90° specimen reaches 0.37 MPa, which is dramatically increased with 184.6% from β=45°.

For trace length l, the σt exhibited a strong negative correlation especially when β=0°. As β increases from l=5mm to l=20mm, the σt decreases in a negatively correlated manner from 0.25 MPa to 0.11 MPa, with a drop of 56%. The weakening of l to σt will decrease with the increase in β until β=90°, and the trace length l has almost no influence on σt.

Figure 10 shows the displacement distribution with respect to various l, and it can be seen that the increase in l will affect the uniformity of particle displacement near joints. In addition, it can be found from Figure 10 that the trace length l will affect the initiate position of crack as the β increases. Take β=90° as the example; cracks initiate from the tip of the joint at l=5mm, l=10mm, and l=20mm. But when l=15mm, the crack will initiate from the middle of the joint, and it did not vary with d.

3.4. The Effect of n on the Tensile Mechanical Properties of Rock

Figure 11 shows the relationship between σt and β under different n. It can be seen that both β and n have a significant effect on σt. When β=0°, the σt is positively correlated with joint density n, and the σt increases by 27.8% from 0.18 MPa to 0.23 MPa as n increases from 1 to 4. But with the increase in the joint angle β, the relationship between σt and n gradually changed from positive correlation to negative correlation. Take the β=45° as an example; the σt of n=1 specimen is 0.21 MPa, with increase from 1 to 4; the σt decreases in a negatively correlated manner to 0.13 MPa, with a drop of 38.1%. However, when β reaches 90°, the σt of n=3 specimen is increased in a positive correlation manner from 0.26 MPa to 0.38 MPa, which is dramatically increased with 46.2% and does not satisfy the regularity.

In terms of the failure pattern, significant effect of both β and n can be clearly found from Figure 12. For the effect of β with n=1 as an example, the presence of a joint can significantly affect the deformation inside the specimen, especially near the joint. But the effect becomes weaker when the β increases; the similar affecting patterns can also be found with different d from n=2 to n=4. For the effect of n with β=90° as an example, variation of n affects the deformation on both sides of the joint. In particular when n=2, the deformation on both sides of the joint is quite different, which is different from other n. And it can be seen from Figure 12 that when n=2, the σt is smaller and the inflection point was observed on the β-σc curve. Notably, the crack initiation is also affected by β and n. When β=0°, the cracks initiate at the tip of a joint and expand along the joint until the specimen fails. As the β increases, cracks usually initiate at the tip of multiple joints, and the cracks will expand and connect the joint.

In this part, the mathematical analysis will be used to study the associations between σt and joint geometric parameters, and the result will be compared with the numerical test to verify the reliability of numerical simulating jointed rock. Furthermore, the relative importance (RI) of β, d, l, and n was quantified.

4.1. Verification of the Results

According to the above study, it is found that the β has the most obvious effect on the mechanical properties of the rock mass. In particular, the σt and the β have an obvious relationship that changes from being negatively correlated to positively correlated, which will be analyzed in this part.

The geometric relationship between the joint and specimens is analyzed as shown in Figure 13, where l1 is the horizontal projection length of the l, l2 is the horizontal projection length of the d, and lmin is the minimum cross-sectional length at which the load can be transferred on both sides of the joint (lmin is distributed symmetrically on both sides of the joint).

According to Figure 13, the l1 and l2 can be calculated as follows:
(1)l1=l·cosβ,l2=d·cos90β=d·sinβ.
Then, the calculations of the lmin in Figure 13 are as follows:
(2)lmin=Dl·cosβd·sinβ2.
It can be seen from equation (2) that lmin is inversely proportional to l and d and proportional to β when the effect of d is ignored. The uniaxial strength of the specimen is calculated as follows [37]:
(3)σ=PA,
where σ is the compressive strength of the jointed specimen, P is the peak load, and A is the bottom area of the specimen.
For this analysis, the subsequent calculation will be based on the following underlying assumptions: (1) the internal material of the specimen satisfies the homogenization assumption and isotropy assumption [38] and (2) the load is uniformly transmitted inside the specimen. On the basis of the above assumptions, the value of P is positively correlated with lmin, and the strength relationship between the jointed specimen and the intact specimen can be obtained as follows:
(4)σ=2lminDσ,
where σ is the compressive strength of the intact specimen.
Obviously, the σ is also positively correlated with lmin; then, equation (2) is incorporated into equation (4), and the strength relationship, with the change of β, between the jointed specimen and the intact specimen can be obtained as follows:
(5)σ=1l·cosβ+d·sinβDσ.

From equation (5), it is clear that the σ is positively correlated with β but negatively correlated with l when d is ignored. Moreover, the curve of σ varies with the β obtained as shown in Figure 14 (taking l=15mm and d=1mm as an example). According to the curve, it can be seen that σ showed a decreasing tendency first and then an increase with the increase in β. And the lowest point of the curve is at β=4°, when the JISR is 0.617.

4.2. Analysis of Relative Importance

From the above analysis, it is seen that β, d, l, and n were shown to have a regular effect on σc, but which joint parameter affects σt most still needs further investigation to provide an insight for the study of surrounding rock control, while the relative importance (RI) analysis is a good candidate.

The linear regression can be expressed as follows:
(6)y=a+j=1Jbjxj+e,
where y is the dependent variable (σc) and x is the explanatory variable (β, d, l, and n).
According to Fields [39], the contribution of each explanatory variable towards the R2 can be isolated as R2y:
(7)R2y=j=1JbjCovxj,yVary=1Cove,yVary,
where Vary is the total sum of squares, Covxj,y is the regression squares, and Cove,y is the error squares.
There is one way to represent the contribution: when the subset model contains only variable k, the contribution of k can be displayed as follows:
(8)Rk2y=Rk2a+xkxk+e.
The preceding analysis indicates that the effect of one joint parameter on σc significant correlations exists with other joint parameters. In this regard, Shorrocks [40] argues that the contribution of the explanatory variable should be equal to the marginal effect (M) of the explanatory variable to R2, and the M is calculated as follows:
(9)Mk=R2y=a+bkxk+jSbjxj+eR2y=a+jSbjxj+e,
where S is the other explanatory variable that does not contain the variable k. Since the regression coefficients usually change when one explanatory variable is removed, the regression coefficients that do not include variable k are marked with .

Due to three explanatory variables in the model (β, d, l, and n), the original model corresponds to seven subset models. Regression has to be done on all subset models to calculate the relative contribution of each variable using the average value as the final calculated degree of contribution of the variable.

To facilitate analysis, the method of Krasikova et al. [41] and Tonidandel and LeBreton [42] is used to standardize the values of RI reported below. Specifically, the RI of all variables is aggregated into RI total. Then, the ratio of the RI of each variable (dominance stat, Ck) to RI total is calculated to get the standardized degree of contribution (standardized dominance stat, RIk), and the equation is shown below:
(10)RIk=Ckj=1JCj.

This standardization has the advantage of making the sum of the standardized degrees of contribution of all explanatory variables equal to one, making the relative importance of each variable becomes easily comparable with others [43]. The Stata command “domin” is used for RIk analysis in this paper, and the results are displayed in Table 5.

According to the above analysis, it can be seen that the RI of β is 65.54 times higher than that of d, 9.24 times higher than that of l, and 11.97 times higher than that of n. It was indicated that the β plays a decisive role on the mechanical behavior of the specimen among the joint geometric parameters, and it should be firstly considered avoiding the joint of the weakest angle (the angle is 30° in this study) in rock mass to improve the bearing capacity.

In present paper, a particle flow model was established by utilizing the flat-joint contact model (FJM) to represent the intact rock materials. And the relationship between the main microparameters and overall mechanical properties was investigated, which is used for rapid and effective calibration of the numerical model. 53 sets of the numerical model were used for investigating the relationship between jointed geometric parameters and uniaxial tensile strength σt, which was confirmed by mathematical calculations. In addition, the relative importance (RI) of β, d, l, and n was analyzed and quantified.

A control variable experiment was performed by numerical simulation to investigate the effect of jointed geometric parameters (β, d, l, and n) on the σt and failure pattern. There is an obvious relationship between the σt and β showing first a weak negative correlation and then a positive correlation, consistent with mathematical calculations. For the failure pattern, it can be seen that the tensile-shear crack will appear as β increases, and the crack initiation is related to l and β. Typically, cracks are initiated at the tip of the joint, but when β increases to 90° and l=15mm, the crack will initiate from the middle of the joint.

According to calculations, it clearly indicated that lmin is a primary factor for σc of the jointed specimen. And it is found from the mathematical calculations that lmin is closely related to the β. In addition, the calculated trend of σc versus β is consistent with numerical simulation results, which yielded very strong evidence in favor of the practicability of the present research. Furthermore, the relative importance analysis (RI) is conducted to quantify the effect of jointed geometric parameters on σt. The results show that the β plays a decisive role among the joint geometric parameters, and the effect factors of σt were joint angle β, trace length l, density n, and aperture d in that order. The present study provides fundamental results for the further study of jointed rocks and serves as an insight for further study of surrounding rock control, and the RI analysis results can provide a basis for estimating the variable values in the inversion analysis in practical engineering projects.

Data used to support the results of this study can be found in this manuscript text.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This study was financially supported by the National Natural Science Foundation of China (52074166) and National Natural Science Foundation of Shandong Province (ZR2021YQ38) and the Open Grant of State Key Laboratory of Mining Response and Disaster Prevention and Control in Deep Coal Mines (SKLMRDPC20KF02).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).