Abstract
Time-dependent behavior of subcritical crack growth is one of the main characteristics in rocks. The double-torsion test is commonly used to study the slow crack growth behavior of brittle and quasi-brittle materials. However, double-torsion specimen is difficult to processing, the process of the laboratory test is irreversible, and the current numerical simulation is difficult to consider the time-dependent behavior, and so on. In view of all above problems, an idealized particle model was built, and the crack was identified in this article, based on the theory of particle flow. The numerical model was built using Particle Flow Code in 3 Dimensions, and the macromechanical and micromechanical parameters of the model were calibrated. The process of the macroscopic crack propagation and its evolution were analyzed. The intrinsic relations with the load, the displacement, and the time were established. The results show that the Particle Flow Code in 3 Dimensions can reproduce the time-dependent behavior of subcritical crack growth in double-torsion test. And, the peak and the law of the curves are in good agreement with the laboratory test results. Therefore, the Particle Flow Code in 3 Dimensions numerical simulation can be used as a new effective method to reveal the slow crack growth behavior, to get the relevant parameters such as V, KI, and KIC, and to build the relationship between V and KI. The results of this article will have some reference value for the simulation and application of double-torsion test.
Keywords
Introduction
The double-torsion (DT) test is commonly used to characterize the slow crack growth behavior of brittle materials. DT specimen was adopted first by Evans 1 and Williams and Evans. 2 It can be regarded as a symmetrical system of two independent elastic torsions. This test consists in bending one of the edges of a “V” notched plate put on four rollers at each corner (Figure 1). In such setup, one face is loaded in tension whereas the other one is loaded in compression. After a pre-cracking procedure, the DT test can be used for the determination of the fracture toughness or the analysis of slow crack growth. The toughness measurement is conducted at very high loading rate and the record of the maximum load before failure and pre-crack size allows obtaining KIC. Slow crack growth curves as a function of the mode I stress intensity factor can be measured with a single specimen using different methods: either by a load relaxation test via a compliance calibration curve or by loading the specimen at different loads and measuring the crack growth rates with an optical microscope.

(a) Specimen shape, (b) cross-section in width, and (c) cross-section in length.
According to the relation between the stress intensity factor and the strain-energy release rate, and the applied load on torsion, and the relation of Poisson’s ratio and shear modulus of the rock, the stress intensity factor can be derived 3
where KI is the stress intensity factor, N m−3/2; E is Young’s modulus, Pa; g is the strain-energy release rate, N m−1; F is the applied load, N; W is the width of specimen, m; Wm is the moment arm of the torsion, m;
According to the relation between the compliance of specimen and the length of crack, the relation between the change rate of compliance and subcritical propagation of crack can be set up. Test indicates that the change rate of compliance can be obtained through the change rate of displacement under invariable loads or the change rate of loads under invariable displacement. Under the invariable displacement condition, the subcritical crack velocity of rock equals
where a is the crack length, m; y is the constant displacement imposed on the loading point, m.
Equation (2) shows that if the size of specimen and the displacement on loaded point are known, then the subcritical crack velocity is relate to the load relaxation rate in current loads under the given displacement state. Therefore, we can establish the relation between stress intensity factor KI and subcritical crack velocity v through the test.
Subcritical crack growth plays an important role in evaluating the long-term stability of structures in rocks. Wu et al. 4 and Zhao et al.5–8 carried out a series of theories for crack development. Henry et al., 9 Atkinson, 10 and Swanson 11 found that the characteristics of subcritical crack growth can be described by a relationship between the stress intensity factor and the crack velocity. At present, the KI-V relationship and KIC of rock are measured by the method of constant displacement relaxation of DT specimen. Many scholars12–15 investigated experimentally subcritical crack growth in granite using the DT test. TY Ko and J Kemeny 16 present the results of studies conducted to validate the constant stress-rate test for determining subcritical crack growth parameters in rocks, compared with the conventional testing method, the DT test. Nara et al. 17 studied the influence of relative humidity on fracture toughness of rock using DT test. P Leplay et al. 18 designed a dedicated experimental setup to perform in situ tests in an X-ray computed tomography scanner on a porous brittle ceramic and to analyze the three-dimensional images with digital volume correlation. Attia et al. 19 investigated the fracture behavior of DT specimens using the finite element method taking no account of the time effect. Zhang 20 and Gao and Li 21 simulated the crack propagation with equivalent cohesive zone method based on virtual internal bond theory. Shyam and Lara-curzio 22 summarized the analytical compliance and finite element stress analyses of the DT test specimen. And, they concluded based on the review that standardization of this test method is required in order to make it more practicable.
Beyond its wide use in the literature, the DT test presents some questionings: (1) the brittle rock is processed into an extremely thin plate with a thickness of 3–5 mm, and the process is extremely difficult. (2) The heterogeneity of the rock results in a larger discrepancy. (3) The test loading method is complex and irreversible. (4) Subcritical crack growth is one of the main causes of time-dependent behavior in rocks. The existing numerical simulation results did not involve time.
Fortunately, with the rapid development of numerical calculation method, numerical simulation is solving more and more complex rock mechanics and engineering problems.23–26 Particle Flow Code in 3 Dimensions (PFC3D) was established, to research the mechanical characteristic of rock material through the movement and interaction of the particulate media, which is widely used in the study of discontinuous deformation and failure characteristics of rock materials. In this article, based on the theory of particle flow, the command flow of FISH language in PFC3D will be used to simulate the DT test. And, the feasibility of DT numerical test by PFC3D will be verified through comparison with indoor experiment.
Crack idealization and identification
Parallel bond behavior
A parallel bond approximates the physical behavior of a cement-like substance lying between and joining the two bonded particles. Parallel bonds establish an elastic interaction between particles that acts in parallel with the slip or contact-bond constitutive models. Parallel bonds can transmit both force and moment between particles, while contact bonds can only transmit force acting at the contact point.
A parallel bond can be envisioned as a set of elastic springs uniformly distributed over either a rectangular (in two-dimensional (2D)) or a circular (in three-dimensional (3D)) cross section lying on the contact plane and centered at the contact point, as shown in Figure 2. These springs act in parallel with the point-contact springs (which come into play when two particles overlap). The behavior of the parallel bond springs is similar to that of a beam whose length, L in Figure 2, approaches zero. Relative motion at the parallel bonded contact (occurring after the parallel bond has been created) causes axial- and shear-directed forces (T and V, respectively) as well as a moment (M) to develop. In PFC3D, a twisting moment (Mt) may also develop. The maximum normal and shear stresses (σ and τ, respectively) carried by the bonding material can be written as
where A and I are the area and moment of inertia, respectively, of the parallel bond cross section; J is the polar moment of inertia of the parallel bond cross section; and positive T indicates tension.

Parallel bond idealization (a), forces carried in the 2D bond material (b), and 3D bond material (c).
Crack identification
Microcracks in a PFC3D synthetic specimen may only form between bonded particles. Particle contact parallel bond model assumes that there is an elastic cementing material between the two particles that connects the two particles together. The relative movement of the contact position results in the bond material generating force and moment.
One mechanism (Cundall et al. 27 ) for the formation of compression-induced cracks is shown in Figure 3(a), in which a group of four circular particles is forced apart by axial load, causing the restraining bond to experience tension. When the tension is large enough, the bond between the particles is destroyed, which indicates the formation of cracks (see Figure 3(b) and (c)). Similar crack inducing mechanisms occur even when different conceptual models for rock microstructure are used. The success of the Particle Flow Code (PFC) material at reproducing the behavior of hard rock can be attributed directly to its ability to generate compression-induced tensile cracks.

(a) Contact form the parallel bond, (b) the parallel bond is broken, and (c) crack formation.
With the destroyed bond between the particles increasing, the material particles move far away from each other, and their interactions become weaker and weaker, as shown in Figure 4. The bond becomes weaker and weaker. Therefore, under the condition that the displacement loading on the specimen remains constant, the required external forces that keep the cracks continue to grow will become smaller and smaller.

(a) Prefabricated crack, (b) crack propagation, and (c) crack penetration.
PFC3D implementation
Geometric model
The model is generated as follows: first, the basic dimensions of the model (150.0 mm × 60.0 mm × 4.80 mm) are defined by the six-sided wall (constrained particle displacement) and then the particles of a given radius are randomly generated within the size range. The total number of particles in this model is 104,138, and the particle radius is Gaussian distribution of Rmin ~ Rmax. And, in order to reduce the computational complexity and make the show of crack growth more accurate, the block modeling method is adopted in numerical model as shown in Figure 5. Zone 1 is the main zone in which crack will grow, hence the small particles are adopted (Rmin = 0.10 mm, Rmax = 0.25 mm). Zones 2 and 3 are the zones of little influence, hence the bigger particles are adopted (Rmin = 0.02 mm, Rmax = 0.40 mm for zone 2, Rmin = 0.30 mm, R max = 0.60 mm for zone 3). The cycle command was written to eliminate the internal non-uniform stress of model, in order to make the model reach to a stable state.

Calculating model of double-torsion specimen: (a) top view and (b) bottom view.
After generating a stable numerical model, the six confined walls were removed; then, axial groove and V-type notch of specimen were cut by the range and delete command as shown in Figure 5; then, the fix command was used to constrain the bottom displacement at the four corners of the specimen, in order to achieve the purpose of forming point support. Finally, at the top of the model, two square load plates (walls) with a length of 1 mm are defined, and the model is finished.
Model parameter calibration
Micromechanical parameters in PFC3D are used to characterize the mechanical properties of particles and bonding. Before the numerical model begins to calculate, the micromechanical properties of the model must be given. However, it is difficult to obtain the micromechanical parameters of the model directly from the laboratory test. So, the corresponding parameter calibration is needed.
The parallel bond parameters can more reasonably simulate the destruction of brittle rocks by cementing between particles. So, the parallel bond parameters were used in this article. And, the uniaxial compression numerical test was carried out. Comparing with the macromechanical parameters obtained by the laboratory test, the final micromechanical parameters (in Table 1) are obtained by continuously adjusting the numerical parameters in the numerical test, and the results are consistent with the laboratory test and applied to the actual calculation model. The macromechanical parameters from the numerical test, including elastic modulus E, uniaxial compressive strength, and Poisson’s ratio, are similar to those obtained from the laboratory test (Figure 6).
Micromechanical parameters.

Uniaxial compression curve in numerical simulation.
Loading method
In the laboratory DT test, the four corners of the bottom support platform are set with four steel balls of diameter 5 mm, see Figure 1(a). The specimen is placed on the four steel balls. The front two balls are subjected to the pressure of the loading head, and the tail two pieces only serve to support the specimen. Loading head is also in contact with specimen through two steel balls. The diameter of the two balls is 4 mm. All above steel balls are painted butter. The role of ball bearing support is to ensure the force for the point load, and to reduce the friction between the bearing and the specimen. The test was carried out on an Instron material servo MTS. And, the loading displacement was measured with a high-precision axial extensometer as shown in Figure 7.

DT test loading system.
Before the DT test is started, the specimen should be pre-cracked. The method is as follows: after the servo loading, head is adjusted to just touch specimen, test starts to load with a constant displacement rate of 0.0008 mm/s. 28 And, the loading is recorded at the same time. When the load does not rise with time, pre-crack is completed. At this moment, the loading should be stopped, and the pressure head will be lifted. In PFC3D, the HISTORY command was used to monitor the loading force. When the loading force does not increase with time, the cracks as shown in Figure 8, appearing on the model, will be thought as pre-crack.

Loading–time curve of pre-crack.
After the completion of the pre-crack test, the load will be unloaded first, and reloaded to 95% of the pre-crack load. 29 Then, the corresponding displacement will be keep as a constant. And, the relaxation test begins. The loading rate is shown in Figure 9. The final displacement of the laboratory test and numerical simulation was 0.2 and 0.215 mm, respectively.

DT loading displacement–time curve.
The model is applied by the loading plate shown as Figure 5(a). The loading plate is a rigid body, which is brought into contact with the corresponding particles, when the loading plate is pressed against the model at a fixed speed. When the contact occurs, according to the force–displacement law, the contact force will produce between the particles and the wall. Then, according to Newton’s law of motion, the particles will be displaced under contact force, and new contact will occur with the surrounding particles and the walls, resulting in new contact force. This cycle (Figure 10) realizes the loading for the model. And, the force–displacement law can be written as
where

Calculating cycle of model.
The law of motion can be written as
where
Results and discussion
Figure 11 shows the results of the DT laboratory test, and Table 2 shows the details about the DT specimen in this experiments. Figure 12 shows the results of the DT numerical test, in which Figure 12(a) shows the pre-crack diagram and Figure 12(b)–(f) shows simulation of crack propagation diagram from 0.6 to 9 min. Figure 13 shows the crack length–time curve and crack growth velocity–time curve in the DT test with constant displacement relaxation. It can be seen from the figures that the crack propagates along the guide groove in the DT numerical test, and the crack growth velocity gradually decreases as time increases. The phenomenon of crack propagation in the numerical test is consistent with that in the case of laboratory test, which proves that PFC3D is feasible in simulating DT test and crack propagation.

Photo of the damaged specimens in laboratory test.
The property information of double-torsion specimen.

Crack propagation process in the DT numerical test: (a) 0.3 min (pre-crack), (b) 0.6 min, (c) 1.5 min, (d) 3 min, (e) 6 min, and (f) 9 min.

Crack length–time curve, and crack growth velocity–time curve.
Various engineering examples related to rock mass show that the deformation and failure mechanism of rock mass are closely related to the initiation, propagation, and penetration of cracks in rock mass. The rapid propagation of cracks in rock and the fracture of structure usually occur after the subcritical crack growth to a certain extent, which leads to the time-dependent relationship between the instability of geotechnical engineering and the propagation of cracks in rock.
Figure 14 shows the load–time curve in the DT test with constant displacement relaxation. It is not difficult to find that the load decreases gradually with increase in time and then gradually stabilizes. The peak load of the numerical simulation curve is 165 N, while the laboratory test is 158 N, and the peak difference is 4.2%. The stable load of the numerical simulation curve is 116 N, while the laboratory test is 110 N, and the stable load difference is 5.1%. In the constant displacement relaxation test, the relationship between load and time is in accordance with the following equation
where a, b, c, and I are constants associated with material and loading method.

Load–time curve in the DT test with constant displacement relaxation.
Compared with the laboratory test results, the numerical simulation is in agreement with the laboratory test results, including the crack propagation law, the displacement–time curve, and the load–time curve. It is shown that the numerical test can better reflect the failure process and the microscopic fracture mechanism of the DT specimen. The reason for the difference between the numerical simulation and the laboratory test may be the mineral composition of the rock, the particle size, and the mechanical properties of the uneven.
We can find from Figure 15 that the curves of lgKI-lgV can be well linear, but the data measured in the relaxation test keep high discretion. After researching carefully, we think that the high discretion of data is caused by mineral composition, particle size, and mechanical property of the rock itself. In order to reduce the dispersed influence to results, the tests of the same kind of rock should be done many times using the constant-displacement load-relaxation method and return the average values of the regressed coefficients as the characteristic parameters of the subcritical crack of the rock.

Double logarithmic space lgKI-lgV.
Conclusion
The DT technique was successfully applied to rock for the study of crack propagation. The technique is able to generate stable crack growth at constant driving displacement. It is suitable for the measurement of fracture characteristics of brittle and quasi-brittle materials, which experience non-linearity due to microcracking in the fracture process zone and crack bridging in the wake.
When the tension is large enough, the bond between the particles is destroyed, which indicates the formation of cracks. Similar crack inducing mechanisms occur even when different conceptual models for rock microstructure are used. The success of the PFC material at reproducing the behavior of hard rock can be attributed directly to its ability to generate compression-induced tensile cracks. With the destroyed bond between the particles increasing, the material particles move far away from each other, and their interactions become weaker and weaker. Therefore, under the condition that the displacement loading on the specimen remains constant, the required external forces that keep the cracks continue to grow will become smaller and smaller.
The PFC3D can reproduce the time-dependent behavior of subcritical crack growth in DT test. And, the peak and the law of the curves are in good agreement with the laboratory test results. Therefore, PFC3D numerical simulation can be used as a new and effective method to study the slow crack growth behavior of rock and to get the relevant parameters, such as subcritical crack growth velocity V, stress intensity factor KI, and fracture toughness KIC, and to set up the relationship between V and KI.
Footnotes
Handling Editor: Kenneth Loh
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research work was funded by the National Natural Science Foundation of China (51874112 and 51774107); the State Engineering Laboratory of Highway Maintenance Technology, Changsha University of Science and Technology (kfj170108); the Opening Project of State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology (KFJJ17-12M); the Open Program of State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology (SKLGDUEK1406); the Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences (Z013010).
