When inducing cracks, soundless cracking agents (SCAs) do not generate vibration, harmful gas, dust, nor flying rock fragment, making them suitable for hard rock roof breaking, rock burst prevention, oil or gas reservoir stimulation, and building demolition. In this study, SCA-induced crack initiation and propagation in different stress conditions were modelled using a modified cohesive element method. A new traction-separation law for describing rock compressive shear strength was proposed. The crack length and direction in bidirectional isobaric and unequal stress fields were analyzed in detail. The crack initiation pressure and the incremental ratio of crack length to SCA expansion pressure were proposed as two indicators to evaluate the difficulty in rock breaking in deep underground. Results indicate that (1) the modified cohesive element method used in this study is feasible to model crack propagation in deep rocks; (2) the maximum expansion pressure of SCAs depends on rock elastic modulus and geostress field and should be measured under a condition similar to deep underground prior to SCA borehole spacing design; when using the SCAs with a maximum expansion pressure of 100 MPa in 600 m underground, the suggested borehole spacing is less than 220 mm; and (3) σ3 dominates the crack initiation pressure while the principal stress ratio σ3/σ1 and notch direction control the direction of crack propagation.

1. Introduction

The morphology, direction, and propagation length of artificially induced cracks are the important parameters for rock or concrete fracturing projects [13]. The common methods of rock and concrete breaking include rock blasting [4, 5], diamond wire cutting [6], controlled blasting [7], hydraulic fracturing [8], and SCAs [9]. SCAs are hopeful to be used in hard rock roof breaking in coal mines (Figure 1(a)), rock burst prevention in metal mines, oil or gas reservoir stimulation (Figure 1(b)) [10], and building demolition [11, 12]. In addition, SCAs are good prospects for rock breaking in subsea tunnel projects as they do not generate vibration, harmful gas, dust, and flying rock fragments when inducing cracks [11, 13].

SCAs primarily consist of calcium oxide (CaO) and some other additives such as magnesium (MgO), aluminum oxide (Al2O3), and silicon dioxide (SiO2). CaO is the main constituent of SCA. When adding water into SCAs, the hydration reaction between water and CaO generates calcium hydroxide (Ca(OH)2) and expands in volume. The main role of the additives is to control the hydration process. In general, SCAs are filled into boreholes to expand and push the borehole wall [14, 15]. The expansion induces radial compressive stress and tangential tensile stress in the surrounding rock. It is commonly believed that tangential tensile stress is the main factor of crack initiation and propagation. Thus, different multihole arrangement patterns create different rock crack morphologies to meet various engineering requirements.

SCAs originated in the 1970s, and since then, many researchers have studied its rock breaking mechanisms. SCA-induced stress distributions under different borehole sizes, rock strengths, and temperatures were studied by Laefer et al. [16] and Mateus and Araújo [17] independently. Dar et al. [18] proposed the optimal spacing of the SCA borehole for breaking marble. De Silva et al. [19] summarized SCA utilities in actual rock breaking projects and put forward its feasible usage in oil formation stimulation. Hanif and Al-Maghrabi [20] developed a SCA application for the granite and marble processing industry. Arshadnejad et al. [21] investigated the crack propagation pathway induced by multihole filled with SCAs. Zhongzhe et al. [22] and Jana [23] proposed a method to determine optimal SCA borehole spacing. Tang et al. [11] analyzed SCA-induced crack initiation and propagation in heterogeneous materials using the RFPA software and proposed that SCA expansion pressure depends on stiffness and confining pressure of the surrounding rock. Tang et al. [24] investigated its feasibility to improve coal permeability. Temperature effect on SCA expansion pressure was also studied numerous times [15, 2527]. For deep rock projects, the initiation, direction, and length of SCA-induced cracks are influenced by stress field. However, very few studies focused on the effects of geostress field on SCA-induced crack initiation and propagation.

Crack propagation modelling issues in rock have been addressed by many studies. Many numerical methods have been performed to model fracture propagation in rock, such as the solid element deletion method [28], the extended finite element method (XFEM) [29], the virtual crack closure technique (VCCT) [30], the cohesive element method [31], the discrete element method (DEM) [32], the smoothed particle hydrodynamics method (SPH) [33], and the discontinuous deformation analysis (DDA) method [34]. The solid element deletion method, as its name suggests, simulates crack propagation by deleting solid elements, so it cannot accurately model contact and slip behaviors between crack surfaces. XFEM is usually used to model brittle crack propagation, but it is difficult to simulate coalescence and branching between a large number of cracks in rock mass. When using VCCT, the location of initial crack is usually required, and besides, it is only able to simulate crack propagation along a prespecified cracking path. Therefore, VCCT is not suitable for modelling rock fracturing. The cohesive element method avoids stress singularity in crack tip by adopting a nonlinear damage zone around the crack tip, resulting in good numerical convergence. Compared to the other methods mentioned above, the cohesive element method can also simulate complex contacts and interactions between cracks after rock fracturing.

This study aims to investigate the characteristics of SCA-induced crack propagation in rock under different stress conditions. As shown in Figure 2, the notch direction is a significant parameter for controlling the location of crack initiation and direction of crack propagation, so the notch direction was also considered in this study. This paper developed a modified cohesive element method that uses rock compressive shear strength to model SCA-induced crack initiation and propagation. This study is of significance to deepen the understandings of the characteristics and mechanism of SCA-induced crack propagation in deep underground.

2. Methodology for Modelling Crack Propagation in Rock

2.1. Cohesive Element Method

Based on the cohesive element method, crack initiation and propagation can be simulated using the finite element method. As shown in Figure 3(a), the 2D zero-thickness cohesive element has 4 nodes and 2 integration points. Its geometric thickness is zero, but its constitutive thickness, used to calculate its deformation stiffness, is not equal to zero. The traction forces on the top and bottom faces induce the relative separation displacement between them, and the relative separation displacement components in local 1 direction and 2 direction (Figure 3(a)) represent crack opening and slipping, respectively.

Figure 3(b) illustrates the process of generating zero-thickness cohesive elements in a solid element mesh. The algorithm used was written in the Python programming language and implemented as an ABAQUS plugin. By executing the plugin, cohesive elements are generated at every solid element interface; therefore, the final model consists of solid elements (the yellow elements) and cohesive elements (the blue elements). A solid element and its neighboring cohesive element share two common nodes, and the two adjacent cohesive elements share one node. Crack initiation and propagation are represented by the damage and failure of cohesive elements. Rock elastic-plastic deformation is dominated by solid elements.

The mesh should be fine and irregular to avoid mesh structure effect on crack propagation paths and to improve crack propagation randomness. In this study, the irregular mesh was generated using the advancing front technique. The final mesh model is shown in Section 4.3.

2.2. Initial Stiffness of Cohesive Element

As shown in Figure 4(a), when two solid elements are in series, their total stiffness can be calculated bywhere is the initial stiffness of the solid element.

When embedding one cohesive element between two solid elements, the new total stiffness can be calculated bywhere is the initial stiffness of the cohesive element.

According to Equations (1) and (2), the total stiffness of the finite element model is affected by embedding cohesive elements. To avoid this effect, must be much greater than in Equation (2).

The relationship between element stiffness K and elastic modulus E iswhere h is the constitutive element thickness. For a solid element, h is the same as its geometric thickness. For a zero-thickness cohesive element, h is not equal to its geometric thickness (that is, zero). Equation (3) indicates that the constitutive thickness of a cohesive element must be much smaller than the geometric thickness of a solid element to not influence the global mesh stiffness. In this study, the elastic modulus of cohesive element and solid element are the same, and the initial constitutive thickness of cohesive element hoc is one thousandths of the constitutive thickness of solid element hs.

3. Damage Initiation and Evolution of Cohesive Element

The damage initiation and damage evolution of cohesive elements are described by the traction-separation law [31]. Crack opening and slipping are caused by tractions on the top and bottom faces of cohesive elements including tensile and shear stresses. In this study, tensile stress and tensile displacement are represented by positive values, whereas compressive stress and compressive displacement are flagged as negative values. As shown in Figure 5, once the traction exceeds the strength, the cohesive element begins to damage, and its stiffness decreases as separation displacement increases. When the cohesive element is completely damaged, its stiffness reduces to zero, indicating forward crack propagation. The traction-separation law can be divided into three stages: elastic stage, damage initiation stage, and damage evolution stage.

3.1. Linear Elastic Traction-Separation Equation

It is assumed that traction-separation behavior in the elastic stage obeys the linear elastic law which is the relationship between nominal stresses and nominal strains. Nominal stresses are force components divided by original area at each integration point of a cohesive element, whereas nominal strains are separations divided by initial constitutive thickness, given bywhere and are the nominal tensile and shear strain components, respectively, and and are the corresponding separation displacement components, respectively.

The nominal traction stress vector t consists of two components in 2D problems: tn and ts, which represent normal traction (along the local 2 direction in Figure 3(a)) and shear traction (along the local 1 direction in Figure 3(a)), respectively. The elastic relationship between nominal stresses and nominal strains is represented aswhere and are the tensile and shear elastic moduli, respectively.

3.2. Damage Initiation

In deep underground surrounding a SCA borehole, stress field is affected by the SCA’s expansion pressure and crustal stress field. Rock shear strength usually increases with confining pressure; therefore, compressive shear strength characteristics are crucial for crack propagation modelling.

3.2.1. Compressive Shear Damage Criterion (tn < 0)

When tn < 0, the damage initiation of cohesive elements is assumed to obey the Mohr–Coulomb criterion. The damage criterion function is given bywhere S is the pure shear strength (i.e., the cohesion) of rock and φ is the friction angle. In the tn-ts plane, Equation (6) is plotted as two radials, as illustrated in Figure 6.

3.2.2. Tensile-Shear Damage Criterion (tn ≥ 0)

The cohesive element begins to damage when the nominal stress ratios in the quadratic interaction function below reaches 1. This criterion is represented aswhere T is the tensile strength of cohesive elements. For rock materials, T is usually smaller than S. Therefore, Equation (7) represents a semiellipse in the tn-ts plane, as shown in Figure 6. Equation (7) describes the mixed mode of tensile and shear stresses that causes damage initiation.

3.3. Damage Evolution

It is assumed that a linear softening evolution appears after damage initiation. A scalar damage variable, D, represents the overall damage in one cohesive element and captures the combined effects of all the active damage mechanisms. D is zero initially and then gradually evolves from 0 to 1 upon further loading. The stress components of the traction-separation law are affected by the damage according to the following equation:where and are the stress components predicted by the elastic law at the current strains, as illustrated in Figure 7. Note that the tensile stiffness is degraded when but unchanged when .

The mixed mode happens when cohesive elements undergo normal and shear tractions simultaneously. To describe damage evolution of this mode, an effective displacement δm [35] is defined aswhere the operator is defined as

Therefore, is equal to when (i.e., ).

According to Equations (4) and (5), tractions at the damage initiation point are calculated bywhere and are the normal and shear stress components at the damage initiation point, respectively, and and are normal and shear displacement components at the damage initiation point, respectively (Figure 7).

The mode mixity ratio β is defined as

According to Equation (9), the effective displacement at the damage initiation point, , is given by

According to Equations (12) and (13), and can be represented by β and aswhere can be obtained by substituting Equations (11) and (14) into (6) and (7), and solving for gives

Equation (15) indicates that changes with β, which is indirectly determined from the damage initiation criteria.

For cohesive elements under tensile-shear loading, the most used damage evolution law is the power law criterion [35], given bywhere GI and GII are tension and shear energy components, respectively, and GIC and GIIC are I-mode and II-mode fracture energies, respectively. When GI and GII satisfy Equation (16), the damage variable D reaches 1. GI and GII are calculated aswhere and are the tensile and shear displacement components at the complete damage point.

GIC and GIIC are the area of triangle O-- and triangle O--, respectively (Figure 7), and are calculated by

According to Equations (9) and (12), and can be represented by β and :where is the effective displacement at the complete damage point. When α = 1, is obtained by substituting Equations (18) and (19) into (16), given by

Under compressive shear loading, according to Equation (12), β tends to infinity. Thus, can be solved by letting β → ∞ in Equation (20), which is given by

For the linear softening evolution, D is expressed aswhere is the maximum effective displacement value attained during loading. As shown in Figure 8, the loading history refers to the loading path (Path 1) or the reloading path (Path 2).

D is calculated using Equation (22), and then the stresses are updated by substituting D into Equation (8). Figure 9 illustrates the traction-separation law when considering compressive shear strength characteristics. The damage initiation and evolution laws were written in the Fortran programming language as a UMAT subroutine.

4. Numerical Simulation Procedure

4.1. Mechanical Properties of Rock

The mechanical properties used in this study were obtained from laboratory tests. The rock samples are medium-grained sandstone obtained from the roof of coal seam No. 22 in the Burtai coal mine in the Ordos Basin, China. The elastic modulus and Poisson’s ratio were measured by uniaxial compression tests, the tensile strength was obtained from Brazilian split tests, the cohesion and friction angle were tested by triaxial compression tests, and the mode-I and mode-II fracture energies were determined by three-point bending tests and four-point shear tests, respectively. Table 1 lists the mechanical properties. Note that the elastic modulus Eo, cohesion c, friction angle φo, and tensile strength st of the solid elements are equal to the elastic modulus En, pure shear strength S, friction angle φ, and tensile strength T of the cohesive elements.

4.2. SCA Expansion Property

The variation of SCA expansion pressure over time can be measured in a confined state. As shown in Figure 10, SCAs are put into a thick-walled steel cylinder tube, and then the SCA expansion property will be reflected on the outer surface deformation of the cylinder tube. According to the literature [9] and [24], for any point (ρ, θ, and z) in the cylinder tube, its stress state is calculated bywhere σρ, σθ, and σz are the stress components along the radical, tangential, and axial directions, respectively; R and r are the inner and outer radii of the steel cylinder tube; and Pin and Pex are the pressures of the inner surface and outside surface.

According to Equation (23), the stress components on the outside surface of the steel cylinder tube are calculated by

According to the generalized Hooke law, the tangential strain component iswhere Ep is the elastic modulus of the surrounding materials.

By substituting Equation (24) into Equation (25), Pin can be obtained as

Equation (26) shows that Pin is influenced by both Ep and Pex. Borehole expansion pressure changes with Pex, so the SCA expansion property measured under Pex = 0 is not suitable for numerical simulation when ground stress field is taken into account.

According to Tang et al. [11] and De Silva et al. [9], the SCA expansion process is assumed to be a static process because SCAs usually take 20 h or longer to reach its maximum expansion pressure. Therefore, the loading rate of expansion pressure has very few effects on the mechanical response of its surrounding materials. In this study, the expansion pressure was applied to the borehole wall using a static loading process from 0 MPa to 150 MPa.

4.3. Mesh Model

In fact, SCA boreholes must have an opening that allows it to expand freely along the axial direction. The free expansion behavior along the axial direction is not what mainly causes rock fracturing, so the interaction between SCAs and surrounding rock along the axial direction was neglected in this study. Additionally, in order to reduce computation complexity in ABAQUS, a 2D model was used, which is equivalent to a 3D model by setting an out-of-plane dimension. As shown in Figure 11, the mesh is a 500 mm × 500 mm square. The SCA borehole diameter is 30 mm, and the distance between the two notch tips is 56 mm. The mesh, initially made of 193413 triangular elements, was generated using the advancing front technique, and then 273441 zero-thickness cohesive elements were embedded into the mesh.

The bottom of the model was fixed in the vertical direction. The maximum principal stress σ1 and the minimum principal stress σ3 were applied in both the vertical and horizontal directions, respectively. The numerical simulation consists of two steps: first, σ1 and σ3 were applied to the corresponding model boundaries, and the geostatic stress was equilibrated. Then, SCA expansion pressure was gradually applied to the borehole wall from 0 MPa to 150 MPa.

4.4. Test Schemes

A series of numerical tests were carried out to study the characteristics of SCA-induced crack initiation and propagation under different stress conditions and notch directions. Stress conditions were categorized into σ1 = σ3 = 0 MPa (free boundary conditions where the rock sample is not fixed and not loaded with any external forces), σ1 = σ3 ≠ 0 MPa (bidirectional isobaric stress field), and σ1σ3 = 5 MPa (bidirectional unequal stress field), as shown in Table 2, where θ is the angle between notch direction and σ1 direction (Figure 2).

5. Numerical Simulation Results and Discussion

5.1. Crack Characteristics under the Free Boundary Conditions

To predict crack propagation using the cohesive element method, the characteristics of SCA-induced crack propagation under the free boundary conditions were computed first. The numerical results were compared with the laboratory results of Tang et al. [11], as shown in Figure 12. Figures 12(a) and 12(b) show the test results from two different samples by Tang et al. [11]. The cracks inside the samples are similar to those on the surfaces of the intact rock samples, which is confirmed by Tang et al.’s results [24]. This is because in each of their tests [11, 24], at least one end of the SCA borehole was unrestricted, so the SCA expansion pressure along the axial direction was released and unable to affect the intact rock samples. However, if a rock sample contains obvious natural fractures, the SCA-induced cracks inside it may be different from those on its surface. The effect of natural fractures on SCA-induced cracks is not considered in this study and will be addressed in the future. Figure 12 shows that the crack angles are significantly different for the two experiments but similar in the numerical results. The crack angle difference in Figures 12(a) and 12(b) is due to different mineral grain and natural microfracture distribution. The simulations in this study did not focus on the mesoscopic heterogeneity of rock materials. The difference between the results of the simulation and that of the experiment is due to the rock parameters used, but both results show three main cracks in each sample. In ABAQUS, tensile stress is positive. In Figure 12(c), tensile stress concentrates around the crack tips, and the cracks propagate gradually as expansion pressure increases. Generally, the crack propagation characteristics simulated by the cohesive element method is similar to those in the experimental results in the literature [11], which suggests that the cohesive element method based on the modified traction-separation law can model SCA-induced crack initiation and propagation.

5.2. Crack Characteristics in Bidirectional Isobaric Stress Fields

The numerical results of crack propagation under σ1 = σ3 = 10 MPa conditions are illustrated in Figure 13. There are three main cracks under the free boundary conditions, but only two main cracks in the bidirectional isobaric stress field where the SCA-induced cracks are initiated at the notch tips and the crack directions are roughly parallel to the notch direction. Crack length gradually increases with the SCA expansion pressure P. According to Figure 13, the extent of tensile stress concentrations around the crack tips also increase as the cracks propagate. At the beginning of crack propagation, the tensile stress concentrations only appear around the crack tips. However, as P increases, the crack openings and borehole diameter expand, causing bulge and tensile stress concentration in the middle of the boundaries.

To apply SCAs in deep underground, the spacing of SCA boreholes, which depends on the crack length l, is crucial, so further analysis of the relationship between P and l under different stress conditions should be undertaken. P-l linear fit in different bidirectional isobaric stress fields is plotted in Figure 14, and the corresponding equations are listed in Table 3. According to Figure 14, it is obvious that in each stress field, SCA-induced crack length increases with its expansion pressure. The SCA crack initiation pressure P0 is a proper parameter for evaluating SCAs’ rock breaking ability (the higher the P0, the more difficult it is to break rock). Under the free boundary conditions, P0 is only approximately 4–6 MPa [11]. Under σ1 = σ3 = 5 MPa, P0 is 16.6236 MPa. However, when σ1 = σ3 = 20 MPa, P0 increases to 37.8456 MPa, which is approximately 6.3–9.5 times that under the free boundary conditions.

The slope of the P-l linear fit (unit: mm/MPa) al is the incremental ratio of crack length to SCA expansion pressure. As shown in Figure 14(b), al decreases as confining stress increases. The lower the al, the more difficult it is to break rock. The reduction rate of al gradually slows down as confining stress increases. l is approximately 500 mm under P = 60 MPa and σ1 = σ3 = 5 MPa, but 50 mm under the same P but σ1 = σ3 = 20 MPa. Overall, the results suggest that the bidirectional isobaric stress field does not significantly affect the crack direction, but can significantly affect crack initiation pressure P0 and crack length l.

According to Equation (26), the maximum SCA expansion pressure depends on the rock elastic modulus and geostress field, so it should be measured in a condition similar to deep underground before designing borehole spacing. According to Figure 14, at 600 m underground under 100 MPa of maximum expansion pressure (vertical stress of approximately 15 MPa), SCA borehole spacing should be smaller than 220 mm.

5.3. Crack Characteristics in Bidirectional Unequal Stress Fields

The crack characteristics under the conditions of σ1σ3 = 5 MPa are shown in Figure 15. When the notch direction θ = 0° (Figure 15(a)), the cracks are roughly parallel to the σ1 direction in each stress field. However, when θ = 45°, although the stress difference is just 5 MPa in each 2field, the crack direction is significantly different (Figure 15(b)). When σ1 = 5 MPa, the angle between the crack direction and σ3 direction θc is 82.8°, whereas it is 65.5°, 58.5°, and 53.9° when σ1 = 10 MPa, 15 MPa, and 20 MPa respectively. Figure 16 illustrates the variations of θc with respect to the stress field when θ = 45°. Figure 16(a) shows that the crack direction rotates gradually from the σ1 direction to the notch direction as both σ1 and σ3 increase. In Figure 16(b), it is worth noting that θc decreases linearly as the stress ratio σ3/σ1 increases. When σ3/σ1 approaches 1, θc = 45°, indicating that the crack direction tends toward the notch direction. Conversely, when σ3/σ1 approaches 0, the crack direction tends to the σ1 direction. Their linear relationship is

When θ = 90°, four cracks are initiated, two of which (named cracks N) appear at the notch tips, and the other two (named cracks M) appear at the center of the SCA borehole wall, as shown in Figure 15(c). Whether cracks can continuously propagate depends on the stress field. Cracks N only propagate continuously along the σ1 direction when σ3 = 0; however, as σ1 and σ3 increase, only cracks M can grow parallel to the notch direction. According to the above results under the conditions of θ = 45° and θ = 90°, it can be concluded that the crack propagates roughly along the σ1 direction when σ3 is close to zero; however, as σ3/σ1 increases, both the notch direction and stress field affect the direction of crack propagation.

Figure 17 shows the relationship of P-l in different bidirectional unequal stress fields and under different notch conditions. Under the condition of θ = 45°, generating 200 mm cracks requires approximately 13 MPa of expansion pressure when (σ1, σ3) = (5 MPa, 0 MPa), compared to 100 MPa when (σ1, σ3) = (20 MPa, 15 MPa). Figures 15 and 17 show that P gradually increases with respect to σ3 when cracks of the same length are generated. Compared with the numerical results under the σ1 = σ3 conditions, it can be concluded that σ3 is the main factor affecting the crack initiation pressure P0, while both the stress ratio σ3/σ1 and notch direction control crack direction.

The corresponding P-l linear fit equations are listed in Table 4. Figure 18 shows the variations of P0 and al versus σ3 when σ1 − σ3 is kept constant at 5 MPa. According to Figure 18, P0 increases with σ3 under all notch direction conditions; P0 increases with θ under all stress conditions. The SCA expansion pressure required for crack initiation that lowers as the notch direction tends toward the σ1 direction. When θ = 0° and θ = 45°, the trends of P0-σ3 curves are approximately linear, but when θ = 90°, the trend flatlines as σ3 > 12.5 MPa. Under all notch conditions, al is inversely proportional to σ3. The results suggest that crack generation becomes more difficult as σ3 and θ increase.

6. Conclusions

To better understand the mechanism of SCA-induced crack propagation in deep underground, crack initiation and propagation under different stress conditions were modelled using a modified cohesive element method. A new traction-separation law for describing rock compressive shear strength was proposed. Crack length and crack direction in bidirectional isobaric and unequal stress fields were analyzed in detail. According to the numerical results, the following conclusions are drawn:(1)The crack morphology modelled using the modified cohesive element method is similar to that of the experimental results, indicating that the method proposed in this study can model SCA-induced crack initiation and propagation.(2)The crack initiation pressure P0 and the incremental ratio of crack length to SCA expansion pressure al can be used as useful indicators to evaluate the difficulty to break rock in deep underground. High P0 and low al make rock breaking more difficult.(3)The maximum expansion pressure of SCAs depends on the elastic modulus and geostress field of the rock and therefore should be measured under conditions similar to deep underground before designing SCA borehole spacing, which should be less than 220 mm at 600 m deep with a maximum SCA expansion pressure of 100 MPa. To break rock more effectively, borehole spacing should be reduced under higher stress condition.(4)σ3 is the main factor that affects the crack initiation pressure P0. Both the stress ratio σ3/σ1 and notch direction control the direction of crack propagation. When the notch direction is oblique with the σ1 direction at 45°, the angle between the crack and the σ3 direction θc decreases linearly as σ3/σ1 increases, and their linear relationship is

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declared that no conflicts of interest exit in this work.


This paper was financially supported by the National Natural Science Foundation of China (51604093 and 51474096), the Program for Innovative Research Team at the University of Ministry of Education of China (IRT_16R22), Program of China Scholarship Council (201708410138), the Doctoral Research Fund Project of Henan Polytechnic University, China (B2017-42), and the Scientific and Technological Key Project of Henan Province (172107000016).