Numerical study on tensile failures of heterogeneous rocks

Ying Xu, Wei Yao, Kaiwen Xia

a State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin, 300072, China

b Department of Civil and Mineral Engineering, University of Toronto, Toronto, Ontario, M5S 1A4, Canada

Keywords:Tensile strength Flexural tensile strength Realistic failure process analysis (RFPA)Homogeneity index Non-local failure Characteristic length

A B S T R A C T Many experimental results have demonstrated the apparent discrepancy of a rock material between its flexural tensile strength measured using various bending methods and its tensile strength measured using direct tension method or Brazil disc(BD)method.To understand the physical mechanism for such discrepancy,numerical simulation using the realistic failure process analysis(RFPA)is carried out in this work to simulate the tensile failure of heterogeneous rocks. Direct tension and semi-circular bend(SCB)tests are simulated using RFPA for rock materials with different levels of inhomogeneity, which is characterized by the homogeneity index of the Weibull distribution used in RFPA. The numerical results show that the discrepancy in the tensile strength values is caused by the inhomogeneity of the rock material. Furthermore, non-local failure criterion is adopted to calculate the characteristic length of the rock materials used in the simulation. It is shown that below a certain value of the homogeneity index,both the characteristic length and discrepancy between two types of tensile strengths of rock decrease with increase of the homogeneity index up to a critical value, at which the discrepancy disappears and the rock material is essentially homogeneous.

1. Introduction

As a typical brittle material, rock has a much lower tensile strength in comparison with its compressive strength (Sundaram and Corrales, 1980; Cao et al., 2017; Gao et al., 2017; Liu et al.,2017). The tensile strength of rock is therefore a significant material parameter with applications in many engineering and geological fields. However, in many engineering practices, the tensile strength of rock is mostly neglected in design or construction stage(Chen and Hsu, 2001). Interestingly in many numerical analyses,the tensile strength of rock was assumed to be 10% of the compressive strength. As a matter of fact, the tensile behavior of rocks varies due to formations, loading rates, locations, and microscopic structures of rocks. The inaccuracy of tensile strength value may lead to inappropriate design or safety assessment of structures and rock masses.

As an important mechanical property, the tensile strength of rocks can be measured in laboratory by various direct and indirect methods(ISRM,1978).Since the tensile strength of rock is defined as the failure stress of rock in pure uniaxial tensile loading (Dai et al., 2010a), direct tension method is most suited for measuring the rock tensile strength. However, in practice, the ideal uniform stress state in the rock specimen is rarely provided in the direct tension test, and stress concentration induced by instrumental misalignment may also lead to the doubtful tensile strength value. Several indirect experimental methods thus have been proposed to measure the tensile strength of rocks, such as the Brazilian disc(BD)test(Dai et al.,2010b;Zhou et al.,2014),the ring test (Chen and Hsu, 2001; Chen et al., 2008), and the semicircular bend (SCB) test (Dai et al., 2008). In addition, dynamic tensile strength of rock was also measured through the same type of specimen as in the static tensile strength methods. Dai et al.(2010a) measured the dynamic flexural tensile strength of rock and other brittle solids with the SCB method and a modified split Hopkinson pressure bar (SHPB) system. Zhao and Li (2000)measured the dynamic tensile properties of granite with the BD and three-point bend methods. These indirect methods can provide efficient and convenient ways for determining the quasistatic and dynamic tensile strengths of rocks.

It is noted that the failure of specimens for SCB and direct tension tests is both caused by the tensile stress. However, the flexural tensile strength from the SCB method is higher than the tensile strength from the direct tension and BD methods for the same material under the same loading rate(Dai et al.,2010a).This phenomenon was also observed in the static tensile strength tests(Lajtai, 1972; Carter, 1992; Van de Steen and Vervoort, 2001). In these tests, the effect of stress gradient on the tensile fracture initiation under static bending condition was studied using different specimen shapes and sizes. The results from these tests revealed that the stress for tensile fracture initiation depends on both the maximum tensile stress and the stress gradient along the fracture path. In such case, the equality between the maximum tensile stress concentration and the uniaxial tensile strength is a necessary but not an adequate condition for fracture initiation under static bending condition.In other words,these static tensile strength tests indicated that there exists a discrepancy between the flexural tensile strength and the direct tensile strength. In addition, Zhao and Li (2000) also reported a smaller value obtained from the BD test,compared with that from the three-point bending test.

To understand the discrepancy between the rock tensile strengths from direct and indirect tension tests, the non-local failure model has been used (Lajtai, 1972; Carter, 1992; Van de Steen and Vervoort, 2001). The non-local failure model is an extension of the traditional failure model (Pinsan,1991; Seweryn and Mróz,1995; Isupov and Mikhailov,1998). The difference between the dynamic flexural tensile strength and the tensile strength of Laurentian granite(LG)was successfully explained by the non-local failure approach(Dai et al.,2010a).To date,there are three non-local failure criteria, i.e. minimum stress fracture criterion (MSFC), average stress fracture criterion (ASFC), and fictitious crack fracture criterion (FCFC). The results showed that in most cases,ASFC is the most suitable(Isupov and Mikhailov,1998;Van de Steen and Vervoort,2001),which is actually the commonly used non-local failure theory (Lajtai, 1972; Carter, 1992; Van de Steen and Vervoort, 2001). Although there are some attempts to explain the discrepancy between the tensile strengths from indirect and direct tests using the non-local failure theory,the physical meaning of the characteristic length in the theory and thus the physical mechanisms of the discrepancy are still unclear.

Rock has been recognized as a heterogeneous material for a long term. To analyze the statistical variation in the bulk failure strength of heterogeneous materials, Weibull (1951) adopted the statistics of extreme to characterize the local failure strength using the probability distribution function. This concept has been adopted by the realistic failure process analysis (RFPA) method (Chen et al.,1997; Tang and Kaiser, 1998; Tang et al.,1998; Tham et al.,2001). As shown by Wong et al. (2014), RFPA can reproduce the rock strength, given properly chosen Weibull distribution parameters. Therefore, in this work, we use RFPA to investigate the physical reason behind the discrepancy between the flexural tensile strength from the SCB method and the direct tensile strength of rock materials. The results also shed light on the physical basis of the non-local failure theory.

2. Overview of realistic failure process analysis (RFPA)

Briefly, RFPA is a finite element code used to model the observed evolution of damage and induced seismicity by the progressive failure.In RFPA,the collapse in quasi-brittle materials due to the progressive failure allows linear elastic elements to fail in a brittle manner as long as their peak strengths are reached(Tang and Kaiser, 1998; Zhu and Tang, 2006; Zhu et al., 2012).Many researchers have employed RFPA to investigate the compressive and tensile failures of rock material under static and dynamic loadings (Chen et al.,1997; Wong et al., 2006; Zhu and Tang, 2006; Tang et al., 2007; Zhu et al., 2012, 2018; Dai et al.,2015; Meng and Liu, 2018; Yu et al., 2018). With the digital image process (DIP) techniques, the RFPA has also been used to interpret the effect of grain size on the uniaxial compressive strength (Yu et al., 2018).

Microscopically, rock is a grain-scale heterogeneous material with various minerals.In RFPA,the solid or structure is assumed to be composed of many mesoscopic elements with the same size.In order that the continuum mechanics and linear elastic finite element method(FEM)can be used to simulate the failure process of rocks involved in nonlinear and discontinuum mechanics, the Weibull probability function is introduced to characterize the local failure strength in RFPA (Chen et al.,1997; Tang and Kaiser,1998;Tang et al.,1998; Tham et al., 2001), which can be written as follows (Chen et al.,1997):

where A is a property parameter (e.g. density, Young's modulus,Poisson's ratio,etc.), and Acis the average value of A.In this study,the density,Young's modulus and Poisson's ratio were specified to each element according to the Weibull distribution in the simulation (Zhu and Tang, 2004; Lei and Gao, 2019).

Mathematically, the parameter m represents the shape of the distribution function and thus controls the distribution of material properties as the average value is fixed in the Weibull distribution.Consequently, the parameter m defines the degree of material homogeneity. It is thus referred to as the homogeneity index. Fig.1 illustrates that the higher the value of m, the more concentrated the distribution of A. Therefore, it is apparent that the larger the value of m, the more homogeneous the material (Tang and Kaiser,1998; Tang et al.,1998).

In RFPA,it is assumed that the material constitutive response is brittle and nonlinear. The brittle response here means that the material breaks without significant plastic deformation at failure.The element fails when the stress exceeds its strength,resulting in a microcrack in the material.The constitutive model used in RFPA is illustrated in Fig. 2. One can find the details of damage constitutive relation in the literature(Tang et al.,2002;Zhu et al.,2005,2006; Li et al., 2009). It should be emphasized that the non-local failure model is not used as a failure criterion in RFPA and it will be shown later that the non-local failure emerges from the simulator naturally using this simple damage constitutive relation shown in Fig. 2.

When the damage accumulates to a certain value, rock failure occurs (Chen et al.,1997; Tham et al., 2001; Hao et al., 2002; Liu et al., 2014, 2015). According to the Lemaitre theory of strain equivalence, one has the following formula:

Fig.1. Distribution of material parameter A with different homogeneity indices m.

Fig. 2. The elastic damage constitutive law used in RFPA. σ1 and σ3 are the uniaxial compressive stress and tensile stress in an element, respectively; ft0 and ftr are the uniaxial tensile strength and residual uniaxial tensile strength,respectively;fc0 and fcr are the uniaxial compressive strength and residual uniaxial compressive strength,respectively.

σ = E(1 D)ε (2)

where σ is the element stress, E is the elastic modulus, ε is the strain, and D is the damage parameter. It is evident that D = 0 corresponds to undamaged state in the rock and D=1 to complete failure of the rock. When the mesoscopic element is in multiaxial stress states,the constitutive relation of the element is illustrated in Fig.2.Initially,no.damage exists,i.e.D=0＞εt0，whereis the equivalent principal strain, and=in which ε1, ε2and ε3are the three principal strains andis a function defined as〈x〉 = x if x 0 and〈x〉 = 0 if x ＜0).When the failure criterion is satisfiedεtu)，damage occurs in the element(D = 1). Between these two conditions, the damage parameter D evolves following the equation D = 1 λεt0/ε (εtu＜εt0)(where λ is the residual tensile strength coefficient)in RFPA model.The principle of the damage parameter evolution was detailed and discussed in Zhu et al.(2005).Furthermore,the calculation process of RFPA method is illustrated in Fig. 3 (Yang et al., 2012).

3. Simulation of direct tension test

In numerical simulation of direct tension test using RFPA, a cylinder rock specimen with 50 mm in diameter and 100 mm in length was employed (Fig. 4). The ratio of specimen length to diameter is 2:1, which is the same as that used in experiments(Okubo and Fukui, 1996; Lin and Takahashi, 2008). Square elements are utilized in the simulation. The element size for all models is 0.25 mm and the number of elements is 200 400 = 80,000. Table 1 lists the basic physico-mechanical properties of the rock material (Dai et al., 2010a). These properties were experimentally measured from LG, which is a wellstudied brittle rock material (Douglass and Voight,1969; Nasseri and Mohanty, 2008; Dai et al., 2010a; Lan et al., 2010; Li et al.,2018a, b). Here, the density, Young's modulus, Poisson's ratio and ratio of tensile to compressive strength were used in the numerical simulation.

Fig.5 shows the stress-displacement curves obtained from the direct tension tests. Damage occurs in the rock before the stress reaches the ultimate tensile strength.From Fig.5,we can obtain the macroscopic tensile strengths of rock specimens with different values of m.It is noted here that even though the basic parameters are the same, we could have different rock materials by changing the value of m in RFPA.

Fig. 3. The calculation flowchart of RFPA.

Fig. 4. Schematic diagram of direct tension test.

Fig. 5. The stress-displacement curves from direct tension simulation.

As seen in Table 2, the direct tensile strength of material increases with m (the mean elastic modulus and strength of the material remain the same),and when the value of m is larger than 33, the direct tensile strength reaches its maximum value of 12.78 MPa, matching the experimental result (12.8 MPa) (see Table 1).

4. Simulation of semi-circular bend test

SCB test was conducted to obtain the flexural tensile strength of various materials (Huang et al., 2006; Dai et al., 2008). The dimension of the SCB specimen in the modeling here is the same as that used in the experiments(Dai et al.,2008).The square element is also employed in the SCB simulation. In the SCB model, the element size is 0.25 mm and the number of elements is 20,150.The density, Young's modulus, Poisson's ratio and ratio of tensile to compressive strength of LG in Table 1 are also used in the SCB numerical simulation. As shown in Fig. 6, the span between two supporting pins is S = 21.8 mm.

The direction of load P along radius is perpendicular to the bottom of the half disk. During loading, failure starts at the failure spot O on the specimen due to bending. According to the dimensional analysis,the tensile stress in the horizontal direction at O can generally be expressed by the following equation(Dai et al.,2008):

where P(t) is the time-varying load, B is the thickness of the specimen,and R is the radius of the disc.The dimensionless stress Y(S/(2R)) is a function of dimensionless distance S/(2R) and can be calibrated using the finite element analysis (Adamson et al.,1996;Huang et al., 2006; Dai et al., 2010a) as follow:

In this study,Y is equal to 5.14 based on the configuration of the SCB specimen (R = 20 mm, S = 21.8 mm, and B = 10 mm). The flexural tensile strength σfis taken as the maximum tensile stress in the history of σ(t).

As shown in Fig.7,the maximum vertical load increases with the increasing m value, indicating that the maximum tensile stress in the SCB specimen enhances with the value of m. Here, in order to compare the stresses in both direct tension and SCB tests, the displacement (i.e.the absolute value of the displacement between the upper and bottom loading points along the loading direction in both direct tension and SCB tests)instead of strain is used in the xaxis of both Figs. 5 and 7. The damage (e.g. microcracks) accumulates around the failure spot O with the increase of the load.As the load increases, the crack reaches a critical length and its growth becomes unstable, and then the rock specimen fails.

Fig.8 shows the distribution of the stress along the loading axis in the specimen when m = 3. It can be seen that the maximum tensile stress is reached around the failure spot O(Fig.6)despite of the significant variation of the stress as a result of rock inhomogeneity (Fig. 8). We denote x as the distance between the point along the loading axis and the failure spot O. Based on the simulations, when 0 mm ＜x ＜15 mm, the specimen is under tensile stress along the loading axis,but when 15 mm ＜x ＜25 mm,the specimen is under compressive stress.In two-dimensional(2D)view, there is a tensile zone enclosing the failure point O. Outsidethis zone, the rock is under compression. The boundary between the two zones goes through the two supporting points and a point along the loading axis (Fig. 8). According to Eq. (3), the flexural tensile strength can be calculated and the results are listed in Table 3. It can be seen that when m 33, the flexural tensile strength reaches its maximum value of 12.83 MPa, which matches well with the experimental tensile strength of the rock.In this case,the rock material is considered as a homogenous material(Table 1).

Table 1 Basic physico-mechanical properties of the rock material (Dai et al., 2010a).

Table 2 Direct tensile strengths calculated using RFPA.

Fig. 6. Schematics of the SCB test.

Fig. 7. The stress-displacement curves from SCB simulation.

Fig. 8. Distribution of stress along the loading axis in the SCB specimen (m = 3).

5. Analysis of simulation results

The simulation results are plotted in Fig. 9. As shown in the figure, both the direct tensile strength and the flexural tensile strength increase with the increase of m. According to the mathematical meaning of m, the material is very heterogeneous when the m value is small, but when it is large enough(＞33 in the current case), the material is essentially homogenous. For different m values, the model describes different rock materials. It is thus only meaningful to check the strength ratio σf/σtfor comparison purpose. It can be seen that when m = 2,the material is extremely heterogeneous, and the strength ratio is the maximum. This ratio dramatically decreases with the increase of m. When the value of m reaches the critical value of 33, the ratio σf/σtbecomes 1 and the material is basically homogenous.

The above result indicates that the microscopic properties of material (m) significantly influence rock macroscopic properties.Our simulation results can well reveal the discrepancy between the flexural tensile strength and the tensile strength, i.e. the rock inhomogeneity can be characterized by the parameter m. Also,this discrepancy is evident when the material is extremely heterogeneous. This has also been proven by the experiments on three typical granitic rocks featured with different grain sizes(Yao et al., 2018). The experimental results (Yao et al., 2018) revealed that rocks with greater grain size (which can be considered as greater heterogeneity for rocks) have more discrepancy in the tensile strength between the bending and BD tests. This is consistent with the simulation results in this study.

5.1. Material characteristic length in the non-local failure model

The non-local fracture criteria can be generally expressed as follows (Isupov and Mikhailov,1998; Guo et al., 2009):

where f(sij(y))a function of stress tensor sij,and σtcis is a material constant.

Some general forms of the most popular non-local fracture criteria for 2D cases were given by Isupov and Mikhailov (1998).This failure criterion can be expressed as

Table 3 Flexural tensile strengths calculated using RFPA.

Fig. 9. The flexural tensile strength compared with the direct tensile strength.

where the center of a local polar coordinate system(ρ,θ)of a body is at the point y; η(θ) is a unit vector with angle θ to the coordinate axis; σθis the stress component in this coordinate system; the constant σtchere denotes the strength of a body under uniform traction,i.e.intrinsic tensile strength;and d1is a material constant with length measure, i.e. the material characteristic length.

If the maximum value of θ in Eq.(6)is θ0,the ASFC equation can be simplified as follows:

This criterion has been used to reconcile the difference between static flexural tensile strength and uniaxial tensile tests (Lajtai,1972; Carter,1992; Van de Steen and Vervoort, 2001) and the one between the dynamic flexural tensile strength and dynamic tensile strength (Dai et al., 2010a).

According to Eq. (7), the failure of material occurs when the average local stress σ over a distance d1along the potential fracture path reaches the tensile strength(Dai et al., 2010a):

where σ is the tensile stress along the potential fracture path.

We also calculated the flexural tensile test using a linear elastic model as we did in the analysis for real laboratory tests.The tensile stress gradient along the potential fracture path (i.e. the loading axis through O)of the SCB specimen was calculated with the finite element analysis(as shown in Fig.10).Exponential fitting equation is used to describe the normalized stress gradient as follows:

where σmis the tensile stress at the failure spot (the maximum tensile stress in the specimen). Sinceequals the direct tensile strength σt,and σmequals the flexural tensile strength σfat failure,d1can be calculated by substituting Eq. (9) into Eq. (8).

Fig. 10. Normalized tensile stress along the prospective fracture path in the SCB sample (Dai et al., 2008). x is the distance from the failure point O along the prospective fracture path (see insert).

The values of σtand σfare taken from Tables 2 and 3, respectively.The results of d1are listed in Table 4.It can be seen that d1has a one-to-one relation to the parameter m and it decreases with the increase of m until 33,after which d1reaches its minimum value of 0.019. This result clearly verifies that the material characteristic length d1indirectly measures the material heterogeneity.

5.2. Discussion

In order to further understand the correlation between m and d1as shown in Table 4, a parameter X is defined as follows:

where X is called the inhomogeneity index,and M is the minimum value of m when the material becomes homogeneous (m = 33).According to the definition of X, the domain of X is ( ∞， +∞).

Fig.11 shows the relationship between the inhomogeneity index and the material characteristic length d1. The curves of d1can be divided into three stages and each stage can also be divided into two parts. With the increase of X, the material changes from homogeneous to heterogeneous. In stage I, the material is relatively homogeneous; from stage II to III, the material is more and more heterogeneous. In each stage, with increase of the inhomogeneity index, the characteristic length of material is first constant (called‘stable state') and then increases.

Table 4 Characteristic length for rocks based on the ASFC criterion.

Fig.11. Relationship between the inhomogeneity index and the characteristic length of material.

The relationship between the inhomogeneity index and the material characteristic length (Fig.11) can be written as

d1= k1(1 k2)f1(X)+k2(1 k1)f2(X)+(1 k1)(1 k2)f3(X)(11)

The three terms in Eq. (11) represent the three stages of d1(stages I,II and III).The prediction of Eq.(11)is plotted in Fig.11.It can be seen that the proposed equation represents well the trend of the data.

According to Eqs. (10)-(12), the relationship between the homogeneity index and the characteristic length of material can be written as follows:

d1= k1(1 k2)f1(m)+k2(1 k1)f2(m)+(1 k1)·(1 k2)f3(m) (m 3) (13)

Fig.12. Relationship between the homogeneity index and the characteristic length of material.

Fig.12 shows the relationship between the homogeneity index and the characteristic length of material d1 as predicted by Eqs.(13)and (14). With the increase of the homogeneity index, the characteristic length of material decreases monotonically till it roughly reaches a constant value (called the ‘stable state'), and then it decreases again. There are three such stable states in Fig.12 and the last one corresponds to a homogenous material. Although empirical in nature, Eq. (13) well describes the data as shown in the figure.6. ConclusionsIn this paper, RFPA was used to investigate the discrepancy between the rock flexural and direct tensile strengths. The Weibull distribution was applied in the numerical simulation and the parameter m in the Weibull distribution was varied to obtain different rock materials with different levels of inhomogeneity.The direct tension and SCB tests were simulated by RFPA,and the direct tensile strength and the flexural tensile strength were obtained.The numerical results showed that when the material is inhomogeneous or heterogeneous, the experimentally observed discrepancy between the flexural and the direct tensile strength was reproduced. Therefore, this discrepancy was induced by the heterogeneity of the rock material. This discrepancy disappeared with the homogeneity index m reaching the critical value of 33,when the material was essentially homogenous.

Furthermore,the non-local fracture criterion used to explain the difference between these two tensile strengths was studied using the data from the numerical experiments. The relationship between the homogeneous index and characteristic length was established and three stages were proposed for such relationship.The material characteristic length can be viewed as an indirect measure of the material inhomogeneity according to this one-toone relationship. This relationship thus provided a physical meaning for the homogenous index in the Weibull distribution for heterogeneous materials.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This work was supported by the Natural Science Foundation of China(Grant Nos.51704211 and 11602165).The authors would like to acknowledge Changyi Yu of Tianjin University for his technical assistance in using software. The authors gratefully acknowledge anonymous reviewers for their valuable suggestions and comments.