I. INTRODUCTION

The mother board assembly constitutes morphologically the brain of an engineering device. The integrated circuits that are photolithographically etched onto silicon chips constitutes morphologically the brain cells. However, the brain cells on individual silicon chip are isolated from other chips; and they need to be electrically interconnected. The brain cells are very delicate and even the interconnects are mechanically fragile; and they need to be mechanically protected. Lastly, the brain cells generate high intensity of heat when running, and the heat needs to be dissipated before the brain cells would be burned. The jobs of electrical interconnection, mechanical protection, and thermal dissipation rest on the design and engineering of electronic and microelectronic assemblies [1].

Figure 1 shows the schematic of a silicon chip that is electrically interconnected to a substrate, which is in turn interconnected to a printed circuit board, through arrays of solder joints. The former assembly is frequently referred to as an integrated circuit component/assembly while the latter is frequently referred to as a printed circuit board assembly. We shall refer to both assemblies as simply electronic assemblies. The solder joints are formed by melting solder balls/paste to form metallurgical bonds with the metal pads on the chip and the substrate, and with the metal pads on the substrate and the printed circuit board.

Electronic assemblies are susceptible to undue stress during manufacturing and while in service resulting in functional failure. In general, electronic assemblies may experience three main physics of damage: the violent vaporisation of the ingress moisture in the integrated circuit assembly during the solder joint forming process leading to cracking and/or delamination of the assembly [2][3]; fracturing of interconnection caused by drop-shock of mobile electronic products [2][4]; and lastly, creep-fatigue damage of interconnecting solder joints after repeated cycles of powering on/off of engineering devices [2][5]. The last is especially critical and has attracted maximum interest in electronic assemblies. In essence, the differential thermal expansions between the silicon chip and the substrate and between the integrated circuit assembly and the printed circuit board give rise to cyclical deformation of the solder joints at temperature above their homologous temperature, driving them towards failure by the mechanism of creep fatigue. The propensity for creep fatigue failure of solder joints is aggravated by the trend towards increasing functionality of consumer electronic products, which is driving increase size of integrated circuit assembly and reduce size of solder joints.

The failure of a nuclear plant or a commercial aircraft is accompanied by unacceptable catastrophic consequences. The structural integrity and reliability of such products are therefore the paramount design considerations. In contrast, the two paramount design considerations for electronic assemblies are electrical performance and space. The former to support the ever-increasing performance of electronic products and the latter to support the increasing functionality of electronic products. Failure of consumer electronic products in service though annoying is acceptable to most users. This has encouraged a relatively relax attitude towards the structural integrity of electronic assemblies; and this is further encouraged by the relatively short life cycle of consumer electronic products. Nevertheless, there is a positive aspect of this more relax attitude towards structural analysis. The structural design of electronic assemblies is not bounded by a design protocol or a design code. Electronic assembly engineers are free to use any analysis method so long as the designed electronic assemblies will meet the integrity and reliability test requirements. It is inherently easy to monitor the structural integrity of an electronic assembly through monitoring the electrical connectivity of the assembly [6]. If necessary, the growth of damage in an electronic assembly can be tracked in real time through monitoring the changing electrical impedance of the assembly [6][7]. The absence of a strict design protocol and the ease of validating an analysis with tests has encouraged the exploration and adoption of new analysis methods, notwithstanding some of these methods may not be robust.

This manuscript gives a condensed presentation of the recent advances in the creep fatigue analysis of solder joints in electronic assemblies. This comprises two major analyses: thermomechanical analysis of solder joints in electronic assemblies; and creep fatigue life modelling of solder joints. It is believed that these advanced analyse techniques are fundamentally robust and they can be adopted to similar applications in other engineering field.

II. THERMOMECHANICAL ANALYSIS OF SOLDER JOINTS

Electronic assembly engineers routinely performed thermomechanical analysis of electronic assemblies using finite element analysis software in which the solder joints are modelled using solid finite elements. This has inadvertently led to singularity of stress/strain at discontinuities of geometry and materials giving rise to inconsistent analysis - because the magnitude of the stress/strain is dependent on the size and shape of the finite element at the site of singularity. To circumvent such singularity, electronic assembly engineers have adopted the practice of volume-averaging the stress/strain over a selected volume of solder joints [8][9]. Unfortunately, such arbitrary volume averaging act is equivalent to smearing the geometry of solder joints, effectively denying the engineers the ability to analyse the geometrical effects of solder joints on the magnitude of stress/strain. The issue of stress/strain singularity can be addressed by modeling the components of electronic assemblies as shells and beams. Maximum insights into the mechanics of the subject matter can be achieved through analytical modeling.

Analytical models of varied sophistication have been reported in the electronic assembling community [10][11][12][13]. In essence, an electronic assembly is treated as a sandwich structure constituting of array of solder joints sandwiched between two outer members. The simplest model treats the outer members as being infinitely rigid and the solder joints as being infinitely compliant such that the solder joints experience only shear strain whose magnitude increases linearly with distance from the neutral plane of an electronic assembly and is given by $\gamma =\mathrm{\Delta }\alpha T\cdot x/h$ (Figure 2a) [10], where $\mathrm{\Delta }\alpha T$ is the differential thermal strain between the two outer members, $x$ is the distance of a discrete joint from the neutral plane of the assembly, and $h$ is the height of the discrete joint. This unrealistic model would grossly overestimate the magnitude of shear strain. The more sophisticated models treat the outer members as rectangular beams; and the solder joints as linear springs [11], or as cylindrical beams that are capable of shearing, flexing, and stretching [12][13]. The inclusion of elasticity of the outer members in the model have led to the vital understanding that shear strain in the solder joints does not increase linearly but exponentially with distance from the neutral plane of an electronic assembly. It takes the form $\gamma ={\gamma }_o{e}^{\beta (x-l)}$ (Figure 2b), where ${\gamma }_o$ is approximately the shear strain at the outmost solder joint, $\beta$ is a compliance parameter of the assembly, $l$ is the half-length of the assembly.

a) Analytical Modelling

Evaluating the stresses in the discrete joints of an assembly involves four steps of analysis, starting with smearing the discrete joints into a continuously bonded joint and then evaluating the compliances of the smeared assembly. This is followed by evaluating the smeared stresses in the smeared joints. The third step integrates the stresses into boundary forces and moments acting on individual discrete joint. The last step evaluates the stresses in the discrete joint due to the boundary forces and moments.

The discrete joints are assumed to be of identical shape and size and are distributed at uniform spacing. Figure 3 shows the schematic of discrete joints with a height $h_d$ and spacing at pitches $p_x$ and $p_y$ along the $x$ and the $y$ coordinates, respectively. In the case that individual discrete joint is not of cylindrical shape but one with non-uniform sections along its height, it is represented by a pseudo cylindrical joints with a representative shear area, $A_{d,rep}$, and a representative second moment of area, $I_{d,rep}$, that would return identical shear and flexural stiffnesses as the original discrete joints [14]. Theses representative parameters are given by

Referring to the outer members as member #1 and member #2, and the discrete joint as member $\mathrm{\#}d$, the height, the stretch modulus, the shear modulus, and the flexural stiffness of member $\mathrm{\#}i$ are denoted as $h_i$, $E_i$, $G_i$, and $D_i$, respectively; the shear and the in-plane stretch compliances of the assembly are denoted as ${\kappa }_s$ and ${\lambda }_x$ respectively. For ease of reference, we shall refer to the moduli of the smeared joints as smeared moduli and the compliances as smeared compliances. Those characteristics that are associated with the smeared joints will be marked with an asterisk.

b) Compliances of a Smeared Assembly

The in-plane stretch compliance of the assembly is a function of the outer members and, for the case of plane stress, is given by [12][13]

where in $D_i=E_i{h}_i^3/12$ for plane stress. The smeared shear compliance of the assembly, ${\kappa }_s^{\ast }$, is given by [12][13]

where ${\kappa }_{si}=\frac{h_i}{8G_i},{\kappa }_{sd}^{\ast }=\frac{h_d}{G_d^{\ast }},{\kappa }_{sd}\phi =\frac{h_d^3}{12D_d^{\ast }}.$ (4) wherein ${\kappa }_{si}$ is associated with the shear compliance of member #i; ${\kappa }_{sd}^{\ast }$ is associated with the shear deformation of the smeared joints; while ${\kappa }_{sd\phi }$ is associated with the flexural deformation of the discrete joints - referring to Figure 2b. The smeared shear modulus, $G_d^{\ast }$, and the smeared flexural rigidity, $D_d^{\ast }$, are given by $G_d^{\ast }=G_d{A}_{d,rep}/(p_x{p}_y)$ and $D_d^{\ast }=E_d{I}_{d,rep}/(p_x{p}_y)$.

c) Stresses in Smeared Joints

The shear stress in the smeared joints along the bonded length of a balanced assembly is given by [12][13][15][16]

where

${\beta }^{\ast }=\sqrt{{\lambda }_x/{\kappa }_s^{\ast }};{\epsilon }_I=\mathrm{\Delta }\alpha T$ is the differential thermal strain between the outer members.

d) Sectional Force and Moment in a Discrete Joint

The magnitude of the sectional shear force on a discrete joint that is at a distance $\mathrm{\Gamma }$ from the mid plane of the assembly may be evaluated approximately as

The sectional shear force does not vary along the height of the discrete joint. On the other hand, and referring to Figure 4, rotational equilibrium dictates that the sectional moment varies linearly along the height of the discrete joint and is given by

where $z$ is the local coordinate of a solder joint as shown in Figure 3.

e) Shearing and Bending Stresses in a Discrete Joint

The distribution of shear stress and bending stress along the height of a discrete joint, assuming it being an Euler beam, are simply [12][13]:

where $A_d(z),I_d(z)$, and $r_d(z)$ are the local cross-sectional area, the local second moment of area, and the local outer fibre of the discrete joint. Assuming circular cross-section, as in the case of solder joints, the ratio $r_d/I_d$ is reduced to $4/(\pi r_d^3)$. The largest magnitudes of shear force and bending moment, and hence shearing and the bending stresses, occur at the discrete joint furthest from the mid-plane of the assembly. Assuming ${\mathrm{\Gamma }}_{\max }=l$, these stresses are given by

where $F_{\tau ,max}=\frac{{\epsilon }_T{p}_x{p}_y}{\sqrt{{\lambda }_x{\kappa }_s^{\ast }}},$ (11) and $A_{d,min}$ is the minimum cross-sectional area of the discrete joint.

f) Optimum Shape of Solder Joints

The solder joints are formed through controlled heating of solder into liquid form followed by controlled cooling the assembly to room temperature, forming metallurgical bonds with metal pads at its two ends. A solder joint will take up the natural shape of a spherical barrel, as shown in Figure 5a, that has the minimum surface energy. It is clear from Eq. (10) that a standard solder joint will experience the maximum magnitude of shear stress, ${\tau }_{d,max}$, and the maximum magnitude of bending stress, ${\tau }_{b,max}$, at its ends joining the outer members; and

where $r_{d,\text{end}}$ is the radius of the solder joint joining the outer member. In a standard barrel-shape solder joint, the magnitude of $h_d$ is much larger than that of $r_{d,\text{end}}$. In other words, ${\sigma }_{b,\text{max}}$ is a far dominant stress in a standard barrel-shape solder joint. On paper, the magnitude of the dominant stress, ${\sigma }_{b,\text{max}}$, can be lowered by increasing the end radius of barrel-shape solder joints, which for the same volume of solder joint, will result in flatten barrel-shape solder joints leading to very significant reduction in the magnitude of ${\sigma }_{b,\text{max}}$. In practice, such a manipulation would inevitably raise the risk of bridging between adjacent solder joints, as is illustrated in Figure 5b.

It is clear from the linear distribution of bending moment in solder joints, as depicted in Figure 4, that solder joints should ideally have the shape of an hourglass. Solder joints of progressive hourglass shape can be designed into electronic assembly [14]. Modeling the curvature of an hourglass-shape solder joint as a hyperbolic sine curve and evaluating its representative shear area, $A_{d,rep}$, and representative second moment of area, $I_{d,rep}$, using Eq. (1), a stout hourglass-shape solder joint – similar to that illustrated in Figure 5c – has been found to return the minimum magnitude of stress for the same volume of solder as a standard barrel-shape solder joint. The maximum magnitude of bending stress, ${\sigma }_{b,max}$, in the stout hourglass solder joints is less than $15\mathrm{\%}$ that in a standard barrel-shape solder joints [14]. It is also worth noting that the use of stout hourglass solder joints does not increase the risk of bridging between solder joints.

III. CREEP FATIGUE MODELING

The fatigue life of a metal experiencing pure low cycle fatigue – that is, in the absence of aggravating element – and under a constant amplitude of cyclic stressing has been found to be satisfactorily modelled using the Coffin-Manson equation:

where ${\epsilon }_p$ is the amplitude of the incremental plastic strain in a cycle; $N_f$ is the number of cycles to failure; $C_o$ and ${\beta }_o$ are material dependent fitting constants. While it is tempting to extend the equation to creep fatigue modelling by replacing the plastic strain amplitude, ${\epsilon }_p$, with inelastic strain amplitude, ${\epsilon }_{in}={\epsilon }_p+{\epsilon }_c$, wherein ${\epsilon }_c$ is the incremental creep strain in a single cycle, this simplicity approach of lumping creep strain with fatigue strain contradicts with the different macrostructural damages of creep and fatigue in metals [17][18] and has been convincingly disproved by abundant experimental data [19].

a) A Brief Review of Practising Creep Fatigue Life Prediction Models

The power generation community and the aerospace engineering community have vast knowledge and experience in modelling creep fatigue life of metals. Both the power generation and the aerospace communities have subscribed to the idea that creep fatigue damage can be evaluated by summing independently the damages due to creep, due to fatigue, and due to interaction of these two damages. The summative creep fatigue damage in a single creep fatigue cycle may be expressed mathematically as

Interestingly, the two communities have subscribed to different idea of defining the respective damages, $d_f,d_c$, and $d_{cf}$.

The aerospace community characterises the three cyclic damage indices from the hysteresis loop of a tension-compression creep fatigue experiment. Three components of inelastic strain: ${\epsilon }_{pp}$, ${\epsilon }_{cc}$, ${\epsilon }_{cp}$ (or ${\epsilon }_{pc}$ ) are partitioned and extracted from the hysteresis loop, wherein the first letter of the subscript ( $c$ for creep and $p$ for plastic strain) refers to the type of strain imposed in the tensile portion of the cycle, and the second letter refers to the type of strain imposed during the compressive portion of the cycle. Individual strain component is assumed to follow a power-law relation with the number of hysteresis cycles to failure. That is,

The damage per creep-fatigue cycle due to individual component is then simply

This is known as the strain range partitioning model [19]. While this noble model has served the aerospace community well, its characterisation is inherently challenging.

The power generation community conveniently treats the cyclic fatigue damage, $d_f$, as that due to pure fatigue, which can be evaluated using the Coffin-Manson equation; that is,

and the cyclic creep damage, $d_c$, as that due to pure creep, which may be evaluated using the creep strain exhaustion rule:

wherein ${\dot{\epsilon }}_c(t)$ is the instantaneous creep strain rate and $t_c$ the cyclic period. However, the cyclic creep-fatigue interaction damage, $d_{cf}$, is a fitting index, which can only be established through extensive experimental characterisation [20][21][22].

The electronic packaging professionals have fallen for the creep fatigue equation of Darveaux [9],

where $w_{in}=w_p+w_c$ is the sum of the cumulative plastic work density and the cumulative creep work density in a single creep fatigue cycle. Just like the failed idea of substituting plastic strain with inelastic strain in Eq. (13), the act of lumping the two work densities is clearly against the macrostructural evidence of the two damages. Consequently, and unsurprisingly, the fitting constants, $K_1$ and $K_2$, are found to be dependent on the size and shape of individual electronic assembly [9], in other words, on the magnitude of the inelastic work density, $w_{in}$. Nevertheless, the electronic packaging community have stubbornly stuck with the model.

b) Creep Integrated Fatigue Equation

In the case of fatigue being the dominant mechanism in creep fatigue failure, the role of creep may be treated as one to lower the material capacity in fatigue. This has led to the idea of creep integrated fatigue equation [23][24][25]:

where $c({\epsilon }_p,T,t_c)$ is a function. Expressing the fatigue capacity in Eq. (13) for the case of pure fatigue as ${\epsilon }_{p,ref}$ and it becomes clear that $c({\epsilon }_p,T,t_c)={\epsilon }_p/{\epsilon }_{p,ref}$ describes the fractional fatigue capacity of a subject in the presence of creep. Its magnitude ranges from zero to unity - a zero magnitude corresponds to the case of pure creep while a magnitude of unity corresponds to the case of pure fatigue. The function 1- $c({\epsilon }_p,T,t_c)$ describes the fractional creep damage acting on the subject.

i. The Fundamental Equations of Pure Creep In Metals

The strain rate in a uniaxial steady-stress creep rupture experiment may be described in the form of Sherby-Dorn equation, ${\dot{\epsilon }}_{SD}$, or Larson-Miller equation, ${\dot{\epsilon }}_{LM}$, or Manson-Haferd equation, ${\dot{\epsilon }}_{MH}$:

wherein $f({\sigma }_s),H({\sigma }_s)$, and $r({\sigma }_s)$ are functions of the applied tensile stress, ${\sigma }_s$; $k$ is the Boltzmann's constant; $H,B,$ and $D$ are material dependent constants; and $T_{ref}$ is the temperature below which the mechanism of creep is assumed to be dormant. Assuming the dominance of the secondary stage of creep, the eventual creep rupture strain, ${\epsilon }_R$, is given by

where $t_R$ is the time to creep rupture. Assuming ${\epsilon }_R$ to be independent of the applied stress and temperature, Eq. (21) may be rearranged into:

wherein $B_R=B/{\epsilon }_R$ and $t_{\mathrm{\infty }}={\epsilon }_R/D$. These functions are known respectively as the Sherby-Dorn parameter, the Larson-Miller parameter, and the Manson-Haferd parameter. These parameters describe the relations between the applied stress, ${\sigma }_s$, and the gradient of the respective time-temperature function. These relations can be readily characterized and are used extensively in the creep rupture design of metal structures.

The corresponding stress parameter for a single stressing cycle may be expressed as [25]:

where $\sigma$ is the amplitude of the cyclic stress; and $B_c=B/{\epsilon }_c$, $t_{c\mathrm{\infty }}={\epsilon }_c/D$, and ${\epsilon }_c$ is the cumulative creep strain in a single cycle. Assuming (i) identical creep damage due to tensile and compressive stresses, and (ii) linear cumulation of creep strain over varied magnitudes of stress, then the cumulative creep strain in a single cycle may be expressed as ${\epsilon }_c={\int }_0^{t_c}k(|\sigma (t)|,T)dt$, where the function $k(|\sigma (t)|,T)$ represents one of the rate equations of the $SD$, $LM$ and $MH$; and $|\sigma (t)|$ is the instantaneous magnitude of the cyclic stress. The cyclic stress parameter function may then be evaluated mathematically from the steady-stress parameter function as [25]

where ${\varphi }_{\epsilon }={\epsilon }_c/{\epsilon }_R$

ii. Fractional Fatigue Capacity Function

Let the fractional fatigue capacity function $c({\epsilon }_p,T,t_c)$ takes the form:

Herein $\chi ({\epsilon }_p)\eta (T,t_c)=\epsilon$ is the fractional creep damage function. It is intuitive that the function $\eta (T,t_c)$ shall take the form of the rate equation of creep; that is,

It is worth noting that $\eta (T,t_c)$ vanishes at $T\le 0$ for ${\eta }_{SD}$ and ${\eta }_{LM}$ and at $T\le T_{ref}$ for ${\eta }_{MH}$ when the mechanism of creep becomes dormant; the condition of pure fatigue prevails and the fractional fatigue capacity function $c({\epsilon }_p,T,t_c)$ acquires the maximum magnitude of unity. Similarly, the condition of pure creep requires that the magnitude of fractional fatigue capacity function vanishes to nil. This implies, from Eq. (24), that

wherein the subscript $k$ signifies $SD$, $LM$, and $MH$, respectively. The cyclic parameter function, $P_{c\mathrm{\_}k}(\sigma )$, can be evaluated from the steady-stress parameter function, $P_k({\sigma }_s)$, using Eq. (25). Using the Ramberg-Osgood relation, $\sigma =$

$\overline{K}{\epsilon }_p^{\overline{n}}$, where $\overline{K}$ and $\overline{n}$ are assumed to be the representative material constants over the range of temperature of interest, the cyclic parameters may then be expressed as a function of plastic strain; that is,

The material dependent fitting constants, $C_o,{\beta }_o,\overline{K},\overline{n}$, and ${\varphi }_{\epsilon }$ (the constant $H$ may be evaluated from the steady-stress creep rupture test data) or $B_c$ or $t_{c\mathrm{\infty }}$, may be established through regressing the experimental creep fatigue data with Eq.(20). It is worth mentioning that the constants $\overline{K},\overline{n}$, and ${\varphi }_{\epsilon }$ or $B_c$ or $t_{c\mathrm{\infty }}$ represent the damages due to pure creep and also creep fatigue interaction.

In practice, one can do away with the experimental characterisation of the steady-stress parameter function, $P_k({\sigma }_s)$, and the subsequent mathematical evaluation of the cyclic parameter function, $P_{c\mathrm{\_}k}(\sigma )$. By expressing the steady stress parameter function as a power law in the form $P_k({\sigma }_s)=p{\sigma }_s^q$, where $p$ and $q$ are fitting constants, it can be shown that the cyclic stress parameter function will be reduced to a power law function, $P_{c\mathrm{\_}k}(\sigma )=\hat{p}{\sigma }^{\hat{q}}$, where $\hat{p}$ and $\hat{q}$ are constants, for the sinusoidal and the triangular stress-time profile. This could then be transformed to $P_{c\mathrm{\_}k}({\epsilon }_p)$ using the Ramberg-Osgood relation. In other words, the function $\chi ({\epsilon }_p)$ may simply be assumed to be a power-law function,

where $a$ and $b$ are arbitrary constants. These two arbitrary constants can be established together with three material dependent fitting constants $C_o$, ${\beta }_o$, and $H$ or $B_c$ or $t_{c\mathrm{\infty }}$, through regressing the experimental creep fatigue data with Eq. (20). This shall be illustrated in the following section.

iii. Illustration and Validation

The experimental creep fatigue data of Sn37Pb under cyclic triangular stress-time stressing generated by Shi et al. [26] has been analysed and the fitting constants, $C_b$ and ${\beta }_b$, for nine sets of $({\epsilon }_p,N_{cf})$ data have been extracted and these are tabulated in Table 1 [23]. Using these fitting constants, the $({\epsilon }_p,N_{cf})$ data were regenerated and these are depicted in Figure 6.

It is clear from Eq. (20) that a set of creep fatigue data, $({\epsilon }_p,N_{cf})$, can be transformed into a set of pure fatigue data, $({\epsilon }_{p,ref},N_f)$, by transforming the magnitude of ${\epsilon }_p$ into ${\epsilon }_{p,ref}$ using the relation,