Fracture parameter determination for the orthotropic interface crack with friction

Engineering Fracture Mechanics 204 (2018) 542–556

Contents lists available at ScienceDirect

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

Fracture parameter determination for the orthotropic interface crack with friction

T

F. Ernesto Penado Department of Mechanical Engineering, Northern Arizona University, 15600 S. McConnell Dr., Flagstaff, AZ 86011-5600, United States

ARTICLE INFO

ABSTRACT

Keywords: Stress intensity factors Frictional interfacial crack Orthotropic materials Composite materials Bonded joints

Many important engineering problems involve cracks at bimaterial interfaces. Some examples of these problems include crack propagation and delamination at composite ply interfaces, the fracture of adhesively bonded joints, the analysis of test specimens to determine fracture properties, the interfacial failure of rock and concrete, and others. Herein, an approach that combines local analytical with global numerical solutions for the analysis of the frictional interface crack is proposed and demonstrated. The method allows the determination of near-tip stresses and stress intensity factors of complex structures with a moderate computational effort and requires no special finite element codes or elements. As an illustration, an infinite plate with a central bimaterial crack under a combination of shear and compressive loads is considered. Results have been presented for a variety of isotropic and orthotropic material combinations and coefficients of friction, and comparisons with available results in the literature show very close agreement. Overall, the method provides a practical way to analyze a variety of engineering structures with frictional interfacial cracks including isotropic and orthotropic materials.

1. Introduction The study of interfacial cracks is of considerable importance in fracture mechanics. Many engineering problems of interest involve cracks at bimaterial interfaces; for example, crack propagation and delamination at composite ply interfaces [1,2], the fracture of adhesively bonded joints [3,4], the analysis of test specimens to determine fracture properties [5,6], the interfacial fracture of rock and concrete [7], and others [8]. The traditional fracture mechanics approach to crack analysis assuming stress free faces leads to the physically inadmissible condition of crack face interpenetration, which occurs even when a crack tends to open. In order to address this ambiguity, Comninou [9,10] proposed a crack contact model that eliminates the oscillatory behavior and results in physically realistic behavior. Furthermore, for structures subjected to combined shear and compressive loads, sizeable contact zones can develop and the effects of friction cannot be neglected [11]. Particularly, frictional crack contact has a considerable effect on the magnitude of the stress intensity factor and the customary inverse square root singularity becomes r , where λ can be smaller or larger than ½, leading to weak or strong singularities. In this case, the energy release rate becomes zero or infinite for the weak or strong singularities, respectively, and requires special considerations if it is to be used as a fracture criterion [11,12]. Additionally, the type of singularity that leads to the correct stress intensity factors must be selected based on knowledge of the direction of slip of the crack surfaces. The use of a fracture criterion requires the determination of fracture parameters such as stresses and stress intensity factors. Methods currently available to determine these parameters include exact analytical solutions and numerical approximations.

E-mail address: [email protected]. https://doi.org/10.1016/j.engfracmech.2018.10.038 Received 17 July 2018; Received in revised form 28 October 2018; Accepted 30 October 2018 Available online 01 November 2018 0013-7944/ © 2018 Elsevier Ltd. All rights reserved.

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Analytical solutions for interfacial cracks with frictional sliding are remarkably difficult to achieve and quickly become intractable. As a consequence, only relatively simple configurations and loadings can be considered. Very few analytical solutions are available for a frictional contact crack, including prominently that of Comninou and Dundurs [13] for an infinite plate with a central crack loaded in shear. They formulated the problem by considering the crack as an array of distributed edge dislocations and eventually reduced the result to a singular integral equation with a Cauchy-type kernel which was solved numerically. They reported a discrepancy with previous results for the frictionless case [14], which they attributed to the slow convergence of some of the hypergeometric functions used. However, a more recent investigation by Dorogoy and Banks-Sills [15] using finite differences has suggested that some of their results for the frictional case contain numerical errors. Overall, for complex, practical applications, numerical methods such as finite elements, finite differences, or boundary elements become necessary. Although stress intensity factors in principle can be determined directly from the stresses and displacements produced by these approximate methods, the process can be computationally intensive and the results may not have the best accuracy. A better option available when considering fracture mechanics problems is to calculate the stress intensity factors indirectly using an energy-based approach such as the Jintegral, M-integral, or crack closure technique. Using these methods, the frictional contact crack has been successfully examined by a number of investigators. Dorogoy and Banks-Sills [15] used an extended J-integral in conjunction with a finite difference technique to investigate the effect of friction on a shear loaded interface crack. The same authors later used the M-integral and finite differences to study the effect of contact and friction on bimaterial Brazilian disk specimens containing a crack [5]. Sun and Qian [16] used a finite element based crack closure method to calculate strain energy release rates from which accurate stress intensity factors for isotropic materials were obtained. They subsequently used this method to analyze problems with strong [11] and weak [12] singularities and establish a fracture criterion based on finite extension strain energy release rate [12]. In this paper, an alternative to commonly used numerical methods for the analysis of the frictional interface crack is proposed and demonstrated. The method consists of combining local analytical with global numerical solutions and offers a number of advantages such as the accurate and direct determination of near-tip stresses and stress intensity factors and a relatively modest computational effort. Analysis of both isotropic and orthotropic materials can be performed while requiring only an unmodified commercially available finite element program. Subsequently, the method is used to study an infinite plate with a central bimaterial crack under a combination of shear and compressive loads. Weak and strong singularities are examined using both singular and non-singular terms to determine the shear stresses ahead of both crack tips and validate the proper set of eigenvalues that must be used. Additionally, stress intensity factors for both isotropic and orthotropic materials for a variety of material combinations and coefficients of friction are presented. Finally, the present results for the stress intensity factor are compared with available solutions in the literature and very close agreement is found. 2. Method of analysis In the present work, a method previously demonstrated for frictionless interfacial cracks in contact [17] is extended to include friction at the crack faces. The method consists of using a numerical procedure such as finite elements to determine the displacements away from the crack tip. Concurrently, an elasticity solution based on eigenfunction expansions of the stresses and displacements is used in a region neighboring the crack tip (local region) where the numerical solution is not valid anymore. Frictional contact at the crack faces is considered in both the numerical solution and the formulation of the eigenfunction expansion. The resulting unknown constants in the eigenfunction expansion are found by matching displacements at a suitable interface contour with the numerical model, leading to a set of linear equations in these arbitrary constants. The method overcomes the limitations of the numerical solution near the crack tip, and provides stress intensity factors at the point of singularity (crack tip) as well as stress values at any point near the singularity. In this work, the numerical results used with the elasticity solution were obtained by means of the finite element method. A graphical representation of the method is given in Fig. 1.

Fig. 1. Schematic illustration of the method of analysis used. 543

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

3. Local eigenfunction solution 3.1. Basic relationships For the local elasticity solution, we assume linear elastic, orthotropic materials. The usual limitations inherent in the assumption of linear elastic material behavior apply, and plasticity effects at the crack tip are not taken into consideration. In rectangular coordinates, with {σ} = stresses and {ε} = strains, the stiffness relation (generalized Hooke’s law) becomes:

E1 (1 E2 ( 12 + 1 E3 ( 13 + = 0 0 0

x y z yz xz xy

23 32 )

E1 ( 21 + E2 (1 23 ) E3 ( 23 + 0 0 0

13 32 ) 12

31 23 )

E1 ( 31 + E2 ( 32 + E3 (1 13 ) 0 0 0

13 31 ) 21

21 32 ) 12 31 ) 12 21 )

0 0 0 G23 0 0

0 0 0 0 G13 0

0 0 0 0 0 G12

x y z yz xz xy

(1)

where Ei, νij and Gij (i, j = 1, 2, 3) are the usual engineering constants, and

=1

12 21

13 31

23 32

12 31 23

(2)

21 13 32

For the special case of isotropic materials, E1 = E2 = E3 = E, ν23 = ν13 = ν12 = ν, and G23 = G13 = G12 = G = E/[2(1+ν)]. Assuming plane strain conditions, εz = 0, the strain-displacement and equilibrium equations become: x

u ; x

= x

xy

+

x

y

v ; y

=

y

xy

= 0;

= 0;

z

y

+

x

yz

w ; y

=

xz

= 0;

y

x

xz

+

=

yz

y

w ; x

xy

=

u v + y x

=0

(3) (4)

Now, replacing Eqs. (3) into (1) and defining cij (i, j = 1,..,6) to be the entries of the stiffness matrix in Eq. (1), we can write the equilibrium Eqs. (4) in terms of displacements: 2u

c11 c66

(

2v x y

+ c12

x2

2u

x y

+

2v

x2

c55

+ c66

)+c

2w

x2

(

2u

y2

2u

+ c22

12 x y

+ c44

2w

y2

2v x y

+

)=0

2v

y2

=0

=0

(5)

It can be seen in Eqs. (5) that the equations for u and v are decoupled from that for w. The first two of (5) correspond to the inplane problem, while the third to the out-of-plane problem. In the present work, only the in-plane problem is considered. Based on the results given by Penado [17] and Ting and Chou [18], we can assume that the solution to the in-plane problem in the vicinity of the crack tip is given by: 1+ p

u=

x

Zp 1+

1+ p

,

v=

Zp 1+

y

p

(6)

p

where Zp = x + μpy, μp represents the eigenvalues of the elasticity constants to be determined from the equilibrium Eqs. (5), ζx and ζy are functions of μp, and λp are eigenvalues to be determined from the boundary conditions near the crack tip. Substituting Eqs. (6) into the first two of Eqs. (5) leads to the following homogeneous system of equations:

c11 + c66 µp2

(c12 + c66) µp

x

(c12 + c66) µp

c66 + c22 µp2

y

{}

= 0 0

(7)

For a non-trivial solution of these equations, their determinant of coefficients must be zero. This leads to the following quartic equation in μp:

c22 c66 µp4 + (c11 c22

2c12 c66

2 c12 ) µp2 + c11 c66 = 0

(8)

It can be shown [18] that μp can only be complex or purely imaginary. Therefore, we have two complex conjugate pairs for μp. Once the μp’s are known, ζx and ζy can be found from (7) up to arbitrary multiplicative constants. Now, substituting (6) into (3) and (1) gives:

c11

x y xy

=

c12

x

+ µp c12

x

+ µp c22

y

+ c66

y

µp c66

x

y

x

Zp p

y

Zp p (9)

xy

In the special case of isotropic materials, the conjugate pair µp = ±i = ± 544

1 is a double root of (8), and the new solution cannot

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

be obtained by simply inserting isotropic material properties into the equations. In this case, the in-plane displacements given by Eqs. (6) become [18]: 1+ p

u=

x

Zp 1+

+ p

1 1+

p

1+ p

d ( Z1 + p) , dµp x

v=

Zp 1+

y

+ p

1 1+

p

d ( Z1 + p) dµp y

(10)

3.2. Asymptotic displacements and stresses The complete asymptotic solutions for displacements and stresses near the crack tip can be obtained by taking linear combinations of Eqs. (10) corresponding to the different eigenvalues μp (summations over k), and λp (summations over m), where k and m are defined below. This leads to the following expressions for the displacements and stresses converted to cylindrical coordinates with x = r cos and y = r sin :

ur(i ) =

M

am m=1

(i ) r

=

r 1 + ( p)m 1 + ( p)m

M a [r ( p)m m=1 m

4

(i ) Akm (Fr(i ) )km ,

M

u (i ) =

k=1

r 1+ ( p)m 1 + ( p )m

am m =1

4 M (i ) A (i ) (Gr(i) )km], = m = 1 am [r ( p)m k = 1 km M 4 (i ) (i ) a [r ( p)m k = 1 Akm (Gr(i ) )km] r = m=1 m

4

(i ) Akm (F (i) )km

(11)

k=1

4 A (i) (G (i ) )km], k = 1 km

(12)

where the superscript i = 1, 2 denotes each material region; k and m are dummy indices of summation associated with the eigenvalues of the problem; M is the number of terms in the truncated infinite series; am’s are multiplicative constants associated with the in-plane eigenvectors; Fij and Gij are defined in Table 1; and (λp)m and Akm will be defined later in relation to the boundary conditions at the crack tip. 3.3. Stress intensity factors The mode I and II stress intensity factors can be calculated from the corresponding singular stresses (those with a negative λp = (λp)1) as follows:

KI = lim(2 r ) r

( p )1

(r , 0),

0

KII = lim(2 r ) r

( p)1

0

r

(r , 0)

(13)

In terms of (λp)1, Ak1, a1, as well as Gθ and Grθ defined in Table 1, we have: 4

KI = (2 )

( p )1 · a 1

4

Ak1 (G |

= 0 )k1,

KII = (2 )

( p )1 · a 1

k=1

Ak1 (Gr |

= 0 )k1

(14)

k=1

4. Boundary conditions for local region 4.1. Boundary conditions near the crack tip Two types of boundary conditions must be considered to find all of the unknown parameters: near the crack tip and at the Table 1 Definitions Used in Eqs. (11)–(12). Variable

Orthotropic materialsa

(Fr)k

[( x )kcos

(Fθ)k

[ ( x )ksin

(Gr)k

[( x )k cos2

(Gθ)k

[( x )k sin2

(Grθ)k

a b

p

Isotropic materialsb 1+ p p) k

+ ( y )k sin ](

+ ( y )k cos ](

+ ( y)ksin2

{[ ( x )k + ( y)k] cos sin

= cos

1+ p p) k

+ 2(

+ ( y )k cos2

1 [cos(2 2G 1 [2(1 2G

2(

xy )kcos xy )kcos

+(

+

1 [ sin(2 2G 1 [4(1 2G

sin ](

p) k

sin ](

xy )k cos(2

)}(

p) k p) k

p

(2 p

p

+ µp sin ; (µp )1 = (µ¯ p )3, (µp )2 = (µ¯ p )4 ; Im[(µp )1, (µp)2] > 0 .

G = shear modulus = E/2(1 + ) .

545

cos(2 +

p)

p](cos

+

p)

)+

p)

cos(2 + sin(2 + p (sin p

± i·cos(2 +

p](sin

± i ·sin

p)

p )(cos p p)

p)

p p)

p

);

;

);

p

p)

);

p

± i ·cos(2 +

i·cos

p)

i·cos

i·sin(2 + ± i ·sin

];

k = 1, 3

± i·sin

p

± i·sin(2 +

p )(cos p

(2 +

± i·sin(2 +

2 )

p)

p

);

]; p

k = 2, 4

k = 1, 3 );

k = 2, 4

k = 1, 3 k = 2, 4

;

k = 1, 3

k = 2, 4 ];

k = 2, 4

k = 1, 3

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Table 2 Appropriate near tip boundary conditions for a contact crack with friction. At the bimaterial interface (θ = 0°):

ur(1) (r , 0) = ur(2) (r , 0); (1)

(r , 0) =

(2)

(r , 0);

At the crack faces (θ = ± 180°):

u (1) (r , 0) = u (2) (r , 0) (1) r (r ,

0) =

(2) r (r ,

u (1) (r , 180°) = u (2) (r , (1)

0)

(r , 180°) =

(1) r (r , (1) r (r ,

180°) =

f

180°) = (2)

(r ,

(2)

180°)

(r ,

180°)

(2) r (r ,

180°)

180°) for

(2)

(r ,

180°)

0

Fig. 2. Near tip boundary conditions and coordinate system definition. (a) Crack faces to the right of the crack tip (left tip). (b) Crack faces to the left of the crack tip (right tip).

interface with the numerical model. The boundary conditions near the crack tip provide some of the unknown constants; specifically (λp)m and Akm. The appropriate near tip boundary conditions used in the definition of the contact problem with friction are shown in Table 2 and Fig. 2. In Table 2, the coefficient of friction is denoted by f, which can assume positive or negative values depending on the direction of slip. So, assuming that the upper material (material 1) is stiffer than the lower material (material 2), positive values of f correspond to the upper material slipping with respect to the lower material in the positive x-direction, while negative values of f denote slip in the opposite direction. Stated alternatively, the sign of f is consistent with the sign of slip, so that:

sgn[f ] = sgn[ur(1) (r , 180°)

ur(2) (r ,

(15)

180°)]

In cases where the crack faces are located to the right of the crack tip (e.g., the left tip of an internal crack), the local coordinate system (r, θ) must be rotated by 180° with respect to the global x-axis and the material numbers get reversed, see Fig. 2(a). Using the new material number definitions shown in Fig. 2(a), the sign of f remains consistent with Eq. (15). Substitution of Eqs. (11)–(12) into the appropriate equations from Table 2 leads to a homogeneous system of linear equations:

r

p

(16)

[Kp ( p)]{A} = {0}

where [Kp(λp)] is an 8 × 8 matrix with entries that depend on λp, and {A} is an 8 × 1 matrix defined as follows: (17)

{A} = [A1(1) A2(1) A3(1) A4(1) A1(2) A2(2) A3(2) A4(2) ]T For non-trivial solutions to Eq. (16), the determinant of its matrix of coefficients must be equated to zero:

(18)

|Kp ( p)| = 0

The solution of Eq. (18) leads to an infinite number of roots λp. These roots were found by truncating the system and using a nonlinear equation solver based on Muller’s method [19]. Then, with the roots known, the eigenvector {A} in Eq. (16) was found by using singular value decomposition [20,21]. This process, however, only defines the eigenvectors up to arbitrary multiplicative constants, denoted by am in Eqs. (11)–(12). These unknown constants were determined from the boundary conditions at the interface with the remote region, as discussed next. 4.2. Boundary conditions at the interface with remote region The determination of the unknown constants am requires enforcement of displacement continuity at the interface between the elasticity and numerical solutions (see Fig. 1). In the numerical solution, a nonlinear finite element model with frictional contact elements was used to avoid interpenetration of the crack faces and determine the extent of the contact region. Furthermore, the best results can be achieved by using more boundary stations than unknown constants [22–24]. The resulting overdetermined system of 546

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 3. Percent error in satisfying the remote boundary conditions as a function of material properties for various friction coefficients and σ/ τ = −1.5. (a) Isotropic materials, left tip. (b) Isotropic materials, right tip. (c) Orthotropic materials, left tip. (d) Orthotropic materials, right tip.

equations was satisfied in the present work in a least squares sense by using singular value decomposition [20,21]. A total of M = 24 terms were used in the eigenfunction expansion, which typically resulted in an error of roughly 2% in satisfying the remote boundary conditions. The error was calculated as the largest element of the residual divided by the maximum boundary displacement, where the residual is the mismatch between the truncated eigenfunction series and the remote traction boundary conditions, that is

% error =

max

× 100

(19)

max

where Δmax is the largest element of the residual and δmax is the maximum boundary displacement. This error gives an indication of the effectiveness of the least squares approximation in satisfying the remote boundary conditions. As an illustration, Fig. 3 shows the percent error for various materials and coefficients of friction for the plate configuration shown in Fig. 4 and a load ratio σ/τ = −1.5. Note that the Dundurs parameter, β, and layup angle, θ, used in Fig. 3 are defined in the next section. In the majority of cases, this error was in the neighborhood of 2%, with a few exceptions having an error as high as 4%. These exceptions occurred in the case of a frictionless crack with orthotropic materials. It should be noted that the resulting error in the stress intensity factor is generally lower than this error. Overall, errors this low in satisfying the boundary conditions at the local/remote interface result in excellent values of the stress intensity factors as determined originally by Penado [17] for a frictionless crack and confirmed in the course of the present investigation for a crack with friction. 5. Analysis of remote region 5.1. Cracked bimaterial plate under shear and axial loading The problem analyzed consists of an infinite plate with a central crack at the interface of two isotropic or orthotropic materials (or 547

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 4. Plate configuration and loading used.

a combination of both, such as a bonded joint with composite adherends and isotropic adhesive) and subjected remotely to a combination of axial and shear loading. This configuration was chosen in order to be able to consider the effect of friction and compare directly with results in the literature for the limiting case of isotropic materials. The plate configuration used is shown in Fig. 4. In order to ensure that the dimensions approximate those of an infinite plate, overall dimensions of 254 × 254 mm with a crack length of 2a = 25.4 mm were chosen, which produced a plate-width to crack-length ratio of 10. In order to get a rough idea of the effect of finite dimensions, for tensile loading the ratio of finite-width to infinite-width stress intensity factor is approximately KI / KI sec( a/ W ) [25] where a = crack half-length and W = plate width. In the present case with a/W = 0.05, the ratio KI/ KI∞ = 1.006, and the difference due to the finite size effect is on the order of 0.6%. 5.2. Finite element model The displacement field of the plate providing the boundary conditions for the local region was found by using the finite element model shown in Fig. 5. The analysis was conducted using the commercially available finite element code ANSYS [26], which includes frictional contact element capabilities. The model consisted of 21,600 8-noded quadrilateral plane strain elements (PLANE 183) and 65,560 nodes. In addition, 160 2D 3-node surface-to-surface contact elements (TARGE169 and CONTA172) were used at the crack faces to simulate the frictional contact problem and avoid crack face interpenetration. A static (steady-state) analysis type incorporating nonlinear geometric effects was used. The element size in the crack region in the x- and y- directions was, respectively, 0.3175 and 0.635 mm. The local region size used for both crack tips was 7.62 × 7.62 mm (Fig. 5). Axial loading was applied by

τ

CRACK

TYPICAL LOCAL REGION SIZE USED FOR LEFT AND RIGHT CRACK TIPS IS 7.62 mm x 7.62 mm

y x

Fig. 5. Mesh, loading and boundary conditions used in the finite element analysis. 548

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Table 3 Effective orthotropic material properties used in the analysis. Material Number

Rotation about y-axis

Ex (GPa)

Ey (GPa)

Gxy (GPa)

νxy

1 (fixed) 2 (rotated) ↓

0-deg 0-deg 15-deg 30-deg 45-deg 60-deg 75-deg 90-deg

141.00 141.00 52.29 20.87 12.89 10.40 9.74 9.65

9.65 9.65 9.65 9.65 9.65 9.65 9.65 9.65

4.90 4.90 4.80 4.54 4.22 3.95 3.77 3.71

0.300 0.300 0.213 0.195 0.214 0.248 0.284 0.300

constraining the bottom edge of the plate from displacements in the y-direction, while stresses were applied on the top edge. Shear loading was achieved by tangential and self-equilibrating stresses on opposite edges, as shown in Fig. 5. For isotropic materials, the properties for materials 1 and 2 were varied to produce a range of the Dundurs parameter, β. The Dundurs parameters, α and β, used in this investigation to describe the isotropic material property mismatch, are defined as [27]:

=

G1 ( G1 (

2 2

+ 1) G2 ( + 1) + G2 (

+ 1) , 1 + 1) 1

=

G1 ( G1 (

2 2

1) G2 ( + 1) + G2 (

1) + 1) 1 1

(20)

where G = shear modulus = E/[2(1 + ν)], E = Young’s modulus, ν = Poisson’s ratio, κ=3 − 4ν for plane strain, and the subscripts 1 and 2 refer to the material number. Note that β > 0 if material 1 is stiffer than material 2, and the sign of β is reversed when materials 1 and 2 are exchanged. In all cases except β = 0.5, the properties for material 1 were kept constant at E1 = 68.95 GPa and ν1=0.3, while the properties of material 2 were varied. For β=0.5, which implies that material 1 is rigid with respect to material 2, the properties used were E1 = 6.89 × 105 GPa, E2 = 3.45 × 10−2 GPa, and ν1=ν2=0, giving a value of β very close to 0.5. For orthotropic materials, the material variation used corresponds to the angle of rotation of a unidirectional composite ply about the y-axis. It should be noted, however, that for angles other than 0° and 90° coupling of in-plane and out-of-plane effects occurs. However, the out-of-plane coupling effects caused by off-axis plies may be small in many cases, and effective ply moduli can be used as an approximation, as shown by Penado [28]. The resulting effective orthotropic properties used are shown in Table 3. 5.3. Contact and normal displacement of the crack faces The implementation of the method of solution requires knowledge of the remote displacement field, as well as the appropriate frictional contact conditions at the crack faces. The extent of the contact region was found directly from the finite element results and is defined in Fig. 6 by the parameter b, which can be positive or negative. Note that for a completely open crack b/a = 1, while for a completely closed crack b/a = −1. Generally, unless the crack is completely open or closed, for a given value of the ratio of axial to shear stress, σ/τ, a large and a small contact region can be observed. For a positive applied shear, on the right end of the crack a large contact region of length a − b occurs, while at the other end a very small contact region on the order of 10−8 of the crack length is present [29]. Following [30], the dimensionless normal gap is defined as

V =

2G2 (1 + ) (1

2 )(

2

+ 1) a

2

+

2

|v (1) (x , 0+)

v (2) (x , 0 )|

(21)

where v(1)(x, 0+) and v(2)(x, 0−) are, respectively, the vertical displacements of the upper and lower crack faces, and the coordinate origin is at the center of the crack as shown in Fig. 6. Typical results for the normal gap of isotropic materials are given in Fig. 7 for σ/

Fig. 6. Nomenclature used to define the partially closed crack. 549

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 7. Dimensionless normal gap comparison with [15] for σ/τ = 0 (pure shear) and a variety of Dundurs parameters, β, and coefficients of friction, f.

τ = 0 (pure shear) and several values of β and f. The present results are also compared with the results given by Dorogoy and BanksSills [15], which they obtained by using a solution based on finite differences and validated by comparison with analytical and finite element solutions. Excellent agreement between the present results and [15] is observed, providing confirmation of the validity of the finite element model used. As anticipated, the size of the gap decreases with increasing f or decreasing β. 6. Results and discussion 6.1. Strong and weak singularities at crack tips The customary inverse square root singularity, r 1/2 , that is observed in homogeneous materials (either isotropic or anisotropic) is also present in bimaterial cracks when the type of contact is frictionless. The r 1/2 singularity can be observed also in homogeneous materials even when the crack contact includes friction. However, for a crack between two dissimilar materials in frictional contact, the singularity becomes r , where the order of the singularity, λ, can be smaller (weak singularity) or larger (strong singularity) than 1/2. This has a number of implications with respect to the solution of a fracture mechanics problem. First, as pointed out by Qian and Sun [11], the classical definition of strain-energy release rate becomes zero for the weak singularity and unbounded for the strong singularity, and therefore cannot be used directly as a fracture criterion. Second, the selection of the suitable set of roots (λp)m that lead to the correct definition of the displacements, stresses, and stress intensity factors in Eqs. (11)–(14) must be based on the proper evaluation of the relative tangential slip of the crack faces. In order to illustrate the occurrence of weak and strong singularities, values for the singular eigenvalues at the right and left tips obtained from Eq. (18) for a variety of material combinations and coefficients of friction are shown in Fig. 8. It can be seen that, while the isotropic eigenvalues are linear with respect to the Dundurs parameter, β, the corresponding orthotropic values are nonlinear with respect to the layup angle, θ. In all cases the eigenvalues are real and the right and left singular eigenvalues are seen to be related by

Fig. 8. Singular eigenvalues (λp)1 at the right and left tips for a variety of material combinations and coefficients of friction. (a) Isotropic materials. (b) Orthotropic materials. 550

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Table 4 Eigenvalues from Eq. (18) up to M = 24 for two material and friction coefficient combinations. Root

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Isotropic materials, β=0.3

Orthotropic materials, θ=90°

f=0

f=0

−0.500000 0.000000 0.500000 1.000000 1.500000 2.000000 2.500000 3.000000 3.500000 4.000000 4.500000 5.000000 5.500000 6.000000 6.500000 7.000000 7.500000 8.000000 8.500000 9.000000 9.500000 10.000000 10.500000 11.000000

f = 0.6 Strong

Weak

−0.556689 0.000000 0.443311 1.000000 1.443311 2.000000 2.443311 3.000000 3.443311 4.000000 4.443311 5.000000 5.443311 6.000000 6.443311 7.000000 7.443311 8.000000 8.443311 9.000000 9.443311 10.000000 10.443311 11.000000

−0.443311 0.000000 0.556689 1.000000 1.556689 2.000000 2.556689 3.000000 3.556689 4.000000 4.556689 5.000000 5.556689 6.000000 6.556689 7.000000 7.556689 8.000000 8.556689 9.000000 9.556689 10.000000 10.556689 11.000000

−0.500000 0.000000 0.500000 1.000000 1.500000 2.000000 2.500000 3.000000 3.500000 4.000000 4.500000 5.000000 5.500000 6.000000 6.500000 7.000000 7.500000 8.000000 8.500000 9.000000 9.500000 10.000000 10.500000 11.000000

f = 0.6 Strong

Weak

−0.525567 0.000000 0.474433 1.000000 1.474433 2.000000 2.474433 3.000000 3.474433 4.000000 4.474433 5.000000 5.474433 6.000000 6.474433 7.000000 7.474433 8.000000 8.474433 9.000000 9.474433 10.000000 10.474433 11.000000

−0.474433 0.000000 0.525567 1.000000 1.525567 2.000000 2.525567 3.000000 3.525567 4.000000 4.525567 5.000000 5.525567 6.000000 6.525567 7.000000 7.525567 8.000000 8.525567 9.000000 9.525567 10.000000 10.525567 11.000000

(22)

( p)1, left = 1 + ( p)1, right

Using the material definition in Fig. 2(b) in conjunction with Eq. (18) produced weak singularities, |(λp)1| < 1/2, while definition as in Fig. 2(a) provided strong singularities |(λp)1| > 1/2. The isotropic singular eigenvalues also satisfy Comninou’s determinantal equation given by [10]:

cot(

(23)

) = ±f

where the “+” sign leads to weak, and the “−” sign to strong, singularities when β > 0. In order to further elucidate the solutions of Eq. (18) for the eigenvalues, solutions up to M = 24 roots including nonsingular values are shown in Table 4 for some sample material combinations. 6.2. Bimaterial central crack under combined shear and compression For a central crack subjected to a combination of axial loading and shear, several different contact regions can be observed depending on the loading and material properties. For axial loading consisting of tensile or small compressive values, a large contact region at the right tip with a significant open region at the left tip develops, similar to the situation illustrated in Fig. 6 assuming that material 1 is stiffer than material 2. Upon close examination, the open region is not totally open but has, in fact, a very small contact region at the left tip, of the order of 10−8 of the crack length [29]. Due to the extremely small size of the contact zone at the left tip, the results for the stress intensity factor can be approximated by considering the results assuming stress free crack faces [16,17,31]. However, as the compressive load on the plate is increased, the left tip opening becomes increasingly smaller but can never completely close, except for β = 0, for finite values of compressive loading [11,32]. This behavior is illustrated in Fig. 9(a), which shows the nonlinear variation of the open contact zone as a function of the ratio σ/τ for negative values of σ (compression) for β = 0.49 and f = 0.5. The present results are also compared with those of Qian and Sun [11] as a check and good agreement is found. It should be noted that Qian and Sun [11] used 4-noded quadrilateral finite elements, while higher order 8-noded elements were used in the present investigation, which may account for the small discrepancies. Since the size of the open zone can be very small for low values of σ/τ, the mesh in Fig. 5 was refined near the left crack tip in order to achieve an element size of 0.0353 × 0.0706 mm, or 1/9 the size of the original mesh, allowing for very small gaps to be detected. The refined mesh is shown in Fig. 10. Additionally, the gap along crack faces is shown in Fig. 9(b). It can be seen that both the gap opening and open contact zone become progressively smaller as the compressive load increases, but never go away completely. A similar trend has been reported by Comninou and Schmueser [32] and Comninou and Dundurs [33] for the size of the open contact zone when f = 0; however, they reported an increasing gap magnitude with increasing compression. Results for f = 0 reported by Van der Zande and Grootenboer [34] show a decreasing gap 551

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 9. Variation of open contact zone at left tip with load ratio σ/τ for β = 0.49 and f = 0.5. (a) Size of open contact zone, b/a. (b) Dimensionless gap along crack faces, V .

CRACK 6.35 mm

6.35 mm Fig. 10. Refined finite element mesh near left crack tip. A typical local region is shown by the dotted contour.

Contact region, f=0

Contact region, f=0.6

Fig. 11. Normal stresses along the crack faces for isotropic materials for σ/τ = −1.5, β = 0.3, and f = 0 and 0.6. The inserts show a close up of the crack tips for the regions enclosed by the rectangles.

magnitude with increasing compression, consistent with the present trends (Fig. 9b), and a disagreement with [32] was pointed out in their paper [34]. The distribution of normal stresses along the crack faces is illustrated in Fig. 11 for the special case of σ/τ = −1.5, β = 0.3, and two values of the coefficient of friction including the frictionless case, f = 0. It can be seen that increasing friction leads to a larger and flatter contact region, to the extent that the normal stress distribution appears nearly constant for most of the crack length when f = 0.6. Regarding the behavior at the crack tips, some comments are now pertinent. At the left tip, as discussed previously, there is 552

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 12. Comparison of shear stresses ahead of both crack tips for isotropic materials and σ/τ = −1.5, β = 0.3, f = 0.6. (a) Left crack tip. (b) Right crack tip.

an extremely small contact zone followed by an open region prior to reaching a large contact region. The open and large contact regions have zero and negative contact stresses, respectively, and can be seen clearly in Fig. 11 for the finite element results. However, these approximate results do not have the resolution to detect the small contact region, which is of the order of 10−8 of the crack length [29]. At the same time, the analytical results from the local eigenfunction solution are also approximate at the left tip because they are based on the assumption of a completely closed crack with no open region. Although this assumption provides a very good approximation when the open region is small, it can be seen that it is not strictly correct. The open region should have zero normal stresses but the analytical results show the normal stress becoming positive. This behavior is the result of the condition of continuity of normal stresses at the crack faces in Table 2. On the other hand, at the right tip, the crack is closed with no open region and the normal stresses become singular. It can be seen that the local analytical results with a closed crack provide an excellent representation of the singular stresses, but the finite element results do not. The latter display oscillatory behavior near the right crack tip for both f = 0 and f = 0.6. Such irregular behavior has also been observed by others [15] and is believed to be caused by the inability of the numerical results to model accurately the singular field. 6.3. Shear stress distribution ahead of crack tips In order to assess the validity of the proper set of eigenvalues that must be used to calculate the stresses when weak and strong singularities are present, the shear stresses from the eigenfunction expansion ahead of both crack tips were evaluated and compared with finite element results. The results are shown in Figs. 12 and 13, which correspond to isotropic and orthotropic materials, respectively. The finite element stress results were calculated using a refined mesh at both tips similar to the one shown in Fig. 10. The analytical results included both singular and non-singular terms from the third of Eqs. (12) with M = 24 roots, and were

Fig. 13. Comparison of shear stresses ahead of both crack tips for orthotropic materials and σ/τ = −1.5, θ = 90°, f = 0.6. (a) Left crack tip. (b) Right crack tip. 553

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Fig. 14. Normalized mode II stress intensity factor, KII , at both crack tips for isotropic materials, σ/τ = −1.5, and all β. (a) Left crack tip. (b) Right crack tip.

calculated using the regular mesh shown in Fig. 5 to determine the displacement boundary conditions at the interface with the remote region. These results confirm the anticipated behavior: on the left tip the analytical results from the eigenfunction expansion and strong singularity (see Table 4) closely agree with the finite element results, while the weak singularity terms do not properly describe the near-tip stress distribution. Conversely, at the right tip it is the weak singularity terms that properly describe the near-tip stress distribution. A similar result for the strong singularity at the left tip of an isotropic bimaterial crack was reported by Qian and Sun [11] by using an approximate one-term analytical solution based on the singular term only and a stress intensity factor calculated from the finite crack extension strain energy release rate. 6.4. Isotropic stress intensity factors and comparisons with published results Results for the stress intensity factor for isotropic materials are now presented. At a closed crack tip, only the mode II stress intensity factor KII is present. All displacements at the local/remote interface were obtained using the baseline finite element mesh in Fig. 5, with the exception of the convergence evaluation discussed later in relation to Fig. 14. The latter necessitated the use of the refined mesh shown in Fig. 10. Also, in order to gain confidence in the validity of the results, comparisons were made with published values whenever possible. For ease of comparison, we now define the normalized stress intensity factor as:

KII =

KII a

( p)1·

2

+

(24)

2

Table 5 presents the comparison of the present results with the analytical solution of Nemat-Nasser and Horii [35] for a homogeneous, isotropic plate under combined shear and compression:

KII =

a(

(25)

f )

It can be seen that there is excellent agreement for a range of coefficients of friction for a representative loading ratio σ/τ = −1, with an error below 2% in all cases. For dissimilar materials, Table 6 displays the comparison with the analytical results of Gautesen Table 5 Comparison between analytical and present results at the right tip for σ/τ = −1 and a homogeneous, isotropic plate (β = 0). f

0.0 0.1 0.2 0.3 0.4 0.5 0.6

Normalized stress intensity factor, KII Nemat-Nasser and Horii [35]

Present

% error

1.253 1.128 1.003 0.877 0.752 0.627 0.501

1.236 1.112 0.988 0.864 0.740 0.616 0.492

−1.36 −1.42 −1.50 −1.48 −1.60 −1.75 −1.80

554

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

Table 6 Comparison between published and present results at the right tip for σ/τ= 0 and various β and f. f

Normalized stress intensity factor, KII β = 0.25

0 0.5 1.0

β = 0.5

Literature

Present

% error

Literature

Present

% error

1.79 [29] 1.74 [15] 1.68 [15]

1.78 1.73 1.69

−0.56 −0.57 0.60

1.83 [29] 1.72 [15] 1.62 [15]

1.84 1.74 1.65

0.55 1.16 1.85

and Dundurs [29] for the special case of f = 0 and σ/τ = 0, as well as results from Dorogoy and Banks-Sills [15] from a finite difference solution for various β and f. Again, very close agreement is observed with the error under 2% in all cases and under 1% in most cases. The results for the normalized stress intensity factor, KII , at the left and right crack tips for the full range of Dundurs parameter 0≤β≤0.5 and a loading ratio of σ/τ = −1.5 are shown in Fig. 14. Note that, at the extremes, β = 0 corresponds to equal materials (homogeneous crack), while β = 0.5 if one of the materials is rigid. Available results in the literature for β = 0 and any f [35], as well as all β but f = 0 [29] are superimposed in Fig. 14 for comparison. It should be observed that Gautesen and Dundurs [29] results are valid for σ/τ = −0.4; however, they are compared with the present results under the assumption that, for frictionless contact, the value of compression has a negligible effect on KII when b/a → − 1. Additionally, it was assumed that Nemat-Nasser and Horii results [35] have the same magnitude at the right and left crack tips. Close agreement for all limiting cases examined is observed. Although for the frictionless crack, f = 0, we expect the values at the right and left tips to be equal in magnitude, this is no longer the case when β > 0 and friction is present, as a result of the unequal singularity strengths at the two tips (see Fig. 8(a)). In order to check the convergence of the results, the refined mesh in Fig. 10 (with a similar configuration for the right tip), was used in the analysis to calculate another set of values of KII for f = 0.6. For the refined mesh, a smaller local region of 3.67 × 3.67 mm was used around the crack tips. The results are plotted in Fig. 14 and show excellent agreement with those from the original mesh, indicating that mesh convergence has occurred. Additionally, the number of terms in the truncated infinite series, M, in Eqs. (11) was varied from 12 to 72 (the baseline is M = 24) and the results were virtually unchanged. This indicates that the mesh refinement used in the baseline mesh in Fig. 5 provides satisfactory convergence. As an additional check, the following relationship given by Comninou [9,36] for the frictionless crack (f = 0) was checked and found to be exactly satisfied in all cases at both the left and right crack tips:

K I(

)

(26)

= ± KII

where K I( ) is a stress intensity factor associated with compression, which is always negative, and the “+” or “−” sign applies when the bonded interface is to the left or the right of the crack tip, respectively. The relationship given by Eq. (26) implies that KII has opposite signs at the two crack tips. Both K I( ) and KII were calculated from Eqs. (14).

Fig. 15. Normalized mode II stress intensity factor, KII , at both crack tips for orthotropic materials and σ/τ = −1.5 as a function of layup angle, θ. (a) Left crack tip. (b) Right crack tip. 555

Engineering Fracture Mechanics 204 (2018) 542–556

F.E. Penado

6.5. Orthotropic materials Next, the normalized stress intensity factor was calculated for a cracked bimaterial plate with orthotropic materials and a loading ratio σ/τ = −1.5. The results for layup angles in the range of 0–90° as defined in Table 3 are presented in Fig. 15. No results are found in the literature for the orthotropic case to allow for a comparison of results as was possible in the isotropic case. As expected, increasing the coefficient of friction, f, results in a decrease in the stress intensity factor. The values at the right and left tips are practically equal at the right and left tips for the special case of f = 0, when the conventional r−1/2 singularity prevails. For f > 0, there are small differences between the right and left tip as we go from 0° to 90°, with the KII values being virtually identical at both tips for θ = 0° regardless of f, tend to be lower at the right tip (weak singularity) as the layup angle increases, and become increasingly lower for larger f. The largest decrease, 24%, occurs for the largest coefficient of friction and layup angle (f = 0.6 and θ = 90°). Note that a coefficient of friction of 0.6 is typical for delamination surfaces of graphite fiber/epoxy matrix composites [37]. 7. Conclusions An approach for computing directly the stress intensity factors of a crack in a bimaterial interface under combined loading has been demonstrated. The approach extends existing capabilities to enable consideration of crack surfaces in frictional contact, allowing the solution of problems with complex geometries, materials, loading, and boundary conditions. The method is portable, efficient, and requires no special finite element codes or elements. As a specific example, a bimaterial plate with a central crack under combined shear and compression loading was considered. In this case, the existence of friction alters the r−1/2 singularity to become either weak or strong and causes the conventionally defined strain energy release to become either zero or infinity. Direct calculation of mode I and II stress intensity factors as well as near-tip stresses with singular and non-singular terms allows the present approach to provide a powerful tool for analyzing fracture problems. Results have been presented for a variety of isotropic and orthotropic material combinations and coefficients of friction. In order to validate the method, the present results were compared with a variety of closed form and numerical results found in the literature, producing excellent agreement. Overall, the method provides a practical way to analyze a variety of engineering structures with frictional interfacial cracks and including isotropic or orthotropic materials. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

Rogel L, Banks-Sills L. A through interface crack between a transversely isotropic pair of materials (+30/-60,-30/+60). Eng Fract Mech 2010;77:3261–91. Banks-Sills L. Review on interface fracture and delamination of composites. Strain 2014;50:98–110. Penado FE. Determination of stress intensity factors in cracked bonded composite joints. J Thermoplast Compos Mater 2000;13:174–89. Haraga S, Okazaki N, Oda K, Tsutsumi N. Evaluation of adhesive scarf joint strength by using singular stress field of small interface crack. IOP Conf Ser: Mater Sci Eng 2018;372:012015. Dorogoy A, Banks-Sills L. Effect of crack face contact and friction on Brazilian disk specimens—A finite difference solution. Eng Fract Mech 2005;72:2758–73. Kessler H, Schuller T, Beckert W, Lauke B. A fracture-mechanics model of the microbond test with interface friction. Compos Sci Technol 1999;59:2231–42. Dong W, Yang D, Zhang B, Wu Z. Rock-concrete interfacial crack propagation under mixed mode I-II fracture. J Eng Mech 2018;144:04018039. Govorukha V, Kamlah M, Loboda V, Lapusta Y. Interface cracks in piezoelectric materials. Smart Mater Struct 2016;25:023001. Comninou M. The interface crack. J Appl Mech 1977;44:631–6. Comninou M. Interface crack with friction in the contact zone. J Appl Mech 1977;44:780–1. Qian W, Sun CT. A frictional interfacial crack under combined shear and compression. Compos Sci Technol 1998;58:1753–61. Sun CT, Qian W. A treatment of interfacial cracks in the presence of friction. Int J Fract 1998;94:371–82. Comninou M, Dundurs J. Effect of friction on the interface crack loaded in shear. J Elast 1980;10:203–12. Comninou M. The interface crack in a shear field. J Appl Mech 1978;45:287–90. Dorogoy A, Banks-Sills L. Shear loaded interface crack under the influence of friction: a finite difference solution. Int J Numer Methods Eng 2004;59:1749–80. Sun CT, Qian W. The use of finite extension strain energy release rates in fracture of interfacial cracks. Int J Solids Struct 1997;34:2595–609. Penado FE. Analysis of the partially closed interface crack in orthotropic materials. J Thermoplast Compos Mater 2018;31:820–36. Ting TCT, Chou SC. Edge singularities in anisotropic composites. Int J Solids Struct 1981;17:1057–68. Muller DE. A method for solving algebraic equations using an automatic computer. Math Comput 1956;10:208–15. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in Fortran. Cambridge: Cambridge University Press; 1992. Dongarra JJ, Bunch JR, Moler CB, Stewart GW. LINPACK Users’ Guide. Philadelphia: Society for Industrial and Applied Mathematics; 1979. Penado FE. Analysis of singular regions in bonded joints. Int J Fract 2000;105:1–25. Gross B, Mendelson A. Plane elastostatic analysis of V-notched plates. Int J Fract Mech 1972;8:267–76. Carpenter WC. A collocation procedure for determining fracture mechanics parameters at a corner. Int J Fract 1984;24:255–66. Feddersen CE. Discussion. ASTM STP 1967;410:77–9. ANSYS, Inc. ANSYS Mechanical APDL Release 19.0. Canonsburg, PA; 2017. Dundurs J. Discussion: edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J Appl Mech 1969;36:650–2. Penado FE. Out-of-plane coupling effects on stress intensity factors at interface cracks in composite lap joints. J Thermoplast Compos Mater 2001;14:178–200. Gautesen AK, Dundurs J. The interface crack under combined loading. J Appl Mech 1988;55:580–6. Gautesen AK. The interface crack under combined loading: an eigenvalue problem for the gap. Int J Fract 1993;60:349–61. Penado FE. A computational procedure for the determination of fracture parameters at interfacial cracks in composite laminates. J Compos Mater 2006;40:701–15. Comninou M, Schmueser D. The interface crack in a combined tension-compression and shear field. J Appl Mech 1979;46:345–8. Comninou M, Dundurs J. On the behaviour of interface cracks. Res Mechanica 1980;1:249–64. Van der Zande HD, Grootenboer HJ. A finite element approach to interface cracks. J Appl Mech 1986;53:573–8. Nemat-Nasser S, Horii H. Compression-induced nonplanar crack extension with application to splitting, exfoliation and rockburst. J Geophys Res 1982;87:6805–21. Comninou M. An overview of interface cracks. Eng Fract Mech 1990;37:197–208. Schön J. Coefficient of friction of composite delamination surfaces. Wear 2000;237:77–89.

556