Lifetime prediction of a viscoplastic lead-free solder in power electronics modules under passive temperature cycling

– It is well-established now that solder joint in power electronics undergoes a thermal-mechanical fatigue in service. In order to provide a predictive tool to assess fatigue lifetime of solder joint subjected to passive temperature cycling, a cohesive zone method (CZM) based fatigue model, accounting for mixed-mode loadings, is presented in this paper. A cohesive zone model with a speciﬁc fatigue traction-separation law is used in conjunction with the Anand viscoplastic behavior for the bulk solder for modelling the interfacial fatigue crack growth. Temperature dependence of both stiﬀness and fracture energy is incorporated in the CZM based fatigue model. The CZM based fatigue model as well as the Anand viscoplastic constitutive law is implemented in the Abaqus ﬁnite element code. Numerical simulations are carried out under passive temperature cycling. A particular attention is given to information in the region ahead of the crack tip such as stress concentration, process zone size, damage distribution. The evolution of the crack extension is used to estimate the number of cycles to failure at which the device becomes faulty.


Introduction
Thermomechanical reliability of power electronics components is an important issue for hybrid and electric vehicle industry.In service, electronic components are subjected to temperature fluctuations caused by both heat dissipation of silicon chips and environmental changes.Combined with the thermal expansion mismatch between the different materials of the assembly (silicon chip, solder joint and copper substrate), cyclic thermal loadings result in stress reversals and accumulation of inelastic strain in the solder joint.This inelastic strain accumulates with repeated cycling and ultimately causes solder joint cracking and interconnect failure [1].During more than 50 years, SnPb alloys are widely used for electronic interconnections, and so far simulation of SnPb solder joint reliability is quite well developed and applied.The transition from SnPb to a new lead-free technology is a challenging and demanding task for the companies.It has an impact on material supply, process equipment and conditions but also the reliability will change.Several issues are showing up which are typically for leada Corresponding author: zhidan.sun@utt.frfree solder joints and make simulation and also reliability testing much more complicated [2].
In the case of IGBT (Insulated-Gate Bipolar Transistor) modules (Fig. 1), solder joint delamination is one of the main failure modes along with wire bond cracking.The mechanical behavior of these systems under fatigue has been widely studied experimentally and numerically.However, prototyping and testing of electronic components devices are often expensive and time consuming while numerical modeling can be more efficient in the design and optimization steps.Various methods, based on empirical laws, have been used for predictions of solder joint reliability.For example, in Coffin-Manson type of fatigue laws, the width of the hysteresis loop, which is an estimate of the inelastic strain range that the solder joint experiences, is used to obtain the number of cycles to failure.In Darveaux's model [3], the damage evolution is based on the inelastic strain energy density per cycle.Although these models provide a valuable engineering predictive framework, they do not give a clear description of the physical process of crack development in the solder joint.
The cohesive zone method (CZM), originally introduced by Dugdale [4] and Barenblatt [5], has been used as an efficient tool to analyze delamination problems under monotonic loadings [6,7].In the basic concepts for cohesive zone models, the fracture process is lumped into the crack line and is characterized by a traction-separation law that relates cohesive stress and crack opening across the cohesive surface.The CZM emerges as an attractive approach for fatigue loading.Some models were presented in the literature [8][9][10][11][12][13] to analyze different fatigue failure problems.Under cyclic loadings, the situation is rather different.Roe and Siegmund [9] introduced a damage evolution law to account for the irreversible accumulation of damage through the cycling process.The inclusion of unloading-reloading hysteresis into the cohesive law is done by employing a condition based on an analogy to an elastic-plastic material under the consideration of damage.The same concepts were applied in a mode II (pure shear) version of the CZM for fatigue life prediction of solder joints under isothermal conditions [10].The application of the CZM for simulating fatigue cracking of solder joint is still in its early stage.Yang et al. [10] considered rate-independent plasticity for representing solder behavior and implemented a pure-shear (mode II) cohesive law to simulate the damage evolution.In fact, the normal loading (mode I) also plays a role in the damage evolution process and its contribution may be significant, especially in the case of IGBT modules.In addition, in the case of solder alloys, the temperature dependent properties should be taken into account.Indeed, the effect of temperature changes during cyclic thermal loadings may be relatively important for solders compared with their very low melting temperature.Accordingly, the temperature would affect both the cohesive strength and the damage accumulation process.
In this paper, a CZM based fatigue model derived from energy potentials is used to analyze the thermomechanical fatigue cracking of solder joint.This fatigue model includes the mixed-mode effect, as well as the temperature dependence of the model's parameters such as energy to fracture.The mixed-mode effect is considered by introducing equivalent parameters which combine pure tensile and shear mode loadings.A cumulative damage variable is introduced in order to simulate the cracking under cyclic loading.The traction-separation equations are irreversible and history dependent.The fracture energy associated to each mode is taken as a function of both temperature and separation rate.The structure of this paper is as follows: the CZM based fatigue model is firstly presented.Then, the parameters for the fatigue model and the constitutive law of a SnAgCu lead-free solder are detailed.Finally, the numerical simulations are presented and analyzed, and lifetime predictions are proposed for a system subjected to passive temperature cycling.

CZM based fatigue model
A mode dependent CZM based fatigue model is presented in this section.Modes I and II fracture parameters are combined with a mixed-mode failure criterion based on energy.A global damage variable is incorporated in the cohesive zone constitutive law in order to account for the degradation process.The damage variable is supplemented by an evolution law to take into account the damage evolution.
Let's consider the current thermodynamic potential of the cohesive zone as follows: where k t , k n and k comp are the stiffnesses respectively under shear, traction and compression.D is a global damage variable.The symbol • is defined as x = 1 2 (x + |x|).The compression stiffness k comp is generally taken very high in order to penalize any interpenetration of matter.The cohesive stress vector T is the partial derivative of the thermodynamic potential with respect to the displacement jump vector Δ: ( The effective displacement jump Δ and the equivalent cohesive stress T in the mixed mode are calculated respectively by: and with α = kt kn .α measures contribution of displacement jumps in modes I and II to the equivalent displacement jump.The current damage evolution threshold for mixed mode, f = 0, is chosen as follows: where Y m is the mixed mode thermodynamic force, conjugate to the damage variable D; G i m is the equivalent energy to damage initiation, and G c m represents the equivalent total energy to fracture.G i m ensures the existence of a damage initiation threshold under which fracture can never initiate.Y m , G i m and G c m are calculated from the two pure modes accounting for the mode mixity which depends on the parameter α and the loading angle.The damage evolution Ḋ is then obtained by applying the consistency condition ḟ = 0: The second equation states that the damage variable does not evolve during unloading.The detailed formulation of the model is presented in reference [14].
3 Material parameters

Parameters of the fatigue model
As described above, the set of input parameters for the CZM based fatigue model includes the stiffness k x , the damage threshold energy G i x and the fracture energy G c x for each mode x (in tensile and shear loadings).The tensile and shear stiffnesses can be derived from the elastic modulus E and the shear modulus G of solder alloys by introducing a constitutive thickness for the cohesive zone [14].For solder alloys, the elastic modulus E is both temperature and strain rate dependent.A formula for calculating the elastic modulus of a Sn-3.8Ag-0.7Cuwas experimentally established by Pang and Xiong [15]: As for the fracture energies G c x , they can be estimated from the stress-strain curves under pure tensile (mode I) and shear loadings (mode II).For the simulations carried out in this work, the fracture energy densities in mode I and in mode II are derived from the work of Pang and Xiong [15].In their work, tensile and shear stress-strain curves are presented as function of temperature and strain rate for a Sn-3.8Ag-0.7Cusolder alloy.For each condition of temperature and strain rate, the volume density of fracture energy is obtained by estimating the surface beneath the corresponding curve.Formulas of volume density of fracture energy dependent on temperature and strain rate are thus derived respectively for tensile loading condition: 704 ln ( ε) + 5.591Q 0.0684 + 6.998 (8) and for shear loading condition: The bonds of validity for these functions correspond to a temperature range from 25 • C to 125 • C, and to a strain rate range from 2.5 × 10 −4 s −1 to 2.5 × 10 −2 s −1 .For use in the cohesive zone formulation, the volume density of fracture energies are then transformed into surface density of fracture energies using the following relations: where h CZ denotes the constitutive thickness of the cohesive layer.As regards the damage threshold energy density G i x , the value is expected to be very small compared to G c x in the case of solder alloys.In this work its value is taken as 2% of the fracture energy density G c x .

Constitutive behavior of the bulk solder
The constitutive behavior of the bulk solder is modeled with Anand viscoplastic model.This model was proposed initially by [16], and completed later by Brown et al. [17] The set of Anand constitutive equations accounts for the phenomena of strain-rate and temperature dependence, strain rate history effects, etc.In this constitutive model, the following functional form of the flow equation was used to accommodate the strain rate dependence on the stress at constant temperature: The internal state variable s enters into the flow equation as a ratio with the equivalent stress σ eq .The temperature dependence is incorporated via a classical Arrhenius term, while the stress and state dependence are a modification of the hyperbolic sine form.Figure 2 illustrates the influence of temperature and equivalent stress on the strain rate, which represents the steady-state creep behavior of the solder alloy.The evolution equation for the internal variable s is assumed to be of the form: where the function h (σ, s, T ) is associated with dynamic process.A simple form of the evolution relation was given as follows: with where h 0 is the hardening/softening constant, a is the strain rate sensitivity of hardening/softening.The scalar s * represents a saturation value of s associated with a set of given temperature and strain rate, while ŝ is a coefficient, and n is the strain rate sensitivity for the saturation value of deformation resistance.In the above viscoplastic Anand model, nine material parameters are required.In this work, the values of these parameters are taken from [18] for a Sn-3.8Ag-0.7Cusolder alloy, as shown in Table 1.Note that s 0 is the initial value of the deformation resistance s.
It should be noted that the Anand model is a unified viscoplastic model and was originally formulated to be used more for hot working process rather than for thermal cyclic loading.Two features for modeling thermal cyclic loading are not considered in the Anand model: (1) kinematic hardening and (2) aging term in the evolution equation.According to [19], the importance of kinematic hardening is somewhat reduced at high homologous temperature and relatively low strain rate regimes.As for the aging of lead-free solder alloys, it is still not presently sufficiently studied.Some further studies on the aging of lead-free solder alloys are needed in the future.
As the Anand viscoplastic law is not available in the Abaqus finite element code, it is defined with the user subroutine CREEP which provides a capability for implementing viscoplastic models.This subroutine can be used to define time-dependent, viscoplastic deformation in a material.The deformation is divided into deviatoric behavior (creep) and volumetric behavior (swelling).The strain rate can be written as a function of Mises equivalent deviatoric stress, and any number of solution-dependent state variables.The values of the solution-dependent state variables evolve with the solution and can be defined in the subroutine.Nonlinear equations must be solved at each time step and the variations of different variables should be defined in the subroutine.To obtain good convergence during implicit integration, it is essential to define these quantities accurately.The subroutine must use the corresponding values of time, field variables, and solution-dependent state variables in the calculation of the creep strain increment.
In Figure 3, the constitutive behavior of the Sn-3.8Ag-0.7Cualloy is presented as function of temperature according to the prediction from Anand model described above.The steady-state plastic flow occurs when the deformation resistance s equals to the saturation value s * .
The other components such as the die (made up of silicon) and the substrate (made up of Cu-Mo) are assumed to behave as linear isotropic elastic materials, with temperature independent properties.Finally note that the assembly is assumed to be initially stress-free at a reference temperature of 20 • C.

Simulation and analysis of results
In this section, the simulation with the simplified assembly composed of the three components die-soldersubstrate is carried out.The assembly is subjected to a cyclic temperature change from 20 to 80 • C with a period of 200 s (Fig. 4a), which is representative of the service condition.The temperature is applied uniformly to the entire assembly.This type of loading permits to simulate the thermomechanical fatigue under passive temperature cycling conditions.Only one half of the assembly is modelled due to the symmetry existing in the geometry and the corresponding boundary conditions are applied (Fig. 4b).Table 1.Values of parameters of Anand viscoplastic model for a Sn-3.8Ag-0.7Cusolder alloy [18].broken, i.e. a crack is formed.The details of stress concentration can be seen on the zoomed view of the crack tip region.Along the crack path, one can identify the point A before which the cohesive elements are intact and the point B after which the cohesive elements are totally broken.In the region between these two points, the cohesive elements are partially damaged.The damaged zone is in fact known as the process zone wherein the cohesive elements are in an intermediate state between the intact state (point A) and the fully damaged one (point B).In Figure 6 the distribution of damage value is plotted along the process zone at cycles 1500.In this process zone, the damage level varies from 0 (intact element) to 1 (fully damaged element).Location of the real crack tip is given by the cohesive element which is fully damaged (corresponding to the point B in Fig. 5).From Figure 6, it can be seen that the length of the process zone is approximately 0.35 mm.Note that the length of process zone depends essentially on the material properties, the geometry/size of the structure, and the loading mode.

Information at crack tip
Figure 7 shows the length of the process zone as function of the number of cycles, which corresponds to the evolution of the process zone during the crack propagation.In fact, the intact elements behave elastically whereas complete loss of traction occurs in the completely damaged elements.Thus, energy is dissipated only in the moving process zone as the crack front advances.At the beginning of loading, the process zone length increases up to 0.45 mm.Then, it drops suddenly, which corresponds to crack initiation (corresponding to the breaking of the first cohesive element at the end of the structure) at approximately 100 cycles.This is followed by a slight fluctuation of the process zone length until stabilization at 200 cycles.Finally the process zone length remains quite constant during the stable crack propagation regime and is about 0.35 mm.
Figure 8 shows stress-strain hysteresis loops for an element nearby the stress concentration region.These loops were obtained during the temperature cycling.It can be seen that the stress level increases with the temperature cycling and it tends to be stabilized.This stress-strain hysteresis behavior reflects also the behavior of the solder alloy modeled by Anand viscoplastic law.

Fatigue lifetime prediction
Figure 9 shows the crack evolution as function of the number of cycles.The crack initiation occurs at about 100 cycles, which corresponds to the inflection point on the crack evolution curve.Note that this number of cycles to initiation is consistent with the number estimated from the curve of the process zone size evolution (Fig. 7).Prior to the point of inflection, some cohesive elements were degraded to some extent but the crack was not yet formed.Once the crack is initiated, the crack evolution becomes faster.Then, its acceleration is observed between 120 and 600 cycles.Beyond 600 cycles, the crack propagation is   almost stable, which corresponds typically to a stable stage of crack propagation.
The lifetime of the system can be obtained by considering that the component fails when the crack extension covers all the solder width.In reality, the solder joint may be considered faulty long before the crack crosses the solder from end to end.In semi-empirical studies of solder reliability, such as in the Coffin-Manson law, the number of cycles to failure is defined with the 50% load drop criterion.Hence, the number of cycles at half of the initial load is taken practically as the number of cycles to failure.The system lifetime, thus defined, seems more relevant in the context of industrial specifications related to the reliability of electronic components.Therefore, according to this criterion, the computed lifetime is obtained when the crack extends to nearly half of the solder width (load carrying capacity decreased by half).Based on this failure criterion, the number of cycles to failure of the assembly under passive temperature cycling conditions is about 2000 cycles as obtained from Figure 9.It is noteworthy that the order of magnitude of this number of cycles is consistent with the literature about fatigue lifetime predictions for BGA (Ball Grid Array) solder joints under harsh thermal cycling conditions [3].The lifetime estimation, albeit approximate due to material data values collected from different sources, would logically become more accurate, provided that the materials are fully characterized through appropriate experimental tests.However, in 306-page 7  the present study, the objective was primarily to propose a solder joint reliability model incorporating physicallybased degradation mechanism.

Concluding remarks
In this work, a simplified IGBT module has been modeled for the prediction of solder reliability under passive temperature cycling condition.A CZM based fatigue model was used in conjunction with a viscoplastic behavior for the bulk solder.This approach gives the damage evolution in interface elements based on their tractionseparation law.The fatigue model presented in this paper incorporates features such as: (i) introduction of an energy density based damage variable for describing the damage state of the interface; (ii) treatment of the mixed mode situation based on the equivalent energy to fracture; (iii) consideration of temperature dependence of the material properties.
To estimate the fatigue lifetime of the solder joint, a failure criterion based on a critical value of the crack length is used.The number of cycles to failure was obtained when the crack length reaches half of the total structure length, which is considered as being a critical level at which the system becomes faulty.Although the numerical simulations have not been directly validated with experimental data obtained from real temperature cycling tests, the results appear to be in qualitative agreement with the data provided by many authors.
The present model offers good future development perspectives for modelling the failure problem often seen in the solder layers of power modules.For example, the extension of the cohesive zone approach to 3D should be done as the cracking mechanism in the solder layer is strongly influenced by the real 3D geometry of the die.Indeed, the crack in the case of 3D analysis is expected to initiate from the four corners of the die where stress concentrations exist, resulting in modified crack growth kinetics.

Fig. 2 .
Fig. 2. Steady strain rate as function of equivalent stress and temperature.

Figure 5
Figure 5 illustrates an in-plane shear stress field obtained after 1500 cycles.From the global view of the inplane shear stress field, one can easily identify the position of the crack tip where the concentration of shear stress occurs.The cohesive elements preceding the crack tip are

Fig. 4 .
Fig. 4. (a) Schematic representation of the applied passive thermal cycling, (b) the boundary conditions applied to the assembly, and (c) the mesh and a zoom in the region of the die-solder interface.

Fig. 5 .Fig. 6 .
Fig. 5. Example of in-plane shear stress field showing a stress concentration at the crack tip (obtained after 1500 cycles).

Fig. 7 .
Fig. 7. Evolution of the process zone length as function of the number of cycles.

ZFig. 8 .
Fig. 8. Stress-strain hysteresis loops for an element in the solder (modeled by Anand constitutive law) nearby the crack tip.

Fig. 9 .
Fig. 9. Crack length extension with respect to the number of cycles.