A technique for studying interacting cracks of complex geometry in 2D

Engineering Fracture Mechanics 73 (2006) 1086–1114 www.elsevier.com/locate/engfracmech

A technique for studying interacting cracks of complex geometry in 2D S.C. TerMaath a b

a,*

, S.L. Phoenix b, C.-Y. Hui

b

Exponent Failure Analysis Associates, 149 Commonwealth Drive, Menlo Park, CA 94025, USA Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14853, USA

Received 5 December 2003; received in revised form 17 September 2004; accepted 19 September 2004 Available online 3 February 2006

Abstract A method for studying brittle fracture in an infinite plate containing interacting cracks of complex shape under general loading conditions is developed and studied for accuracy and potential applications. This technique is based on superposition and dislocation theory and can be used to determine the full stress and displacement fields in a cracked body. In addition, stress intensity factors at both crack tips and wedges, created by crack kinking and branching, are calculated so that crack growth and initiation can be analyzed at these locations of possible crack propagation. Such information can then be used to study damage accumulation in structures containing a large number of interacting cracks.
2004 Published by Elsevier Ltd. Keywords: Superposition; Dislocation distribution; Interacting cracks; Stress intensity factor; Brittle fracture

1. Introduction Cracks may change direction or branch as they grow due to many factors including fatigue loading conditions, preferred crack paths such as grain boundaries, material imperfections, or environmental effects (i.e., corrosion). Damage zones containing many of these randomly shaped cracks of varying size are common in aging structures. During early formation, these crack arrays, as illustrated in Fig. 1, typically consist of small cracks that may not adversely affect structural performance. However, subsequent loading could cause these cracks to propagate and coalesce into a larger failure-inducing crack. Components plagued by crack arrays include a broad range of materials and applications, and staggering replacement costs often *

Corresponding author. Tel.: +1 817 737 3131. E-mail address: [email protected] (S.C. TerMaath).

0013-7944/$ - see front matter
2004 Published by Elsevier Ltd. doi:10.1016/j.engfracmech.2004.09.009

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1087

Fig. 1. Random array of cracks of complex shape.

dictate life extension over part replacement. To ensure public safety and preserve structural integrity when exercising a life extension option, it is imperative to accurately predict the behavior of existing crack arrays to assess structural reliability, safety, and performance. Evaluating such damage zones poses a challenging problem. To predict damage propagation and structural failure in cracked regions, the stress state in the material and stress intensity factors must be determined at all locations and in all directions of possible crack growth. The stress state and stress intensity factors must be calculated at crack tips, wedges formed by crack kinking (or turning), and small-scale material flaws (which may be justifiably ignored in initial modeling, but may grow and become important if located in an area with a high stress state). Limiting the scope of this problem to two-dimensional components modeled with isotropic, linear elastic material properties still leaves complex mathematical barriers to overcome. To study interacting crack behavior under the aforementioned conditions, an analytical weighted superposition technique based on the dislocation distribution method [1] is further developed. This method is applicable to arrays of cracks of virtually any configuration with an arbitrarily selected number of straight segment lengths. Cracks need not be symmetric, periodically placed, nor limited to a certain number of crack tips or growth directions. Furthermore, loading is not restricted by type and can be any combination of shear and normal loading modes. With this technique, key equations are solved analytically, and, only minimal numerical approximations are necessary when evaluating opening displacements of cracks and stress singularities at wedge locations. Since both approximations have a negligible effect on final results, continued development of this method is a solid step forward towards accurately predicting the evolution of crack arrays. Applications for three different types of crack shapes, V-shaped, multiply kinked, and branched, are presented. Convergence and parameter studies provide insight into efficient implementation of the method, while verification cases are included to assess accuracy. In addition, a propagation example demonstrates the use of this method to study growth of damage zones and the link-up of small cracks into a possible failure-inducing crack. 1.1. Literature review Brittle rock fracture, cracking in composites, and fatigue cracking in aging aircraft are a few of the many fracture examples that may consist of crack arrays or distributed fractures. Attention in the literature has focused on many diverse solution methods for widespread applications involving this type of damage. For

1088

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

example, dislocation-based analytical and hybrid-analytical methods to analyze arrays of multiple cracks are abundant in the literature (e.g. [2–11] as a representative sampling). To solve the singular integral equations involved in these methods, many researchers follow the numerical solution schemes of Erdogan et al. [12], or others [13–16]. An overview of applications and developments of the dislocation approach is provided by Hills et al. [17]. Unfortunately, few exact solutions are available for interacting cracks and cracks of complex shape; therefore, accuracy comparisons for new approaches are limited to specific crack configurations [18–22]. To evaluate cracks of complex shape, singularities occurring at kinks and branches must be considered. Williams [23,24] performed some of the first work on the classical elasticity problem of singularities at material wedges. Timoshenko and Goodier [25] and Barber [26] both included review chapters on wedge theory. Meanwhile, multiple authors have studied kinked and branched cracks [27–36]. A more comprehensive literature review can be found in TerMaath [37]. Since most researchers evaluate stress singularities and crack propagation only at crack tips, the possibility of crack development and growth is neglected elsewhere in the material. When included at wedge locations, these singularities are usually approximated, and most methods utilize numerical techniques for integration, leading to additional approximations and inaccuracies. Moreover, the application of many methods is restricted to straight cracks, and often only a small number of cracks can be evaluated in a reasonable amount of time due to computational requirements. To overcome these limitations, the following technique was developed for determining stress fields and stress intensity factors in a cracked, two-dimensional, isotropic, linearly elastic material under conditions of brittle fracture.

2. Mathematical formulation and concepts The mathematical formulation is based on the dislocation distribution method, coupled to the principle of superposition applied at both the global and local levels. Additional reviews of dislocation theory and the dislocation distribution method are available [17,26,38,39]. 2.1. Global superposition At the global level, the cracked plate solution is the sum of two separate solutions, trivial and auxiliary (Fig. 2). The trivial problem is the given plate under the specified far field loading but without the cracks. The auxiliary problem is the given cracked plate without the far field loading, but instead, prescribed tractions applied to the crack faces. To fulfill the original boundary condition of traction-free crack faces, these tractions are equal and opposite to those induced by the stresses in the uncracked material at the location of the crack faces. Obtaining the solution to the auxiliary problem comprises the bulk of the analytical effort. Satisfying the prescribed tractions necessitates the calculation of the corresponding opening displacement profile (shape of the deformed crack) for every crack in the array. 2.2. Local superposition Solving the auxiliary problem for the stress field and stress intensity factors [40] requires the development and superposition of detailed crack geometric features. To determine the opening displacements, cracks are subdivided into a series of straight segments. Each crack segment is treated as though it alone is open in an infinite plate but with an opening displacement form in terms of certain shapes that anticipate those needed to model the full crack system. A local coordinate system (xi, yi) is aligned along each crack segment i. For example, each arm of a V-shaped crack represents a crack segment of length ai, i = 1, 2, with a unique local coordinate system (Fig. 3). To simplify notation in equations, the arbitrary parameter t is aligned with xi. For the V-shaped crack case, each segment is closed at one end (a crack tip at t = ai). The other end (t = 0)

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1089

Fig. 2. Global superposition.

Fig. 3. Coordinate systems and notation for a V-shaped crack.

is open, since it must connect to the opposite segment and satisfy continuity. Superposition of the local solutions for all crack segments yields the full solution to the auxiliary problem.

1090

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

2.3. The distributed dislocation approach Dislocations are utilized to create the prescribed crack face tractions of the auxiliary problem. Glide dislocations will produce relative tangential sliding between opposing crack faces, while climb dislocations will induce relative normal displacements. For crack segments of finite length, opening displacements are generally not equal at both ends, and the displacement between crack faces along the crack segment is not constant. Therefore, a single climb or glide dislocation is insufficient, and an array of dislocations must be applied. For crack segment, i, these dislocation distributions are (l1i(t), l2i(t)) for the tangential and normal directions respectively. In local coordinate system (xi, yi) the stress induced at an arbitrary point (x, y) caused by crack segment i acting alone (as though all other crack segments are closed) is written as Z ai 
 i 2G 1 h
2 2 ðiÞ 2 2 sxy ðx; yÞ ¼
ðx
tÞ
y þ ðx
tÞ
y yl ðx
tÞl dt 2i 1i pð1 þ jÞ 0 w4 Z ai 
 i 2G 1 h
2 2 sðiÞ ðx
tÞ þ 3y 2 ðx
tÞl2i þ ðx
tÞ
y 2 yl1i dt ð1Þ yy ðx; yÞ ¼
4 pð1 þ jÞ 0 w Z ai 
 i 2G 1 h
2 2 2 2 sðiÞ ðx
tÞ
y þ 3ðx
tÞ
y ðx
tÞl yl dt 2i 1i xx ðx; yÞ ¼
pð1 þ jÞ 0 w4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 where w = ðx
tÞ þ y 2 . The shear modulus of the material is G, and m is Poisson’s ratio. Kosolov’s constant is denoted by j ((3
4m), for plane strain, and (3
m)/(1 + m), for plane stress). These equations can be recast with a complex variable formulation, z = x + iy (Eqs. (2)).  2G  y ReðZ 22i Þ þ ReðZ 11i Þ þ y ImðZ 12i Þ sðiÞ xy ðzÞ ¼
pð1 þ jÞ  2G  ReðZ 21i Þ
y ImðZ 22i Þ þ y ReðZ 12i Þ sðiÞ ð2Þ yy ðzÞ ¼
pð1 þ jÞ  2G  ReðZ 21i Þ þ y ImðZ 22i Þ þ 2 ImðZ 21i Þ
y ReðZ 12i Þ sðiÞ xx ðzÞ ¼
pð1 þ jÞ Z g1i and Z g2i , are Cauchy singular integrals to be evaluated in closed form in terms of the dislocation distributions (Eqs. (3)). When the point (x, y) falls on a crack segment, these integrals are solved for the Cauchy principal value. Z ai lgi ðtÞ dt Z g1i ¼ z
t Z0 ai ð3Þ l d g gi ðtÞ dt Z ¼
Z g2i ¼ dz 1i ðz
tÞ2 0 The total stress at a given point in the auxiliary problem, r, is the sum of the contributions, s (Eqs. (2)), from all crack segments after they are converted to the global coordinate system (identified by capital letters). With c as the number of crack segments, c c c X X X ðiÞ ðiÞ ðiÞ rXY ¼ sXY rYY ¼ sYY rXX ¼ sXX ð4Þ i¼1

i¼1

i¼1

2.4. Enforcing the traction-free condition To obtain traction-free crack faces in the full problem, series of equations in both the tangential and normal directions are enforced simultaneously at a given set of points along each crack segment. At point (X, Y) in the global coordinate system, these equations take the form

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114 1 r1 XY nY þ rXX nX ¼
nY rXY
nX rXX 1 r1 YY nY þ rXY nX ¼
nY rYY
nX rXY

1091

ð5Þ

The left sides of Eqs. (5) are the prescribed tractions induced at the crack faces by the loading conditions, while the right sides represent the tractions caused by the opening displacements (dislocation distributions) of the crack segments (Eqs. (4)). Also, nX and nY are the normals to the bottom crack face. The r1 are the far field stresses applied to the plate in the directions denoted by their subscripts. Dislocation distributions, l1i and l2i, must be determined that satisfy these equations for every crack segment. 2.5. Solving for dislocation distributions Enforcing the prescribed tractions (Eqs. (5)) at a sufficient number of points along each crack segment in an array of cracks results in a large system of equations. Furthermore, constraints applied to adjoining crack segments (required to enforce displacement continuity and to eliminate mathematically induced but non-physical singularities) produce additional algebraic equations that must be satisfied explicitly. For a single crack with a complex geometry (or an array of cracks with complex geometries, where the crack face tractions are influenced by the opening displacements of the other cracks), exact solutions for the dislocation distributions typically cannot be obtained and must be approximated. A library of dislocation distributions, as described shortly, was created for the standard crack segment shapes required to analyze V-shaped, multiply kinked, and branched cracks. Dislocation distributions will be calculated as the derivatives of certain opening displacement profiles (opening shapes of the crack segments). Development of these opening displacement profiles is based on solid and fracture mechanics fundamentals, resulting in different types of series, each representing a specific characteristic of crack segment behavior. Every individual term in a series is multiplied by a weighting coefficient (or degree of freedom) whose value is determined by the solution of Eqs. (5), and each type of series is used independently in both the tangential and normal modes of deformation. To maximize computational efficiency, points where the traction conditions are enforced are allocated along crack segments according to the leverage effects of an equal and opposite opening point force pair on stress intensities (i.e., at tips, kinks, and branches). This scheme (presented in Appendix A) allows for a greater density of points to be allocated near the ends of crack segments, where accuracy is more difficult to achieve and traction errors more significantly affect results. Points must be assigned in sufficient numbers to fully capture crack behavior, and increasing the number of points is computationally less demanding than adding additional degrees of freedom. Thus, the minimum number of degrees of freedom (columns in the matrix), necessary for convergence and accuracy will be included in an analysis, while adequate points (rows in the matrix) will be assigned to capture traction behavior. This selection of points and number of terms produces an over-determined matrix that is solved by a least squares fit. This solution method has the advantage of yielding an error estimate based on the fit of the tractions induced by the computed opening displacement profile, as compared to the prescribed tractions.

3. Opening displacement series Different types of series (wedge, polynomial, and tip) comprise the library of opening displacement profiles (or dislocation distributions). Poor convergence and increased computation time arise when terms of different series compete to represent similar behavior, as demonstrated by near cancellation of extremely large equal and opposite weighting coefficients. To prevent this difficulty and to minimize computation time through efficient use of degrees of freedom, each term from each type of series must independently represent a specific physical characteristic. And, to prevent non-physical singularities, series are constrained for zero

1092

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 4. Multiply kinked crack.

opening displacement and zero slope (first derivative) at their end points. In exploratory work, additional constraints applied to the third and fourth derivatives were found to impede convergence [37]. Therefore, while series must be prohibited from non-physical singular behavior, they must retain enough ability to perform minor shaping adjustments. Details for evaluating the Cauchy singular integrals for each type of term are presented in Appendix B. Crack segments are defined as interior or exterior depending on whether they include a crack tip (i.e. Fig. 4 for a multiply kinked crack). 3.1. Wedge series A power series based on Williams’ [23] displacement equation (Eq. (6)) for a traction free wedge (Fig. 3) is developed to include the effects of wedge distortions at kinks and material singularities and has the form X ODðr; hÞ  cq ðhÞrq ð6Þ q

where the q are eigenvalues that satisfy either of the two boundary conditions (for the mating wedges) sinðqxÞ ¼ q sinðxÞ sinðqð2p
xÞÞ ¼ q sinðxÞ

ð7Þ

Only eigenvalues greater than zero automatically satisfy the opening displacement continuity condition across kinks, and only an eigenvalue whose real part is less than one will cause a stress singularity. At wedge locations, there is always a dominant q1 eigenvalue less than unity representing a circumferentially symmetric, singular stress field in the lower wedge, Mode I. There is a second singular eigenvalue, q2, when the lower wedge angle is greater than 257.4, accounting for a circumferentially anti-symmetric singular stress field, Mode II. Eigenvalues that do not cause singularities are not explicitly included in the analysis, since their effects are approximated by the higher order terms of the series. In addition, the eigenvalues q = 0 and 1 that model relative translation and rotation of upper and lower wedges also satisfy Eqs. (7) but are included in the polynomial series considered shortly. Though a wedge may be flanked by crack segments of different lengths, a1 and a2, the generic symbol a is used when writing series equations to simplify notation. The general form of a wedge series, W(t), for a crack segment emanating from a kink is given by W ðtÞ ¼ cq0


t q a

þ cq1


t qþ1 a

þ cq2


t qþ2 a

þ    þ cqn


t qþn a

ð8Þ

where n is fixed, depending on the desired accuracy, and q is either q1 or q2. The cqj are the weighting coefficients. Applying conditions to prohibit (from this series) a displacement jump and non-zero slope of the opening profile at the opposite end of the crack segment reduces the number of weighting coefficients to n
1 resulting in the final form of the wedge series (Eq. (9)).

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1093

Fig. 5. Wedge opening displacement shapes for n = 3 and q = 0.7.

W ðtÞ ¼

n
2 X

cqj


qþj t

j¼0

a


ðn
jÞ


t qþn
1 a


t qþn  þ ðn
j
1Þ ; a

nP2

ð9Þ

When both q1 and q2 apply to a wedge angle, a separate wedge series must be included for each eigenvalue, and both series must be assigned separately to all adjoining crack segments to the wedge. Therefore, an interior crack segment with a kink at each end could conceivably be assigned four wedge series, as two eigenvalues could exist for the wedge at each kink. In this case, two series, W(t), originate from the local segment coordinate origin, and two additional series written as W(a
t), would represent the wedge at the opposite end of the crack segment. If a crack is branched, wedge terms are included only when the junction contains a wedge angle greater than 180, the limit to produce singular eigenvalues. Shapes for individual components of Eq. (9) follow the two general forms depicted in Fig. 5. When j = 0, the shape has an infinite slope at the kink end, since q + j < 1 (representing the singularity as t approaches 0). Therefore, wedge series for each crack segment emanating from a kink, and thus sharing the same q value, require constraints relating their cq0 coefficients to eliminate unwanted non-physical singularities in the tractions at the kink. The final form of the constraints is given in Appendix C, and a detailed derivation is contained in [37]. The shapes for j P 1 (e.g., right example for q + j = 1.7) approximate the influence of higher order eigenvalues and do not generate a singularity at either end. 3.2. Polynomial series Polynomial series account for the q = 0 and 1 eigenvalue effects (opening and rotation) and also perform the task of manipulating the overall shape of the crack segment opening displacement profile by allowing greater flexibility in midspan deformations. A polynomial series of degree n is simply a wedge series with q = 0; therefore, the constrained polynomial series, P(t), is written as 
j n
2
t n
1
t n  X t P ðtÞ ¼ c0j
ðn
jÞ þ ðn
j
1Þ ; nP2 ð10Þ a a a j¼0 As with the wedge series, polynomial series can also be applied as P(a
t) depending on local coordinate system orientation relative to the crack segment. When applied to interior crack segments, however, one end should be restricted to n = 3 (two degrees of freedom, c00 and c01), controlling only the slope and opening displacement at that end. This restriction avoids redundant higher order terms emanating across a single crack segment from both ends (resulting in a rank deficiency in the matrix inversions) but retains the necessary shaping for the opening displacements and slopes at both ends. Fig. 6 shows the shapes of components associated with c00, c01 and c02 for n = 5. The j = 0 shape has a jump opening, associated with (t/a)0 = 1, that allows translation of an upper wedge relative to its lower counterpart. The non-physical singularity induced by this shape is removed by applying constraints to enforce continuity of opening displacements between adjoining crack segments (Appendix C). For j = 1,

1094

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 6. Polynomial opening displacement shapes for n = 5.

the slope discontinuity, (t/a)1, is another source of a non-physical singularity that is eliminated once constraints on adjoining crack segments are applied. This shape allows for relative rotation of upper and lower wedges in addition to other possible linear distortions in shear. The shapes for the higher order terms (e.g., j = 2) manipulate the overall shape of opening displacement profiles. These terms have zero displacement and slope at both ends and require no additional constraints. 3.3. Tip series Tip series incorporate crack tip behavior including the square-root singularity and higher order behavior [41]. Therefore, they are only included in opening displacement profiles of exterior crack segments. Tip series are constructed as series of powers to the order (2j + 1)/2 to exclude redundant integer terms that are already included in the polynomial series. Tip series are subject to the constraints of zero slope and displacement at the kink or branch end of a crack segment, since they always originate from crack tips (Eq. (11)). T ðtÞ ¼ 2

n
2 X j¼0

cð1=2Þj


a
t2jþ1 2 a


ðn
jÞ


a
t2ðn
1Þþ1 2 a

!
a
t2nþ1 2 þ ðn
j
1Þ ; a

nP2

ð11Þ

Eq. (11) applies to an origin located at the kink or branch end of an exterior crack segment. The factor 2 in front of the series provides coefficient consistency with the wedge series, since a crack tip is a wedge with a collapsed mating wedge (x ! 0). No additional constraints are required, since non-physical singularities are not induced by these terms. The shape for j = 0 represents the square-root singularity (Fig. 7). The shape of the j = 1 term with leading exponent 1.5 is similar to the shape depicted in Fig. 5 for the leading exponent 1.7.

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1095

Fig. 7. Tip opening displacement shape for n = 2, j = 0.

4. Calculation of stress intensity factors Stress intensity factors are calculated at both crack tip and wedge locations using the weighting coefficients corresponding to the terms inducing singular behavior. At the crack tips, the equations for the Mode I and Mode II stress intensity factors, KI and KII respectively, are rffiffiffiffiffiffi0 ½2 1   KI 2G 2p@ cð1=2Þ0 A ð12Þ ¼ 1 þ j a c½1 K II ð1=2Þ0 ½1

½2

where cð1=2Þ0 and cð1=2Þ0 are the weighting coefficients for the j = 0 tip terms representing the square-root singularity. The superscript notations, [1] and [2], denote the tangential and normal directions respectively. At the kink, generalized stress intensity factor equations are   ½2 cq1 0 =aq1 2G cos½ðq1
1Þx0 =2 ð1
q Þ ð2pÞ 1 q1 ðq1 þ 1Þ 1
KI ¼
1þj cos½ðq1 þ 1Þx0 =2 sin½ðq1
1Þx0 =2 ð13Þ   ½1 q2 0 c =a 2G sin½ðq
1Þx =2 q
1 q2 0 2 ð2pÞð1
q2 Þ q2 ðq2 þ 1Þ K II ¼

2 1þj sin½ðq2 þ 1Þx0 =2 q2 þ 1 sin½ðq2
1Þx0 =2 where x 0 = 360
x. When q2 < 1 does not exist for certain wedge angles, KII at the kink is irrelevant. A detailed derivation of Eqs. (13) is provided in Appendix D. For x = 0 and q = 1/2, these equations collapse to those for a crack tip, Eq. (12).

5. V-shaped cracks Two crack segments forming the shape of a V with angle x define a V-shaped crack (Fig. 8). In addition to verifying the accuracy of the method, this configuration was used to investigate convergence, point allocation, and the effects of eigenvalue approximation. Eigenvalues and their rational number approximations, a/b, corresponding to various values of h (x = 180
h) are listed in Table 1. These values will be used in all subsequent calculations unless specified otherwise. To verify accuracy, results from this method were compared with results from other methods [37]. Since crack tip stress intensity factors are typically the only values specified in the literature, this parameter formed the basis for comparison. However, for all examples addressed in this paper, tractions on all crack segments were evaluated and matched those prescribed in the auxiliary problem with high accuracy, and generalized stress intensity factors computed at wedges displayed rapid convergence. First, a straight crack (h = 0) divided into two crack segments with varying segment lengths and loading conditions was studied, and the results from all cases matched Irwin’s [42] exact results. And, multiple configurations of a single Vshaped crack were analyzed, including varying crack segment lengths, loading conditions (tension, shear,

1096

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 8. V-Shaped crack in an infinite plate.

Table 1 Wedge eigenvalues for varying values of h = 180
x h

q

qapprox = a/b

Relative error (%)

30 45 60 90

0.7519745 0.6735834 0.6157311 0.5444837 0.9085292 0.5122214 0.7309007

70/93 31/46 8/13 43/79 10/11 21/41 19/26

0.095 0.049 0.056 0.033 0.062 0.005 0.018

120

and biaxial), and deviation angles. Results agreed with those of other researchers who utilized a variety of methods [31,43–46]. 5.1. Effect of point allocation on convergence Enforcing tractions at only a few points will lead to nearly exact correlation in those regions, but accuracy will be compromised elsewhere. Therefore, a sufficient number of points must be assigned along crack faces to sufficiently capture traction behavior. Since computation time is proportional to the number of points, it is imperative to optimize point allocation. An error estimate is quantified through the fit of the tractions induced by the computed opening displacement profiles as compared to the prescribed tractions. This error estimate (relative root mean square

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1097

Table 2 Convergence results for varying number of points for 120 V-shaped crack with unit crack segment lengths loaded in unit tension Number of points

KI

KII

KI,kink

RRMS error Normal

Shear

25 30 40 50 100

1.2548 1.2548 1.2548 1.2548 1.2548

0.8456 0.8456 0.8456 0.8456 0.8456


0.4494
0.4494
0.4493
0.4492
0.4490

1.50 · 10
5 1.52 · 10
5 1.44 · 10
5 1.34 · 10
5 9.74 · 10
6

1.46 · 10
5 2.52 · 10
5 4.46 · 10
5 6.08 · 10
5 1.10 · 10
4

error, RRMS) is calculated by dividing the sum of the squared errors of the fit at each point by the sum of the squared applied tractions, dividing by the number of points, and then taking the square root of this quantity. Table 2 lists results obtained for varying number of points for a V-shaped crack (Fig. 8) with unit crack 1 1 segment lengths (a1 = a2 = 1) under far field unit tension ðr1 y ¼ 1; rx ¼ 0; and sxy ¼ 0Þ and x = 120. The eigenvalue required for the solution corresponds to h = 60 from Table 1, and the number of degrees of freedom (DOF) included in the analysis was 48. This number was chosen in order to compare this point allocation scheme to that of previous work [47]. This scheme achieves adequate convergence with 25 points. For 48 DOF, no less than 24 points can be placed along each crack segment, since the use of fewer points generates an under-determined matrix with more columns than rows. (Both normal and tangential tractions are calculated, resulting in a doubling factor.) To reach a compromise between time efficiency and fully capturing behavior at crack segment ends, 50 points will be used along crack segments for all subsequent examples. 5.2. Effect of number of degrees of freedom on convergence Table 3 lists results for a V-shaped crack (Fig. 8) with unit crack segment lengths under far field unit 1 1 tension and shear ðr1 y ¼ 1; rx ¼ 0; and sxy ¼ 1Þ and x = 60. The two eigenvalues required for the solution correspond to wedge angle h = 120 from Table 1. Fig. 9 shows normal traction plots along crack segment 1 as representative examples of the tractions calculated along each crack segment. The value of the prescribed traction for this case is
0.6160254. Reasonable results were obtained with only 9 DOF when only the basic crack behaviors (tip and wedge singularities and jump opening at the kink) are included with no fine-tuning. As more DOFs (up to series order n = 5 or 53 DOF) are added to allow greater flexibility in shaping the crack opening displacement profiles, convergence occurs rapidly with extremely low errors, as noted in the dramatic change of scale for the normal tractions axis. Based on these results, all subsequent analyses will use series of order n = 5 per crack segment, in both normal and tangential modes.

Table 3 Convergence results for 60 V-shaped crack with unit crack segment lengths loaded in unit tension and shear DOF

9 21 37 53

KI

KII

KI,kink

Segment 1

Segment 2

Segment 1

Segment 2


1.1652
1.1549
1.1554
1.1554

1.7196 1.8208 1.8206 1.8206

0.9087 0.9033 0.9037 0.9037


0.6459
0.7567
0.7554
0.7554


0.0146
0.0561
0.0565
0.0564

KII,kink
2.2846
2.4146
2.4062
2.4051

RRMS error Normal

Shear

0.0317 1.98 · 10
3 5.26 · 10
5 1.16 · 10
6

0.0472 1.41 · 10
3 2.50 · 10
5 2.25 · 10
6

1098

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 9. Plots of normal tractions for varying numbers of degrees of freedom along crack segment 1 for a 60 V-shaped crack with unit crack segment lengths loaded in unit tension and shear.

5.3. Effects of wedge eigenvalue approximation While most methods in the literature fix the singular exponent at kinks as q = 1/2, this method closely approximates the exact singular eigenvalues calculated from a wedge analysis. To more fully understand the effects of approximation on solution accuracy and convergence, two examples are studied. The first example is the x = 120, V-shaped crack in an infinite plate loaded in far field unit tension from the point allocation study. This problem is symmetric with only one singular eigenvalue, q1; therefore, the effects of independently varying eigenvalue approximations, a/b, can be isolated. The commonly used value of 1/2

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1099

is compared to the values 8/13 (which closely approximates the exact value, 0.61573, with a relative error of 0.06%) and 274/445 (with a relative error of only 5.5 · 10
5%). As shown in Table 4 and Fig. 10, substantial error results from using q = 1/2 in place of a more accurate approximation. Large deviations occur in the tractions, particularly at the kink, and the generalized Mode I stress intensity factor computed at the kink is inaccurate. However, crack tip stress intensity factors appear unaffected. At the other extreme, the traction plots for 8/13 and 274/445 are nearly identical, and the results with 8/13 exhibit sufficient convergence. Therefore, the additional accuracy is not worth the substantially increased computation time. (Since a sum is performed over the denominator integer value of the approximation when integrating the corresponding dislocation distribution term, the large values associated with more accurate approximations dramatically increase computation time.) The second eigenvalue study was performed using the previously described x = 60, V-shaped crack loaded in unit tension and shear from the DOF convergence study. Since this angle has two singular eigenvalues, as seen in Table 1, (q1 = 0.5122, q2 = 0.7309), the effects of inaccurately approximating, as well as ignoring, the second eigenvalue can be evaluated. Results from three cases are studied: Case 1 (as presented in the DOF study, Table 3) correctly accounts for both singular eigenvalues; Case 2 approximates both eigenvalues as q1 = q2 = 1/2; and Case 3 correctly includes the first singular eigenvalue but ignores the second. As shown in Table 5, relatively small errors occur in the Modes I and II stress intensity factors at the crack tips when the analysis is performed with incorrect eigenvalues. However, the effects on the Modes I and II generalized stress intensity factors at the kink are more noticeable. For Case 2, since the correct value of q1 is close to the value of 1/2, little deviation occurred in the Mode I generalized stress intensity factor. However, the Mode II generalized stress intensity factor is inaccurate, and the traction plots show large error at the kink (Fig. 11). When only q1 was included in the analysis (Case 3), even the K1 value at the kink was affected. In addition, no Mode II generalized stress intensity factor can be computed, and the tractions at the kink deviate to an even greater extent. It is therefore concluded that all singular eigenvalues are required in an analysis and must be close approximations of their true values. Moreover, these results further substantiate the need to analytically evaluate the Cauchy singular integrals to obtain accurate stress fields at wedges, since it is clear that results are sensitive to large approximations such as those that can occur from numerical integration. 5.4. Interacting V-shaped cracks Before attempting the case of two interacting V-shaped cracks (Fig. 12), two collinear straight cracks were analyzed for varying separation distances. The results [37] agreed with Erdogan’s exact solution [20]. Having verified accuracy for two straight cracks, two interacting V-shaped cracks with angles h = 60 were analyzed under unit tension loading. The crack segment lengths of a1 = a4 = 2 and a2 = a3 = 1 were held constant while the value of d (the distance between the cracks) was varied.

Table 4 Convergence results for varying eigenvalue of x = 120 V-shaped crack with unit crack segment lengths loaded in unit tension q

KI

KII

KI,kink

RRMS error Normal

Shear

1/2 8/13 274/445

1.2550 1.2548 1.2548

0.8452 0.8456 0.8456


0.2505
0.4492
0.4498

1.03 · 10
4 9.08 · 10
6 1.06 · 10
5

2.19 · 10
3 5.77 · 10
5 8.71 · 10
6

1100

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 10. Effects of varying eigenvalue on normal tractions of crack segment 1 for a 120 V-shaped crack with unit crack segment lengths loaded in unit tension.

Niu and Wu [36] studied interacting kinked and branched cracks by modeling deviations from main cracks as dislocation distributions. The interpolations of their graphically presented results are listed in Table 6 and may suffer inaccuracy in the second decimal place; however, the results from the present method match within this tolerance. As will be demonstrated in the coalescence example, convergence can be achieved for much smaller distances between the cracks than those presented in this example. Moreover, this method is not limited to collinear main cracks nor to equal deviations of the two cracks or segment lengths.

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1101

Table 5 Convergence results for varying eigenvalue of x = 60 V-shaped crack with unit crack segment lengths loaded in unit tension and shear Case

1 2 3

KI

KII

Segment 1

Segment 2

Segment 1

Segment 2


1.1554
1.1495
1.1519

1.8206 1.8147 1.8171

0.9037 0.9134 0.9591


0.7554
0.7457
0.6999

KI,kink

KII,kink


0.0564
0.0556
0.0565


2.4051
3.4794 N/A

RRMS Error Normal

Shear

1.16 · 10
6 0.0147 0.0344

2.25 · 10
6 0.0264 0.0392

Fig. 11. Effects of varying eigenvalue on normal tractions of crack segment 1 for a 60 V-shaped crack with unit crack segment lengths loaded in unit tension and shear.

1102

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 12. Two interacting V-shaped cracks.

Table 6 Results comparison of stress intensity factors at crack tip 2 for varying separation distance of two interacting V-shaped cracks (Fig. 12) under unit tension (a1 = a4 = 2, a2 = a3 = 1 and h = 60) d/a2

KI

KII

[36]

Present work

[36]

Present work

1.05 1.10 1.25 1.50

2.48 2.30 2.13 1.86

2.49 2.34 2.10 1.93

1.19 1.14 1.10 1.06

1.19 1.16 1.11 1.07

6. Multiply kinked and branched cracks Using the same general principles successfully applied to V-shaped cracks, the method is now applied to cracks with more than one kink and branched cracks. 6.1. Examples for multiply kinked cracks To demonstrate the accuracy of the method for multiply kinked cracks, an antisymmetric crack with two kinks (Fig. 13) was analyzed, and the results compared with accepted results from Kitagawa and Yuuki [46] and Vitek [33]. As with previous comparisons, stress intensity factors at crack tips form the basis of the comparisons, although crack face tractions matched those prescribed and wedge singularities demonstrated convergence. The loading was unit tension (/ = 90), and the central crack segment length was a2 = 2. Table 7 presents results for two values of a1 = a3 while h is varied. Next, the same configuration was analyzed under pure shear loading, (/ = 0). For comparison, the central crack segment is again a2 = 2, but the kink angle is now held constant at h = 60. Results showing agreement are given in Table 8. 6.2. Examples for branched cracks Though not presented here, results were obtained for a single branched crack and compared to those from other researchers [16,46,48], including experimental results [49]. Again, only stress intensity factors at crack tips were available. As with previous examples, tractions matched those prescribed, and when applicable, generalized wedge stress intensity factors exhibited rapid convergence. Details appear in TerMaath [37], where agreement is seen between the results of the above researchers and the present method. Next, the problem of two interacting branched cracks was studied. Comparisons were made with the results of Niu and Wu [36] under unit tension loading (Fig. 14). The first example was a symmetric case where all lengths and branch angles were held fixed (a1 = a4 = 2, a2 = a3 = a5 = a6 = 1, and b1 = b2 = 30) while the separation distance, d, was varied. Stress intensity factor results shown in Table 9 for the interacting crack tips are in agreement. The second example was non-symmetric. Lengths were equal to the previous values with the exception that now a5 = a6 = 0.1. One branch angle was fixed at

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1103

Fig. 13. Multiply kinked crack in an infinite plate (two kinks).

Table 7 Comparison of crack tip stress intensity factors for a crack with two kinks (Fig. 13) under unit tension a1/a2

h

[33]

[46]

Present work

[33]

[46]

Present work

0.1

30 45 60

1.570 1.208 N/A

1.605 1.237 0.800

1.604 1.237 0.810

0.705 0.909 0.972

0.714 0.926 0.995

0.716 0.926 0.995

0.2

30 45 60

1.686 1.283 0.806

1.720 1.311 0.763

1.695 1.260 0.762

0.829 1.085 1.099

0.839 1.102 1.124

0.827 1.061 1.126

KI

KII

Table 8 Comparison of crack tip stress intensity factors for a crack with two kinks (Fig. 13) under unit shear a3 0.10 0.20 0.60 1.0

KI

KII

[46]

Present work

[46]

Present work

0.3972 0.5681 1.0111 1.3212

0.3991 0.5687 1.0120 1.3223


0.2465
0.3497
0.6123
0.7954


0.2464
0.3498
0.6127
0.7959

b1 = 30, and the separation distance was constant at d = 0.3. The other angle, b2, was varied. The stress intensity results at the interacting crack tips are again in agreement (Table 10).

1104

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 14. Two interacting branched cracks in an infinite plate.

Table 9 Stress intensity factor comparison for two symmetric interacting branched cracks d

KI [36]

Present work

[36]

Present work

1.05 1.10 1.25 1.50

2.04 1.95 1.84 1.68

2.00 1.93 1.81 1.68

0.94 0.89 0.78 0.69

0.94 0.89 0.78 0.70

KII

Table 10 Stress intensity factor comparison for two interacting branched cracks b2 30 60

K 2I

K 2II

K 6I

K 6II

[36]

Present work

[36]

Present work

[36]

Present work

[36]

Present work

0.53 0.62

0.53 0.61

1.12 1.15

1.12 1.15

1.65 1.10

1.66 1.10


0.34
0.96


0.31
0.95

7. Crack propagation and coalescence To demonstrate the ability of the method to treat local behavior of interacting cracks, an example problem is presented that propagates and joins two initially straight cracks (Fig. 15) in an infinite plate. The loading is unit tension applied perpendicular to the cracks. To solve for the stress intensity factors, the two straight cracks are divided into two crack segments each of unit length and solved in the same manner

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1105

Fig. 15. Localized coalescence of two interacting cracks.

as V-shaped cracks. For this particular problem, the cracks are forced to propagate along the shortest path between the cracks, and possible growth at the extremes is not considered, so that the local joining growth is isolated. As the cracks begin to grow towards each other, they become V-shaped cracks. Convergence is achieved for very small ligaments (as short as 0.0774) between the two V-shaped cracks. (The main cracks are now modeled as entire crack segments, while the deviations are also modeled as individual crack segments.) Note, that for crack tip 1 (and tip 3 by symmetry arguments), KII grows faster than KI, so in reality, the crack paths may curve as the two cracks tips attempt to avoid each other. When the two V-shaped cracks coalesce, they form a multiply kinked crack. The shear stress intensity factor, KII, at the ends of the main cracks, represented by crack tip 2, is negligible, and therefore not tabulated. As this example demonstrates, the method is capable of providing converged values of stress intensity factors and opening displacement shapes for propagating cracks with very small ligaments of material between them. With this information, growth of larger crack arrays can be studied.

1106

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

8. Future work A straightforward extension of this research is to study cracks in a finite plate. This problem would be modeled as a polygon-shaped plate created by crack segments embedded in an infinite plate. For example, a rectangular plate would be modeled as four connected crack segments with four wedges of 90. In addition, interpenetration of crack surfaces is not currently prevented and should be eliminated, and friction laws governing crack face sliding could be implemented through traction criteria. Also, small-scale yielding at crack tips could be readily incorporated providing a more realistic model of crack behavior for materials exhibiting limited ductile behavior. Moreover, this method is conducive to programming in a parallel environment, and if transferred to a supercomputer, large arrays of cracks could be processed in a feasible amount of time, allowing for Monte Carlo simulation of fracture processes. For varying probabilistic characteristics (for example, crack orientations, growth properties, or material defects), statistical properties could be determined by investigating randomly assigned crack patterns. This method could be applied to curved cracks by modeling the crack as a series of straight crack segments and using the solution method for a multiply kinked crack. For propagating cracks where there is a substantial mode II component, crack growth will not be linear. This case could be handled by modeling the curved crack growth as a new crack segment attached to the initial crack with a newly formed kink. By assuming a brittle matrix material, this method could be applied to cases of cracking in composite materials. The problem of crack-bridging by fibers could be studied by modeling the effects of the fibers as point forces acting on the crack faces. T- and H-shaped cracks, that induce debond between the fiber and the matrix, could be investigated once this method is modified for bimaterial interfaces. 9. Conclusions Superposition methods based on a dislocation distribution approach require the determination of crack opening displacement profiles that satisfy the traction-free condition on crack faces. Even for single cracks of complex shape, this solution approach is mathematically cumbersome creating a barrier to the development and application of methods to analyze large arrays of interacting cracks. To overcome this difficulty, the method presented in this paper creates these displacement profiles from sets of series, each representing a physically realistic behavior of the cracks. Such an approach proves accurate and efficient. Solutions exhibited rapid convergence (with respect to prescribed tractions and stress intensity factors) using few degrees of freedom, and stability of the weighting coefficients was achieved regardless of loading conditions or crack configuration. This research shows promising applicability to a wide variety of important problems of practical interest and begins to fill the need for a means to efficiently evaluate general configurations of interacting cracks under arbitrary loading conditions. Acknowledgement S.L.P. and C.Y.H. acknowledge support under the Institute for Future Space Transport (IFST), a NASA University Institute that is monitored through NASA Glenn Research Center. The authors would also like to thank Dr. Brun Hilbert for his time in reviewing this paper.

Appendix A. Point allocation scheme For exterior crack segments, the effects of an opening point force pair on the stress intensity factors at kinks and tips are assumed to be proportional to 1=xqj 1 and 1/(1
xj)1/2 respectively. (xj is the distance mea-

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1107

sured from the kink or branch point when the crack segment length is normalized by a to unit length.) These relationships are multiplied by constants and integrated from the origin of the crack segment to the point xj. Enforcing continuity across the crack segment leads to the solution for the constants. The final point allocation equations are   1 ð3
2q1 Þ j 1
q1 1 ; 0 < xj 6 xj ¼ q1 2ð1=2Þ n 2 !2 ðA:1Þ pffiffiffiffiffiffiffiffi ð3
2q1 Þðj=nÞ
1 1 pffiffiffi 6 xj < 1 xj ¼ 1
1=2
; 2 2 2ð1
q1 Þ For interior crack segments, both ends are bounded by kinks so the effects of a point force pair on wedge q stress intensities at both ends must be included and are assumed proportional to 1=xj 1;1 and 1=ð1
xj Þq1;2 . The second subscripts on the eigenvalues designate the origin, 1, or the end of the crack segment, 2. Following the previous derivation, final equations for point allocation along interior crack segments are !1
q1 1;1 2
q1;1
q1;2 j 1 xj ¼ ; 0 6 xj 6 q1;1
1 n 2 ð1=2Þ ð1
q1;2 Þ ðA:2Þ  1
q11;2 j ð2
q1;1
q1;2 Þ n þ ðq1;2
1Þ 1 1 1
6 xj 6 1 xj ¼ 1
; ð1
q1;1 Þ 2 2

Appendix B. Integration of dislocation distributions Dislocation distributions, l, derived from the terms of the various series (Eqs. (9)–(11)), are analytically integrated according to Eqs. (3). All terms have the form (t/a)k, where the exponent k P 0 is a real number. One of five different solutions applies for Z1 and Z2 depending on the particular value of k. B.1. Solution 1: k ¼ q where 0 < q < 1 For this case, opening displacement terms of the form (t/a)q are evaluated. These terms are first converted to dislocation distribution terms by taking the derivative as  d
t q q
t q
1 lðtÞ ¼ ðH ðtÞ
H ðt
aÞÞ ¼ ðH ðtÞ
H ðt
aÞÞ
dðt
aÞ ðB:1Þ dt a a a where H(t) is the Heaviside unit step function, and d(t) is the Dirac delta function. Use of these functions allows crack opening displacements to be turned on and off at the beginning and end of each crack segment. Eqs. (3) can be solved exactly if q is a rational number, if not, q can be approximated to any desired degree of accuracy by choosing integers a and b such that q  a/b. In terms of this approximation, the final forms of the integral solutions are q
z a=b
1 X 2pika=b ðz=aÞ e2pik=b
1 Z1 ¼
e ln a a ðz=aÞ1=b e2pik=b k¼0 b
1

1=b

!


1 z
a

b
1 1=b q ð1
qÞ q
z a=b
1 X 2pika=b ðz=aÞ e2pik=b
1 Z2 ¼
e ln 1=b zðz
aÞ z a a ðz=aÞ e2pik=b k¼0

ðB:2Þ

!


1 2

ðz
aÞ

1108

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

B.2. Solution 2: k ¼ q + j where q + j > 1 and q + j is not an integer In this case, q is again approximated as a rational number a/b. A recursive procedure using the identity (t/a)(q+j)
1 =
((z
t)/a)(t/a)(q+j)
2 + (z/a)(t/a)(q+j)
2 is applied successively until the results from Eqs. (B.2) are applicable. B.3. Solution 3: k ¼ q ¼ 0 For this exponent, the corresponding dislocation distribution term is lðtÞ ¼ dðtÞ
dðt
aÞ

ðB:3Þ

Using Eqs. (3), it is trivial to evaluate Z1, and Z2 as 1 1 Z1 ¼
z z
a 1 1 Z2 ¼ 2
2 z ðz
aÞ

ðB:4Þ

B.4. Solution 4: k ¼ q ¼ 1 The dislocation distribution term for this exponent is   1 lðtÞ ¼ ðH ðtÞ
H ðt
aÞÞ
dðt
aÞ a

ðB:5Þ

The integration is again trivial and results in the solution 1 1 1

z  1 ¼ ln
Z 1 ¼ ðlnð
zÞ
lnða
zÞÞ
a z
a a a
z z
a 1 1
Z2 ¼ zðz
aÞ ðz
aÞ2

ðB:6Þ

B.5. Solution 5: k ¼ N where N > 1 is an integer For the final case, the dislocation distribution is lðtÞ ¼

N
t N
1 ðH ðtÞ
H ðt
aÞÞ
dðt
aÞ a a

ðB:7Þ

yielding the following integral solutions N Z1 ¼
a N Z2 ¼ 2 a

! N
1 X 1
z N
j
1
z N
1

z
a 1 þ ln
j a a z z
a j¼1

! N
1
z N
2

z
a
z N
1  a2  X N
j
1
z N
j
2 1 þm ln

2 j a a z a zða
zÞ ðz
aÞ j¼1

ðB:8Þ

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1109

Appendix C. Coefficient constraints ½1

½2

To develop the constraints on coefficients at kinks, the notation, cqj;k and cqj;k is used to represent the tangential and normal opening displacement coefficients, respectively, where k = 1, 2 refers to the two adjoining crack segments (Fig. 8). These constraints eliminate mathematically possible but physically inadmissible singularities in the crack face tractions and enforce continuity where two or more crack segments join. C.1. Constraints for the jump openings induced by (t/a)0 polynomial terms for V-shapes To ensure the continuity in crack opening at a kink, coefficient constraints are " ½1 #  " ½1 # c00;2 c00;1 cos x sin x ¼
½2
sin x cos x c½2 c00;2 00;1

ðC:1Þ

C.2. Constraints for the slope discontinuity associated with (t/a)1 polynomial terms for V-shapes The linear opening terms allow one wedge to rotate relative to another as well as providing for possible Poisson effects in connection with the superimposed trivial problem. The constraints for these terms are     a2 ½1 a2 ½2 ½1 ½2 c01;2 ¼
c01;1 ; c01;2 ¼
c ; x¼p a1 a1 01;1     a2 ½1 a2 ½2 p ½1 ½2 ðC:2Þ c01;2 ¼ c01;1 ; c01;2 ¼
c01;1 ; x ¼ 2 a1 a1   a2 ½2 p ½1 ½1 ½2 c01;2 ¼ c01;1 ¼ 0; c01;2 ¼
c01;1 ; x 6¼ ; p 2 a1 C.3. Constraints for traction singularities associated with (t/a)q polynomial terms for V-shapes For the terms based on the wedge eigenvalue, q1, the coefficients must satisfy   q  q1  q1 a2 a2 1 þ B a2 1 ½2 ½2 ½2 ½1 ½1 cq1 0;2 ¼ cq1 0;1 ; cq1 0;2 ¼
cq1 0;1 ¼ cq1 0;1 A a1 a1 a1 For the terms based on the wedge eigenvalue, q2, the coefficients must follow   q  q2  q2 a2 a2 1
B a2 2 ½2 ½2 ½2 ½1 ½1 cq2 0;2 ¼
cq2 0;1 ; cq2 0;2 ¼ cq2 0;1 ¼ cq2 0;1 A a1 a1 a1

ðC:3Þ

ðC:4Þ

where ð1 þ qÞ sin x cosðqðx
pÞÞ cosðpqÞ cos x cosðqðx
pÞÞ þ q sin x sinðqðx
pÞÞ B¼ cosðpqÞ A¼


ðC:5Þ

C.4. Constraints for branched cracks If a branch junction includes a wedge angle greater than 180, constraints will be applied to the two crack segments adjacent to the wedge on the corresponding (t/a)q terms using the equations developed

1110

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

Fig. 16. Branched crack configuration for constraint equations.

for V-shaped cracks. However, constraint equations for the two singular polynomial terms must be modified to include additional crack segments; therefore, the subscript i now denotes crack segment 1, 2 or 3 (Fig. 16). Constraint equations for branches consisting of more than three crack segments follow the same form. For the (t/a)0 terms, the jumps must all be compatible in assembling the three wedges leading to ½1

½1

½2

½1

½2

½2

½1

½2

½1

½2

c00;1 ¼
cos x1 c00;2 þ sin x1 c00;2
cos x2 c00;3 þ sin x2 c00;3

ðC:6Þ

c00;1 ¼
sin x1 c00;2
cos x1 c00;2
sin x2 c00;3
cos x2 c00;3

If the branch has no right angle wedges or angles equal to p (collinear segments), the constraint equations on the (t/a)1 terms are ½1

½1

½1

c01;1 ¼ c01;2 ¼ c01;3 ¼ 0 ½2

½2

ðC:7Þ

½2

c01;1 c01;2 c01;3 ¼

a1 a2 a3

Physically, these constraints specify that the angles must be preserved and that tangential stretching or compression cannot occur. Tangential behavior (associated with Poisson effects) is only allowed at right angles (p/2) or angles equal to p, as will be seen by non-zero tangential coefficients. Constraints for configurations containing such angles are summarized in Table 11. Table 11 Constraints on t1 for branched crack configurations Two right angles

½1

c01,1 a1

½1

¼

c01,2 a2

½2

c01,1 a1

½1




c01,2 a2

½1

c01,3

c01,1

a3

a1

½2

¼


One right angle

One right angle

½1

¼

½1

c01,2

c01,1

a2

a1

Branch off a main crack

½1

¼

½1

c01,3

c01,1

a3

a1

½1

¼


c01,3 a3

½2




c01,3 a3

½1

½1

c01,3 ¼ 0 ½2

c01,1 a1

½2

¼


½1

c01,2 ¼ 0 c01,2 a2

½2




½2

c01,3

c01,1

a3

a1

c01,2 ¼ 0 ½2

¼


c01,2 a2

½2




½2

c01,3

c01,1

a3

a1

½2

¼


c01,2 a2

½2




c01,3 a3

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

1111

Appendix D. Generalized stress intensity factors derivation Based on Williams [23], generalized stress intensity factors at wedges greater than 180 are formulated using polar coordinates (Fig. 3). For simplicity this derivation neglects rigid body motion. Final expressions are functions of the eigenvalues and corresponding weighting coefficients. From Williams, the stress equations (for a given q) are expressed as rhh ðr; hÞ ¼ qðq þ 1Þrq
1 F ðhÞ

srh ðr; hÞ ¼
qrq
1 F 0 ðhÞ

ðD:1Þ

where the prime denotes differentiation with respect to h. From Eqs. (D.1), it is clear that values 0 < q < 1 will cause singular stress behavior as r ! 0. The displacement equations are   1 q 0 N ðhÞ 2G ur ¼ r
ðq þ 1ÞF ðhÞ þ 1þm   ðD:2Þ 1 q 0 2G uh ¼ r
F ðhÞ þ ðq
1ÞN ðhÞ 1þm D.1. Symmetric rhh case (Mode I) The biharmonic and harmonic functions of the stress and displacement equations are F ðhÞ ¼ b2 cosððq1 þ 1ÞhÞ þ b4 cosððq1
1ÞhÞ 4 b4 sinððq1
1ÞhÞ N ðhÞ ¼ q1
1

ðD:3Þ

where b2 ¼ Cb4 ;

C¼


cosððq1
1ÞaÞ cosððq1 þ 1ÞaÞ

ðD:4Þ

and a = (2p
x)/2. The normal stress along the wedge bisector, h = 0 is rhh jh¼0 ¼ q1 ðq1 þ 1Þrq1
1 ð1 þ CÞb4

ðD:5Þ

To obtain the relationship between b4 and the coefficient corresponding to the singular wedge term, set t = r in the opening displacement profile. The equation for the normal displacement along the wedge boundary h = a can then be written as ½2

uh jh¼a ¼
cq1 0 ðr=aÞq1

ðD:6Þ

An important subtlety is that setting the displacement along a wedge boundary equal to this term causes this boundary to absorb all of the displacement. From the wedge theory formulation (and generalizing to plane stress or plane strain), the displacement equation is 2G uh jh¼a ¼ rq1 ð1 þ jÞb4 sinððq1
1ÞaÞ Equating Eqs. (D.6) and (D.7), yields  ½2  cq1 0 2G 1 b4 ¼
q a 1 1 þ j sinððq1
1ÞaÞ

ðD:7Þ

ðD:8Þ

Finally, the stress intensity factor for symmetric behavior is K I ¼ limr!0 rhh jh¼0 ð2prÞ

1
q1

ðD:9Þ

1112

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

and by substitution of the above relations, the final form is !  ½2  cq1 0 2G 1 1
q1 K I ¼ ð2pÞ q1 ðq1 þ 1Þð1 þ CÞ
q a 1 1 þ j sinððq1
1ÞaÞ

ðD:10Þ

D.2. Anti-symmetric rrh case (Mode II) In a similar manner, the generalized stress intensity factor for anti-symmetric, Mode II, behavior is obtained. The biharmonic and harmonic functions are F ðhÞ ¼ b1 sinððq2 þ 1ÞhÞ þ b3 sinððq2
1ÞhÞ
4 b3 cosððq2
1ÞhÞ N ðhÞ ¼ q2
1

ðD:11Þ

where b1 ¼ Xb3 ;

X¼


sinððq2
1ÞaÞ sinððq2 þ 1ÞaÞ

ðD:12Þ

so that the shear stress along h = 0 is srh jh¼0 ¼
q2 rq2
1 ðXðq2 þ 1Þ þ ðq2
1ÞÞb3

ðD:13Þ

To obtain the relationship between b3 and the coefficient for the singular wedge term for tangential deformation, begin with the tangential displacement along h = a and t = r, which is ½1

ur jh¼a ¼
cq2 0 ðr=aÞ

q2

ðD:14Þ

From the wedge theory formulation and generalizing to plane stress or plane strain, 2G ur jh¼a ¼ rq2 ð1 þ jÞb3 sinððq2
1ÞaÞ From these two equations,  ½1  cq 0 2G 1 b3 ¼
q2 a 2 1 þ j sinððq2
1ÞaÞ

ðD:15Þ

ðD:16Þ

Finally the stress intensity factor for anti-symmetric behavior is K II ¼ limr!0 srh jh¼0 ð2prÞ1
q2

ðD:17Þ

and the final form is obtained by substitution as K II ¼ ð2pÞ

1
q2

!  ½1  cq2 0 2G 1 q2 ððq2 þ 1ÞX þ ðq2
1ÞÞ q a 2 1 þ j sinððq2
1ÞaÞ

ðD:18Þ

References [1] Burton Jr JK, Phoenix SL. Superposition method for calculating singular stress fields at kinks, branches and tips in multiple crack arrays. Int J Fract 2000;102:99–139. [2] Bilby B, Eshelby J. Dislocations and the theory of fracture. In: Liebowitz H, editor. Fracture, an advanced treatise, vol. I. New York: Academic Press; 1968. p. 99–182. [3] Lam K, Phua S. Multiple crack interaction and its effect on stress intensity factor. Engng Fract Mech 1991;40(3):585–92.

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114 [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] [38] [39] [40] [41] [42] [43] [44] [45] [46]

1113

Melin S. Why do cracks avoid each other? Int J Fract 1983;23:37–45. Gol’dstein R, Salganik R. Brittle fracture of solids with arbitrary cracks. Int J Fract 1974;10:507–23. Han X, Wang T. Interacting multiple cracks with complicated crack surface conditions. Int J Fract 1996;82:R53–7. Li S, Mear M. Singularity-reduced integral equations for displacement discontinuities in three-dimensional linear elastic media. Int J Fract 1998;93:87–114. Hu K, Chandra A. Interactions among cracks and rigid lines near a free surface. Int J Solids Struct 1993;30(14):1919–37. Hu K, Chandra A. Interactions among general systems of cracks and anticracks: an integral equation approach. J Appl Mech 1993;60:920–8. Benveniste Y, Dvorak G, Zarzour J, Wung E. On interacting cracks and complex crack configurations in linear elastic media. Int J Solids Struct 1989;25(11):1279–93. Kachanov M, Montagut E. A simple analysis of interacting cracks and cracks intersecting a hole. Int J Fract 1989;40(3):R61–5. Erdogan F, Gupta G, Cook T. Numerical solution of singular integral equations. In: Sih G, editor. Methods of analysis and solutions of crack problems. Leyden: Noordhoff; 1973. p. 368–425. Theocaris P, Ioakimidis N. Numerical integration methods for the solution of singular integral equations. Q Appl Math 1977;35:173–83. Gerasoulis A. The use of quadratic polynomials for the solution of singular integrals of Cauchy type. Comput Math Appl 1982;8:15–22. Melin S. On singular integral equations for kinked cracks. Int J Fract 1986;30:57–65. Chen Y, Hasebe N. New integration scheme for the branch crack problem. Engng Fract Mech 1995;52(5):791–801. Hills D, Kelly P, Dai D, Korsunsky A. Solution of crack problems: the distributed dislocation technique. Dordrecht: Kluwer Academic Publishers; 1996. Westergaard H. Bearing pressures and cracks. J Appl Mech 1939;6:49–53. Koiter W. An infinite row of collinear cracks in an infinite elastic sheet. Ing Arch 1959;28:168–72. Erdogan F. On the stress distribution in plates with collinear cuts under arbitrary loads. In: Proceedings, fourth US national congress of applied mechanics, 1962. p. 547–53. Sih G. Stress distribution near internal crack tips for longitudinal shear problems. J Appl Mech 1965;32:51–8. Murakami Y, editorStress intensity factors handbook. 1st ed. New York: Pergamon; 1987. Williams M. Stress singularities resulting from various boundary conditions in angular corners of plates in extension. J Appl Mech 1952;19:526–8. Williams M. On the stress distribution at the base of a stationary crack. J Appl Mech 1957;24:109–14. Timoshenko S, Goodier J. Theory of elasticity. 3rd ed. New York: McGraw-Hill; 1970. Barber J. Elasticity. Boston: Kluwer Academic Publishers; 1992. Cotterell B, Rice J. Slightly curved or kinked cracks. Int J Fract 1980;16(2):155–69. Karihaloo B, Keer L, Nemat-Nasser S. Crack kinking under nonsymmetric loading. Engng Fract Mech 1980;13:879–88. Karihaloo B, Keer L, Nemat-Nasser S, Oranratnachai A. Approximate description of crack kinking and curving. J Appl Mech 1981;48:515–9. Bilby B, Cardew G. The crack with a kinked tip. Int J Fract 1975;11:708–12. Chatterjee S. The stress field in the neighborhood of a branched crack in an infinite sheet. Int J Solids Struct 1975;11:521–38. Lo K. Analysis of branched cracks. J Appl Mech 1978;45:797–802. Vitek V. Plane strain stress intensity factors for branched cracks. Int J Fract 1977;13(4):481–501. Hayashi K, Nemat-Nasser S. Energy release rate and crack kinking under combined loading. J Appl Mech 1981;48:520–4. Blanco C, Martinez-Esnaola J, Atkinson C. Kinked cracks in anisotropic elastic materials. Int J Fract 1998;93:387–407. Niu J, Wu M. Strong interactions of morphologically complex cracks. Engng Fract Mech 1997;57:665–87. TerMaath S. A two-dimensional analytical technique for studying fracture in brittle materials containing interacting kinked and branched cracks. Dissertation, Cornell University, 2000. Lardner R. Mathematical theory of dislocations and fracture. Great Britain: University of Toronto Press; 1974. Hirth J, Lothe J. Theory of dislocations. New York: John Wiley & Sons; 1982. Bueckner H. A novel principle for the computation of stress intensity factors. Z Angew Math Mech 1970;50:529–46. Hui C-Y, Ruina A. Why K? High order singularities and small scale yielding. Int J Fract 1995;72:97–120. Irwin G. Analysis of stresses and strains near the end of a crack transversing in a plate. J Appl Mech 1957;24:361–4. Abe H, Hayashi K, Yamamoto T. Growth path of a crack in earth’s crust. Trans Jpn Soc Mech Engrs 1985;51–465:1359–66. Isida M. Analysis of stress intensity factors of plate with arbitrary array cracks and bent cracks. Trans Jpn Soc Mech Engrs 1978;44–380:1122–33. Isida M, Nishino T. Formulae of stress intensity factors of bent cracks in plane problems. Trans Jpn Soc Mech Engrs 1982;48– 430:729–38. Kitagawa H, Yuuki R. Stress intensity factors for branched cracks in a two-dimensional stress state. Trans Jpn Soc Mech Engrs 1975;41–346:1641–9.

1114

S.C. TerMaath et al. / Engineering Fracture Mechanics 73 (2006) 1086–1114

[47] TerMaath S, Phoenix SL. Investigation of a new analytical method for treating kinked cracks in a plate. In: Fatigue and fracture mechanics. ASTM STP 1389, vol. 31. 2000. p. 331–47. [48] Isida M, Noguchi H. Formulae of stress intensity factors of branched cracks in plane problems. Trans Jpn Soc Mech Engrs 1983;49–440:469–79. [49] Theocaris P. Complex stress-intensity factors at bifurcated cracks. J Mech Phys Solids 1972;20:265–79.