# Acceleration Response First Passage Failure Probability Analysis Method for Nonlinear Package

Authors: Dapeng Zhu
Journal: Shock and Vibration
 (2018)

---

## Abstract

The acceleration response first passage failure problem of the nonlinear package base excited by Gaussian white noise is analyzed. The model correction factor method (MCFM) is implemented in conjunction with the first-order reliability method (FORM) to analyze the first passage failure probability of the nonlinear package. The white noise is discretized in standard normal space, and an iterative algorithm is proposed to find the design point of the packaging system. On the design point, the hypersurface representing the limit-state function of the nonlinear package is replaced approximately by a hyperplane representing the limit-state function of an equivalent linear system, and the FORM is employed to calculate the failure probability of the packaging system. The accuracy of this method is verified by crude Monte Carlo simulations. Numerical simulations are carried out to observe the effects of system parameters variations on failure probability which can be used for the improvement of packaging design.


---

## Body

## 1. Introduction

In the process of product distribution from manufacturers to customers, protective package is designed to withstand the real operational conditions to which the product will be subjected. The main operational conditions include drop shock, random vibration from trucks, ships, or airplanes, and static compression. It is important and necessary to analyze the response of the product in the package under various operational conditions, based on the response analysis, and a proper packaging design may be carried out to increase package safety and reliability and avoid over-package costs. In the package distribution process, if the package is excited by drop shock, the motion state of the package changes greatly in relatively short time; therefore, in the drop shock excitation conditions, the package is more inclined to get failed.In recent years, many researchers have devoted their attentions to the package safety and cushion material dynamic properties under drop impact condition. Wang and Jiang [1, 2] used the two-degree-of-freedom mechanical system to represent the packaging system; the response of the key component in the product was analyzed numerically, the dropping damage boundary curve was developed, and the effects of system parameters variations on the damage boundary curve were also observed. Zhong et al. [3] took the friction force into account and proposed a new friction equivalent drop theory; this theory provided a new theoretical basis for experimental analysis of structural cushioning materials. To model the closed-cell foam dynamic properties efficiently in the finite element method, Piatkowski and Osowski [4] presented a modified method for the nonlinear dynamic stress-strain curve determination with the use of cushion curves. In [5, 6], new methods for evaluating cushioning properties in drop shock conditions were also presented. Lu and Gao [7] proposed the visual mass method and analyzed quasi-static compression and drop impact behaviors of the multilayered corrugated paperboard cushioning structure. In order to evaluate safety of the package excited by shock, the nonlinear packaging system was modeled as a lumped mass system, the shock responses of the key component were investigated numerically, and the effects of system parameters variation were also discussed in [8–10].In the package distribution process, most of the time, it is the random vibration that excites the packaging system; if the acceleration response of the package reaches or exceeds a certain fixed threshold, that is, the fragility of the product in the package, package failure will occur. This package failure mode is the first passage failure, also mentioned as the first excursion failure, which is a major failure mode of a structure. The main goal of this paper is to analyze the first passage failure probability under random excitation conditions, and the influence of system parameters on the failure probability was also discussed. Based on the analytical conclusions, the package can be designed properly to prevent failure or to minimize the likelihood of the first passage failure occurrence. For the nonlinear package excited by random vibration, nonlinear stochastic dynamic analysis methods can be employed to estimate failure probability. A variety of methods for nonlinear stochastic dynamic analysis have been explored and presented. In the structural reliability field, first- and second-order reliability methods (FORM and SORM) [11, 12] have been proven to be effective methods, for their simplicity and efficiency. The mirror image excitation method presented by Koo et al. [13] provided an approximate method for the nonlinear system failure probability estimation, and the solution of this method does not require the finding of the design point and response gradient; however, this method is restricted to the displacement first passage of single-degree-of-freedom problems. In the tail equivalent linearization method (TELM) [14, 15], the nonlinear system was equivalently replaced by a linear system at the design point; this method can be adopted in many nonlinear cases with great accuracy, especially for high response threshold cases. In recent years, a promising method for the structural reliability problems is the response surface methodology [16, 17]; in this method, the discrete points on the hypersurface around the design point were discretized and obtained, and the points were used to reconstruct a surface to replace the actual hypersurface relating to the real system.The most important step in the nonlinear stochastic dynamic analysis is the determination of the design point. In this paper, a correction idealized system is developed to equivalently replace the real system, and an iterative algorithm is proposed for determining the design point of the nonlinear packaging system, based on the model correction factor method (MCFM). A hyperplane representing the limit-state function of a linear system is used to replace the hypersurface representing the limit-state function of the actual nonlinear system. The FORM solution of the packaging system is used to estimate the reliability approximately. Through the comparison analysis between the crude Monte Carlo simulation and the results from MCFM-FORM, it is shown the presented MCFM-FORM approach allows accurate first passage failure probability estimation. Finally, the proposed MCFM-FORM is applied to analyze the effects of parameters variations on nonlinear package reliability.


## 2. Acceleration Response First Passage of the Packaging System

The package is modeled by a single-degree-of-freedom system, which in shown in Figure1. The rigid mass in the model represents the product in the package, and the spring, paralleled with a damping, represents the cushion material. In the process of package distribution, the package is excited by the base random acceleration vibration input b¨. In order to apply structural reliability methods to the package reliability problems, it is necessary to discretize the random input process and represent it in terms of a finite number of random variables. For the zero-mean base input, the acceleration excitation b¨ can be represented as(1)b¨t=∑i=1nsitui=st⋅u,where u=u1,u2,…,unT is an n-dimensional standard normal vector, st=s1t,s2t,…,sntT is a row vector of deterministic basis functions dependent on the covariance structure of the input process, n is a measure of the resolution of the representation, and T denotes the matrix transposition. As presented by (1), the randomness of the acceleration input b¨ is determined by the vector u, and the time-evolution property of b¨ is presented by the vector s. At a given time constant t, the excitation b¨ can also be represented by a specified point in the standard normal space u; on the contrary, a given point in the space u denotes a random input.Figure 1
Simplified single-degree-of-freedom package model.From (1), the input of the packaging system is a function of the standard normal variable u, and the acceleration response of the package can also be represented as a function of u. In order to solve the package acceleration response first passage problem, the limit-state function of this problem should be formulated firstly. For the acceleration response first passage of the package failure event, that is, at a specified time t, the acceleration response of the product in the package exceeds a specified product fragility Gc, and the limit-state function for this problem can be presented by(2)gu,Gc,t=Gc−x¨u,t.The limit-state functiong can also be presented in the standard normal space u. For the nonlinear packaging system, g can be represented by a hypersurface in the space u. Therefore, the space u is divided by the hypersurface representing the limit-state function g into two parts: in the first part, g>0, the package is safe, and in the other part of the space, g≤0, the acceleration response first passage of the package failure event occurs and the product in the package gets failed. The limit-state function, the safety region, and the failure region in the u space are shown in Figure 2. For the case g≤0, in all the realizations of u that give rise to the event x¨t,u≥Gc, the one that has the highest likelihood is the one nearest to the origin (Figure 2). This point known as the design point in the theory of structural reliability is defined as(3)u∗=minu,s.t.gu,Gc,t=Gc−x¨u,t=0,where ⋅ denotes the norm of a vector. In the FORM, a hyperplane is used to replace the actual limit-state hypersurface at the design point; applying the FORM to the acceleration response first passage probability estimation, one obtains(4)Prx¨u,t≥Gc≅Φ−βGc,t,where Φ⋅ denotes the standard normal cumulative probability function and βGc,t is the reliability index and is also the distance between the design point and origin; the reliability index can be presented as(5)βGc,t=αGc,t⋅u∗Gc,t=−∇ugu,Gc,t∇ugu,Gc,t⋅u∗Gc,t,where αGc,t is the negative normalized gradient vector of the limit-state function at the design point (Figure 2).Figure 2
Vibration reliability analysis using FORM.In the first step, we assume the packaging system is linear, which means the stiffness coefficientk and the damping coefficient c are constant; therefore, the motion equation of the package can be given by(6)mz¨+cz˙+kz=mb¨,where z=b−x is the relative displacement of the product in the package. Applying the Duhamel integral, we can represent the acceleration response of the package as(7)x¨t=2ξωnz˙+ωn2z=2ξωn∫0tb¨τhz˙t−τdτ+ωn2∫0tb¨τhzt−τdτ=∫0tb¨τht−τdτ,where ωn and ζ are the natural frequency and damping ratio, respectively, and hzt and hz˙t are the unit impulse displacement response and unit impulse velocity response functions; similarly, ht is the unit impulse acceleration response function of the packaging system and can be given by(8)ht=2ξωnexp−ξωntcosωdt+1−2ξ21−ξ2ωnexp−ξωntsinωdt,where ωd=ωn1−ξ2; from (1) and (7), one obtains(9)x¨t,u=∫0tht−τ∑i=1nsiτuidτ=at⋅u,where at=a1t,a2t,…,antT is a vector, the elements of which are convolutions of the shape functions of the excitation and the unit impulse acceleration response function of the system:(10)ait=∫0tht−τsiτdτ,i=1,2,…,n.From (9), one can obtain the limit-state function of the linear package base excited by random vibration by(11)gu,Gc,t=Gc−at⋅u=0.From (11), the acceleration response of the linear packaging system is a linear function of standard normal random variables u. Equation (11) defines the limit-state function surface of the linear package for the given fragility Gc at time t, which defines a hyperplane in the space of u; therefore, the design point can be expressed in the closed form as(12)u∗Gc,t=Gcatat2.And the solution for the reliability index is(13)βGc,t=Gcat.Supposing we know the design pointu∗Gc,t and solving for at from (12), the gradient vector at is obtained as(14)at=Gcu∗Gc,tu∗Gc,t2.This expression indicates that, for a given vibration reliability problem, if the design point is obtained, one can determine the gradient vectorat. On the contrary, using the components of at in (10), it is possible to obtain a set of n equations for determining the IRF of the linear system at a set of n discrete time points.From the discussion above, one can conclude that, for the linear packaging system excited by zero-mean vibration, its acceleration response can be determined by closed-form expression, its limit-state function is a linear function of standard normal random variablesu, and the limit-state function and reliability index can be expressed in the closed form. However, for the nonlinear packaging system, the acceleration response of the system cannot be obtained by steady closed-form expression; accordingly, the limit-state function and reliability index cannot be obtained in closed-form expression either. An attractive and reasonable approach is to find an equivalent linear packaging system to replace the actual nonlinear system. This method [18, 19] can accurately estimate the second moments of the nonlinear packaging system response. However, the accuracy of this method in computing first passage probabilities is not good. Most commonly, the package response is assumed to have the Gaussian distribution, and a fairly good approximation of the response variance is obtained when the system nonlinearity is not severe. However, the accuracy deteriorates as the degree of system nonlinearity increases, particularly for the first passage probability, because the Gaussian assumption is not valid for such cases.The discussion above highlights the importance of determining the design point of the nonlinear packaging system subjected to stochastic excitation. In this paper, a procedure based on the model correction factor method is formulated for determining the design point of the nonlinear packaging system, which can be adopted to analyze the first passage failure probability of the nonlinear packaging system.


## 3. Model Correction Factor Method for Nonlinear System Dynamic Analysis

The MCFM, originally formulated in [20], in nature, is a special form of the response surface method. Instead of selecting a more or less arbitrary mathematical form, in MCFM, the response surface is selected by an equivalent idealized system. The idealized system contains main features of the real nonlinear system. The hypersurface representing the limit-state function of the real system and hypersurface representing the limit-state function of the idealized system are brought to match at the neighbourhood of the design point in the u space. After some correction factor iterations, we can use the idealized system to replace the actual system approximately, and the idealized system can be used to analyze problems in vibration reliability. For the simplicity of the linear system in vibration reliability problems, we usually use a linear system as the idealized system.Letx¨ru,t denote the acceleration response of the real system at the time constant t to the base random vibration and x¨iu,t denote the acceleration response of the idealized system to the same excitation. Accordingly, the limit-state functions of the real and idealized systems gr and gi can be expressed:(15)grx¨ru,t,Gc=Gc−x¨ru,t=0,gix¨iu,t,Gc=Gc−x¨iu,t=0.The correction factorνu is defined as the ratio of the acceleration response of the real system and idealized system:(16)νu=x¨ru,tx¨iu,t.Based on (16), we define a correction idealized system, whose limit-state function can be expressed as(17)giνux¨iu,t,Gc=Gc−νux¨iu,t=0.Substituting (16) into (17), we have(18)giνux¨iu,t,Gc=grx¨ru,t,Gc,which shows that the correction idealized system is completely identical to the real system. The correction idealized system is also difficult to solve in the direct way. However, because the two models describe the same system, the models behave in the same way; in this paper, by the reasonable approximate of the correction factor in the idealized model, we can obtain a good approximation for the real system.At the design pointu∗ of the real system, the correction factor νu can be expressed as the Taylor-series expansion, and we obtain(19)νu=νu∗+∇uνu∗u−u∗+⋯.By approximating the correction factor with the constantν∗=νu∗ and substituting it in (17), the zero-order correction idealized model is obtained:(20)giν∗x¨iu,t,Gc=Gc−ν∗x¨iu,t.From (20), in order to find a good approximate of the idealized system, the task is then to reasonably find the value of the correction factor ν∗ so that the estimated design point ui∗ of the zero-order correction idealized model is as close as possible to the design point u∗ of the real system. The estimation of the correction factor ν∗ can be obtained by the following iterative algorithm:Step 1. Select an idealized system. For the simplicity of the linear system in solving vibration reliability problems, in this paper, an equivalent linear system is selected as the idealized system, and the design point of this linear system is obtained by (12) and is taken as the initial iteration value ui0 of the real system design point. At the time constant t, the acceleration response of the idealized linear system is Gc, that is, x¨i1=Gc.Step 2. Since each point in the standard normal space u can be considered as random excitation, the acceleration response of the real system x¨r1ui0,t at the time constant t excited by the excitation corresponding to ui0 can be achieved by using numerical analysis.Step 3. Using (16), calculate the initial estimate of the correction factor ν1u=x¨r1/x¨i1 and solve the corresponding correction idealized model gi1=Gc−ν1x¨iu,t by FORM to obtain its design point ui1 and the corresponding reliability index βi1 through use of (12) and (13), respectively.Step 4. Taking ui1 as the initial point of iteration, calculate the new correction factor ν2: ν2u=x¨rui1,t/x¨iui1,t, and formulate the corresponding limit-state function of the new correction idealized model: gi2=Gc−ν2x¨iu,t; calculate the design point ui2 and the corresponding reliability index βi2 of this new model by FORM.Step 5. Repeating Step 3 and Step 4, a sequence ν1,ui1,ν2,ui2,… is obtained. If the sequence converges, then νp→1 and uip+1−uip<ε, where ε is a set threshold. Then, uip can be considered as the design point of the correction idealized model. Furthermore, because νp=1, we can conclude that uip is also located on the limit-state surface of the real system. If the acceleration response of the correction idealized model x¨iu,t is not too far from the real model acceleration response, it is reasonable to expect the iteration of uip will lead to the convergence of the design point. On the contrary, if the sequence does not converge, then it suggests that the correction idealized model is too far from the real model. Therefore, it is necessary to select a reasonable correction idealized model.From the steps above, one can obtain the design point of the nonlinear system by the MCFM iteration and FORM, and the above procedure normally converges if the correction idealized model is properly selected.


## 4. Choice of the Correction Idealized Model and Optimization

As described above, the most important step for the design point of the nonlinear system using MCFM is the proper selection of the correction idealized model. The linear system is an attractive selection, whose acceleration response, design point, and reliability index can all be obtained in the closed form, which have been discussed in the preceding part. The drawback of the linear idealized system is that as shown in Figure3, the MCFM simply moves the limit-state hyperplane of the correction idealized model without changing its gradient. A linear system is selected as the correction idealized model; assume the design point of the system is ui1, the corresponding limit-state function of the model is gi1, which is indicated by the dashed line in Figure 3, and the reliability index is βi1. In MCFM, the limit-state function of the correction idealized model can be obtained from (9) and (17):(21)giν∗x¨iu,t,Gc=Gc−ν∗at⋅u=0.Figure 3
Design point determining using the linear correction idealized model.From (21), we can observe that, for the linear correction idealized model, with the variation of ν∗, the hyperplane Gc/ν∗−at⋅u=0 moves in a direction parallel to at until the design point ui∗=Gc/ν∗at/at2 falls onto the limit-state hypersurface of the real model. It is obvious that the perfect estimation of the nonlinear system design point can be obtained only if the fixed vector at is the same as the gradient vector of the real limit-state hypersurface at its design point. Since the design point ui∗ is located on the limit-state surfaces of both the correction idealized and real models, the estimated reliability index βi is an upper bound to the reliability index of the real system; therefore, we obtain(22)β≤βi,where βi is the reliability index corresponding to the design point ui∗ in the iterations and β is the reliability index corresponding to the real design point u∗. From (22), one can conclude that β=minβi. Suppose we define the correction idealized linear system in terms of a set of j parameters, θ=θ1,θ2,…,θj. For a given product fragility Gc, taking into account the significance of (22), we can obtain the optimal correction idealized linear system by solving the optimization problem:(23)θoptGc=argminβiθ,where βiθ is the reliability index of the idealized linear system after achieving convergence according to the proposed algorithm. The optimization problem in (23) can be solved by the simplex method [20]. However, it is necessary to have a good starting point. For this purpose, we select a grid of potential starting points θm,m=1,2,…, compute the corresponding objective values βiθ, and select the starting point with the minimal β value. From [20], we concluded that two points per parameter are sufficient. From the solution of (23), the optimized parameters θopt for the correction idealized linear system can be obtained; the design point u∗ corresponding to this optimized linear model is the design point of the real model.


## 5. Example

The package is modeled by a cubic nonlinear single-degree-of-freedom system, shown in Figure1, and the motion equation of the system is(24)mz¨+cz˙+kz+k3z3=mb¨,where z is the compression of the cushion material, z=b−x; the basic parameters of the system are m = 10 kg, c = 491 Ns/m, k = 930100 N/m, and k3 = 3 × 1010 N/m3, and b¨ is the base acceleration excitation of the system. Assume b¨ is Gaussian white noise with the power spectral density Wf = 0.1 g2/Hz, where g = 9.8 m/s2. Let ti=Δt×i,i=0,1,…,n denotes equally spaced discrete time points. The white noise can be approximately represented by a rectangular wave process defined by(25)b¨t≅b¨i=1Δt∫ti−1tib¨τdτ,ti−1<t≤ti,where Δt is the time interval (in this paper, Δt = 0.001 s) and b¨i denotes the intensity of the rectangular pulse in the time interval ti−1<t≤ti. It is easy to show that b¨i are statistically independent normal random variables having a zero mean and the variance σ=Wf/2Δt = 69.3 m/s2. Introducing the standard normal random variables ui=b¨i/σ, the above discrete representation of the white noise process is found to be of the form (1) with(26)sit=σ,ti−1<t≤ti,0,otherwise.Note that this representation approximates the white noise by a band-limited broadband process with the upper frequency cutoff at approximately 1/(2Δt) Hz.Assume the fragility of the product in the package is 20 g, that is,Gc = 20 g. The equivalent correction idealized linear system in this case can be(27)mz¨+ceqz˙+keqz=mb¨,where the free parameters are θ=keq,ceq. One advantage of using this equivalent idealized linear system is that it is possible to construct correspondence between the properties of the original nonlinear system and the correction idealized linear system. For the cubic nonlinear packaging system in (24), we can expect that ceq≈c and keq>k. Therefore, only one parameter keq needs to be optimized. Let the start point of keq=k; in each iteration, the value of k increases, and the reliability index βi is obtained by MCFM. In the iterations, when the minimum value of βi is obtained, the corresponding k is keq. In this example, the optimization process of k is shown in Figure 4.Figure 4
Optimization process ofkeq.From Figure4, the optimized keq in this example is keq = 955700 N/m, the corresponding reliability index β = 2.6599, and the correction factor ν∗ = 1.0643.Obtaining the correction idealized linear system, in the standard normal spaceu, the hypersurface representing the limit-state function of the real nonlinear system can be replaced by the hyperplane representing the limit-state function of the correction idealized linear system at the design point (Figure 3). Therefore, the limit-state function of the system can be obtained by (21), and the FORM can be used to analyze the acceleration response reliability of the package. The reliability index is obtained by (13), the vector at is obtained by (10), ht is obtained by (8), the vector st can be obtained by (26), and the acceleration response first passage probability of the package is obtained by (4). In this example, the first passage failure probability of the package is calculated and shown in Figure 5. From Figure 5, the first passage failure probability of the package increases approximately in an exponential way with time to a steady state, and the convergence time is less than 0.2 s. To verify the effectiveness and accuracy of the method proposed in this paper, the crude Monte Carlo simulation is carried out to analyze the acceleration response failure probability of the package for comparison. In the Monte Carlo simulation, the number of simulations N = 106, and the simulation of the base acceleration random vibration is realized by the harmonic superposition method [21]. From the comparison between the Monte Carlo simulation result and the analysis result using the MCFM-FORM, we can conclude the nonlinear package reliability analysis method presented in this paper has a good accuracy. The failure probability analysis result obtained by using the ELM method is also shown in Figure 5. The error of the analysis by the ELM is that the response of the nonlinear package excited by Gaussian excitation is non-Gaussian, while the model adopted in ELM can only estimate the second moments of the package response accurately in a general manner and cannot reproduce higher moments of the response. The error of ELM is determined by the extent of the nonlinearity of the package and the value of threshold.Figure 5
Acceleration response first passage failure probability of the nonlinear package.In order to find the effects of system parameters on package reliability, numerical simulations were carried out to observe the effects of parameters variations on the first passage failure probability. The parameters variation ranges werek = 4 × 105∼14 × 105 N/m, k3 = 0.5 × 1010∼6 × 1010 N/m3, c = 100∼1000 Ns/m, and Gc = 10 g∼30 g. The basic parameters are c = 491 Ns/m, k = 930100 N/m, k3 = 3 × 1010 N/m3, and Gc = 20 g. While one parameter varies in the set range, the other parameters remain constant as the basic parameters. We used the MCFM-FORM to analyze the failure probability of the system, and the effects of parameters variations on failure probability are shown in Figure 6.Figure 6
Effects of parameters variations on system failure probability.
(a)
(b)
(c)
(d)From Figure6, we can conclude the following:(1)
With the increase of the stiffness parametersk and k3, the first passage probability of the package also increases. This indicates that, for the package excited by white noise, a softer cushion in the package can lead to smaller package failure probability.(2)
Slight variation of damping can lead to a great change in failure probability, which means the failure probability is very sensitive to damping. We can conclude that the most effective method to improve vibration reliability of the package is to select a cushion material with good damping behavior in the package.(3)
The product fragilityGc also affects the package vibration reliability, and higher product fragility can also improve the package reliability. Therefore, a durable product means good vibration reliability of the package.


## 6. Conclusions

In the package distribution process, most of the time, the package is base excited by random vibration, and it is necessary to analyze the package failure probability under random loads. The theory and method in the structure reliability field are implemented to analyze the failure probability of the nonlinear package.The random loads are assumed to be white noise and are discretized in the standard normal space. A model correction factor method is employed to replace the real nonlinear model by an equivalent idealized model. Because of the simplicity and efficiency of the linear system, an equivalent linear package model is used to replace the nonlinear package model. An iterative algorithm is proposed to find the design point of the nonlinear system. At the design point, a hyperplane representing a linear system is used to replace the real hypersurface representing the nonlinear system. The FORM is used to obtain failure probability of the system. The accuracy and efficiency of the MCFM-FORM presented in this paper are verified by the crude Monte Carlo simulation. Numerical simulations are carried out to analyze the effects of each parameter on the failure probability. It can be concluded from simulations that the cushion material with softer stiffness and higher damping can improve vibration reliability of the package effectively, and the increase of product fragility can also lead to better vibration reliability behavior. The failure probability of the package is very sensitive to damping, and the most effective method to improve the vibration reliability behavior is the selection of cushion materials with higher damping property. The softer cushion material and higher product fragility can also improve the package vibration behavior.In order to improve the accuracy of the package in the distribution process, the following topics are recommended for further research:(1)
Development of the package parameters sensitivity method. From parameters variation analysis (Figure6), it is clear that the failure probability is sensitive to the variation of damping. Sensitivity measures are useful for many purposes, including identification of important sources of uncertainty and optimal package design.(2)
Development of the discrete representation for acceleration excitation with specified PSD characteristics. In order to improve the package failure probability analysis accuracy, consideration of the PSD characteristics of the acceleration excitation is inevitable. There can be many other approaches for representation of a continuous random process. For example, one can use the frequency domain analysis to represent a continuous random process. More investigation on the effective representation of a continuous random process is desirable.



---

*Source: 7513971-2018-07-30.xml*