Intergranular cracking simulation of the intermetallic compound layer in solder joints

Computational Materials Science 79 (2013) 1–14

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Intergranular cracking simulation of the intermetallic compound layer in solder joints Tong An, Fei Qin ⇑ College of Mechanical Engineering and Applied Electronics Technology, Beijing University of Technology, Beijing 100124, PR China

a r t i c l e

i n f o

Article history: Received 11 March 2013 Received in revised form 5 May 2013 Accepted 27 May 2013 Available online 28 June 2013 Keywords: Intergranular crack Cohesive interface element Finite element Intermetallic compounds Polycrystalline materials

a b s t r a c t This paper presents a grain level finite element model to simulate the cracking behavior of the intermetallic compound (IMC) layer in solder joints. The grain microstructure of the IMC layer is explicitly included in the model by Voronoi tessellations. Cohesive interface elements with a coupled cohesive law are embedded along the grain boundaries to simulate microcrack initiation, propagation and coalescence in the IMC layer. A model with a Weibull distributed grain interfacial strength is adopted to account for randomly distributed grain boundary defects. The average thickness of the IMC layer and the wavelength and the roughness of the waved solder/IMC interface are used to characterize the IMC microstructure. Using the numerical approach developed, the effects of the grain shape, the randomly distributed grain boundary defects, the thickness of the IMC layer and the morphology of the solder/IMC interface on the microcrack patterns and on the overall response of solder joints are investigated. The results indicate that the overall mechanical strength is not sensitive to the grain shape, but the microcrack pattern and the crack path depend heavily on the grain shape. In the model containing randomly distributed grain boundary defects, the weak grain interface plays a critical role in the overall strength and the crack path of the model. The average thickness of the IMC layer has the greatest impact on the overall strength and the failure mode of the solder joint. The wavelength and the roughness of the solder/IMC interface have little impact on the overall strength but do have an impact on the failure mode of the solder joint. The predicted failure mode agrees well with the experimental observation in solder joints. The presented approach is feasible for simulating microcracking and the failure behavior of the IMC layer in solder joints and other quasi-brittle polycrystalline materials. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Solder joints serve as electrical connections and mechanical supports in electronic packages. The failure of the solder joint can cause electronic devices to malfunction. The solder joint is usually fabricated by the reflow process, in which a liquid solder alloy reacts with a solid metal substrate, such as a Cu pad, to form a thin layer of intermetallic compounds (IMCs) at the bonding interface. The IMC layer is essential to the mechanical strength of the metallurgical bond in the solder joint [1]. However, an excessively grown IMC layer, which is thick and has a rough solder/IMC interface [2], becomes vulnerable and severely degrades the mechanical strength of the solder joint [1,3,4]. To investigate the influence of the IMC layer on the mechanical behavior of solder joints, one approach is to treat the IMC layer as an elastic continuum layer that has a thickness and different mechanical properties compared with the solder [5,6]. Another approach is to treat the IMC layer as an interface without thickness, ⇑ Corresponding author. Tel.: +86 1067392760. E-mail address: [email protected] (F. Qin). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.05.044

and the behavior of the interface, such as crack initiation and propagation, is modeled by using the cohesive zone model (CZM) [7]. Although these approaches are capable of analyzing the stress and the strain in solder joints with a needle shape, hemisphere shape [8] or scallop shape solder/IMC interface [9], they do not consider the grain microstructure and randomly distributed grain boundary defects in the IMC layer. As a result, they cannot model the microcracking of the IMC layer, and their results do not agree well with the experimental findings. The IMC layer is a polycrystalline quasi-brittle material, and its mechanical behavior depends on the microstructural characteristics, such as the grain shape and the grain orientation, and strongly depends on grain boundary defects, such as second phases, voids and microcracks, which are randomly distributed at grain boundaries before the material is subjected to loads. These defects are important sources of damage and cracks, and thus, the intergranular fracturing is the main failure mechanism in this type of material. Therefore, it is necessary to develop a numerical approach that can consider the microstructure and the defects, simulate microcracking in the IMC layer and simulate the macro failure behavior of the solder joint. This type of approach is important

2

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

not only for understanding the underlying mechanisms that govern the micro-mechanical behavior of the IMC layer but also for understanding the failure mode of solder joints. Two important aspects of this approach are the geometric representation of the grain microstructure of the IMC and the numerical modeling of the grain interfacial behavior, such as sliding, separating and intergranular cracking. 1.1. Geometric representation of grain microstructure Random polygonal grains generated by Voronoi tessellation have been used extensively in computational material science. Ghosh and Liu [10] developed a Voronoi cell finite element model (VCFEM) to simulate arbitrary microstructures of a multi-phase material. Two-dimensional Voronoi tessellation was used to analyze the brittle compressive fracture behavior of polycrystalline ice [11]. Three-dimensional Voronoi tessellation was used to study the fracture behavior of ceramic materials [12] and was also applied to model the heterogeneity of polycrystalline structures, in which different clusters of Voronoi tessellations were assigned to different crystallographic orientations and mechanical properties to represent different material phases [13]. Voronoi tessellations were applied to the damage analysis of aluminum based on arbitrary distribution functions of grain sizes [14]. 1.2. Numerical modeling of grain interfacial behavior Although there are some numerical methods, such as the multiscale boundary element method [15], Monte-Carlo simulation [16] and the lattice or spring network model [17], that are being developed to simulate cracking in polycrystalline materials, they cannot explicitly capture the microcracking behavior of the materials on the grain level. Based on the finite element method, Zavattieri et al. [18] presented an approach to simulate intergranular cracking in polycrystalline brittle materials. In their work, the grain microstructure was retained explicitly, and the initiation and propagation of microcracks at grain boundaries were simulated by interface elements that were placed at the grain boundaries. The interface element is one type of finite element that is embedded between two volume elements to model the debonding behavior of interfaces. One type of interface element is the thin conventional continuum finite element with damage effect, and the damage is described in terms of the stress and the strain [19]. The second type is the discrete interface element or the linkage element, in which only the opposite nodes of the two volume elements neighboring the interface are linked by rod or spring type elements [20]. The third type is the zero thickness cohesive interface element, which is a continuum type element, and its constitutive behavior is described by the cohesive zone model. The cohesive zone model was originally proposed by Barenblatt [21] to model the fracture behavior of quasi-brittle materials. Then, Dugdale [22] applied it on the fracturing process zone at the crack tips. In the cohesive zone model, the cohesive law or the relationship between the traction and the separation of the interface is such that, as the interfacial separation increases, the traction across the interface initially increases until achieving a maximum value, and then it decreases and eventually goes to zero when a complete separation occurs. The cohesive zone model mimics the physical failure process of two atoms that are initially bonded, and the cohesive law resembles the relationship between the interatomic force and the interatomic separation. As a result, the parameters in the cohesive law depend only on the materials, rather than on the geometry and loading condition of the structures. This feature is very different from the fracture mechanics

used in macrocrack analysis, and the cohesive zone model is applicable to problems with various dimensional scales and with finite inelastic deformation. The cohesive laws used in the cohesive interface element can be categorized into three groups based on the shape of the traction– separation curves: the exponential potential-based cohesive law originally proposed by Xu and Needleman [23], the idealized trapezoidal traction–separation law developed by Tvergaard and Hutchinson [24], and the bilinear type proposed by Geubelle and Baylor [25]. The cohesive law can be uncoupled or coupled. In the uncoupled cohesive law, the normal traction is independent of the tangential sliding displacement, and the tangential traction is independent of the normal opening displacement. The uncoupled cohesive law is applied in pure mode conditions, such as pure opening (mode-I) or pure shearing (mode-II) [26]. For mixed modes, the coupled cohesive law is widely used. In the coupled cohesive law, both the normal and tangential tractions depend on the normal and tangential separating displacements, and they are coupled by using an effective separating displacement [27] or coupling parameters [23]. This paper presents a grain level model to study the intergranular cracking of the IMC layer in solder joints. The polycrystalline microstructure of the IMC layer is explicitly included in the model by Voronoi tessellations. The cohesive interface elements with the coupled cohesive law are embedded along grain boundaries to simulate microcrack initiation and propagation in the IMC layer. The use of the cohesive interface element is advantageous because the initiation, growth and coalescence of microcracks are natural outcomes of the mechanical response of the IMC layer under the applied loadings. The cohesive interface element is numerically implemented by incorporating a user-defined subroutine (UEL) into the finite element code ABAQUS [28]. Then, based on the numerical approach, the effects of the grain shape, the randomly distributed grain boundary defects, the thickness of the IMC layer and the morphology of the solder/IMC interface on the microcrack patterns and the overall response of solder joints are investigated. The paper is organized as follows. In Section 2, the experimentally observed microstructures and the failure modes of the IMC layer in solder joints are introduced briefly. In Section 3, the procedure to generate two-dimensional Voronoi tessellations of the IMC grains is presented, and then the cohesive law and the numerical implementation of the cohesive interface element used in this paper are presented in detail. In Section 4, the numerical results on the micro and macro behavior of the IMC layer in the solder joint are presented and discussed.

2. Microstructure and failure mode of the IMC layer 2.1. Microstructure A schematic cross section and a scanning electron microscopy (SEM) observation of a solder joint are shown in Fig. 1a and b, respectively. As depicted, Cu6Sn5 and Cu3Sn are the two major phases in the IMC layer of the solder joint. The Cu6Sn5 is formed first as the Cu atoms diffuse from the Cu pad into the molten Sn3.0Ag0.5Cu solder during the liquid-state soldering process, which is dominated by a chemical reaction [29]. After the soldering, the Cu6Sn5 continues to grow due to the effect of volume diffusion during the solid-state isothermal aging. When the Cu6Sn5 grows into a certain thickness, the Cu3Sn begins to form between the Cu pad and the Cu6Sn5. The Cu3Sn layer is a thin and uniform layer, and the Cu6Sn5 layer is thicker and has a scallop-like waved interface with the solder [30].

3

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

Solder (Sn3.0Ag0.5Cu)

Cu6Sn5

IMC layer Cu pad

(a) (a)

Solder (Sn3.0Ag0.5Cu)

Cu6Sn5

IMC layer

Cu3Sn Cu pad

(b)

(b) Fig. 1. (a) Schematic cross-section of a solder joint; and (b) SEM observation of the IMC layer.

2.2. Evolution of IMC microstructure in solder joints In the liquid-state soldering process, the solder/IMC interface shows a scallop-like morphology and is extremely rough [31], as shown in Fig. 2a. During the solid-state isothermal aging, the thickness of the IMC layer and the average grain size of the Cu6Sn5 increase with increasing aging time [32]. Because the solder in the valley of the IMC layer is closer to the Cu pad, the Sn atoms from the solder and the Cu atoms from the Cu pad reach each other faster in the valley, which leads to a higher growth velocity in the valley than that in other parts of the IMC layer. As a result, the roughness of the scallop-like interface is gradually reduced with increasing aging time, as shown in Fig. 2b. The grain microstructure of the Cu6Sn5 is hexagonal-shaped, as shown in Fig. 2c. 2.3. Failure modes The microstructure has a significant influence on the fracture toughness and the failure mode of the IMC layer. With the increasing IMC thickness, the failure mode migrates from ductile failure in the solder bulk to brittle cleavage within the IMC layer, which degrades the strength of the solder joints [3]. There are three failure modes observed in solder joints under tensile or shear loading conditions [3,4], as shown in Fig. 3.

Cu6Sn5

(c) Fig. 2. SEM observations of the IMC layer after isothermal aging at 150 °C for (a) 0 h, (b) 500 h and (c) 72 h, showing the Cu6Sn5 grains.

Mode 1. When the thickness of the IMC layer is approximately 1.0 lm, the failure is dominated by the fracturing of the solder close to the Cu6Sn5. Intact Cu6Sn5 grains can be observed in the fractography of this failure mode. Fracturing usually does not occur within the IMC layer when the IMC layer is very thin. Mode 2. When the thickness of the IMC layer is between 1 and 10 lm, the failure is dominated by the mixed fracturing in both the solder and the Cu6Sn5 grains. The fracture path is along the waist of the protruding portion of the IMC layer. Brittle intergranular cracking of the Cu6Sn5 grains can be observed in the fractography of this failure mode.

4

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

Solder Crack path

Mode 1 Mode 2 Mode 3 Cu6Sn5

Cu3Sn

Cu pad Fig. 3. Three failure modes of solder joints.

m = 18, n = 10 α =β = 0

Mode 3. When the thickness of the IMC layer is greater than 10 lm, the failure is dominated by the fracturing of the Cu6Sn5 grains, which occurs completely in the IMC layer. Brittle cleavage and broken Cu6Sn5 grains can be observed in the fractography of this failure mode.

(a)

Fracturing in the Cu3Sn layer is seldom observed; therefore, more attention is given to the Cu6Sn5 layer in the following discussion.

3. Grain level model of the IMC layer 3.1. Geometric model of IMC grains Voronoi tessellations were used to generate the geometry of the IMC grains. For a two-dimensional problem, let P denote a set of n points (called the nucleus) in a plane. The Voronoi diagram of P is defined as the subdivision of the plane into n cells, which is

 Vðpi Þ ¼ fq 2 Planedðq; pi Þ 6 dðq; pj Þ for each pi 2 P; j–ig

m = 18, n = 10 α = β = 0.5

ð1Þ

where pi represents a nucleus, q is a point lying in the cell, and d(q, pi) denotes the Euclidean distance between q and pi. Each cell represents an individual grain. A two-dimensional finite element model including randomly shaped grains was constructed based on Voronoi tessellations. The nuclei define the shape and arrangement of grains; therefore, their distribution should be predefined carefully. After isothermal aging, the crystalline structure of Cu6Sn5 is typically hexagonal with a relatively regular grain size distribution, as shown in Fig. 2c. Thus, the distance between the neighboring nuclei should not be excessively short or long. First, the plane was divided into m  n regular hexagons, where m and n were the number of regular hexagons in the x and y directions in a plane coordinate system, respectively. Then, the center points of each regular hexagon (xi, yi) (i = 1, 2, . . ., m  n), were obtained. The coordinates of the nuclei (ai, bi) were randomly selected within a range around the center points, which are (xi  axi, xi + axi) and (yi  byi, yi + byi). The coefficients a and b, which are in the range of (0, 1), were used to control the location of the nuclei within the hexagon. When a = b = 0, the nucleus was located in the hexagonal center, and the collection of regular hexagons was generated (see Fig. 4a). With increasing a and b, the irregularity of the Voronoi tessellations increases (see Fig. 4b). The coordinates of the Voronoi tessellations generated by MATLAB were saved as a data file; then, the vertex coordinates were read into a file using the Python scripting language, and the vertices were connected sequentially. Finally, the geometric model of the polycrystalline microstructures was constructed by ABAQUS/CAE and the Python interface.

(b) Fig. 4. Voronoi tessellations of polycrystalline materials with (a) regular hexagonally shaped grains and (b) irregular polygonally shaped grains.

3.2. Constitutive behavior of the cohesive interface element 3.2.1. Cohesive law The behavior of the cohesive interface element is described by a constitutive relationship, or a cohesive law, which is specified by the relationship between the separating displacement d and the corresponding traction T on the interface. One of the most popular cohesive laws is the exponential cohesive law [23], which is expressed as

!  ( /n dn dn d2t Tn ¼ exp  exp  2 dcr;n dcr;n dcr;n dcr;t " !# ) 1q d2t dn r 1  exp  2 þ r1 dcr;n dcr;t Tt ¼ 2

ð2Þ

!     r  q d  /n dt dn d2 n qþ exp  exp  2t r  1 dcr;n dcr;t dcr;t dcr;n dcr;t ð3Þ

where Tn and Tt are the normal and tangential traction components, respectively; dn and dt are the normal and tangential separating displacements, respectively; dcr,n and dcr,t are the critical separating

5

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

displacements in the normal and tangential directions, respectively; and q and r are two coupling parameters, defined as q = /t//n and r = dn /dcr,n, where dn is the value of dn after complete shear separation with Tn = 0. Here, /n and /t are the energy release rates under the pure mode-I and pure mode-II, respectively, and the values of /n and /t are equal to the areas below the traction–separation curves. The exponential cohesive law above describes the coupling between the opening and shearing modes only in the case of /n = /t, which limits its application in mixed modes. In this paper, a modified exponential cohesive law proposed by van den Bosch et al. [33] was used, as depicted in Fig. 5. The modified cohesive law preserves all features of the original law but has a better, controllable coupling effect. In the modified cohesive law, the two energy release rates corresponding to the opening and shearing modes are treated as independent parameters. Thus, Tn and Tt are expressed as

!     / dn dn d2 exp  exp  2t ; Tn ¼ n dcr;n dcr;n dcr;n dcr;t

dn P 0; Ddn P 0 ð4Þ

!      /t dt dn dn d2t 1þ exp  exp  2 ; Tt ¼ 2 dcr;t dcr;t dcr;n dcr;n dcr;t

Dd t P 0 ð5Þ

1.2 0.8

Elastic unload

where Ddn and Ddt are the increments of the separating displacement. According to the cohesive law, for normal separation, the normal traction increases with the increasing of the opening displacement until the maximum value, Tn,max, is reached, and then the normal traction gradually decreases to zero. Under the pure shearing mode, the tangential traction increases with the increase of the tangential displacement until the maximum value, Tt,max, is reached, and then the tangential traction gradually decreases to zero. Thus, the maximum values of normal and tangential tractions can be derived as

T n; max ¼

/n ; dcr;n expð1Þ

T t; max ¼

/t qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 dcr;t 2 expð1Þ

ð6Þ

By observing Eqs. (4)–(6), we found that /n, /t, dcr,n (or Tn,max) and dcr,t (or Tt,max) are four independent parameters in the cohesive law. 3.2.2. Unloading and reloading Unloading and reloading occur when external loadings change or the internal stresses in the neighboring elements are redistributed. In the cohesive interface element, before the maximum traction is reached, the traction for unloading and reloading obeys the same law as the loading process. After the maximum traction is reached, the traction for unloading and reloading is along an elastic path that connects the current point and the origin of the traction– separation curve, as shown in Fig. 5, i.e.,

DT n ¼ k n Dd n ; Dd n < 0

ð7Þ

DT t ¼ kt Ddt ; Ddt < 0

ð8Þ

Tn / Tn,max

0.4

where DTn and DTt are the increments of the normal and tangential tractions and kn and kt are the unloading/reloading stiffness in the normal and tangential directions, respectively.

Reload

0.0 -0.4

Compression -0.8 -1.2 -1

0

1

2

4

3

5

6

δ n / δ n,cr

(a)

1.2

 2 dn T c ¼ a ; dtol

0.8

Elastic unload Tt / Tt,max

Reload Reload

-0.4 -0.8

-1.6

-0.8

ð9Þ

3.3. Finite element formulas of the cohesive interface element

Elastic unload

-1.2 -2.4

dn < 0

in which a is the contact stiffness and dtol is a specific penetration tolerance. In this way, the resistance of the interface element against the compression is increased, and the penetration can be prevented.

0.4 0.0

3.2.3. Compression When the cohesive interface element is in compression, it is necessary to prevent the two volume elements connected with the interface element from penetrating into each other. There are several ways to address this issue, including implementing contact algorithms at the interfaces [34], carefully choosing the value of normal stiffness [14], and modifying the constitutive equation in compression to guarantee that the penetration is tolerable [35]. In this paper, an extra compressive traction Tc was introduced if the normal opening displacement became negative,

0.0

0.8

1.6

2.4

δ t / δ t,cr

(b) Fig. 5. Cohesive laws: (a) the normal traction–separation relationship defined by Eq. (4) with the maximum traction Tn,max at dn/dn,cr = 1 and dt  0; and (b) the tangential traction–separation pffiffiffi relationship defined by Eq. (5) with the maximum traction Tt,max at dt/dt,cr = 1= 2 and dn  0.

Fig. 6 shows a two-dimensional quadrilateral cohesive interface element with four nodes, two integration points, and a local element coordinate n (1 6 n 6 1) defined along the mid-line of the cohesive interface element [36]. The initial thickness of the cohesive interface element before deformation is zero, as shown in Fig. 6a. As the adjacent volume elements deform, the two surfaces of the cohesive interface element separate, and the mid-line of the interface element deforms simultaneously. To describe the traction and separation in the normal and tangential directions, another local coordinate system (n, t) connected with the deformed mid-line was defined, as shown in Fig. 6b.

6

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

xRI ¼ Continuum element 4 Mid-plane 1

ξ

2

Continuum element

0.0 Cohesive interface element

xR ðnÞ ¼

4 Mid-plane 1 y

Continuum element

×

xR ðnÞ yR ðnÞ

¼ HðnÞxRI

ð18Þ

0 @N 2 ðnÞ @n

ð19Þ

# xRI

The 2  2 transformation matrix H is defined by the direction cosines of the local frame to the global frame

3

 H¼

2 Continuum element

cos h  sin h sin h cos h

 ð20Þ

with @yR

@xR

(b)

@n cos h ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     ; 2

@xR @n

0



0 R1 @x @xR ðnÞ @ @n A @HðnÞ R ¼ ¼ x @yR @n @n I @n " @N1 ðnÞ 2 ðnÞ 0 @N@n @n ¼ @N 1 ðnÞ 0 0 @n

×

t



and its derivatives with respect to n can be written as

×: Gaussian integration points

(a)

n

ð17Þ

Therefore, analogous to Eq. (13), the coordinates of any point on the mid-line can be derived by 3

×

×

1 ½I44 jI44 ðXI þ uI Þ 2

þ

@yR @n

2

@n sin h ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     @xR @n

2

þ

@yR @n

2

ð21Þ

x

Fig. 6. The cohesive interface element in (a) the initial configuration and (b) the deformed configuration.

Thus, the local separating displacements can be obtained by the global displacements and the transformation matrix H as

Duloc ¼ HT Du

ð22Þ

The four-node cohesive interface element has 8 (2  4) degrees of freedom, and its nodal displacement vector (8  1) in the global coordinate system is written as

Then, the 2  1 local traction vector Tloc is determined by the constitutive law for the cohesive interface element, which was discussed in Section 3.2

uI ¼ ðu1 ; m1 ; u2 ; m2 ; u4 ; m4 ; u3 ; m3 ÞT

Tloc ¼

ð10Þ



utop ¼ ðu4 ; m4 ; u3 ; m3 ÞT I

ð11Þ

The separating displacements of the cohesive interface element are the relative displacements between each pair of nodes:

DuI ¼ utop  ubottom I I



DuðnÞ ¼

Dux ðnÞ Duy ðnÞ

 ¼ HðnÞDuI

ð13Þ

¼ HðnÞutop  HðnÞubottom ¼ BðnÞuI I I

ð23Þ

where

T loc;n ¼

!     Du2loc;t /n Duloc;n Duloc;n exp  exp  2 dcr;n dcr;n dcr;n dcr;t

ð12Þ

The continuous separating displacement field within the element is then expressed as



T loc;t

The nodal displacements of the bottom and top nodes are

ubottom ¼ ðu1 ; m1 ; u2 ; m2 ÞT ; I

T loc;n

T loc;t ¼ 2

ð24Þ

!      Du2loc;t /t Duloc;t Duloc;n Duloc;n 1þ exp  exp  2 dcr;t dcr;t dcr;n dcr;n dcr;t ð25Þ

The global tractions can be calculated by the local tractions and the transformation matrix H as

T ¼ HTloc

ð26Þ

with

BðnÞ ¼ ½HðnÞjHðnÞ

ð14Þ

3.4. Numerical implementation of the cohesive interface element

ð15Þ

The cohesive interface element was incorporated into the finite element code ABAQUS through a user-defined subroutine UEL. In the UEL subroutine, the nodal force vector and the element stiffness were computed and output to the main routine of ABAQUS. The nodal force vector was computed by

ð16Þ

Fel ¼

where H(n) is a 2  4 matrix of the shape function given by

 HðnÞ ¼

N1 ðnÞ

0

N2 ðnÞ

0

0

N1 ðnÞ

0

N2 ðnÞ



and the shape functions are

N1 ðnÞ ¼

1 1 ð1  nÞ N2 ðnÞ ¼ ð1 þ nÞ 2 2

The separating displacements in the global coordinate system can be transformed into the local coordinate system (n, t) by a transformation matrix H. Let XI be the initial coordinate of nodes on the mid-line; after deformation, its new coordinate xRI can be calculated by the displacement of the top and bottom nodes in the interface element as

Z Ael

BT TdA ¼ W

Z

1

BT HTloc detðJÞdn

ð27Þ

1

where Ael and W are the area and the width of the cohesive interface element, J is the Jacobian matrix, and

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  R 2  R 2 @x @y detðJÞ ¼ þ @n @n

ð28Þ

7

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

The tangential stiffness of the interface element was computed by el

Kel ¼

@F ¼W @uI

Z

1

BT HDloc HT B detðJÞdn

ð29Þ

1

where Dloc is the cohesive constitutive matrix derived by

Dloc ¼

 Dnn @Tloc ¼ @ Duloc Dtn

Dnt

 ð30Þ

Dtt

where Dnn, Dnt Dtn and Dtt are the derivatives of the tractions with respect to the separations. Two numerical integration schemes, Newton–Cotes and Gaussian integration, can be used in Eqs. (27) and (29). The difference between Newton–Cotes and Gaussian integration is the location of the integration points. In the four-node interface element, the location of the integration point for the Newton–Cotes scheme is at the

Displacement U2

mid-point of the line connecting the paired bottom and top nodes, and thus, only the relative displacement of one nodal pair is involved in the integrating calculation, while the relative displacement of the other nodes in the element has no influence on the integration. In contrast, for the Gaussian integration scheme, the relative displacements of all the nodes in the element are involved in the integration at each integration point. Some studies have suggested that the Newton–Cotes integration scheme should be used for the interface elements because the application of Gaussian integration leads to oscillatory traction profiles in the interface elements [37]. However, a thorough investigation for the impact of the Gaussian integration showed that the oscillating stress was caused by using a large element size in the regions with a great stress gradient, not by the integration scheme, and thus, the stress oscillation can be overcome by carefully setting the size of the elements adjacent to the interface [38]. Because the Newton–Cotes integration scheme has no benefit over Gaussian integration, the Gaussian integration scheme was used in this paper.

4. Result and discussion

54 μm

tm = 9.5 μm

Solder

Cu6Sn5

24 μm

Cu pad Interface elements 46 μm

(a) Displacement U2

4.1. Effect of grain shape To study the effect of grain shapes, two models with regular hexagonal and randomly irregular polygonal grain shapes generated by Voronoi tessellations were constructed, as shown in Fig. 7a and b, respectively. The two models had a size of 54  46 lm2, and the thickness of the Cu pad was 24 lm. The thickness of the IMC (Cu6Sn5) layer was measured as an average thickness tm, which was defined by dividing the area A of the IMC layer by its width L, i.e., tm = A/L. The average thicknesses of the Cu6Sn5 layers in the regular and irregular grain shape models were 9.5 and 9.3 lm, respectively. The two models contained the same number of grains and grain boundaries. There were 48 grains and 138 grain boundaries in each model, and the average length of the grain facets was approximately 2 lm. The four-node bilinear plane strain elements, type CPE4 in ABAQUS, were used to mesh the Cu pad, the solder part and the interior of the Cu6Sn5 grain. The four-node cohesive interface elements described in Section 3 were embedded along all of the grain boundaries to simulate the initiation, propagation and coalescence of microcracks. There were 15,791 elements and 15,948 nodes in total in the regular grain shape model and 13,850 elements and 14,001 nodes in the irregular grain shape model. The two models were subjected to the same vertical displacement loading on the

Solder

0.12

A

Cu6Sn5

24 μm

Cu pad Interface elements

A' Reaction force / N

54 μm

tm = 9.3 μm

0.10

0.08

B' B

0.06

C' C

0.04

Regular grain shape model

0.02

46 μm

Irregular grain shape model 0.00 0.0

(b) Fig. 7. Models with (a) a regular grain shape and (b) an irregular grain shape.

0.2

0.4

0.6

0.8

1.0

Displacement / μm Fig. 8. Effect of the grain shape on the overall force–displacement response.

8

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

top side of the model. The boundary conditions of the model are shown in Fig. 7. The Young’s modulus E of the Cu6Sn5 was 119 GPa, the Poisson’s ratio m was 0.31, and the stress intensity factor KI was 2.73 ± 0.63 MPa m1/2 [39]. Thus, the critical energy release rate Gc of the Cu6Sn5 can be calculated by

Gc ¼

ð1  m2 ÞK 2I E

ð31Þ

The normal and tangential critical energy release rates used in the simulation were /n = /t = Gc = 33.52 N/m, and the maximum tractions in the normal and tangential directions were Tn,max = 68 MPa [40] and Tt,max = 28 MPa [41], respectively. The Young’s modulus of the Cu pad was 117 GPa, and the Poisson’s ratio was 0.34; the Young’s modulus of the Sn3.0Ag0.5Cu was 54 GPa, and the Poisson’s ratio was 0.36. The Cu pad, the solder and the Cu6Sn5 grain in the models were treated as isotropic elastic materials.

The overall force were obtained by averaging the reaction forces at the nodes on the top side of the model, and the force–displacement response curves are shown in Fig. 8. The overall responses of the two models are almost the same, which suggests that the overall response is not sensitive to the grain shape. To examine the difference in the microcrack patterns in the two models, three load steps labeled as A, B and C for the regular grain shape model in Fig. 8 were selected, and their corresponding microcrack patterns are shown in Fig. 9. Similarly, the microcrack patterns corresponding to load steps A0 , B0 and C0 for the irregular grain shape model in Fig. 8 are shown in Fig. 10. Fig. 9 shows the microcrack patterns and the vertical normal stress ryy in the regular grain shape model. At load step A, microcracks initiate at the top left and bottom right of the Cu6Sn5 layer, as shown in Fig. 9a. In the central region of the Cu6Sn5 layer, microcracks initiate almost simultaneously because the grain boundaries perpendicular to the loading direction are subjected to almost the

Solder A starting position of cracking

A starting position of cracking Cu pad

(a)

Solder

These grain boundaries debond at the same time.

Solder

Cu pad

Cu pad

(b)

(c)

Fig. 9. Microcrack patterns and stress ryy of the regular grain shape model at (a) load step A, (b) load step B and (c) load step C (unit: MPa).

9

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

Solder Microcracks are stopped and deviated by grain boundaries.

Triple-grain junction Cu pad

(a) Solder

Solder

Coalescence of microcracks

New microcracks form.

Cu pad

Cu pad

(b)

(c)

Fig. 10. Microcrack patterns and stress ryy of the irregular grain shape model at (a) load step A0 , (b) load step B0 and (c) load step C0 (unit: MPa).

0.12 30

Weibull parameters: κ = 5, λ = 0.59 0.10

Reaction force / N

Number of interfaces

25

20

15

10

0.08

0.06

0.04

0.02 5

0 0.4

Uniform grain interfacial strength model Weibull grain interfacial strength model

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

i Tmax Tmax

Fig. 11. Weibull distribution of the grain interfacial strength.

1.6

0.00 0.0

0.2

0.4

0.6

0.8

1.0

Displacement / μm Fig. 12. Effect of Weibull distributed grain interfacial strength on the overall force– displacement response.

10

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

same stress. At load step B, microcracks propagate, and more grain interfaces open; they coalesce to form a main crack across the grain layer, as shown in Fig. 9b. The path of the main crack has a zigzag shape, but it tends to remain in the horizontal direction, which is perpendicular to the loading direction. Fig. 9c shows that the stress contour is substantially symmetrical because the orientation and starting time of the microcracks are almost the same in the regular grain shape model. Fig. 10 shows the microcrack patterns and the vertical normal stress ryy in the irregular grain shape model. At load step A0 , microcracks initiate at locations similar to those in the regular grain shape model, as shown in Fig. 10a, but microcrack deviation is observed at the triple-grain junctions. Some cracks are stopped when they meet a grain boundary, which gives the other grain interfaces in the model an opportunity to develop new microcracks, as shown in Fig. 10b. As a result, the crack path in the irregular grain shape model is more zigzag, and the stress contour is asymmetrical, as shown in Fig. 10c.

4.2. Effect of grain interfacial strength Polycrystalline materials usually contain various defects, such as voids, pre-existing cracks and secondary phases at grain boundaries. These defects should be taken into account when predicting the micro and macro-mechanical behavior of the materials [42]. In this paper, a Weibull distribution [43] of grain interfacial strengths was adopted to model the effect of these randomly distributed defects. The four parameters in the cohesive law, both of the energy release rates /n and /t, and the maximum tractions Tn,max and Tt,max were assigned values randomly according to the following probability density functions:

f ð/in Þ

¼

f ð/it Þ ¼

!j # /in ; exp  /n /n !j1 " !j # /it /it exp  /t /t

j /in

kj

j kj

!j1

"

Solder

Cu pad

(a)

Solder

Cu pad

(b) Fig. 13. Microcrack patterns of models with (a) a uniform grain interfacial strength and (b) a Weibull distributed grain interfacial strength (unit: mm).

ð32Þ

11

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

Displacement U2

0.105 0.090

Solder

0.075

Cu6Sn5

Reaction force / N

12 µm

56 µm

λ =12 µm Ra=8 µm

24 µm

Cu pad L = 112 µm

0.060 0.045 0.030

Model I: tm=8 μm, λ=12 μm, Ra=8 μm Model II: tm=9.5 μm, λ =24 μm, Ra=8 μm

0.015

(a)

Model III: tm=13 μm,λ =24 μm, Ra=4 μm

Displacement U2

0.000 0.00

0.15

0.30

0.45

0.60

0.75

Displacement / μm Solder

Cu6Sn5

Fig. 15. Effect of IMC microstructure on the overall force–displacement response.

13.5 µm

58 µm

24µm λ ==24

Ra = 8 µm

parameters were j = 5 and k = 0.59, and then all four interfacial strength parameters changed in a range from 50% to 150% of the average values used in the uniform grain interfacial strength model. There were 138 grain interfaces in this model, and the distribution of interfaces with different strengths is shown in Fig. 11. The overall force–displacement responses of the Weibull distributed grain interfacial strength model are compared with that of the uniform grain interfacial strength model in Fig. 12. The overall strength of the Weibull distributed grain interfacial strength model is approximately 15% lower than that of the uniform grain interfacial strength model, even though half of the total interfaces are stronger in the Weibull distributed grain interfacial strength model than in the uniform grain interfacial strength model. This result implies that the weaker grain interfaces play a critical role in the overall strength of the model. Microcrack patterns at the final load step of the two models are shown in Fig. 13. Although the same grain shape and configuration were used in the two models, their microcrack patterns were different. In the uniform grain interfacial strength model, the crack path is mainly governed by the shape and the orientation of grains, while in the Weibull distributed grain interfacial strength model, the crack path is controlled not only by the shape and the orientation of grains but also by the weak grain interfaces.

24 µm

Cu pad L = 112 µm

(b) Displacement U2

Solder Ra = 4 µm

Cu6Sn5

15 µm

61.5 µm

λ = 24 µm

24 µm

Cu pad L = 112 µm

(c)

4.3. Effect of the IMC microstructure

Fig. 14. Models with different IMC microstructures: (a) Model I: tm = 8 lm, k = 12 lm, Ra = 8 lm; (b) Model II: tm = 9.5 lm, k = 24 lm, Ra = 8 lm; and (c) Model III: tm = 13 lm, k = 24 lm, Ra = 4 lm.

!j # T in; max ; T n; max T n; max !j1 " !j # T it; max T it; max exp  T t; max T t; max

j T in; max

f ðT in; max Þ ¼

kj

f ðT it; max Þ ¼

kj

j

!j1

"

exp 

ð33Þ

where j and k are the Weibull parameters. For comparison, a model with uniform grain interfacial strength was built, and its model parameters were /n = /t = 33.52 N/m, Tn,max = 68 MPa, and Tt,max = 28 MPa, which were the statistical averages of the Weibull distribution. In the model with the Weibull distributed grain interfacial strength, every grain interface was randomly assigned a set of the four parameters according to Eqs. (32) and (33). The Weibull

In solder joints, the IMC microstructure can be characterized by the thickness of the IMC layer and the morphology of the solder/ IMC interface. In this paper, three geometric parameters, the average thickness tm of the IMC layer, the wavelength k and the roughness Ra of the solder/IMC interface, were used to describe the IMC microstructure. The average thickness tm is defined in the same way as in Section 4.1. The wavelength k is defined as the average spacing between two adjacent valleys of the waved solder/IMC interface. The roughness Ra is equal to the distance between the peak and the valley of the waved solder/IMC interface. To investigate the effect of the IMC microstructure on the failure behavior of the solder joints, three models with different IMC microstructures, Model I, Model II and Model III, were constructed, as shown in Fig. 14a–c, respectively. The microstructural parameters of the IMC (Cu6Sn5) layer in Model I were tm = 8.0 lm, k = 12.0 lm and Ra = 8.0 lm, which characterized a thin, rough and scallop-like IMC layer just after the soldering process. For Model II, the parameters were tm = 9.5 lm, k = 24.0 lm and Ra = 8.0 lm. For Model III, the parameters were tm = 13.0 lm,

12

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

k = 24.0 lm and Ra = 4.0 lm. Model III characterized a thick and relatively flat IMC layer after isothermal aging. There were 151 grains and 394 grain interfaces in Model I, 148 grains and 397 grain interfaces in Model II, and 184 grains and 502 grain interfaces in Model III. In all three models, the interfacial strength parameters were set to be the same as in the Weibull distributed interfacial strength model in Section 4.2. Fig. 15 shows the overall responses of the three models. Model I has the greatest overall strength, while Model III has the lowest one. The overall strength of Model II is slightly lower than that of Model I. Compared to Model I, Model II has a slightly thicker IMC layer, twice the wavelength and the same roughness, which seems to imply that the wavelength has little impact on the overall strength, but the thickness of the IMC layer plays a key role. This conclusion can be confirmed by further examining the response of Model III. Model III has the thickest IMC layer, the same wavelength as Model II and a smaller roughness, but it has the lowest overall strength. One explanation is that the model with the thicker IMC layer contains more grain interfaces, and these interfaces are vulnerable sites to initiate cracks. The results also indicate that

the roughness of the solder/IMC interface has little impact on the overall strength of the model. Fig. 16 shows the final failure modes of the three models. In Model I, all of the microcracks start at the valley bottom of the waved solder/IMC interface because severe stress concentration occurs there, and then they propagate across the root of the protruding Cu6Sn5 grains, as shown in Fig. 16a. The failure mode of Model II is presented in Fig. 16b. In Model II, some microcracks start at the valley bottom of the waved solder/ IMC interface, which happens in Model I, but more microcracks start within the IMC layer. It is interesting that the microcracks starting in the IMC layer propagate toward the waist of the protruding Cu6Sn5 grains. This result coincides with the experimental observation [4] and can explain why broken Cu6Sn5 grains are observed in the fractography. The failure mode in Model II is similar to the observed failure Mode 2 depicted in Fig. 3. The difference between the predicted failure mode and the experimental observation is that the numerically simulated crack does not propagate into the solder because the failure in the solder bulk was not considered in this model.

Solder

Solder

A starting position of cracking is at the solder/IMC layer interface.

A starting position of cracking is at the waist of the protrudent Cu6Sn5 grains.

A starting position of cracking is at the solder/IMC interface.

Microracks propagate across the root of the protruding Cu6Sn5 grains.

Microcracks propagate toward the waist of the protruding Cu6Sn5 grains.

A starting position of cracking is within the IMC layer.

Cu pad

Cu pad

(a)

(b)

Solder

Microcracks start and propagate completely within the IMC layer. Cu pad

(c) Fig. 16. Effect of IMC microstructure on failure mode: (a) Model I: tm = 8 lm, k = 12 lm, Ra = 8 lm; (b) Model II: tm = 9.5 lm, k = 24 lm, Ra = 8 lm; and (c) Model III: tm = 13 lm, k = 24 lm, Ra = 4 lm (unit: mm).

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14

The failure mode of Model III is presented in Fig. 16c. In Model III, almost all of the microcracks start and propagate within the IMC layer. The failure mode of Model III is a typical failure Mode 3, as depicted in Fig. 3. The results indicate that, from Model I to Model III, the overall strength decreases, and the failure mode transforms from Mode 2 to Mode 3. These results can be interpreted as follows. As discussed in Section 4.2, the weak grain interfaces play a critical role in the overall strength. Model III has the thickest IMC layer and thus contains the most weak grain interfaces. As a result, the overall strength of Model III is the lowest. Furthermore, a thicker IMC layer contains more grain interfaces and thus increases the opportunity of cracking in the IMC layer, which is one reason why the failure mode transforms from Mode 2 in Model II to Mode 3 in Model III. The other reason is that Model III has the smallest relative roughness among the three models if we define a relative roughness Rr as Rr = Ra/k to measure the morphology of the solder/IMC interface. In other words, the solder/IMC interface in Model III is the flattest one, which greatly reduces the stress concentration at the interface and thus reduces the opportunity for cracking to start at the solder/IMC interface.

5. Conclusions This paper presents a grain level numerical model to simulate the cracking behavior of the IMC layer in solder joints. The grain microstructure of the IMC layer was explicitly included in the model by Voronoi tessellations. The cohesive interface elements with the coupled cohesive law were embedded along grain boundaries to model the initiation, growth and coalescence of microcracks in the IMC layer. The cohesive interface element was numerically implemented by incorporating a user-defined subroutine UEL into the finite element code ABAQUS. Using the numerical approach developed, the effects of the grain shape, the randomly distributed grain boundary defects, the thickness of the IMC layer and the morphology of the solder/IMC interface on the microcrack patterns and on the overall response of solder joints were investigated. (1) Two models with regular hexagonal and randomly irregular polygonal grain shapes were constructed to investigate the effect of the grain shape. The results show that the overall mechanical strength is not sensitive to the grain shape, but the microcrack pattern and the crack path heavily depend on the grain shape. (2) A model with a Weibull distributed grain interfacial strength was adopted to account for the randomly distributed grain boundary defects. The numerical results show that the overall strength of the model with the Weibull distributed grain interfacial strength is significantly lower than that of the model with a uniform grain interfacial strength. The microcrack pattern and the crack path are different from those in the model with a uniform grain interfacial strength. The weak grain interface plays a critical role in the overall strength and the crack path of the model. (3) The average thickness of the IMC layer and the wavelength and roughness of the solder/IMC interface were used to characterize the IMC microstructure. Three models with different IMC microstructures were constructed to investigate the effect of the IMC microstructure on the failure behavior of solder joints. The results indicate that the average thickness of the IMC layer has the greatest influence on the overall strength and the failure mode of the solder joint. The thicker the IMC layer is, the lower the overall strength and the greater the number of microcracks starting within the IMC layer. The wavelength and the roughness of the

13

solder/IMC interface have little influence on the overall strength but affect the failure mode of the solder joint. When the solder/IMC interface has a rough and scallop-like morphology, the microcrack prefers to start at the valley bottom of the waved solder/IMC interface and to propagate across the root of the protruding Cu6Sn5 grains. When the solder/ IMC interface is relatively smooth or flat, the microcrack prefers to start and propagate within the IMC layers. The predicted failure mode agrees well with the experimental observation.

Acknowledgment The authors would like to acknowledge financial support from the National Natural Science Foundation of China (NSFC) under Grant No. 10972012. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.commatsci.2013. 05.044. References [1] H.T. Lee, M.H. Chen, Mater. Sci. Eng. A – Struct. Mater. Prop. Microstruct. Process. 333 (2002) 24–34. [2] T.Y. Lee, W.J. Choi, K.N. Tu, J.W. Jang, S.M. Kuo, J.K. Lin, D.R. Frear, K. Zeng, J.K. Kivilahti, J. Mater. Res. 17 (2) (2002) 291–301. [3] S.M. Hayes, N. Chawla, D.R. Frear, Microelectron. Reliab. 49 (2009) 269–287. [4] H.T. Lee, M.H. Chen, H.M. Jao, T.L. Liao, Mater. Sci. Eng. A – Struct. Mater. Prop. Microstruct. Process. 358 (2003) 134–141. [5] John H.L. Pang, F.X. Che, Drop impact analysis of Sn–Ag–Cu solder joints using dynamic high strain rate plastic strain as impact damage driving force, in: Proceedings of 56th Electronic Components and Technology Conference, San Diego, CA, United States, 2006. Piscataway, NJ, United States, IEEE, 2006, pp. 49–54. [6] M.O. Alam, H. Lu, Chris Bailey, Y.C. Chan, Comput. Mater. Sci. 45 (2009) 576– 583. [7] T. An, F. Qin, Trans. ASME J. Electron. Packag. 133 (2011) 031004-1–031004-7. [8] W. Ronnie Teo, Y.F. Sun, Acta Mater. 56 (2008) 242–249. [9] W.H. Müller, T. Hannach, H.J. Albrecht, FE-investigation of the stress/strain and fracture mechanics properties of intermetallic phase regions in leadfree solder interconnects, in: Proceedings of 8th Electronic Packaging Technology Conference, Singapore, 2006. Piscataway, NJ, United States, IEEE, 2006, pp. 114–120. [10] S. Ghosh, Y. Liu, Int. J. Numer. Methods Eng. 38 (1995) 1361–1398. [11] M.S. Wu, Jian Niu, Mech. Mater. 20 (1995) 9–32. [12] Y. Toi, T. Kiyosue, Eng. Fract. Mech. 50 (1995) 11–27. [13] F. Barbe, L. Decker, D. Jeulin, G. Cailletaud, Int. J. Plasticity 17 (2001) 513–536. [14] T. Luther, C. Könke, Eng. Fract. Mech. 76 (2009) 2332–2343. [15] G.K. Sfantos, M.H. Aliabadi, Comput. Method Appl. Mech. Eng. 196 (2007) 1310–1329. [16] L. Chen, R. Ballarini, M. Grigoriu, Int. J. Fract. 125 (2004) 353–369. [17] Martin Ostoja-Starzewski, Appl. Mech. Rev. 55 (2002) 35–59. [18] P.D. Zavattieri, P.V. Raghuram, H.D. Espinosa, J. Mech. Phys. Solids 49 (2001) 27–68. [19] A. Musienko, G. Cailletaud, Acta Mater. 57 (2009) 3840–3855. [20] D. Xie, A.M. Waas, Eng. Fract. Mech. 73 (2006) 1783–1796. [21] G.I. Barenblatt, Adv. Appl. Mech. 7 (1962) 56–129. [22] D.S. Dugdale, J. Mech. Phys. Solids 8 (1960) 100–104. [23] X.-P. Xu, A. Needleman, Model. Simul. Mater. Sci. Eng. 1 (1993) 111–132. [24] V. Tvergaard, J.W. Hutchinson, J. Mech. Phys. Solids 40 (1992) 1377–1397. [25] P.H. Geubelle, J.S. Baylor, Compos. Part B-Eng. 29 (1998) 589–602. [26] S. Li, M.D. Thouless, A.M. Waas, J.A. Schroeder, P.D. Zavattieri, Compos. Sci. Technol. 65 (2005) 281–293. [27] P.P. Camanho, C.G. Dávila, NASA/TM-2002-211737 (2002) 1–37. [28] ABAQUS Version 6.10 Documentation, ABAQUS Version 6.10 Documentation, Dassault Systèmes Simulia Corp., Providence, RI, USA, 2010. [29] J. Görlich, G. Schmitz, K.N. Tu, Appl. Phys. Lett. 86 (2005) 1–3. [30] K.N. Tu, T.Y. Lee, J.W. Jang, L. Li, D.R. Frear, K. Zeng, J.K. Kivilahti, J. Appl. Phys. 89 (2001) 4843–4849. [31] Jicheng Gong, Changqing Liu, Paul P. Conway, Vadim V. Silberschmidt, Acta Mater. 56 (2008) 4291–4297. [32] Fubin Song, Experimental investigation on testing conditions of solder ball shear and pull tests and the correlation with board level mechanical drop test,

14

[33] [34] [35] [36]

[37]

T. An, F. Qin / Computational Materials Science 79 (2013) 1–14 Ph.D. Thesis, The Hong Kong University of Science and Technology, Hong Kong, China, 2007. M.J. van den Bosch, P.J.G. Schreurs, M.G.D. Geers, Eng. Fract. Mech. 73 (2006) 1220–1234. H.D. Espinosa, P.D. Zavattieri, G.L. Emore, Mech. Mater. 29 (1998) 275–305. M.J. van den Bosch, P.J.G. Schreurs, M.G.D. Geers, Eur. J. Mech. A-Solids 26 (2007) 1–19. S. Feih, Development of a user element in ABAQUS for modelling of cohesive laws in composite structures, Technical Report Risø-R-1463 (EN), Roskilde, Denmark, Risø National Laboratory, 2004. J.C.J. Schellekens, R. De Borst, Int. J. Numer. Methods Eng. 36 (1993) 43–66.

[38] R.A. Day, D.M. Potts, Int. J. Numer. Anal. Methods Geomech. 18 (1994) 689– 708. [39] G. Ghosh, J. Mater. Res. 19 (2004) 1439–1454. [40] K.H. Prakash, T. Sritharan, Mater. Sci. Eng. A – Struct. Mater. Prop. Microstruct. Process. 379 (2004) 277–285. [41] Jie Zhao, Cong-qian Cheng, Lin Qi, Cheng-yu Chi, J. Alloys Compd. 473 (2009) 382–388. [42] P.D. Zavattieri, Computational modeling for bridging size scales in the failure of solids, Ph.D. Thesis, Purdue University, West Lafayette, IN, USA, 2000. [43] W. Weibull, S. Sweden, J. Appl. Mech. 18 (1951) 293–297.