A Model of BGA Thermal Fatigue Life Prediction Considering Load Sequence Effects

Accurate testing history data is necessary for all fatigue life prediction approaches, but such data is always deficient especially for the microelectronic devices. Additionally, the sequence of the individual load cycle plays an important role in physical fatigue damage. However, most of the existing models based on the linear damage accumulation rule ignore the sequence effects. This paper proposes a thermal fatigue life prediction model for ball grid array (BGA) packages to take into consideration the load sequence effects. For the purpose of improving the availability and accessibility of testing data, a new failure criterion is discussed and verified by simulation and experimentation. The consequences for the fatigue underlying sequence load conditions are shown.


Introduction
Solder joint interconnects serve as electrical connections and mechanical bonds between components and the substrate [1]. Small package size trends dictate leadless solder ball grid arrays (BGA) becoming the mainstream, as a result BGA package reliability has been put at the forefront of research issue.
Fatigue life prediction models have been investigated based on a diverse range of solder joints, all the models require some specific geometry and material related information to establish constitutive equations [2]. Some experts concentrate on the damage mechanisms and conditions, for example: electronic products may experience thermal cycle, power, shock, and vibration during the whole life cycle. Literature [3] shows three quarters of BGA fatigue failures are the result of thermal and vibration load, and researchers developed models by associating the stress-strain data under certain load with damage situation. For now, most of the models applied in engineering follow the Miner linear damage accumulation rule for its convenience. However, some researchers believe that the variation of load conditions also has an effect on the BGA products failure time [4], study [5] indicates that the linear cumulative damage law cannot conform to the experiments well. So they focus on the sequence influence of individual cycles in physical fatigue damage accumulation. A study by Okuyama et al. [6] demonstrates that the load conditions of thermal to vibration are harsher than in reverse. Fatemi and Yang [7] provided an overview on cumulative fatigue damage and life prediction, but these models do not have good practicability.
This paper tries to propose a life prediction model for BGA packages that takes load sequence influence into consideration that is considered more appropriate for application. Actually, it is extremely difficult to monitor and track the initiation and propagation of fatigue cracks of BGA packages, consequently we define electrical resistance strain, whose value is equal to ratio of the resistance variation to initial electrical resistance, and can reflect the quantitative relationship between the extent of solder joint resistance variation and crack expansion as an alternative failure criteria for measuring the physical fatigue damage.

Model Assumption
Failure occurs when the fatigue crack of BGA solders reaches a critical length C a ; the fatigue failure process of BGA solder joints can be divided into initial crack stage and crack propagation stage, the fatigue life of solder joints can be described as Equation (1) shows, where f N means the whole fatigue lifetime of solder joints, o N represents the lifetime of initial cracks, namely the cycles to crack initiation, and p N represents the lifetime of the propagation stage, respectively [8]; assuming that cyclic thermal load conditions in the initial crack stage may affect the life of solder joints. Several studies [9,10] indicate that cyclic thermal load conditions result in material hardening or softening, which may change the material hardness of solder joints. Bao [11,12] points out that when the recovery resistance of material is the same, the elastic modulus is proportional to the square root of the material hardness, therefore the elastic modulus of solder joints may be influenced and subsequently the stress situation is affected. In addition, loads with different amplitudes have different effects on inner defect of solder joints [13], Figure 1 shows that low loading amplitude can help release the inner stress of defects and extend solder joints fatigue life. On the contrary, high loading conditions may promote the expansion of defects and shorten their lifetime. 1 2 3 , , L L L in this figure represent three loading conditions and 1 2 3 , , N N N are their corresponding fatigue cycles, respectively. Assuming that the load condition effects only exist in initial crack stage. Abundant tests [14] show that cyclic softening or hardening situation will be stable within a short time, and the lifetime of the initial crack stage is much shorter than that of propagation stage (nearly 10%-20%). Therefore, we suppose that the load condition effects only exist in initial crack stage to simplify the model.

Base Model for BGA Fatigue Life
According to the Darveaux model [15], the relation between crack growth and final failure time can be described as follow, Assuming that the load condition effects only exist in initial crack stage. Abundant tests [14] show that cyclic softening or hardening situation will be stable within a short time, and the lifetime of the initial crack stage is much shorter than that of propagation stage (nearly 10%-20%). Therefore, we suppose that the load condition effects only exist in initial crack stage to simplify the model.

Base Model for BGA Fatigue Life
According to the Darveaux model [15], the relation between crack growth and final failure time can be described as follow, in which ∆W ave is the incremental viscoplastic strain energy density, da dN is the crack growth rate and k 1 ∼ k 4 are correlation coefficients of material properties.
The Darveaux model has been proven to be available in engineering applications, Equation (3) indicates this model suppose that the crack growth rate is constant, which can be linear fitted by test data. However, research [16] demonstrates that crack growth rate increases with the increasing number of cycles in propagation stage, and the increasing trend is similar to exponential growth, as Figure 2 shows.  k k are correlation coefficients of material properties.
The Darveaux model has been proven to be available in engineering applications, Equation (3) indicates this model suppose that the crack growth rate is constant, which can be linear fitted by test data. However, research [16] demonstrates that crack growth rate increases with the increasing number of cycles in propagation stage, and the increasing trend is similar to exponential growth, as Figure 2 shows. Hence we improve the Darveaux model with the Paris formula, which is widely used in describing the crack propagation stage, as Equation (5) shows.
where C and n are material constants, K  means amplitude of stress intensity factor,   represents the differences between the maximum stress value and minimum stress value under certain load conditions, Y means the geometry influence function whose value is related to the sample size and load form, in our model we set Y = 1 for simplification [17].
Transforming Equation (5) and integrating from an initial crack length to a variable crack length gives ( take a derivative with respect to cycle number,  Hence we improve the Darveaux model with the Paris formula, which is widely used in describing the crack propagation stage, as Equation (5) shows.
where C and n are material constants, ∆K means amplitude of stress intensity factor, ∆σ represents the differences between the maximum stress value and minimum stress value under certain load conditions, Y means the geometry influence function whose value is related to the sample size and load form, in our model we set Y = 1 for simplification [17]. Transforming Equation (5) and integrating from an initial crack length to a variable crack length gives (N = 0, a 0 = 0) take a derivative with respect to cycle number,  where α is related to thermal cyclic stress condition ∆σ in crack propagation stage, thus it can be called in-time load factor. β is denoted as the exponential increasing trend of crack growth ratio, according to the assumption 4 the value of β is related to the stress condition in initial crack stage. When β = 0 and da/dN = α is constant, our model coincides with the Darveaux model. Integrating Equation (8) from zero to critical crack length a C gives equation of N p , thus N f can be refined as

Load Sequence Effects in Model
In daily engineering practice, load sequence effects are often ignored and variable amplitude fatigue life are calculated using the linear damage accumulation rule, but the results rarely coincide with experimentally determined data [18]. With the help of α and β in the base model, we are able to consider the load sequence effects in thermal fatigue life prediction of BGA solder joints.
There are two kinds of thermal cyclic load conditions: Ω 1 and Ω 2 , based on Equation (7) the instantaneous crack growth rates of the sample are, when the crack length a = a 1 = a 2 , integrating Equation (11) gives combine these two equations together, relation formula between N 1 and N 2 is given Assuming that the BGA solder joint sample has experienced N 1 = n 1 cycles since passed initial crack stage in Ω 1 , then the thermal cyclic condition is transformed rapidly into Ω 2 , the instantaneous crack growth rate gives let A = α 2 ) β 1 +1 and the modified crack growth equation of BGA sample can be deduced as follows, finally, we work out the thermal fatigue lifetime prediction model:

Parameters Determination Based on Historical Crack Length Data
It is convenient that determining the values of α and β if the crack length data of a specimen in each load cycle is available. Take the logarithm of Equation (10), a relation between crack length and cycle number can be derived, we can work out the least square estimates of α and β by However, these history data are deficient because these micro defects like crack and void are unobservable without destructive physical analysis (DPA) method, not to mention that this method cannot take sequence effects into consideration. Therefore, improved methods are urgently needed.

Discussion of Numerical Computation Method for Model Parameter Determination
In this paper, we propose a numerical computation method for determining parameters of a BGA life prediction model considering load sequence effects. This method is available without keeping a record of crack length, as long as the cycle numbers in several conditions are provided. It will be more convenient for forecasting the lifetime of BGA packaged products.

Parameters Determining Equations
In order to reflect the load sequence effects, more relations between load condition and life cycle need to be established. At least four equations are needed to determine all the parameters in Ω 1 and Ω 2 as Equation (19) shows. where , the simplification of these simultaneous equations gives relations between β 1 and β 2 , as Equation (20) shows, the approximate values of β 1 and β 2 can be worked out by Matlab software [19] with the test records of cycle numbers in different conditions. Afterwards, α 1 and α 2 can be calculated with Equation (21).

Resistance Strain Theory
By nature of metal conductor resistance, according to solid fracture theory, the quantitative relationship can be revealed for the solder joint resistance strain and the crack expansion. The electrical resistivity of metal material increases as fatigue load cycles. Resistivity changes slightly in the initial crack stage, in the crack propagation stage resistivity varies relatively faster, when the material fails, fatigue damage mutates significantly as well as the resistivity [20].
Based on the Griffith fracture theory, the resistance strain of solder joints can be solved by Equation (22), where R 0 is the solder joints resistances at time point t 0 , and R t is the solder joints resistances at time point t.
the approximate values of 1  and 2  can be worked out by Matlab software [19] with the test records of cycle numbers in different conditions. Afterwards, 1  and 2  can be calculated with Equation (21).

Resistance Strain Theory
By nature of metal conductor resistance, according to solid fracture theory, the quantitative relationship can be revealed for the solder joint resistance strain and the crack expansion. The electrical resistivity of metal material increases as fatigue load cycles. Resistivity changes slightly in the initial crack stage, in the crack propagation stage resistivity varies relatively faster, when the material fails, fatigue damage mutates significantly as well as the resistivity [20].
Based on the Griffith fracture theory, the resistance strain of solder joints can be solved by Equation (22), where 0 R is the solder joints resistances at time point 0 t , and t R is the solder joints resistances at time point t. The geometry model of BGA solder joint is shown as Figure 3, I represents current direction passing the solder joint,  represents the stress in the welding points. Supposing that only part of the damage crack defects, which are perpendicular to the current direction, have effects on the specimen resistance. When the solder joint is undamaged, its resistance can be solved by The geometry model of BGA solder joint is shown as Figure 3, I represents current direction passing the solder joint, τ represents the stress in the welding points. Supposing that only part of the damage crack defects, which are perpendicular to the current direction, have effects on the specimen resistance. When the solder joint is undamaged, its resistance can be solved by ρ is material resistivity, H is solder joint thickness, S is the cross-sectional area that perpendicular to the current, and S = πr 2 , r is the radius of welding points between solder point and PCB. When crack defects grow under load the resistance changes by S(a) is a function of crack length, which represents the equivalent area of crack defects. Equation (24) indicates that when S(a) = 0, the solder joint is undamaged and R = R 0 . When S(a) << S, means that crack is in initial stage, the resistance increases slowly and linearly dependent on the crack growth. The solder joint fractures when S(a) reaches the maximum value and resistance value tends to infinity nonlinearly ( R → ∞ ).
Differentiating Equation (24) with respect to all variables we obtain, In this paper, it is the thermal load conditions that dominate electrical resistance changes, and the influence of dr and dH can be neglected. So Equation (25) can be transformed into Assume that the equivalent area is a circle with radius a, thus S(a t ) = πa t 2 , a t is the crack length at time point t, S (a) = 2π∆a 2 and ∆a is the increment of crack length in ∆t. [21] points out that ∆ρ = ρ 0 × k × ∆T, ρ 0 is the resistivity of material at a certain ambient temperature, k is the temperature coefficient of resistivity, ∆T is temperature range. Thus we have where T 1 is the initial temperature of the thermal cycle load. Crack growth rate [22] can be expressed by da dt = m × ω γ , in which m = 8.34 × 10 −8 , γ = 1.08, ω is related to the thermal stress inside of solder joint, in this paper we assume ω = 0.003 which also follows the reference because of similar geometry and thermal conditions of solder joints. Combined with Equation (5) we have With rem(t, t c ) means taking reminders of t/t c , t c is the time that a single cycle may costs. Besides, parameters in the equation that C = 2.6 × 10 −13 , n = 2.1, Y = 1, and ∆σ can be given by simulation which will be shown in following Section 3.3. As a result, the electrical resistance of solder joint in the (N + 1)th cycle is let t r represents the time of heating-cooling process, t s represents the holding time, T L represents the minimum temperature while T H means the maximum temperature, and the temperature changing rate is b. Equation (30) can be refined into in which

Stress Range Simulation
Mechanical behavior of solder joints under certain load profile can be computed by ANSYS software [23]. In order to analyze the stress and strain situation of BGA under different loading conditions and obtain the equivalent stress range mentioned in the previous section, simulation models of BGA specimens are constructed in software. To ensure the accuracy of simulation results, simulation models need to be in accord with the real specimen.

Geometry Construction
In this paper MROM TF-BGA48 is selected as the research sample, the structure of which is as Figure 4 shows, detail size information is listed in Table 1.

Stress Range Simulation
Mechanical behavior of solder joints under certain load profile can be computed by ANSYS software [23]. In order to analyze the stress and strain situation of BGA under different loading conditions and obtain the equivalent stress range mentioned in the previous section, simulation models of BGA specimens are constructed in software. To ensure the accuracy of simulation results, simulation models need to be in accord with the real specimen.

Geometry Construction
In this paper MROM TF-BGA48 is selected as the research sample, the structure of which is as Figure 4 shows, detail size information is listed in Table 1.    According to the size information, a simulation sample can be constructed in software, as Figure 5 shows. According to the size information, a simulation sample can be constructed in software, as Figure 5 shows.

Initial Simulation Setting
For the purpose of describing the relationship between stress tensor and strain tensor, we can directly invoke a classical Anand constitutive model, which has been embedded into ANSYS. Furthermore, we chose Solid185 instead of Visoco107 for material unit in ANSYS, since the Solid185 can be used in simulating the deformation of more complex materials, such as incompressible elasticplastic material or incompressible hyperelastic material, other material settings in the simulation are listed in Table 2. The full transient dynamic analysis method in software is used for large displacement analysis and the time increment of load substep values 10, other simulation settings like calculating times and convergence condition remain at defaults to ensure the convergence and accuracy of results.

Mesh and Boundary Condition
Both the quantity and quality of mesh are significant to guarantee the accuracy of simulating results, especially the meshes on the welded ball of the specimen. Firstly, we divided a single welded ball model into quarters along two mutually perpendicular planes. Then we trisected the radius of a

Initial Simulation Setting
For the purpose of describing the relationship between stress tensor and strain tensor, we can directly invoke a classical Anand constitutive model, which has been embedded into ANSYS. Furthermore, we chose Solid185 instead of Visoco107 for material unit in ANSYS, since the Solid185 can be used in simulating the deformation of more complex materials, such as incompressible elastic-plastic material or incompressible hyperelastic material, other material settings in the simulation are listed in Table 2. The full transient dynamic analysis method in software is used for large displacement analysis and the time increment of load substep values 10, other simulation settings like calculating times and convergence condition remain at defaults to ensure the convergence and accuracy of results.

Mesh and Boundary Condition
Both the quantity and quality of mesh are significant to guarantee the accuracy of simulating results, especially the meshes on the welded ball of the specimen. Firstly, we divided a single welded ball model into quarters along two mutually perpendicular planes. Then we trisected the radius of a quarter of a welded ball, simultaneously quartering the remaining ball arc. After that, we executed the 'sweep' command to generate sufficient meshes, repeating the operation until the whole welded ball is meshed, as the Figure 6a shows.
To reflect the actual boundary condition of BGA specimens, symmetry plane constraint command is attached to the edge surfaces for prohibiting their displacements and rotary movements in simulation. Additionally, we set the PCB bottom point's DOF (degree of freedom) as zero to avoid rigid body motion as Figure 6b shows. the 'sweep' command to generate sufficient meshes, repeating the operation until the whole welded ball is meshed, as the Figure 6a shows.
To reflect the actual boundary condition of BGA specimens, symmetry plane constraint command is attached to the edge surfaces for prohibiting their displacements and rotary movements in simulation. Additionally, we set the PCB bottom point's DOF (degree of freedom) as zero to avoid rigid body motion as Figure 6b shows.

Thermal Load Profile
The simulation thermal load profiles are set as Figure 7 shows. In 1  , the maximum temperature is 100 °C and the minimum temperature is −60 °C, holding times are 10 min, the rate of temperature change is 1 °C/min. In 2  , the maximum temperature is 80 °C and the minimum temperature is −40 °C, holding time and temperature changing rate are the same.

Thermal Load Profile
The simulation thermal load profiles are set as Figure 7 shows. In Ω 1 , the maximum temperature is 100 • C and the minimum temperature is −60 • C, holding times are 10 min, the rate of temperature change is 1 • C/min. In Ω 2 , the maximum temperature is 80 • C and the minimum temperature is −40 • C, holding time and temperature changing rate are the same.
the 'sweep' command to generate sufficient meshes, repeating the operation until the whole welded ball is meshed, as the Figure 6a shows.
To reflect the actual boundary condition of BGA specimens, symmetry plane constraint command is attached to the edge surfaces for prohibiting their displacements and rotary movements in simulation. Additionally, we set the PCB bottom point's DOF (degree of freedom) as zero to avoid rigid body motion as Figure 6b shows.

Thermal Load Profile
The simulation thermal load profiles are set as Figure 7 shows. In 1  , the maximum temperature is 100 °C and the minimum temperature is −60 °C, holding times are 10 min, the rate of temperature change is 1 °C/min. In 2  , the maximum temperature is 80 °C and the minimum temperature is −40 °C, holding time and temperature changing rate are the same.   Figure 8 is the deformation nephogram of solder joints in BGA package structure under one of our thermal load conditions. It indicates the solder joints that are closer to the edge of package the larger plastic strain they will suffer, which coincides with the experiment results in [24]. Table 3 shows the simulation results of equivalent stress (Von Mises stress) range ∆σ, plastic shear strain range ∆γ, and fatigue life N f . Default Coffin-Manson model is selected in ANSYS to give a preliminary fatigue life in Ω 1 and Ω 2 , respectively. Figure 8 is the deformation nephogram of solder joints in BGA package structure under one of our thermal load conditions. It indicates the solder joints that are closer to the edge of package the larger plastic strain they will suffer, which coincides with the experiment results in [24]. Table 3 shows

Thermal Fatigue Life Computation
According to Equation (31), resistance strain curve under certain load conditions can be graphed as Figure 9 shows. The figure indicates that the variation of resistance strain is reversible, it increases or decreases as the temperature rises or falls during the thermal cycle. The amplitude of resistance strain variation depends on the temperature differences, the values of resistance strain increment are 0.335 and 0.243 in 1  and 2  , respectively.

Thermal Fatigue Life Computation
According to Equation (31), resistance strain curve under certain load conditions can be graphed as Figure 9 shows. The figure indicates that the variation of resistance strain is reversible, it increases or decreases as the temperature rises or falls during the thermal cycle. The amplitude of resistance strain variation depends on the temperature differences, the values of resistance strain increment are 0.335 and 0.243 in Ω 1 and Ω 2 , respectively.
The relationship between resistance strain and number of cycles can be obtained by extracting the information of resistance strain at certain temperature in each cycle. The relation curve under different thermal load conditions Ω 1 , Ω 2 , Ω 1→2 , Ω 2→1 are graphed by programming in Matlab according to the theory described in Section 3.2 and the result are shown in Figure 10. As (a) and (b) describe that there are some irreversible resistance strains during thermal cycle, this can be the results of irrecoverable defects like crack and void. The irrecoverable defects are negligibly small with respect to the recoverable defects in Figure 10, however, when these kinds of As (a) and (b) describe that there are some irreversible resistance strains during thermal cycle, this can be the results of irrecoverable defects like crack and void. The irrecoverable defects are negligibly small with respect to the recoverable defects in Figure 10, however, when these kinds of As (a) and (b) describe that there are some irreversible resistance strains during thermal cycle, this can be the results of irrecoverable defects like crack and void. The irrecoverable defects are negligibly small with respect to the recoverable defects in Figure 10, however, when these kinds of defects accumulate to some extent, a surge may occur in resistance strain. It turns out to be the consequence of crack length reaches the critical length a C (radius of solder points) that surge occurs at the same time.
The simulation of thermal load sequence conditions Ω 1→2 and Ω 2→1 are based on the previous work. Solder joints spend their initial crack phase, which occupies 10% of the whole fatigue lifetime according to the fourth model assumption, under the first condition. We hypothesize that the specimen remains in the first condition until the next 1/10 of fatigue lifetime has passed, the load condition turns into the other condition afterwards. The simulation results showed in (c) and (d) indicate that a crack grows linearly before the surges. Moreover, the crack growth rate varies with the thermal load condition.
According to the simulation analysis, all information of fatigue life needed for calculating the life prediction model parameters are shown in Table 4. Plug these fatigue data into Equation (19), the values of β 1 and β 2 can be solved by Matlab as Figure 11 shows, β 1 = 1.136 and β 2 = 1.052. So we can work out α 1 = 7.8192 × 10 −9 and α 2 = 4.4132 × 10 −9 .
Materials 2016, 9,860 13 of 18 defects accumulate to some extent, a surge may occur in resistance strain. It turns out to be the consequence of crack length reaches the critical length C a (radius of solder points) that surge occurs at the same time.
The simulation of thermal load sequence conditions 1 2   and 2 1   are based on the previous work. Solder joints spend their initial crack phase, which occupies 10% of the whole fatigue lifetime according to the fourth model assumption, under the first condition. We hypothesize that the specimen remains in the first condition until the next 1/10 of fatigue lifetime has passed, the load condition turns into the other condition afterwards. The simulation results showed in (c) and (d) indicate that a crack grows linearly before the surges. Moreover, the crack growth rate varies with the thermal load condition. According to the simulation analysis, all information of fatigue life needed for calculating the life prediction model parameters are shown in Table 4.

Experimental Verification
Studies [25][26][27] demonstrate that daisy chain structure is an effective way to reduce the difficulty of monitoring a batch of solder joints. In daisy chain structure, severe damage in one single solder joint, especially the corners, often leads to the destruction of the entire package. The BGA package sample is designed as Figure 12 shows.

Experimental Verification
Studies [25][26][27] demonstrate that daisy chain structure is an effective way to reduce the difficulty of monitoring a batch of solder joints. In daisy chain structure, severe damage in one single solder joint, especially the corners, often leads to the destruction of the entire package. The BGA package sample is designed as Figure 12 shows. Temperature cycle test is conducted following the 1 2   and 2 1   conditions, we record the resistances of solder joints in daisy chain at T = 298 K every five cycles. The relation curve between resistance strain and cycle numbers can be fitted by the records as shown in Figure 13a, b. After thermal cycle test, the SEM (scanning electron microscope) is used for location the crack defection region, and we observed that cracks arose a lot, especially at the corner solder joints as reference mentioned, as shown in Figure 13c. According to Figure 13, we can see that in the early stage each condition the resistance strain grows slowly and linearly. The values of resistance strain at the condition changing time is similar to the simulation results, and the tendencies of R  approximate to the results in Figure 10. In real experimental environment, fatigue is not the only factor that leads to damage, creep and other mechanisms may also exist. As a result, the experimental value of R  is basically bigger than simulation results, and experiment life is shorter. Temperature cycle test is conducted following the Ω 1→2 and Ω 2→1 conditions, we record the resistances of solder joints in daisy chain at T = 298 K every five cycles. The relation curve between resistance strain and cycle numbers can be fitted by the records as shown in Figure 13a,b. After thermal cycle test, the SEM (scanning electron microscope) is used for location the crack defection region, and we observed that cracks arose a lot, especially at the corner solder joints as reference mentioned, as shown in Figure 13c.
According to Figure 13, we can see that in the early stage each condition the resistance strain grows slowly and linearly. The values of resistance strain at the condition changing time is similar to the simulation results, and the tendencies of ε R approximate to the results in Figure 10. In real experimental environment, fatigue is not the only factor that leads to damage, creep and other mechanisms may also exist. As a result, the experimental value of ε R is basically bigger than simulation results, and experiment life is shorter. Temperature cycle test is conducted following the 1 2   and 2 1   conditions, we record the resistances of solder joints in daisy chain at T = 298 K every five cycles. The relation curve between resistance strain and cycle numbers can be fitted by the records as shown in Figure 13a, b. After thermal cycle test, the SEM (scanning electron microscope) is used for location the crack defection region, and we observed that cracks arose a lot, especially at the corner solder joints as reference mentioned, as shown in Figure 13c. According to Figure 13, we can see that in the early stage each condition the resistance strain grows slowly and linearly. The values of resistance strain at the condition changing time is similar to the simulation results, and the tendencies of R  approximate to the results in Figure 10. In real experimental environment, fatigue is not the only factor that leads to damage, creep and other mechanisms may also exist. As a result, the experimental value of R  is basically bigger than simulation results, and experiment life is shorter.

Conclusions
Traditional fatigue prediction approaches are based on isothermal fatigue test data, but basically they cannot be appropriately used for engineering applications. Actually, most of the approaches ignore the load sequence influence and fatigue test data for mini-sized devices is hard to obtain, as a result of that, errors may be introduced if the data is not accurate.
In this article, we develop a new BGA life prediction model considering load sequence effects. Specifically, the analytic mathematic model is derived based on Darveaux model and Paris formula, whose parameters are determined by numerical simulation approaches. A resistance strain theory is proposed to reflect the crack growth progress by more convenient resistance data from fatigue test, and conducted by Mat lab software; the stress and strain situation that is needed for the theory is analyzed by modeling the BGA in ANSYS software. This work can help decrease the difficulty of getting fatigue test data effectively. Moreover, the resistance strain is discussed as microstructure damage related failure criteria and verified by experiment. However, the infinite element simulation of our model in ANSYS was able to analyze the stress and strain situation under different loading conditions, but cannot take the crack influence into consideration. Further research will be focused on improving the model, and our life prediction model may be a promising tool for future research in electronic device monitoring and life prediction.
Author Contributions: Weiwei Hu and Yaqiu Li conceived and designed theory and experiments, then performed the experiments; Yaqiu Li analyzed the data; Yufeng Sun contributed experiment and material tools, Ali Mosleh provided feedback on the described experiments and the writing of the manuscript; Yaqiu Li wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
Traditional fatigue prediction approaches are based on isothermal fatigue test data, but basically they cannot be appropriately used for engineering applications. Actually, most of the approaches ignore the load sequence influence and fatigue test data for mini-sized devices is hard to obtain, as a result of that, errors may be introduced if the data is not accurate.
In this article, we develop a new BGA life prediction model considering load sequence effects. Specifically, the analytic mathematic model is derived based on Darveaux model and Paris formula, whose parameters are determined by numerical simulation approaches. A resistance strain theory is proposed to reflect the crack growth progress by more convenient resistance data from fatigue test, and conducted by Mat lab software; the stress and strain situation that is needed for the theory is analyzed by modeling the BGA in ANSYS software. This work can help decrease the difficulty of getting fatigue test data effectively. Moreover, the resistance strain is discussed as microstructure damage related failure criteria and verified by experiment. However, the infinite element simulation of our model in ANSYS was able to analyze the stress and strain situation under different loading conditions, but cannot take the crack influence into consideration. Further research will be focused on improving the model, and our life prediction model may be a promising tool for future research in electronic device monitoring and life prediction.
Author Contributions: Weiwei Hu and Yaqiu Li conceived and designed theory and experiments, then performed the experiments; Yaqiu Li analyzed the data; Yufeng Sun contributed experiment and material tools, Ali Mosleh provided feedback on the described experiments and the writing of the manuscript; Yaqiu Li wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

N f
the whole fatigue lifetime of solder joints, also the stress cycles of a specified character that a specimen sustains before failure of a specified nature occurs N f i the total fatigue life cycles of ball grid arrays under the ith thermal load condition N f i→j the total fatigue life of ball grid arrays consists of certain fatigue cycles under the ith thermal condition and the rest cycles under the jth condition N o the initial crack lifetime of solder joints, also the cycles to crack initiation