1. Introduction
Electronic devices appear almost everywhere, including telecommunications and the vehicle industry. These are technically printed circuit boards (PCB) on which electronic components are mounted by using solder joints. Solder joint failure is one of the most critical issue in electronics’ reliability [
1], since they are quite vulnerable. The miniaturization of electronic devices, the application of surface-mounted devices (SMD) and the development of increasingly dense integrated circuit (IC) packaging types, such as the ball grid array (BGA) and quad flat package (QFP), make the solder joint even more vulnerable [
2,
3]. SMD components in general are more critical than through-hole components when the lifetime is considered. The main loads on PCBs and solder joints include mechanical shocks [
4], mechanical vibration [
3,
5], thermal shocks [
6] and cyclic thermal expansion-induced stress [
1,
6]. In this study, thermal cycling is the focus. A change in the temperature causes thermal stress in the microelectronic assembly due to the mismatch in the coefficient of thermal expansion (CTE) between the materials, usually denoted by
, and the dimension is 1/K.
To improve the reliability of electronic equipment subjected to cyclic thermal loads, methods to predict the expectable lifetimes of solder joints have been intensively investigated and developed in the last few decades, but there are still open questions. Physical thermal cycling tests require a long time, even though these tests are usually accelerated [
7,
8]. Besides physical testing [
8,
9], which uses actual temperature cycling profiles [
1,
8], simulation-based lifetime prediction methodologies have been continuously under development [
10,
11]. There is difficulty in estimating the lifetime as the variations in the package assembly and component vendor variations strongly affect the true lifetime [
12].
The introduction of lead-free solder alloys introduced challenges due to the lack of experience related to manufacturing, lifetime properties [
9,
13,
14] and mechanical properties [
15]. After many studies, in the industry, the Sn-Ag-Cu (SAC) alloy material is the most widely used among lead-free alloys [
16]. Regarding the thermal–fatigue properties, SAC solders are sensitive to the thermal cycle profile to a greater extent than Sn-Pb alloys [
13]. Nevertheless, the SAC solder material is widely used [
17], even with different additives [
2]. It is emphasized in [
18] that slight changes in the solder alloy composition and microstructure can result in large changes in the temperature-dependent properties.
It has been found that the microstructure of a lead-free solder plays a crucial role, and it is highly sensitive, not only to the alloy applied but to the conditions of molten solder solidification [
19,
20,
21,
22,
23]. In fact, intermetallic compounds (IMCs) can be formed at the interface between the solder and the copper substrate [
24]. The IMC grain size has a significant effect on the lifetime as solder joints are becoming smaller. Results regarding microstructural crack propagation are summarized in [
13], where it is concluded that Sn-Pb and SAC have very different behaviors: SAC recrystallizes in the plastic zone under temperature cycling, and this finer crystal structure provides an easy path for cracks. The assumption of a homogeneous material leads to errors in lifetime estimation, since the microstructural effects become stronger with miniaturization [
24]. The recent paper [
25] reviewed the thermal fatigue failure mechanisms of solder joints, including microstructural changes during thermal cycling, factors affecting the fatigue life and simulation methods for thermal fatigue life prediction. It is pointed out based on the reviewed literature that the excessive growth of the IMC layer contributes greatly to solder fatigue by allowing crack growth along the IMC and the solder. Still, most of the solder joint fatigue life models assume a homogeneous distribution of the material properties throughout the solder volume.
The most reliable and widely used fatigue models since the 1960s are reviewed in [
26,
27]. The authors of [
27] classify the fatigue life failure models into the following four fundamental categories:
plastic strain-based,
creep strain-based,
energy-based and
damage accumulation-based models. There is also a distinct group for other types of models. Stress-based models are also reviewed in [
26]. It is important to note that some of the reviewed fatigue models are not particularly related to soldering. The reviewed fatigue models are generally applicable for low cycle fatigue (LCF), which is also the focus in the present paper. In the LCF of solder joints, typically, temperature variation cause plastic deformation and ultimately lead to failure.
1.1. Methodology of Virtual Thermal Fatigue Lifetime Estimation
In the case of the FEM-based virtual estimation of thermal fatigue lifetimes, knowledge of the geometry and material behavior is a key aspect, such as the elastic and creep parameters and the thermal expansion of the different materials in the assembly [
28,
29,
30]. The determination of the accurate values of the material parameters for the FEM simulation is often not possible from measurements, as [
31] shows. Since the material parameters deviate, the outputs from the FEM simulation will also be affected by this deviation. With the output values from the FEM simulation, which could be the different stress, strain or energy values calculated from the components, the fatigue models predict the lifetimes of the investigated components [
26,
27]. However, there is a very large deviation between the predicted lifetimes according to different fatigue models; moreover, the tuning of these fatigue models is not prescribed and can deviate in the literature [
27,
28,
29,
30].
1.2. Material Models in Fatigue Life Estimation
Lead-free solder alloys have large variety, which necessitates the assessment of the material parameters involved in the constitutive model [
13,
19]. The constitutive model has a large effect on the simulation-based prediction of the solder lifetime [
32]. Solder alloys are generally modeled by elastic, elastoplastic and viscoplastic models [
27]. Viscoplastic models are capable of modeling the creep phenomena. It is not completely clear which model should be chosen in a certain situation; therefore, they are benchmarked in many papers. Furthermore, the material parameters are under research [
31].
Sn-Pb and SAC alloys are compared in [
6] using FE simulations. It has been reported that crack growth in the solder joints of chip resistors is faster for SAC alloys than for Sn-Pb alloys. Interestingly, reference [
25] concludes that the crack growth is faster in lead-containing solders than lead-free solders such as SAC alloys. The Sn-Bi and SAC305 materials were compared in [
7] by means of accelerated thermal cycling tests, which allowed the parameter identification of the Sn-Bi constitutive model.
Viscoplastic models are capable of modeling the creep behavior of the material, and several variants can be found in the literature, referred as the Garofalo model (also referred to as the hyperbolic sine or Garofalo–Arrhenius model) [
6,
10,
11,
15,
28,
32,
33,
34], the Wiese model or double-power model [
34] and the Anand model [
35,
36]. The Garofalo, Wiese and Anand models are benchmarked in [
16] for different variations of SAC alloys. The Garofalo and hyperbolic sine creep constitutive models are benchmarked to simulate damage in the SAC solder joints in flip chip and resistor components in [
32]. The accumulated creep strain and strain energy density are generated from simulations using various models and then compared in terms of the solder fatigue life. The creep parameters, damage index, creep model and solder geometry indeed have significant effects on the lifetime. A similar study in [
11] compares the effects of the creep parameters. The mentioned paper suggests a specific constitutive model based on the damage indication (creep strain or energy density) and component type (resistor or BGA flip chip).
The Anand material model [
35,
36,
37] was developed in response to the need related to the high-temperature plastic forming of metal products. In addition to the temperature dependence of the material behavior, the model takes into account strain hardening, recrystallization and the strain rate. The model has been made coherent with solder joint alloys and geometry by conducting creep tests and constant strain rate tests in [
38]. Since then, the Anand model has been widely used, e.g., in [
1,
7,
12,
16,
17,
29,
30,
31,
39]. Two material parameter sets were introduced for the Anand model for the SAC material in. The Anand model parameters are assessed/collected in [
16,
29,
30,
31].
1.3. Fatigue Models
Fatigue models can be divided into two groups. Some analytical models, such as the Engelmaier model [
8] and the Blattau model and its variations [
11,
18,
40,
41], assume a rectangular-shaped solder geometry and require geometric data, the temperature range and some material properties as input data. Other fatigue models use inputs from FE simulations, where the strains, energies and other types of inputs are calculated numerically [
10,
29,
30]. In this research, FE simulation-aided fatigue models are the focus. In
plastic strain-based models, the plastic strain is assumed as a measurable factor in the predicted fatigue life [
26,
27]. This means that any finite element (FE) analysis is required to calculate the plastic strain. The accuracy of the calculated plastic strain might depend on the FE model, including the elastoplastic material model. Creep strain is typically not considered, but the thermal cycle acceleration factor can be considered. Plastic strain-based models include those in [
26,
27].
Creep damage-based models are generally founded on viscoplastic material models and suppose that deformations develop in time, even when the load is time-invariant and smaller or larger than the yield strength [
26,
27]. For steel, creep appears at high temperatures. However, for some materials, such as lead-free solder materials, creep occurs even at room temperature due to the material’s highly homologous temperature [
42]. Three stages are observable for most solder materials. In the secondary stage, the deformation grows proportionally to the elapsed time. It is therefore called steady-state creep. Most of the creep models assume steady-state creep, while the first and third stages are omitted due to computational challenges caused by nonlinearities. High-temperature creep and fatigue damage are handled separately or combined in the models. Ultimately, the creep–fatigue interaction is considered. Models that combine both creep effects and fatigue effects are presented in [
26,
27].
Energy-based models rely on the non-zero-area stress–strain loop in cycling tests [
26,
27]. The hysteresis energy is absorbed due to plastic deformation and the absorbed energy accumulates until fracture or crack. The absorbed energy is therefore a damage indicator. Several models have been developed in response to the need to better calculate the area within the stress–strain loop and to more accurately link the energy with the fatigue life. Energy is related to both stress and strain and seems more accurate than plastic strain-based or creep-based models. On the other hand, energy-based models are capable of only calculating crack initiation and cannot predict the crack propagation inside the solder joint. Energy-based models include those in [
26,
27].
Damage accumulation-based models, also known as
fatigue cumulative damage theories, assume that fatigue damage accumulates in the material based on two different concepts: damage accumulation and crack propagation. Such models can be found in the literature [
26,
27]. Linear damage accumulation theories assume that damage is independent of the load sequence and the ratio of damage accumulation is independent of the stress level. Nonlinear theories have the ability to unify the damage originating in various types of load by capturing their interplay. Crack initiation and propagation approaches are fundamentally damage tolerance methods based on fracture mechanics. Damage accumulation models based on cumulative damage are given in [
26,
27].
The authors of [
2,
25] also consider these fatigue models, but with no additional information. However, several fatigue life models are available, and [
29] concludes that the predictions are strongly inconsistent regarding the predicted fatigue life, even in the case of the very careful selection of the model in an actual situation. This conclusion was achieved after carrying out FE temperature cycling simulations and then benchmarking different fatigue models: Coffin–Manson, Engelmaier, Solomon and Syed. It is also important to note that most of the models are developed for uniaxial loads; therefore, their application to multi-axial loads [
43] raises some questions and might be a crucial reason for inaccurately estimated lifetimes.
1.4. Solder Geometry
The focus of scientific studies is mostly the BGA solder type [
7,
8,
9,
13,
16,
23,
28,
29,
44]. The reason is that it is widely used in the industry, and its modeling is relatively simple compared to that of more complicated solder geometries in other SMD devices, such as gull-wing packages or SMD capacitors and resistors with three or five wetted sides. SMD capacitors with five wetted sides are the focus of this study.
The manufacturing process generally affects the geometry of the solder material, and the solder shape has a direct relationship with the lifetime [
9,
30,
32,
45,
46]. For instance, hourglass- and barrel-shaped BGA joints are compared in [
9].
Note that different component types require the re-tuning of the lifetime estimation model. A model tuned for BGA is not directly applicable to a resistor or capacitor component. For instance, the Syed lifetime model [
47,
48] is applied in [
11] for resistors with proper tuning.
Cyclic thermal load-induced fatigue phenomena can be analyzed only via sophisticated numerical simulations, as [
18] states, in which semi-analytical fatigue models are developed based on the strain energy density for light-emitting diode (LED) components soldered on metal core PCBs. Design guideline instructions related to solder geometries are formulated in [
16]. The optimal solder shape is addressed in [
30,
45] for a resistor component: concave, straight and convex shapes lead to different lifetimes. Similar results are achieved for SMD capacitors in [
49]. It is reported in [
6] that the estimated lifetime clearly grows with a decreasing size in SMD resistor components mounted with Sn-Pb alloys, and a less clear but similar tendency is reported for SAC solder alloys. Indeed, the applied solder volume directly affects the solder shape and therefore the lifetime of the resistor [
17,
46] and capacitor [
50] in solder joints. The optimal solder volume is calculated in [
17] for resistor components. In [
24], tensile strength results for different ratios of the BGA solder volume to substrate area are reviewed. Another important factor that affects the lifetime is the standoff height [
51] of the SMD component. The effect of the standoff height was assessed in [
44] for BGA joints, and it was found that reduced stresses developed in response to temperature loads in solder joints with elevated standoff heights, referred to as a column grid array. In [
52], the effect of the standoff height on the fatigue lifetime was analyzed in parallel to the benchmarking of the Coffin–Manson model, the Shi frequency-modified Coffin–Manson model, the Engelmaier model and the Solomon fatigue model. Similarly to the effect of the real solder shape, the geometry of the solder joint in the simulation model also plays a crucial role in the predicted lifetime. The modeling of the solder geometry is therefore an essential aspect when solder joint fatigue theories are applied in industrial complexity problems. Manually generated geometries are used in many research works [
6,
7,
10,
11,
28,
29,
32,
38,
39]. However, simple models from computer-aided design (CAD) software tools are effective in industrial applications regarding human effort, although the achievable accuracy of the fatigue life estimation when using them is questionable. Another concept is to create the solder geometry by means of physics-based simulations, such as the widely used
Surface Evolver computational tool [
53] applied in [
17,
30,
45,
49,
50,
51,
52,
54]. Indeed, the FE results and therefore the fatigue life estimation are more accurate in general compared to manually created models [
51,
52,
55]. On the other hand, the creation of such models is time-consuming when the industrial complexity is considered, since, prior to the FE modeling, each different electronic component requires a simulation only for geometry creation. Related the solder geometry issues, References [
50,
54] investigate the phenomenon of component self-alignment during reflow soldering and solder material solidification. Information related to the effect of misalignment [
51] on the lifetime can be found in [
17,
49].
In solder joints, the strains and stresses are indeed multi-axial due to the highly complex geometry. The challenges in relation to fatigue caused by multi-axial stresses are the focus of [
43]. Although the presented application is not related to soldering and is related to high cycle fatigue, the issue of multi-axial fatigue must be kept in mind.
The authors in [
6] estimated the lifetime based on the creep strain energy density, which served as input for Miner’s law. The concept of the elastoplastic strain energy density (SED) is applied in [
56] for fatigue life estimation in the presence of surface imperfections on structural components. The results indicate that the elastoplastic model provides considerably higher accuracy than the linear elastic models. Meanwhile, [
56] uses the idea from [
57], according to which fatigue failure theoretically occurs when a critical value is reached by the SED averaged over the control volume. The idea of a properly chosen control volume is widely used in FE-based solder joint fatigue analysis too—for instance, in [
1,
6,
17,
29,
30,
45,
46,
58]. However, the choice of the control volume is not obvious.
1.5. Hypotheses and Goals
This paper benchmarks ten fatigue lifetime models on two types of SMD capacitor components, C0603m and C1005m, with four solder geometry modeling approaches and three different standoff heights. The first main goal of the present paper is to benchmark a variety of solder geometry modeling methods. We hypothesize that
(1) the FE results from simplified, manually created geometries offer good approximations of the results obtained from a physics-based geometry. The physics-based geometries, which are used as a reference, are generated by the widely used numerical software called Surface Evolver [
53]. Secondly, this work aims at the comparison of the available fatigue models, which are applicable in the post-processing phase of solder joint FE analysis. Different fatigue models require different inputs and rely on different material models. An elastoplastic material model and two viscoplastic models are included in this study: the Garofalo and the Anand viscoplastic models. We hypothesize
(2) that the more advanced the material model is, the more robust the results are. Hypothesis
(3) is that by reducing the SMD component size, the predicted lifetime is improved. We also hypothesize that
(4) increasing the solder volume increases the predicted lifetime. Our last hypothesis is that
(5) increasing the standoff height also improves the predicted lifetime.
2. Materials and Methods
Ten different fatigue models are applied in this paper to the outputs from the FE simulation. Each geometrical case considers three different materials defined for the solder joint part of the assembly: • elastoplastic material model; • Garofalo viscoplastic material model; • Anand viscoplastic mateiral model. Four of the fatigue models are based on the results from the elastoplastic simulation: • Coffin–Manson, • Shi frequency-modified Coffin–Manson, • Solomon–Tolksdorf, • Morrow model. Five of the fatigue models are based on the results from the viscoplastic simulation in the case of the Garofalo and Anand material models: • Syed-I, • Darveaux, • Akay-II, • Joseph–Jerries, • Gustaffson model. One of the fatigue models is based on the results from both the elastoplastic material and the viscoplastic material models: • the Pan model. The inputs of the fatigue models are either the accumulated strain outputs or the strain energy density outputs, which are calculated for the geometry set from the solder joint’s lower portion using the volume-weighted average function of the values. The investigation is carried out for industrially used SMD capacitor components of two packaging types: C0603m and C1005m, from the manufacturer TDK (Suzhou) Co., Ltd., China.
Two assumptions are made in the evaluation process: (1) the Surface Evolver geometry is the closest to the reality; therefore, it provides the most accurate result; (2) the Syed fatigue life model is the most reliable for resistor and capacitor components, since the reliable parameters of other models are not validated for resistor-type components. Therefore, the Surface Evolver geometry, together with the Garofalo material model and Syed lifetime model, is considered as a reference.
2.1. Material Models
2.1.1. Elastoplastic Material Model
The elastoplastic model was defined to be linear–elastic below the strain value of 0.01. In this region, the elastic modulus was set to
GPa. Over the strain value of 0.01, the behavior of the material was plastic, defined as reaching the strain of 0.03, at which the stress in the material is 50 MPa. The stress–strain characteristics are taken from [
59].
2.1.2. Garofalo Viscoplastic Material Model
The Garofalo viscoplastic material model was introduced in [
33] and is still used, e.g., in [
6,
10,
11,
28,
32]. The Garofalo model gives the relation between the strain rate and the stress in the following way:
where
is the strain rate,
A is the material constant,
is the hyperbolic law multiplier,
is the stress,
n is the stress exponent,
Q is the activation energy,
R is the universal gas constant, and
T is the absolute temperature. The Garofalo model parameters are taken from [
28] and collected in
Table 1.
2.1.3. Anand Viscoplastic Material Model
The Anand viscoplastic material model can be described by the following equations:
where
is the viscoplastic strain rate,
A is the pre-exponential factor,
is the applied stress,
s is the deformation resistance,
m is the strain rate sensitivity of stress,
Q is the activation energy,
k is the Boltzmann constant, and
T is the absolute temperature. The evolution of the deformation resistance is governed by
where
is the hardening/softening constant,
a is the strain rate sensitivity of hardening,
is the saturation value of the deformation resistance, and
is the coefficient for deformation resistance saturation. The parameters in
Table 2 are taken from [
30].
2.2. Fatigue Lifetime Estimation Models
The Syed lifetime model [
47,
48] is proven to be applicable for resistors in [
11]. The Syed model was tuned for the Surface Evolver geometry in such a way that the predicted lifetime was N(50%) = 10,000. Note that there is only one free parameter in the Syed model, but there are generally two parameters in other lifetime models. Therefore, the other lifetime models were tuned in such a way that the load obtained from the FE resulted in N(50%) = 10,000. Furthermore, the tuning also guaranteed that the lifetime at a double load was N(50%) = 5000, similarly to the Syed model. This consideration is physically not necessarily reasonable, but a basis for comparison must be chosen. In the further subsections, the lifetime models and their parameters are presented.
2.2.1. Coffin–Manson Fatigue Model
The Coffin–Manson model estimates the lifetime with the following equation:
where
m is the fatigue exponent and
C is the ductility coefficient. The chosen parameters for the Coffin–Manson model are
,
.
2.2.2. Frequency-Modified Coffin–Manson Model by Shi
The frequency-modified Coffin–Manson model by Shi estimates the lifetime with the following equation:
where
v is the frequency and
k is the frequency exponent. The chosen parameters for the Shi model are
,
,
,
.
2.2.3. Morrow Fatigue Model
The Morrow model estimates the lifetime with the following equation:
where
m is the fatigue exponent,
C is the material ductility coefficient, and
is the plastic strain energy in a steady-state loop. The chosen parameters for the Morrow model are
,
.
2.2.4. Solomon–Tolksdorf Fatigue Model
The model by Solomon and Tolksdorf uses the following equation for fatigue lifetime estimation:
where
v is the frequency,
k is the frequency exponent determined from the connection of the fatigue life and frequency,
n is the frequency exponent determined from the relationship between the strain energy density and the frequency,
is the strain energy density, and
C and
m are the same as in the previous Morrow fatigue model. The chosen parameters for the Solomon–Tolksdorf model are
,
,
,
,
.
2.2.5. Pan Fatigue Model
The Pan model uses the following equation for fatigue lifetime estimation:
where
C is the critical strain energy density, which is geometry-dependent.
is the number of cycles to failure,
is the time-dependent plastic energy,
is the time-dependent creep energy, and
a and
b are the weighting factors of the corresponding energies. The chosen parameters for the Pan fatigue model are
,
,
.
2.2.6. Joseph and Jerries Fatigue Model
The Joseph and Jerries model estimates the lifetime with the following equation:
where
is the intrinsic material property, and
W is the energy dissipated in one cycle. The chosen parameter for the Joseph and Jerries fatigue model is
.
2.2.7. Syed-I
The Syed lifetime model [
47,
48] is proven to be applicable for resistors in [
11]. The Syed-I model estimates the lifetime with the following equation:
where
is the accumulated creep strain for the whole cycle, and
C is the creep ductility for the creep mechanism. The chosen parameter for the Syed-I fatigue model is
.
2.2.8. Akay-II Fatigue Model
The Akay-II model uses the following equation for fatigue lifetime estimation:
where
is the total strain energy calculated form the volume-weighted average over a cycle in a steady-state region;
and
k are load-independent material constants. The chosen parameters for the Akay-II fatigue model are
and
.
2.2.9. Darveaux Fatigue Model
The Darveaux model was originally introduced in [
60,
61] and is still used by many researchers [
1,
12,
14,
16,
17,
30,
39,
62,
63]. The model uses the following equation for fatigue lifetime estimation:
where
is the characteristic life for the solder joint and
a is the final crack length.
is the fatigue life at the initiation of the crack.
The equation for the calculation of the cycles for crack initiation is the following:
The equation for the calculation of crack growth is
where
,
,
, and
are life correlation constants. The chosen parameters for the Darveaux fatigue model are
,
,
, and
, with
in the case of the C0603m capacitor and
in the case of the C1005m capacitor.
2.2.10. Gustafsson Fatigue Model
The Gustafsson model estimates the lifetime with the following equation:
where
is the characteristic life of the solder joint and
a is the total possible crack length. The Gustafsson model is a modification of the Darveaux model. This model views both primary and secondary cracks with initiation as growing towards each other. The primary crack is indicated with the subscript
p and the secondary crack is indicated with the subscript
s.
is the secondary crack energy-based term and
is the primary crack energy-based term.
is the secondary crack propagation and
is the primary crack propagation. The chosen parameters for the Gustafsson model are the same as in the case of the Darveaux model, so the primary and secondary crack initiation and crack propagation rates are the same in both cases:
,
,
,
, with
m in the case of the C0603m capacitor and
m in the case of the C1005m capacitor.
2.3. FE Models with a Variety of Geometries
A variety of geometry models are compared. The first concept in
Figure 1a is a simple extrusion geometry on the basis of a copper pad. This geometry can be quickly and easily generated by engineers. More realistic solder geometries can be achieved by using the ’loft’ feature of any CAD software. Here, the rectangular surface of the copper pad is connected to another rectangle placed in a tilted position inside the component volume, as
Figure 1b shows. A further enhanced modeling concept is based on the loft connection of ellipse shapes, as presented in
Figure 1c. Although the latter two geometries are not very complicated to create in industrial circumstances, they are closer to the physics-based reference solder shape shown in
Figure 1d. The meshed models for the four different geometric input methods can be seen in
Figure A1 The Surface Evolver numerical solver was used [
53] to generate the reference geometry.
The geometric data of the two capacitors are collected in
Table A3. The solder volumes and the lower solder volume values are, respectively, collected in
Table A4 and in
Table A5. The Young’s modulus, Poisson’s ratio and thermal expansion coefficient values for each material are presented in
Table A6. The meshing of the assembly is illustrated in
Figure 2. For the solder joint and termination parts, a tetrahedral mesh was used; for the other parts in the assembly, a hexahedral mesh was used. The mesh sizes were different for the C0603m and C1005m capacitor cases for the assembly parts listed in
Table 3. Mesh verification tests were carried out for the mesh sizes listed in
Table A1 and
Table A2. The meshed models during the verification test are given in
Figure A2. In the FE analysis, the “static general” for the elastoplastic cases and the “visco general” solver from Abaqus were utilized, allowing thermal expansion over time. The strains and stresses were induced by the thermal expansion differences in the different materials.
2.4. Definition of Loads
A temperature cycling profile similar to those in [
1,
30] was applied. One time period of the temperature cycle was 600 s. The temperature cycle was symmetric, so the ramp-up time and ramp-down time were both 100 s, and the ramp rate was constant. The idle time at high temperatures and at low temperatures was defined as 200 s. The low temperature was 223 °K and the high temperature was 393 °K. The starting temperature of the test was set to 291 °K. It is important to mention that the chosen temperature shock cycle profile period was shorter, similarly as in [
30], while the JEDEC standard for temperature cycles [
64] recommends longer period times.
Because of the CAD models’ symmetry to their center points on both the X plane and Y plane, the FE models can be reduced to quarters when adding the symmetry constraints. Moreover, the quarter models are fixed in the vertical position at the intersection of the two symmetry planes at the bottom point of the PCB.
3. Results and Discussion
In this section, the FE simulation results and predicted lifetimes are presented. All of the different geometry cases used the symmetry constraints in order to reduce the computational effort.
Figure 2 shows the reduced meshed quarter assembly model on the left side and the results from the FE simulation with the creep energy density variable on the right side. Due to the symmetry constraints, as the left side of
Figure 2 shows, the meshed assembly is cut to a quarter, while the solder joint geometry is cut to a half.
The outputs from the FE simulation in the elastoplastic material model case for the solder joint were the equivalent plastic strain, referred to as PEEQ, and the energy dissipated by rate-independent and rate-dependent plasticity, per unit volume, referred to as PENER. In the case of the Anand and Garofalo viscoplastic material models, the two outputs were the equivalent creep strain, referred to as CEEQ, and the energy dissipated by creep, swelling and viscoelasticity, per unit volume, referred to as CENER.
To evaluate the FE simulation outputs, the volume-weighted average was used regarding the volume of the lower portion of the solder joint for every output. The volume-weighted average of the outputs from the FE simulations was the input for the 10 fatigue models, which gave the resulting lifetime predictions.
3.1. Standoff Height and Solder Volume
In
Figure 3, the estimated lifetimes from the 10 different lifetime estimation methods are presented in the case of the elastoplastic and Anand viscoplastic material models (see
Table A8 and
Table A10 in
Appendix B). The bar plot groups with 3, 6, and 9
m standoff heights are from the simulation of the C0603m capacitor, while the bar plot groups with 5, 10, and 15
m standoff heights belong to the simulation results of the larger-sized C1005m capacitor.
Figure 3a shows the lifetimes calculated from the geometry input of extruded rectangle-shaped solder joints (EXT),
Figure 3b shows the solder joint shaped from the rectangular loft (REC),
Figure 3c shows the ellipse loft-shaped solder joint (ELL), and
Figure 3d shows the resulting lifetimes for the Surface Evolver-generated solder geometry (SE).
Four out of the the ten fatigue models used elastoplastic material models: the Coffin–Manson, frequency-modified Coffin–Manson (labeled as Shi Freqmod-CM), Morrow, and Solomon–Tolksdorf fatigue models. Meanwhile, the other six fatigue models relied on the Anand viscoplastic material model: the Pan, Joseph and Jerries, Syed-I, Akay-II, Darveaux and Gustafsson fatigue models. These six fatigue models were applied also for the Garofalo viscoplastic material model. All of the lifetime estimation methods were tuned to predict 10,000 cycles as the mean lifetime in the case of the Surface Evolver geometry in the elastoplastic and viscoplastic Anand material models with a standoff height value of 6
m, as the barplot in
Figure 3d shows.
Figure 4 shows the estimated lifetimes in the case of the Garofalo material model in the same bar groups, categorized by the standoff height value of the assembly, as in
Figure 3, in each case of the geometry inputs.
Figure 4a is for the EXT,
Figure 4b for the REC,
Figure 4c for the ELL and
Figure 4d for the SE geometry input cases. In this case, only six of the ten fatigue models could be applied: the Pan, Joseph and Jerries, Syed-I, Akay-II, Darveaux and Gustafsson fatigue models. The related lifetime data are presented in
Table A12 in
Appendix B.
In
Figure 5, the PEEQ, PENER, CEEQ and CENER output variables can be seen as a function of the component standoff height. The points marked with ‘
x’ are the results for the C0603m capacitor, while the points marked with ‘
o’ are the results for the case of the C1005m capacitor. The different colors indicate the geometry input methods. The points for the same geometry inputs are connected with dashed lines to improve the visualization of the tendencies. In
Figure 5a,b, the plastic material case outputs are presented, while, in
Figure 5c,d, the outputs from the Anand viscoplastic material case are presented. The related FE output data are presented in
Table A7 and
Table A9.
Figure 6 presents the PEEQ, PENER, CEEQ and CENER output data in the same order as in
Figure 5, but with the difference that the data are plotted as a function of the solder volume in the lower portion of the capacitor (for the lower volume data, see
Table A5). In
Figure 6a,b, the plastic material case outputs are presented, while, in
Figure 6c,d, the outputs for the Anand viscoplastic material case are presented. The related FE output data are presented again in
Table A7 and
Table A9.
Figure 7 shows the CEEQ and CENER outputs from the Garofalo viscoplastic material models; in
Figure 7a,b, they are shown as a function of the component standoff height, while, in
Figure 7c,d, they are given as a function of the solder volume in the lower portion. The related FE output data are collected in
Table A11 and the lower portion volumes are given in
Table A5.
3.2. Geometry Input Type
In
Figure 8, the lifetime estimation bars are grouped according to the geometry input methods in the case of the elastoplastic and viscoplastic Anand material models. Each subplot in
Figure 8 only shows the predicted lifetimes in the case of the given standoff height. For the numerical data, see
Table A8 and
Table A10.
Figure 8a shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C0603m capacitor;
Figure 8b shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C0603m capacitor;
Figure 8c shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C1005m capacitor; and
Figure 8d shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C1005m capacitor.
In
Figure 9, the lifetimes are grouped in the same way as in
Figure 8, but, in this case, the lifetimes were predicted from the outputs of the simulation results, where the Garofalo viscoplastic material model was applied. In this case, only the six fatigue models based on the viscoplastic outputs are plotted. For the numerical data, see
Table A12.
Figure 9a shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C0603m capacitor;
Figure 9b shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C0603m capacitor;
Figure 9c shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C1005m capacitor; and
Figure 9d shows the predicted lifetimes in the case of different geometry input methods at a
m component standoff height in the case of the C1005m capacitor.
In
Figure 10, the FE outputs are plotted as a function of the lower portion solder volume (see
Table A5), in the same way as in
Figure 6. The numerical values are shown in
Table A7 and
Table A9. However, in this case, the color grouping and the marker types are intended to highlight the cases where the standoff heights have the same value in the model and the capacitor is the same size; only the geometry input method is different. For every four points in a group, a fitted line is plotted to visualize the tendency of the outputs as a function of the solder volume in the lower portion of the solder joint.
Figure 10a shows the tendencies in the case of the PEEQ output, and
Figure 10b shows the tendencies for the PENER output, when the elastoplastic material model was applied.
Figure 10c shows the tendencies in the case of the CEEQ output, and
Figure 10d shows the tendencies for the CENER output, when the Anand viscoplastic material model was used.
Figure 11 shows the tendencies in the same order as in
Figure 10c,d, but in the case of the Garofalo viscoplastic material model:
Figure 11a shows the tendencies in the case of the CEEQ output, and
Figure 11b shows the tendencies in the case of the CENER output.
4. Discussion
The predicted lifetimes in the plots in
Figure 3 and
Figure 4 show the same tendency, ranging from
m to
m for the C0603m capacitor and from
m to
m for the C1005m case. In particular, with an increase in the component standoff height, the predicted lifetime will also increase in 9 of the 10 fatigue models; the Akay-II fatigue model predicts the lifetime with an inverse tendency. In the case of
Figure 8 and
Figure 9, this tendency can also be seen.
In the case of the Syed model, in the case of the Anand viscoplastic model for m with the Surface Evolver geometry input, the predicted lifetime is 6309 cycles; for m with the same geometry input and same material model, the predicted lifetime is 10,000 cycles; and with the same parameters but with m, the predicted lifetime is 13,206 cycles. However, in the case of the Akay-II model, with the same material model and geometry input method, for m, it is 11,827 cycles; for m, it is 10,000 cycles; and, for m, the predicted lifetime is 9035 cycles.
This could be due to the equation of the Akay-II fatigue model, where the input is the total creep strain energy calculated for the investigated volume set. This means that, in this case, with an increase in the standoff height, the investigated volume will also increase. If the volume-weighted average strain energy does not decrease at a higher rate, this volume set increases the rate, and the predicted lifetime will decrease due to the higher input value. This prediction is feasible for passive SMD components, particularly when comparing it to the experimental results in the literature [
30].
Investigating the outputs from
Figure 5 and
Figure 7a,b, in all four cases of output variables, in all four cases of geometry inputs, and for all three different material models applied, the tendency indicates a decrease in the volume-weighted average output variables with an increase in the standoff height for one type of capacitor. For the elastoplastic case, regarding the PENER values in the case of the C1005m capacitor and in the case of the REC geometry input, shown by the green ‘
o’ marker in
Figure 5b, the following results are found: for a
m component standoff height, the PENER value is 0.0812; for
m, it is 0.0358; and, for
m, the PENER value is 0.0257.
Regarding the outputs as a function of the solder volume in the lower portion, as
Figure 6 and
Figure 7c,d show, the tendency is the same. For the Garofalo viscoplastic case, the CENER values in the case of the C0603m capacitor and in the case of the ELL geometry input, shown by the blue ‘
x’ in
Figure 7d, are as follows: when the lower portion volume of
and the standoff height is
m, the CENER value is 0.7418; when the lower portion volume is
, the CENER value is 0.2634. It is remarkable that not only does the solder volume increase, but the standoff height also increases.
From
Figure 5,
Figure 6 and
Figure 7, it can be seen that the outputs are larger for the C1005m capacitor than the C0603m capacitor. Overall,
Figure 3 and
Figure 4 show that the lifetime is smaller in the case of a larger capacitor; in 8 of the 10 fatigue model cases, the two exceptions are the Darveaux and the Gustafsson models. The explanation is that, in both fatigue model equations, the total crack length is required as a parameter. If the capacitor size is larger than the total crack length, the
a variable is also larger, which causes an increase in the lifetime.
In the comparison of the geometry input methods, in
Figure 10 and
Figure 11, the fitted lines with higher standoff height cases (green and blue dashed lines) show a decreasing tendency for the outputs, if the standoff height is fixed but the solder joint volume in the lower portion increases. For C1005m, with a
m standoff height and with the Anand viscoplastic model applied, as shown in
Figure 10c with the blue ‘
o’ marks, the following results are found: in the case of the ELL geometry input, the CENER output variable is 0.2010 with the lower portion volume of
, while, in the case of the EXT geometry input, the CENER output variable is 0.1799 with the solder volume in the lower portion of
.
In
Figure 10 and
Figure 11a,b, depicting the fitted lines for the points, in the case of the
m and
m standoff heights for C0603m, and in the case of the
m and
m standoff heights for the C1005m capacitor, the approximated line is fit to the corresponding points of the output variables calculated from different geometry inputs. • In the case of CENER with the Anand viscoplastic material model, for the C1005m capacitor, the corresponding points shown in
Figure 10d with blue ‘
o’ markers indicate the following. • In the case of the EXT geometry input method, when the lower portion volume of the solder is 0.001500 and the CENER variable is 0.1799, the fitted linear CENER value here is 0.1721. The error of the fitting at this point is
. • In the case of the REC geometry input method, when the lower portion volume of the solder is 0.001404 and the CENER variable is 0.1776, the fitted linear CENER value here is 0.1791. The error of the fitting at this point is
. • In the case of the ELL geometry input method, when the lower portion volume of the solder is 0.001143 and the CENER variable is 0.2010, the fitted linear CENER value here is 0.1980. The error of the fitting at this point is
. • In the case of the SE geometry input method, when the lower portion volume of the solder is 0.001400 and the CENER variable is 0.1701, the fitted linear CENER value here is 0.1794. The error of the fitting at this point is
.
However, when the component standoff height is
m in the case of C1005m or
m in the case of C0603m, the errors are much larger, as
Figure 10 and
Figure 11 show.
From the simulation results for all four different geometry input methods and for all three material models, an increase in the solder volume and standoff height has a positive effect on the predicted lifetime of the solder joint in the case of thermocyclic loading, as can be seen in
Figure 5,
Figure 6 and
Figure 7. The results in
Figure 10 and
Figure 11 show that, in the case of higher component standoff heights, the sensitivity to the geometry input method type is lower. Moreover, if the component standoff height is small, the sensitivity to the geometry input method type is higher; if the component standoff height and the volume are properly set, this gives the opportunity to use a substitute geometry, instead of creating a geometry with Surface Evolver, which requires more effort. The larger-sized capacitor C1005m has slightly higher values for the volume-weighted output variables as calculated from the corresponding volume set than the C0603m capacitor, as shown in
Figure 5,
Figure 6 and
Figure 7. In contrast, the predicted lifetime is smaller in the case of a smaller-sized capacitor, as
Figure 3 and
Figure 4 show.
In the studied cases, the stresses were approximately multi-axial and proportional [
43], which simplified the evaluation process: none of the 10 fatigue life models was developed considering multi-axial non-proportional stresses. On the contrary, they are derived for uni-axial cases. In further research, the possibility of multi-axial non-proportional fatigue models will be addressed.
5. Conclusions
In this study, two different types of capacitors were investigated, with six different component standoff heights; in each case, four different approaches to the geometry modeling of the solder joint were applied for FE analysis. The investigated outputs in the case of the plastic model for the solder were the accumulated plastic strain and the plastic strain energy density, while, in the case of the viscoplastic material model, the accumulated creep strain and the creep strain energy density were calculated. Three different material models were defined for the solder in all geometry cases. The outputs of the FE analysis were compared to each other as a function of the solder joint volume, component standoff height and geometry input method. Moreover, 10 different lifetime estimation methods were applied to the FE outputs to determine the impacts of the different cases on the lifetime.
The results from the benchmark analysis show that the solder volume and the standoff height have a positive effect on the lifetime, and the average strain and strain energy values decrease with an increase in the solder volume and component standoff height. Therefore, our fourth and fifth hypotheses are supported.
The component size had a slightly negative effect on the lifetime in this comparative study. Consequently, Hypothesis 3 is supported.
Regarding the different geometry input methods, the results given in the strain and strain energy plots as a function of the standoff height and solder volume indicate that two important parameters impact the strain, strain energy and lifetime: (1) the solder volume in the critical portion of the joint and (2) the component standoff height. With a minimum value for the standoff height and solder volume in the critical lower portion of the solder joint, a CAD-based substitute geometry can be used instead of Surface Evolver. This substitution works well if the standoff height is equal to or greater than m in the case of C0603m and if the standoff height is greater than m in the case of the C1005m capacitor. If the standoff height is too small, the simplified geometry cannot be applied as a substitution for the Surface Evolver model. For small standoff heights, the sensitivity to the geometry modeling type is large. Consequently, the first hypothesis is supported with restrictions.
Overall, in regard to the geometry input methods and the predicted lifetimes, the modeling of the REC geometry input overestimates the lifetime. Relying on overestimated lifetimes can create risks regarding the derived conclusions and can lead to premature failures. The ELL geometry input underestimates the lifetime, since, in this case, the ellipsoid geometry constraint does not allow the solder to fill the volume in the lower portion of the component. In this case, the solder volume will always be smaller, and the predicted lifetime will be also smaller, making the model inaccurate. The EXT geometry is easy to generate in CAD software, since it is a simple extrude command, and it is applicable in most of the cases investigated in this paper.
Hypothesis 2 is not supported, which suggests that more advanced material models result in less accurate modeling of the geometry.
Lastly, it is important to note that the accuracy of the material and fatigue models can be only validated by comparing the obtained results to experimental results.