Issue 
Int. J. Simul. Multidisci. Des. Optim.
Volume 13, 2022
Simulation and Optimization for Industry 4.0



Article Number  7  
Number of page(s)  19  
DOI  https://doi.org/10.1051/smdo/2021038  
Published online  06 January 2022 
Research Article
Advanced Reliability Analysis of Mechatronic Packagings coupling ANSYS^{©} and R
^{1}
Normandie Univ, INSA Rouen, LMN, 76000 Rouen, France
^{2}
Hassan First University, Faculté des Sciences et Techniques, LIMII, BP: 577, Settat, Morocco
^{*} email: hamidelhamdani@gmail.com
Received:
31
August
2021
Accepted:
8
November
2021
The complexity challenges of mechatronic systems justify the need of numerical simulation to efficiently assess their reliability. In the case of solder joints in electronic packages, finite element methods (FEM) are commonly used to evaluate their fatigue response under thermal loading. Nevertheless, Experience shows that the prediction quality is always affected by the variability of the design variables. This paper aims to benefit from the statistical power of the R software and the efficiency of the finite element software ANSYS^{©}, to develop a probabilistic approach to predicting the solder joint reliability in Mechatronic Packaging taking into account the uncertainties in material properties. The coupling of the two software proved an effective evaluation of the reliability of the TCSP using the proposed method.
Key words: R / ANSYS^{©} / chipscale packages / solder joint / kriging metamodel / finite element analysis / MonteCarlo Simulation
© H. Hamdani et al., Published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 chipscale package technology (TCSP) 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 TCSPs 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 TCSP configurations. Consequently, several studies are carried out in order to characterise the impact of TCSP 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 TCSP will be very expensive or even impractical. In this respect, a metamodelbased approach is proposed to effectively overcame the high computation time and assess precisely the reliability of TCSP.
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 finiteelement 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_{0}” 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_{0} (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 temperaturedependent and timedependent 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 preexponential 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_{0} 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, highorder 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 noisefree 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_{1}, … , 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^{(i)}, x^{(j)}) is the correlation function Gaussian process between any two samples x^{(j)} and x^{(j)} and σ^{2} 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^{2}(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), β = (β_{1}, ..., β_{k}) the associated regression coefficients, and r(x) = [R(x, x_{1}), … , R(x, x_{k})]^{T} is the vector of correlation functions between the untried point x and the k samples x^{(1)}, ..., 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^{2}. it gives information on how the model reproduces the observed results. Otherwise, if , y^{(i)}, are respectively the predicted, experimental and the mean of the response. R^{2} expression is given as follow:(13)where 0 ≤ R^{2} ≤ 1. An R^{2} value close to 1 means that the constructed metamodel is more accurate.
2.2.4 Crossvalidation
The main concept of Crossvalidation 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 “leavekout” approach is a kind of crossvalidation where all exluded subsets are size k. In the case where k = 1, the crossvalidation is called “Leaveoneout” [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^{(i)} is the exact response of the finite element model at x^{(i)} and is the prediction of the built metamodel at x^{(i)} based on all the set excluding (x^{(i)}, y^{(i)}).
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 TCSP.
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 Krigingbased 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 TCSP
The studied package is a 13 × 13 mm TCSP. This device is composed of 15 × 15 full ball matrix (225ball), 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 TCSP 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 TCSP 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 TCSP. 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 rateindependent plasticity phenomenon. However, Darveaux [14] has modified Anand™ model in order to consider timedependent and timeindependent phenomena. The modified Anand constants (cf. Tab. 5) are activated for the solder joint material in the studied CSP.
Model material properties.
Tape material properties.
Die attach material properties.
Mold compound material properties.
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 TCSP 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 worstcase 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 (xdimension) was taken at 1.5 of the slice xdimension and the slice model width at onehalf 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 TCSP 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 nonadhesion 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 ChipScale 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 TCSP 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. 
Simulation results of chipscale package.
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 spaceflling DOE, including quasirandom lowdiscrepancy sequences, orthogonal array sampling, and pseudorandom 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 socalled 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 spacefilling 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^{2} = 0.89 for Ball/Test Board Solder Joint and RMSE_{CV} = 3.3, MAE_{CV} = 2.42 and R^{2} = 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 TCSP. 
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/SubstrateSolder Joint and the Ball/PCBSolder 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 QQ plot representing the empirical quantiles (yaxis) against the theoretical quantiles (xaxis).
a PP plot which represents the empirical distribution function computed at each point (yaxis) against the fitted distribution function (xaxis).
In the deterministic simulation using the TCSP FE model, the required computation time is about 4 min (240 s). Therefore, To perform an MCS with 10^{6} samples using the FE model, the computation time reaches 4 × 10^{6} 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 (characteristicfatigue life) of the solder joint presented in Figures 9 and 10 is estimated with 10^{5} sampling points using Monte Carlo simulations. The maximum likelihood estimator is used to find the best probability distribution of fatigue life adjustment. The lognormal 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 Goodnessoffit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Test Board Solder Joint. 
Fig. 10 Four Goodnessoffit plots for various distributions fitted to continuous data (Weibull, gamma and lognormal distributions) in Ball/Substrate Solder Joint. 
Lognormal distribution parameters.
Failure probabilities.
6 Conclusion
In the present paper, we have shown, by means of a probabilistic methodology to assess the reliability of the TCSP 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 TCSP 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
 H. Hamdani, B. Radi, A. El Hami, Metamodel assisted evolution strategies for global optimization of solder joints reliability in embedded mechatronic devices, Microsyst. Technolog. 25, 3801–3812 (2019) [Google Scholar]
 A. Babu Seelam, N. Ahmed Zakir Hussain, S. Hassan Krishanmurthy, Design and analysis of disc brake system in high speed vehicles, Int. J. Simul. Multidisci. Des. Optim. 12, 19 (2021) [Google Scholar]
 Y.H. Pao, A fracture mechanics approach to thermal fatigue life prediction of solder joints, IEEE Trans. Compon. Hybrids Manufactur. Technol. 15, 559–570 (1992) [Google Scholar]
 R. Darveaux, K. Banerji, A. Mawer, G. Dody, Reliability of plastic ball grid array assembly (1995) [Google Scholar]
 A. Makhloufi, Y. Aoues, A.E. Hami, B. Radi, P. Pougnet, D. Delaux, 10study on the thermomechanical fatigue of electronic power modules for traction applications in electric and hybrid vehicles (igbt), in Reliability of HighPower Mechatronic Systems 1, edited by A.E. Hami, D. Delaux, H. Grzeskowiak (Elsevier, 2017), pp. 213–251 [Google Scholar]
 H. Hamdani, A. El Hami, B. Radi, Reliability analysis of tape based chipscale packages based metamodel, Microelectr. Reliab. 102, 113445 (2019) [Google Scholar]
 B.A. Zahn, Impact of ball via configurations on solder joint reliability in tapebased, chipscale packages, in Electronic Components and Technology Conference, 2002. Proceedings. 52nd (IEEE, 2002), pp. 1475–1483 [Google Scholar]
 S. Jayesh, J. Elias, Finite element modeling and random vibration analysis of bga electronic package soldered using lead free solder alloy − Sn1Cu1Ni1Ag, Int. J. Simul. Multidisci. Des. Optim. 10, A11 (2019) [Google Scholar]
 I. Moukhtar, A.Z. El Dein, A.A. Elbaset, Y. Mitani, Introduction and Literature Review, Springer International Publishing, Cham (2021), pp. 1–27 [Google Scholar]
 S. Yu, Z. Wang, J. Fan, C. Qian, Z. Deng, D. Gui, High temperature performance evaluation and life prediction for titanium modified silicone used in lightemitting diodes chip scale packages, J. Electr. Packa. 142 (2020) [Google Scholar]
 J.A. Depiver, S. Mallik, D. Harmanto, Solder joint failures under thermomechanical loading conditions − a review, Adv. Mater. Process. Technolog. 7, 1–26 (2021) [Google Scholar]
 H. Hamdani, Metamodels for the reliability study of mechatronic systems. Phd thesis, Normandie University, 2019. https://tel.archivesouvertes.fr/tel02923958 [Google Scholar]
 M. Jannoun, Y. Aoues, E. Pagnacco, P. Pougnet, A. ElHami, Probabilistic fatigue damage estimation of embedded electronic solder joints under random vibration, Microelectr. Reliab. 78, 249–257 (2017) [Google Scholar]
 R. Darveaux, Effect of simulation methodology on solder joint crack growth correlation, in Electronic Components & Technology Conference, 2000. 2000 Proceedings. 50th (IEEE, 2000), pp. 1048–1058 [Google Scholar]
 L. Anand, Constitutive equations for the ratedependent deformation of metals at elevated temperatures, J. Eng. Mater. Technol. 104, 12–17 (1982) [Google Scholar]
 J. Lin, Modeling test responses by multivariable polynomials of higher degrees, SIAM J. Sci. Comput. 28, 832–867 (2006) [Google Scholar]
 J. Wu, Z. Luo, Y. Zhang, N. Zhang, An interval uncertain optimization method for vehicle suspensions using chebyshev metamodels, Appl. Math. Modell. 38, 3706–3723 (2014) [Google Scholar]
 J. Wu, Z. Luo, J. Zheng, C. Jiang, Incremental modeling of a new highorder polynomial surrogate model, Appl. Math. Modell. 40, 4681–4699 (2016) [Google Scholar]
 C. Huang, B. Radi, A. El Hami, H. Bai, Cma evolution strategy assisted by kriging model and approximate ranking, Appl. Intell. 48, 4288–4304 (2018) [Google Scholar]
 G. Matheron, Principles of geostatistics, Econ. Geol. 58, 1246–1266 (1963) [Google Scholar]
 J. Sacks, W.J. Welch, T.J. Mitchell, H.P. Wynn, Design and analysis of computer experiments, Stat. Sci. 409–423 (1989) [Google Scholar]
 C.E. Rasmussen, C.K. Williams, Gaussian Processes for Machine Learning, MIT Press, Cambridge (2006), vol. 1 [Google Scholar]
 B. Sudret, Metamodels for structural reliability and uncertainty quantification, arXiv preprint arXiv:1203.2062 (2012) [Google Scholar]
 H. Hamdani, B. Radi, A. El Hami, Métamodélisation pour une conception robuste des systèmes mécatroniques, Incertitudes et fiabilité des systèmes multiphysiques 2, 10 (2017) [Google Scholar]
 A. Forrester, A. Keane et al., Engineering Design via Surrogate Modelling: A Practical Guide, John Wiley & Sons (2008) [Google Scholar]
 C. Huang, B. Radi, A. El Hami, Uncertainty analysis of deep drawing using surrogate model based probabilistic method, Int. J. Adv. Manufactur. Technol. 86, 3229–3240 (2016) [Google Scholar]
 C. Ling, Z. Lu, B. Sun, M. Wang, An effcient method combining active learning kriging and monte carlo simulation for profust failure probability, Fuzzy Sets Syst. 387, 89–107 (2020) [Google Scholar]
 M.K. Thompson, J.M. Thompson, ANSYS mechanical APDL for finite element analysis, ButterworthHeinemann (2017) [Google Scholar]
 R.C. Team, R: A language and environment for statistical computing. r foundation for statistical computing, Vienna, Austria, 2013 (2014) [Google Scholar]
 B.A. Zahn, Comprehensive solder fatigue and thermal characterization of a silicon based multichip module package utilizing finite element analysis methodologies, in Proceedings of the 9th International ANSYS Conference and Exhibition (2000), pp. 1–15 [Google Scholar]
 J. Standard, Jesd22a104b, Temperature Cycling, July, 2000 [Google Scholar]
 A. Makhloufi, Y. Aoues, A. El Hami, Reliability based design optimization of wire bonding in power microelectronic devices, Microsyst. Technolog. 22, 2737–2748 (2016) [Google Scholar]
 A.A. Giunta, S.F. Wojtkiewicz, M.S. Eldred et al., Overview of modern design of experiments methods for computational simulations, in Proceedings of the 41st AIAA aerospace sciences meeting and exhibit, AIAA20030649 (2003) [Google Scholar]
 A. Dean, M. Morris, J. Stufken, D. Bingham, Handbook of design and analysis of experiments, CRC Press (2015), Vol. 7 [Google Scholar]
 A. Ben Abdessalem, A. ElHami, Global sensitivity analysis and multiobjective optimisation of loading path in tube hydroforming process based on metamodelling techniques, Int. J. Adv. Manufactur. Technol. 71, 753–773 (2014) [Google Scholar]
 D. Dupuy, C. Helbert, J. Franco et al., Dicedesign and diceeval: Two r packages for design and analysis of computer experiments, J. Stat. Softw. 65, 1–38 (2015) [Google Scholar]
Cite this article as: Hamid Hamdani, Bouchaïb Radi, Abdelkhalak El Hami, Advanced Reliability Analysis of Mechatronic Packagings coupling ANSYS^{©} and R, Int. J. Simul. Multidisci. Des. Optim. 13, 7 (2022)
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 ChipScale 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 TCSP. 

In the text 
Fig. 9 Four Goodnessoffit 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 Goodnessoffit 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.