Computation of T-stresses for multiple-branched and intersecting cracks with the numerical manifold method

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Computation of T-stresses for multiple-branched and intersecting cracks with the numerical manifold method H.H. Zhang a, S.M. Liu a, S.Y. Han a, L.F. Fan b,∗ a b

School of Civil Engineering and Architecture, Nanchang Hangkong University, Nanchang, Jiangxi 330063, China College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100084, China

a r t i c l e

i n f o

Keywords: T-stress Multiple-branched and intersecting cracks Numerical manifold method Interaction integral

a b s t r a c t Due to its dual cover systems (i.e., the mathematical cover and the physical cover), the numerical manifold method (NMM) is outstanding in crack modeling, especially in the simulation of complex cracks. In this paper, the NMM is extended to calculate the T-stresses for two-dimensional multiple-branched and intersecting cracks. The displacement jump across crack surface is essentially portrayed by the NMM, and the stress singularity at the crack tip is represented through the incorporation of related asymptotic basis in the NMM local functions. The T-stresses are extracted with the domain-form interaction integral in the NMM postprocessing. To validate the proposed scheme, three numerical examples are tested on uniform mathematical covers nonconforming with all domain boundaries, and the computed results are in nice consistence with the existing solutions. Besides, the effects of crack configurations and boundary conditions on the T-stresses are also revealed.

1. Introduction In the classical linear elastic fracture mechanics, the stress intensity factors (SIFs) play a dominant role in the evaluation of the crack tip fields and were even supposed to be adequate to characterize the associated fracture behaviors [1]. However, lots of studies (e.g., [2–6]) have shown that the influences of another fracture parameter, i.e., the T-stress (the first non-singular term in the William’s expansion of the normal stress parallel to the crack and ahead of a crack tip [1,7]) cannot be overlooked. The T-stress was found to affect the size and shape of the plastic zone (e.g., see [2,3,8]), crack path stability (e.g., see [4,9,10]), the crack tip constraint (e.g., see [11–13]) and also the fracture toughness (e.g., see [14–16]). Accordingly, much work declares that the crack tip behaviors are governed by both SIFs and T-stress [17–19]. However, compared with the intensively studied SIFs, the T-stress is far less understood. Therefore, to better represent the fracture features, it is of significant importance to accurately compute the T-stress. Among various tools for the calculation of T-stress (e.g., see [7]), numerical methods such as the finite element method (FEM), the boundary element method (BEM), the extended finite element method (XFEM) and the scaled boundary finite element method (SBFEM) are more attractive because of their wider applicability. Chen et al. [20] evaluated the Tstress of mode-I crack using the hierarchical p-version FEM and two path-independent integrals. By virtue of the interaction integral and the FEM, Kim and Paulino [21] obtained the T-stresses, mixed mode SIFs

∗

and crack initiation angles in functionally graded materials (FGMs). Jogdand and Murthy [22] combined an interior collocation method with the FEM to analyze the mixed-mode SIFs and T-stress. Bouchard et al. [23] studied the effect of SIFs and T-stress on the crack propagation in the kerf-less spalling process of single crystal silicon foils with the FEM. Yu et al. [24] reported the SIFs and elastic T-stress results for through-wall-cracked pipes under internal pressures by 3D FEM. Sladek and Sladek [25] investigated the T-stress of interface cracks between dissimilar materials through the BEM and path-independent integral. Tan and Wang [26] incorporated quarter-point crack-tip elements into the BEM to compute the T-stress for mode-I cracks. Guduru et al. [27] took the symmetrical Galerkin BEM to assess the dynamic SIFs and T-stress for cracked particulate composite materials. Chen et al. [28] summarized the SIFs and T-stress for multi-crack problems using a spline fictitious boundary element alternating method. Yu et al. examined the T-stress of FGMs under mechanical loading [29] and thermal loading [30] through the XFEM and interaction integral. Song and his coauthors presented the T-stress of in-plane cracks under static loading [31] and dynamic loading [32] with the SBFEM; besides, they further coupled the XFEM and SBFEM to measure the SIFs and T-stress for interfacial cracks and cracks terminating at the interface [33]. In recent decade, the numerical manifold method (NMM) gradually comes into our sights due to its predominant advantages in discontinuous analysis, especially in crack modeling, benefiting from its two cover systems (i.e., the mathematical cover (MC) and the physical cover (PC))

Corresponding author. E-mail address: [email protected] (L.F. Fan).

https://doi.org/10.1016/j.enganabound.2019.07.011 Received 11 April 2019; Received in revised form 11 July 2019; Accepted 11 July 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

) 𝐾 ( 𝐾 𝜃 3𝜃 𝜃 𝜃 𝜃 3𝜃 𝜎22 = √ I 1 + sin sin + 𝑂(𝑟1∕2 ) cos + √ II sin cos cos 2 2 2 2 2 2 2𝜋𝑟 2𝜋𝑟 (2) ) 𝐾 𝐾 ( 𝜃 3𝜃 𝜃 𝜃 3𝜃 𝜃 𝜎12 = √ I sin cos cos + √ II 1 − sin sin cos + 𝑂(𝑟1∕2 ) 2 2 2 2 2 2 2𝜋𝑟 2𝜋𝑟 (3) where 𝜎 ij (i, j = 1, 2) represent the stress components defined in the crack tip coordinate system x1 ox2 plotted in Fig. 1, and the subscripts 1 and 2 refer to, respectively, x1 and x2 axis, with x1 axis parallel and x2 axis perpendicular to the crack surface. (r, 𝜃) are the polar coordinates as indicated in Fig. 1. KI and KII are, respectively, the mode-I and mode-II SIF. The r-independent term T in Eq. (1) has dimension of stress and is frequently named T-stress [3]. The corresponding displacement fields are expressed as [59] √ [ ] 𝐾 𝑟 𝜃 3𝜃 𝑢1 = I (2𝜅 − 1) cos − cos 4𝜇 2𝜋 2 2 √ [ ] 𝐾 𝑟 𝜃 3𝜃 + II (2𝜅 + 3) sin + sin 4𝜇 2𝜋 2 2

Fig. 1. A linear elastic isotropic body with a traction-free crack and coordinate systems.

[34]. The superiorities of the NMM for fracture simulation are primarily reflected in: (a) the MC (similar to the mesh in the FEM) can be inconsistent with all domain boundaries including crack surfaces; (b) the discontinuity of physical fields (e.g., displacements) can be naturally captured without any special treatments, which is very suitable for complex problems containing multi-branched and intersecting cracks; (c) the singularity of field gradients (e.g., stresses) at the crack tip can be accurately manifested through the appropriate definition of local functions on associated physical patches (the subset of the PC). Since its advent, the NMM has been extensively developed to study a variety of crack problems. Chiou et al. [35] utilized the NMM and the virtual crack extension method to study the mixed-mode fracture propagation. Li et al. [36] linked the meshless concept, the NMM cover systems and crack-tip enrichment idea to simulate 2D crack problems. Terada et al. [37] traced the progressive and cohesive failure of heterogeneous solids and structures with the NMM. Wu et al. considered the frictional crack initiation and propagation process [38], the effect of water flow in fractures on the stability of rock structures [39] and the cohesive fracture of rocks [40] using the NMM. Zheng and his coauthors [41–43] established a moving least square NMM for static and dynamic crack problems. Li et al. [44] proposed a particle-based NMM to investigate rock fracturing behavior under dynamic loading. Yang et al. [45] modeled the fluid driven fracturing in rocks by an enriched NMM. Yang et al. [46] employed the NMM and a non-local crack tracking method to capture the 3D crack growth paths. Zhang and his coauthors focused on 2D crack problems in homogeneous and nonhomogeneous materials under mechanical or thermo-mechanical loadings [47–55]. However, to date, the report of T-stress results in the framework of NMM is still an open issue, although some related schemes were proposed [56,57]. In this paper, seeing the importance of T-stress and the advantages of the NMM for complex crack simulation, we further extend the NMM to compute the T-stresses for multi-branched and intersecting cracks. To this end, the remaining contents are arranged as follows. In Section 2, 2D crack tip asymptotic fields regarding the contribution of T-stress term are presented. Following, some basic concepts and equations of the NMM for crack modeling, and the extract scheme for T-stress are provided in Section 3. Next, detailed computer implementations on the calculation of T-stress are discussed in Section 4. Then, in Section 5, three numerical examples are tested. Finally, the corresponding concluding remarks are drawn in Section 6.

+

𝑢2 =

𝑇 𝑟(𝜅 + 1) cos 𝜃 + 𝑂(𝑟3∕2 ) 8𝜇

(4)

√

] [ 𝑟 𝜃 3𝜃 (2𝜅 + 1) sin − sin 2𝜋 2 2 √ [ ] 𝐾 𝑟 𝜃 3𝜃 − II (2𝜅 − 3) cos + cos 4𝜇 2𝜋 2 2 𝐾I 4𝜇

+

𝑇 𝑟(𝜅 − 3) sin 𝜃 + 𝑂(𝑟3∕2 ) 8𝜇

(5)

where 𝜇 = 0.5𝐸∕(1 + 𝜈) is the shear modulus, with E the Young’s modulus and 𝜈 the Poisson’s ratio. 𝜅 equals (3 − 𝜈)∕(1 + 𝜈) in plane stress problem and (3 − 4𝜈) for plane strain case. 3. The NMM for the computation of T-stress 3.1. Key terminology of the NMM The NMM, originally proposed for the efficient and accurate deformation analysis of rock mass containing a series of joints and cracks [34], has been proved to be a very powerful tool for fracture modeling (e.g., see [47,48]). Compared with other numerical methods (e.g., the FEM and BEM), the leading highlight of the NMM is its dual cover systems (i.e., the MC and PC), associated with which another three basic concepts, namely, the mathematical patch (MP), the physical patch (PP) and the manifold element (ME) are also unique. Herein, an illustrated example will be used to clarify the above mentioned terminology. Consider a rectangular plate with a multi-branched crack shown in Fig. 2(a). As a primary step to model this problem with the NMM, domain discretization is firstly carried out. For this purpose, a large enough MC, that is, the union of all overlap-permitted mathematical patches (MPs), should be constructed to cover the whole physical domain. In theory, the MC is allowed to be nonconforming with all domain boundaries, and the shape and number of MP can be arbitrary (note that the MP geometry and amount may affect the solution accuracy and efficiency). Hereby, an MC in Fig. 2(b), formed by MPs consisting of square mathematical elements (e.g., the MP MI in Fig. 2(b) is made of four squares, and each square is just a mathematical element), and consistent with neither external boundary nor internal crack surfaces, is built to cover the rectangular domain (see Fig. 2(c) for the interaction of the MC and the physical domain). Following, through the intersection of the MPs and the physical domain, the physical patches (PPs) can be obtained. Generally, if an MP is divided by the domain boundaries into

2. Crack tip asymptotic fields Consider a traction-free crack with sharp crack tips in a planar linear elastic isotropic domain as shown in Fig. 1. The associated near-tip stresses can be written as [58] ) ) 𝐾 ( 𝐾 ( 𝜃 3𝜃 𝜃 𝜃 3𝜃 𝜃 cos − √ II 2 + cos cos sin 𝜎11 = √ I 1 − sin sin 2 2 2 2 2 2 2𝜋𝑟 2𝜋𝑟 ( 1∕2 ) +𝑇 +𝑂 𝑟 (1) 150

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Fig. 2. Illustration of some basic concepts in the NMM: (a) physical domain, (b) MP and MC, (c) interaction of physical domain and the MC, (d) typical PPs and (e) typical MEs.

n isolated parts, then n PPs are generated. For example, regarding the MP MA (representative of the MPs intersecting with none of the boundaries), MB (representative of the MPs only intersecting with external boundary) and MD (representative of the MPs partially separated by the crack surface) in Fig. 2(c), each gives 1 PP, namely PA , PB and PD in Fig. 2(d); as for MC in Fig. 2(c), since it is completely cut by the branched crack into three segments, 3 PPs are made, i.e., PC1 , PC2 and PC3 in Fig. 2(d). The set of all PPs then forms the PC. Subsequently, each common zone of as many as possible PPs gives a manifold element (ME); accordingly, in Fig. 2(e), the overlapped area of 2 PPs from MA -type MP and 2 from MB -type generates the ME Ea , 4 from MB -type makes Eb , 2

from MC -type and 2 from MD -type establishes Ec , and 4 from MD -type determines Ed . 3.2. Crack modeling by the NMM After the definition of the core ingredients of the NMM in Section 3.1, the displacement approximation on a certain ME E can be generally expressed as [48] 𝐮ℎ (𝐗) =

𝑁 ∑ 𝐽 =1

151

𝑤𝐽 (𝐗)𝐮𝐽 (𝐗)

(6)

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

where N denotes the number of PPs forming the ME E (e.g., for the marked MEs in Fig. 2(e), N = 4). X = (X1 , X2 ) is a point in this ME, with X1 and X2 the global coordinates as shown in Fig. 1. wJ (X) and uJ (X) are, respectively, the partition of unity weight function and the local function defined on the Jth PP. wJ (X) relates to the shape of MP containing the Jth PP, and in practice frequently takes the form of the shape function for finite elements with similar configuration to the mathematical element [47,48]. uJ (X) is associated with the local physical property in the Jth PP. For PPs containing no crack tip (e.g., PPs from MA -type, MB -type and MC -type MPs in Fig. 2(d)), continuous functions are often adopted, i.e., 𝐮𝐽 (𝐗) = 𝐏(𝐗)𝐚 𝐽 (𝐗)

(7)

where aJ (X) is the vector of conventional unknowns defined on the Jth PP, and P(X) is the polynomial basis being [ ] 1 0 𝑋1 0 𝑋2 0 ⋅⋅⋅ 𝐏(𝐗) = (8) 0 𝑋1 0 𝑋2 ⋅⋅⋅ 0 1

Fig. 3. An illustration of the interaction integral domain.

In Eq. (12), 𝜎𝑖𝑗aux and 𝑢aux represent, respectively, the stresses and 𝑖 displacements in the auxiliary state and they are usually taken as [61]

For PPs containing the crack tip (e.g., PPs from MD -type MPs in Fig. 2(d)), to better capture the near-tip properties, local functions involving crack tip asymptotic basis deduced from Eqs. (4) and (5) are used. For the present study, they are formulated as [48] 𝐮𝐽 (𝐗) = 𝐏(𝐗)𝐚𝐽 (𝐗) + 𝚽(𝐗)𝐛𝐽 (𝐗)

(9)

where bJ (X) is the vector of additional unknowns also defined on the Jth PP. 𝚽(X) is the asymptotic basis being [ ] Φ1 0 Φ2 0 Φ3 0 Φ4 0 𝚽(𝐗) = (10) 0 Φ1 0 Φ2 0 Φ3 0 Φ4 where ] [ ] [√ 𝜃 √ 𝜃 √ 𝜃 √ 𝜃 𝑟 sin , 𝑟 cos , 𝑟 sin 𝜃 sin , 𝑟 sin 𝜃 cos (11) Φ1 , Φ2 , Φ3 , Φ4 = 2 2 2 2 Therefore, through Eq. (9), the crack tip singularity can be well manifested. As for the displacement discontinuity across the crack surface, since the NMM is able to represent the jump of physical field naturally, no additional operations are needed (e.g., the use of extra functions such as the Heaviside functions or jump functions as in the XFEM is avoided), which means that arbitrarily branched and intersecting cracks can be modeled in the same manner as a single crack [47,48]. Upon the displacement approximation in Eq. (6), the corresponding governing equation and boundary conditions (see [48]), the NMM discrete formulations for crack simulation can be derived by, e.g., variational principle, and the associated results can be found in our previous papers (e.g., see [47,48]). Herein, for simplicity, we will not go into more details.

aux 𝜎11 =−

𝐹 cos3 𝜃 𝜋𝑟

(13)

aux 𝜎22 =−

𝐹 cos 𝜃sin2 𝜃 𝜋𝑟

(14)

aux 𝜎12 =−

𝐹 cos2 𝜃 sin 𝜃 𝜋𝑟

(15)

𝑢aux =− 1

𝐹 (𝜅 + 1) 𝑟 𝐹 ln − sin2 𝜃 8𝜋𝜇 𝑑 4𝜋𝜇

(16)

𝑢aux =− 2

𝐹 (𝜅 − 1) 𝐹 𝜃+ sin 𝜃 cos 𝜃 8𝜋𝜇 4𝜋𝜇

(17)

where F is a point force acting at the crack tip in the x1 direction and d is the coordinate of a fixed point on the x1 axis. The interaction integral I in Eq. (12) is related to the T-stress by [61] 𝑇 =

𝐸0 𝐼 𝐹

(18)

where E0 = E for plane stress problem and 𝐸0 = 𝐸∕(1 − 𝜈 2 ) for plane strain case. 4. Numerical implementations for interaction integral

3.3. The interaction integral for the T-stress

4.1. Derivatives of the auxiliary displacements

So far, many techniques have been proposed to compute the T-stress in the background of numerical simulations, e.g., stress extrapolation [28], displacement extrapolation [60] and the interaction integral method [21]. In this work, the domain-independent interaction integral [21] is taken for multi-branched and intersecting crack analysis. In the interaction integral method, two admissible states of the cracked body are considered, that is, the actual state and the auxiliary state, and the interaction integral can be obtained by the linear superposition of the J-integral of the two states as [21] [ ( ) ] 1 𝐼= + 𝑢aux (12) 𝜎 𝑢aux + 𝜎𝑖𝑗aux 𝑢𝑖,1 − 𝜎𝑖𝑘 𝑢aux 𝛿1𝑗 𝑞,𝑗 𝑑𝐴 𝑖,𝑘 𝑘,𝑖 ∫𝐴 𝑖𝑗 𝑖,1 2

To compute the interaction integral in Eq. (12), the derivatives of the auxiliary displacements should be figured out. From Eqs. (16) and (17), they are derived as 𝜕𝑢aux 1 𝜕 𝑥1 𝜕𝑢aux 1 𝜕 𝑥2 𝜕𝑢aux 2 𝜕 𝑥1

where A is the closed domain bounded by the contour C = C0 ∪C1 ∪IJ∪MN as illustrated in Fig. 3. 𝜎 ij and ui are, respectively, the stress and displacement of the true state expressed in the crack tip coordinate system x1 ox2 as shown in Fig. 3. The subscript , j (j = 1, 2) symbolizes the partial derivative with respect to xj . 𝛿 ij is the Kronecker delta being unity while i = j and zero otherwise. q(x1 , x2 ) is a sufficiently smooth function taking a value of unity on C0 , zero on C1 and arbitrary elsewhere.

𝜕𝑢aux 2 𝜕 𝑥2

=−

𝐹 cos 𝜃 (𝜅 − 1 + 2 cos 2𝜃) 8𝜋𝜇𝑟

(19)

=−

𝐹 sin 𝜃 (𝜅 + 3 + 2 cos 2𝜃) 8𝜋𝜇𝑟

(20)

=

𝐹 sin 𝜃 (𝜅 − 1 − 2 cos 2𝜃) 8𝜋𝜇𝑟

=−

𝐹 cos 𝜃 (𝜅 − 1 − 2 cos 2𝜃) 8𝜋𝜇𝑟

(21)

(22)

According to Eqs. (12), (18) and (19)–(22), we can observe that the value of F and d in Eq. (16) has no influence on the result of T-stress. 152

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Fig. 4. Determination of the interaction integral domain A and the smooth function q.

4.2. Determination of interaction integral domain To obtain the interaction integral in Eq. (12), the proper choice of integral domain A in Fig. 3 in the setting of the NMM is critical. Although the interaction integral is theoretically domain-independent, its numerical evaluation may depend on the size of the integral domain, as revealed in [49,51,52]. Take the leftmost crack tip in Fig. 2(a) for example, to determine the associated T-stress, herein a circle of radius rI centered at the very tip (see Fig. 4) is drawn and all the MEs interior or intersecting with this circle then comprise the domain A. Generally, the radius rI is taken as √ 𝑟𝐼 = 𝛼 𝐴t (23)

Fig. 5. A three-branched crack in a rectangular plate under uniaxial tension.

used to discretize the physical domains. Accordingly, the weight function wJ (X) in Eq. (6) takes the shape function of the four-nodded rectangular finite element [51]. Moreover, the polynomial basis in Eq. (8) is chosen to be constant and the material constants are fixed at E = 1.0 and 𝜈 = 0.3. Besides, the displacement boundary conditions are prescribed through the penalty method [47].

where At is the area of the ME containing the considered crack tip and 𝛼 is a scaling factor which will be further addressed in the first numerical example. 4.3. Determination of q(x1 , x2 ) and its derivatives

5.1. A three-branched crack under uniaxial tension

In the NMM, the smooth function q(x1 , x2 ) in Eq. (12) is such chosen that its value is 1.0 at the MP stars inside the circle in Fig. 4 and 0.0 at the stars outside it (the star of an MP is the point where the weight function wJ (X) in Eq. (6) is unity. For example, the square and the diamond in Fig. 4 denote, respectively, the star of the MP MA and MD ). Once the values on the associated stars are established, the corresponding result at any point interior an ME involved in the interaction integral is calculated by interpolation, that is, 𝑞(𝑥1 , 𝑥2 ) =

𝑁 ∑ 𝐽 =1

𝑤𝐽 𝑞𝐽

In the first example, as shown in Fig. 5, a crack with three branches in a rectangular plate subjected to uniform tension 𝜎 is tested. The dimension of the plate is 2W × 2H. The crack length of the horizontal branch is a and the two inclined segments are of the same size b. Besides, the included angles between the oblique branches and the X1 axis are both 𝜃. The SIFs for similar problem was studied by the NMM and the results agree well with the reference solutions [47]. Herein, we further examine the associated T-stress. Moreover, since the two inclined branches are symmetrical, we mainly care about crack tip A and B in the simulation. In the NMM modeling, the following parameters are adopted: W = 5.0, H = 4.0, 𝜎 = 1.0. Firstly, the convergence of the method and the domain-independence of the interaction integral are inspected. For this part, 𝑎∕𝑊 = 𝑏∕𝑊 = 0.1 and 𝜃 = 45°. Four MCs of size h = 0.48, 0.24, 0.12, 0.09 are, respectively, employed to cover the plate (h is the edge length of the square mathematical element). The resulting numbers of PPs and MEs are listed in Table 1 and the discretized domain when h = 0.12 is plotted in Fig. 6. For each MC, the scaling factor 𝛼 in Eq. (23) is, respectively, taken as 1.0, 1.5, 2.0, 2.5… and 6.0. The computed T-stresses are normalized by 𝜎 and the results for tip A and B are summarized, respectively, in Tables 2 and 3, together with the analytical values given by Chen and Lin [62] for the same cracks in an infinite plate (seeing that the crack lengths a and b are much smaller than the plate width in our simulation, it is reasonable to choose the results in infinite plate as reference solutions).

(24)

where qJ is the corresponding value at the star of the MP generating the PP PJ . From Eq. (24), the derivatives of the q(x1 , x2 ) can be further deduced as 𝑞,𝑗 =

𝑁 ∑ 𝐽 =1

𝑤𝐽 ,𝑗 𝑞𝐽

(25)

Therefore, according to the assignment of q(x1 , x2 ) and combining with Eqs. (12) and (25), we can readily see that only the MEs intersecting with the circle really contribute to the interaction integral, since the term q, j in Eq. (12) is identically vanishing for MEs interior the circle. 5. Numerical examples In this section, three numerical cases are conducted in plane strain state. At first, a crack with three branches is considered to verify the convergency of the proposed method and the domain-independence of the interaction integral. Then, another two problems with a four-branched and a six-branched crack are, respectively, studied to inspect the influences of crack geometry and boundary conditions on the T-stress. Throughout the simulation, uniform MCs formed by square mathematical elements and being inconsistent with all domain boundaries are

Table 1 Number of PPs and MEs at different MC size.

153

MC size h

0.48

0.24

0.12

0.09

Number of PPs Number of MEs

422 381

1517 1436

5806 5643

10,319 10,097

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Fig. 6. Discretized plate when h = 0.12: (a) global view and (b) enlarged zone around the crack.

Fig. 7. A four-branched crack in a square plate under biaxial tension: (a) physical model and (b) discretization of the plate when a = b = 0.5 and h = 0.03.

Table 2 Normalized T-stress for tip A at different MC size and scaling parameter.

Table 3 Normalized T-stress for tip B at different MC size and scaling parameter.

𝛼

h = 0.48

h = 0.24

h = 0.12

h = 0.09

𝛼

h = 0.48

h = 0.24

h = 0.12

h = 0.09

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0

−1.0241 −0.9794 −0.9864 −0.9856 −1.0053 −1.0069 −1.0068 −1.0072 −1.0076 −1.0076 −1.0077

−0.9067 −0.9098 −0.9188 −0.9159 −0.9155 −0.9748 −0.9930 −1.0067 −1.0093 −1.0096 −1.0103

−0.9803 −0.9046 −0.8925 −0.8891 −0.8876 −0.8875 −0.8862 −0.8942 −0.8945 −0.9112 −0.9202

−1.0161 −0.9220 −0.8898 −0.8916 −0.8879 −0.8881 −0.8870 −0.8872 −0.8860 −0.8882 −0.8885

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0

−0.0193 −0.0879 −0.0408 0.0050 −0.0004 0.0022 0.0032 0.0029 0.0030 0.0030 −0.0022

−0.2913 −0.3119 −0.2890 −0.2180 −0.0794 −0.0568 −0.0080 −0.0053 −0.0028 −0.0029 −0.0024

−0.3545 −0.3536 −0.3313 −0.3266 −0.3260 −0.3264 −0.3267 −0.3274 −0.2741 −0.2632 −0.0997

−0.3960 −0.3228 −0.3316 −0.3317 −0.3308 −0.3309 −0.3318 −0.3312 −0.3325 −0.3324 −0.3301

The reference solution reported by Chen and Lin in [62] is −0.8741.

The reference solution reported by Chen and Lin in [62] is −0.3291.

As we can observe from Tables 2 and 3, for coarser MCs (e.g., h = 0.48), the results at all 𝛼 deviate much from the reference ones; while for finer MCs (e.g., h = 0.09), most of the values are very close to the analytical ones. Such kind of outcome is rational. As we know, the solution accuracy of crack problems strongly depends on the mesh density around the crack tip; therefore the use of very coarse MC cannot well reflect the local property in the near-tip zone. With the refinement of the MC, the T-stresses then gradually tend to the exact ones, which numerically demonstrates the convergence of the proposed method. However, as we can note, the results at some 𝛼 are not so satisfactory even for finer MC (e.g., refer to the values at 𝛼 = 1.0, 1.5 when h = 0.09, and 𝛼 = 1.0,

1.5, 4.5, 5.0...when h = 0.12); in fact, these phenomena are also well founded. For 𝛼 around the lower end (e.g., 𝛼 = 1.0), the MEs involved in the interaction integral are such close to the crack tip that the errors of the mechanical fields are larger; while for those near the upper end (e.g., 𝛼 = 6.0), the integral domain is so big that it may cover the MEs around other branch tips, which will also bring down the accuracy. Among the cases when 𝛼 is far from the two extremes (e.g., 𝛼 = 2.0 ∼ 4.0), at finer MC (e.g., h = 0.12, 0.09), the simulated T-stress is very close to each other, which numerically infers the domain-independence of the interaction integral in Eq. (12). Accordingly, for MC with proper density (e.g., each crack branch is discretized into more than 5 parts as in the case 154

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Table 4 Normalized T-stresses for tip A and B at different inclined angle. Crack tip

Method

𝜃 = 15°

𝜃 = 30°

𝜃 = 45°

𝜃 = 60°

𝜃 = 75°

Tip A Tip B

NMM Chen and Lin [62] NMM Chen and Lin [62]

−0.9927 −0.9759 −1.2484 −1.2319

−0.9531 −0.9356 −0.8495 −0.8382

−0.8916 −0.8741 −0.3317 −0.3291

−0.8130 −0.7936 0.1833 0.1782

−0.7235 −0.7024 0.5546 0.5444

Fig. 10. Normalized T-stress for tip B at different branched length b.

Fig. 8. Normalized T-stresses for tip A, C and D at different branched length b.

Fig. 11. Normalized T-stress for tip C at different branch length a. Fig. 9. Normalized T-stress for tip A at different branch length a.

5.2. A four-branched crack under biaxial tension when h = 0.12, 0.09), the choice of 𝛼 = 2.0 ∼ 4.0 is preferred and similar findings were also obtained in the computation of SIFs [49,51,52,54]. Further considering that larger 𝛼 means the participation of more MEs in the interaction integral and thus causes higher modeling cost, hereafter we take 𝛼 = 2.5. Next, the accuracy of the proposed method is further tested by calculating the T-stress at various branched angles. Therefore, different inclinations with 𝜃 = 15°, 30°, 45°, 60°, 75° are, respectively, inspected at 𝑎∕𝑊 = 𝑏∕𝑊 = 0.1 and h = 0.09. The normalized T-stresses by the NMM together with the analytical solutions from [62] are provided in Table 4, from which we can catch that the present results agree well with the theoretical ones. Furthermore, it is found that the T-stress at both crack tips increases with the inclination.

As illustrated in Fig. 7(a), two intersecting cracks with four branches in a square plate loaded by biaxial tension 𝜎 = 1.0 are studied. The width of the plate is 2W = 2.0, and the size of the horizontal and inclined branches, is respectively, a and b. Additionally, the included angle between the oblique branch and the horizontal one is 𝜃. In addition, the plate is fixed on the bottom side. Firstly, the influences of the branched length b on the T-stress are concerned. For this, several cases with b = 0.2, 0.3, 0.4, ..., 0.8 are separately considered with a = 0.5 and 𝜃 = 90°. Uniform MC of size h = 0.03 is constructed to cover the plate. The associated discretization parameters are reported in Table 5 and the discretized domain when b = 0.5 is plotted in Fig. 7(b). Because of symmetry, we mainly focus on crack tip A, C and D. The corresponding T-stresses normalized by 𝜎 are 155

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Fig. 14. Normalized T-stress for tip A at different branch length and loading ratio.

Fig. 12. Normalized T-stress for tip D at different inclined angle.

Table 5 Number of PPs and MEs at different branch length. b

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Number of PPs Number of MEs

4440 4268

4456 4276

4468 4282

4480 4288

4492 4294

4508 4302

4520 4308

tips move closer to the boundaries, which is due to stronger impact of the boundaries. 5.3. A six-branched crack under biaxial tension In the last example, as shown in Fig. 13(a), a star-shaped crack with six branches in a square plate under biaxial tension is investigated. The width of the plate is 2W. All branches are of the same length a and the included angles between the adjacent segments are all 𝜃. The boundary condition is such prescribed that the uniform tension along the vertical and horizontal direction is, respectively, 𝜎 and k𝜎 (k is a loading ratio), and the bottom of the plate is clamped. In the NMM modeling, W = 1.0, 𝜃 = 60°, 𝜎 = 1.0. Herein, we mainly inspect the effects of the branch length and the loading ratio on the Tstress. For this purpose, the T-stresses under different combinations of branch size and loading ratio, with a = 0.2, 0.3, 0.4, ..., 0.8 and k = 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0 are, respectively, calculated at h = 0.03. The discretization data is summarized in Table 6 and the discretized plate when a = 0.5 is illustrated in Fig. 13(b). Seeing the symmetry, we just plot the results of T/𝜎 at tip A, C and E in Figs. 14–16 (theoretically, the values at tip B, D and F are, respectively, equal to those at tip A, E and C), from which we can tell that the influences of the loading ratio on the T-stress at a given branch size are similar for all three tips, while those

presented in Fig. 8. From this figure, it is easy to find that the T-stress at tip A decreases with the increased b; the value of tip C first reduces then increases, while the result for tip D undergoes a trough, then reaches a peak and finally dramatically falls. Then, by assuming that the four branches are equal-sized, the effects of crack length and oblique angle on the T-stress are analyzed. To this end, different crack lengths with a = b = 0.2, 0.3, 0.4, ..., 0.8 are, respectively, considered, and for each branch size, selected inclinations with 𝜃 = 15°, 30°, 45°, 60°, 75°, 90° are further examined on the same MC of size h = 0.03. The associated results are plotted in Figs. 9–12, from which we can detect that for a given branch length, the T-stress increases as the inclination angle grows larger, except individual cases for crack tip D when a = 0.7, 0.8. While for a specified inclined angle, the T-stress is reduced with the raising crack length in all the cases but the situations when 𝜃 = 60°, 75°, 90° for tip C, where the result increases with a; what is more, we can also view that the T-stress varies faster when the branch

Fig. 13. A six-branched crack in a square plate under biaxial tension: (a) physical model and (b) discretization of the plate when a = 0.5 and h = 0.03.

156

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

Table 6 Number of PPs and MEs at different branch length. a

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Number of PPs Number of MEs

4444 4272

4492 4296

4534 4316

4586 4342

4637 4367

4685 4391

4729 4413

gency, accuracy and also the domain-independence of the interaction integral were numerically tested through a typical numerical example, and then the effects of crack configurations and boundary conditions on the T-stress were analyzed via another two cases. Good accuracy of the proposed method was displayed through the comparison with published solutions in the literature. The present work can be further developed to study the T-stress for cracked nonhomogeneous materials under both mechanical and thermo-mechanical loadings. Acknowledgments The present work was supported by the National Natural Science Foundation of China (Grant nos. 11462014, 51778021 and 11572282), the Provincial Natural Science Foundation of Jiangxi, China (Grant no. 20151BAB202003), the Science and Technology Program of Educational Committee of Jiangxi Province of China (Grant no. GJJ180533) and the Graduate Innovation Foundation of Nanchang Hangkong University (Grant no. YC2018069). References [1] Phan AV. A non-singular boundary integral formula for determining the T-stress for cracks of arbitrary geometry. Eng Fract Mech 2011;78(11):2273–85. [2] Larsson SG, Carlsson AJ. Influence of non-singular stress terms and specimen geometry on small-scale yielding at crack tips in elastic-plastic materials. J Mech Phys Solids 1973;21(4):263–77. [3] Rice JR. Limitations to small-scale yielding approximation for crack tip plasticity. J Mech Phys Solids 1974;22(1):17–26. [4] Cotterell B, Rice JR. Slightly curved or kinked cracks. Int J Fract 1980;16(2):155–69. [5] Ueda Y, Ikeda K, Yao T, Aoki M. Characteristics of brittle-fracture under general combined modes including those under bi-axial tensile loads. Eng Fract Mech 1983;18(6):1131–58. [6] 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 Mech 2001;24(2):137–50. [7] Gupta M, Alderliesten RC, Benedictus R. A review of T-stress and its effects in fracture mechanics. Eng Fract Mech 2015;134:218–41. [8] Huang XL, Liu YH, Dai YW. Characteristics and effects of T-stresses in central-cracked unstiffened and stiffened plates under mode I loading. Eng Fract Mech 2018;188:393–415. [9] Hutchinson JW, Suo Z. Mixed-mode cracking in layered materials. In: Hutchinson JW, Wu TY, editors. Advances in Applied Mechanics, 29. San Diego: Elsevier Academic Press Inc; 1992. p. 63–191. [10] Yang B, Ravi-Chandar K. Crack path instabilities in a quenched glass plate. J Mech Phys Solids 2001;49(1):91–130. [11] Tong J. T-stress and its implications for crack growth. Eng Fract Mech 2002;69(12):1325–37. [12] Sobotka JC, Dodds RH. T-stress effects on steady crack growth in a thin, ductile plate under small-scale yielding conditions: three-dimensional modeling. Eng Fract Mech 2011;78(6):1182–200. [13] Ayatollahi MR, Sedighiani K. A T-stress controlled specimen for mixed mode fracture experiments on brittle materials. Eur J Mech A – Solid 2012;36:83–93. [14] Dyskin AV. Crack growth criteria incorporating non-singular stresses: size effect in apparent fracture toughness. Int J Fract 1997;83(2):191–206. [15] Saghafi H, Monemian S. A new fracture toughness test covering mixed-mode conditions and positive and negative T-stresses. Int J Fract 2010;165(1):135–8. [16] Zhao YH, Dong W, Xu BH, Liu J. Effect of T-stress on the initial fracture toughness of concrete under I/II mixed-mode loading. Theor Appl Fract Mech 2018;96:699–706. [17] Berg E, Skallerud B, Thaulow C. Two-parameter fracture mechanics and circumferential crack growth in surface cracked pipelines using line-spring elements. Eng Fract Mech 2008;75(1):17–30. [18] Seitl S, Vesely V, Routil L. Two-parameter fracture mechanical analysis of a near-crack-tip stress field in wedge splitting test specimens. Comput Struct 2011;89(21–22):1852–8. [19] Kravchenko SG, Kravchenko OG, Sun CT. A two-parameter fracture mechanics model for fatigue crack growth in brittle materials. Eng Fract Mech 2014;119:132–47. [20] Chen CS, Krause R, Pettit RG, Banks-Sills L, Ingraffea AR. Numerical assessment of T-stress computation using a p-version finite element method. Int J Fract 2001;107(2):177–99. [21] Kim JH, Paulino GH. T-stress, mixed-mode stress intensity factors, and crack initiation angles in functionally graded materials: a unified approach using the interaction integral method. Comput Method Appl Mech 2003;192(11–12):1463–94. [22] Jogdand PV, Murthy KSRK. A finite element based interior collocation method for the computation of stress intensity factors and T-stresses. Eng Fract Mech 2010;77(7):1116–27. [23] Bouchard PO, Bernacki M, Parks DM. Analysis of stress intensity factors and T-stress to control crack propagation for kerf-less spalling of single crystal silicon foils. Comp Mater Sci 2013;69:243–50. [24] Yu PS, Wang Q, Zhang CF, Zhao JH. Elastic T-stress and I-II mixed mode stress intensity factors for a through-wall crack in an inner-pressured pipe. Int J Pres Ves Pip 2018;159:67–72.

Fig. 15. Normalized T-stress for tip C at different branch length and loading ratio.

Fig. 16. Normalized T-stress for tip E at different branch length and loading ratio.

of the crack length at a fixed loading ratio is more complex. In detail, the T-stress decreases with increasing k for a constant branch length in all the cases. For an assumed k, the values for crack tip A is reduced with the increased a but some special situations when k = 10.0; as for tip C, the T-stress slowly increases then decreases for k < 1.0, while first decreases then increases when k ≥ 1.0; concerning tip E, the monotonic reduction is found for k ≤ 2.0, while it experiences a trough and then raises in the remaining; moreover, we can also note that under most of the situations, the T-stress is not so sensitive to branch length when a ≤ 0.5. 6. Concluding remarks The elastic T-stresses for multiple-branched and intersecting cracks have been evaluated by the numerical manifold method. The displacement discontinuity across crack surface was described in nature by the NMM, and the near-tip stress singularity was addressed with the use of the crack tip asymptotic basis in the corresponding local functions. The domain-form interaction integral for T-stress computation was introduced and its numerical implementations were provided. The conver157

H.H. Zhang, S.M. Liu and S.Y. Han et al.

Engineering Analysis with Boundary Elements 107 (2019) 149–158

[25] Sladek J, Sladek V. Evaluations of the T-stress for interface cracks by the boundary element method. Eng Fract Mech 1997;56(6):813–25. [26] Tan CL, Wang X. The use of quarter-point crack-tip elements for T-stress determination in boundary element method analysis. Eng Fract Mech 2003;70(15):2247–52. [27] Guduru V, Phan AV, Tippur HV. Transient analysis of the DSIFs and dynamic T-stress for particulate composite materials-Numerical vs. experimental results. Eng Anal Bound Elem 2010;34(11):963–70. [28] Chen M, Xu Z, Fan XM. Evaluation of the T-stress and stress intensity factor for multi-crack problem using spline fictitious boundary element alternating method. Eng Anal Bound Elem 2018;94:69–78. [29] Yu HJ, Wu LZ, Guo LC, Li H, Du SY. T-stress evaluations for nonhomogeneous materials using an interaction integral method. Int J Numer Methods Eng 2012;90(11):1393–413. [30] Guo FN, Guo LC, Huang K, Bai XM, Zhong SY, Yu HJ. An interaction energy integral method for T-stress evaluation in nonhomogeneous materials under thermal loading. Mech Mater 2015;83:30–9. [31] Song CM. Evaluation of power-logarithmic singularities, T-stresses and higher order terms of in-plane singular stress fields at cracks and multi-material corners. Eng Fract Mech 2005;72(10):1498–530. [32] Song CM, Vrcelj Z. Evaluation of dynamic stress intensity factors and T-stress using the scaled boundary finite-element method. Eng Fract Mech 2008;75(8):1960–80. [33] Natarajan S, Song CM, Belouettar S. Numerical evaluation of stress intensity factors and T-stress for interfacial cracks and cracks terminating at the interface without asymptotic enrichment. Comput Method Appl Mech 2014;279:86–112. [34] Shi GH. Manifold method of material analysis. In: Proceedings of the transactions of the ninth army conference on applied mathematics and computing; 1991. p. 57–76. [35] Chiou YJ, Lee YM, Tsay RJ. Mixed mode fracture propagation by manifold method. Int J Fract 2002;114(4):327–47. [36] Li SC, Li SC, Cheng YM. Enriched meshless manifold method for two-dimensional crack modeling. Theor Appl Fract Mech 2005;44(3):234–48. [37] Terada K, Ishii T, Kyoya T, Kishino Y. Finite cover method for progressive failure with cohesive zone fracture in heterogeneous solids and structures. Comput Mech 2007;39(2):191–210. [38] Wu ZJ, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39:38–53. [39] Wu ZJ, Wong LNY. Extension of numerical manifold method for coupled fluid flow and fracturing problems. Int J Numer Anal Met 2014;38(18):1990–2008. [40] Wu ZJ, Xu XY, Liu QS, Yang YT. A zero-thickness cohesive element-based numerical manifold method for rock mechanical behavior with micro-Voronoi grains. Eng Anal Bound Elem 2018;96:94–108. [41] Zheng H, Liu F, Li CG. The MLS-based numerical manifold method with applications to crack analysis. Int J Fract 2014;190(1–2):147–66. [42] Zheng H, Liu F, Du XL. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Comput Method Appl Mech 2015;295:150–71. [43] Li W, Zheng H, Sun GH. The moving least squares based numerical manifold method for vibration and impact analysis of cracked bodies. Eng Fract Mech 2018;190:410–34.

[44] Li X, Zhang QB, He L, Zhao J. Particle-based numerical manifold method to model dynamic fracture process in rock blasting. Int J Geomech 2017;17(5). [45] Yang YT, Tang XH, Zheng H, Liu QS, Liu ZJ. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl Math Model 2018;53:462–86. [46] Yang SK, Cao MS, Ren XH, Ma GW, Zhang JX, Wang HJ. 3D crack propagation by the numerical manifold method. Comput Struct 2018;194:116–29. [47] Ma GW, An XM, Zhang HH, Li LX. Modeling complex crack problems using the numerical manifold method. Int J Fracture 2009;156(1):21–35. [48] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation problems using the numerical manifold method. Eng Anal Bound Elem 2010;34(1):41–50. [49] Zhang HH, Zhang SQ. Extract of stress intensity factors on honeycomb elements by the numerical manifold method. Finite Elem Anal Des 2012;59:55–65. [50] An XM, Zhao ZY, Zhang HH, He L. Modeling bimaterial interface cracks using the numerical manifold method. Eng Anal Bound Elem 2013;37(2):464–74. [51] Zhang HH, Ma GW. Fracture modeling of isotropic functionally graded materials by the numerical manifold method. Eng Anal Bound Elem 2014;38:61–71. [52] Zhang HH, Ma GW, Ren F. Implementation of the numerical manifold method for thermo-mechanical fracture of planar solids. Eng Anal Bound Elem 2014;44:45– 54. [53] Zhang HH, Ma GW, Fan LF. Thermal shock analysis of 2D cracked solids using the numerical manifold method and precise time integration. Eng Anal Bound Elem 2017;75:46–56. [54] Zhang HH, Liu SM, Han SY, Fan LF. Modeling of 2D cracked FGMs under thermo-mechanical loadings with the numerical manifold method. Int J Mech Sci 2018;148:103–17. [55] Zhang HH, Liu SM, Han SY, Fan LF. The numerical manifold method for crack modeling of two-dimensional functionally graded materials under thermal shocks. Eng Fract Mech 2019;208:90–106. [56] Liu F, Zheng H, Du XL. Hybrid analytical and MLS-Based NMM for the determination of generalized stress intensity factors. Math Probl Eng 2015;2015. [57] He J, Liu QS, Ma GW, Zeng B. An improved numerical manifold method incorporating hybrid crack element for crack propagation simulation. Int J Fract 2016;199(1):21–38. [58] Sladek J, Sladek V. Evaluation of T-stresses and stress intensity factors in stationary thermoelasticity by the conservation integral method. Int J Fract 1997;86(3):199–219. [59] Yates JR, Zanganeh M, Tai YH. Quantifying crack tip displacement fields with DIC. Eng Fract Mech 2010;77(11):2063–76. [60] Chen CH, Wang CL. Stress intensity factors and T-stresses for offset double edge-cracked plates under mixed-mode loadings. Int J Fract 2008;152(2):149–62. [61] Sutradhar A, Paulino GH. Symmetric Galerkin boundary element computation of T-stress and stress intensity factors for mixed-mode cracks by the interaction integral method. Eng Anal Bound Elem 2004;28(11):1335–50. [62] Chen YZ, Lin XY. Evaluation of the T-stress in branch crack problem. Int J Fract 2010;161(2):175–85.

158