Theoretical and numerical prediction of crack path in the material with anisotropic fracture toughness

Engineering Fracture Mechanics 180 (2017) 330–347

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Theoretical and numerical prediction of crack path in the material with anisotropic fracture toughness Yue Gao, Zhanli Liu ⇑, Qinglei Zeng, Tao Wang, Zhuo Zhuang, Keh-Chih Hwang AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China Center for Mechanics and Materials, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 1 May 2017 Received in revised form 8 June 2017 Accepted 8 June 2017 Available online 16 June 2017 Keywords: Anisotropic material Crack growth Energy release rate Fracture toughness Extended finite element method (XFEM)

a b s t r a c t The crack path in the anisotropic medium is studied theoretically and numerically in this paper, with focusing on the effects of the anisotropic fracture toughness. A weak plane model is adopted to characterize the anisotropic fracture toughness and the maximum energy release rate criterion (MERR) is chosen to predict the crack path. We prove that the crack deflecting direction in a weak plane material, theoretically, only relates to the weak plane direction, the weakness of the plane, and the ratio of stress intensity factors before crack extending. Two critical weakness ratios are found: one is for that the weak plane captures all cracks for any loading state; the other is for that the weak plane traps the crack once it is deflected to the plane. To further numerically study the complex crack path in an anisotropic medium, the extended finite element method (XFEM) embedded MERR criterion and weak plane model is developed, in which a mesh independent piecewise linear crack formulation is proposed to capture the curved crack path. By modeling the three point bending loading test, the influences of the weak plane angles and the weaknesses on crack path deflection are studied, and periodically oscillatory crack path behaviors are found as observed in experiments. Further more, the crack extending through multiple layers of weak plane material and isotropic material is numerically studied, which shows that the crack diffracts at the material interface, and the crack path in the anisotropic toughness material is rougher than in the isotropic material. The presented work in this paper will be helpful to understand and control the crack path in the rocks as well as other anisotropic materials. Ó 2017 Published by Elsevier Ltd.

1. Introduction Anisotropic materials are commonly used in the world, such as rocks, crystals, and many other nature materials. However, their fracture behaviors have not been fully understood. In this article, we will focus on the crack path in the anisotropic material, and demonstrate a theoretical prediction of it with numerical verifications, then study the complex fracture behavior in anisotropic material numerically. When we state that a medium is anisotropic, we usually mean that the constitutive model in the medium is anisotropic. However, there is another type of anisotropy, the anisotropic fracture toughness. Nasseri and Mohanty [1] measured fracture

⇑ Corresponding author at: Center for Mechanics and Materials, Tsinghua University, Beijing 100084, China; AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China. E-mail address: [email protected] (Z. Liu). http://dx.doi.org/10.1016/j.engfracmech.2017.06.013 0013-7944/Ó 2017 Published by Elsevier Ltd.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

331

toughness for granitic rocks in different directions. They found that orientations of the grain axis and microcracks might be various in different faces of a rock specimen. Kataoka et al. [2] measured two types of anisotropic rocks, African granodiorite and Korean granite. They found that, the fracture toughness might be 29% larger in the direction perpendicular to the bedding plane than parallel to it. Chandler et al. [3] measured Mancos shale, and they found that the fracture toughness is 0:72 MPa m1=2 when the crack is normal to the bed plane of shale, but only 0:21 MPa m1=2 when parallel to the bed plane in their specimen. Besides rocks, the anisotropic fracture toughness are also widely observed in other materials, e.g., NiAl single crystals [4], graphene aerogels infiltrated with epoxy resin [5], microcellular foamed ABS polymer [6], polyurethane and polyisocyanurate foams [7]. Compared with isotropic material, the fracture paths in the material or composite structures with anisotropic fracture toughness are much more complicated, but also more interesting. Takei et al. [8] performed tearing experiments on polypropylene sheets, and validated the anisotropic fracture toughness model. They pointed out that for a well-designed fracture toughness profile, there would be forbidden directions for crack extending. Audoly et al. [9] and Reis et al. [10] observed an oscillatory fracture when a cylinder cuts a brittle elastic film in experiments. Nam et al. [11] showed results of controlled crack paths in a multilayer medium, which is combined with a silicon nitride Si3N4 film, a single-crystalline silicon substrate, and a silicon wafer. Three different crack path patterns were found: straight, oscillatory, and orderly bifurcated, by controlling the system parameters. Nasseri and Mohanty [1] and Kataoka et al. [2] also found that the cracks would extend more roughly in the larger fracture toughness direction, when they measured the fracture toughness of granitic rocks in laboratory. There have already been some theoretical and numerical works to predict the crack path in anisotropic materials. Theoretically, plenty of criteria are developed to predict crack extending direction, such as (1) the principle of local symmetry (K II ¼ 0), (2) maximum K I criterion, (3) maximum hoop stress criterion, (4) maximum energy release rate (MERR) criterion. These criteria have been studied well in the isotropic constitutive model, and they usually give the same prediction in an isotropic medium. However, in an anisotropic medium, the crack extending directions predicted by them are usually different [12]. There are many people who endorsed the MERR criterion, and warned that other criteria may give incorrect predictions even in isotropic materials [e.g., 8,12–15]. Zeng and Wei [16] used MERR to analyze the initiation of hydraulic fracture deflection when it crosses a given nature fracture with weak strength, based on the study of He and Hutchinson [17] and Hutchinson and Suo [18]. The MERR criterion is generally accepted in anisotropic material fracture study [3,19], and is adopted in this paper. Recently, Li et al. [20] also developed a microplane constitutive model to consider the anisotropic mechanics behavior of shale rocks. From the aspect of numerical simulation, phase field fracture model is often used on the anisotropic fracture toughness topic. Hakim and Karma [12,14] demonstrated a force-balance principle (FL) and compared the results with maximum energy release rate criterion (MERR). They found out that the kink angle in their phase field simulation result was between the predictions of FL and MERR. Li et al. [19] focused on the anisotropic fracture toughness that has rotational symmetry of order four, which looks like a flower type profile. A zigzag crack extending between two constrained edges is simulated. Shanthraj et al. [21] developed a finite-strain anisotropic phase field model to study the kinking and branching of the crack path resulting from the crystallographic misorientation across the laminate boundary and the grain boundaries. Though phase field method is widely used in modeling complex fracture phenomenon recently, it is hard to extract the fracture parameters from the phase field modeling. Li et al. [22] formulated a discrete micromechanical approach based on the lattice discrete particle model to simulate the anisotropic mechanical behavior of shale. The extended finite element method (XFEM) is a powerful tool in numerical fracture modeling [23–27]. With this method, mesh refinement is not required when crack extends, and stress intensity factors are easy to be obtained with interaction integrals [28]. To the best knowledge of the authors, there is yet no work dealing with the material with anisotropic fracture toughness using XFEM. The weak plane model is introduced in this paper as a special case of the anisotropic fracture toughness models. The model assumes that the fracture toughness is the same in almost every direction as an isotropic material, but is lower in a specific direction, which is the direction of the weak plane. The sedimentary rock, like shale, is the typical material which obeys the weak plane model and will be of particular interest in this paper. In this study, the MERR criterion is adopted to predict the crack path in an anisotropic medium, with the weak plane model to simulate the anisotropic fracture toughness. A simple crack deflection criterion is demonstrated from the weak plane model, and the forbidden areas phenomenon of crack deflection is discussed. An XFEM based numerical modeling technique, which embeds the MERR criterion and the weak plane model to handle the anisotropic fracture toughness, is developed to study the complex crack path in anisotropic material, and the results between theoretical prediction and numerical simulation are compared. 2. Theoretical model 2.1. Maximum energy release rate criterion The MERR criterion is adopted in this paper to predict crack extending direction, and is briefly introduced in this section.

332

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

We would like to focus on the fracture toughness anisotropy. Therefore, the fracture toughness Gc ðhÞ is assumed to be a function of the crack extending direction h. The original MERR criterion in the isotropic material states as follows [13]. Let the crack driving force G be a function, which only relates to the crack extending direction h in the next moment. Then the crack path after extending would satisfy

maxfGðhÞg:

ð1Þ

h

In the material with anisotropic fracture energy Gc ðhÞ, the MERR criterion1 is extended to that the crack path after extending should satisfy [12,17]


 GðhÞ : max h Gc ðhÞ

ð2Þ

If GðhÞ and Gc ðhÞ are smooth enough, a necessary condition of Eq. (2) could be obtained as


 d GðhÞ ¼ 0: dh Gc ðhÞ

ð3Þ

However, the weak plane model formulated in this work does not satisfy the continuity requirement in Eq. (3), see Section 2.2, thus only Eq. (2) is used in this paper as the crack direction criterion. To use the MERR criterion, we must figure out this problem: assume there is a crack whose stress intensity factors (SIFs) ~ I ðhÞ; K ~ II ðhÞg after the in mode I and mode II are known as, separately, K I and K II , then we want to know what are the SIFs fK crack propagating in the direction of h, with an small enough length. Nuismer [29] studied this problem in an isotropic med~ I ðhÞ; K ~ II ðhÞg should be a linear combination of fK I ; K II g for a given h. He and Hutchinson [30] studium, and found out that fK ied this problem when a crack penetrates the interface of two different isotropic material. He and Hutchinson [31] numerically calculated this problem, and a table of more accurate coefficients in the linear relations has been obtained. Amestoy and Leblond [32] used a complex analysis method to study this problem, and obtained series expansion results of these coefficients related to h. With the results obtained by Amestoy and Leblond [32], the crack driving force could be demonstrated as

Gðh; K I ; K II Þ ¼

i 1  m2 h ~ 2 ~ 2 ðhÞ ; K I ðhÞ þ K II E

ð4Þ

~ I ðhÞ; K ~ II ðhÞg and fK I ; K II g are shown in Appendix A. and the full relations between fK The results obtained from Nuismer [29], He and Hutchinson [30,31], and Amestoy and Leblond [32] are compared in Appendix B, from which it could be found that the results in Amestoy and Leblond [32] are accurate enough, and these expressions will be adopted in this article, instead of the tabular data from He and Hutchinson [31], or the results from Nuismer [29], which are too simple and not appropriate to be used in crack with a large deflection angle. In the isotropic material, the crack deflection angle could be obtained by

dGðh; K I ; K II Þ ¼ 0; dh

ð5Þ

and by using Eqs. (4), (A.1) and (A.2), it could be proved that the obtained deflection angle h0 , in isotropic medium, is only related to the ratio of SIFs (K II =K I ). Thus, the function could be written as

h0 ¼ h0

  K II ; KI

ð6Þ

which is also proved by Hutchinson and Suo [18] based on the results of He and Hutchinson [30,31]. 2.2. Weak plane model Hoek [33] proposed a classical image of rock material as follows. There are two types of cracks in a rock: one is the primary cracks, which usually lie along the bedding plane; the other is the secondary cracks, which are actually randomly oriented grain boundaries; see Fig. 1a for details. This image inspired us an anisotropic fracture toughness model, which has different fracture energies in different directions. In most of the directions, the rock has a usual fracture toughness Gc0 , which is related to the randomly oriented grain boundaries, therefore has the same value, independent of h. In the direction of bedding plane, the rock is weaker and has a lower toughness Gcw . The weak plane fracture toughness model is discontinuous, and fracture energy Gc ðhÞ could be describes as




Gc ðhÞ ¼

1

Gcw

h ¼ a or

Gc0

other

pþa

;

In order to clarify, the phrase MERR criterion in this paper below indicates the MERR criterion with anisotropic fracture toughness Gc ðhÞ, i.e. Eq. (2).

ð7Þ

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

333

Fig. 1. (a) A schematic diagram of fracture inside a rock with weak plane. The hatched diagonal lines denote the weak plane direction in the medium, but not real discrete weak plane positions. (b) Fracture energy Gc ðhÞ of the weak plane model relates to crack extending direction h. a denotes the weak plane direction, where counterclockwise is positive. Gc ðhÞ has the value Gc0 almost in every direction except the weak plane direction a, where Gc ðhÞ has a weaker value Gcw .

where a is the weak plane direction, see Fig. 1b. It is noted that, the weak plane model does not implicate multiple discrete weak planes in the medium, but means that there is a weak plane direction in everywhere of the material. This simple weak plane fracture toughness model is not only proper in rocks, but also in many crystal materials, which have a similar fracture behavior. There are some advantages in the proposed weak plane model. Firstly, the fracture toughness parameters of the material are easy for measurement. The weak plane direction is the same as the bedding plane direction, and standard three point bending test could be applied in the bedding plane direction and corresponding perpendicular direction, separately, to measure Gcw and Gc0 . Secondly, although the discontinuity makes the model inconvenient in using some mathematical techniques, e.g. Eq. (3), it is still easy to be embedded in numerical modeling with MERR criterion, by following three steps. Step 1. For any given crack with known fK I ; K II g, we have Gðh; K I ; K II Þ for any crack extending direction h, which could be obtained from Eqs. (4), (A.1) and (A.2). Step 2. Assume the crack is extending in an isotropic fracture toughness material, find h0 ¼ arg maxh GðhÞ, which could be obtained from Eq. (5). Step 3. Compare Gðh0 Þ=Gc0 with GðaÞ=Gcw . If GðaÞ=Gcw > Gðh0 Þ=Gc0 , then the crack will deflect to h ¼ a, otherwise it still goes to h0 . The values compared in the step 3 could be rewritten as GðaÞ=Gðh0 Þ and Gcw =Gc0 , and we have Theorem 1. The crack deflection direction could be obtained by

Gða; K I ; K II Þ Gcw ) crack deflects to h ¼ a; > Gðh0 ; K I ; K II Þ Gc0 Gða; K I ; K II Þ Gcw ) crack deflects as in an isotropic material; to h0 ¼ arg maxh GðhÞ; < Gðh0 ; K I ; K II Þ Gc0

ð8Þ

where h0 could be obtained from Eq. (5). A similar result has been found by He and Hutchinson [17], when they dealt with the problem whether the crack would penetrate or deflect to the interface of two different elastic materials. Then the value compared in the left hand of Eq. (8), i.e. ½Gða; K I ; K II Þ=Gðh0 ; K I ; K II Þ, is only related to the angle between the crack and the weak plane, as well as the load configurations, while the right hand Gcw =Gc0 is only related to the material’s properties. 2.2.1. Crack path prediction in pure Mode I and Mode II In this section, two examples are given to show the convenience of using weak plane model. Firstly, assume there is a pure mode I type crack, and in the isotropic material, the crack would extend directly, therefore h0 ¼ 0. Then the left hand value in Eq. (8) could be simplified by using Eqs. (4) and (A.1) as

    ~ 2 ðaÞ ~ 2 ðaÞ þ K F 2 pa þ F 221 pa Gða; K I Þ K II ; ¼ I ¼ 11 ~ I ðh0 Þ þ K ~ 2 ðh0 Þ F 2 ð0Þ þ F 2 ð0Þ Gðh0 ; K I Þ K 11 21 II

ð9Þ

334

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

where F ij ðmÞ’s expressions are shown in Eq. (A.2). Equation (9) is a function only related to the weak plane direction a, since K I is eliminated. Therefore, the mode I crack path is only determined by the direction of weak plane a, and the weakness of the weak plane Gcw =Gc0 . For any boundary loading state resulting to pure mode I crack, a figure of energy release rate ratio GðaÞ=Gðh0 Þ in Eq. (9) is shown in Fig. 2, where h0 ¼ 0. A horizontal line denoting the weak plane property Gcw =Gc0 is drawn in Fig. 2 to show an example of crack path prediction. Once we have a weak plane direction a in the real case, we can find out the corresponding point on the mode I curve in Fig. 2. And according to Eq. (8), if this point is above the horizontal line Gcw =Gc0 , then the crack will extend to the weak plane direction a, otherwise extend to h0 ¼ 0. In the example of Gcw =Gc0 ¼ 0:6 shown in Fig. 2, if the weak plane direction jaj 6 0:317p  56:97 , then the crack is attracted to the weak plane. On the other hand, if jaj > 56:97 , the crack path will just extend straightly instead. There is a minimum value GðaÞ=Gðh0 Þ  0:2594 for mode I crack on Fig. 2, when a ¼ p=2. Thus for a material in which the weak plane is weak enough such that Gcw =Gc0 < 0:2594, the crack would always extend to the weak plane for any material direction, in the mode I loading case. For mode II cracks, Eq. (5) finds out that the maximized energy release rate direction is h0  75:7566 , which is a bit pffiffiffi

larger in value than the result from maximum hoop stress principle h ¼ 2 arctan 2=2  70:5288 . The result is different from the local symmetry principle either, which predicts a kink angle h  77:3301 . This difference has been discussed by Amestoy and Leblond [32] and Hutchinson and Suo [18]. As in mode I crack, a pure mode II crack could also be found that the left part of Eq. (8) has nothing to do with K II . By letting K I ¼ 0, the ratio of energy release rate could be obtained as

    F 212 pa þ F 222 pa Gða; K II Þ ¼    : Gðh0 ; K II Þ F 212 hp0 þ F 222 hp0

ð10Þ

The result of mode II crack is also shown in Fig. 2. In the example of Gcw =Gc0 ¼ 0:6, the entire curve is above the horizontal line, thus crack would extend to weak plane direction for any a. In fact, the minimum value is 0:6621 at a ¼ 0. 2.2.2. Crack path prediction in any loading case A key technique we described in mode I and mode II crack path prediction is the elimination of K I and K II , so that the curve of GðaÞ=Gðh0 Þ is irrelevant to an actual loading case, which makes crack path prediction only relate to material direction a and weakness property Gcw =Gc0 . For any other loading cases, i.e. K I – 0 and K II – 0, it could be demonstrated that

h     i2 h  a  K   i2 F 11 pa þ KKIII F 12 pa þ F 21 p þ KIII F 22 pa Gða; K I ; K II Þ ¼h     i2 h  h  K   i2 ; Gðh0 ; K I ; K II Þ F 11 hp0 þ KKIII F 12 hp0 þ F 21 p0 þ KIII F 22 hp0

ð11Þ

where h0 is obtained from Eq. (5) and only related to K II =K I . Thus we have Theorem 2. The crack extending direction in the weak plane model is only related to the following three parameters: (1) weak plane direction a, (2) weakness ratio of the weak plane Gcw =Gc0 , (3) ratio of SIFs before crack extending K II =K I . A surface of GðaÞ=Gðh0 Þ varying with a and K II =K I is shown in Fig. 3. The figure only shows the K II =K I > 0 part, and the other part is symmetrical to the z-axis since

    K II K II ¼ G a; ; G a;  KI KI

ð12Þ

because of the odd function property of F 12 ðmÞ and F 21 ðmÞ. For a given material with known Gcw =Gc0 , we could draw a horizontal plane in Fig. 3. Then all the information we need to determine a crack extending direction are in the figure. Any cases in which the loading state K II =K I and weak plane direction a leads to the GðaÞ=Gðh0 Þ surface above the horizontal plane Gcw =Gc0 , would let the crack extend to weak plane direction a, otherwise to the isotropic direction h0 obtained by Eq. (5). 2.2.3. Crack path deflection and forbidden area If we know the weak plane direction a but do not know the loading details on the boundaries or the SIFs {K I ; K II } at crack tip, then vertically slicing the Fig. 3 with the known a is a way for crack path prediction.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

335

Fig. 2. The energy release rate ratio GðaÞ=Gðh0 Þ varies by the weak plane direction a, for pure mode I and mode II cracks separately. The horizontal line denotes a known Gcw =Gc0 ¼ 0:6 of given material, and any part of the curve above the horizontal line means corresponding crack is attracted by the weak plane and would extend to a, otherwise to h0 . For mode I crack, h0 ¼ 0; for mode II, h0  75:76 .

Fig. 3. Energy release rate ratio GðaÞ=Gðh0 Þ for any weak plane direction a and loading state K II =K I .

Figure 4 shows five sliced curve for a ¼ fp=2; p=4; 0; p=4; p=2g separately. For a given material direction a, cases for both K II =K I > 0 and K II =K I < 0 should be considered. Therefore, we need to take care of Fig. 4, and watch both a and a curves for a given a. For a known material direction a, we could find the easiest loading case that leads the crack to deflect to weak plane in Fig. 4, as well as the most difficult case. For example, when a ¼ p=4, the highest point on the curve a ¼ p=4 is at K II =K I ¼ 0:5734, therefore loading cases nearby K II =K I  0:5734 would help the weak plane attract crack deflection. On the other hand, the lowest point on the curve a ¼ p=4 indicates that K II =K I  0:7578 would make the crack just extend to isotropic direction unless Gcw =Gc0 < 0:2376. In the isotropic material, crack extending direction h0 is determined by K II =K I , and a monotonic function h0 ðK II =K I Þ could be obtained from Eq. (5) with the MERR criterion. However, in the material with a weak plane, the energy release ratio GðaÞ=Gðh0 Þ in Fig. 4 is cut by the weakness Gcw =Gc0 horizontal line, and is divided into multiple parts. A similar cutting is performed in Fig. 2 for the pure mode I crack situation. All the parts that above this horizontal line, which has a corresponding original deflection direction in the isotropic material, will deflect to the weak plane direction now. For instance, the dashed line part at a ¼ p=4 in Fig. 4 is cut by Gcw =Gc0 ¼ 0:9, where 0:1736 < K II =K I < 4:789 is cut down. The crack path will never

336

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

Fig. 4. Energy release rate ratio GðaÞ=Gðh0 Þ varies by loading case K II =K I , for given material weak plane directions a ¼ fp=2; p=4; 0; p=4; p=2g. The red dashed line denotes the toughness ratio Gcw =Gc0 of a given material. As an example, the curve a ¼ p=4 is cut by Gcw =Gc0 ¼ 0:9, and the cut part above the straight line is plotted in dashed blue line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

deflect to directions in these cut area but to weak plane direction, and we would call it the crack path forbidden area phenomenon. For a given material with certain weak plane direction a and the weakness Gcw =Gc0 , the critical values of K II =K I for the forbidden area boundaries could be obtained theoretically by solving equation Eq. (11) with the given a and Gcw =Gc0 . And the corresponding critical deflection angles could be obtained by Eq. (5). For instance, the deflection angles for a material with weak plane at a ¼ p=12 ¼ 15 and weakness Gcw =Gc0 ¼ 0:9 are shown in Fig. 5. Nine contours K II =K I ¼ f2; 1; 0:5; 0:25; 0; 0:25; 0:5; 1; 2g are marked on the figure with corresponding deflection angles. The forbidden area is plotted as an orange area in this figure, which shows the minimum and maximum boundaries that a crack would deflect to the weak plane. Fig. 5 shows that the loading cases with 0:0988 6 K II =K I 6 0:472 will lead the crack to deflect to the weak plane, which is obtained from Fig. 3 with a vertical slice plane in a ¼ p=12. And any deflecting angle within 38:95 6 h 6 11:07 in an isotropic medium, will now deflect to a ¼ 15 in the weak plane material. The two critical boundary angles are obtained from Eq. (5) with corresponding K II =K I . Therefore, we find out that in this weak plane material, only crack deflecting angle that satisfies h 2 ½75:7566 ; 38:95 Þ [ f75:7566 g[ ð11:07 ; 75:7566  is allowable, and all the loading cases leading to h 2 ½38:95 ; 11:07 are now captured by the weak plane a ¼ 15 . Thus we call the deflecting area ½38:95 ; 15 Þ [ ð15 ; 11:07  is a forbidden area for crack deflecting. The minimum energy release rate ratio value with respect to K II =K I


 Gða; K II =K I Þ min K II =K I Gðh0 ; K II =K I Þ

ð13Þ

is a function only related to a, and is plotted in Fig. 6. This means that for a given weak plane direction a, if the weakness Gcw =Gc0 is lower than the corresponding minimum value in Fig. 6, the crack would always extend to the weak plane for any loading boundaries for any K II =K I . The smaller the angle between the crack and the weak plane, the higher is critical value. Once the weakness Gcw =Gc0 in the weak plane is lower than the value given in Fig. 6 for the corresponding angle a, the forbidden area would extend to the whole possible directions and only remains the weak plane direction. In this case, the orange area in Fig. 5 would cover the entire contour plot. The minimal value in Fig. 6 is

 GðaÞ  ¼ 0:0580346; Gðh0 Þ a ¼ p=2 K II ¼ 0:39K I

ð14Þ

and any weak plane with Gcw =Gc0 less than this value will capture the crack in any loading state. The maximum value in Fig. 6 is

 GðaÞ  Gðh0 Þ

¼ 0:553617; a¼0 K II ¼ 2:02K I

ð15Þ

which means, once a crack is trapped in a weak plane with weakness Gcw =Gc0 < 0:554, then the crack could never escape from the controlling of the weak plane.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

337

Fig. 5. Polar contours of crack deflection angle for given K II =K I , which are outlined on the figure with corresponding angles. The weak plane parameters are set to a ¼ p=12 and Gcw =Gc0 ¼ 0:9. Orange area in the figure shows a forbidden area for crack deflection angle, because the deflection angles near a are attracted by the weak plane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3. Numerical techniques In this section, the MERR criterion and the weak plane model is composed into the extended finite element method (XFEM) to numerically study the complex crack path in an anisotropic material. In our modeling, full Heaviside enrichment functions and crack tip enrichment functions are used to obtain the fracture parameters as accurate as possible. To precisely capture the kink crack extending, a piecewise linear crack formulation is introduced into XFEM modeling, which separates crack into multiple segments, and every segment’s position, direction, and length is independent of element meshes. To use the piecewise linear crack formulation efficiently, the grid ray tracing method and horizontal ray method is introduced to our modeling, which are originally used in computer graphics. The details of XFEM and piecewise linear crack formulation are discussed in Sections 3.1 and 3.2. 3.1. Extended finite element method The extended finite element method is based on partitions of unity, and using additional enrichment functions with the basic shape functions of classical finite element method. Four-node elements of plane strain model are used in the program, and the trial displacement field in any direction could be expanded as

uðxÞ ¼

4 X X X X NI ðxÞuI þ N J ðxÞHðxÞaJ þ NK ðxÞ Ul ðxÞbKl ; I2S

J2Sh

K2Sc

l¼1

Fig. 6. The minimal energy release rate ratio GðaÞ=Gðh0 Þ for any loading case K II =K I , varies by the weak plane direction a.

ð16Þ

338

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

where S denotes the set of whole nodes, and Sh and Sc denotes Heaviside function enriched nodes and crack tip shape function enriched nodes separately. N I ; N J , and N K are basic shape functions of four-node element. HðxÞ denotes a Heaviside-like function, which could distinguish that a point with given coordinate x locates at which side of a crack. Since piecewise linear crack is formulated and the crack is not one straight line in every element, see Section 3.2, the often used Heaviside function with the level set Hðf ðxÞÞ could not be adopted here. The details of HðxÞ is shown in Appendix C.

Ul ðxÞ denotes the crack tip enrichment functions

UðxÞ ¼


 pffiffiffi pffiffiffi h pffiffiffi h pffiffiffi h h r sin ; r cos ; r sin sin h; r cos sin h ; 2 2 2 2

ð17Þ

where r and h denote the polar coordinates referenced to the crack tip’s location and direction. Besides uI , there are additional variables in the displacement fields aJ and bKl in Eq. (16), related to the Heaviside part and crack tip part separately. They are only activated in sets Sh and Sc . To calculate the SIFs K I and K II for a given displacement and stress field obtained by the XFEM modeling, the traditional interaction integral method is used [23,28,34], which states that

I

where

ð1;2Þ

n

Z "

# ð2Þ ð1Þ @ui @p ð2Þ @ui ð2Þ ð1Þ ¼ r þ rij  rmn emn d1j dA @xj @x1 @x1 A

2 ð1Þ ð2Þ ð2Þ ð1Þ ¼ 0 K I K II þ K I K II ; E ð1Þ ij

o

n

ð18Þ ð19Þ o

ð1Þ ð1Þ ð1Þ ð1Þ ð2Þ ð2Þ ð2Þ ð2Þ ð2Þ and rij ; eij ; ui ; K I ; K II denote any two different stress, strain, displacement, and correrð1Þ ij ; eij ; ui ; K I ; K II

sponding SIFs states separately. p denotes any scalar field which equals 1 near the tip, but equals 0 far away, thus the integral domain A could be chosen as only including elements where @p=@xj – 0. The x1 direction in Eq. (18) must be the crack propagation direction. In practice, to obtain SIFs K I and K II , let the state (1) be the solution obtained from XFEM modeling, and let the state (2) be a specific solution with known SIFs. Then integrate the expression in Eq. (18) with these two solutions, and obtain a equation between Eqs. (18) and (19). After changing the solution using in state (2), two different equations are obtained, and SIFs of state (1) could be solved. After obtaining SIFs K I and K II by interaction integrals, ½Gða; K I ; K II Þ=Gðh0 ; K I ; K II Þ can be calculated by Eqs. (11) and (A.2). Then the crack deflection angle in the next time step can be predicted by using Section 1. The details of the process has been described in the steps 1–3 in Section 2.2.

3.2. Mesh independent piecewise linear crack description The traditional XFEM for simplicity, usually letting the crack extend straightly to element boundary for every propagating step, and start from this boundary to another boundary in next step. Thus, the crack is constructed by multiple linear crack parts, which are straight in every element, and only deflect at the element boundaries. Obviously, this implementation limits the crack path, and binds it with element mesh. A mesh independent propagable crack is developed in our XFEM modeling by the piecewise linear crack description. Fig. 7 shows a schematic diagram of the piecewise linear crack. A square grid mesh is used in this example, though any types of mesh element are acceptable. The crack path is divided into multiple crack parts, whose position, length, and direction of every segment are totally independent to the mesh. The split point positions marked with crosses in Fig. 7 are the only information need to be stored to reconstruct the piecewise crack. The length of every crack part could be different, e.g., in Section 4.1 and Fig. 9, the initial crack is a single crack part with length 10, but every extended crack part with lengths Dlc ¼ 0:5Dx. In traditional XFEM, the level set function is used to obtain the Heaviside value. The level set function works by using a directional distance between the probe point and the straight crack segment to distinguish a point in which side of a crack, thus this method requires the cracks are straight in every single element. However, the level set function is not suitable for the mesh independent piecewise crack, since there could be multiple crack parts in one element. To resolve this problem, an effective mesh independent horizontal ray method is developed to obtain the Heaviside value, which could be summarized as: (1) for a given probe point, imagine a horizontal ray to the left starting from it; (2) count intersections between the horizontal ray and every crack part segment; (3) increase the count by one if the probe point is above the crack starting point; (4) the Heaviside value of the probe point is 1 if the above-mentioned count number is odd, otherwise is 1. The details and effectiveness of the horizontal ray method are described and discussed in Appendix C.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

339

Fig. 7. A schematic diagram of mesh independent piecewise linear crack. A grid mesh is used, and the crack is divided to multiple straight short parts, which split points are marked with a cross.

3.3. Verification of XFEM numerical results In Section 2.2.1, for a given weak plane direction a in the pure mode I crack type, the critical weakness Gcw =Gc0 to determine a crack deflect or not is obtained. For example, when a ¼ p=6, the critical value

Gcw ¼ 0:870851; Gc0

ð20Þ

is obtained by Eq. (9). Therefore, the weak plane model predicts any material with weakness parameter Gcw =Gc0 < 0:870851 will let the crack to deflect to the weak plane direction; otherwise the crack will extend straightly in K I loading direction. Two numerical examples with Gcw =Gc0 ¼ f0:87; 0:88g separately are chosen as verifications of the weak plane model and the MERR criterion embedded in the XFEM. To compose a pure mode I crack loading case, the following numerical modeling case is set up: a square plane with side length 200, and an initial crack which is horizontally placed in the middle point of left boundary, with length 10. A square grid mesh with 99  99 elements is used. Tensile loadings are applied normally to the upper and lower square boundaries. The two nodes on the middle of right boundary are fixed, i.e. ux ¼ 0 and uy ¼ 0. There is a weak plane in the material, with direction a ¼ p=6. Fig. 8 shows the numerical results for crack path, where the toughnesses of weak planes are set to (a) Gcw =Gc0 ¼ 0:87, and (b) Gcw =Gc0 ¼ 0:88. We find a great match between the numerical results and the critical weakness in Eq. (20), since the crack deflects to the weak plane in Fig. 8a but extends straightly in Fig. 8b. 4. Numerical modeling results 4.1. Numerical results for three point bending tests In this section, three point bending (TPB) test is modeled to show the effect of weak plane model. As in Section 3.3, a square plate with side length 200 is modeled, and an initial horizontal crack with length 10 is horizontally placed in the middle point of the left boundary. Fig. 9 shows the loading details. Two points at a quarter and three quarters of the left side boundary are fixed on displacement on x-direction, or ux ¼ 0, separately. At the middle point of the right boundary, the displacement on y-direction is fixed (uy ¼ 0), and the loading is applied to the x-direction. A square grid mesh with 99  99 elements is used, where every element is also a square with side length about 2:0202. There are 100  100 nodes in total. The diagonal lines hatched in Fig. 9 indicate the weak plane in the test medium. 4.1.1. Crack path affected by the weakness ratio of weak plane Gcw =Gc0 With a given weak plane direction a ¼ p=6, Fig. 10 shows the crack path of XFEM modeling results for five different weak plane weakness, Gcw =Gc0 ¼ f0:51, 0:61, 0:71, 0:81, 0:91g separately. It could be found that, there is a competition between weak plane and mode I loading, and the crack path is controlled by the stronger one. The weak plane is more attractive to the crack path when the weakness ratio Gcw =Gc0 is lower. The limitation obtained in Eq. (20) still works in Fig. 10, and any weak plane with Gcw =Gc0 > 0:870851 for a ¼ p=6 leads to a straight crack, as shown for Gcw =Gc0 ¼ 0:91, since the pure mode I loading state holds for every time step after crack extending. Fig. 10 shows that the weakness Gcw =Gc0 ¼ 0:51 will lead the crack to weak plane, which is a combined effect of Eqs. (15) and (20). Firstly, Eq. (20) shows that the crack will turn to weak plane at the very first step if Gcw =Gc0 < 0:870851 for a pure mode I loading. Secondly, Eq. (15) shows that, once the crack direction turns to weak plane direction, then the weakness Gcw =Gc0 ¼ 0:51 < 0:553617 is strong enough to trap the crack.

340

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

a)

Gcw = 0.87 Gc 0

b)

Gcw = 0.88 Gc 0

Weak plane direction

Crack path Initial crack

Fig. 8. Crack paths obtained by the XFEM numerical modeling in a square plane with a pure mode I loading as a verification, where the weak plane direction a ¼ p=6. The weak plane model theoretically predicts the crack path would deflect to weak plane direction if Gcw =Gc0 < 0:870851, otherwise the crack would extend straightly. Different weakness is used as (a) Gcw =Gc0 ¼ 0:87, and (b) Gcw =Gc0 ¼ 0:88. The crack paths show great matches between XFEM numerical results and theoretically predictions.

50

Initial crack

100

l0 = 10

F

50

200

Fig. 9. A schematic diagram of three point bending test.

For the case Gcw =Gc0 ¼ 0:81 and a ¼ p=6, Fig. 11 shows three contour plots of ryy as a group of examples to view crack path history and stress distribution. An oscillatory crack path could also be observed for this case both in Figs. 10 and 11, which will be discussed in Section 4.3. 4.1.2. Crack path affected by the weak plane direction a The crack paths obtained from XFEM modeling with different weak plane directions a ¼ fp=16, p=8, 3p=16, p=4g, are shown in Fig. 12. The weakness is set to Gcw =Gc0 ¼ 0:8 for the four examples. We could find that, the weak plane is more attractive to crack if its direction is closer to the mode I loading direction. The weak plane of a ¼ p=4 is too far away from the mode I loading direction, and does not capture the crack at all. In fact, as the critical weakness obtained in Eq. (20), the corresponding critical weak plane direction for the weakness ratio Gcw =Gc0 ¼ 0:8 could also be obtained as

a ¼ 0:211144p ¼ 38:0 ;

ð21Þ

by using Eq. (9). Therefore, all the weak plane whose direction angle is above this critical value would not affect the crack deflection, and leads to a straight path. The cracks in other three cases all extend to the weak plane direction firstly, and then leave the initial deflected path. The weak plane a ¼ p=16, which is closer to the mode I loading direction, almost traps the crack path in its direction, but the crack path of a ¼ 3p=16 is just oscillating near the middle line of the plane. This will be further discussed in Section 4.3.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

341

Fig. 10. Crack paths obtained from XFEM modeling for 5 different weakness Gcw =Gc0 separately. The weak plane direction is denoted as a dashed line, where a ¼ p=6. It could be observed that weak plane is more attractive to the crack path with a lower weakness ratio Gcw =Gc0 .

Fig. 11. Contour plots of ryy , which are loaded by three point bending, show the crack extending history. The weak plane properties are set to Gcw =Gc0 ¼ 0:81 and a ¼ p=6.

We also try to study the effects of initial crack direction, which is denoted by the angle between initial crack direction and x-axis hinit . However, it is found that hinit only affects the initial part of the crack path. 4.2. Modeling of crack extending in a multilayered medium To show the difference of crack behavior between the isotropic material and the weak plane material, a numerical example for a crack penetrating multiple layers of materials is modeled. The loading state holds as TPB, which is the same as in Section 4.1, but the square plane is divided horizontally into five parts equally, see Fig. 13, where the first, third, and fifth parts are weak plane material with properties a ¼ p=6 and Gcw =Gc0 ¼ 0:82, and the second and fourth parts are isotropic material. The crack path shows that the crack is willing to deflect to weak plane direction in every weak plane material part, but turns back to the horizontal middle line, i.e. the mode I loading direction. We can call this ‘crack diffraction’ at the material interface like the traveling light. The crack shape looks much rougher in the weak plane parts than in the isotropic parts. In the last part of weak plane material, the crack extends straightly to the weak direction, which might be a result of the compress stress loaded at the middle point at the right boundary. 4.3. Discussion of the oscillation cracks Periodically oscillatory crack paths are observed in Sections 4.1.1 and 4.1.2. For example, the configuration set {Gcw =Gc0 ¼ 0:81; a ¼ p=6} in Fig. 10, and the configuration set {Gcw =Gc0 ¼ 0:8; a ¼ 3p=16} in Fig. 12. In fact, the oscillatory

342

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

Fig. 12. Crack paths obtained from XFEM modeling for 4 different a, where Gcw =Gc0 ¼ 0:8. The weak plane is more attractive to crack if its direction is closer to the mode I loading direction. The example of a ¼ p=4, which is above the critical angle in Eq. (21), shows that weak plane direction too far away could not change crack direction.

Fig. 13. Crack path obtained by XFEM modeling. The crack, which is driven by the TPB test loading, starts from left boundary and extends through a 5layered medium, which consists weak plane material and isotropic material in turn. The weak plane material’s properties are set to a ¼ p=6 and Gcw =Gc0 ¼ 0:82.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

343

Fig. 14. The oscillation cracks have been found in many other experiments with different materials, such as (a) a Si3N4 film deposited on a silicon substrate [11], (b) a Bigwood granite rock [1], and (c) an African granodiorite rock [2].

crack path could often be reproduced with material properties in the range Gcw =Gc0 2 ½0:7; 0:9 and a 2 ½p=8; 3p=16. The oscillation cracks have been found in many experiments from past researches, see Fig. 14. In Theorem 2, we have already proved that the crack paths are only related to the angle between crack and weak plane directions a, the weakness Gcw =Gc0 , and the SIFs ratio K II =K I . Since the weakness Gcw =Gc0 is a certain value for one simulation, only the variations of a and K II =K I could lead to the oscillatory cracks. The three point bending loading boundaries are set to create a mode I loading state, but it is only effective in the area near the middle horizontal line of the square. When the crack path leaves the middle line, the TPB loadings will try to bend the crack path back. On the other hand, the weak plane direction is attractive to the crack path because of inherent property. For a proper weak plane material property pair fa; Gcw =Gc0 g, the crack path would turn to weak plane direction firstly, but the bending force of TPB loadings would quickly increase when the crack extends close to the boundary of effective TPB area. This competition between weak plane and TPB loading might be a main reason for the periodically oscillatory crack path. This oscillatory crack path phenomenon might be an inherent crack pattern of weak plane material. Details of the phenomenon, such as relations of period and amplitude with the material properties and loading state, still need future study.

5. Conclusion A weak plane model is demonstrated in this article, as an example to show the effects of anisotropic fracture toughness on crack extending. The MERR criterion is used to obtain the crack deflecting direction in an anisotropic fracture toughness material. It is proved that the crack deflecting direction in a 2D problem only relates to three values in the weak plane model: (1) the weak plane direction a, (2) the weakness ratio of the weak plane Gcw =Gc0 , (3) the stress intensity factor ratio K II =K I at the crack tip induced by loading. And the crack deflecting direction for any loading states and material constants could be theoretically obtained. The crack deflecting forbidden area phenomenon is discussed. Two critical weakness ratios of material are found, which are, separately, the critical weakness ratio that weak plane captures all cracks for any loading state, and the critical weakness ratio that the crack is trapped by the weak plane once it deflects to the plane. A numerical method is developed to simulate the crack propagation in the anisotropic material, by embedding the MERR criterion and the weak plane model into XFEM. It is verified by a pure mode I loading state, which finds a pretty good match between numerical results and theoretical prediction. By changing the material properties Gcw =Gc0 and a, groups of crack paths have been compared, where a three point bending test loading state is applied. Periodically oscillatory crack path behaviors are found in the numerical modeling, which is attributed to the competition between weak plane and loading state. Further works are needed to find out the relations between the period and amplitude of the oscillation, and the material properties. An example of crack extending through multiple layers of weak plane material and isotropic material is studied, which shows that the crack diffracts at the material interface, and the crack path in the anisotropic toughness material is rougher than in the isotropic material. The weak plane model is useful for crack behavior study in rocks, crystals, and some kinds of composite materials. This model is easy to be measured and verified in laboratory, and is convenient to be used in the numerical simulation. The presented work in this paper will be helpful to understand the fracture laws in the rocks as well as other anisotropic materials. The crack propagation path could be controlled by designing the anisotropic properties in the medium. In addition, the material properties could be inferred from the crack propagation path.

344

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

Acknowledgements This work is supported by the National Natural Science Foundation of China under Grant No. 11532008, Research Institute of Petroleum Exploration and Development of China National Petroleum Corporation, Drilling Research Institute of China National Petroleum Corporation, Tsinghua University Initiative Scientific Research Program, Chinese 1000-talents Plan for Young Researchers, and Science Challenge Program No. JCKY2016212A502. Appendix A. SIFs after crack deflection ~ I ðhÞ; K ~ II ðhÞg and This appendix shows the results from Amestoy and Leblond [32], which gives the relations between fK fK I ; K II g for a crack extended in the direction of h, with a enough small length. Amestoy and Leblond proved that the SIFs after crack deflecting with direction h ¼ m p are linear combinations of SIFs before crack extending,

~ I ðmpÞ ¼ F 11 ðmÞK I þ F 12 ðmÞK II ; K ~ II ðmpÞ ¼ F 21 ðmÞK I þ F 22 ðmÞK II ; K

ðA:1Þ

and coefficients in Eq. (A.1) expanded up to 20th order series are obtained,

F 11 ðmÞ ¼ 1 

3m2 p2 þ 8



p2 

  2  5p4 p 11p4 119p6 6 m4 þ m þ 5:0779m8  2:88312m10  0:0925m12  þ 128 9 72 15360

þ 2:996m14  4:059m16 þ 1:63m18 þ 4:1m20 ; F 12 ðmÞ ¼ 

    3mp 10p p3 133p3 59p5 m3 þ 2p  m5 þ 12:3139m7  7:32433m9 þ 1:5793m11 þ þ þ 2 3 16 180 1280

þ 4:0216m13  6:915m15 þ 4:21m17 þ 4:56m19 ; F 21 ðmÞ ¼

ðA:2aÞ

ðA:2bÞ

    mp 4p p 3 2p 13p3 59p5  þ m3 þ  þ  m5  6:17602m7 þ 4:44112m9  1:534m11 2 3 48 3 30 3840  2:07m13 þ 4:684m15  3:95m17  1:32m19 ;

ðA:2cÞ

      3p 2 8 29p2 5p4 32 4p2 1159p4 119p6 þ  F 22 ðmÞ ¼ 1  4 þ m2 þ  m4 þ   þ m6 þ 10:5825m8 3 15 8 18 128 9 7200 15360  4:78511m10  1:8804m12 þ 7:28m14  7:591m16 þ 0:25m18 þ 26:6681m20 :

ðA:2dÞ

It is noted that F 11 ðmÞ and F 22 ðmÞ are even functions; F 12 ðmÞ and F 21 ðmÞ are odd functions. Appendix B. A comparison of kinked SIFs results among Nuismer, He and Hutchinson, and Amestoy and Leblond Nuismer [29], He and Hutchinson [30,31] and Amestoy and Leblond [32] all found that the relations of crack SIFs after ~ I ðhÞ; K ~ II ðhÞg and before deflection fK I ; K II g are linear as shown in Eq. (A.1). deflection fK ~ II ðhÞg and the SIFs before deflection fK I ; K II g in the ~ I ðhÞ; K Nuismer [29] found a simple relation between the kinked SIFs fK isotropic material,

~ I ðhÞ ¼ 1 cos h ½K I ð1 þ cos hÞ  3K II sin h; K 2 2 ~ II ðhÞ ¼ 1 cos h ½K I sin h þ K II ð3 cos h  1Þ; K 2 2

ðB:1Þ

which is also used by Zeng and Wei [16]. He and Hutchinson [30] faced this problem when they studied the interface crack kink behavior on the boundary of two different elastic media. They gave the results in a table when the material constants in the two media are the same [31, p. 2– 5] for a deflection angle 6 6 h 6 135 , which could be used in this comparison. The comparison results are shown in Fig. B.1. It could be found that the results between He and Hutchinson [30] and Amestoy and Leblond [32] match perfectly, particularly in the range h < 90 , which is needed for this article. But the results of Nuismer [29], i.e. Eq. (B.1), deviate from the others, especially in F 12 and F 22 . For instance, when h ¼ p=6, the value of F 22 from Nuismer [29] is about 10:6% smaller than Amestoy and Leblond [32]. Therefore, we think Nuismer [29]’s results Eq. (B.1) are not accurate enough, and are not appropriate to be used for crack deflection problem if the deflection angle h > p=6, or the problem whose crack direction could deviate from the weak plane by p=6 or more, when K II could not be ignored before the crack deflection.

345

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

b) 0.0

0.8

0.2 0.4

0.6

F 12

F 11

a) 1.0

0.4

0.8

0.2

1.0 1.2

0.0 0

0.4

0.6

8

4

3 8

2

5 8

3 4

0

m

c)

8

4

3 8

2

5 8

3 4

2

5 8

3 4

m

d) 1.0

0.3

F 22

F 21

0.5 0.2

0.0

0.1 0.0

0.5 0

8

4

3 8

2

5 8

3 4

0

8

3 8

m

m Amestoy & Leblond (1992)

4

He & Hutchinson (1989)

Nuismer (1975)

Fig. B.1. A comparison among the kink SIFs ratio results obtained separately by Nuismer [29], He and Hutchinson [30], and Amestoy and Leblond [32]. The relations between SIFs before and after a crack deflection are shown in Eq. (A.1), which is controlled by four ratios, and these ratios varying with kink angle h ¼ mp are shown in this figure in the order of (a) F 11 , (b) F 12 , (c) F 21 , (d) F 22 , separately.

Therefore, we find that the results obtained by Amestoy and Leblond [32] are accurate enough and easy for using, and rest the article uses their formulations Eq. (4), (A.1) and (A.2) only. Appendix C. Horizontal ray method To determine whether a point is inside a polygon, a classical algorithm is introduced by Shimrat [35], which uses a ray, usually horizontal, starting from the probe point, and ending to infinity in any direction. It could be found that, the probe point locates inside the polygon if and only if there are odd intersections between the ray and polygon boundaries. This algorithm is easy to be extended to use in XFEM with piecewise linear cracks, to determine the relative position between a probe point and the crack path, and obtain the Heaviside value of the point. The original version of the horizontal ray method could only determine whether the point is in a polygon or not. But in the crack Heaviside problem, although the crack path is composed of linear segments, the domain boundaries might be curved, and the shape built up with crack path and domain boundaries is usually not closed. Therefore, a modified algorithm is developed by introducing a imaginary close shape with three additional virtual boundaries. Fig. C.1 shows an example of the modified horizontal ray method. Consider a domain X containing a crack from the red solid point to the green point, with the blue crack path. Thus, every point above the crack path should has a Heaviside value 1, otherwise 1. Therefore, The Heaviside function HðxÞ in Eq. (16) could be formulated as


HðxÞ ¼

1;

point x above the crack path;

1; otherwise:

ðC:1Þ

This criterion could be converted to that, the Heaviside value should be 1 in the gray domain in Fig. C.1, otherwise 1. The modified method uses a horizontal ray to the left, which is shown as a dash-dotted line for every probe point {A, B, C} in Fig. C.1, with the following three virtual boundaries: (1) The green horizontal dashed line is a virtual boundary of the Heaviside value, since the Heaviside values of the points ahead the crack path does not affect the results of XFEM. The green dashed line would never be intersected, since it is horizontal.

346

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

A

B C

Ω

Fig. C.1. A schematic diagram of the mesh independent horizontal ray method to obtain the Heaviside value of a probe point. Count the number of intersections nI between the horizontal ray from the probe point and the crack path as well as the red dashed line on the left, which starts from the crack starting point. The probe point is on upper side of the crack if nI is odd, otherwise ni is even. It could be counted that nI ðAÞ ¼ 1; nI ðBÞ ¼ 3, and nI ðCÞ ¼ 4.

(2) The red vertical dashed line on the left indicates the left limitation of the gray domain, which should be a vertical line placed on the left of any points in the domain, and the starting dashed circle point is a horizontal translation of the red solid point. The intersection counting rule for the red dashed line could be simplified as, counting number nI increases by one if the probe point locates above the crack starting point. (3) The purple dashed curve out of the domain X indicates a virtual boundary for the gray domain. It should be located on the above and right of the domain X, and any probe point in X should not gain intersections amount nI from the purple curve line. This horizontal ray method to obtain Heaviside value is also totally independent to the mesh. All we need to obtain Heaviside value for a given probe point, are the coordinates of crack parts vertexes and crack starting point. And the computational complexity does not increase when mesh is refined. This algorithm could be further optimized, e.g., there would be no intersections if the two vertexes of the crack part segment are both above or below the probe point, which could be determined by the y-coordinate of these three points.

References [1] Nasseri M, Mohanty B. Fracture toughness anisotropy in granitic rocks. Int J Rock Mech Min Sci 2008;45(2):167–93. doi: http://dx.doi.org/10.1016/j. ijrmms.2007.04.005. [2] Kataoka M, Obara Y, Kuruppu M. Estimation of fracture toughness of anisotropic rocks by semi-circular bend (SCB) tests under water vapor pressure. Rock Mech Rock Eng 2015;48(4):1353–67. doi: http://dx.doi.org/10.1007/s00603-014-0665-y. [3] Chandler MR, Meredith PG, Brantut N, Crawford BR. Fracture toughness anisotropy in shale. J Geophys Res: Solid Earth 2016;121(3):1706–29. doi: http://dx.doi.org/10.1002/2015JB012756. [4] Ast J, Przybilla T, Maier V, Durst K, Göken M. Microcantilever bending experiments in NiAl evaluation, size effects, and crack tip plasticity. J Mater Res 2014;29(18):2129–40. doi: http://dx.doi.org/10.1557/jmr.2014.240. [5] Wang Z, Shen X, Akbari Garakani M, Lin X, Wu Y, Liu X, Sun X, Kim J-K. Graphene aerogel/epoxy composites with exceptional anisotropic structure and properties. ACS Appl Mater Interf 2015;7(9):5538–49. doi: http://dx.doi.org/10.1021/acsami.5b00146. [6] Gómez-Monterde J, Schulte M, Ilijevic S, Hain J, Sánchez-Soto M, Santana O.O., et al. Effect of microcellular foaming on the fracture behavior of ABS polymer. J Appl Polym Sci 133 (7). http://dx.doi.org/10.1002/app.43010. [7] Andersons J, Cbulis U, Stiebra L, Kirpluks M, Sprniš E. Modeling the mode I fracture toughness of anisotropic low-density rigid PUR and PIR foams. Int J Fract 2017;205(1):111–8. doi: http://dx.doi.org/10.1007/s10704-017-0194-2. [8] Takei A, Roman B, Bico J, Hamm E, Melo F. Forbidden directions for the fracture of thin anisotropic sheets: an analogy with the Wulff plot. Phys Rev Lett 2013;110(14):1–5. doi: http://dx.doi.org/10.1103/PhysRevLett.110.144301. [9] Audoly B, Reis PM, Roman B. Cracks in thin sheets: when geometry rules the fracture path. Phys Rev Lett 2005;95(2):025502. doi: http://dx.doi.org/ 10.1103/PhysRevLett.95.025502. [10] Reis PM, Audoly B, Roman B. Cracking sheets: oscillatory fracture paths in thin elastic sheets. Chaos 2008;18(4):041108. doi: http://dx.doi.org/ 10.1063/1.2997335. [11] Nam KH, Park IH, Ko SH. Patterning by controlled cracking. Nature 2012;485(7397):221–4. doi: http://dx.doi.org/10.1038/nature11002. [12] Hakim V, Karma A. Crack path prediction in anisotropic brittle materials. Phys Rev Lett 2005;95(23):235501. doi: http://dx.doi.org/10.1103/ PhysRevLett.95.235501. [13] Gao H, Chiu C-H. Slightly curved or kinked cracks in anisotropic elastic solids. Int J Solids Struct 1992;29(8):947–72. doi: http://dx.doi.org/10.1016/ 0020-7683(92)90068-5. [14] Hakim V, Karma A. Laws of crack motion and phase-field models of fracture. J Mech Phys Solids 2009;57(2):342–68. doi: http://dx.doi.org/10.1016/j. jmps.2008.10.012. [15] Chambolle A, Francfort G, Marigo J-J. When and how do cracks propagate? J Mech Phys Solids 2009;57(9):1614–22. doi: http://dx.doi.org/10.1016/j. jmps.2009.05.009.

Y. Gao et al. / Engineering Fracture Mechanics 180 (2017) 330–347

347

[16] Zeng X, Wei Y. Crack deflection in brittle media with heterogeneous interfaces and its application in shale fracking. J Mech Phys Solids 2017;101:235–49. doi: http://dx.doi.org/10.1016/j.jmps.2016.12.012. [17] He M-Y, Hutchinson JW. Crack deflection at an interface between dissimilar elastic materials. Int J Solids Struct 1989;25(9):1053–67. doi: http://dx.doi. org/10.1016/0020-7683(89)90021-8. [18] Hutchinson J, Suo Z. Mixed mode cracking in layered materials. In: Advances in applied mechanics, vol. 29; 1991. p. 63–191. http://dx.doi.org/10.1016/ S0065-2156(08)70164-9. [19] Li B, Peco C, Millán D, Arias I, Arroyo M. Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy. Int J Numer Methods Eng 2015;102(3–4):711–27. doi: http://dx.doi.org/10.1002/nme.4726. [20] Li C, Caner FC, Chau VT, Bazˇant ZP. Spherocylindrical microplane constitutive model for shale and other anisotropic rocks. J Mech Phys Solids 2017;103:155–78. doi: http://dx.doi.org/10.1016/j.jmps.2017.03.006. [21] Shanthraj P, Svendsen B, Sharma L, Roters F, Raabe D. Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture. J Mech Phys Solids 2017;99(November 2016):19–34. doi: http://dx.doi.org/10.1016/j.jmps.2016.10.012. [22] Li W, Rezakhani R, Jin C, Zhou X, Cusatis G. A multiscale framework for the simulation of the anisotropic mechanical behavior of shale. URL . [23] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46(1):131–50. doi:10.1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J.. [24] Liu ZL, Menouillard T, Belytschko T. An XFEM/spectral element method for dynamic crack propagation. Int J Fract 2011;169(2):183–98. doi: http://dx. doi.org/10.1007/s10704-011-9593-y. [25] Gordeliy E, Peirce A. Coupling schemes for modeling hydraulic fracture propagation using the XFEM. Comput Methods Appl Mech Eng 2013;253:305–22. doi: http://dx.doi.org/10.1016/j.cma.2012.08.017. [26] Zhao J, Bessa MA, Oswald J, Liu Z, Belytschko T. A method for modeling the transition of weak discontinuities to strong discontinuities: from interfaces to cracks. Int J Numer Methods Eng 2016;105(11):834–54. doi: http://dx.doi.org/10.1002/nme.4995. [27] Zeng Q, Liu Z, Wang T, Gao Y, Zhuang Z. Fully coupled simulation of multiple hydraulic fractures to propagate simultaneously from a perforated horizontal wellbore. Comput Mech. http://dx.doi.org/10.1007/s00466-017-1412-5. [28] Yau JF, Wang SS, Corten HT. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. J Appl Mech 1980;47(2):335. doi: http://dx.doi.org/10.1115/1.3153665. [29] Nuismer RJ. An energy release rate criterion for mixed mode fracture. Int J Fract 1975;11(2):245–50. doi: http://dx.doi.org/10.1007/BF00038891. [30] He M-Y, Hutchinson JW. Kinking of a crack out of an interface. J Appl Mech 1989;56(2):270. doi: http://dx.doi.org/10.1115/1.3176078. [31] He M-Y, Hutchinson JW. Kinking of a crack out of an interface: tabulated solution coefficients. Harvard University Report MECH-113A. URL . [32] Amestoy M, Leblond J. Crack paths in plane situations – II. Detailed form of the expansion of the stress intensity factors. Int J Solids Struct 1992;29 (4):465–501. doi: http://dx.doi.org/10.1016/0020-7683(92)90210-K. [33] Hoek E. Fracture of anisotropic rock. J S Afr Inst Min Metall 1964;64(10):501–18. [34] Bao J, Fathi E, Ameri S. Uniform investigation of hydraulic fracturing propagation regimes in the plane strain model. Int J Numer Anal Methods Geomech 2015;39(5):507–23. doi: http://dx.doi.org/10.1002/nag.2320. [35] Shimrat M. Algorithm 112: position of point relative to polygon. Commun ACM 1962;5(8):434. doi: http://dx.doi.org/10.1145/368637.368653.