A general methodology for calculating mixed mode stress intensity factors and fracture toughness of solder joints with interfacial cracks

Engineering Fracture Mechanics 131 (2014) 9–25

Contents lists available at ScienceDirect

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

A general methodology for calculating mixed mode stress intensity factors and fracture toughness of solder joints with interfacial cracks Z. Huang a, P. Kumar b, I. Dutta a,⇑, J.H.L. Pang c, R. Sidhu d a

The School of Mechanical and Materials Engineering, Washington State University, Pullman, WA 99164, USA Department of Materials Engineering, Indian Institute of Science, Bangalore 560012, India c School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore 639798, Singapore d Intel Corp., Assembly Technology Development, 5000 W. Chandler Blvd., Chandler, AZ 85226, USA b

a r t i c l e

i n f o

Article history: Received 27 January 2014 Received in revised form 30 September 2014 Accepted 7 October 2014 Available online 14 October 2014 Keywords: Joint fracture Interfacial crack Mixed-mode fracture

a b s t r a c t Solder joints in electronic packages undergo thermo-mechanical cycling, resulting in nucleation of micro-cracks, especially at the solder/bond-pad interface, which may lead to fracture of the joints. The fracture toughness of a solder joint depends on material properties, process conditions and service history, as well as strain rate and mode-mixity. This paper reports on a methodology for determining the mixed-mode fracture toughness of solder joints with an interfacial starter-crack, using a modified compact mixed mode (CMM) specimen containing an adhesive joint. Expressions for stress intensity factor (K) and strain energy release rate (G) are developed, using a combination of experiments and finite element (FE) analysis. In this methodology, crack length dependent geometry factors to convert for the modified CMM sample are first obtained via the crack-tip opening displacement (CTOD)-based linear extrapolation method to calculate the under far-field mode I and II conditions (f1a and f2a), (ii) generation of a master-plot to determine ac, and (iii) computation of K and G to analyze the fracture behavior of joints. The developed methodology was verified using J-integral calculations, and was also used to calculate experimental fracture toughness values of a few lead-free solder-Cu joints. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Most handheld electronic devices undergo thermo-mechanical cycling during service and therefore, tiny low-cycle fatigue cracks form early in the service life of the solder joints. These pre-existing cracks may propagate under a combination of tensile and shear loading (i.e. mixed-mode loading) when the package sustains an impact during a drop. In addition, the scale of plasticity at the crack tip is not the same for solder joints processed with different solder composition or joint-reflow conditions or the post-reflow thermo-mechanical excursions, which results in varied effective crack growth as well as fracture toughness [1–5]. Thus, it is necessary to develop a general methodology to evaluate the fracture resistance capability of solder joints after duly addressing the varied effective crack growth. Critical stress intensity factor, KC, and strain energy release rate, GC, are the most common parameters used to characterize the fracture toughness of materials. Unlike a crack embedded in a homogeneous body (Fig. 1a), stress-field ⇑ Corresponding author. Tel.: +1 509 335 8354; fax: +1 509 335 4662. E-mail address: [email protected] (I. Dutta). http://dx.doi.org/10.1016/j.engfracmech.2014.10.003 0013-7944/Ó 2014 Elsevier Ltd. All rights reserved.

10

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Nomenclature a C E E0 f G H J k K K0 l P r s t* T * u u UX UY w W

crack length (m) compliance (m/N) Young’s modulus (MPa) plane strain Young’s modulus (MPa) geometry factor (dimensionless) strain energy release rate (J/m2) substrate height (m) J integral value (J/m2) tangent modulus (MPa) stress intensity factor (MPa m0.5) normalized stress intensity factor (MPa m0.5) characteristic length (m) force (N) distance from the crack tip (lm) path along a J-integral contour (m) thickness of the specimen (m) stress vector acting on a contour (MPa) displacement vector (m) node displacement of finite element model after loading (m) displacement in direction x, applied as boundary condition in the finite element model (m) displacement in direction y, applied as boundary condition in the finite element model (m) strain energy density (J/m3) width of the joint (mm)

Greek symbols a Dundurs parameters (dimensionless) b Dundurs parameters (dimensionless) d relative crack-tip opening displacement (m) e Dundurs parameters (dimensionless) e_ strain strain rate (s1) / loading angle (°) C integral path for J-integral (m) l shear modulus (MPa) m Poisson’s ratio (dimensionless) h angle originating from the crack tip of finite element model (°) r normal stress (MPa) s shear stress (MPa) w mode mixity (°) Subscripts 0 initial value 1 material 1 or node 1 2 materials 2 or node 2 1a mode I component for interfacial crack 2a mode II component for interfacial crack a combined effect of adhesive and substrate c critical value I mode I component x x direction of the coordinate at the crack tip y y direction of the coordinate at the crack tip

heterogeneities in the presence of interfacial cracks (Fig. 1b and c) are also caused by discontinuity in the elastic properties of the two adjoining materials composing the interface. This elastic (or materials) mismatch is contained in the three Dundurs parameters, namely e, a and b, which are given as follows [6]:

e¼

1 1b ln 2p 1 þ b

a ¼ ðE01  E02 Þ=ðE01 þ E02 Þ

ð1aÞ ð1bÞ

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

11

Fig. 1. Schematics of a (a) crack embedded in a homogeneous body; (b) bi-material interfacial crack; and (c) interfacial crack in a sandwiched structure.

b¼

1 2

l1 ð1  2m2 Þ  l2 ð1  2m1 Þ l1 ð1  m2 Þ þ l2 ð1  m1 Þ

ð1cÞ

where E0 is the plane strain Young’s modulus [=E/(1  m2), here E is the Young’s modulus]; l and m are the shear modulus and Poison’s ratio of a material, and subscripts 1 and 2 represent the two materials forming the interface. In this report, subscripts 1 and 2 represent substrate and solder, respectively. These Dundurs parameters may be employed to theoretically determine the stress field around an interfacial crack [7,8]. Sandwich structure systems, as shown in Fig. 1c and akin to solder joints, containing interfacial cracks have been studied and formulations utilizing Dundurs parameters have been proposed to calculate stress intensity factors [9–19]. Early on, a universal relationship between the stress intensity factor(s) of the sandwich joints and the homogeneous specimens with the same geometry was developed for the condition when the thickness of the adhesive layer is negligible as compared to the other geometry length scales of the specimen [9,10]. Therefore, the stress intensity factor(s) of such adhesive joints can be derived by solving for the stress field at the crack tip of the corresponding homogeneous specimen, which is relatively easier to solve. However, such models are not appropriate for solder joints, as the thickness of a solder joint may not be negligible relative to the other dimensions; for example, thickness of a ball grid array (BGA) solder joint is comparable to the diameter of solder balls. Besides the above method, a few methodologies based on stress or displacement extrapolation method [11–13], J integral [14,15], stiffness derivative technique [16,17], etc. and utilizing finite element method (FEM) were also developed to determine the stress intensity factors. Among these, the extrapolation method is widely used due to its simplicity. Since stress values, being derived variables, are less accurate than the displacement values at the finite element nodal points, the displacement extrapolation method is relatively more accurate, and hence preferred [12]. Pang and Shi [12] utilized the displacement extrapolation method and proposed a scheme to calculate the stress intensity factor(s) as the product of the applied stress, square root of the crack length and the geometry factor(s).1 The geometry factor(s) depend only on the sample geometry, such as the thickness of the joint, width of the sample, and crack length. So, the stress intensity factor(s) of a joint can be conveniently calculated provided these geometry factors are known. However, most of the studies calculated the fracture toughness of joints based on ASTM E399 [20], ignoring the fact that the adhesive, which is solder in this case, may be an elastic–plastic material, and hence the critical crack length may be affected by effective crack growth. In addition, these studies focused on joints containing cohesive cracks in the solder rather than adhesive cracks at the interface. It is noted that the latter is more likely due to the brittle intermetallic compounds formed between the solder and substrate, which behave as potential sites for pre-cracks formation. Meanwhile, the constraint on the joint is large enough to limit deformation of the solder and cause transition from cohesive to adhesive fracture. As such, this limited deformation is responsible for the constrained plastic zone, and is a pre-requisite for near-interface failure, which is the fracture mode observed in real joints during drop [21,22]. Thus, the literature lacks a systematic methodology to calculate mode I and II stress intensity factor(s) for solder joints in the presence of interfacial crack. As mentioned earlier, determination of the mixed-mode stress intensity factor in the presence of an interfacial crack is particularly important for the reliability of solder joints due to the commonly observed fracture behavior in these joints during a drop. This study presents a general methodology to calculate stress intensity factors for a CMM solder joint [23] under mixed mode I and mode II loading and for the plane strain condition. Here, the displacement extrapolation method, as proposed by Pang and Seetoh [12], is employed. Using FEM, J-integrals were also calculated for identical conditions and compared with 1

The geometry factors are named as ‘‘non-dimensionalised stress intensity factors’’ in Ref. [18,19].

12

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 2. Schematic of CMM specimen with (a) a single interfacial crack; and (b) double interfacial cracks. In our experiments, W = 9 mm, a/W = 0.5, t = 6.35 mm.

the GC values, as calculated using KC values, to check the validity of the derived geometry factors. Also, a methodology to calculate the effective crack growth and determination of the critical crack length is presented. The developed methodology to calculate critical crack length and the derived geometry factors, f1a and f2a, were then used to calculate effective stress intensity factors in experiments. 2. Experimental and FEM details 2.1. Experimental details The methodology developed in this study primarily targets the fracture of CMM samples as shown in Fig. 2a. The experimental samples used in this study were composed of a thin (500 lm) layer of commercially available Sn–3.8%Ag–0.7%Cu (SAC387) solder between two copper substrates in the CMM configuration. The experimental details are available elsewhere [3] and here it is summarized for continuity. The copper substrates were polished up to 0.05 lm colloidal silica. To form a sharp interfacial crack, a 500 nm thick and 4.5 mm long Al thin film was deposited on Cu surface before soldering. The following reflow parameters were used: (a) reflow temperature: 260 °C, (b) dwell time: 30 s, 90 s, 180 s or 300 s, and (c) cooling rate: 10 °C/s (water quenched, WC) or 3 °C/s (air cooled, AC). Some samples were aged at 150 °C for 48 h after being soldered. Fracture tests were conducted by a servo-hydraulic test frame at ram speeds of 0.5–50 mm/s (corresponding to joint strain rates of 1–100s1), where the crack tip opening displacement (CTOD) was measured using a clip gauge attached to a pair of knife edges glued to the Cu substrate within 1–2 mm of the interfaces. The loading angle applied in this study is 0°, 30°, 45° and 75°. 2.2. Determination of stress state in fracture samples It is necessary to determine the stress-state in the sample during a typical fracture test prior to the development of FEM models. Therefore, a few double crack samples,2 as shown in Fig. 2b, were also prepared to observe the crack propagation 2

The details of double crack samples were illustrated in Ref. [3].

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

13

Fig. 3. Crack profile of a sample with the reflow condition as (i) cooling rate = 10 s1, (ii) dwell time: 30 s and reflow temperature: 260 °C. The sample was tested in as-reflowed condition at 1 s1 under far-field Mode I loading. The arrows point out some of the surface undulations. The crack tip plastic zone size (2ry) is 160 lm. The inset pictures show magnified views of the regions and the arrow indicates the boundary where the surface undulations cease.

behavior and to ascertain the role of plasticity during fracture through observation of the regions near the tip of the partially propagated crack (also termed as secondary crack). One surface of the double crack samples was mechanically polished up to 0.05 lm colloidal silica before a fracture test in order to enable observation of the crack profile. After the test, the region ahead of the partially propagated crack was observed under scanning electron microscope (SEM). The region ahead of the crack tip where plastic deformation takes place should show apparent surface undulations including short slip traces, whereas the regions without plastic deformation should remain flat [3]. Fig. 3 shows the crack profile of a sample processed with typical reflow parameters (WC, dwell time 30 s, as-reflowed condition) and tested at relatively low strain rate of 1 s1. Plasticity ahead of the crack tip resulted in slip bands and surface undulations on the polished surface of the solder joint, thus delineating the plastic zone. By measuring the regions up to which the plasticity induced surface undulations extended, the size of the plastic zone (2ry) was estimated to be 160 lm for this condition. Similarly, plastic zones for other cases were also measured, and it was concluded that the plastic zone extended 100–200 lm ahead of the tip for all test conditions. Hence, the plastic zone radius (ry) was always close to or smaller than 1/50th of both the crack length (a0 = 4.5 mm), and the specimen thickness (t = 6.35 mm). This suggests that the lateral constraints imposed by the massive copper substrates on the solder joint limits the size of the plastic zone, including a plane strain condition for all tests in this study. Two-dimensional plane strain finite element (FE) models with same dimensions as the real joints, as shown in Fig. 2, were built using the commercial software, ANSYSÒ.3 Fig. 4 shows the strain field distribution in the solder joint at an equivalent strain rate of 100 s1 at an applied stress of 70 MPa.4 Fig. 4 reveals that although the solder is inherently ductile, it experiences only limited plasticity ahead of the crack tip due to the stress triaxiality associated with the joint between two relatively stiff substrates. This FEM result is consistent with the experimental observation that the stress-state in these samples within the strain rate range of 1 to 100 s1 is plane strain. 2.3. FEM modeling details Consistent with the stress state as observed experimentally and in FE models, 2-D plane strain FE models of CMM joints were constructed. For modeling expediency, only the part between the lines PP0 and QQ0 that pass through the centers of the loading holes (Fig. 2) were built in ANSYSÒ. 6-noded triangular elements were utilized to mesh the model. The Cu substrates were modeled as elastic material for its large Young’s modulus and yield strength as compared to the solder. The solder adhesive was modeled as iso-strain hardening, bi-linear, elastic–plastic material with the following room temperature yield strength (rYS) and tangent modulus (ktangent) [24]:

rYS ¼ 76:54e_ 0:07 strain

ð2aÞ

ktangent ¼ 201:53e_ 0:14 strain

ð2bÞ

where e_ strain is the strain rate. In FE models, the effects of strain rate on the solder was represented by using the value of yield strength and tangent modulus at the corresponding strain rate as calculated using Eq. (2). The material constants used in FE modeling are summarized in Table 1. FE models with crack-length over joint-width ratio (a/W) ranging from 0.45 to 0.6, at increments of 0.05, were built. A sharp crack-tip was simulated by constructing a line at a very small angle h (tan (h) = 1.6  102) originating from the crack 3 4

The development and solutions of FE models will be discussed in detail in the next sub-section. It was observed from the experiments that most of the samples failed well below a stress value of 70 MPa.

14

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 4. FE model with the same dimension as the real joint showing very limited plastic strain ahead of the crack tip. Some plastic strain is also produced at the end of the unbroken ligament of the solder joint. The equivalent strain rate for the simulation was 100 s1 (yield strength of solder = 106 MPa) [3].

Table 1 Various material constants for the model.

Young’s modulus (GPa) Yield strength (MPa) Tangent modulus (GPa) Poison’s ratio

Copper

SAC387

111 [25] No plasticity in Cu

55 [24] 106 at 100 s1 [24] 0.384 [24] 0.4 [25]

0.3 [25]

Fig. 5a. Schematic of the crack tip and the divided areas around the crack tip, where tan h = 1.6  102.

tip at the interface between the adhesive and the substrate (Fig. 5a). Because of the existence of this interfacial crack (where the crack face is not parallel to the ‘‘horizontal’’ x-axis), the area around the crack tip cannot be directly meshed by the conventional singular finite elements. Therefore, the region near the crack tip was divided into several smaller areas, radially distributed around the crack tip (Fig. 5a). Then, each small area was mapped meshed, separately. In order to minimize the total number of degrees of freedom (DOF) to reduce the computation time, the mesh was made coarser as it got further away from the crack tip. Nevertheless, the element sizes were small enough to ensure numerical convergence and mesh-size independence of the results. The boundary conditions are shown in Fig. 5b. Point A was fixed in both x and y directions. Line AD was constrained in y direction, i.e. displacement of this line in y direction, UY, was equal to zero. Displacements in direction x and y (UX and UY, respectively) were applied on Line BC. The loading mode mixity depends on the ratio of the applied UX and UY. For example, UY was set to be a constant value and UX was set to be equal to 0 for Mode I loading; while UY was set to be equal to 0 and

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

15

Fig. 5b. Boundary conditions and displacements applied to FE models.

UX was set to be a constant value for Mode II loading. CMM models containing varied crack lengths were loaded in both Mode I and Mode II. Fig. 6 shows typical examples of the displacement and the von Mises stress fields for one of the FEM runs. The von-Mises stress field, showing kidney shaped stress contours with maximum near the crack tip, as well as the displacement field is consistent with the physics of the problem. 3. Calculation of K and GC and determination of geometry factors Fig. 7 shows force vs. CTOD behavior of two CMM solder joints. Both solder joints were processed with same reflow parameters and have the same post-reflow thermo-mechanical history; nevertheless, the two joints show two different but typical fracture behaviors that are usually observed in these solder joints. Fig. 7a shows linear force vs. CTOD behavior of the sample tested at a high strain rate (100 s1), implying limited ‘‘effective crack growth due to plastic dissipation prior to catastrophic failure’’5 as well as low fracture toughness. On the other hand, when the sample was tested at a low strain rate of 1 s1, as shown in Fig. 7b, significant nonlinearity between the force and the CTOD occurred. This indicates considerable effective crack growth, resulting in relatively higher fracture toughness. However, it should be noted that the fracture was catastrophic (i.e. unstable) in both cases. Fracture toughness can be expressed as the critical strain energy release rate, GC for such catastrophic failures [26]. The strategy utilized for computing GC comprised two steps: (i) estimation of the actual crack length prior to the onset of the unstable failure (i.e. critical crack length, ac) by accounting for the effective crack growth (which is described in Section 4), and (ii) application of the proper equations to calculate K1a and K2a, which are the stress intensity factors for interfacial fracture of adhesive joints in mode I and II, for the CMM joints containing an interfacial crack, respectively. Here, K1a and K2a represent the real and imaginary components of the complex stress intensity factor for an interfacial crack, and are given by [19]: pffiffiffiffiffiffi / pac ac K 1a ðCMMÞ ¼ P cosWt f 1a ðW Þ pffiffiffiffiffiffi P sin / pac ac K 2a ðCMMÞ ¼ f 2a ðW Þ Wt

ð3Þ

where P is the critical load at which the unstable fracture initiates, / is the loading angle, ac is the critical crack length (equal to the initial crack length a0 plus the effective crack growth Da), t is the thickness of the joints, and f1a and f2a are the geometry factors for Mode I and II, respectively, that depend on the sample geometry and the crack length. Following the calculations of K1a and K2a, GC and the mode mixity, w, can be calculated as [10]:

GC ¼

K 21a þ K 22a 2

E0a cos h ðpeÞ

w ¼ tan1

K 2a K 1a

ð4Þ

ð5Þ

Here, e is a Dundurs parameter given by Eq. (1a), and Ea0 is the combined Young’s modulus and is given as:

E0a ¼ 2

5

E01  E02 E01 þ E02

The term ‘‘effective crack growth’’ will be short for ‘‘effective crack growth due to plastic dissipation prior to catastrophic failure’’ in this paper.

ð6Þ

16

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 6. FEM solution of the model with a/W = 0.5, UX = 0 and UY = 3 lm. (a) von-Mises stress field of the complete model; (b) von-Mises stress field near the crack tip; and (c) the displacement field.

where E01 and E02 are the plane strain Young’s moduli [=E/(1  m2), where E is the Young’s modulus] of the substrate and the adhesive, respectively. It is noted that using Eq. (5), in conjunction with Eq. (3), will result in a mode mixity of w = 0° when the loading angle / = 0°. Strictly speaking, interfacial fracture at a loading angle of 0° does not necessarily correspond to pure mode I (i.e. in general, w – /), because the stress field near the tip of an interfacial crack involves an oscillatory singularity (rie) [27]. However, in the present work, the use of Eqs. (3) and (5) can be justified as follows. For joints with very small Dundurs parameters a, b and e (given by Eq. (1)), the oscillatory singularity as well as the deviation from pure mode I (w = 0°) will be very small [28].

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

17

Fig. 7. Example ‘‘force-crack tip opening displacement’’ plots showing (a) limited effective crack growth, and (b) substantial effective crack growth [3].

In that case, Eqs. (3) and (5) can be used to estimate the stress intensity factors and mode mixity [27,29]. In this work the e, a and b of SAC387/Cu joints are calculated as 0.002755, 0.45 and 0.00865, respectively. For such small values of a and b, the deviation from the mode I mode mixity for / = 0° is less than 5° [28], enabling the use of Eqs. (3) and (5). In order to use Eq. (3) for calculating stress intensity factors, it is necessary to calculate the geometry factors. Two different methods, namely CTOD extrapolation and J-integral, were used to calculate the geometry factors and they are explained in the following two subsections. 3.1. CTOD extrapolation method Fig. 8 shows a schematic of a crack before and after loading. For this geometry of the crack, the displacement field at a distance r behind the crack tip, along the crack face is given by [12,27]:

dy þ idx ¼

 r 1=2 r ie 8 ðK 1a þ iK 2a Þ 0 l ð1 þ 2ieÞEa cos hðpeÞ 2p

ð7Þ

where dx and dy are the relative CTODs in the x and y directions, respectively, between two key points, whose distances from the origin O before loading were equal to r, Ea0 is the combined Young’s modulus, K1a and K2a are the real and imaginary components of the stress intensity factors, respectively, l is a characteristic length, which weakly affects the stress field of the

18

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 8. A schematic of crack tip showing the relative displacement between two key points before (solid line) and after loading (broken line). Circular points show the positions of key-points before loading and the square points show their positions after loading. The relative displacement between these keys points (1 and 2) in x and y directions are equal to dx and dy, respectively.

region near the crack tip. The choice of l is somewhat arbitrary. It may be chosen either based on the sample geometrical length parameters in the plane of the crack, such as crack length and joint thickness, or a material determined length parameter, such as the size of plastic zone ahead of the crack tip [27]. Eq. (7) can also be rewritten as [12]:

K 1a ¼ K 2a ¼

A cos ðe ln ðrlÞÞþB sin ðe ln ðrl ÞÞ D

ð8aÞ

B cos ðe ln ðrlÞÞA sin ðe ln ðrl ÞÞ D

where

A ¼ dy  2edx ; B ¼ dx þ 2edy ; pffiffiffiffi r

D ¼ E0 cos8hðpeÞ

ð8bÞ

2p

a

The displacements of node 1 in x and y directions after loading were solved as u1,x and u1,y, respectively and that for node 10 were solved as u10 ,x, and u10 ,y, respectively. Hence, the relative displacements between this pair of nodes in x and y direction, dx and dy can be calculated as:

dx ¼ u1;x  u10 ;x

ð9Þ

dy ¼ u1;y  u10 ;y

For simplicity, the joint thickness (=0.5 mm) was taken as the characteristic length, (l).6 Now, K1a and K2a can be calculated by substituting the nodal displacement values, calculated using FEM, into Eqs. (8) and (9). On the other hand, P is the total reaction force on the body, corresponding to the measured load in Eq. (3), which was calculated as:

P ¼ Py ðMode IÞ or Px ðMode IIÞ

ð10aÞ

where

Px ¼

R R R R syx dxdy rxx dxdy 2H R WðyÞ R WðyÞ R px ¼ t þ 2Ht 2H

X

Whole Body

Py ¼

X

2H

dy

2H

WðyÞ

dxdy

R R R R ryy dxdy sxy dxdy 2H R WðyÞ py ¼ t þ 2Ht R2H RWðyÞ

Whole Body

2H

dy

2H

WðyÞ

ð10bÞ

dxdy

Here, r and s are the nodal stresses on finite element nodes of the model. W(y) and H are illustrated in Fig. 5b. Px and Py are estimated as the summation of average shear force and normal force in x and y directions, respectively. By substituting K1a, K2a and P into Eq. (3), f1a and f2a can be derived. For given crack-tip stress intensity factor values (K1a, K2a), d, r and e give the displacement, stress and strain at a given distance from the crack tip r (Eqs. (7), (8a), (8b), (9)). When the d, r and e values are obtained from FEM, and the K values are back-calculated, they represent ‘nominal’ stress-intensity factor values (K1a0 , K2a0 ). Analogously, using Eq. (3), ‘nominal’ values of the geometry factor (f1a0 and f2a0 ) at a given distance r from the crack-tip. These K0 or f0 values rise slightly with 6 It is noted that when the characteristic length l is changed from 0.5 mm to 4.5 mm (increased by 800%), the f1a and f2a values increased by only 0.04%, indicating that the choice of the characteristic length exerts negligible effect on the values of the geometry factors. This is consistent with other reports [12,27].

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

19

Fig. 9. (a) f1a0 at different r behind of the crack tip for the model a/W = 0.5 under Mode I loading; (b) f2a0 at different r behind of the crack tip for the same model as in (a) but under Mode II loading. The region left to solid vertical line, I, and the data points encircled by the broken rectangles are affected by the plasticity ahead of the crack tip.

decreasing r, and may be extrapolated to r = 0 to find the actual crack-tip stress intensity factor or analogously, the crack-tip geometry factors (f1a and f2a), Ref. [12]. Fig. 9a and b show the f1a0 and f2a0 at different r behind the crack tip for the model with a/W = 0.5 under Mode I loading. Here r/a = 0 represents the crack tip. Within the plastic zone (r/a 6 0.03), the combination of inelastic effects and limited mesh-fineness to adequately handle singularity cause anomalous behavior (Region I in Fig. 9a) [27]. However, away from the plastic zone, a linear relationship between the geometry factors and the distance r was established as r increased. By extrapolating the K–r graph back to r = 0, the geometry factor f1a was calculated to be equal to 0.765 for a/W = 0.5. Since /, and hence sin(/), are equal to zero in mode I loading, f2a cannot be derived by using Eq. (3) from mode I loading results. Therefore, f2a was derived under mode II loading condition for which the result is shown in Fig. 9b. Similarly, application of the above methodology to a whole range of CMM FE models containing different crack lengths gives the following for the geometry factors:

   2  3 f 1a ¼ 0:90881  0:63489 Wa þ 0:46791 Wa þ 0:45339 Wa    2  3 f 2a ¼ 0:07525 þ 2:89730 Wa  5:97820 Wa þ 4:22990 Wa

ð11Þ

Eq. (11) is valid for 0.45 6 a/W 6 0.6. 3.2. J-integral method To ascertain the validity of the CTOD-based extrapolation method, J-integral method was also employed to calculate the geometry factors f1a and f2a for the CMM joints. Along a path enclosing the crack tip, the J-integral was calculated as:

20

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 10. Schematic showing several different paths defined near the crack tip for J-integral calculation.

J¼

Z  C

 @~ u T wdy  ~ ds @x

ð12Þ

where x, y are Cartesian coordinates with origin at the crack tip; x being parallel to the crack face and y being perpendicular to the crack face; s is the path along a contour U; ~ u is the displaceT is the stress vector (or, traction) acting on the contour; ~ ment vector and w is the strain energy density. J-integral is independent of the integral path near the crack tip [26]. As shown in Fig. 10, the J-integral was calculated over several contours for each FE model.7 Fig. 11 is an example plot showing the dependence of J-integral value on the length of the integral path. As mentioned earlier, the stress and displacement values calculated using FEM are prone to errors very close to the crack tip (i.e. the origin O), resulting in the scatter for path lengths less than 9.6  104m (Fig. 11). Beyond that, the J-integral remains relatively constant, and start decreasing slightly when the integral path became larger than 3.5  103m. This is most likely because for longer path lengths, a significant proportion of the paths reside in the copper substrate (as opposed to the solder), complicating the computed value of J integral. For integral path lengths of 1–3.5  103m, J-integral remained roughly constant, displaying very slight fluctuation (±0.001 J/m2) from a constant value of 2.513 J/m2. This path independence suggests that over this range, the computed J-integral value is valid. Therefore, for the test condition shown in Fig. 11, a J-integral value of 2.513 ± 0.001 J/m2 was ascertained. Similarly, J-integral values for all other test conditions were calculated. For small scale yielding, the J-integral is equal to the GC value [26], i.e.:

J  GC ¼

K 21a þ K 22a E0a

2

cos h ðpeÞ

ð13Þ

Under Mode I loading, K2a can be ignored compared with K1a; while K1a can be ignored under Mode II loading condition [18,19]. Thus, Eq. (13) can be simplified as: For Mode I loading:

K 1a >> K 2a ;

K 1a 

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 JE0a cos h ðpeÞ

ð14aÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 JE0a cos h ðpeÞ

ð14bÞ

For Mode II loading:

K 2a >> K 1a ;

K 2a 

Combined with Eq. (3), f1a and f2a can be written as:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pffiffiffiffi : Mode I loading JE0a cos h ðpeÞ P cosWt / pa qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pffiffiffiffi : Mode II loading ¼ JE0a cos h ðpeÞ P sinWt / pa

f 1a ¼ f 2a

ð15Þ

Here, the force P is calculated by integrating the stress at the top edge as given by Eq. (10). The geometry factors derived using J-integral method for the FE model with a/W = 0.5 is shown in Fig. 12. The geometry factor values are very close (within 2.5%) to the ones calculated from the CTOD-extrapolation method. This essentially validates the calculation of the f1a and f2a using Eq. (11). 7 Strictly speaking, the entire contour for J calculation should be drawn through one single material (i.e. a homogeneous medium) [30]. Since this study deals with an interfacial crack, though, it is not possible to draw a contour through one medium. However, plotting the line integral as a function of path length, it was found that J remains roughly constant once the path length exceeds a minimum value, even for a bi-material contour. Although a theoretical justification for this is elusive, since the values remain constant for large contours, they appear to be characteristic of the sample, and hence provide a working estimate of the crack driving force.

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

21

Fig. 11. Variation in J-integral values with the total length of the integral path for the a/W = 0.5 FE model. A displacement of 1 lm was applied in Y-direction to create Mode I loading. The stress and displacement values calculated using FEM are prone to errors very close to the crack tip (i.e., the origin O), resulting in the scatter for path lengths of 9.6  104m.

Fig. 12. A comparison between the geometry factors derived using J-integral method and those derived using CTOD-extrapolation method for the a/W = 0.5 FE model. The two sets of results are close to each other (2.5% difference).

4. Calculation of critical crack length ac It is essential to calculate the critical crack length, accounting for any effective crack growth, in order to use Eq. (11). A methodology based on ASTM standard E561-05 [31], was employed to accomplish this. The methodology accounts for any effective crack growth, which can be significant if, as in the case of solder joints, significant plasticity precedes catastrophic fracture. This involved comparing the compliance (C = d/P) just prior to unstable fracture from the force vs. CTOD plot (P vs. d) with the initial compliance to determine the fractional increase in compliance (DC/C0), as shown in Fig. 7. Therefore, a master plot is needed to correlate the fractional compliance change (DC/C0) and the fractional crack length increase (Da/a0), which is developed as follows:  By FEA, force (P) vs. CTOD (d) was plotted for models with a/W values of 0.5, 0.55, 0.6 and 0.65. The simulation was conducted under far-field mode I loading, at the strain rate of 100 s1 (Fig. 13a).  From these P–d plots, for the models containing different crack lengths, the compliances C at a small-level force, used as 800N8 in this study, were calculated. The fractional changes in the compliance were then plotted against the fractional changes in the crack length (DC/C0 vs. Da/a0), as shown in Fig. 13b, using the following two equations: 8 The ‘‘small force’’ applied is arbitrary. In this study, experiments show that when the force is below 800N, P and d obey linear relationship, implying for this force level the effective crack growth can be ignored. So at this force level, the crack length is very close to the initial crack length.

22

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 13. (a) Force (P) vs. CTOD (d) for models with a/W values of 0.5, 0.55, 0.6 and 0.65. The simulation was conducted under far-field mode I, at the strain rate of 100 s1; and (b) master plot correlating the fractional compliance change (DC/C0) and the fractional crack length increase (Da/a0) [3].

DC=C 0 ¼ ðC  C 0 Þ=C 0 ¼ Da=a0 ¼ ða  a0 Þ=a0

"    #,  d d d  P P a0 P a0

ð16aÞ ð16bÞ

where C0 and a0 are the compliance and crack length of the model containing smallest (or initial) crack size, which is 4.5 mm (a/W = 0.5) in this study; C and a correspond to the model with a/W larger than 0.5 (a/W = 0.55, 0.6, 0.65). Although the master plot was developed under mode I loading, it can be also used to calculate ac for the mixed mode loading. This is because of the fact that the Mode I compliance (dI/PI) only depends on the crack length (a0 + Da) and sample geometry and it does not change with the loading condition. To calculate the mixed mode compliance, the experimental P–d plot needs to be replotted in the form of PI vs. dI plot, where PI and dI are the mode I components of the applied load and displacement:

PI ¼ P measured cos h dI ¼ dmeasured

ð17Þ

where h is the loading angle. Since the clip-gauge was affixed to the joint along the center-line of the specimen, it always measured displacement normal to the joint (i.e. dI). The measured P, on the other hand, is at an angle h to the center-line. Thus, DC/C0 on the PI vs. dI plot can be measured. From the master plot, as shown in Fig. 13b, Da/a0 and hence ac can be determined corresponding to a measured DC/C0 at which unstable fracture initiated. Then, f1a and f2a can be calculated by substituting the value of ac in Eq. (11). By substituting the values of f1a, f2a and ac into Eqs. (3)–(5), K, GC and w can be determined. The procedure for computing Gc is summarized in the form of a flow chart in Fig. 16. 5. Experimental illustration CMM solder joints with a0/W = 0.5 and different reflow, aging and testing conditions were tested. Their GC values were calculated by the methodology explained in the previous sections. The experimental results are shown in Fig. 14. The details

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

23

Fig. 14. GC values of CMM solder joints with different reflow and aging preconditions under various loading conditions [3,4].

Fig. 15. (a) Contours lines of J-integral; and (b) the comparison of GC values for different reflow and aging conditions in Fig. 14 with the FEM calculated Jintegral values for the same critical force. It is shown that these two sets of values are consistent with each other.

of the experimental results are reported elsewhere [3,4], where it was also shown that the limited extent of crack-tip plasticity resulted in small-scale yielding conditions to prevail in all cases. FEM analysis was performed to check the validity of this methodology by comparing the experimental GC data with the FEM calculated J-integral values. J-integral values at var-

24

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

Fig. 16. Flow chart summarizing the methodology to calculate GC.

ied force levels were derived for the models with different crack lengths, and then were plotted as contours lines (Fig. 15a). It was observed that J-Integral depended only weakly on the CTOD as compared to the load. This is attributed to the small scale yielding at the crack tip, implying that force is the primary parameter determining J-integral values. The GC values for different reflow and aging conditions, as shown in Fig. 14, were then compared with the FEM calculated J-integral values of a/W = 0.5 for the same critical force and the results are shown in Fig. 15b. Fig. 15b shows that these two sets of values almost coincide with each other, implying that the above mentioned methodologies can be used to calculate the fracture behavior of solder joints. 6. Conclusions A methodology for studying the fracture behavior of the solder joints with an interfacial (i.e. adhesive) crack and under the plane strain condition was developed in this study. By utilizing the CTOD based linear extrapolation method, the dependence of geometry factors f1a and f2a on crack length was determined; the result of which was then also confirmed by the J integral method. Also, a master plot correlating the change of the fractional compliance and the effective crack extension was developed. Based on this master plot, the critical crack length ac can be calculated for a fracture test. Once the critical crack length was estimated, the corresponding geometry factors f1a and f2a were calculated to determine the stress intensity factors and the critical fracture toughness. Although the geometry factors and the master plot were derived only for the CMM solder joints, the methodology presented in this study can be universally applied to develop the specific geometry factors and master plot, and hence obtain GC values, for elasto-plastic joints of other geometries under small-scale yielding conditions. Acknowledgements This research was supported by National Science Foundation (NSF) Grant # DMR-0939392 and Semiconductor Research Corporation (SRC) contract # 2008-KJ-1855. Partial support from Strategic Research and Development Program (SERDP) contract # W912HQ-10-C-0041) is also acknowledged gratefully. References [1] Kumar P, Dutta I, Sarihan V, Frear DR, Renavikar M. Effects of strain rate and aging on deformation and fracture of Sn–Ag–Cu solder joints. ITHERM 2008. 11th. p. 660–7. [2] Kumar P, Huang Z, Dutta I. Roles of process and microstructural parameters of mixed mode fracture of Sn–Ag–Cu solder joints under dynamic loading conditions. In: Proceedings of InterPACK2009, Paper no. 89205, IEEE/ASME CD-ROM. [3] Huang Z, Kumar P, Dutta I, Pang JHL, Sidhu R, Renavikar M, et al. Fracture of Sn–Ag–Cu solder joints on Cu substrates: I. Effects of loading and processing conditions. J Electron Mater 2012;41:375–89.

Z. Huang et al. / Engineering Fracture Mechanics 131 (2014) 9–25

25

[4] Kumar P, Huang Z, Dutta I, Sidhu R, Renavikar M, Mahajan R. Fracture of Sn–Ag–Cu solder joints on Cu substrates: II. Fracture Mechanism Map. J Electron Mater 2012;41:412–24. [5] Hayes SM, Chawla N, Frear DR. Interfacial fracture toughness of Pb-free solders. Microelectron Reliab 2009;49:269–87. [6] Dundurs J. Edge-bonded dissimilar orthogonal elastic wedges. J Appl Mech 1969;36:650–2. [7] Suo Z, Hutchinson JW. Interface crack between two elastic layers. Int J Fract 1990;43:1–18. [8] Rice JR, Suo Z, Wang JS. Mechanics and thermodynamics of brittle interfacial failure in bimaterial systems. In: Rühle M, Evans AG, Ashby MF, Hirth JP, editors. New York: Pergamon Press; 1990. p. 269–94. [9] Suo Z, Hutchinson JW. On sandwich test specimen for measuring interface crack toughness. Mater Sci Engng A-Struct 1989;107:135–43. [10] Majumdar BS, Ahmad J. Fracture of ceramic-metal joints: zirconia/nodular cast iron system in. Metal-ceramic joining. In: Kumar P, Greenhut VA, editors. TMS; 2009. p. 67–97. [11] Yuuki R, Liu JQ, Xu JQ, Ohira T, ONO T. Mixed mode fracture criteria for an interface crack. Engng Fract Mech 1994;47:367–77. [12] Shi XQ, Zhang XR, Pang JHL. Determination of interface fracture toughness of adhesive joint subjected to mixed-mode loading using finite element method. Int J Adhes Adhes 2006;26:249–60. [13] Yuuki R, Xu JQ. Boundary element analysis of dissimilar materials and interface crack. Comput Mech 1994;14:116–27. [14] Zavattieri P, Hector L, Bower AF. Determination of the effective mode I toughness of a Sinusoidal interface between two elastic solids. Int J Fract 2007;145:167–80. [15] Chow WT, Beom HG, Atluri SN. Calculation of stress intensity factors for an interfacial crack between dissimilar anisotropic media, using a hybrid element method and the mutual integral. Comput Mech 1995;15:546–57. [16] Hamoush SA, Salami MR. A stiffness derivative technique to determine mixed mode stress intensity factors of rectilinear anisotropic solids. Engng Fract Mech 1993;44:297–305. [17] Toparli M, Aksoy T. Fracture toughness determination of composite resin and dentin/composite resin adhesive interfaces by laboratory testing and finite element models. Dent Mater 1998;14:287–93. [18] Pang JHL. Mixed mode fracture analysis and toughness of adhesive joints. Engng Fract Mech 1995;51:575–83. [19] Pang JHL, Seetoh CW. A compact mixed mode (CMM) fracture specimen for adhesive bonded joints. Engng Fract Mech 1997;57:57–65. [20] ASTM E399-08. Standard test method for linear-elastic plane-strain fracture toughness KIC of metallic materials. Annual book of ASTM standards 04.01. p. 491–524. [21] Xu Y, Ou S, Tu KN, Zeng K, Dunne R. Measurement of impact toughness of eutectic SnPb and SnAgCu solder joints in ball grid array by mini-impact tester. J Mater Res 2008;23:1482–7. [22] You T, Kim Y, Jung W, Moon J, Choe H. Effect of surface finish on the fracture behavior of Sn–Ag–Cu solder joints during high-strain rate loading. J Alloys Compd 2009;486:242–5. [23] Arcan M, Hashin Z, Voloshin A. A method to produce uniform plane- stress states with applications to fiber-reinforced materials. Exp Mech 1978;28:141–6. [24] Long X, Guduru R, Dutta I, Sarihan V, Frear DR. Deformation behavior of Sn3.8Ag0.7Cu solders at intermediate strain rates: effect of microstructure and test conditions. J Electron Mater 2008;37:189–200. [25] Nogita K. Stabilisation of Cu6Sn5 by Ni in Sn–0.7Cu–0.05Ni lead-free solder alloys. Intermetallics 2010;18:145–9. [26] Hertzberg RW. Deformation and fracture mechanics of engineering materials. 3rd ed. Wiley; 1989. [27] Hutchinson JW, Suo Z. Mixed mode cracking in layered materials. Adv Appl Mech 1991;29:63–191. [28] Evans AG, Dalgleish BJ, He M, Hutchinson JM. On crack path selection and the interface fracture energy in biomaterial system. Acta Metall 1989;37:3249–54. [29] Hutchinson JW, Mear M, Rice JR. Crack paralleling an interface between dissimilar materials. J Appl Mech 1987;54:828–32. [30] Knowles JK, Sternberg E. On a class of conservation laws in linearized and finite elastostatics. Arch Ration Mech Anal 1972;44:187–211. [31] ASTM E561-05. Standard test method for K–R curve determination. Annual book of ASTM standards 04.01. p. 593–611.