© H. Hamdani et al., Published by EDP Sciences, 2022

1 Introduction

In engineering simulation, various programs, focused on Mechatronics reliability analysis and based on the finite element method (FEM), are available (i.e., Abaqus, ANSYS^©, etc.) [1,2]. In [3], the FEM is used for thermal fatigue life prediction, in [4], to analyse the reliability of plastic grid array, and in [5], to study the thermomechanical fatigue of electronic power modules. Even though some of these programs allow the user to perform probabilistic and reliability algorithms, it does not provide the opportunity to develop more efficient methodology for reliability analysis. However, R is a powerful statistical computing tool that provides not only a rich library of reliability assessment algorithms and statistical tools previously implemented but the possibility to build and define a user algorithm. Therefore, to get the best from the FEM software ((such as ANSYS^©)) and R while allowing full user control, the ideal solution is to couple them [6]. In addition, since the user can adjust the device geometry as well as the reliability parameters, the proposed global loop has proven to be an effective tool in advanced analysis in mechatronics engineering. The proposed coupling method allows the program to run automatically, and the built script does not need user intervention until one of the stopping criteria is attained or the solution is found.

This article describes the different steps to implement this coupling with an application in reliability assessment of Chip Scale Packaging (CSP). The use of CSP in Mechatronic products is increasing due to the growing demand for portable and smaller devices [7]. Over traditional Ball Grid Array (BGA) packaging [8], the main advantage of CSPs technologies is that they save large space. CSP packages are classified into four types [9]: Lead frame (chips are placed on a frame), Rigid (chips rest on a ceramic or laminated substrate), Flex or Tape where chips are placed on a flexible or tape material, and Wafer Level. In this paper, Tape based chip-scale package technology (T-CSP) has been investigated. In operating conditions, Mechatronic packages are subject to large temperature variations [10]. The coefficient of Thermal Expansion (CTE) mismatch between the various materials used in the package construction leads to thermomechanical stresses in solder joints, which can cause the damage in solder joints and consequently the failure of the entire package [11].

To predict reliability in the T-CSPs packages, the most simulation tools are based on deterministic approaches, which do not take into consideration uncertainty of input parameters, environmental and operational conditions. For example, in [7], package assembly and tape vendors diversity leads to variation in T-CSP configurations. Consequently, several studies are carried out in order to characterise the impact of T-CSP configurations on the solder joints reliability. However, the natural variability in the material parameters and manufacturing process are not taken into account in this kind of study [12]. Therefore, it is necessary to take into account the variability of the parameters through a probabilistic method which evaluates the reliability from random inputs [13]. In this work, the material properties are taken as random inputs and their impact on the solder joints reliability is analysed. More precisely, the used nonlinear finite element simulation requires a considerable computation time to calculate the solder joint fatigue response. Therefore, using MCS to predict the reliability of T-CSP will be very expensive or even impractical. In this respect, a metamodel-based approach is proposed to effectively overcame the high computation time and assess precisely the reliability of T-CSP.

The layout of the document is as follows. in the second section we present the set of models and governing equations used to carry out the thermomechanical study and reliability assessment of electronic devices. In Section 3 the proposed probabilistic methodology was presented. The coupling of the finite element code ANSYS^© and R is described in Section 4 as well as how to manage files. Section 5 present the numerical simulation results to turn out the efficiency of the proposed coupling method. Finally, conclusions about the proposed methodology are provided in the last section (Sect. 6).

2 Governing equations

2.1 Life time prediction model subject to thermal cycles

Solder joint fatigue is the primary thermal cycling failure mechanism for most mechatronic packages. Therefore, to evaluate the reliability of an electronic device, it is essential to calculate the fatigue life cycles of its solder joints.

Solder ball fatigue life prediction requires a combination of thermal fatigue methodology with the finite-element simulations. The fatigue model is usually got using accelerated testing and experimental data. This methodology is employed to compute the number of cycles the package can support before damage. One of the proposed approaches, is the Darveaux's model [14] expressed in two equations, where the number of cycles to crack initiation as well as the crack propagation rate per thermal cycle can be obtained using finite element simulation results.

Equations (1) and (2) allow to compute, respectively the thermal cycles to crack initiation “N₀” and the crack propagation rate per thermal cycle “da/dN”.(1) (2)where ΔW_ave is the element volumetric average of the stabilised change in plastic work within the solder element thickness. The characteristic solder joint fatigue life “α” is the number of cycles to 63.2% population failure. It is expressed as a sum of the number of cycles for crack propagation across the solder ball diameter “a” and N₀ (cf. Eq. (3)).(3)

To calculate precisely ΔW_ave in equations (1) and (2), a finite element model with a detailed description, taking into consideration the temperature-dependent and time-dependent deformation behaviour of the solder joints, is essential. The viscoplastic model introduced by Anand is one of the most powerful time and temperature dependent models of solder joints in power modules [15]. Its equation is defined below:(4)where Q is the activation energy, ε_p is the inelastic strain rate, T is the temperature, A is a pre-exponential factor (1/s), R is universal gas constant, m is the strain rate sensitivity of the stress, ξ represents the materials constant, and “s” is an internal variable which expresses the resistance to plastic deformation. Its evolution equation is given as:(5)where:(6)where “a” represents the strain rate sensitivity of hardening or softening, h₀ is the hardening/softening constant, s^* express the saturation value of s associated with a set of given strain rate and temperature, n represents the strain rate sensitivity for the saturation value of deformation resistance, and is a coefficient.

2.2 Metamodeling techniques

A metamodel can be simply defined as a model of model. In literature, metamodels are treated under various names such as approximator, response surface, simplified model, surrogate model, or emulator. Their construction process is summed up in finding an approximation function that sufficiently gives the relationship between the input data , and their corresponding output obtained using the FEM . various metamodelling techniques have been employed in several research, such as polynomial regression, radial basis functions, support vector regression, and Kriging metamodel. Based on a polynomial to mimic the model behaviour, Polynomial regression (PR), is the most used kind of metamodel. The quadratic response surface is the popular form of the polynomial regression model. However, it is not easy to build a complex nonlinear function with high accuracy using lower order polynomials. To overcome this limitation, high-order polynomials [16] such as Gegenbauer functions, Chebyshev polynomials and the Bernstein polynomials [17] have been used for metamodelling [18]. Support vector regression (SVR) and support vector machines (SVMs) are also a powerful methods in machine learning. From all types of metamodel, the major challenge is to choose the appropriate one for FEM approximation. Kriging, which can be seen as the realisation of a Gaussian process, is proved as a powerful tool to approximate complex nonlinear computer code [19]. It is one of the powerful and most used metamodel. The Kriging is used in this paper.

2.2.1 Kriging metamodel

Kriging is a geostatistical method to interpolate deterministic noise-free data. The kriging theory was formalised by Matheron [20], afterwards, Kriging become a standard technique for building metamodels [21]. It is then used to predict the value of an costly finite element model. This technique is also recognised as Gaussian process regression. Let D = { x₁, … , x_n } , a design of experiment (DOE) with dimension n * d and x ∈ R_d, y(x) is a function of x, and y is a vector of n observed values of y(x) on D. Kriging suppose that y(x) is a realisation of a Gaussian process designated by Y(x), and its expression is as follows:(7)where h(x) represent the mean of Y(x), and Z(x) is a Gaussian process whose mean is zero and its covariance is defined as follow: (8)

with R(x⁽ⁱ⁾, x^(j)) is the correlation function Gaussian process between any two samples x^(j) and x^(j) and σ² is its variance.

The selection of the correlation function is a leading element of kriging. Several covariance functions are given in the literature [22], such as spherical correlation, Gaussian, matter, or exponential functions. The Gaussian correlation function, which allows control of both the range of influence and the smoothness of the approximation model, is the most commonly used method.(9)where θ_k(k = 1, 2, ..., d) are unknown parameters of the correlation function, d is the dimension of design space, and respectively are the k^th component of the sampling point x_i and x_j.

The mean and variance of prediction for any new point x [23] can be, respectively, computed by:

(10)(11)where s²(x) and are respectively the variance of and estimated mean value, R represents the correlation matrix of D, h^T = [h_i], i = 1, … , k a set of basis functions (e.g., polynomial functions), H the matrix corresponding to the values of h^T(D), β = (β₁, ..., β_k) the associated regression coefficients, and r(x) = [R(x, x₁), … , R(x, x_k)]^T is the vector of correlation functions between the untried point x and the k samples x⁽¹⁾, ..., x^(k).

The unknown model parameters are determined using the Maximum Likelihood Estimation (MLE) method [23].

2.2.2 Metamodels validation

Switching to the operating phase, requires validation of the metamodel to ensure its quality and its ability to imitate the behaviour of the FE model. The accuracy of the constructed metamodel is highly dependent on its type, its quality, and the quantity of data used in the building process. Several methods, known as metamodel validation, allow to evaluate the precision of the built metamodel. The most used techniques to evaluate the metamodel and comparing it with others are Cross Validation (CV) and model testing.

2.2.3 Metamodel Testing

Examining the residual errors is the simple and the most used way to assess a model performance [24]. Model testing methods aims to determine the difference between the values predicted by the constructed metamodel and the observed values y. The most popular technique is the so called root mean square error (RMSE):(12)

Furthermore, another commonly used error method is the coefficient of determination R². it gives information on how the model reproduces the observed results. Otherwise, if , y⁽ⁱ⁾, are respectively the predicted, experimental and the mean of the response. R² expression is given as follow:(13)where 0 ≤ R² ≤ 1. An R² value close to 1 means that the constructed metamodel is more accurate.

2.2.4 Cross-validation

The main concept of Cross-validation technique (CV) is to use the same set of samples in both building and validation of the metamodel. Precisely, the set of samples generated by the DOE is randomly divided into equal subsets, then in each iteration, one of these subsets is excluded and the metamodel is built based on the remaining subsets. The CV error is the sum of resulting errors of each iteration between excluded subset and the constructed metamodel [25]. The “leave-k-out” approach is a kind of cross-validation where all exluded subsets are size k. In the case where k = 1, the cross-validation is called “Leave-one-out” [26], and CV prediction of the Root Mean Square Error (RMSE) and the Maximum Absolute Error (MAE_CV) are expressed as follow:(14) (15)where y⁽ⁱ⁾ is the exact response of the finite element model at x⁽ⁱ⁾ and is the prediction of the built metamodel at x⁽ⁱ⁾ based on all the set excluding (x⁽ⁱ⁾, y⁽ⁱ⁾).

3 Reliability analysis based metamodel

Monte Carlo simulation is the most commonly used method in the reliability assessment and uncertainty analysis. Its principle is based on using one of the sampling methods to generate the input range from a distribution of input variables, then the sampling results are used to calculate the responses (outputs) through the finite element model. However. The MCS implementation requires a large number of input samples to get more accurate results. Moreover, multiphysics simulations which use the FEM are more expensive in terms of computation time and finally classical MCS becomes unaffordable. For these reasons, in this paper the proposed Kriging based probabilistic method proposed method use the kriging metamodel to approximate the thermomechanical simulation of T-CSP.

The probabilistic method combines MCS with the Kriging to perform efficiently the reliability analysis [27]. More precisely, the proposed method uses the DOE method to obtain the sampling inputs, and the deterministic FE simulation is performed for all the sampling inputs, then Kriging metamodel is constructed based on the sampling inputs and their corresponding outputs. The last step consists in implementing the MCS using the constructed metamodel. The procedure can be summarised in four stages:

• Generate the sampling inputs using the chosen DOE method;

• Compute the response data corresponding to the obtained inputs through FEM simulation;

• Chose, construct and assess the metamodel quality;

• Use the built metamodel to perform the MCS.

Figure 1 presents the flowchart of the Kriging-based probabilistic method. The implementation procedure requires an interconnection between ANSYS^© [28] and R programming [29] in order to combine the metamodel, MCS and finite element analysis.

Fig. 1 Flow chart of the proposed reliability method.

4 Coupling of the finite element code ANSYS^© and R

The coupling of the two software ANSYS^© and R consists in defining a function in a file R. The role of this function is to call ANSYS^© in batch mode. Moreover ANSYS^© runs automatically under the operating system without intervention of the user. Consequently, the defined function “FEModel.R” allows to run ANSYS^© in batch mode, evaluates the finite element model, and store the results in the output files.

Figure A3 presents an example of the “FEModel.R” function which allows you to read the variables to be varied in the finite element “ANSYS^©” code, build the ANSYS^© input file, execute the built file under ANSYS^© with the system command and return finally the results from the finite element simulation. Figure A1 shows an example of the packages needed for our study that should be installed and loaded into R.

As explained in the previous part, to build a metamodel that can emulate the behaviour of the finite element model, it is necessary to create a design of experiment (DOE) and to calculate the response of the samples generated by the DOE through the “FEModel.R” function. The reading commands must be defined as shown in Figure A2.

5 Numerical simulation

5.1 Studied T-CSP

The studied package is a 13 × 13 mm T-CSP. This device is composed of 15 × 15 full ball matrix (225-ball), with a die size which was measured at 8.24 × 9.12 mm and 0.80 mm pitch (cf. Figs. 2 and 3). In [30], the authors study the impact of T-CSP configurations on solder joint reliability, these various configurations are caused by the variations in both tape vendors and package mounting. Figure 2 presents the stackup layer thickness detail for the PCB. The present study is based on a single configuration [30], where we take into consideration the uncertainties resulting from material properties and the thermal expansion mismatch between different materials composing the T-CSP packages.

Fig. 2 Package outline drawing (a) and Layer dimensions of pcb (b).

Fig. 3 Principle of Tape automated bonding.

5.1.1 Material properties

Several material properties are used In the thermomechanical simulation of the studied T-CSP. They vary between linear and nonlinear, plastic and elastic, dependent or independent on time and temperature (Tables 1–4). The Anand constitutive model [15] is used to represent the deformation behaviour of solder joints.

Anand model which expresses the material viscoplastic behaviour, does not take into consideration the rate-independent plasticity phenomenon. However, Darveaux [14] has modified Anand™ model in order to consider time-dependent and time-independent phenomena. The modified Anand constants (cf. Tab. 5) are activated for the solder joint material in the studied CSP.

5.2 Deterministic case

5.2.1 Solder ball finite element model

Before introducing uncertainties, it is necessary to carry out a deterministic simulation in order to characterise the solder joints behaviour. this kind of nonlinear analysis using the global model of the T-CSP package generates tedious calculations. However, to avoid computational cost, only the diagonal slice of the analysed package was developed [30]. Using the diagonal part (cf. Fig. 2), it is ensured that the worst-case situation is simulated. The meshed 3D diagonal part developed by ANSYS^© [28] is shown in Figure 4 and more details of the configuration are provided in Table 6.

The studied slice model (cf. Fig. 2a) contains all major components and a full set of solder joints, crossing the entire thickness of the device. For the boundary constraints, the model plane is neither a free surface a true symmetry plane, nor, i.e. Therefore, the surface of the slice plane has a free displacement in the y direction. The boundary constraint applied to the numerical slice model are presented in Figure 4. In the present modelling, the length PCB (x-dimension) was taken at 1.5 of the slice x-dimension and the slice model width at one-half of the solder ball pitch. In the diagonal slice model, the ball pitch (1.1314 mm) is the hypotenuse of the real ball pitch (0.80 mm).

Fig. 4 Boundary constraints applied to a studied slice model.

5.2.2 Fatigue life prediction model

Darvaux model presented in Section 2.1 is applied to calculate the number of fatigue life cycles of T-CSP package. In the method application, the simulation results are affected by the singularity issues due to the size of the mesh. However, we must take care of the sensitivity of the finite element FE simulation. On one hand, at the limit between the the copper pad and ball solder, the element thickness needs to be well controlled. On the other hand, to perform the simulation, this controlled element thickness was used to establish the element volumetric averaging of the stabilised change in plastic work. For this reason, Darveaux [14] adopts equation constants for various element thicknesses in the interface.

In our case study, we have taken 0.0254 mm (1 mil) in the element thickness [7]. this value is used for the first two layers at both interfaces of the solder joints (see Fig. 5). However, in the finite element model, because of non-adhesion between tape materials and the solder ball, Darveaux's approach uses a 0.5 mil (0.0127 mm) gap between the solder ball and tape material in the finite element simulation (Fig. 5). Table 7 contains the crack growth correlation constants (K1 through K4) for a 1 mil (0.0254 mm) solder ball element thickness. Figure A4 shows an example of the method implementation under ANSYS^©.

Fig. 5 Modeled ball for Chip-Scale package.

5.2.3 Analysis results and Interpretation

After the modelling of the slice under ansys FE software (cf. Fig. 4), the application of boundary conditions and thermal loading (cf. Fig. 1), the FE simulation was performed in order to get the viscoplastic strain energy at the printed circuit board and the package substrate/solder joints. The obtained results are used to compute the solder joint characteristic life through Darveaux's model. The applied thermal loading in finite element simulation is the accelerated thermal cycling (see Fig. 1), advisable by JEDEC [31]. The profile varies from –40 °C to 125 °C. Figure 6 presents the thermal cycle loading whose maximum value is set as the ANSYS^© zero strain reference temperature [7]. Two cycles are taken as a thermal loading condition in the FE simulation. Figure A4 shows the two first steps of the temperature profile implementation under ANSYS^©.

Figure 7 shows the plastic strain distribution in the solder joints at last loading thermal cycle and Table 8 contains the simulation results of the studied model. In the diagonal section of the package (from the centre, containing the central solder ball), the table shows the position of the failing solder joint, the Delta Plastic Work/Cycle which expresses the variation of viscoplastic deformation energy density per loading cycle, the propagation cycles and crack initiation computed respectively through equations (2) and (1), and finally the characteristic life through equation (3). The results show that the solder joint/PCB failure occurred at the end of the T-CSP slice model in the solder ball number 8 (813 cycles), while Ball/Substrate/Solder Joint failure appeared in the 7th (317 cycles).

Fig. 6 Thermal cycle loading.

Fig. 7 Plastic strain distribution in the solder joints at last loading thermal cycle.

5.3 Probabilistic case

Besides the impact of the variation in package assembly and tape vendors on solder joint reliability [7], the material properties making up each configuration show variations in the form of uncertainties [32]. In the present study, the material properties are taken as uncertain inputs. We have chosen the normal distribution, with a mean value presented in Tables 1 and 2 and standard deviation of 0.02, to describe the variations of these properties.

5.3.1 DOE and metamodels building

The process of metamodel building requires a matrix X of inputs and a vector Y containing their corresponding outputs. These outputs are obtained by passing the matrix X through the FE model. However, to built an efficient metamodel the input matrix must be well distributed over the set of each variable. It is the main purpose of the DOE, which consists in maximising the quantity of information obtained from a limited data [33]. Several DOE approaches can be found in the literature, such as space-flling DOE, including quasi-random low-discrepancy sequences, orthogonal array sampling, and pseudo-random sampling [34].

Latin hypercube sampling (LHS) is one of the most efficient DOE methods [35]. Indeed, If we take the example of a square grid containing sample points, this grid is called a Latin square if each row and column contain only one sample. The generalisation of this idea for n dimension is the so-called LHS. In practice, for a function with m variables, to create an LHS sampling we must divide each variable into n equal intervals, then we create n sampling points so that a Latin hypercube is developed.

In the operating phase, to get the the maximum amount of information using a limited number of sampling points, the analyst must verify that the created experimental design properly covers the experimental domain. Therefore, several criteria are used to study the distance between points to evaluate how close the distribution is to a uniform distribution. The used criteria, to describe the distribution of the sampling points in the experimental domain, can be classified in two categories: the discrepancy measures which aim to quantify how the points distribution is different from an uniform distribution and other criteria computed using the distance between pairs of points.

mindist is a distance measurement which aims to measure the minimum distance between two points in the experimental domain. A close pair of points is expressed by a small mindist value, while a good distribution of the points in the experimental domain is expressed by large value. The maximisation of mindist which is called the maximin criterion [36] is the most used method to optimise the LHS design in order to ensure better space-filling properties. In this paper, we have chosen the function maximinLHS [36] as a DOE method.

In the metamodel building, the number of samples in the created DOE affects directly the metamodel quality. Hence, to define correctly the number of samples, several kriging models are constructed and validated using a multiple number of samples (from 60 to 180).

5.3.2 Validation results

Figure 8 shows the metamodels validation results. The sample size can be set as 180, since the kriging performance does not have a great improvements when the size is greater than 160. The validation results based on 180 samples are RMSE_CV = 1.79, MAE_CV = 1.30 and R² = 0.89 for Ball/Test Board Solder Joint and RMSE_CV = 3.3, MAE_CV = 2.42 and R² = 0.9 for Ball/Substrate Solder Joint. Consequently, The obtained results indicates that the constructed metamodel based on 180 samples is sufficiently accurate.

Fig. 8 Kriging metamodels validation for reliability prediction of T-CSP.

5.3.3 Kriging based MCS

After metamodel building and validation, the last stage of the probabilistic method consist in performing the MCS using the constructed metamodel. Figures 9 and 10 present respectively the MCS results of the Ball/Substrate-Solder Joint and the Ball/PCB-Solder Joint. Under R software, the four graphs of each figure are given by the functions denscomp, qqcomp, cdfcomp, and ppcomp of the package “fitdistrplus”. The four graphs contain:

• a density plot which represents the histogram of the empirical distribution and the density function of the fitted distribution,

• a CDF plot of the empirical distribution and the fitted distribution,

• a Q-Q plot representing the empirical quantiles (y-axis) against the theoretical quantiles (x-axis).

• a P-P plot which represents the empirical distribution function computed at each point (y-axis) against the fitted distribution function (x-axis).

In the deterministic simulation using the T-CSP FE model, the required computation time is about 4 min (240 s). Therefore, To perform an MCS with 10⁶ samples using the FE model, the computation time reaches 4 × 10⁶ minutes, while performing the MCS through the kriging metamodel reduces significantly the calculation cost and makes the probabilistic approach more practical and affordable.

The probability density function of the fatigue life (characteristic-fatigue life) of the solder joint presented in Figures 9 and 10 is estimated with 10⁵ sampling points using Monte Carlo simulations. The maximum likelihood estimator is used to find the best probability distribution of fatigue life adjustment. The log-normal distribution with its two parameters (Mean, Standard Deviation) (see Tab. 9) corresponds to the probability distribution of fatigue life. Table 10 presents the calculated probability of failure for target values of the number of fatigue life cycles.

Fig. 9 Four Goodness-of-fit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Test Board Solder Joint.

Fig. 10 Four Goodness-of-fit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Substrate Solder Joint.

6 Conclusion

In the present paper, we have shown, by means of a probabilistic methodology to assess the reliability of the T-CSP package, how to perform an advanced reliability analysis in the Mechatronic field based on the usage of R and ANSYS^©. The coupling method has been presented, and the necessary files to carry out the reliability analysis of T-CSP are given. The obtained results showed the importance of benefiting, on the one hand, from the richness of R as very powerful software in terms of processing statistical data, developing new complex scripts or using available packages to perform numerical studies either for metamodel fitting and validation or MCS implementation. On the other hand, the finite element software ANSYS^© proves its efficiency to perform the finite element simulation to assess the reliability of CSP through accelerated temperature loading. The method presented in this paper is not limited to CSP packages and can be employed for other Mechatronic devices.

Appendix: A

Fig. A1 Packages installation and loading.

Fig. A2 DOE and Metamodel construction.

Fig. A3 Programming of the FEModel function in FEModel.R.

Fig. A4 The temperature profile implementation under ANSYS©.

References

All Tables

All Figures

	Fig. 1 Flow chart of the proposed reliability method.
In the text

	Fig. 2 Package outline drawing (a) and Layer dimensions of pcb (b).
In the text

	Fig. 3 Principle of Tape automated bonding.
In the text

	Fig. 4 Boundary constraints applied to a studied slice model.
In the text

	Fig. 5 Modeled ball for Chip-Scale package.
In the text

	Fig. 6 Thermal cycle loading.
In the text

	Fig. 7 Plastic strain distribution in the solder joints at last loading thermal cycle.
In the text

	Fig. 8 Kriging metamodels validation for reliability prediction of T-CSP.
In the text

	Fig. 9 Four Goodness-of-fit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Test Board Solder Joint.
In the text

	Fig. 10 Four Goodness-of-fit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Substrate Solder Joint.
In the text

	Fig. A1 Packages installation and loading.
In the text

	Fig. A2 DOE and Metamodel construction.
In the text

	Fig. A3 Programming of the FEModel function in FEModel.R.
In the text

	Fig. A4 The temperature profile implementation under ANSYS©.
In the text