Phase field model for simulating the fracture behaviors of some disc-type specimens

Journal Pre-proofs Phase field model for simulating the fracture behaviors of some disc-type specimens Longfei Wang, Xiaoping Zhou PII: DOI: Reference:

S0013-7944(19)31392-X https://doi.org/10.1016/j.engfracmech.2020.106870 EFM 106870

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

9 November 2019 2 January 2020 4 January 2020

Please cite this article as: Wang, L., Zhou, X., Phase field model for simulating the fracture behaviors of some disc-type specimens, Engineering Fracture Mechanics (2020), doi: https://doi.org/10.1016/j.engfracmech. 2020.106870

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Phase field model for simulating the fracture behaviors of some disc-type specimens Longfei Wang1,2 , Xiaoping Zhou 1,2, (1School of Civil Engineering, Chongqing University, Chongqing 400045, P. R. China; 2State

Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400044, P.R. China;)

Abstract: A staggered scheme for phase field modeling of brittle fractures is implemented to simulate mixed mode fracturing of disc-type specimens. The crack propagation process and load-displacement curves of each type of disc specimen are investigated. The effects of specimen geometry and loading mode on the crack propagation process and load-displacement curves are studied. The predicted results for crack initiation angle and fracture toughness obtained by the phase field method are compared with the analytical solutions. The results show that the crack initiation angle and fracture toughness obtained by the phase field method agree well with those determined by the generalized maximum tangential stress (GMTS) criterion. This implies that the phase field method has high accuracy in simulating the problem of crack-weakened Brazilian discs. Keywords: Phase field method; Disc-type specimens; Mixed mode fracture; Crack initiation and propagation; Fracture toughness

Nomenclature:  ∂Ω

An arbitrary domain External boundary

Ro Ri

Γ

Internal discontinuity boundary

a

u u t

Preserved displacement Displacement boundary Preserved traction





P 2S

Corresponding author.

E-mail address: [email protected] (X.P. Zhou). Tel./Fax.+86-23-6512-3511 1

The outer radius of the specimen The inner radius of the specimen Crack length (for CCCD, it is half length) Crack inclined angle Concentrated force Support spacing

t

Traction boundary

d

Phase field



Total potential energy

e()

Elastic strain energy density

gc

Critical energy release rate

l0

Length scale parameter

 (d,d)

Crack surface density

u

Displacements vector

ε



Strain tensor Tensile strain tensor Compressive strain tensor Principal strain Principal strain direction Elastic density caused by tension Elastic density caused by compression Volume force Cauchy stress tensor

H

History variable

+ a na

e+ eb

Nu i Bu i Bd i

t KCD-I, KSCB-I, KHD-I KCD-II, KSCB-II, KHD-II YCD-I, YSCB-I, YHD-I YCD-II, YSCB-II, YHD-II TCD, TSCB, THD TCD, TSCB, THD

Specimens thickness

rr,  , r 0

Stress components in the polar system Crack initiation angle Minimum mesh size Maximum mesh size Mode mixities Crack initiation time Stress intensity factor International Society for Rock Mechanics Short rod Chevron bend Cracked chevron notched Brazilian disc

hmin hmax Me t* SIF ISRM SR CB CCNBD

Displacement shape function associated with the ith node Strain gradient matrix of the node i Cartesian derivative matrices of the node i

Mode II fracture toughness

Mode I dimensionless SIF

Mode II dimensionless SIF

T stress

Dimensionless of T-stress

SCB

Semi-circular bend

CCCD

Centrally cracked circular disc

HCCD

Hollow central cracked disc

Ru

Displacement residue vector

GMTS

Rd

Phase filed residue vector

DDM

Ku

Displacement tangent matrix

DDA 2

Mode I fracture toughness

Generalized maximum tangential stress Displacement discontinuity method Discontinuous deformation analysis

Kd D H (x)

 

LEFM PFC

Phase field tangent matrix Fourth-order elasticity tensor Heaviside function Fourth-order projection tensor Fourth-order tensor Linear elastic fracture mechanics Particle flow code

DEM SPH XFEM GPD PD DLSP

Discrete element method Smooth particle hydrodynamics Extended Finite Element Method General particle dynamics Peridynamics Diagonally loaded square plate

1 Introduction Cracks are one of the most common discontinuities in rock masses and which can significantly affect the mechanical properties of rocks [1-2]. Moreover, under an external load, cracks may initiate, propagate and coalesce, leading to the failure of rock masses and engineering safety accidents. Therefore, understanding crack initiation, propagation and coalescence in rock masses has become a hot topic in the field of rock fracture mechanics in recent years [3-11]. In the framework of linear elastic fracture mechanics (LEFM), the stress intensity factor (SIF) is used to describe the displacement and stress field distribution at the crack tips [12]. Fracture toughness, which is defined as the critical value of the SIF, is an important parameter for crack problems, and it characterizes the ability of materials to prevent crack growth. The International Society for Rock Mechanics (ISRM) proposed four configurations of specimens to test the mode I rock fracture toughness: (a) short rod (SR) specimens [13], (b) chevron bend (CB) specimens [13], (c) cracked chevron notched Brazilian disc (CCNBD) specimens [14] and (d) straight through cracked semicircular bend specimens (SCB) [15]. However, in practical rock engineering, cracks usually present a certain angle between the crack centerline and the external loading direction. Therefore, in most cases, the crack is subjected to various mixed mode (combinations of mode I and mode II) loading. Therefore, rock failure under mixed mode loading should be studied. For mixed mode fracture testing, although the ISRM has not yet provided standard specimens and related methods, disc-type specimens [16-22] and bending beam type specimens [23-28] have been widely used in the study of mixed mode fracture toughness. Compared with bending beam specimens, disc specimens are more popular because of their simple 3

configuration, convenient preparation process and lower requirements for the loading equipment. Disc-type specimens, which are used for mixed mode fracture testing, usually include centrally cracked circular disc (CCCD) specimens [29,30], hollow central cracked disc (HCCD) [31-33] specimens and straight through cracked semicircular bend (SCB) [34-36] specimens. In addition, to better understand the initiation and propagation of mixed mode cracks in brittle and quasi-brittle materials, researchers have proposed many fracture criteria. These fracture criteria can be divided into three categories: stress-based criteria [37-39], strain-based criteria [40,41] and energy-based criteria [42-44]. Because of the simple form and clear physical meaning of each parameter, the generalized maximum tangential stress (GMTS) criterion is widely used in the prediction of mixed mode fractures dominated by tension and shear stresses [45-48]. On the other hand, due to the advantages of cost-savings, clear visualization of the crack growth process and strong repeatability, an increasing number of numerical methods are used to simulate crack growth. Vásárhelyi and Bobet [49] adopted the displacement discontinuity method (DDM) to simulate crack initiation, propagation and coalescence in rocks under uniaxial compression. Ning et al. [50] used an advanced discretization approach to simulate rock failure problems within the framework of the discontinuous deformation analysis (DDA). Bai et al. [51] investigated the fracture mechanism of rock discs containing hole(s) and its influence on tensile strength by using discrete element method (DEM). Based on a new computational approach combining smooth particle hydrodynamics (SPH) and a constitutive model, Wang et al. [52] simulated three mode I rock fracture tests including three-point bending test, Brazilian-disc test and semicircular bending test. Zhou and Yang [53] employed the extended finite element method (XFEM) to model crack initiation, propagation and coalescence in rocks. Zhou et al. [54] presented a numerical study of the cracking behavior of rock-like materials containing multiple pre-existing flaws under compression using general particle dynamics(GPD). Based on peridynamics (PD), Zhou and coworkers [55-56] conducted works on the cracking behaviors of rocks. Yang et al. [57] used a particle flow code (PFC) to investigate the 4

failure mechanical behavior of pre-holed granite specimens after elevated temperature treatment. However, the modeling of crack nucleation, merging, branching and complex crack paths still remains formidable task for discrete approaches and meshless methods. Therefore, a smeared or diffuse approach, known as the phase field method, has been developed as an alternative way to handle complex crack problems. Inspired by the works of Munford and Shah [58] and Ambrosio and Fortorelli [59], Francfort and Marigo [60] approached the fracture problem by the minimizing a potential energy based on Griffith’s theory of brittle fracture. Bourdin et al. [61] defined a scalar variable to regularize a sharp crack topology into a diffuse entity. Amor et al. [62] distinguished the fracture behavior between tension and compression by decomposing the elastic energy density into volumetric and deviatoric portions. Similarly, Miehe [63] proposed the method of spectral decomposition of strain to avoid cracks in the compressed part of the specimen. Moreover, a robust staggered scheme was employed by Miehe et al. [64], in which a local energy history field is introduced as a state variable to capture the irreversible nature of crack growth. After this work, due to simple implementation, no additional fracture criteria and easy handling of multiple cracks and complex crack trajectories, the phase field method have been widely applied to brittle fracture problems [65,66], cohesive fracture problems [67], dynamic fracture problems [68-70], multiphase materials fracture problems [71], hydraulic fracturing problems [72], and interface crack fracture problems [73-75]. Moreover, the phase field method (PFM) is different from the gradient damage model (GED): the PFM behaves like a cohesive zone model (a subclass of fracture mechanics) when the length scale is sufficiently small and the response is insensitive to this length scale, whereas the GED is a nonlocal damage model (a subclass of continuum damage mechanics) in which the response depends heavily on the length scale [76]. The PFM has the following advantages in simulating crack propagation: (1) it can model complex crack behavior without any ad hoc criterion; (2) it allows naturally incorporating multi-field physics owing to the variational principle; (3) it can intrinsically deal with merging and branching of 5

multiple cracks with no additional technique; (4) it can be generalized effortlessly to three dimensions and its computational implementation is straightforward in any dimension. In this paper, the mixed mode fracture behaviors of disc-type specimens are investigated using the phase field method. The innovation of this paper is described as follows: (1) the phase field method is adopted to study the mixed mode fracture behavior of CCCD, SCB and HCCD specimens; (2) the effects of the prefabricated crack inclination angle, specimen geometry and loading mode on the peak load, crack initiation angle and propagation path are analyzed; and (3) the numerical results obtained by the phase field method are compared with the theoretical results. This paper is organized as follows: A brief introduction to the phase field method is given in Section 2. Then, the validity and accuracy of the phase field method are verified to simulate tension-shear fractures in Section 3. In Section 4, the relevant theory of disc-type specimens and the mixed fracture criterion are introduced. The numerical results are analyzed and compared with the theoretical results in Section 5. Finally, conclusions are drawn in Section 6. 2 Introduction to the phase field method 2.1 Phase field approximation for total potential energy According to previous works [61,63], when the energy generated by external forces is neglected, the total potential energy for a cracked body  can be expressed as the sum of the elastic energy and crack surface energy:     e ( )d  +  g c d  



(1)

where    n , n  1, 2,3, represents an elastic n-dimensional body,    is the crack set,  e ( ) and gc represent the elastic energy density and critical energy release rate, respectively. In addition, the strain tensor     u  can be written as a function related to the displacement field u 1 2

  [u  (u)T ] 6

(2)

For an undamaged solid, the elastic energy density is defined as 1 2

2

 o ( )   tr ( )   tr ( ) 2

(3)

where  and  are the Lamé constants. According to the works by Miehe et al. [63,64], to account for only stress degradation in tension, the strain tensor can be decomposed into tensile and compressive parts

   

(4)

n      a na  na    a 1   n     a na  na   a 1 

(5)

In Equation (4),

where + and  - are the tensile and compressive strain tensors, respectively;

a is the

principal strain and na is its normal direction; and the symbol 〈·〉 denotes the Macaulay bracket, which can be defined as follows:

    

 

 (   ) / 2

(6)

 (   ) / 2

Therefore, the positive and negative elastic energy densities are obtained as follows:

      tr       e 2         tr     e 2

2 

2 

  tr   2    tr  

2 



(7)

The phase field method is an energy-based method that belongs to the category of smeared approaches. In the phase field model, the total energy functional (sum of bulk energy and surface energy) is minimized to obtain the crack nucleation and propagation. Moreover, a sharp crack is handled by regularizing over the entire domain using a scalar parameter called a phase field variable. The phase field variable is a function of spatial coordinates and takes the following form: x

d  x   e l0 7

(8)

where l0 is the length parameter used to represent the crack width, as shown in Fig.1. The crack approximation converges to a sharp crack as l0 approaches zero. That is, the larger the value of l0 is, the wider the crack. The smaller the value of l0 is, the narrower the crack.

Fig.1 Phase field parameter as a function of spatial coordinates

To describe the degradation of the stored elastic energy, a monotonically decreasing function

g  d   1  d   k is introduced. Therefore, the elastic strain 2

energy density function can be written as





 e   , d   1  d   k  e     e    2

(9)

To calculate the crack surface energy conveniently, Miehe et al. [63] introduced the crack surface density function

  d , d   Because





1 2 l0 2 d  d 2l0 2

(10)

d     d  , the fracture energy is approximated by 

 d 2 l0 2 g d   g   d   d c c   2l0 2  

(11)

Finally, substituting Equation (9) and Equation (11) into Equation (1) yields





 1 2 l0 2 2    1  d   k  e ( )  e ( ) d    g c  d  d  d      2  2l0 

(12)

2.2 Governing equations The variation of the functional  can be derived and its first variation should be zero, which leads to the following governing equations:

8

  b  0 in    2 1  d    g  d  l d   0 in   e c 0  l0        n = t on  t    u  n = u on  u   d  n = 0 on   where  is the external boundary of domain  ;

(13)

t and u are the traction 



boundary and displacement boundary, respectively; t and u are the preserved traction and displacement, respectively; and  is the Cauchy stress tensor and can be calculated by





      1  d   k   tr   1   2    tr    1  2  2





(14)

For two-dimensional plane problems, 1 is the second-order identity tensor. To ensure the irreversibility of crack growth under unloading conditions, Miehe et al. [64] introduce a strain-history field H (x, t), which can be defined by

H  x , t   max  e  x , t  t[0,T]

(15)

Replacing  e with H  x , t  in Equation (13), the damage evolution function can be rewritten as d  2 1  d  H  x , t   g c   l0 d   0 in   l0 

(16)

2.3 Numerical implementation of the phase field method Now, the numerical model is discretized in a numerical setting. The discretized displacement is given by m

u   N iu ui

(17)

i

where m is the total number of nodes in one element, ui is the displacement of node I; and N iu is the shape function associated with node i. The discretized phase field is given by 9

m

d   N id di

(18)

i 1

where di is the phase field of node i, and N id is also a shape function related to node i, but its form is different from that of N iu . For two-dimensional problems, the shape function N iu and N id take the following forms:  u  N1 0  N m 0   Ni     0 N1  0 N m   d  N i   N1 N 2  N m 

(19)

The strain field and phase field gradient can be discretized into m

   Biu ui

(20)

i

m

d   Bid di

(21)

i 1

where Biu is the strain gradient matrix of node i and Bid is the Cartesian derivative matrix. They can be respectively calculated by

  N m N1 N 2 0 0  0   x x x     N m N1 N 2 u 0 0 Bi    0  y y y    N1 N1 N 2 N 2 N m N m     x y x y x   y

 N1  x Bid    N1  y

N 2 x N 2 y

N m  x   N m   y 

(22)



(23)

The residue vectors Ru and R d at the element level take the following form Ru     ( Biu )T  d 

10

(24)

 T T   g R d     c d  2 1  d  H   N id   g c l0  Bid  d  d    l0 

(25)

The Newton-Raphson algorithm is used to update the tangent matrix and the residue vector at each internal iteration:  K mu um 1   Rmu  d d  K m d m 1  Rm

(26)

where d m 1 and um 1 contain the new phase field values and displacement values at each integration point, respectively. The tangent matrix on the element level is thus calculated by K u    ( Biu )T DBiu d 

  T T g  K d     2 H  c   N id  N id  g c l0  Bid  Bid  d  l0   

(27)

(28)

where D is the elasticity tensor of fourth order, and can be given by D

  2  1  d   k    H  (tr( )) + 2       H  ( tr( )) + 2      

(29)

where H (x) is the Heaviside function with H (x) = 1 for x  0 and H (x) = 0 for x  0, and  is a fourth-order tensor such that if i = k and j = l, the result is unity, otherwise the result is zero.   are the fourth-order projection tensors and can be expressed as 3

3

 Pijkl   H  ( a ) ab nai naj nbk nbl a 1 b 1

1  a     a   nai nbj  nak nbl  nbk nal    a  b a 1 b  a 2 3

3

3

3

(30)

 Pijkl   H  ( a ) ab nai naj nbk nbl a 1 b 1

1  a     a   nai nbj  nak nbl  nbk nal  a  b a 1 b  a 2 3

3

 

(31)

3 A benchmark example A benchmark example is taken to verify the validity and accuracy of the phase 11

field method for simulating crack propagation under tension-shear loads.

Fig.2 A diagonally loaded square plate (DLSP) specimen

As shown in Fig.2, a benchmark example named the diagonally loaded square plate (DLSP) specimen proposed by Ayatollahi and Aliha [77] is applied to investigate the mixed mode fracture. According to the work by [77], the dimensions of the DLSP specimens are listed as follows: 2w = 150 mm and 2a = 45 mm. Thus, the ratio of a / w is equal to 0.3 in the test samples. Two holes with a radius of 4 mm are drilled near two opposite corners of the DLSP specimen. The center of each hole is 25 mm from the specimen corner. The following crack angles are considered for experiments:  = {0o (pure mode I), 15o, 30o, 45o and 62.5o (pure mode II)}.

Fig.3 Comparisons of the specimen failure modes between experimental results and phase field simulation results

Fig.3 shows the comparisons of the specimen failure modes between the experimental results and the phase field simulation results. It can be seen from Fig.3 that the crack propagation path obtained by the phase field method is in good agreement with the experimental observation, implying that the phase field method can accurately obtain the crack propagation path under tension-shear loading. Fig.4 illustrates that the peak load trend obtained by the phase field method agrees well with the experimental results, which implies that the phase field method not only can capture the path of crack propagation, but also can accurately obtain the trend of the peak load when the specimen is subjected to tensile-shear loading.

Fig.4 Comparisons of the peak load between the experimental results and phase field simulation results

12

4 Disc-type specimens and fracture criterion 4.1 Disc-type specimens Typically, disc specimens used for rock mixed mode fracture testing include CCCD, SCB and HCCD, as shown in Fig.5. Fig.5 Three disc-type specimens As seen in Fig.5, the three types of disc-type specimens can be divided into two categories based on the loading mode, namely, CCCD specimens and HCCD specimens belong to the diametric compression type, while SCB specimens belong to the three-point bending type. For CCCD specimens, the radius of the specimen is Ro, and there is an inclined crack with a length of 2a in the center. The crack inclination angle  indicates the angle between the crack centerline and the line of diametric concentrated force P. In SCB specimens, there is a crack with a length of a , whose angle  is with respect to the loading force P, and the support spacing of the specimen is 2S. In HCCD specimens, Ro and Ri stand for the outer and inner radii, respectively; P stands for the diametric concentrated forces, which are applied on the top and bottom of the disc; two diametric cracks with lengths of a are located at the edge of the central hole; and the crack angle with respect to the loading force is still . If  = 0o, the specimens are under pure mode I loading. With increasing , the specimens gradually change from pure mode I loading to mixed mode I/II loading. When  increases to a certain value, the specimens are in a pure mode II loading state, and this value depends on factors such as crack length and geometric configuration of the specimens. According to previous works [15,17,32], the stress intensity factors (KI and KII) and the T-stress in these specimens can be calculated by:

13

 P 2 Ro YCD-I  K CD-I  Rot   P 2 Ro  YCD-II for CCCDspecimen  K CD-II  Rot   4P * TCD  TCD  Rot 

(32)

 P a YSCB-I  KSCB-I  2 Rot   P a YSCB-II for SCBspecimen  KSCB-II  2 Rot   P * TSCB  TSCB  2 Rot 

(33)

 P  aYHD-I  K HD-I   R R t ) ( i o   P  aYHD-II for HCCDspecimen  K HD-II   R R t ) ( i o   P * THD  THD   t R R ) ( i o 

(34)

where YI, YII and T* are the dimensionless forms of KI, KII and T*, respectively, t is the thickness of the specimen, and t = 1 is assumed for a two-dimensional problem. 4.2 Mixed mode fracture criterion According to work by Williams [78], the elastic stress field around the crack tip under mixed mode I/II loading can be written as follows:  1      3 cos  K I 1  sin 2   K II  sin   2 tan     rr  2  2 2  2 r 2    T cos 2   O  r1/2    1   3  cos  K I cos 2  K II sin    T sin 2      2 2 2 2 r     O  r1/2    1    cos   K I sin  K II  3cos   1   T sin  cos   r  2 2 2 r    1/2   O r  

(35)

where r and  are the local crack tip coordinates; rr,  and r are the stress 14

components in the polar system, whose origin is located at the crack tip; KI and KII are the mode I and mode II stress intensity factors (SIFs); and T is a nonsingular and constant stress term, that is independent of the distance with respect to the crack tip, while it is dependent on the specimen geometry and the loading condition. Fig.6 The stress and displacement fields around the crack tip Fig.7

The definition of the crack initiation angle 0

Fig.7 shows the definition of the crack initiation angle; that is, 0 is the angle between the newly generated crack and the pre-existing crack. Erdogan and Sih [37] proposed the traditional maximum tangential stress (MTS) criterion to predict mixed mode crack propagation. According to their work, the crack initiation angle 0 can be obtained as follows: cos

0 2

 K I sin 0  K II (3cos 0  1)  0

(36)

However, some researchers [79-81] showed that there is a large discrepancy between the results predicted using the MTS criterion and the experimental observations because the MTS criterion considers only singular terms (SIFs) in Williams’ expansion and ignores higher order terms. To solve this problem, some improved MTS criteria have been proposed [82]. Among these improved MTS criteria, the generalized maximum tangential stress criterion (GMTS) is widely used, in which the second term of the higher order term (T-stress) is taken into account in the fracture criteria and the higher order terms of the expansion are ignored. In the GMTS, the crack initiation angle 0 satisfies the following relation   

 0

0   K I sin  0  K II  3cos  0  1  

16T 3

2 rc cos  0 sin

0 2

(37)

0

where rc is a critical distance and is considered a material constant. According to work by Taylor et al. [83], rc can be expressed as the following 15

function of the mode I fracture toughness KIc and tensile strength t

1  K Ic  rc    2   t 

2

(38)

Substituting Equations (32) - (34) into Equation (37) yields

 64 CCCD : YI sin  0  YII  3cos  0  1  T *  3   SCB and HCCD :Y sin   Y 3cos   1  16 T *  0 0 I II   3

rc  cos  0 sin 0  0 2 R 2rc  cos  0 sin 0  0 2 a

(39)

Based on the GMTS criterion, the onset of mixed mode I/II fracturing takes place when

K Ic  cos

0 

 3  K I cos 2 0  K II sin  0   2 rc T sin 2  0  2 2 2 

(40)

By substituting Equations (32) - (34) into Equation (40), the following dimensionless expressions can be obtained: K      3Y  I  cos 0  cos 2 0  II sin  0   2 2 2 YI  K Ic   CCCD :   K II   0  YI   3  cos  cos 2 0  sin  0    2  YII 2 2   K Ic 

 rc 4T * sin 2  0  R YI 

1

 rc 4T sin 2  0  R YII 

1

*

1 K   0  2  0 3 YII   2rc T * 2 I    cos  cos sin  0  sin  0   2 2 2 YI a YI  K Ic    SCB and HCCD :  1   K II   0  YI  2rc T * 3 2 0 2  cos  cos  sin  0   sin  0   2  YII 2 2 a YII   K Ic  

(41)

(42)

5 Results and discussion 5.1 The numerical model and experimental comparison In this section, based on the numerical results, the effects of the crack inclination angle on crack propagation and fracture toughness of three disc-type specimens are studied. Furthermore, the numerical results are compared with the theoretical solutions. Finally, the influences of the geometry and loading mode of the specimens are discussed. 16

In this paper, isotropic fine-grained sandstone specimens are used. According to previous works [84], the elastic modulus is E = 20 GPa, Poisson’s ratio is v = 0.25, its critical energy release rate is gc = 15.68 N/m, and the internal length parameter is l0 = 0.5 mm. Moreover, to investigate the influence of mesh size on the simulation results, six different mesh densities are considered as follows: h = l0 /5, h = l0 /4, h = l0 /3, h = l0 /2, h = l0 and h =2l0. The numbers of elements and nodes in the models with different mesh densities are shown in Tab.1. Fig.8 represents the influence of mesh density on the numerical results. Fig.8(a) shows that the mesh density has a certain influence on the load displacement curve. Fig.8(b) illustrates that the peak load of the specimen increases continuously with increasing mesh size. However, when the mesh density is no more than l0 /2, the influence of mesh density on the numerical results can be ignored. Considering the speed and accuracy of the numerical calculation, l0 /2 is chosen as the refinement density of the mesh in the region of possible crack propagation. That is, the minimum mesh size is hmin = 0.25 mm, and the mesh size of the other regions is hmax/hmin=8. According to the works by Zhang et al. [85] and Yin and Zhang [86], the displacement increment is u=110-5 mm. . Tab.1 Numbers of elements and nodes in the models with different mesh densities

Fig.8 Influence of mesh density on numerical results

The geometric parameters of the specimen are listed as follows: the specimen radius is Ro = 25 mm, and the length of prefabricated crack a is equal to 7.5 mm (a/R = 0.3). For SCB specimens, the support spacing is S = 10 mm (S/Ro = 0.4). For HCCD specimens, the inner diameter of the specimen is Ri = 5 mm (Ri /Ro = 0.2). Many specimens can be used for mixed mode fracture tests, such as Brazilian discs, threepoint bending beams and four-point bending beams. For specimens with different configurations, when the transition from pure mode I loading to pure II mode loading is realized, the required range for the crack inclination angle is obviously different. 17

For example, for a centrally cracked circular disc(CCCD) specimen, to achieve mode II loading, the crack inclination angle is generally located in the range from 200 to 300, while for the straight through cracked semicircular bend (SCB) specimen, the pure mode II loading angle is in the range of 350 ~ 600 [2,87]. Therefore, when the influence of specimen configuration on the numerical simulation results is studied, it is not obviously appropriate to adopt the same inclination angle for the prefabricated crack, so we need to introduce a parameter that can describe the degree of mixed mode loading of the specimens. Me is the mixity parameter which shows the ratio of the mode I and mode II stress intensity factors. It can be determined by Me 

K  tan 1  I    K II  2

(44)

Because the calculation formula of the mode I stress intensity factor KI is similar to that of the mode II stress intensity factor KII, the difference is related only to the dimensionless stress intensity factors YI and YII. Therefore, Equation (44) can be expressed as Me 

Y  tan 1  I    YII  2

(45)

In this paper, the mode I and mode II stress intensity factors (SIFs) of three kinds of specimens are obtained by using the interaction integral method. Then, the dimensionless SIFs can be found by dimensionless treatment of the SIFs. Finally, the mixity parameter Me can be obtained by substituting the mode I and mode II dimensionless SIFs into Equation (45). Therefore, three dimensionless crack tip parameters are obtained by the interaction integral method, as shown in Fig.9. Fig.9 Dimensionless SIFs and T-stress of three disc-type specimens All the specimens were simulated at similar mode mixities of Me = 1 (pure mode I), 0.8, 0.6, 0.4, 0.2, and 0 (pure mode II ) . The crack angles corresponding to each mode mixity for each specimen are shown in Tab.2. 18

Tab.2 The values of  for the required mode mixities in three disc-type specimens

The numerical and experimentally observed [21,87-88] crack propagation paths of the three kinds of disc specimens under mixed mode loading are plotted in Fig.10. As shown in Fig.10, the crack paths are reasonably affected by the test configurations, and the simulation results reasonably match the realistic crack pattern.

Fig.10 Crack path comparison between numerical results and experimental observation

5.2 Centrally cracked circular disc (CCCD) specimens

Fig.11 Crack propagation of CCCD specimens with different mode mixities Me

Fig.11 shows the failure modes of CCCD specimens with different mode mixities Me. As shown in Fig.11, when Me is equal to 0-0.8, all cracks initiate from the prefabricated crack tips, and then propagate to the loading end of the specimen along curvilinear paths. When Me = 1, the cracks propagate from the prefabricated crack tips to the loading end along straight paths. In addition, Fig.11 shows that the crack initiation angle is quite different when the prefabricated crack inclination angle changes. To investigate the difference quantitatively, the crack initiation angles of different mode mixities Me are measured. In addition, it is well known that both the MTS criterion and GMTS criterion can theoretically predict the crack initiation angle. Therefore, the results obtained by the phase field method are compared with those obtained by both the MTS criterion and GMTS criterion, as shown in Fig.12.

Fig.12 Comparison between theoretical and numerical crack initiation angles for CCCD specimens 19

Fig.12 shows that for CCCD specimens, the crack initiation angle obtained by MTS (without considering T-stress) is larger than that obtained by GMTS (considering T-stress). Many experiments [81-82,84] have shown that the crack initiation angle predicted by GMTS is more accurate. In addition, the results obtained by the phase field method are in good agreement with those predicted by the GMTS criterion. This result implies that the phase field method can accurately capture the crack initiation angle in CCCD specimens.

Fig.13 Load-displacement curves for CCCD specimens

Fig.13 shows the load-displacement curves of CCCD specimens with different mode mixities Me. As shown in Fig.13, the peak loads and maximum displacements increase continuously with increasing Me. This relation implies that for CCCD specimens, when the prefabricated crack inclination angle is between the pure mode I loading angle and the pure mode II loading angle, the larger the prefabricated crack inclination angle is, the smaller the load-bearing capacity of the specimens. By substituting the peak load obtained by the phase field method and the dimensionless stress intensity factor obtained by the interaction integral method into Equation (32), the mode I and II fracture toughness of CCCD specimens with different Me can be determined. Fig.14 shows a comparison between the numerical mode I and II fracture toughness obtained by the phase field method and theoretical results obtained by the MTS and GMTS criteria in the tested CCCD specimens. As shown in Fig.14, the numerical results obtained by the phase field method are in good agreement with the theoretical predictions of the GMTS criterion.

Fig.14 Mixed mode fracture toughness obtained by the phase field method and theoretical criteria in CCCD specimens

5.3 Straight through cracked semicircular bend (SCB) specimens 20

Fig.15 Failure mode of SCB specimens with different mode mixities Me

Fig.15 shows the failure modes of SCB specimens with different mode mixities Me. It can be seen from Fig.15 that the failure modes of SCB specimens with different prefabricated crack inclination angles are similar to those of CCCD specimens. All new cracks initiate from the prefabricated crack tip, and then propagate to the loading point of specimens along straight line paths (pure mode I) and curvilinear paths (mixed mode I/II or pure mode II). In addition, the numerical crack initiation angles of SCB specimens from the phase field method are determined and compared with the theoretical predictions using the MTS criterion and GMTS criterion, as shown in Fig.16. Fig.16 clearly illustrates that the crack initiation angles obtained by the phase field method are consistent with those predicted by the GMTS criterion, and they are all larger than those predicted by the MTS criterion, which is contrary to the results for CCCD specimens; similar results have been obtained in the literature [87]. This contrast is related to the difference in the T-stress of the specimens with the two configurations. It can be seen from Fig.9 (c) that when the inclination angle of the prefabricated crack is between the pure mode I loading angle and the pure mode II loading angle, the T-stress of the CCCD specimens is negative, while the T-stress of the SCB specimens is positive in most cases. According to previous work [39], the positive T-stress causes an increase in the crack initiation angle, while the negative T-stress leads to a decrease in the crack initiation angle. Therefore, the initiation angle of the SCB specimen predicted by the GMTS criterion is larger than that predicted by the MTS criterion, while the initiation angle of the CCCD specimen predicted by the GMTS criterion is smaller than that predicted by the MTS criterion.

Fig.16 Comparison between theoretical and numerical crack initiation angles for SCB specimens

Fig.17 Load-displacement curves for SCB specimens 21

Fig.17 shows the load-displacement curves for the SCB specimens. After the SCB specimens are loaded, the load increases continually with increasing displacement. When the peak load is reached, the load-carrying capacity of the specimens decreases rapidly with an increase in displacement. Comparing the peak loads of samples with different Me, it can be found that the peak load of specimens decreases with increasing Me, which is contrary to the trend of change for CCCD specimens. Before the peak loads, the load-displacement curves in specimens with different Me are approximately the same, implying that the prefabricated crack inclination angle has little effect on the load and displacement of the specimens. It can be seen from Fig.13 that when the crack propagation is completed, the load of the CCCD specimen does not decrease to zero, which indicates that the CCCD specimen still has a certain bearing capacity after failure. It is not difficult to observe from Fig.17 that the bearing capacity of the SCB specimen is basically close to 0 after failure. The reason for this difference is related to the loading modes of the specimens. A similar phenomenon was observed in the experiment reported in [87].

Fig.18 Mixed mode fracture toughness of SCB specimens obtained by the phase field method and theoretical criteria

Fig.18 shows a comparison of the numerical and theoretical results for fracture toughness of SCB specimens. It can be seen from Fig.18 that the fracture toughness obtained by the phase field method is in good agreement with that predicted by the GMTS criterion, which is smaller than that predicted by the MTS criterion. This difference is related to the positive T-stress in SCB specimens. 5.4 Hollow central cracked disc (HCCD) specimens Fig.19 Failure mode of HCCD specimens with different mode mixities Me

Fig.19 shows the failure mode of HCCD specimens with different prefabricated crack inclination angles. It is found from Fig.19 that the crack propagation paths of 22

HCCD specimens are similar to those of CCCD and SCB specimens, i.e. there are straight line paths under pure I mode loading and curvilinear paths under mixed mode I/II or pure mode II loading. Similarly, the crack initiation angles obtained by the phase field method are compared with those obtained by the fracture criterion, as shown in Fig.20. It can be seen from Fig.20 that when Me is equal to 0 or 1, that is, when the specimen is in the pure mode I or II loading state, the crack initiation angle predicted by the GMTS criterion is the same as that predicted by the MTS criterion. When the specimen is in the mixed mode I/II loading state, the crack initiation angle predicted by the GMTS is smaller than that obtained by MTS. In addition, the crack initiation angle of HCCD specimens obtained by the phase field method is in good agreement with that predicted by GMTS, which validates the accuracy of the phase field method.

Fig.20 Comparison between theoretical and numerical crack initiation angles for HCCD specimens

Fig.21 Load-displacement curve for HCCD specimens

Fig.21 illustrates that as the displacement increases, the load-bearing capacity of HCCD specimens increases first, then decreases, and finally increases. In addition, the second local maximum load in HCCD specimens is lower than the first local maximum load, which is quite different from the trends of CCCD specimens and SCB specimens. In the CCCD specimens and SCB specimens, the load-bearing capacity first increases, and then decreases. The reason for this phenomenon is related to the hollow in the center of the HCCD specimens.

Fig.22 Mixed mode fracture toughness obtained by the phase field method and fracture criteria for HCCD specimens

Fig.22 shows a comparison of the fracture toughness obtained by the numerical 23

method and fracture criteria. It is not difficult to determine that the numerical results obtained by the phase field method are in good agreement with those obtained by the GMTS criterion. Notably, the fracture toughness of HCCD specimens predicted by the GMTS criterion is close to that predicted by the MTS criterion under pure I and II loading conditions. This result occurs because the crack initiation angle 0 is equal to 0o when the crack is in the pure mode I loading state, and the last terms in Equations (41) and (42) become to 0. Therefore, the GMTS criterion prediction formula is the same as that of the MTS criterion. When HCCD specimens are in the pure mode II loading state, it can be seen from Fig.9(c) that the T-stress is basically close to 0, and the last terms in Equations (41) and (42) are also changed to 0. Therefore, the GMTS criterion prediction results are the same as the MTS criterion prediction results. 5.5 Geometry and loading mode effects Geometry and loading mode effects are two main factors influencing the accuracy of fracture toughness testing. Taking disc-type specimens as an example, CCCD and HCCD specimens are subjected to diametric compression loads, and have the same loading mode, therefore we can compare the numerical results of these two types of specimens to study the geometric effects. CCCD specimens and SCB specimens are straight cracked disc specimens, the loading mode of CCCD specimens is diametric compression, while SCB specimens are subjected to three-point bending. Therefore, their numerical results are compared to study the effect of the loading mode. First, it can be seen from Figs.13 and 21 that the load-displacement curves of the CCCD and HCCD specimens are not zero after failure, which indicates that the CCCD and HCCD specimens still have a certain load-bearing capacity after failure. This result occurs because the CCCD specimens and HCCD specimens are split into left and right parts after failure. Due to the diametric loading, the left and right parts are not completely separated. Therefore, the CCCD specimens and HCCD specimens still have load-bearing capacity after failure. As shown in Fig.17, the load-bearing capacity of SCB specimens drops almost to zero after failure. This result occurs because the SCB specimens are subjected to three-point bending loads, and the left 24

and right parts are completely separated after failure. It can be found from Fig.21 that a "local minimum" occurs on the load-displacement curve of HCCD specimens. Some researchers [89-92] used the relation between the local minimum load and the maximum dimensionless stress intensity factor to test the fracture toughness of HCCD specimens, thus avoiding the determination of crack length in the tests. Generally, both the geometry effects and loading modes have effects on the load-displacement curves of specimens. Geometric factors mainly affect the trends of the curves, while loading modes mainly affect the amplitudes of the curves.

Fig.23 Comparison of crack initiation angles of three disc specimens

Fig.23 shows the numerical results for the crack initiation angles of the three disc specimens. It can be seen from Fig.23 that under the same mode mixities Me, the crack initiation angle of SCB specimens is the largest, followed by HCCD specimens and CCCD specimens in turn. These results indicate that the geometry and loading mode of the specimen have great influences on the crack initiation angle. Fig.24 The trend of mode I and mode II fracture toughness of three kinds of disc specimens The mode I and mode II fracture toughness of three kinds of disc specimens are plotted in Fig.24. It can be concluded from Fig.24 that as Me increases, the mode I fracture toughness of all disc specimens increase, while the mode II fracture toughness decreases. When Me increases from 0 to 1, the increase in the mode I fracture toughness of CCCD specimens is the smallest, and the increase in SCB specimens is the largest; the decrease in the mode II fracture toughness of SCB specimens is the smallest, while the decrease in mode II fracture toughness of CCCD specimens is the largest. In addition, the geometry and loading mode affect the crack initiation angle and fracture toughness predicted by the fracture criterion for positive and negative T-stress. 25

For example, when Me varies from 0 to 1, the T-stress of the CCCD and HCCD samples are negative, while the T-stress of the SCB samples is mostly positive. Therefore, the crack initiation angles and fracture toughness for the CCCD and HCCD specimens predicted by the MTS criterion are larger than those predicted by the GMTS criterion, but the crack initiation angle and fracture toughness for the SCB specimens predicted by the MTS criterion are smaller than those predicted by the GMTS criterion.

Fig.25 Comparison of crack initiation time for the three disc specimens

To study the relation between the crack initiation time and the peak load time, the ratio t* of the specimen displacement at the crack initiation to the specimen displacement at the peak load is defined as the dimensionless crack initiation time. When t* > 1, the crack initiates after the peak load; when t* = 1, the crack initiates at the peak load; when t* < 1, the crack initiates before the peak load. The dimensionless crack initiation time t* of different specimens is plotted in Fig.25 (a) and (b). It can be seen from Fig.25 (a) that the dimensionless crack initiation time of all disc specimens is less than 1, which indicates that all specimens initiate before the peak load. Fig.25 (b) shows that the inclination angles of prefabricated cracks have little effect on the initiation time of SCB specimens, while it has great influence on the initiation time of CCCD and HCCD specimens. 6 Conclusions The phase field method is used to simulate the mixed mode fracturing of disc-type specimens. The crack propagation process and load-displacement curves of disc specimens are obtained. The crack initiation angle and fracture toughness predicted by the phase field method are compared with those obtained by mixed mode fracture criteria. Finally, the geometry and loading mode effects are discussed. The following conclusions are drawn: (1) The phase field method can efficiently capture the initiation and propagation process of tensile-shear mixed mode cracks without external criteria, and can 26

accurately determine the fracture toughness of Brazilian disc specimens in combination with the interaction integral method. (2) When the specimen is in the mode I loading state, the cracks initiate at the prefabricated crack tips, and propagate to the loading end of the specimen along straight line paths. When the specimen is in mixed mode or pure mode II loading, the cracks propagate along curvilinear paths. (3) The load-bearing capacity of the SCB specimens basically decreases to zero after failure. The CCCD and HCCD specimens still maintain their load-carrying capacity after failure. (4) As the prefabricated crack inclination angle undergoes the transition from mode I (Me = 1) to pure mode II (Me = 0), the mode I fracture toughness of all disc specimens decreases continuously, while the mode II fracture toughness increases continuously. (5) The crack initiation angle and fracture toughness of the CCCD and HCCD specimens predicted by the MTS criterion are larger than those predicted by the GMTS criterion, while the crack initiation angle and fracture toughness of the SCB specimens predicted by the MTS criterion are smaller than those predicted by the GMTS criterion. The numerical results obtained by the phase field method are in good agreement with those predicted by the GMTS criterion. Data Availability Data supporting this research article are available from the corresponding author via e-mail. Declaration of Competing Interest The authors declare that there is no conflict of interest. Acknowledgements The work is supported by the National Natural Science Foundation of China (Nos. 51809198, 51839009 and 51679017), Fundamental Research Funds for the Central Universities (No. 2042018kf0008). References [1] Aliha MRM, Mahdavi E, Ayatollahi MR. The influence of specimen type on tensile fracture 27

[2]

[3]

[4]

[5] [6]

[7]

[8] [9]

[10]

[11]

[12] [13] [14]

[15]

[16] [17]

toughness of rock materials. Pure Appl Geophys 2016;174:1237-53. Xie YS, Cao P, Jin J, Wang M. Mixed mode fracture analysis of semi-circular bend (SCB) specimen: A numerical study based on extended finite element method. Comput Geotech 2017;82:157-72. Wang YT, Zhou XP, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics. Eng Fract Mech 2016;163:248-73. Haeri H, Sarfarazi V, Yazdani M, Shemirani AB, Hedayat A. Experimental and numerical investigation of the center-cracked horseshoe disk method for determining the mode I fracture toughness of rock-like material. Rock Mech Rock Eng 2018;51:173-85. Cheng Y, Wong LNY, Zou C. Experimental study on the formation of faults from en-echelon fractures in Carrara Marble. Eng Geol 2015;195:312-26. Yang SQ, Tian WL, Jing HW, Huang YH, Yang XX, Meng B. Deformation and damage failure behavior of mudstone specimens under single-stage and multi-stage triaxial compression. Rock Mech Rock Eng 2019;52:673-89. Liu K, Zhang QB, Wu G, Li JC, Zhao J. Dynamic mechanical and fracture behaviour of sandstone under multiaxial loads using a triaxial Hopkinson bar. Rock Mech Rock Eng 2019; 52:2175-95. Wasantha PLP, Ranjith PG, Zhao J, Shao SS, Permata G. Strain rate effect on the mechanical behaviour of sandstones with different grain sizes. Rock Mech Rock Eng 2015; 48:1883-95. Zhu Q, Li D, Han Z, Li X, Zhou Z. Mechanical properties and fracture evolution of sandstone specimens containing different inclusions under uniaxial compression. Int J Rock Mech Min Sci 2019;115:33-47. Yang SQ, Huang YH, Tian WL, Yin PF, Jing HW. Effect of High Temperature on deformation failure behavior of Granite specimen containing a single fissure under uniaxial compression. Rock Mechanics and Rock Engineering 2019; 52: 2087-2107. Haeri H, Shahriar K, Marji MF, Moarefvand P. Modeling the propagation mechanism of two random micro cracks in rock Samples under uniform tensile loading. 13th International Conference on Fracture 2013; ICF 2013 2:1472-1481. Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 1957;24:361-4. Ouchterlony F. ISRM Suggested methods for determining fracture toughness of rocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25:71-96. Fowell RJ. Suggested method for determining mode I fracture toughness using cracked chevron notched Brazilian disc (CCNBD) specimens. Int J Rock Mech Min Sci Geomech Abstr 1995;32:57-64. Kuruppu MD, Obara Y, Ayatollahi MR, Chong KP, Funatsu T. ISRM-suggested method for determining the mode I static fracture toughness using semi-circular bend specimen. Rock Mech Rock Eng 2014;47:267-74. Jiang C , Zhao GF , Khalili N. On crack propagation in brittle material using the Distinct Lattice Spring Model. Int J Solids Struct 2017;118-119:41-57. Hou C, Wang ZY, Liang WG, Li JB, Wang ZH. Determination of fracture parameters in center cracked circular discs of concrete under diametral loading: A numerical analysis and experimental results. Theor Appl Fract Mech 2016;85:355-66. 28

[18] Mirsayar MM. On fracture analysis of dental restorative materials under combined tensile-shear loading. Theor Appl Fract Mech 2017;93:170-6. [19] Cai M. Fracture initiation and propagation in a brazilian disc with a plane interface: a numerical study. Rock Mech Rock Eng 2013;46:289-302. [20] Ayatollahi MR, Sistaninia M. Mode II fracture study of rocks using Brazilian disk specimens. Int J Rock Mech Min Sci 2011;48:819-26. [21] Chen CH, Chen CS, Wu JH. Fracture toughness analysis on cracked ring disks of anisotropic rock. Rock Mech. Rock Eng 2008;41:539-62. [22] Funatsu T, Kuruppu M, Matsui K. Effects of temperature and confining pressure on mixed-mode (I–II) and mode II fracture toughness of Kimachi sandstone. Int J Rock Mech Min Sci 2014;67:1-8. [23] Belli R, Wendler M, Petschelt A, Lohbauer U. Mixed-mode fracture toughness of texturized LS2 glass-ceramics using the three-point bending with eccentric notch test. Dent Mater 2017;33:1473-7. [24] Aliha MRM, Mousavi SS, Ghoreishi SMN. Fracture load prediction under mixed mode I + II using a stress based method for brittle materials tested with the asymmetric four-point bend specimen. Theor Appl Fract Mech 2019;103:102249 [online]. [25] Dilip KM, Girish YS. Experimental investigations on stable crack growth in three-point bend specimens under mode I and mixed mode (I and II) loading. Eng Fract Mech 1992;43:41-53. [26] Haeri H, Sarfarazi V, Hedayat A. Suggesting a new testing device for determination of tensile strength of concrete. Struct Eng Mech 2016;60:939-952. [27] Razavi SMJ, Aliha MRM, Berto F. Application of an average strain energy density criterion to obtain the mixed mode fracture load of granite rock tested with the cracked asymmetric four-point bend specimens. Theor Appl Fract Mech 2018;97:419-25. [28] Xeidakis GS, Samaras IS, Zacharopoulos DA, Papakaliatakis GE. Crack rowth in a mixed-mode loading on marble beams under three point bending. Int J Fract 1996;79:197-208. [29] Lim IL, Johnston IW, Choi SK. Assessment of mixed-mode fracture toughness testing methods for rock. Int J Rock Mech Min Sci Geomech Abstr 1994; 31:265-72. [30] Chang SH, Lee C, Jeon S. Measurement of rock fracture toughness under modes I and II and mixed-mode conditions by using disc-type specimens. Eng Geol 2002; 66:79-97. [31] Amrollahi H, Baghbanan A, Hashemolhosseini H. Measuring fracture toughness of crystalline marbles under modes I and II and mixed mode I–II loading conditions using CCNBD and HCCD specimens. Int J Rock Mech Min Sci 2011; 48:1123-34. [32] Zhou XP, Wang LF, Berto F, Zhou LS. Comprehensive study on the crack tip parameters of two types of disc specimens under combined confining pressure and diametric concentrated forces. Theor Appl Fract Mech 2017;103:102317 [online]. [33] Akbardoost J, Ghadirian HR, Sangsefidi M. Calculation of the crack tip parameters in the holed-cracked flattened Brazilian disk (HCFBD) specimens under wide range of mixed mode I/II loading: The mixed mode crack tip parameters for HCFBD specimens, Fatigue Fract Eng Mater Struct 2017; 40:1416-27. [34] Saha G, Biligiri KP. Fracture properties of asphalt mixtures using semi-circular bending test: A state-of-the-art review and future research. Constr Build Mater 2016;105:103-12. [35] Dai F, Xia K, Zheng H, Wang YX. Determination of dynamic rock Mode-I fracture 29

[36]

[37] [38] [39] [40] [41] [42] [43] [44]

[45] [46]

[47] [48] [49] [50]

[51]

[52]

[53]

parameters using cracked chevron notched semi-circular bend specimen. Eng Fract Mech 2011;78:2633-44. Du HB, Dai F, Xia KW, Xu NW, Xu Y. Numerical investigation on the dynamic progressive fracture mechanism of cracked chevron notched semi-circular bend specimens in split Hopkinson pressure bar tests. Eng Fract Mech 2017;184:202-17. Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J. Basic Eng. ASME 1963;85:519-25. Smith DJ, Ayatollahi MR, Pavier MJ. The role of T-stress in brittle fracture for linear elastic materials under mixed-mode loading. Fatigue Fract Eng Mater Struct 2001;24:137-50. Akbardoost J, Ayatollahi MR. Experimental analysis of mixed mode crack propagation in brittle rocks: the effect of non-singular terms. Eng Fract Mech 2014;129:77-89. Chang KJ. On the maximum strain criterion-a new approach to the angled crack problem. Eng Fract Mech 1981;14:107-24. Mirsayar MM, Razmi A, Aliha MRM, Berto F. EMTSN criterion for evaluating mixed mode I/II crack propagation in rock materials. Eng Fract Mech 2018;190:186-97. Lazzarin P, Zambardi R. A finite-volume-energy based approach to predict the static and fatigue behavior of components with sharp V-shaped notches. Int J Fract 2001;112:275-98. Sih GC. Strain-energy-density factor applied to mixed mode crack problems. Int J Fract 1974; 10: 305-21. Ayatollahi MR, Moghaddam MR, Berto F. A generalized strain energy density criterion for mixed mode fracture analysis in brittle and quasi-brittle materials. Theor Appl Fract Mech 2015;79:70-76. Inoue R, Yang JM, Kakisawa H, Kagawa Y. Mixed-mode fracture criterion of short carbon fiber-dispersed SiC matrix composite. J Ceram Sci Technol 2017;8:223-32. Mohtarami E, Baghbanan A, Hashemolhosseini H, Bordas SPA. Fracture mechanism simulation of inhomogeneous anisotropic rocks by extended finite element method. Theor Appl Fract Mech 2019;104:102359 [online]. Fayed AS. Numerical analysis of crack initiation direction in quasi-brittle materials: effect of T-stress. Arab J Sci Eng 2019;44:7667-76. Liu HY. Wing-crack initiation angle: A new maximum tangential stress criterion by considering T-stress. Eng Fract Mech 2018;199:380-91. Vásárhelyi B, Bobet A. Modeling of crack initiation, propagation and coalescence in uniaxial compression. Rock Mech Rock Eng 2000;33:119-39. Ning Y, Yang J, An XM, Ma GW. Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework. Comput Geotech 2011;38(1):40-49. Bai QS, Tu SH, Zhang C. DEM investigation of the fracture mechanism of rock disc containing hole(s) and its influence on tensile strength. Theor Appl Fract Mech 2016;86:197-216. Wang YN, Ha HB, Giang DN, Ranjith PG. A new SPH-based continuum framework with an embedded fracture process zone for modelling rock fracture. Int J Solids Struct 2019; 159:40-57. Zhou XP, Yang HQ. Multiscale numerical modeling of propagation and coalescence of multiple cracks in rock masses. Int J Rock Mech Min Sci 2012;55:15-27. 30

[54] Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock Eng 2015;48:1097-14. [55] Zhou XP, Wang YT. Numerical simulation of crack propagation and coalescence in pre-cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics. Int J Rock Mech Min Sci 2016;89:235-49. [56] Wang YT, Zhou XP, Shou YD. The modeling of crack propagation and coalescence in rocks under uniaxial compression using the novel conjugated bond-based peridynamics. Int J Mech Sci 2017;128-129:614-43. [57] Yang SQ, Tian WL, Huang YH. Failure mechanical behavior of pre-holed granite specimens after elevated temperature treatment by particle flow code. Geothermics 2018;72: 124-137. [58] Mumford D, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems. Commun Pure Appl Math. 1989;42(5):577-685. [59] Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via t-convergence. Commun Pure Appl Math. 1990;43(8):999-1036. [60] Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46(8):1319-42. [61] Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. J Elast 2008;91:5-148. [62] Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids, 2009;57(8):1209-29. [63] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int J Numer Methods Eng 2010;83(10):1273-1311. [64] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng 2010;199(45):2765-78. [65] Borden MJ, Hughes TJR, Landis CM, Verhoosel CV. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Comput Methods Appl Mech Eng 2014;273:100-18. [66] Kiendl J, Ambati M, De Lorenzis L,Gomez H, Reali A. Phase-field description of brittle fracture in plates and shells. Comput Methods Appl Mech Eng 2016;312:374-94. [67] Zhang P, Hu XF, Wang XY, Yao W. An iteration scheme for phase field model for cohesive fracture and its implementation in Abaqus. Eng Fract Mech 2018;204:268-87. [68] Carlsson J, Isaksson P. Dynamic crack propagation in wood fibre composites analysed by high speed photography and a dynamic phase field model. Int J Solids Struct 2018;144-145:78-85. [69] Zhou SW, Rabczuk T, Zhuang XY. Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies. Adv Eng Softw 2018;122:31-49. [70] Ren HL, Zhuang XY, Anitescu C, Rabczuk T. An explicit phase field method for brittle dynamic fracture. Comput Struct 2019;217:45-56. [71] Zhang P, Hu XF, Yang ST, Yao WA. Modelling progressive failure in multi-phase materials using a phase field method. Eng Fract Mech 2019;209:105-24. [72] Li KC, Zhou SW. Numerical investigation of multizone hydraulic fracture propagation in porous media: New insights from a phase field method. J Nat Gas Sci Eng 2019;66:42-59. 31

[73] Schneider D, Schoof E, Huang Y, Selzer M, Nestler B. Phase-field modeling of crack propagation in multiphase systems. Comput Methods Appl Mech Engrg 2016;312:186-95. [74] Paggi M, Reinoso J. Revisiting the problem of a crack impinging on an interface:A modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model. Comput Methods Appl Mech Engrg 2017;321:145-72. [75] Hansen-Dörr AC, de Borst R, Hennig P, Kästner M. Phase-field modelling of interface failure in brittle materials. Comput Methods Appl Mech Engrg 2019;346: 25-42. [76] Tushar KM, Vinh PN, Amin H. Phase feld and gradient enhanced damage models for quasi-brittle failure: A numerical comparative study. Eng Fract Mech 2019; 207:48-67. [77] Ayatollahi MR , Aliha MRM. Analysis of a new specimen for mixed mode fracture tests on brittle materials. Eng Fract Mech 2009;76:1563-73. [78] Williams ML. On the Stress distribution at the base of a stationary crack. J Appl Mech 1957; [79] [80]

[81]

[82]

[83] [84]

[85] [86] [87]

[88]

[89]

[90]

24:109-14. Erdogan Ayatollahi MR, Aliha MRM. Mixed mode fracture in soda lime glass analyzed by using the generalized MTS criterion. Int J Solids Struct 2009;46:311-21. Mohtarami E, Baghbanan A, Hashemolhosseini H. Prediction of fracture trajectory in anisotropic rocks using modified maximum tangential stress criterion. Comput Geotech 2017; 92:108-20. Seitl S, Miarka P, Bílek V. The mixed-mode fracture resistance of C 50/60 and its suitability for use in precast elements as determined by the Brazilian disc test and three-point bending specimens. Theor Appl Fract Mech 2019;97:108-18. Wei MD, D F, Zhou JW, Liu Y, Luo J. A further improved maximum tangential stress criterion for assessing mode I fracture of rocks considering non-singular stress terms of the Williams expansion. Rock Mech Rock Eng 2018;51:3471-88. Taylor D, Merlo M, Pegley R, Cavatorta MP. The effect of stress concentrations on the fracture strength of polymethyl methacrylate. Mater Sci Eng A 2004;382:288-94. Wei MD, Dai F, Xu NW, Liu Y, Zhao T. Fracture prediction of rocks under mode I and mode II loading using the generalized maximum tangential strain criterion. Eng Fract Mech 2017;186:21-38. Zhang X, Vignes C, Sloan SW, Sheng D. Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale. Comput Mech 2017;59:737-52. Yin BB, Zhang LW. Phase field method for simulating the brittle fracture of fiber reinforced composites. Eng Fract Mech 2019;211:321-40. Aliha MRM, Ayatollahi MR, Smith DJ, Pavier MJ. Geometry and size effects on fracture trajectory in a limestone rock under mixed mode loading. Eng Fract Mech 2010;77:2200-2212. Ban H, Im S, Kim YR. Mixed-mode fracture characterization of fine aggregate mixtures using semicircular bend fracture test and extended finite element modeling. Constr Build Mater 2015; 101:721-729. Thiercelin M, Roegiers JC. Fracture toughness determination with the modified ring test. In: Hartman HL, editor. Proc 27th US Symp on Rock Mech, Littleton, CO., SME; 1986. p:615-22. Thiercelin M. Fracture toughness under confining pressure using the modified ring test. In: Rotterdam, Proc 28th US Sym on Rock Mech, Balkema, CO., SME; p:149-56. 32

[91] Thiercelin M. Fracture toughness and hydraulic fracturing. Int J Rock Mech Min Sci Geomech Abstr 1989;26:177–83. [92] Tutluoglu L, Keles C. Effects of geometric factors on mode I fracture toughness for modified ring tests. Int J Rock Mech Min Sci 2012;51(1):149-61.

33

Fig.1 Phase field parameter as a function of spatial coordinates

Fig.2 A diagonally loaded square plate (DLSP) specimen

34

Fig.3 Comparisons of the specimen failure modes between experimental results and phase field simulation results

Fig.4 Comparisons of the peak load between the experimental results and phase field simulation results

35

(a) CCCD

(b) SCB

(c) HCCD

Fig.5 Three disc-type specimens

Y



r



rr

X

Crack tip Fig.6 The stress and displacement fields around the crack tip

36

Fig.7 The definition of the crack initiation angle 0

(a) Load-displacement curves

(b) Peak load

Fig.8 Influence of mesh density on numerical results

37

(b)

(a)

(c)

Fig.9 Dimensionless SIFs and T-stress of three disc-type specimens

(a) CCCD [87]

(b) HCCD [21]

(c) SCB [88] Fig.10 Crack path comparison between numerical results and experimental 38

observation

u = 0.041mm

u = 0.043mm

u = 0.0451mm

u = 0.0465mm

u = 0.0482mm

u = 0.043mm

u = 0.045mm

u = 0.047mm

u = 0.048mm

u = 0.0495mm

u = 0.0447mm

u = 0.0467mm

u = 0.0484mm

u = 0.0491mm

u = 0.0504mm

u = 0.045mm

u = 0.0469mm

u = 0.0487mm

u = 0.0497mm

u = 0.0519mm

u = 0.0459mm

u = 0.0488mm

u = 0.0497mm

u = 0.0502mm

u = 0.0516mm

u = 0.0458mm

u = 0.0481mm

u = 0.0496mm

u = 0.0503mm

u = 0.0515mm

Me = 0

Me = 0.2

Me = 0.4

Me = 0.6

Me = 0.8

Me = 1.0

Fig.11 Crack propagation of CCCD specimens with different mode mixities Me

39

Fig.12 Comparison between theoretical and numerical crack initiation angles for CCCD specimens

Fig.13 Load-displacement curves for CCCD specimens

40

Fig.14 Mixed mode fracture toughness obtained by the phase field method and theoretical criteria in CCCD specimens

41

u = 0.0219mm

u = 0.0232mm

u = 0.0242mm

u = 0.025mm

u = 0.0268mm

u = 0.0204mm

u = 0.0218mm

u = 0.0228mm

u = 0.0236mm

u = 0.0279mm

u = 0.0202mm

u = 0.0218mm

u = 0.0228mm

u = 0.0236mm

u = 0.0263mm

u = 0.0173mm

u = 0.0189mm

u = 0.0202mm

u = 0.0212mm

u = 0.0233mm

u = 0.0170mm

u = 0.0186mm

u = 0.0198mm

u = 0.0208mm

u = 0.022mm

u = 0.0162mm

u = 0.018mm

u = 0.0193mm

u = 0.0205mm

u = 0.0230mm

Me = 0

Me = 0.2

Me = 0.4

Me = 0.6

Me = 0.8

Me = 1

Fig.15 Failure mode of SCB specimens with different mode mixities Me

42

Fig.16 Comparison between theoretical and numerical crack initiation angles for SCB specimens

Fig.17 Load-displacement curves for SCB specimens

43

Fig.18 Mixed mode fracture toughness of SCB specimens obtained by the phase field method and theoretical criteria

44

u = 0.0219mm

u = 0.0232mm

u = 0.0242mm

u = 0.025mm

u = 0.0268mm

u = 0.0204mm

u = 0.0218mm

u = 0.0228mm

u = 0.0236mm

u = 0.0279mm

u = 0.0202mm

u = 0.0218mm

u = 0.0228mm

u = 0.0236mm

u = 0.0263mm

u = 0.0173mm

u = 0.0189mm

u = 0.0202mm

u = 0.0212mm

u = 0.0233mm

u = 0.0170mm

u = 0.0186mm

u = 0.0198mm

u = 0.0208mm

u = 0.022mm

u = 0.0162mm

u = 0.018mm

u = 0.0193mm

u = 0.0205mm

u = 0.0230mm

Me = 0

Me = 0.2

Me = 0.4

Me = 0.6

Me = 0.8

Me = 1

Fig.19 Failure mode of HCCD specimens with different mode mixities Me

45

Fig.20 Comparison between theoretical and numerical crack initiation angles for HCCD specimens

Fig.21 Load-displacement curve for HCCD specimens

46

Fig.22 Mixed mode fracture toughness obtained by the phase field method and fracture criteria for HCCD specimens

Fig.23 Comparison of crack initiation angles of three disc specimens

47

(a)

(b)

Fig.24 The trend of mode I and mode II fracture toughness of three kinds of disc specimens

(a)

(b)

Fig.25 Comparison of crack initiation time for the three disc specimens

48

Table 1 Numbers of elements and nodes in the models with different mesh densities Mesh size h

l0 / 5

l0 / 4

l0 / 3

l0 / 2

l0

2 l0

Number of elements

36877

24721

15321

8043

3200

1278

Number of nodes

37017

24843

15429

7953

3276

1332

Table 2 The values of  for the required mode mixities in three disc-type specimens  (degree)

specimen Me = 1

Me = 0.8

Me = 0.6

Me = 0.4

Me = 0.2

Me = 0

CCCD

0

4.4

8.8

14.2

19.8

27

SCB

0

12.2

23.1

30.9

37.2

42.5

HCCD

0

5.6

10.3

15.2

20.6

26.5

49

Highlights (1) Phase field method is used to study the mixed mode fracture behaviors of disc-type specimens without external criterion. (2) As the crack inclination angle transits from mode I to pure mode II, mode I fracture toughness decreases, while mode II fracture toughness increases. (3) The numerical results obtained by phase field method are in good agreement with those predicted by GMTS criterion.

50

Highlights (1) Phase field method is used to study the mixed mode fracture behaviors of disc-type specimens without external criterion. (2) As the crack inclination angle transits from mode I to pure mode II, mode I fracture toughness decreases, while mode II fracture toughness increases. (3) The numerical results obtained by phase field method are in good agreement with those predicted by GMTS criterion.

51