mechanical loading

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

Crack kinking out of interface of two orthotropic materials under combined thermal/mechanical loading

T

Oldřich Ševeček, Michal Kotoul , Tomáš Profant, Miroslav Hrstka ⁎

Brno University of Technology, Institute of Solid Mechanics, Mechatronics and Biomechanics, Technická 2896/2, 616 69 Brno, Czech Republic

ARTICLE INFO

ABSTRACT

Keywords: Interface crack Finite fracture mechanics Crack kinking Fracture criterion FEM T-stress

The problem of crack path stability along the interface between two orthotropic elastically dissimilar materials under the presence of in-plane residual stresses is analyzed using the concept of Finite Fracture Mechanics and matched asymptotic procedure. The complex stress intensity factor and the T-stress characterizing the stress state at the crack tip are calculated both for the thermal (residual stresses) and mechanical loading using the two-state integral. The matched asymptotic procedure together with FEM is used to derive the change of the potential energy induced by the incremental crack growth (and correspondingly the incremental energy release rate of such a crack), which is needed for the fracture criterion formulation. An energy based fracture criterion is introduced for this problem and it is investigated whether and how is the criterion for the prediction of crack kinking from the interface affected by residual stresses.

1. Introduction Kinking of a crack out of a bi-material interface in layered systems considering the influence of residual stresses was investigated in the classical paper by He et al. [1]. They represented the effect of residual stresses by the non-singular term known as T-stress in the series expansion of the stress components near the tip of the parent crack. They were able to show that compressive T-stresses are beneficial to permit higher interface fracture energy to achieve crack deflection and to prevent crack kinking, respectively. The validity of this approach requires that the residual stresses do not contribute to the singular stress field. As pointed out by Leguillon [2], the main role of residual thermal stresses in fracture problems of composite materials consists in modifying the stress intensity factor, which by the principle of superposition, can be expressed as the sum of mechanical loading contribution and thermal stresses contribution. The criterion governing the crack path selection remains unchanged if in the asymptotic analysis only K-dominant field is considered. However, if higher order terms of the asymptotic expansion are needed (e.g. when the fracture process is not confined well within the boundary of the K-dominant field), then under the combined influence of mechanical loading and thermal stresses the criterion evolves from a local formulation (i.e. depending only on local properties around the primary crack tip) to non-local one in which the geometry of the whole body plays an important role. On the contrary, if mechanical and thermal loadings are taken separately,

⁎

the criterion still preserves the local character. Lee at al. [3] reconsidered the kinking out of the interface problem in layered systems loaded in four-point bending under the presence of in-plane residual stress using the Finite Element Method (FEM) and assuming non-oscillatory character of the singular field of the parent interface crack. It was shown that the residual stress develops its own singular field near the tip of the parent interface crack thereby causing the transition of the local stress field from one state to the other which is reflected in the change of the phase angle and the strain energy release rate. Assuming that the K-dominance holds, they showed that regardless of the magnitude of the residual stress, there is a unique relation between the maximum permissible level of interface fracture energy relative to the bulk fracture energy and the phase angle of loading for a given material combination. However, since the total phase angle depends on the global distribution of the residual stresses, the fracture process within a K-dominant field can no longer be treated as local problem alone. The singularity of the stress field at the interface crack tip has the form r1/2+iε, with r being the radial coordinate emanating from the crack tip. Analysis of the interface crack kinking usually brings difficulties with extension of the classical criteria for fracture in homogeneous materials in situations when the oscillatory index ε defined by

=

1 1 ln 2 1+

,

(1)

where β is the second Dundurs parameter, different from zero. These

Corresponding author. E-mail address: [email protected] (M. Kotoul).

https://doi.org/10.1016/j.tafmec.2019.102397 Received 4 September 2019; Received in revised form 21 October 2019; Accepted 22 October 2019 Available online 07 November 2019 0167-8442/ © 2019 Elsevier Ltd. All rights reserved.

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

Nomenclature

y1, y2 z α, β αij, αijeff Γd, Γi

scaled up coordinates (for zoomed-in view) argument of the complex potential function first and second Dundurs’ parameter coefficients of thermal expansion (CTE), effective CTE toughness of the matrix material of intact layer, toughness of the interface I , II I , II 1 ( ), i ( ) angular function of the singular stress term and its auxiliary counterpart I , II ( ), 1I , II ( ) angular function of the displacement field of the 1 singular term and its auxiliary counterpart I , II ( ), iI , II ( ), i = 2, 3 angular functions of the displacement i field of the T-stress term and their auxiliary counterparts Δa kinked crack length δW change of the potential energy ΔT change of the temperature εij, εel components of the strain tensor, elastic strain ε0 initial strain ε oscillation index θ angular polar coordinate I , II , I , II characteristics of the material orthotropy (layer I and II) μi material eigenvalues (complex numbers) νI, νII, ν́ Poisson’s ratio for material I and II, in-plane Poisson’s ratio ξ residual stress magnitude parameter ρ zoomed polar coordinate r with factor τ ρI, ρII composite parameters for material I and II I , II components of the singular field stresses 1, ij I , II non-zero component of the stress field corresponding to T2,11 stress term σkl(u ) stress tensor, auxiliary solution for stresses σ0,ij initial residual stress tensor components in the uncracked specimen σr,ij residual stress tensor components in the cracked specimen τ perturbation parameter ϕjJ vector of complex coefficients ψ, ψ́ local phase angles of SIF associated with different reference length l or ĺ appl , res , tot local phase angles of applied, residual and total SIF

a length of the interfacial parent crack A1, A2, A3, A4 auxiliary parameters introduced in Eq. (6) matrix of material characteristic eigenvalues and compoA, Aij nents sij A(ω),B(ω),C(ω),D(ω) functions in the expression for the change of potential energy induced by a finite crack increment B Hermitian matrix c,d dimensionless complex parameters in Eq. (2) C, C1, C2 integration contours É , E Young’s modulus in longitudinal and transversal directions D integration domain fi(τ) coefficients of the outer asymptotic expansion f(zj), fj(zj) complex potential and its components P applied mechanical force coefficients of the inner asymptotic expansion Fi(τ) G, Ǵ out of plane and in plane shear modulus G( ) energy release rate in case of kinked crack Gi = G (0) energy release rate in case of an interface crack Gmax maximal value of ERR H, H11, H22 Bimaterial Hermitian matrix, components of the H matrix i imaginary unit I, II material I, Material II indices I(.,.) path independent integral kappl, kres normalized stress intensity factors Kij coefficients of the perturbed solution K,K SIF, Complex interfacial SIF Ktot, Kres, Kappl total SIF, thermal loading SIF, mechanical loading SIF l reference length L characteristic length L, Liα matrix containing material characteristic eigenvalues. nl normal unit vector components r polar coordinate R(θ),Ψ(θ) polar coordinates in the complex variable formalism p2I , II , p3I , II real coefficients in T-stress expressions sij components of compliance matrix S compliance matrix Si, So inner/outer span of the supports upon the 4 point bending thickness of material layer I and II tI, tII Ttot, Tres, Tappl total T-stress, T-stress from thermal and mechanical loading j j tappl , tres normalized T-stresses initial displacement u0 u displacement vector auxiliary displacement vector u wi eigenvector components W potential energy W0, Wτ potential energies of unperturbed and perturbed state x1, x2, x3 Cartesian coordinates

appl ,

res ,

tot modified local phase angles in Eq. (21) of applied, residual and total SIF ω crack kink angle Ω, Ωτ, Ωin perturbation domain, Inner domain

Abbreviations: CTE FEM FFM ERR GSIF LEFM SIF

difficulties are associated with the oscillatory character of the singular elastic solution for the parent interface crack. Due to this oscillatory character, the fracture criteria based upon maximum of the energy release rate (ERR), local symmetry (KII = 0), maximum KI or maximum circumferential stress do not predict a unique value of kink angle ω without specifying a fixed kink crack length. Different kink lengths imply different kink angles values predicted by these criteria. Applying dimensional analysis, He and Hutchinson [4] derived the expressions for Stress Intensity Factors (SIFs) at a straight kink crack of length Δa much smaller in comparison to the parent crack length 2a as functions

Coefficient of Thermal Expansion Finite Element Method Finite Fracture Mechanics Energy Release Rate Generalized Stress Intensity Factor Linear Elastic Fracture Mechanics Stress Intensity Factor

of the complex interfacial SIF K = Kli :

kI ( , a) + ikII ( , a) = c ( , , ) K

a l

i

+ d ( , , )K

a l

i

(2) where α is the first Dundurs parameter and c and d are two dimensionless complex parameters and l is the reference length introduced by Rice [5]. According to Eq. (2), local SIFs are oscillating with kink length Δa and consequently no limit exists for Δa → 0 when ε ≠ 0. They also suggested to overcome the difficulties caused by the 2

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

oscillatory character of the expression by arbitrarily taking β = 0, especially when β is small. They showed that such an approach is reasonable because the fracture variables of interest depend weakly on β. Leguillon [6] suggested different approach based upon the concept of the Finite fracture mechanics (FFM). This allows getting rid of the problem brought by oscillations. The essence of FFM consists in the key assumption that crack propagation is a discontinuous process occurring in finite steps, rather than continuously and smoothly as in the traditional Linear Elastic Fracture Mechanics (LEFM) theory. Mathematically, instead of using the differential form of the Griffith energy balance, the integral formulation of the Griffith criterion is applied. Besides the energy balance also a stress condition is required to be simultaneously fulfilled. Consequently, the crack advance becomes a structural parameter which allows taking into account the interaction between the finite crack advance and the geometry of the specimen. We will extend this concept to orthotropic bi-materials including an effect of the T-stress generated due to applied loading as well as due to residual stresses. The matched asymptotic expansion method is applied to evaluate the change of potential energy connected with a short crack kink. An expression for the change of potential energy and/or the energy release rate (ERR), which does not suffer from the oscillations mentioned above, is obtained.

modulus is G′ and in-plane Poisson’s ratio is ν′. The compliance matrix S then reads 1 E E

S=

0

0

0

0

0

0

1 E

E

E

1 E

E

0

0

0

0

0

0

1 G

0

0

0

0

0

0

1 G

0

0

0

0

0

0

1 G

E

E

.

(3)

The theory used for the stress field description is derived purely for the plane strain or plane stress conditions. The 2D FEM mesh was generated in the finite element system ANSYS 18.1 [7]. Combined loading (thermal + mechanical) was simulated using the principle of superposition. The thermal loading due to cooling down from the processing temperature to room temperature induces residual stresses in the specimen. In case of the mechanical loading, the results obtained from 2D plane strain calculations entirely correspond to those obtained from the full 3D analysis in the centre of the specimen. However, in case of the thermal loading of laminate, certain disproportions are found between the results obtained with 2D and 3D model, respectively. The reason is that the plane strain conditions cannot be fulfilled in case of the thermal loading (because of the material extension in all directions) and the biaxial thermal residual stress σ11 ≠ 0, σ33 ≠ 0, σ22 = 0 is not reproduced. To overcome this problem while still taking an advantage of 2D calculations also for the thermal problems, certain modification of the material coefficients of thermal expansion (CTEs), in all directions, has to be performed. In the next, the initial residual stresses, strains, and displacements respectively in the uncracked specimen will be denoted by the subscript 0 separated by comma from coordinate indices. Assuming σ22 = 0 and considering the compliances in Eq. (3), one obtains:

2. Computational model A symmetric laminate made of three orthotropic layers is considered in our study in order to prevent bending moments after the cooling process. The outer layers are denoted by the index II, the inner layer is denoted by the index I. It is assumed that the principal axes of all layers are aligned and the bi-material interfaces are parallel or perpendicular to the principal axes. Since all of the materials' principal axes in the out of plane direction are assumed to be parallel with the x3-axis, the antiplane deformation can be decoupled from the in-plane deformation. Four-point bending delamination test geometry was adopted for the analysis, see Fig. 1. This is a convenient test method to characterise both the interface fracture property and residual stress at the same time. Two-dimensional plane strain geometry was assumed considering the beam-type structure of the laminates. Thermo-elastic properties of the layers I and II possess the tetragonal symmetry with elastic modulus E′ in the direction of x1-axis and elastic moduli E2 = E3 = E in x2 and x3-axis. The in-plane shear

j 0,11 j 0,33

=

=

Ej Ej

2 j Ej

E j Ej Ej

2 j Ej

(

(E j

el, j 0,33

el, j 0,11

+

+ Ej

el, j j 0,11),

el, j j 0,33),

j = I , II .

(4)

The residual stresses 0,11, 0,33 have to satisfy the following static equilibrium conditions derived for the three-layer geometry shown in

Fig. 1. Scheme of the analyzed four-point bending specimen, detail of the deformed FE mesh and considered tetragonal material properties of both layers. Additional j

elastic moduli EII , EI , EI are calculated for selected values of the Dundurs bi-material parameter α and the orthotropic parameters , j = I , II defined in Appendix A. Shear moduli GI , GII are then calculated for selected values of the orthotropic parameters j , j = I , II also defined in Appendix A. Negative thermal expansion coefficient is typical e.g. for some epoxy/graphene composites. 3

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

Fig. 1.

2 2

II 0,11 tII II 0,33 tII

I 0,11 tI I 0,33 tI

+ +

beginning. In case of the interface crack in a finite bi-material anisotropic body, FEM is mostly an indispensable tool to extract the stress intensity factors. The possible oscillatory behaviour of the near tip solution makes difficult to extract the stress intensity factors directly from the near tip field, hence several path independent integrals were proposed to compute the stress intensity factor which allow to avoid the complexity of stresses around the crack tip. The complex interfacial SIFs of the parent crack, Kapp, Kres, where Kapp is due to the applied load and Kres is due to residual stresses are calculated using the Betti reciprocal theorem [18–22]. The reciprocal theorem states for two arbitrary admissible fields u and u

= 0, (5)

= 0.

Substituting from Eq. (4) into Eq. (5) and making use of the relael, II el, I el, II el, I = 0,11 + 11, 0,33 = 0,33 + 33 , where 11, 33 are detions 0,11 I II I II fined as 11 = ( 11 33 = ( 33 11 ) T , 33 ) T , α11, α33 stand for thermal expansion coefficients and ΔT is an uniform cooling range from the processing temperature, one obtains el, I 0,11 (A1 el, I 0,11 (A3

+ A2 ) +

el, I 0,33 (A3

+ A4 ) +

el, I 0,33

(

A3 II

2tII EII2 , 2 EII II EII

where A1 =

el, I 0,33

+

el, I 0,11

11 [A1 (A3 I + A 4 II )

=

A4 I

A1

)= tI EI

A2 =

Solving Eq. (6) for el, I 0,11

+ A4 ) =

EI

and

11

A3

2

2 , I EI

el, I 0,33

A3

33,

A3

33,

11

A3 =

we get

II

(6)

2tII EII EII II , 2 EII II EII

A3 I II (A3 + A 4 )] + 33 A3 A4 ( II 2 I II (A3 + A 4 )

I)

(A1 + A2 )(A3 I + A4 II ) I[

=

33 A3 ( II (A3 + A 4 )

(A1 + A2 )) +

(A1 + A2 )(A3 I + A4 II )

11 II (A1 A 4 2

A2 A3 )]

I II (A3 + A 4 )

D

EI

,

.

kl (u ) kl (u )

kl (u ) kl (u )]dS

(7)

Consider now plane strain conditions 33 = 0 and seek the effective eff eff el, I , 33 CTEs 11 which lead to the value of 0,11 given by Eq. (7)1. The thermo-elastic constitutive law under plane strain conditions reads

j r

= C j( ) ·( r (urj )

eff , j

n·

22

(

=

E

+

(

2

)

E

E

+

)

0 11

2

1 E

0

E

0

E

22

11

(

12

1 G

0

+

(

+

22

E

E

+

33

33

)

I , eff 33

due el 0,33 +

eff 33

el, I 0,33

T=

to tetragonal T = 0 we have I[

= II , eff 33

I[

=

33 A3 ( II (A3 + A 4 )

,

) T

T=

33 A3 ( II (A3 + A 4 )

(A1 + A2 )) +

(A1 + A2 )(A3 I + A 4 II )

α22 = α33.

(A1 + A2 )) +

(A1 + A2 )(A3 I + A 4 II ) el, II el, I 0,33 = 0,33 11 II (A1 A 4 2 I II (A3 + A 4 )

11 II (A1 A 4 2

A2 A3 )]

I II (A3 + A 4 )

D = D I + D II

Since

j 11

+

j

Ej

Ej

eff , j 33

(14)

33.

where j = I , II

[(

j kl (u )

j j r , kl (u r )) kl (u )

j j kl (u )( kl (u )

j r , kl (u r ))]dS

(15) II

Clearly, the integration over the subdomains D and D is performed using pertinent stress and strain field components marked by superscripts I and II, respectively. Applying the divergence theorem and using FE approximation to the field u denoted by the superscript h we get

,

(9)

=

n1 0 n2 0 n1 n2

I

The effective CTE in the direction x1 is obtained from Eq. (8) if is substituted for α33: eff , j 11

, j = I , II , n =

= 0, j = I , II .

33= A2 A3 )]

= 0 on

satisfies compatibility. The components of the stiffness matrix C j( ) are given in Eq. (17) and the components of the effective thermal expansion vector eff , j are given in Eqs. (10) and (11). The subscript r refers to a thermal contribution and nl denotes a component of the unit normal vector. Then Eq. (12) can be rewritten as

T

0

symmetry

j r

j r (u r )

(8) where 0,33 =

(13)

and homogeneous boundary conditions

2 12 E E 2

(12)

T ), j = I , II

11

1 E

= 0,

where the domain D is any subset of the original domain obtained by excluding the crack tip. The boundary of D is made up of arbitrary contours circumventing the crack tip and connecting the traction free crack faces. Consider now initial thermal stresses r , kl (ur ) and strains I , II developed due to a uniform r , kl (u r ) in the body. Thermal stress r temperature change ΔT satisfies equilibrium, the thermo-elastic constitutive law

tI EI EI I 2 . I EI

A4 =

[

eff 33

(10)

Similarly, the effective CTE in the direction x2 is eff , j 22

=

eff , j 33 (1

+ j ), where j = I , II .

(11)

The above derived effective CTEs were subsequently employed in 2D modelling of the cracked four-point bending specimen by means of j , j = I , II as a function FEM. Fig. 2 shows the initial residual stress 0,11 of the Dundurs parameter α for a special choice of values of the orj thotropic parameter , j = I , II . 3. Determination and representation of the complex stress intensity factor and the T-stress for primary interface crack Starting with Suo classical paper [8], a great number of papers were devoted to the problem of interface crack between dissimilar anisotropic media, see e.g. [9–17] and citations herein. Lekhnitskii, Eshelby, and Stroh (LES) formalism, summarized in [16], was broadly applied to solve an anisotropic elasticity problem. To obtain SIFs as useful design parameters, a singular stress field needs to be developed at the

Fig. 2. Residual stress

j 0,11 ,

j = I, II as a function of the Dundurs parameter α for

several values of the orthotropic parameter

4

I

=

II

=

.

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

C = C1 C2

= +

(

h, j j kl (u ) nl uk

eff , j ) j [C j(,11 T 11 (u ) 11 D eff , j ( ) j C j,22 22 (u ) 22 T ]dS,

+

j h, j kl (u ) nl uk )ds=

eff , j ) C j(,12 ( 11j (u ) 22

=

j eff , j 22 (u ) 11

T+

T) (16)

Ej 2 ( j

1) , 1) + 2Ej j 2

Ej ( j

) C j(,22 =

E j Ej j

) C j(,12 =

Ej (Ej j 2 E j ) (1 + j )[2Ej j 2 + E j ( j

, 1) + 2Ej j 2

Ej ( j

1)]

, j = I , II .

(17)

The domain D is enclosed between the contours C1 and C2 which are arbitrary. Observe that the term on the right-hand side of Eq. (16) disappears in case of pure mechanical loading. Reversing the flow direction of the contour C1 we obtain: C1

h, j j kl (u ) nl uk

( +

D

+

j h, j kl (u ) nl uk )ds

) j [C j(,11 11 (u )

eff , j 11

eff , j ) j C j(,22 22 (u ) 22

=

) T + C j(,12 (

h, j j kl (u ) nl uk

(

C2

j eff , j 11 (u ) 22

j h, j kl (u ) nl uk )ds+ j eff , j 22 (u ) 11

T+

T ]dS, (18)

^1j (r , ) = r

) = Kr 1

3 i ^ j 2 1(

1 +i 2

j 1(

) + Kr

1 2

i

¯ 1j ( ),

2+ i

j 1 (

) + Kr 1

2 i

¯1j ( ),

j ), u^1 (r , ) = r

1 2 i

^ j ( ), j = I , II . 1

I , II 2,12 I , II 2,22

Kres = +

D

C2

h, j ( kl nl u1,I ,kj

j h, j 1, kl nl uk )ds j j ( ) n l 1, k ( ))d 1, kl

( 1,j kl ( ) nl 1,j k ( ) ( h, j n u I , II C2 kl l 1, k ( 1,j kl ( ) nl 1,j k ( )

j h, j 1, kl nl uk )ds j j ( ) n l 1, k ( ))d 1, kl

Eq. (24) shows that and the elements of well, are real or pure imaginary quantities. Choosing

w2I = [1 0]T , (19)

I , II 2,11 I , II 2,12

= p2I , II

= i(P¯ I , II

I , II 2,11

0

+ p3I , II

I , II 3,11

0

,

(25) (26)

P I , II )[ wiI , II 0 ], i = 2, 3,

P I , II = LI , II µI , II (LI , II ) 1,

µI , II = diag(µ1I , II , µ2I , II ).

(27)

Eqs. (24)–(27) result in:

,

TI

= 2 Im(µ1I + µ 2I ) p2I

T II = 2 Im(µ1II + µ2II ) p2II

2 Im(µ1I . µ2I ) p3I , 2 Im(µ1II . µ2II ) p3II ,

(28)

where

where the subscripts “appl” and “res” refer to applied mechanical loading and residual stress loading, respectively. The expression in Eq. (20)2 corresponds to the formula derived for isotropic bimaterial notches in [23]. In the next, a reference length l [5] is introduced which defines the new complex stress intensity factor K = Kli whose physical dimension does not depend on ε and the local phase angle arg(K ) does not depend on the used length unit. The crucial parameter is the local phase angle defined by

H11 H22

as

w3I = [0 1]T

I , II = T 0

I , II i

(20)

= arg(K1 + i H11 H22 K2) = arg

w3I , II

where the matrix B is defined in Appendix A and

j = I , II ,

1,12 (l , 0) , 1,22 (l , 0)

(24)

w2I , II ,

(BI + B I )(p2I w2I + p3I w3I ) = (BII + B II )(p2II w2II + p3II w3II ),

+

j j 1, kl ( ) nl 1, k ( ))d

p3I , II

and taking into account the stress distribution and the continuity of the displacement at the interface one obtains

,

eff , j j eff , j j eff , j eff , j ) j ) ) j [C (j,11 T + C (j,12 ( 11 (u1 ) 22 T + 22 (u1 ) 11 T ) + C (j,22 T ]dS 11 (u1 ) 11 22 (u1 ) 22

( 1,j kl ( ) nl 1,j k ( )

p3I , II

= 0 = 2i(I {p2I , II w2I , II } + I {p3I , II w3I , II }). 0 p2I , II ,

The singular solution and the corresponding auxiliary solution are denoted by the subscript 1 separated by comma from coordinate indices, circumflex ^ denotes the appropriate auxiliary solution. If the contour C1 shrinks to the crack tip, then the FE approximation can be replaced by the singular solution and the following formulae are obtained for Kappl and Kres respectively:

K appl =

(23)

to the right-hand side of Eq. (A5), where are arbitrary constants and w2I , w3I and w2II , w3II respectively, are pairs of linearly independent vectors. The stress free conditions along the crack faces and the continuity conditions along the interface require

p2I , II ,

If proper auxiliary solutions ij (u ) , u are chosen then Eq. (18) provides sufficient information for determining the amplitude for each term in the asymptotic crack-tip fields. Consider first the singular term and the appropriate auxiliary solution, see Appendix A:

) = Kr

(22)

i z I (LI ) 1 (p2I w2I + p3I w3I ), i z II (LII ) 1 (p2II w2II + p3II w3II )

T)

j = I , II .

j 1 (r , u1j (r ,

l , l

where the local phase angles , are associated to two different reference lengths l′ and l respectively. In contrast to the evaluation of SIFs for an interface crack between dissimilar anisotropic materials, studies dealing with the determination of T-stress for an interface crack in anisotropic bimaterials are very scarce. We are aware of the studies [11] and [24]. The authors of [11] derived a relation between the J-based mutual integral (the M-integral) and the T-stress and illustrated it for an interface crack in an infinite anisotropic bimaterial. In the paper [11] the M-integral approach implemented in [25] for obtaining the T-stress for isotropic bimaterial interface cracks using boundary element method was extended to generally anisotropic case and applied to a centre cracked bimaterial plate and centrally cracked interface disc. To our knowledge there are no studies concerning the evaluation of the T-stresses for an interface crack in a dissimilar orthotropic bodies subjected to combination of applied and residual stresses. In this paper, the T-stress is included by adding the terms

where C2 is a remote path whereas C1 is a path very close to the crack tip and the apparent elastic stiffnesses under plane strain conditions ) C j(, ab are defined as follows ) C j(,11 =

+ ln

p2II = µ1II + µ2II

p2I

II I s22 s11 Im

+ p3I

II I s22 s11 Im

µ1II µ2II

µ1II + µ2II µ1II µ2II

I Im(µ1I + µ2I ) + s11II s22 Im(µ1II µ2II ) Im I Im(µ1I µ2I ) + s11II s22 Im(µ1II µ2II ) Im

p3II =

(21)

II I p2I s22 s11 Im

where H11, H22 are the components of the Hermitian bimaterial matrix H, see Appendix A. The scaling factor H11 H22 was introduced in order to maintain the phase shift rule

II I + p3I s22 s11 Im

5

1 µ1II µ2II 1 µ1II µ2II

2 × det(B II + B¯ II )

2 det(B II + B¯ II )

Im(µ1I + µ 2I )

1 µ1I µ2I

µ1I + µ2I µ1I µ2I

,

×

I s11II s22 Im(µ1II + µ 2II ) Im

I Im(µ1I µ2I ) + s11II s22 Im(µ1II + µ2II ) Im

1 µ1I µ2I

µ1I + µ2I µ1I µ2I

.

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

The corresponding displacements are derived from Eq. (A1)1. After some algebraic manipulations they can be reduced to the following form:

u 2I

= r [p2I 2I ( u 2II = r [p2II 2II (

)+ )+

p3I 3I ( p3II 3II

i i

= [(BI + B I ) cos

= [(BII + B II ) cos wiII

)], ( )],

II

I

=

(BII

+ ( µI BI + µ I B I ) sin ] wiI , i = 2, 3 + ( µII BII + µ II B II ) sin ] wiII , i = 2, 3 + B II ) 1 (BI + B I ) wiI , i = 2, 3.

Remark: For clarity, it is convenient to rewrite the expressions in Eq. (29) to include the T-stress:

(29)

Fig. 3. Comparison of the stress and displacement field in the vicinity of the interface crack tip (on the circular path of r = 0.004 mm) obtained from the asymptotic

and the complete numerical solution for the bi-material combination where 104N/m, (b) pure thermal loading with ΔT = −100 K.

I

=

6

II

= 0.1, ρI = ρII = 2 and α = −0.8: (a) pure mechanical loading with P = 2.5

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

u 2I = r [p2I u 2II = r [p2II where

j 2 (

I 2 ( II 2 (

)=

) + p3I

I I I 3 ( )] = T r 2 ( p3II 3II ( )] = T II r 2II

),

)+

( ),

p2j 2j ( ) + p3j 3j ( ) j 2 Im(µ1 + µ2j ) p2j 2 Im(µ1j . µ2j ) p3j

can be written as follows: (30)

I , II i, kl

, j = I , II .

=r

2(

I , II i, kl (

)+

I , II i, kl (

)), i = 2, 3

(31)

where

To determine the T-stress using the Betti theorem, an auxiliary 0 has to be known. It elastic field with the stress singularity ~r 2 as r

Fig. 4. Comparison of the stress and displacement field in the vicinity of the interface crack tip (on the circular path of r = 0.004 mm) obtained from the asymptotic

and the complete numerical solution for the bi-material combination where m, (b) pure thermal loading with ΔT = −100 K.

I

=

II

7

= 0.1, ρI = ρII = 2 and α = 0.8: (a) pure mechanical loading with P = 2.5 104 N/

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al. I , II i (

I

I , II

µI , II (RI , II ( )e i

= LI , II

( )) 2

(RI , II ( )ei

I , II

( )) 2

1

LI , II wiI , II

,

1

LI , II wiI , II

i = 2, 3 The auxiliary displacement fields follow as

uiI , II

=r

1 ( I , II i

I , II

( )+

i

( )), i = 2, 3,

(32)

I , II

(33)

where i

I , II

( ) = AI , II (RI , II ( )e i

( )) 1

(LI , II ) 1wiI , II , i = 2, 3.

The functions RI , II ( ), I , II ( ) are given at the end of Appendix A. Finally, the values of the above mentioned parameters p2I , p3I can be found using the Betti reciprocal theorem written in the form of path independent integral: C2

piI =

(

j i, kl (

[

h, j j kl nl ui, k

) nl

j i, k (

j h, j i, kl nl uk ]ds j i, kl (

)

) nl

j i, k (

))d

, i = 2, 3, j = I , II . (34)

Determination of the T-stress for a thermal problem proceeds in a similar way as the determination of SIF Kres in Eq. (20)2. Namely, the area integral must be included as [ h, j n u j C2 kl l i, k j ( i, kl ( ) nl i,jk ( )

piI = +

D

j h, j i, kl nl uk ]ds j j i, kl ( ) nl i, k ( ))d

( ij, kl ( ) nl i,jk ( )

j j i, kl ( ) nl i, k ( ))d

Note that K appl is a linear function of the line load P [N/m], see Fig. 1, and Kres is a linear function of the thermo-elastic coefficients mismatch and the temperature change ΔT. It allows to define normalized stress intensity factors kappl and kres , which depend only on the geometry of the specimen and boundary conditions, as follows:

,

i = 2, 3, j = I , II .

K appl = kappl P ,

(35) However, in contrast to Eq. (20)2, the integral over the domain D is highly singular because klj (ui )~r 2 . It can be evaluated in terms of Hadamard’s concept. In accord with Hadamard’s concept the finite part of the following divergent integral (denoted by F.P.) is defined as

F.P.

dr

r2 r1

r

r1

= lim c

r1

dr

r2 c

r

r1

+ log(c

r1) = log(r2

Kres = kres

where r1, r2 are the radii of the integration contours C1 and C2, rer 2 dr 0 r

= lim F.P . r1 0

r 2 dr r1 r r1

eff 11 1 1 tI EI

T

EI 1 1 2 I + 2tII EII EI 2

EII 2 II EII2

, (37)

where the brackets . denote a jump through the interface between the eff eff , II eff , I = 11 layers I and II 11 11 . Observe that the second term on the right-hand side of Eq. (37)2 has also the physical dimension of a line load, i.e. N/m. In order to measure a relative importance of thermal stresses to mechanical ones the following dimensionless parameter ξ is introduced

r1),

r1 < c < r2, spectively. Hence F.P .

(36)

Ktot = K appl + Kres .

+

eff , j j eff , j j eff , j eff , j ) j ) ) j [C (j,11 T + C (j,12 ( 11 (ui ) 22 T + 22 (ui ) 11 T ) + C (j,22 T ]dS 11 (ui ) 11 22 (ui ) 22

II

, ρ II defined in Appendix characterized by the parameters α, , ρ I, A. The Dundurs parameter α was changed in the interval {−0.9…0.9} – I with the step 0.1 and other parameters were chosen as follows = I {0.1; 0.5; 1}, ρ = {0; 1; 2; 3; 4; 5}. Poisson’s ratios were considered to I II = = for simplicity. take a single value ν I = ν′ I = 0.3 and The asymptotic stresses and displacements calculated along the circular path with radius r = 0.004 mm encircling the crack tip together with results obtained by FEM are shown in Fig. 3 and Fig. 4. The plots show a very good agreement of the asymptotic solution with the complete FEM solution obtained using a very fine mesh, which also demonstrates the accuracy of the SIF and the T-stress calculations. Fig. 5 shows the phase angle for the case of pure mechanical loading (a) and pure residual stress loading (b) as a function of the Dundurs parameter α for several combinations of the parameter 1 1 = 2 (2~ s12 + ~ s66 )(~ s11~ s22) 2 and the ratio . It is interesting to notice in Fig. 5(a) that for α≈0.1 and a specified value of the ratio the phase angle in case of pure mechanical loading does not depend on the parameter ρ. A similar behaviour can be observed in case of pure thermal loading, see Fig. 5(b). If the mechanical and residual loading are superimposed, the total intensity factor is the sum of the two contributions

)

= log(r2).

Fracture-mechanical parameters Kappl, Kres, Tappl and Tres were calculated for a wide range of combinations of orthotropic materials

Fig. 5. The phase angle of SIF K as functions of the Dundurs parameter α and parameters ρ and 8

for: (a) pure mechanical loading; (b) for pure thermal loading.

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

=±

I , II I , II Tappl = tappl P,

Gres |K | = ± res Gappl |K appl |

|k | = ± res |kappl |

1 P

I , II I , II Tres = tres

eff 11 1 tI

1 E'I

2

EI

ÂI +

2 E'I

T 1 2tII

EII

1

2 E ' II

E ' II

= |Ktot| exp(i

appl )

+ |Kres| exp(i

Consider a perturbation of the domain Ω with interface crack between materials I and II as shown in Fig. 9. The perturbation is a deflected crack extension of length Δa into the material I or II making an angle ω with the x1-axis, or a crack extension along the interface (ω = 0), respectively. Detailed FEM calculations for the specimen specified in Fig. 1 revealed that the sum of first two asymptotic terms differs by less than 5% from the complete (FEM) solution in the region around the crack tip the size of which is approximately 30 μm. The validity of the applied perturbation approach to determination of the incremental ERR requires the length of finite crack increment to be at least three times smaller, i.e. with the upper bound about 10 μm. The 1, where L is small perturbation parameter τ is defined as = a L the characteristic length of the domain Ω, e.g. the crack half-length a or the layer thickness t. In the next the characteristic length L is selected to be the same as the reference length l introduced in Eq. (2) and both are equal to the crack half-length a. It means that the perturbation parameter τ is approximately 1/200. The change of the potential energy δW between the unperturbed state u0 (without the crack extension) and perturbed state uτ (with the small finite crack extension) is given by the path independent integral:

The amplitude of the total SIF and the total phase angle respectively follow as

tan(

tot )

=

1

2 cos(

appl )

res

sin(

appl )

sin(

res )

cos(

appl )

cos(

res )

.

+

(41)

4. Determination of the incremental ERR for a finite crack increment

H11 tan( ). H22

|Ktot| = |K appl |

.

res )

where the phase angle appl and/or res respectively are connected to the phase angle appl and/or res respectively specified in Eq. (21) by a simple relation

tan( ) =

EII 2 II EII2

Fig. 8 shows the normalised T-stress for the case of pure mechanical loading (a) and the normalized T-stress for pure residual stress loading (b) as a function of the Dundurs parameter α for several values of the ratio and the parameter ρ = 2. Black lines refer to the upper material I, red lines refer to the lower material II, see Fig. 1. As expected, the Tstress has the same value in both materials for zero elastic mismatch, i.e. for the Dundurs parameter α = 0 for both cases of loading. Observe that for given applied loading and residual stress loading respectively the corresponding T-stresses have opposite signs. Compare also with results in Figs. 3 and 4.

(38)

(39)

tot ),

T

EI 2 1 1 I + 2tII EII EI 2

2

ÂII

which is obtained by normalizing the strain energy release rate of the interface crack under residual stresses Gres with respect to the strain energy release rate of the interface crack under given applied loading, Gappl [3]. By a suitable choice of the temperature change ΔT, the thermoelastic coefficients mismatch, and the line load P is it always possible to let the parameter ξ change within the interval 〈−1;1〉. The total SIF can be written then as

Ktot = K appl + Kres = |K appl| exp(i

eff 11 1 1 tI EI

2,

(40)

The total phase angle as a function of the residual stress magnitude parameter ξ is shown in Fig. 6 for several values of the Dundurs parameter α. Observe that the total phase angle changes sign for ξ > −0.9 depending upon the Dundurs parameter α and the ratio . The amplitude of the total SIF as a function of the residual stress magnitude parameter ξ is shown in Fig. 7 for several values of the Dundurs parameter α. I , II I , II and tres Similarly to SIFs, also normalised T-stresses tappl can be introduced as follows:

Fig. 6. The total phase angle as a function of the residual stress magnitude parameter ξ for several values of the Dundurs parameter α and two values of parameter : (a) = 0.1, (b) = 1.

9

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

Fig. 7. The amplitude of the total SIF as a function of the residual stress magnitude parameter ξ for several values of the Dundurs parameter α and two values of parameter : (a) = 0.1, (b) = 1.

Fig. 8. Normalized T-stress for the case of: (a) pure mechanical loading tappl and (b) pure residual stress loading tres, both as a function of the Dundurs parameter α for several values of the ratio and the parameter ρ = 2.

Ωτ

y2

Ωin

x2 MI MII

MI

Ωin

x1 ∆a

MII

MI MII

ω

y1

∆a=1

τ=∆a/L yi = xi / τ

Fig. 9. Outer and inner domain used in the matched asymptotic analysis and coordinate systems of the outer and inner domain.

10

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

W = W0

1 2

W =

(

0 ij (u ) ni uj

ij (u

) ni uj0) ds =

1 I (u0, u ), 2 (42)

where I (u0, u ) denotes the two states path independent integral [2,6,19,26]. Matched asymptotic procedure is used to derive the change of the potential energy induced by a finite crack increment which introduces a perturbation into the crack tip field [2,6]. The following result is obtained for the incremental energy release rate G , for details see Appendix B:

=

W0

W

G( ) =

a

=

1 I (u0 ( , ), v ( , )) 2L

j j (Tappl + Tres ) |Ktot |2 max(A ( ), 0) + L |Ktot |

+

1 2

j j 2 1 (Tappl + Tres ) D ( ) , j = I , II , 2 |Ktot |2

[B ( ) + C ( )]

(43)

where ω stands for the angle of crack deflection, A(ω), B(ω), C(ω), and D(ω) are complicated functions of ω calculated in the inner domain and given in the Appendix B. The functions A(ω) [m/Pa], B(ω) + C(ω) [m3/ 2 /Pa], and D(ω) [m2/Pa] are shown in Fig. 10 for several values of the dimensionless parameter ξ, which measures a relative importance of residual stresses with respect to applied loading, and specified values of the parameters α, λ and ρ. Making use of the relations – and taking I II j j + Tres ) |Ktot | of the functions B = , the coefficients (Tappl (ω) + C(ω) and D(ω) in the expression for the incremental energy release rate can be rewritten as follows j j Tres + Tappl = |Ktot | |kappl |

j tres

1

|kappl | |kres |

j + tappl

2 cos( res

2 appl ) +

, j = I , II .

(44)

Notice that the courses of A(ω) and B(ω) + C(ω) change rapidly with the residual stress magnitude parameter ξ within the interval −1 ≤ ξ ≤ 0 due to the dramatic change of the total phase angle for this range of ξ, see Fig. 6 and the boundary condition on the local domain defined in Eq. (B8). However, for 0 < ξ ≤ 1 the courses of A(ω) and B (ω) + C(ω) change only slightly with increasing ξ. In contrast, D(ω) does not depend on the residual stress magnitude parameter ξ, which is a consequence of the definition of the parameter K22 in Eq. (B10) from which the function D(ω) is inferred. Fig. 11 shows angular distribution of the normalized ERR G (ω)/ G (0) calculated using the relation for the combined loading characterized by several values of the residual stress magnitude parameter ξ, the crack kink length characterized by the perturbation parameter τ = Δa/2a = 1/100, and for two values of the Dundurs parameter α = −0.6 and α = 0.6. The values of orthotropic parameters = 0.1 and ρ = 2 were used in both layers I and II. The pure mechanical loading corresponds to ξ = 0. Observe that the first term in Eq. (43) deduced from the singular field does not depend on the crack kink length and does not exhibit any oscillatory character, contrary to the previous results in literature. Black lines show the normalized ERR G (ω)/G (0) calculated using only the first term in Eq. (43) which corresponds to the situation of kink length approaching to zero, i.e. τ → 0. Red lines show the normalized ERR G (ω)/G (0) calculated with higher order terms in Eq. (43) involved. Bullets on diagrams mark a maximum of G (ω)/G (0) and indicate the direction of the crack deflection which is plotted versus the residual stress magnitude parameter ξ in Fig. 12. 5. Discussion

Fig. 10. The courses of: (a) A(ω), (b) B(ω) + C(ω), and (c) D(ω) for several values of residual stress magnitude parameter ξ and specified values of the parameters α, λ and ρ.

Before discussing the obtained results, it should be noted that the energy release rate depends on the overall anisotropic elastic response of the layers and the fracture toughness, in general, is also orientation dependent and anisotropic, which makes the study of crack kinking very complex. The fracture anisotropy is usually larger than the elastic 11

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

Fig. 11. Angular distribution of the normalized ERR G (ω)/G (0). Bullets mark maxima of the normalized ERR G (ω)/G (0) and indicate crack deflection directions for α = ± 0.6, = 0.1, ρ = 2, τ = 1/100 and various loading parameters ξ: (a) ξ = −1, (b) ξ = −0.5, (c) ξ = 0, (d) ξ = 0.5, (e) ξ = 1.

anisotropy, so that the direction of the crack extension will mainly be determined by planes having lower fracture toughness than other planes of the solid. For example, crack in a fibre reinforced composite material is usually expected to grow parallel to the stiffer material direction. Obviously, only experiments can answer the question concerning the prediction of the fracture path in anisotropic bimaterial solids. Nevertheless, the goal of this study is not to predict the fracture path, a focus is rather on the estimate of tendency of a crack to deflect out of the interface and to quantify the effect of a thermal loading contribution. It is important to note that the calculations in the preceding section concerning the determination of the incremental ERR =~ s11 ~ s22 = 0.1 which means that the material were performed for direction parallel to the interface is 10 times stiffer than the material direction perpendicular to it. Physically, the stiffer reinforcing fibres are considered parallel to the interface. Hence, it is reasonable to expect that the crack remains after kinking out of the interface in the brittle isotropic matrix of the composite layer. Since only short crack kinks are considered (about 4–10 μm long), the simplifying assumption that the fracture toughness for crack kinking is uniform seems to be acceptable. A necessary condition for a crack to deflect out of the interface is based upon the competition between the crack deflection and crack propagation along the interface [1,6,27–32]

G( ) > G (0)

d i(

)

G( )

d

presented in Figs. 3b and 4b that the residual stresses develop a singular field of the parent crack and, consequently, modify the resulting stress intensity factor under combined thermal and mechanical loading. A variable contribution of the residual stresses to the normalized ERR is characterized by the residual stress magnitude parameter ξ. Note that negative values of ξ correspond to the situation when the residual SIF has a negative real part, i.e. the parent interfacial crack closes. Positive values of ξ then tell us that the residual stresses contribute to the opening of the parent crack. The angular distribution of normalized ERR calculated for the Dundurs parameter α = −0.6 (i.e. the upper material I is more compliant) shows that the kink angle ω closely relates to the total phase angle and this relation is one-to-one. For the residual stress magnitude parameter ξ = −1 is the total phase angle positive, see Fig. 6, and correspondingly, the interface crack tends to kink into the lower material II, of course under the condition, that the criterion is fulfilled. However, when ξ = −0.5, the total phase angle is negative and the interface crack prefers to kink into the upper material I. It is also seen that the ERR is larger when the cracks kinks into the more compliant material I. Different results are obtained for the Dundurs parameter α = 0.6 (i.e. the lower material II is more compliant). In this case, the interface crack also tends to kink into the lower material II for the residual stress magnitude parameter ξ = −1, however both the predicted kink angle size and the ERR are larger in comparison with the previous results obtained for α = −0.6. For ξ = −0.5 and larger, the crack kink remains in the lower material II very close to the interface. Summarizing, if only the first term in Eq. (43) is considered, the criterion still preserves its mathematically local form, i.e. it does not explicitly contain terms which relates to the geometry of the whole structure and loading. The non-local character of the deflection criterion is manifested only by the dependence of the total phase angle on the global distribution of the residual stresses. Consider now that higher order terms in Eq. (43) are involved. In

(45)

where i ( ) is the interface toughness at a phase angle of loading, ψ, Γd is the fracture toughness of the material of the intact layer towards which the kink is directed, G (0) is the ERR for the interface crack. Consider first the normalized ERR G (ω)/G (0) calculated using only the first term in Eq. (43) which does not depend on the crack extension length, see dotted lines in Fig. 11 and dashed lines in Fig. 12, and assume that the K-dominance holds. It is seen according to the results 12

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

Fig. 12. (a)–(c) Dependence of the crack kink angle on the loading parameter ξ and Dundurs parameter α (upon consideration of just singular term of ERR G – dashed curve and full expansion of ERR G – solid line curve); (d)–(f) dependence of the normalized ERR G (ωk)/G (0) (calculated for crack kink angle ωk) on the parameter ξ and Dundurs parameter α.

this case, the criterion depends on the loading parameter j j (Tres + Tappl ) |Ktot | specified in Eq. (44), thus it loses its mathematically local form and evolves to a non-local type of criterion. Let us start with the angular distributions of the normalized ERR calculated for the Dundurs parameter α = −0.6 and a crack extension length characterized by the perturbation parameter τ = 1/500. Fig. 12a shows (full black lines) that there is no distinctive change in the prediction of both the kink angle and the normalized ERR in comparison with the prediction based only on the first term in Eq. (43). Similar predictions are also obtained for a crack extension length characterized by the perturbation parameter τ = 1/200 (see Fig. 12b) and τ = 1/100 (see Fig. 12c). This is due to the low value of the loading parameter j j (Tres + Tappl ) |Ktot | in Eq. (44), see Fig. 13. The angular distributions of the normalized ERR calculated for the Dundurs parameter α = 0.6 indicate that the effect of higher order terms is stronger than for α = −0.6, see Fig. 12. The cause of such behaviour is given by (i) the nature of dependence of the mechanical and residual T-stresses on the Dundurs parameter α, see Fig. 8, with a j j + Tappl ) |Ktot |,see direct impact on values of the loading parameter (Tres Fig. 13 and (ii) the nature of the functions B(ω) + C(ω) and D(ω), see Fig. 10. First it should be noted that for α = 0.6 and ξ = −1, the loading parameter in Eq. (44) takes considerably more negative values

Fig. 13. Loading parameter in Eq. (44) as a function of the residual stress magnitude parameter ξ in the material I and II for α = ± 0.6. 13

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

in both materials than for α = −0.6, see Fig. 13. As a result, the normalized ERR and the absolute value of the kink angle are significantly reduced with respect to the corresponding values obtained using only on the first term in Eq. (43). When the residual stress magnitude parameter ξ is increased to −0.5, the loading parameter j j (Tres + Tappl ) |Ktot | in Eq. (44) takes positive values in both materials, however its value in the lower material II is very small. Also observe that the function B(ω) + C(ω) is negative in the upper material I. Consequently, the effect of higher order terms for ξ = −0.5 is very small and practically independent of the crack extension length characterized by the perturbation parameter τ. An interesting situation occurs for ξ = 0, i.e. for zero residual stresses. While for shorter crack extensions characterized by τ = 1/500, 1/200 the crack kink remains in the lower more compliant material II very close to the interface, as also predicted by considering only the first term in Eq. (43), for a longer crack extension characterized by τ = 1/100 the maximum of the angular distribution of the normalized ERR shifts to the upper stiffer material I at ω = 55°. For higher values of the residual stress magnitude parameter ξ the loading parameter steadily increases and, as a results, the crack kink switches from the lower more compliant material II to the upper stiffer material I also for shorter crack extensions, see Fig. 12b. Apparently, with higher order terms involved, the criterion also depends on the crack extension length since the effects of the Tstresses are then more pronounced.

•

•

6. Summary

• Using the LES formalism, the Betti reciprocal theorem and FEM, the

• •

•

complex SIF and the T-stress of the parent interfacial crack in a symmetric laminate made of three dissimilar orthotropic layers were calculated for combinations of applied four-point bending and residual stresses. The calculations were performed for a wide range of the Dundurs parameter α and several values of orthotropic parameters ρ and defined in Appendix A. The accuracy of the SIF and the T-stress calculations was demonstrated by comparison of the asymptotic solution with the complete FEM solution obtained using a very fine mesh. A relative importance of thermal stresses to mechanical ones was 1, 1 characterized by means of the dimensionless parameter which is obtained by normalizing the strain energy release rate of the interface crack under residual stresses Gres with respect to the strain energy release rate of the interface crack under given applied loading Gappl . The pure mechanical loading corresponds to ξ = 0. The residual stress magnitude parameter ξ allows to analyze a contribution of residual stresses to the phase angle and the amplitude of the resultant SIF of the parent interfacial crack if the

mechanical and residual loading are superimposed. It was found that the phase angle of the resultant SIF changes sign for ξ > −0.9 depending upon the Dundurs parameter α and the ratio . Matched asymptotic procedure was used to derive the change of the potential energy induced by a finite crack increment which introduces a perturbation into the crack tip field. Consequently, the incremental ERR, which does not exhibit any oscillatory character, was evaluated. The technique enabled analyzing a wide range of finite crack increment geometries under various mechanical versus residual loading combinations. The incremental ERR was evaluated up to second order terms thus including an effect of the T-stress of the parent interfacial crack. The competition between the crack deflection and crack extension along the interface was analyzed by virtue of the incremental ERR. If the effect of the T-stress was not considered, it was found that the crack kink angle ω closely relates to the total phase angle, which strongly depends on the relative contribution of residual stresses. If the residual stress magnitude parameter ξ increases, the crack deflection either switches from the stiffer material to the more compliant one or it remains in the more compliant material but the crack kink closely approaches the interface. In this case, the non-local character of the deflection criterion is manifested only by the dependence of the total phase angle on the global distribution of residual stresses. If the size of the region of K-dominance is not by order greater than the crack kink length, higher order terms depending on the T-stress have to be involved in the deflection criterion. In this case, the crack deflection depends on a relative importance of the T-stress due to pure mechanical loading and the Tstress due to pure residual stress loading, and on the Dundurs parameter α. It is important of noting that the deflection criterion with higher order terms involved depends on the loading, thus it loses its mathematically local form and evolves to a non-local type of criterion.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The support of the Czech Science Foundation through a grant number 17-08153S is gratefully acknowledged.

Appendix A Assuming two aligned orthotropic materials, a general solution for the displacements that satisfy the equation of equilibrium and corresponding stress components may be written in terms of two analytic complex functions as

ui = 2Re{Aij f j (zj )},

2i

= 2Re{Lij f j (zj )},

1i

= 2Re{Lij µj f j (zj )},

(A1)

where i, j = 1,2, f j (zj ) are analytic in their arguments, z j = x1 + µj x2 , µj = µj + iµj are distinct complex numbers with positive imaginary parts representing eigenvalues of a material, Aij and Lij are matrices depending on elastic constants

~ ~ s11 µ12 + ~ s12 s11 µ22 + ~ s12 A= ~ , L= ~ ~ s12 µ1 + s22 µ1 s12 µ 2 + ~ s22 µ 2

µ1 1

µ2 , 1

(A2)

sij are defined as where for plane strain conditions the compliances ~ 1 ~ s11 = E

E E2

2

, ~ s12 =

E

+

E

1 , ~ s22 = E

2

E

.

A positive definite Hermitian matrix B and, for the bimaterial case, a positive definite Hermitian matrix H are defined as follows (A3)

B = iAL 1 , H = BI + BII ,

14

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

where the superscripts I and II refer to the materials I and II respectively. According to Suo [8] it holds 1 4

2n

B=

(

1 (~ s11~ s22 ) 2

1 i (~ s11 ~ s22 ) 2 + ~ s12

(

1 i (~ s11 ~ s22 ) 2 + ~ s12

)

1 4

2n

)

1 (~ s11 ~ s22 ) 2

1

H11

, H=

i (H11 H22 ) 2 1

i (H11 H22 ) 2

H22

(A4)

where 1 4

H11 = 2n

1 (~ s11 ~ s22 ) 2

1 4

+ 2n

1 (~ s11 ~ s22 ) 2

I

1 1 (H11 H22 ) 2 = (~ s11~ s22 ) 2 + ~ s12

II

(

n=

1 (~ s11~ s22 ) 2

1 4

+ 2n

1 (~ s11 ~ s22 ) 2

I

1 (~ s11~ s22 ) 2 + ~ s12

II

1 4

, H22 = 2n

1 (1 2

I

+ )

,

1 = 2 (2~ s12 + ~ s66 )(~ s11 ~ s22)

s11 , s22

=

II 1 2,

1 2

).

Two generalized Dundurs bi-materials parameters α, β are defined as 1 (~ s11 ~ s22 ) 2

=

1 (~ s11~ s22 ) 2

II

1 (~ s11 ~ s22 ) 2

1 s11~ s22 ) 2 + (~

II

I f1,1 (z1 )

f1I (z ) =

(

= K

I f1,2 (z 2 )

1 2

1 i, 2

w1 =

= diag 1 (H22 2

)

1

K¯

1 +i z2

(

,

, z2

1

1 T

H11 ) 2

1

(

1 2

)

1

=

1 2

i

1 +i 2

I

1

z I 2 + i (LI ) 1w1 + e

1 2

1 (~ s11 ~ s22 ) 2 + ~ s12

II

(H11 H22 ) 2

= K

II f1,2 (z2 ) 2

1 (~ s11 ~ s22 ) 2 + ~ s12

=

I

+i

+e 1

,

II f1,1 (z1 )

f1II (z ) =

where z 2 + i

I

)

+i

1

z II 2 ln

K¯

2

(

1 2

i

)

1

1

zI 2

1

1 z II 2 + i (LII ) 1H¯ 1Hw1+

i

(LII ) 1H¯ 1Hw¯1

( ) and w 1 1+

1

i

(LI ) 1w¯1

(A5)

is an eigenvector defined as

.

(A6)

The analytic complex functions which correspond to the auxiliary solution are as follows: I

f^1 (z ) = k^1

(

1 2

i

)

1

1 2

zI

II f^1 (z ) = k^1

¯ + e 2 k^1

(

(

^ 1 + e 2 k^¯1 (LI ) 1w

i 1 2

1 2

)

i

+i

)

1

1

1 2

z II z II

(

1 2

+i

)

1

zI

^¯ 1 (LI ) 1w

1 +i 2

(LII ) 1H¯ 1Hw^1 +

i

(LII ) 1H¯ 1Hw^¯1

1 +i 2

(A7)

where the eigenvector w1 is defined by 1 i, 2

w1 =

1 (H11 2

T

1

H22 ) 2

(A8)

The angular functions in Eq. (19)1 then follow from Eqs. (A1), (A5) and (A7) as I 1

( )=

(

1 2

+i

)

1

1

AI e 2 lnR

+i

( ln R + 12 ) II 1

1

× AII e 2 lnR I 1(

+e II 1 (

+e

+i

) = LI e 2

L¯ I e

) = LII e 2

L¯ II e

1 lnR 2

+i

(LII ) 1H¯ 1H + e

( ln R

1 2

)(

sin

1 2

2

)

×

+i 2

1

1

A¯ I e 2 lnR +

1

A¯ II e 2 lnR +

+i

( ln R

( ln R

sin

+ µ cos ) (LII ) 1H¯ 1H +

+i

1 2

( ln R+ 12 ) (

)( sin

+i

( ln R

1 2

)

1 2

)

(L¯ I )

1

w1

(LII ) 1H 1H¯ w1

(A9)

+ µ cos ) (LI ) 1+

+ µ¯ cos ) (L¯ I )

1 lnR 2

(cos

(

+e

sin

+i

+i

( )=

1

( ln R + 12 ) (

1 lnR + 2

1 lnR + 2

where RI , II ( ) =

( ln R+ 12 )

(LI )

1

w1,

+ µ¯ cos ) (L¯ II ) 1H 1H¯ w1,

(A10)

+ Re(µI , II ) sin )2 + (Im(µI , II ) sin )2 , θ is the angle of the polar coordinate system with origin in the crack tip.

15

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

= 0,

0 for I 1 cos + Re(µ ) sin Im(µ I ) sin

cot

( )=

II 1 cos + Re(µ ) sin Im(µ II ) sin

cot

for

-

(0, ),

for

(

, 0).

The angular functions of the auxiliary solution in Eq. (19)2 are derived in a similar fashion. Appendix B The change of the potential energy δW between the unperturbed state u0 and perturbed state uτ is given by the relation. A second scale to the problem can be introduced, represented by the scaled-up coordinates

yi =

xi

x1 x2 , ,

, or (y1 , y2 ) =

(B1)

which provides a zoomed-in view into the perturbed region surrounding the crack. The displacement u of the perturbed elasticity problem due to the crack extension can be expressed in terms of the regular coordinate x and the scaled-up coordinate y as u (x ) = u ( y ) = v (y ) . The unperturbed state u0 has the following expansion form τ

u0 (r , ) = Kr 1

2+ i

1(

) + Kr 1

2 i

¯1 ( ) + Tr

2(

(B2)

) + ...higher order terms,

or in terms of the local polar coordinate ρ = r/τ

u0 ( , ) = K 1 2

=K

1 2+i

()

i

1 2+i

1 2+i

l

1(

1(

1 2 i

)+K 1 2

)+K

()

i

l

1 2 i

¯1 ( ) + T

2(

) + ..higher order terms=

1 2 i

¯1 ( ) + T

2(

) + ..higher order terms,

(B3)

where K, K = are the complex stress intensity factors and T is the T-stress which takes the different value in the materials I and II. Consider now the asymptotic outer expansion for uτ and the inner expansion for vτ

Kli

u (r , ) = u0 (r , ) + f1 ( ) r

^ ( ) + f ( )r 1 1

1 2 i

where limfi + 1 ( ) fi ( ) = 0,

i = 1, 2, ... and

0

1 2+i

i(

^¯ ( ) + f ( ) r 1 ^ ( ). .., 1 2 2

(B4)

) form a set of linearly independent basis functions.

Observe that 1 ( ) is the auxiliary solution which is also used in the calculation of the complex stress intensity factor. The inner asymptotic expansion is possible to write in the following form: 1 2

v ( , ) = F1 ( ) 1 2

+ F1 ( )

()

i

()

i

1(

l

) + K^11

¯1 ( ) + K^11

l

+ F2 ( )[ 0

i(

)+

()

i

l

l

i

^¯ ( ) + K^12 1

^ ( ) + K^12 1

1^ ( 2

Kij

j

j

1^ ( 2

1^ ( 2

) + ... +

) + ... +

) + ...] + ...,

(B5)

1 +i 2

, F2 ( ) = T and the expressions in brackets form a set of linearly independent basis

i = 1, 2, ..., F1 ( ) = K

where limFi + 1 ( ) Fi ( ) = 0,

i

()

) + K^21

2(

functions vi(ρ,θ)

vi ( , ) =

1 2

1 2

( ) + ...,

(B6)

j=1

where δi are generally complex eigenvalues with positive real part pertaining to the Williams-like asymptotic expansion and responding eigenfunctions. This allows writing the inner asymptotic expansions in the form

v ( , )=

[Fi ( ) vi ( , ) + Fi ( ) vi ( , )].

i(

),

i(

) are cor-

(B7)

i=1

The outer expansion is valid in the whole domain except near the former main crack tip where the geometry is perturbed. The inner expansion 0 becomes the unbounded ‘inner’ domain in spanned by the stretched describes the near field in the stretched ( ×1 ) domain which as 0, variables y1 and y2. The inner expansion is valid for y→∞. The coefficients Kij are computed on the inner domain in which is unbounded for but in the employed numerical model, particularly developed using the FEM, the in is approximated by a circular region with a radius Rin much larger than the deflected crack extension Δa. Specifically, the coefficients K11 and K12 are calculated in the inner domain whose remote boundary ∂Ωin is subjected to the following boundary condition

u1 |

in

= exp(i

=

Ktot

1 2

|Ktot |

tot )

1 2

()

()

i

l

i

l

1(

1(

)+

Ktot

) + exp( i

( ) ¯ ( )= ) () ¯()

1 2

|Ktot | tot

i

1

l

1 2

i

l

(B8)

1

The coefficients K11 and K12 are then calculated as

K^11 =

I (v1h, I(

1 2(

l)

1 2( i

^ ( ), 1

l)i

1( 1 2(

)) l)i

1(

))

, K^12 =

I (v1h, 2 ( )) , I ( 1 ^2 ( ), 2 ( ))

(B9)

where is FE approximation to v1 and the two states path independent integral I (.,.) is specified in Eq. (42). The coefficients K21, K22 are calculated in the inner domain whose remote boundary ∂Ωin is subjected to the boundary condition u2 | in = 2 ( ) and they follow as

v1h

16

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al.

K^21 =

I (v2h, 1 2(

I(

l)

1 2( i

l)i

1( 1 2(

^ ( ), 1

)) l)i

1(

))

I (v2h, 2 ( )) . I ( 1 ^2 ( ), 2 ( ))

, K^22 =

(B10)

The higher order coefficients Kij can be calculated accordingly. Remark. It is important to notice that the coefficients Kij in the inner asymptotic expansion (B5) depend only on the local geometry and they are independent of the geometry of the whole structure and the applied loads. They are computed once for all use with different loading conditions. As a result, a number of necessary computations of the incremental energy release rate (ERR) within a wide range of loading conditions is significantly reduced in comparison with a standard finite element analysis. Another advantage is that the semi-analytical expression for the incremental ERR, cf. Eq. (B14), is obtained which makes it easier to study the effect of loading parameters such as the T-stress. Making use the property of the orthogonality with respect to the path independent integral I:

I(

i

i(

),

j(

j

0 for i j const 0 for i = j ,

)) =

(B11)

one obtains to the second order the change of the potential energy δW from Eqs. (42), (B3) and (B5) as follows 1 I (u0, 2

W=

1 I 2

u)=

1 +i 2

K K

1 2

i

1 2(

[

[

21

1 2(

l)i

1 +i 2

1(

)+K

^ 2 ( l ) i ^¯ ( ) + K 12 1 ^ 2 ( l) i ^ ( ) + K

1

1(

i

l)

(K

) + K^11 ¯1 ( ) + K^11

l)i

1 2(

T [ 2 ( ) + K^21 1 2 ( l) i K K Re[K^11 I ( 1 2 ( l) i 1 ( 3 T 2 Re[K K^ i I ( 1 2 ( l) i

=

1 I 2

(u0 ( , ), v ( , )) =

1

1^ ( 2 1^ ( 2

12

1

^¯ ( ) + K^21 1

1 2(

l)

i

1 2

i

¯1 ( ) + T

2(

) + ...,

) + ...]+ ) + ...]+

l) ^1 ( ) + K^22 1 ^2 ( ) + ...] + ...)= 3 ), 1 2 ( l) i ^1 ( ))] T 2 K^12 Re(K i ) I ( 2 ( ), 1 ^2 ( )) 1 2 2^ 1 2 ( l ) i ^ ( ))] T K22 I ( 2 ( ), 1 ^2 ( )). 1 ( ), 1 2 1 2(

i

(B12)

Using the definitions of the coefficients Kij in (B9) and (B10) and considering the total SIF specified in Eq. (39), the result for the change of the potential energy δW can be finally written in the following form:

|Ktot |2 { Re[I (v1h,

W=

T |Ktot |

3 2

i

Re[

exp(i

1 2( tot ) I

l)i

1(

(v2h,

T

))]

1 2(

|Ktot |

l)i

1(

3 2

i

Re[

exp(i

T2

1 2I 2 |Ktot |2

))]

tot )] I

(v2h,

(v1h, 2(

2(

))

}

)) .

(B13)

Correspondingly, the incremental ERR G ( ) reads

G( ) =

W0

W

1 |K |2 T I (u0 ( , ), v ( , )) = tot max(A ( ), 0) + 2L L |Ktot |

=

a

1 2

(B ( ) + C ( )) +

1 T2 D( ) , 2 |Ktot |2

(B14)

where

Re[I (v1h,

A( ) = B( ) = C( ) =

Re[

Re[

i

i

exp(i

exp(i

D( ) =

1 2(

l)i

h tot )] I (v1 ,

h tot ) I (v2 ,

I (v2h,

1 2( 2(

)),

1(

))],

2(

l)i m2 Pa

)), 1(

m Pa

,

m3 2 Pa

))],

, m3 2 Pa

,

,

and max(A ( ), 0) selects such a range of crack deflection angles that the first term in the incremental ERR G ( ) deduced from the singular field is positive.

[9] L. Banks-Sills and V. Boniface, in: T.-. Chuang and J.W. Rudnicki (Eds.), Springer Netherlands, 2002, pp. 183-204. [10] H.G. Beom, C.B. Cui, H.S. Jang, Dependence of stress intensity factors on elastic constants for cracks in an orthotropic bimaterial with a thin film, Int. J. Solids Struct. 49 (2012) 3461–3471. [11] J.H. Kim, H.J. Moon, Y.Y. Earmme, Inplane and antiplane T-stresses for an interface crack in anisotropic bimaterial, Mech. Mater. 33 (2001) 21–32. [12] M. Nagai, T. Ikeda, N. Miyazaki, Stress intensity factor analysis of a three-dimensional interface crack between dissimilar anisotropic materials under thermal stress, Eng. Fract. Mech. 91 (2012) 14–36. [13] T. Profant, J. Klusák, O. Ševeček, M. Hrstka, M. Kotoul, An energetic criterion for a micro-crack of finite length initiated in orthotropic bi-material notches, Eng. Fract. Mech. 110 (2013) 396–409. [14] C. Hwu, T.L. Kuo, A unified definition for stress intensity factors of interface corners and cracks, Int. J. Solids Struct. 44 (2007) 6340–6359. [15] X. Deng, Mechanics of debonding and delamination in composites: asymptotic studies, Compos. Eng. 5 (1995) 1299–1315. [16] T.C.T. Ting, Anisotropic Elasticity: Theory and Applications, Oxford University Press, USA, 1996. [17] R. Desmorat, F.A. Leckie, Singularities in bi-materials: parametric study of an isotropic/anisotropic joint, Europ. J. Mech., A/Solids 17 (1998) 33–52. [18] K.C. Shin, W.S. Kim, J.J. Lee, Application of stress intensity to design of anisotropic/

References [1] Ming Yuan He, A.G. Evans, J.W. Hutchinson, Crack deflection at an interface between dissimilar elastic materials: role of residual stresses, Int. J. Solids Struct. 31 (1994) 3443–3455. [2] D. Leguillon, C. Lacroix, E. Martin, Crack deflection by an interface - asymptotics of the residual thermal stresses, Int. J. Solids Structures 38 (2001) 7423–7445. [3] W. Lee, H. Shin, Y. Yoo, Revisited criterion to prevent kinking of a crack out of a bimaterial interface under the presence of in-plane residual stresses - a global approach, Int. J. Eng. Sci. 42 (2004) 257–270. [4] M. He, J.W. Hutchinson, Kinking of a crack out of an interface, Trans. ASME 56 (1989) 270–278. [5] J.R. Rice, Elastic fracture mechanics concepts for interfacial cracks, J. Appl. Mech. 55 (1988) 98–103. [6] D. Leguillon, S. Murer, A criterion for crack kinking out of an interface, Key Eng. Mater. 385–387 (2008) 9–12. [7] ANSYS Inc, ANSYS Release 18.1 Useŕs manual, Swanson Analysis Sys. Inc, Pensylvania, 2016. [8] Z. Suo, Singularities, Interfaces and Cracks in Dissimilar Anisotropic Media, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 427 (1990) 331-358.

17

Theoretical and Applied Fracture Mechanics 105 (2020) 102397

O. Ševeček, et al. isotropic bi-materials with a wedge, Int. J. Solids Struct. 44 (2007) 7748–7766. [19] D. Leguillon, E. Sanchez-Palencia, Computation of Singular Solutions in Elliptic Problems and Elasticity, John Willey & Son, New York, 1987. [20] L. Banks-Sills, A conservative integral for determining stress intensity factors of a bimaterial strip, Int. J. Fract. 86 (1997) 385–398. [21] G.B. Sinclair, M. Okajima, J.H. Griffin, Path independent integrals for computing stress intensity factors at sharp notches in elastic plates, Int. J. Numer. Methods Eng. 20 (1984) 999–1008. [22] P.H. Wen, M.H. Aliabadi, D.P. Rooke, A contour integral for the evaluation of stress intensity factors, Appl. Math. Model. 19 (1995) 450–455. [23] L. Banks-Sills, C. Ishbir, A conservative integral for bimaterial notches subjected to thermal stresses, Int. J. Numer. Methods Eng. 60 (2004) 1075–1102. [24] P.D. Shah, C.L. Tan, X. Wang, Evaluation of T-stress for an interface crack between dissimilar anisotropic materials using the boundary element method, CMES – Comput. Model. Eng. Sci. 13 (2006) 185–197. [25] J. Sladek, V. Sladek, Evaluations of the T-stress for interface cracks by the boundary

element method, Eng. Fract. Mech. 56 (1997) 813–825. [26] O. Ševeček, M. Kotoul, T. Profant, Effect of higher order asymptotic terms on the competition between crack penetration and debond at a bimaterial interface between aligned orthotropic materials, Eng. Fract. Mech. 80 (2012) 28–51. [27] M. He, J.W. Hutchinson, Crack deflection at an interface between dissimilar elastic materials, Int. J. Solids Struct. 25 (1989) 1053–1067. [28] M.Y. He, A. Bartlett, A.G. Evans, J.W. Hutchinson, Kinking of a crack out of an interface: role of in-plane stress, J. Am. Ceram. Soc. 74 (1991) 767–771. [29] D. Leguillon, C. Lacroix, E. Martin, Interface debonding ahead of a primary crack, J. Mech. Phys. Solids 48 (2000) 2137–2161. [30] D. Leguillon, C. Lacroix, E. Martin, Crack deflection by an interface – asymptotics of the residual thermal stresses, Int. J. Solids Structures 38 (2001) 7423–7445. [31] E. Martin, D. Leguillon, C. Lacroix, A revisited criterion for crack deflection at an interface in a brittle bimaterial, Compos. Sci. Technol. 61 (2001) 1671–1679. [32] D. Martínez, V. Gupta, Energy criterion for crack deflection at an interface between two orthotropic media, J. Mech. Phys. Solids 42 (1994) 1247–1271.

18