A critical review of constitutive models for solders in electronic packaging

Abstract

Owing to their superior electrical, thermal, and mechanical properties, solder joints are the most widely used interconnection materials in electronic product packaging. Because the failure of the whole electronic packaging is often induced by the failure of solders, modeling and simulation of solder joint performance are quite important in ensuring the quality and reliability of electronic products. In fact, the accuracy of reliability design of electronic packaging is dependent on accelerated reliability tests, material constitutive models, finite element modeling, and fatigue-life prediction models. Here, the material constitutive model plays the most important role in the development of electronic reliability design because the stress, strain, and strain energy density adopted in finite element modeling and fatigue-life prediction are all computed and derived from the constitutive model. In this work, the constitutive models developed for simulating the complex stress and strain response of solders are critically reviewed. To make the description concise and comprehensive, the constitutive models were categorized into two types: the constitutive models with and without the definition of yield surface. The prediction capabilities, application scope, merits and shortcomings, and objective suggestions for the further development of constitutive models for solders in electronic packaging are proposed.

Introduction

Solder is a fundamental material used in the assembly of electronic devices to printed circuit boards (PCBs). There are three important functions of the solder interconnections: mechanical support, electrical connection, and thermal dissipation. During the service of electronic devices, the package/PCB assembly is subject to changes in temperature owing to cyclical power on/off status and computational activity. In such temperature variations, thermally induced stresses, arising from the dissimilar coefficients of thermal expansion (CTE) of packaging materials, could cause the solder joints to fail, and eventually compromise the functionality of the whole packaging assembly. Along with the progress of science and technology, the ever-increasing demands of high performance, multi-function, and miniaturization for electronic devices have introduced a considerable challenge in reliability design of solder joints. Therefore, the thermomechanical reliability of solder joint in microelectronic packaging assembly needs to be carefully studied.

The design of thermomechanically reliable microelectronic packages is a time-consuming and costly task. The accelerated thermal cycling and power cycling tests are elaborate. In general, modeling methods have become a useful tool to replace the accelerated tests and predict the thermomechanical reliability of solder interconnections. The modeling methods basically include three steps: (1) establish a constitutive model to predict the complex stress–strain response of the key materials in electronic packaging assembly, for example, solder, plastic module, and substrate; (2) implement the constitutive models into finite element software, create the model of electronic packaging assembly, and predict the mechanical response of the solder joints and surrounding assembly under thermomechanical cycling; (3) apply an appropriate fatigue model to predict the thermomechanical reliability of the assembly. Irrespective of the kind of fatigue models involved, the necessary input information such as stress, strain, strain energy density, or plastic work, can be entirely extracted from the finite element analysis (FEA) simulation results. In summary, the precision of numerical methods depends on the quality of constitutive model, finite element method, and fatigue model (failure criterion). In order to make a reliable computational simulation, it is necessary and important to develop an accurate constitutive model to simulate the material behavior. ¹

The following is a critical review on constitutive models for solders in electronic packaging, which were reported in international publications in recent years. Here, to make the description concise and comprehensive, we reclassify the models into two types:

1. Constitutive models without the definition of yield surface. In general, the constitutive models without the definition of yield surface can be categorized into two types: unified constitutive models and uncoupled constitutive models. The unified models consider the sum of plastic strain and creep strain as the unified inelastic strain, which results in simplicity, whereas the uncoupled models consider the plastic strain and creep strain separately. For the unified constitutive models, the notable feature is that rate dependence and temperature dependence of solders can be well captured using a complicated flow rule; therefore, the models are usually used to simulate the uniaxial static test, including tension, creep, and stress relaxation. For the uncoupled models, the main merit of the classified models is straightforward physical-based significance. However, the distinction of plastic strain and creep strain is difficult to be determined during the test, especially for complex loading conditions. Thus, the uncoupled models discussed here are only related to creep simulation.

2. Constitutive models with the definition of yield surface. Most of these types of constitutive models are involved in the framework of classical plastic theory, that is, yield criterion, plastic flow rule, and hardening rule. Loading and unloading criteria can be determined according to Drucker’s rule and yield surface definition. The flow rule defines the evolution of the inelastic strain and follows the consistency condition. The hardening rule includes isotropic and dynamic hardening rule, which can be used to describe the cyclic hardening/softening or the Bauschinger effect. Comparing with the constitutive models without the definition of yield surface, these models allow for an accurate definition of material behavior under wide ranges of loading conditions by some physical assumptions, but require a larger number of material constants.

Models without the definition of yield surface

Unified constitutive models

Anand model

The Anand model is one of the most widely used constitutive models in electronic packaging to predict the complex stress–strain relationship of solders under different temperatures and strain rates. There are two specific features in the Anand model. First, this model needs no explicit yield condition and no loading and unloading criteria. Second, the model employs a single scalar as an internal variable s which expresses the averaged isotropic deformation resistance to macroscopic plastic flow.^2,3 The Anand model can consider the temperature sensitivity, strain rate sensitivity, temperature and strain rate history effects, strain hardening, and the dynamic recovery. The functional form of the flow equation, which accommodates the strain rate and temperature dependence on the stress, can be expressed as follows

$${\stackrel{·}{\epsilon }}^p=A\exp (-\frac{Q}{\mathrm{R}T}){\left[\sinh \left(\xi \frac{\tilde{\sigma }}{s}\right)\right]}^{1/m}$$

(1)

where ${\stackrel{·}{\epsilon }}^p$ is the inelastic strain rate, A is the pre-exponential factor, m is the strain rate sensitivity, Q is the activation energy, R is the universal gas constant, T is the absolute temperature, $\xi$ is the material constant, and s is the internal state variable. The evolution equation for the internal variable s is as follows

$$\stackrel{·}{s}=\left\{h_0{|1-\frac{s}{s*}|}^a\mathrm{sign}(1-\frac{s}{s*})\right\}{\stackrel{·}{\epsilon }}^p$$

(2)

where h₀ and a are material parameters representing the strain hardening. s* represents the saturation values of s and can be expressed as follows

$$s*=\hat{s}{\left[\frac{{\stackrel{·}{\epsilon }}^p}{A}\exp \left(\frac{Q}{\mathrm{R}T}\right)\right]}^n$$

(3)

where $\hat{s}$ is a coefficient of saturation, and n is the strain rate sensitivity for the saturation value of deformation resistance. The nine material parameters, such as A, Q, m, n, $\xi$ , $\hat{s}$ , a, h₀, and s₀ can be easily determined from a limited number of experimental tests.

Lots of scholars implemented the Anand model to simulate the stress–strain relationship of a variety of solders successfully.^4–13 However, Chen et al. ³ found that the Anand model cannot accurately predict the response of material which has the specification of strong strain-hardening effect at low temperature. Thus, they proposed a modified Anand model by correlating h₀ with temperature and strain rate

$$h_0=a_0+a_{1}T+a_2{T}^2+a_3*{\stackrel{·}{\epsilon }}_p+a_4*({\stackrel{·}{\epsilon }}_p{)}^2$$

(4)

where a₀, a₁, a₂, a₃, and a₄ are material constants and can be determined by the nonlinear fitting method. ¹⁴ From Figure 1, one can see that the simple amendment can effectively improve the prediction ability of Anand model. It is worth noting that the modified Anand model is readily available in ABAQUS as a user routine. The modified model and its parameters are widely adopted to predict the SAC305 solder.^15–17

Figure 1.

Comparison of experiments and predictions for lead-free solder Sn–3.5Ag: (a) original Anand model and (b) modified Anand model. ³

Although the Anand model is suitable for predicting the stress–strain relationship of solders under uniaxial monotonic loading at different temperatures and strain rates, the model is not feasible for predicting the response of solders subjected to complex loading conditions, such as multiaxial loading or cyclic loading. ¹⁸ This is because, in the Anand model, only one single scalar is employed as an internal variable to express the averaged isotropic deformation resistance to macroscopic plastic flow. Therefore, the Bauschinger effect and multiaxial stress–strain response cannot be characterized by the Anand model due to the lack of definition of dynamic hardening rule.

Bodner–Partom model

Bodner and Partom proposed a constitutive model, the Bodner–Partom model (BP model), to describe elastic-visco-plastic strain-hardening material behavior for large deformation and arbitrary loading histories. ¹⁹ An essential feature of BP model formulation is that the total deformation rate is considered to be separated into elastic and inelastic components. The state variables in the BP model are used to describe the mechanical behavior of the materials during loading and unloading. There is no yield criterion in the BP model. The plastic strain flow rule is given as follows

$${\stackrel{·}{\epsilon }}_p=\frac{2}{\sqrt{3}}\left(\frac{\sigma }{\left|\sigma \right|}\right)D_0\exp [-\frac{n+1}{2n}{\left(\frac{Z}{\sigma }\right)}^{2n}]$$

(5)

where $\sigma$ is stress, D₀, and n are the material parameters and D₀ is the limiting value of the strain rate based on the maximum loading rate. It was recommended that for quasi-static modeling, D₀ = 10⁴ s⁻¹ and for dynamic modeling, D₀ = 10⁶ s⁻¹. n is the strain-sensitive coefficient. Z is the internal invariable which accounts for the history effect and can be presented as

$$Z=Z_1+(Z_0-Z_1)\exp (-m\int \sigma d\frac{{\epsilon }^p}{Z_0})$$

(6)

where Z₀ is the initial value of Z, Z₁ is the maximum or saturation value of Z, and m is the material parameter. Methods to obtain the BP material parameters D₀, n, Z₀, Z₁, and m have previously been discussed by various authors, for example, Bodner and Partom. ¹⁹

Skipor et al. ²⁰ described the uniaxial tensile curves of 63Sn–37Pb by BP model under a wide range of strain rates and temperatures. These are shown in Figure 2. The results indicated that the predicted results only agreed well with the experimental results at higher strains or higher stress levels. The degree of strain rate sensitivity exhibited by the model was approximately fairly correct. However, the BP model failed to represent the stress and strain response of the solder at lower strain levels or at lower temperatures. The poor fits are undesirable, but perhaps unavoidable, results from highly strain rate sensitivity at the lower strain levels.

Figure 2.

Experimental and computed curves of true stress versus true strain for the four applied strain rates at (a) −40°C, (b) −10°C, (c) 20°C, (d) 60°C, and (e) 100°C. ²⁰

Whitelaw et al. ²¹ modified the BP model by adding a dynamic recovery term, a thermal recovery term, and a temperature-rate-dependent term. It was found that a reasonable representation of the time and temperature-dependent nature could be estimated using the modified model. In addition, the modified model was shown to predict the cyclic behavior of the alloys well, especially at the elevated temperatures. However, the BP model could not readily predict the stress–strain relationships of the solders at low temperature or low strain rate. In addition, the determination of the parameters based on microscopic observation hardly yielded a repeatable and consistent stress–strain relationship of the materials because the microstructure possessed some degree of uncertainty.

VBO model

The “unified” viscoplasticity theory based on overstress (VBO), developed by Ho and Krempl, ²² is used to simulate strain rate dependence, creep, relaxation, cyclic hardening/softening, as well as ratcheting (namely, plastic strain accumulation due to cyclic stressing with nonzero mean stress).^23–27 The VBO model does not use a yield surface and loading/unloading criterion. For isothermal conditions, the flow rule can be written as follows

$$\stackrel{·}{\epsilon }={\stackrel{·}{\epsilon }}^{\mathrm{el}}+{\stackrel{·}{\epsilon }}^{\mathrm{in}}=\frac{1+v}{\mathrm{E}}\stackrel{·}{s}+\frac{3}{2}\frac{s-g}{\mathrm{E}k(\Gamma )}$$

(7)

$${\Gamma }^2=\frac{3}{2}(s-g):(s-g)$$

(8)

$$k=k_1{\{1+\frac{\Gamma }{k_2}\}}^{-k_3}$$

(9)

where $\stackrel{·}{\epsilon }$ is the total strain rate, ${\stackrel{·}{\epsilon }}^{\mathrm{el}}$ is the elastic strain rate, ${\stackrel{·}{\epsilon }}^{\mathrm{in}}$ is the inelastic strain rate, v is Poisson’s coefficient, E is the elastic modulus, s is the deviatoric stress, g is the equilibrium stress, and k₁, k₂, k₃ are material constants.

Maciucescu and colleagues^28,29 proposed a minimal version of the VBO to predict the behavior in tensile and cyclic strain-controlled loadings, ratcheting, and creep of Sn–Pb solder alloy (shown in Figure 3). The comparison between experimental and predicted results indicated that the VBO model can predict the inelastic response of solder alloy very well. However, because one of the main differences between the VBO model and the general constitutive models is that when formulating the back-stress evolution in terms of total strain rate instead of viscoplastic strain rate, it is difficult to incorporate the VBO model into a standard thermodynamic framework. ³⁰

Figure 3.

Comparison of experimental results (a) with VBO predictions and (b) of cyclically stable hysteresis loops in torsion for several strain amplitudes at 0.01 Hz and 25°C. ²⁸

In summary, it can be seen that the Anand model is suitable for describing the temperature dependence and rate dependence of solders in a wide range due to its powerful flow equation in hyperbolic sine form. In addition, it is easy to determine the parameters in the Anand model from a series of uniaxial tensile tests. However, the Anand model cannot accurately represent the cyclic stress–strain response of solders. For the BP model, it can be seen that the model is capable of describing the uniaxial monotonic and cyclic stress–strain response at high temperature and high strain rate, but the invalidity at low temperature and low strain rate limits its applications. On the other hand, because the parameters of BP model are determined based on microscopic observation, it is difficult to obtain a repeatable result due to the deviation of microstructural observation. For the VBO model, it can be used to predict the complex stress–strain response of solders very well under uniaxial and multiaxial loading conditions, by adding more material parameters. The effect of drag stress in VBO model is quite similar to that of back stress defined in classical plastic theory, by which the Bauschinger effect and dynamic hardening effect can be well described. However, as mentioned above, the fatal shortcoming of VBO model makes it difficult to be incorporated into a standard thermodynamic framework due to the replacement of total strain rate to plastic strain rate in the formulation of drag stress evolution.

Uncoupled models for creep simulation

A typical conventional creep curve includes three stages of creep, namely, primary creep, steady-state creep, and tertiary creep.^31,32 Steady-state creep, also referred to as secondary creep, takes up a relatively longer period of creep time than with primary and tertiary creep because the creep rate is constant. Thus, constant creep rate for the secondary creep, resulting from a balance between the competing processes of strain hardening and recovery, is generally referred to as the creep rate. ⁴

Dorn model (Norton power law)

The steady-state creep rate is generally estimated by constitutive models such as those generally described by Weertman–Dorn equation ²¹ as follows

$$\stackrel{·}{\epsilon }=\frac{A\mathrm{G}b}{\mathrm{k}T}{\left(\frac{b}{d}\right)}^p{\left(\frac{\sigma }{\mathrm{G}}\right)}^n{D}_0\exp (-\frac{Q}{\mathrm{k}T})$$

(10)

where $\epsilon$ is the steady-state creep strain, $\stackrel{·}{\epsilon }$ is the creep strain rate, A is a constant, G is the shear modulus, b is the Burgers vector, k is Boltzmann’s constant, T is the absolute temperature (in Kelvin), d is the grain size, p is the material’s grain size exponent, σ is the applied stress, n is the stress exponent, D₀ is the frequency factor, Q is the activation energy of the diffusion, and C is a constant. The above equation can be simplified as the Dorn model (Norton power law)^21,33–36 given by

$$\stackrel{·}{\epsilon }=C_1{\sigma }^n\exp (-\frac{Q}{\mathrm{k}T})$$

(11)

where $C_1$ is another constant. Such a model provides a useful description of the steady creep rate as a function of temperature, but is seriously limited in applicability because the creep strain rate is an exponential function of stress rather than a power law function for high values of stress. Creep behavior in this high stress region is expressed as follows

$$\stackrel{·}{\epsilon }=C_1\exp (C_2\sigma )$$

(12)

where $C_1$ and $C_2$ are material constants.

The constitutive relations are of practical interest for modeling of solder joints. However, STM/BGA/CSP joints are all small enough that solidification of each nucleated in only one point. The joint therefore has cyclic twin boundaries, but no real grain boundaries. Thus, lots of constitutive relations have been based on measurements on much larger, polycrystalline samples, leading to grain size effects that are not relevant to joints. In addition, the inelastic deformation properties are strongly dependent on the distribution of secondary precipitates and precipitate coarsening for thermal cycling where precipitate coarsening is important. Constitutive relations that ignore the microstructural evolution of precipitate distribution and size cannot accurately predict the stress–strain response of solder joints under combined effect of thermal and mechanical stress. In order to consider the microstructural coarsening effect, Dutta and colleagues^37–43 proposed a microstructurally adaptive creep model by incorporating Ansell–Weertman model as follows.

For low stress regime

$${\stackrel{·}{\epsilon }}_c\propto \left(\frac{\mathrm{G}b}{\mathrm{k}T}\right)r{\left(\frac{\sigma }{\mathrm{G}}\right)}^n\exp (-\frac{Q}{\mathrm{k}T})$$

(13)

For high stress regime

$${\stackrel{·}{\epsilon }}_{\mathrm{plc}}\propto \left(\frac{\mathrm{G}b}{\mathrm{k}T}\right){\left(\frac{\sigma -{\sigma }_{\mathrm{th}}}{\mathrm{G}}\right)}^n\exp (-\frac{Q}{\mathrm{k}T})$$

(14a)

$${\sigma }_{\mathrm{th}}=\frac{\sqrt{2(1-{\mathrm{k}}^2)}}{\sqrt[3]{\pi V_f/6}}\frac{\mathrm{G}b}{2r}$$

(14b)

where ${\stackrel{·}{\epsilon }}_c$ is the creep rate, ${\stackrel{·}{\epsilon }}_{\mathrm{plc}}$ is the creep rate for particle-limited climb, G is the shear modulus, b is the Burgers vector, k is Boltzmann’s constant, T is the absolute temperature, r is the particle size, Q is the activation energy of the diffusion, ${\sigma }_{\mathrm{th}}$ is the threshold stress, and $V_f$ is the volume fraction of the intermetallic particles. It can be seen that the creep rate is proportional to the particle size in the low stress regime, and is proportional to r in the high stress regime; thus, the effect of particle size on creep in both the low and high stress regimes can be captured by equations (13) and (14a). The deformation mechanism for the microstructurally adaptive creep model is clear. The most attractive merit of model is that only the initial particle size in typical joints needs to be measured, and the variation of the particle size in cycling and aging can be reasonably predicted as proposed by Kutta and his coworkers.^41–43

Garofalo model (hyperbolic sine-law)

In order to capture both the power law behavior and the exponential behavior, Darveaux and Banerji ⁴ proposed a hyperbolic sine creep equation as follows

$$\frac{d{\gamma }_s}{\mathrm{dt}}=C_s\frac{\mathrm{G}}{T}{\left[\sinh \left(\alpha \frac{\tau }{\mathrm{G}}\right)\right]}^n\exp \left(\frac{-Q}{\mathrm{k}T}\right)$$

(15)

where $C_s$ is a material constant, $\alpha$ represents the stress level at which the power law dependence no longer holds. The equation was later modified by Darveaux ⁴⁴ to include the transient creep strain that occurs immediately after the stress is applied and before a steady-state creep stage is reached. Darveaux’s constitutive model has been implemented into the commercial finite element software, for example, ANSYS, ABAQUS, and is widely used in the analyses of modeling solder joints in electronic packages, that is

$$\frac{d{\gamma }_c}{\mathrm{dt}}=\frac{d{\gamma }_s}{\mathrm{dt}}(1+B{\gamma }_T\exp (-B\frac{d{\gamma }_s}{\mathrm{dt}}))$$

(16)

Zhang et al. ⁴⁵ predicted the creep behavior of 80Au/20Sn solder alloy using the Dorn model and Garofalo model. Figure 4 shows that both models are able to provide an acceptable description of eutectic 80Au/20Sn solder alloy experimental results over the present experimental stress-temperature range. However, the Garofalo model exhibits a lower estimated variance of error terms than the Dorn model.

Figure 4.

Comparison of steady-state creep rate obtained from experimental results and fit models: (a) Dorn model and (b) Garofalo model. ⁴⁵

Similar to the models proposed by Dutta and colleagues,^37–40 Tang et al. ⁴⁶ proposed a modified Garofalo model by incorporating the internal stress to predict the creep behavior of 95.8Sn–3.5Ag–0.7Cu shear lap. The internal stress was correlated to the shear stress in the low-stress range, and to the particle size and volume fraction of the intermetallic particles in the high-stress range, which makes the modified Garofalo model simulate the creep behavior of 95.8Sn–3.5Ag–0.7Cu solder joint at all different temperatures and stress levels reasonably well (shown in Figure 5). However, the authors did not present how to determine the microscopic parameters during the creep process, and the sensitivity analysis of the parameters on simulation results was lack.

Figure 5.

Measured shear stress as a function of creep strain rate behavior compared with modified creep constitutive model predictions at different test temperatures. ⁴⁶

Lee and Basaran ⁴⁷ proposed a modified creep model to account for shear modulus temperature dependency and grain size effect. The model is denoted by the following formula

$$\frac{d{\gamma }_s}{\mathrm{dt}}=A{\left(\frac{b}{d}\right)}^p{\left[\sinh \left(\alpha \frac{\tau }{\mathrm{G}}\right)\right]}^n\exp \left(\frac{-Q}{\mathrm{k}T}\right)$$

(17)

The proposed model described the stress versus strain rate curves very well as shown in Figure 6.

Figure 6.

Creep model for different grain sizes: (a) on grain size 9.7 µm and (b) on grain size 28.4 µm. ⁴⁷

Amalu and Ekere ⁴⁸ evaluated the effect of the use of different parameters in the Garofalo-Arrhenius constitutive creep relation on the predicted magnitude of creep strain and number of cycles to failure (fatigue life). Four sets of parameters, proposed by Lau, ⁴⁹ Pang et al., ⁵⁰ Schubert et al., ⁵¹ and Zhang et al., ⁵² were employed in their work. The results indicated that different values of parameters produced distinctive magnitudes and history of stress, strain, hysteresis loop, accumulated strain, accumulated plastic strain, and fatigue life (in Figure 7), even if the same set of solder joints in electronic assemblies and the same creep model are used.

Figure 7.

Equivalent creep strain in solder joints showing: (a) schematic of strain damage in resistor solder joint, (b) schematic of strain damage in flip chip solder joint, (c) plot of equivalent creep strain against temperature cycle time for resistor assembly, and (d) plot of equivalent creep strain against temperature cycle time for flip chip assembly. ⁴⁸

Creep-mechanism-based model

Creep mechanisms can be categorized into three types: dislocation glide, dislocation creep, and diffusion creep. Unlike the above Dorn and Garofalo models, which are phenomenologically based, the creep-mechanism-based model is physically based. Compared to the phenomenological descriptions, they allow for an accurate definition of material behavior under wide ranges of loading conditions by some physical assumptions and a larger number of material constants:

• 1. Knecht and Fox model

Knecht and Fox ⁵³ used the following creep form

$${\stackrel{·}{\gamma }}_{\mathrm{cr}}=C_1{\tau }^{n_{\mathrm{GB}}}\exp (-\frac{\Delta H_{\mathrm{GB}}}{\mathrm{k}T})+C_2{\tau }^{n_{\mathrm{MC}}}\exp (-\frac{\Delta H_{\mathrm{MC}}}{\mathrm{k}T})$$

(18)

where k is Boltzmann’s constant; GB and MC denote grain boundary and matrix creep, respectively; and $\Delta H$ refers to the activation energies for the corresponding creep mechanism.

Knecht and Fox ⁵³ used representative values of the activation energies and exponents found in the literature to predict the creep data. The results showed that the model was insensitive to the selection of the parameters, all values assumed for $\Delta H_{\mathrm{GB}}$ , $\Delta H_{\mathrm{MC}}$ , $n_{\mathrm{GB}}$ , $n_{\mathrm{MC}}$ fit the creep data quite well (see Figure 8). However, the difficulties of determining the parameters based on creep mechanism limit the application of the model.

Figure 8.

Shine and Fox steady-state creep data. ⁵³

• 2. Heiduschke model

Heiduschke ⁵⁴ used the creep part as the inelastic strain, and proposed the following model

$${\stackrel{·}{\epsilon }}^c=f^c(\nu )n$$

(19)

$$n=\frac{{\epsilon }_e\prime }{\Vert {\epsilon }_e\prime \Vert },\Vert {\epsilon }_e\prime \Vert =\sqrt{{\epsilon }_{\mathrm{ij}}{\prime }^e{\epsilon }_{\mathrm{ij}}{\prime }^e}$$

(20)

$$f^c(\nu )=f^{\mathrm{DC}}(\nu )+f^{\mathrm{GBS}}(\nu )$$

(21)

$$f^{\mathrm{DC}}(\nu )=A_{\mathrm{DC}}{(\Vert {\epsilon \prime }^e\Vert )}^{n_{\mathrm{DC}}}\exp (-\frac{Q_{\mathrm{DC}}}{\mathrm{R}T})$$

(22)

$$f^{\mathrm{GBS}}(\nu )=A_{\mathrm{GBS}}{(\Vert {\epsilon \prime }^e\Vert )}^{n_{\mathrm{GBS}}}\exp (-\frac{Q_{\mathrm{GBS}}}{\mathrm{R}T})$$

(23)

where $\epsilon {\prime }^e$ is the deviatory elastic strain tensor; DC and GBS denote dislocation climb and grain boundary sliding, respectively; $Q_{\mathrm{DC}}$ and $Q_{\mathrm{GBS}}$ denote corresponding activation energies; R is the universal gas constant; and T is the absolute temperature. The model was incorporated into a FEA program and used to simulate the eutectic Sn–Pb solder joints of quad flat packs with gull wing under thermal cyclic loading.

• 3. Shi model

Shi et al. ⁵⁵ proposed a creep model as follows:

$${\stackrel{·}{\gamma }}_{\mathrm{cr}}=C_1{\left(\frac{\tau }{\mathrm{G}}\right)}^{n_1}\exp (-\frac{Q_1}{\mathrm{R}T})+C_2{\left(\frac{\tau }{\mathrm{G}}\right)}^{n_2}\exp (-\frac{Q_2}{\mathrm{R}T})$$

(24)

This model assumes that there are two regimes of steady-state creep, and each regime has a power law dependence of creep strain rate on stress. The results have shown that the model provides a reasonable prediction of creep deformation at the intermediate temperature regime, but underestimates the creep deformation of the solder at low temperature (−40°C) and high temperatures (125°C and 150°C). In addition, the determination of the stress exponent and activation energy used in the second term is elaborate because they are not obtained from the same experimental data.

Based on their own experimental results, the same researchers proposed a modified model as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{\gamma }}_{\mathrm{cr}}=}& C_1\frac{\mathrm{G}}{T}{\left[\sinh \left(\beta \frac{\tau }{\mathrm{G}}\right)\right]}^{n_1}\exp (-\frac{Q_1}{\mathrm{R}T})\hfill \\ \multicolumn{1}{c}{}& +C_2\frac{G}{T}{\left[\sinh \left(\beta \frac{\tau }{\mathrm{G}}\right)\right]}^{n_2}\exp (-\frac{Q_2}{\mathrm{R}T})\hfill \end{array}$$

(25)

The first term in both models is for low temperature dislocation-glide-dominated creep, while the second term is for high temperature dislocation-climb-dominated creep. With this unified dislocation-controlled creep constitutive model, one can predict very well the creep deformation of solder alloys at different temperatures for the entire stress range, as shown in Figure 9.

Figure 9.

Comparison of the creep deformation obtained from the experiments with that obtained from unified dislocation-controlled creep model. ⁵⁵

Models with the definition of yield surface

In this section, the constitutive models with the definition of yield surface are discussed. The models listed here focus on describing the cyclic behavior of solders under complex loading conditions, such as Bauschinger effect, cyclic hardening or softening, thermomechanical fatigue (TMF), and damage evolution during fatigue process. In fact, these constitutive models are also capable of simulating the monotonic tensile test and creep test as well.

McDowell model

McDowell creep plasticity unified constitutive model,^56,57 developed based on the phenomenological theory, is capable of simulating the rate dependence and temperature dependence of metal alloys. This is especially true when describing the creep fatigue interactions, strain-induced dynamic recovery, static thermal recovery, and cyclic aging of tin-based solders used in electronic packaging. ⁵⁶ A two back-stress version, that is, short-range back stress characterizes the slow evolution of initial yielding behavior and long-range back stress that contributes in the asymptotic plastic flow regime, is enough to satisfy the demands of prediction accuracy. ⁵⁸

The temperature-dependent viscoplastic flow rule has the form

$${\stackrel{·}{\epsilon }}^{\mathrm{in}}=A{\left(\frac{<S_v>}{D}\right)}^n\exp \left[B{\left(\frac{<S_v>}{D}\right)}^{n+1}\right]\theta N$$

(26a)

$$\theta =\exp \left(\frac{-Q}{\mathrm{R}T}\right),\mathrm{for}T\ge \frac{T_m}{2}$$

(26b)

$$\theta =\exp \left\{\frac{-2Q}{\mathrm{R}{T}_m}[\ln \left(\frac{T_m}{2T}\right)+1]\right\},\mathrm{for}T\le \frac{T_m}{2}$$

(26c)

where A, D, n, and B are all material constants; S_v is over stress; N is the normal direction of yield surface; Q is the energy of thermal activation; and R is the gas constant; the Macauley brackets $<>$ are defined as $<x>=((x+|x|)/2)$ .

Many scholars used the McDowell model to describe the stress–strain relationship of tin-based solders.^58–60 Fu et al. ⁵⁸ presented numerical simulations of different test cases, which included monotonic tensile test, creep tests, and two-step ratcheting tests, for 62Sn–36Pb–2Ag solder. The predicted results by McDowell model agree well with the experimental results if a semi-implicit time-integration scheme was used in the constitutive level of computation (Figure 10). Due to lack of the experimental data of tin-based solder, only the numerical simulations for oxygen-free high thermal conductivity (OFHC) copper under a non-proportional and an in-phase TMF history were presented to verify the validity of McDowell constitutive model and the integration scheme.

Figure 10.

Second sequence of a two-step ratcheting test at 20°C for 62Sn–36Pb–2Ag and the FE simulations. ⁵⁸

Neu et al. ⁵⁹ proposed a new McDowell unified creep-plasticity (UCP) model through modifying the short-range back-stress evolution equations (Figure 11). They accurately predicted the deformation phenomena associated with the solders’ monotonic tension, strain jump tests, steady-state creep, and cyclic responses by taking the stress-dependence of the activation energy and the strong Bauschinger effect into account. ⁵⁹ In addition, Busso et al. ⁶¹ presented an effective method to measure the back stress through the strain transient tests and the rapid loading/unloading tests, and the latter method provided a better measurement because the time elapsed to determine that the measurement was short enough to avoid the influence of the creep or thermal recovery.

Figure 11.

Comparison of (a) experimental hysteresis loops at a strain rate of 10⁻³ with those (b) predicted by the McDowell UCP model. ⁵⁹

Wen et al. ⁶⁰ proposed a damage-coupled McDowell UCP model to simulate the cyclic hardening and softening, damage evolution during fatigue process, and TMF of high lead 96.5Pb–3.5Sn solder (Figures 12 and 13). It can be seen that the model can predict the isothermal and thermomechanical stress–strain response and the damage accumulation of the solder very well.

Figure 12.

Model predictions of the stress–strain curve for the first two cycles (a) and (b) under various loading and temperature condition for solder 96.5Pb–3.5Sn. ⁶⁰

Figure 13.

Model results from UCP with damage for 96.5Pb–3.5Sn solder under TMF tests: (a) experimental data and (b) model comparison. ⁶⁰

More recently, Yao and Keer ⁶² adopted the dual-back-stress model based on McDowell model to predict the stress and strain response of solder materials through a finite element software. In their work, thermal strain as well as the isotropic hardening was all considered to predict the solder interconnections under coupled thermomechanical loading. The proposed model was extended to predict the dynamic condition under high strain rate. It was found that the inelastic strain predicted using the proposed flow rule is much higher than elastic and thermal strain because of the large deviatory stress; thus, only inelastic strain was mainly considered in relative ductile Sn–3.8Ag–0.7Cu eutectic solder. The nominalized maximum von Mises stress in the solders at different rows was plotted in Figure 14. The figure indicates that rate dependence of solder interconnections can be captured well by the McDowell model under drop impact.

Figure 14.

Nominalized maximum von Mises stress in the solders at different rows. ⁶²

Busso model

Based on the continuum theory, Busso et al. ⁶¹ presented a visco-plastic model which accounted for the stress-dependence of the activation energy and the strong Bauschinger effect for Sn/40Pb. The equations are as follows

$${\stackrel{·}{\epsilon }}^p={\stackrel{·}{\epsilon }}_0\exp \{-\frac{F_0}{\mathrm{R}\theta }{[1-{\left(\frac{|\sigma -B|}{{\sigma }_0\mathrm{G}/{\mathrm{G}}_0}\right)}^p]}^q\}\mathrm{sign}(\sigma -B)$$

(27)

where $F_0$ is the total free energy of activation at a vanishingly small stress, B is back stress, ${\sigma }_0$ is yield stress at 0 K, G is the shear modulus, ${\mathrm{G}}_0$ is the shear modulus at zero Kelvin degree, and p, q are material constants.

Busso et al. ⁶¹ properly predicted the temperature and strain range dependence of the shear stress range of 60/40 tin lead solder, as shown in Figure 15. The Bauschinger effect, thermal mechanical behavior, and ratcheting can be also simulated very well with the model. ⁶³ In addition, they proposed a self-correcting forward gradient time integration scheme to implement the Busso model into finite element software, for example, ABAQUS. However, the Busso model does not exhibit a smooth elastic–plastic transition after yielding, and the isotropic hardening and static recovery effect were also ignored to keep the model simple.

Figure 15.

Model prediction of cyclic stress–strain behavior in shear. ⁶¹

Ohno–Wang model

Ohno and colleagues^64–66 used a multilinear model in which several kinematic hardening rules of the Armstrong–Frederick type were superposed. They assumed that each component of back stress ${\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}$ had a critical state for its dynamic recovery term to be activated. The Ohno–Wang model is proposed in the following form

$$\mathbf{\alpha }\mathbf{=}\sum _1^{\mathbf{M}}{\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}$$

(28)

$$\mathrm{Model}(\mathrm{I}){\stackrel{·}{\mathbf{\alpha }}}^{(i)}=h^{(i)}[\frac{2}{3}{\stackrel{·}{\epsilon }}^p-H(f^{(i)})〈{\stackrel{·}{\epsilon }}^p:\frac{{\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}}{{\overline{\mathbf{\alpha }}}^{(i)}}〉\frac{{\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}}{r^{(i)}}]$$

(29a)

$$f^{(i)}={\overline{\mathbf{\alpha }}}^{(i)}-{r^{(i)}}^2$$

(29b)

$$\mathrm{Model}(\mathrm{II}){\stackrel{·}{\mathbf{\alpha }}}^{(i)}=h^{(i)}[\frac{2}{3}{\stackrel{·}{\epsilon }}^p-{\left(\frac{{\overline{\mathbf{\alpha }}}^{(i)}}{r^{(i)}}\right)}^{m^{(i)}}〈{\stackrel{·}{\epsilon }}^p:\frac{{\mathbf{\alpha }}^{(i)}}{{\overline{\mathbf{\alpha }}}^{(i)}}〉\frac{{\mathbf{\alpha }}^{(i)}}{r^{(i)}}]$$

(30)

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{Model}(\mathrm{III}){\stackrel{·}{\mathbf{\alpha }}}^{(i)}=}& {\zeta }^{(i)}\{\frac{2}{3}{r}^{(i)}{\stackrel{·}{\epsilon }}^p-[(1-{\mu }_i)H(f^{(i)})+{\mu }_i]\hfill \\ \multicolumn{1}{c}{}& 〈{\stackrel{·}{\epsilon }}^p\mathbf{:}\frac{{\mathbf{\alpha }}^{(i)}}{{\overline{\mathbf{\alpha }}}^{(i)}}〉{\mathbf{\alpha }}^{(i)}\}\hfill \end{array}$$

(31)

where M denotes the number of parts of back stress; ${\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}$ is the ith component of the deviatory back stress $\mathbf{\alpha }$ ; ${\overline{\mathbf{\alpha }}}^{(i)}$ is the magnitude of ${\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}$ , ${\overline{\mathbf{\alpha }}}^{(i)}=\sqrt{3/2{\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}:{\mathbf{\alpha }}^{\mathbf{(}\mathbf{i}\mathbf{)}}}$ ; $h^{(i)}$ and $r^{(i)}$ are material constants; H stands for the heaviside step function; <> is the Macauley bracket, that is, $<x>=(x+|x|)/2$ . ${\zeta }^{(i)}=h^{(i)}/r^{(i)}$ , ${\mu }_i$ is the material constant. The parameters in the Ohno–Wang model can be easily determined by a uniaxial tensile stress–strain curve ⁶⁷ and the model is not sensitive to yield stress value and the number of back-stress components. ⁶⁸ Therefore, the Ohno–Wang model was widely used to simulate the uniaxial tension, strain rate dependence, cyclic hardening/softening, and uniaxial and multiaxial ratcheting.

Chen and colleagues^69–71 employed a constitutive model with the Ohno–Wang kinematic hardening rule to simulate the isothermal cyclic behavior of Sn–Pb solder under uniaxial, torsional, and multiaxial loading. They also simulated the rate-dependent multiaxial ratcheting of the 63Sn–37Pb by correlating the material parameter ${\mu }_i$ to the strain rate. Comparisons of the experimental and simulated results verify that the constitutive model with the Ohno–Wang kinematic hardening rule is able to predict the complicated deformation of 63Sn–37Pb very well (Figures 16 and 17).

Figure 16.

Equivalent stress–strain loop with damage under uniaxial and torsional loading: (a) axial and (b) torsional. ⁷⁰

Figure 17.

Axial strain versus cycles for specimen SN21: (a) experimental and (b) modified model. ⁷¹

Based on the radial return method and backward Euler integration, Chen et al. ¹⁸ implemented the viscoplastic Ohno–Wang constitutive model into finite element software ABAQUS to simulate the shear and the ratcheting behavior of the sintered nanosilver joint. The simulated results obtained by Ohno–Wang model were also compared with the Anand model, as shown in Figure 18. It can be concluded that Ohno–Wang type model could predict the ratcheting behavior of the sintered nanosilver joint better than the Anand model, especially at high temperatures. However, the Ohno–Wang model is incapable of accounting for the temperature dependence of the materials. ⁷² Therefore, it is necessary to determine the parameters from uniaxial tension curve at each temperature. As far as the parameter determination is concerned, the Anand model is far more convenient compared to the Ohno–Wang model because one set of parameters is enough for the Anand model to simulate the temperature and strain-rate sensitivity of the materials.

Figure 18.

Prediction relationship between ratcheting displacement and number of cycles of sintered lap-shear joint with different dwell time by different models: (a) 225°C and (b) 325°C. ¹⁸

Damage-coupled constitutive model

Solder joints are assumed to be homogeneous and defect-free in most conventional mechanics approach. However, the presence of micro-defects, which are unavoidable during the manufacturing process, are the main cause of crack initiation, growth, and fracture when the materials subjected to static or cyclic loading. The progressive material degradation has been identified as the main cause of solder failure. ⁷³ Therefore, the damage models, taking into account the material degradation and independence from the constitutive model, have been developed based on the theory of damage mechanics by many scholars.^74–83 Stolkarts et al. ⁷⁵ proposed a typical damage evolution law as follows

$$D=\frac{16(1-{(\overline{v})}^2)(10-3\overline{v})}{45(2-\overline{v})}\omega$$

(32)

$$\overline{v}=-0.579\omega +0.37$$

(33)

where D is the damage parameter, $\overline{v}$ is Poisson’s coefficient, and $\omega$ is the microcrack density. As pointed out by Stolkarts et al., ⁷⁵ the number of cycles to failure in the nucleation governed damage process was determined as follows

$$N_f=\frac{{\overline{\omega }}_c}{{(\Delta {\gamma }^{\mathrm{cycle}})}^{\chi }}$$

(34)

where ${\overline{\omega }}_c$ can be taken as 0.093 for a relatively large specimen corresponding to its characteristic length, $\Delta {\gamma }_i^{\mathrm{cycle}}$ is an average shear strain range at a slip plane along a grain or phase boundary per cycle, and $\chi$ is a material constant. It is assumed that the microcrack density $\omega$ increases at a constant rate during steady state of damage process, that is

$$\omega =(\Delta {\gamma }^{\mathrm{cycle}}{)}^{\chi }N$$

(35)

where N is the cyclic number. Stolkarts et al. ⁷⁵ arbitrarily correlated $\Delta {\gamma }^{\mathrm{cycle}}$ to the total strain range used in the uniaxial test. It was assumed that for continuous fatigue loading, $\Delta {\gamma }^{\mathrm{cycle}}=k\Delta \epsilon$ , where k was the proportionality constant. Figure 19 shows the experimental and predicted results of 63Sn–37Pb by the damage model proposed by Stolkart et al. under different strain amplitudes. However, the model can be only put into application for uniaxial loading condition.

Figure 19.

Modeling of fatigue for various loading conditions. ⁷⁵

In order to expand the application scope of the damage model proposed by Stolkart, Chen et al. ⁷⁰ proposed a modified damage model, which replaced the uniaxial strain range with the sum of maximum shear strain range and normal strain on the critical plane, to predict the damage evolution process of tin–lead solder under multiaxial low cycle fatigue. Based on the critical plane concept, Chen et al. believed that the normal strain will assist microcrack growth, and the expression of $\Delta {\gamma }^{\mathrm{cycle}}=k\Delta \epsilon$ was changed into

$$\Delta {\gamma }_i^{\mathrm{cycle}}=k(\Delta {\gamma }_{max}+S\Delta {\epsilon }_n)$$

(36)

where k is the proportionality constant, $\Delta {\gamma }_{max}$ is the maximum shear strain range, and $\Delta {\epsilon }_n$ is the normal strain range on the maximum shear strain critical plane. S is a material constant and can be determined as 0.5 from the experiments. With this modification, they accurately predicted the peak stress evolution with the number of cycles under various loading paths as shown in Figure 20.

Figure 20.

Equivalent stress amplitude evolution with the number of cycles under various loading conditions: (a) axial and (b) torsional. ⁷⁰

The damage evolution equations brought forward above are independent of the constitutive model, and any solder constitutive model can be used in conjunction with them. However, the numerous parameters in those models tend to restrict the application of the model only to engineers, and many tests have to be performed for the procedure of parameter determination.

In summary, all the constitutive models with the definition of yield surface are established based on the framework of classical plastic theory. Therefore, these models can be easily implemented into finite element software for further structural stress analysis. For the McDowell model, the flow rule is taken as the hyperbolic sine form, which makes it applicable for describing the stress–strain response of solders in a wide range of temperatures and strain rates. For the Busso model and Ohno–Wang model, the flow rule is formulated in the exponential and power form, which means that the viscosity of solders cannot be accounted for well, especially for a high temperature and low strain rate. When it comes to the dynamic hardening rule, only one or two components are included in the back-stress formulation for the Busso model and the McDowell model, whereas at least four components of back stress are included in the Ohno–Wang model. This means that the Ohno–Wang model has the highest capability of predicting the nonlinear stress–strain behavior of solders among the three mentioned models, especially for the case of non-proportional multiaxial loading conditions.

Conclusion

The goal of constitutive modeling for solders in electronic packaging is to provide the precise information of stress, strain, and strain energy density for electronic reliability evaluation and design. The major constitutive models practiced in the literature have been reviewed and these could be categorized into two types:

1. Models without the definition of yield surface. The constitutive models without definition of yield surface are suitable for simulating the static tests of solders, namely, the tensile test and creep test. The strain-rate sensitivity and temperature sensitivity of solders can be well simulated by these models, namely, the Anand model, BP model, Dorn model and VBO model. Most of the constitutive models listed here have been already implemented into commercial finite element software and are freely available as user subroutines. Therefore, it is convenient to employ these models into applications for evaluating the reliability of electronic devices using finite element method. However, the thing that needs to be noted is that most of the first type of constitutive models here are unable to accurately predict the cyclic response of solders because of the lack of yield surface definition and dynamic hardening rule. Most of the first type of models does not incorporate a hardening rule, that is, isotropic or dynamic hardening rule, except for the VBO model.

2. Models with the definition of yield surface. The second type of constitutive models, such as the McDowell model, Busso model, and Ohno–Wang model, all focus on describing the cyclic behavior of solders under complex loading conditions, such as the Bauschinger effect, cyclic hardening or softening, TMF, ratcheting, and damage evolution during fatigue process. The isotropic hardening rule and dynamic hardening rule along with a complex flow rule are developed to simulate the cyclic hardening effect, the Bauschinger effect, and the viscous effect of solders. It is noted that the flow rule included in the first type of constitutive models can also be used here to predict the rate dependency and temperature dependency of solders. With the combination of yield criterion, flow rule and hardening rule, the second type of constitutive models are capable of representing the stress–strain relationship of solders very well, even under complex loading conditions, which is similar to the real stress state in service. However, because the viscous part is too weak to reflect the complex inelastic behavior of solders, the second type constitutive models are difficult to accurately describe the strain rate and temperature sensitivity of solders simultaneously, especially for low strain rates and high temperatures.

There are many aspects that were incompletely developed in this study, with many good references and other original models. An impression at the end is that these current models have made a significant contribution toward predicting the inelastic behaviors of solders for reliability design. Further developments rely on a combination of phenomenological-based models and physical-based models, that is, multiscale constitutive models, to predict the complex deformation of solder joints accurately.