The determination of the activation energy for a vibro-impact system under multiple excitations

In this paper, the stochastic stability of a vibro-impact system with multiple excitation forces is studied. Due to the multiple external excitations and the coexistence of metastable states, the solution of each attractor’s activation energy, which is specifically used to characterize the attractor’s stochastic stability, is much more suitable for the stability analysis rather than the solution of the probability density function. Based on the large deviation theory, the asymptotic analysis is carried out, and a time-varying Hamilton’s equation for the quasi-potential of the vibro-impact system is derived. To verify the effectiveness of the theoretical analysis, two detailed examples, where an impact attractor and a non-impactor coexist in the system, are conducted. By the application of the action plot method, the activation energies and the most probable exit paths for each attractor are derived. Compared with the numerical simulation, the results show very good agreement. Moreover, it shows that the existence of transient chaos near the attractor could seriously deteriorate the attractor’s stability.


Introduction
Vibro-impact systems are usually of many engineering applications, such as vibratory pile drivers, pie placers, heat exchange tubes in nuclear reactors [1]. Because of its importance in the industry, many studies concerning the vibro-impact system's dynamical mechanisms have been made during the past decades. Some studies focus on the modeling of the vibro-impact system, analyzing the bifurcation and the chaotic behavior, and studying the global dynamical behaviors [2][3][4][5][6][7][8].
Because of the existence of impact, some interesting phenomena such as the corner bifurcation [9], the grazing bifurcation [10], the chatter, sticking, and chaotic impacting motion [11,12] are observed. For a vibro-impact system driven by a harmonic excitation, multistability has been a predominant feature [13,14]. Some studies are devoted to controlling the vibroimpact occurrence, where the strategies of the fast base displacement and the fast time-varying stiffness are introduced [15,16], while some other researchers are tried to control it into a monostable system [14]. Those analytical works are excellent. However, they didn't consider the influence of random perturbation.
As is known, random noise exists widely in the real life and has a powerful effect on dynamic behaviors. Thus, the stochastic stability analysis has become a hot research topic for the vibro-impact system. Up to now, the main part of the stochastic researches of the vibroimpact system is devoted to obtaining the response probability density function (PDF), and various methods are tried [17]. One of the most useful methods is the stochastic averaging method, which has been extensively adopted for the vibro-impact systems with various random excitations for the past years [18][19][20][21][22]. The path integration method [23], the multiple scales method [24], the exponential polynomial closure method [25], and the finite element method [26] are also recently proved useful to predict the response PDF of the vibro-impact system. However, after a review of those papers, it can be found that most of the vibro-impact systems in those works are only subjected to random excitations. This contrasts with the fact that the usual vibro-impact systems in the real life are not only perturbed by random perturbations but also propelled by some external deterministic excitations. When the system is under multiple excitations, it will usually present multistability and complex dynamic behavior, and the above-mentioned methods will be extremely difficult to derive the response PDFs [27,28].
On the other hand, when the system is under multiple excitations and has multiple coexisting metastable states, the solution of the PDF usually will not meet the stability analysis anymore. In this case, the activation energy, quasi-potential, and mean first exit time (MFET) are more suitable to characterize stochastic stability [29][30][31]. In the field of biology and chemistry, when facing such a complex multistable system, scientists prefer to work out the activation energy or quasi-potential to have a stochastic stability analysis or thermodynamic stability analysis [32][33][34][35]. While, in engineering, the relevant research is really scarce, especially for the vibroimpact system [36].
In the present paper, based on the large deviation theory, the stochastic stability analysis for a vibroimpact system with multiple excitations is carried out and the activation energies for the attractors are obtained. The main contribution of this paper is twofold: first, through the asymptotic analysis, a timevarying Hamilton's equation for the vibro-impact system with multiple excitations is derived, and it successfully predicted the activation energies and the MPEPs which haven't been done before. Second, through the comparison of two similar examples, it is found that the existence of transient chaos near the attractor could seriously lower the attractor's stability.
The paper is organized as follows: First, in Sect. 2, the schematic model of the vibro-impact system with multiple excitations is given. In Sect. 3.1, based on the large deviation theory, the asymptotic analysis for the vibro-impact system under weak noise is carried out and a time-varying Hamiltonian equation is obtained. In Sect. 3.2, two similar examples, both have the coexistence of the impact attractor and the nonimpactor attractor while one has a transient chaos near the impact attractor, are picked to test the theoretic analysis. With the help of the action plot method [37], each attractor's activation energy and most probable exit path (MPEP) are derived. In Sect. 4, numerical simulations are conducted, and the theoretical results are verified. At last, conclusions are drawn in Sect. 5.

The vibro-impact system
Vibro-impact systems are encountered in many engineering applications such as rolling contact mechanism, percussive drilling, for which impacting behavior is a part of the original design, or gearboxes, bearings, and rotor systems. Therefore, it is essential to study the dynamics of the vibro-impact system. Here, we consider a vibro-impact system, which is a generic mechanical system with a clearance type nonlinearity, as shown in Fig. 1, and the motion of the system obeys the equation as follows where m, c and k denote the mass, damping coefficient, and linear stiffness. F 0 and x are the amplitude and the frequency of the harmonic excitation, respectively. F 1 is the constant excitation, and FðxÞ is the collision force. Systems governed by equations of this type have been studied extensively, theoretically, experimentally, and numerically [8,28,38]. Dividing both sides of the original system (1) by m, the dimensionless equation of motion is derived: where x n ¼ ffiffiffiffiffiffiffiffi ffi k=m p is the natural frequency and f ¼ are the corresponding forces after dimensionless. Assuming the collision force F is governed by the Hertz contact law [19]. Thus, Here, b l and b r are the left and right distances from barriers, respectively; k l and k r are the left and right contact stiffnesses, which are functions of geometries and material properties of mass and its both side constraints. A detailed description of the contact stiffness can be found in the paper [39]. For convenience, we rewrite the system (2) as follows The paper [14,40] shows when the damping ratio f is low, the existence of multiple coexisting metastable states will be a predominant feature for the system.

Noise perturbed system
In engineering, such a situation is often encountered that a machine that runs smoothly suddenly becomes noisy; after a while, it again runs smoothly, the motion state is alternatively between smoothy and noisy, and there are no changes in the system's parameters in the whole procedure. This is usually the case that the noise induces transitions between the metastable non-impact and metastable-impact states. In a real working environment, the perturbation of the noise is always unavoidable. As is seen, for the deterministic system (4), the existence of multiple metastable states is a predominant feature. If the system is perturbed by random noise, the noise will induce transitions between the attractors. Hence, it is very necessary to study each attractor's stability. Traditionally, efforts are made to obtain the stationary PDF of the system's response process. However, for system (4), due to the coexistence of multiple stable states, the solution of the FPK equation is very difficult. Furthermore, even the PDF is obtained, it still could not depict the stabilities of the attractors well. Therefore, it is very important to determine the activation energy, quasipotential, or the MFET of each attractor in this case. Reviewing from large deviation theory [30], the asymptotic analysis for our system's response in weak noise limit is carried out, and the activation energy and the most probable transition path for each attractor are obtained theoretically.

Theoretical analysis
Motions of the randomly perturbed system are usually governed by stochastic differential equations of the type Here, x is an n-dimensional random process, e is a small parameter used to characterize the noise's weakness, r x ð Þ is a matrix with n rows and m columns, and W t ð Þ stands for an m-dimensional standard Wiener process.
If the stationary probability density p is assumed to exist, then p is governed by the stationary FPK equation As e ! 0, the quasi-potential W associated with the stationary probability density p can be shown to be given by [30,41,42]: C A ð Þ is an arbitrary constant. S T 1 T 2x ð Þ is the action functional of fluctuation along with the trajectoryx, defined by This has a form of Lagrangian L for a classic mechanical system.
In the case of multiple attractors A i coexisting in the deterministic system, a global potential can be constructed by taking the minimum of all local pieces W i at each point x: and the constant C A i ð Þ should be appropriately adjusted relative to each other before the minimum is evaluated.
The minimization in Eq. (8) is taken over all paths x, starting from any point in the attractor A and ending at a point x. The endpoints T 1 and T 2 are left free. Actually, the minimum only is achieved for the pathsx taken over infinite time intervals such thatx À1 ð Þ2 A andx T ð Þ = x, and they correspond to the most probable path of Eq. (5) fluctuating from the attractor A to point x. By the method of steepest descents, several extreme paths correspond to dS=dx¼ 0 can be derived, and the least of the action of the extreme paths corresponds to the minimum in Eq. (8). This leads to the Euler-Lagrange equation for the extreme paths, a n second-order nonlinear partial differential equation given by However, the Euler-Lagrange equation (11) is complicated to be solved and usually is transformed into the Hamiltonian equation by using the Legendre transform. By introducing a generalized coordinate vector,q spanning an n-dimensional space with components q 1 ; :::; q n f gsuch that A Hamiltonian structure can be constructed: Here, the generalized momentum p i is defined by Using Eq. (14) to solve for _ q as a function of p and then substituting the result for _ q into Eq. (13) result in a function depending only on q, p. Differentiating Eq. (14) with respect to time and using Eq. (11) give This immediately leads to the Hamilton's equation, and the extreme paths governed by Eq. (11) can be thus found by solving the 2n first-order differential equations: Supposing our system (4) is perturbed by an additive white Gaussian noise, and the motion of the system follows as Here, the random process n t ð Þ has a zero mean, and an autocorrelation function n t ð Þn t þ s ð Þ h i¼ Dd s ð Þ, and D is the intensity of the noise.
Rewritten the system (17) as type (5) follows: For the drift coefficient in system (18) is periodic with respect to time t, the stationary probability density p x; t ð Þassumed to exist should also be periodic for time t, so as the associated quasi-potential W x; t ð Þ.
For the system excited by a periodic excitation, the study is fixed on the Poincaré section plane R t 0 .
Hence, the action functional in Eq. (9) should be changed into the form: Here, T¼2p=x, and the integer N 1 and N 2 are left free.
Furthermore, for the diffusion matrix (a ij ) in our system is singular, the Lagrangian L should be reformed according to the paper [43]: Appling the coordinates transformation Eq. (12), the Lagrangian L would be rewritten into the generalized coordinate, and replacing it into Eq. (14) obtains: Hence, the Hamiltonian for (18) is found to be The extreme path is then the solution of Hamilton's equations (Freidlin equation) and the corresponding action of the extreme path is given by

Examples
Now, applying the above analysis to concrete examples and setting the parameters of the system are as follows [40] f ¼ 0:01; x n ¼ 1; and the parameters in the contact force f ðxÞ as [19] b Varying the force-frequency x, a bifurcation diagram of the system is obtained, as shown in Fig. 2. As is seen, the existence of multiple coexisting metastable states is a predominant feature of the system. It can be checked that the up branch in Fig. 2a is a non-impact branch, and the low branch an impact branch. As the forcing frequency x increases, the motion of the up branch will collide with the barriers eventually, resulting in all the metastable states becoming impact solutions in Fig. 2b. Without loss of generality, the systems under x ¼ 0:525 and x ¼ 0:57, in which both the impact attractor and the non-impactor attractor have coexisted, are selected as examples. For the other cases, such as the coexistence of multiple impactor attractors, the asymptotic analysis is also applicable without any technical problems. Through the GCM method [44], the global structures of the deterministic systems are shown in Fig. 3. And the phase projections of the corresponding metastable states are shown in Fig. 4. Attractor I is used to denoting the metastable-impact state and attractor N the non-impact metastable state, for brevity. As is seen, the basins of the attractor N and attractor I are separated by the stable invariant manifolds of the saddle point. For, the attraction basin of attractor I is larger than that of attractor N in both cases, it would be estimated that the attractor N should be less stable than attractor I for both cases. Also, it is The phase projections of the metastable states: a x ¼ 0:525: a stable periodic-two impact orbit coexists with a periodic-one nonimpact orbit. b x ¼ 0:570: a stable periodic-one impact orbit, and a stable periodic-one non-impact orbit observed that the attraction basin of attractor N in x ¼ 0:525 is slightly larger than in x ¼ 0:57, so attractor N of x ¼ 0:525 should be expected to be more stable than attractor N of x ¼ 0:57, and inverse for attractor I.
Applying the action plot method [37] to Hamilton's equation (22), trajectories initialed from points (100,000 points) on a small circle (radius = 0.000001) centered at each of the attractors are tracked until their first exit from the attraction domain. Then the actions of the sets of the trajectories are parameterized according to the angular position, and the trajectory with the minimum action is considered to be the most probable exit path. The MPEPs escaping from attractor I to attractor N and in reverse found by the action plot are plotted on the Poincaré section, as in Fig. 5. Moreover, the time histories of the momentum p 2 and the action s are shown in Fig. 6.
As shown in Figs. 5 and 6, as the derived MPEPs of the attractor N spiral out from the basin of attraction hitting exactly near the saddle point, the momentums p 2 fall nearly to zero, in full agreement with the large deviation theory. For attractor I, although the trajectory with the minimum action does not hit on the saddle point so exactly and its momentum p 2 does not drop to zero, it is still viewed as the MPEP because it is usually very close to the theoretic MPEP. In addition, from Fig. 6, it can be seen that the farther the MPEP leaves away from the attractor, the much larger the momentum p 2 is, and the much more action is required. That is why the size of the basin of attraction or the distance between the attractor and the saddle point on the boundary is usually expected as a criticism of the stability of the metastable state [45,46]. For x ¼ 0:525, the corresponding activation energies for attractor I and attractor N are 0.000102 and 0.000097, respectively. For x ¼ 0:570, the corresponding activation energies are 0.00107 and 0.000082 for attractor I and attractor N, respectively. It is consistent with the expectation that attractor N's stability should be recessed because the size of the attraction basin of attractor N decreases as the frequency increases. However, it is noted that the stabilities of attractor I differ very far for both the parameters. That is mainly attributed to the transient chaos existing in x ¼ 0:525, which seriously deteriorates the attractor I's stability. In the transient chaos, there is a dense orbit that can carry the particle to any point in the transient chaos without any external forces. Moreover, for the quasi-potential is Lipschitz continuous [30], hence in the transient chaos, the quasi-potential is an equipotential set (as seen in Fig. 6b, the costed action is very tiny before the MPEP leaves the transient chaos). This will greatly decay the attractor's active attraction basin and make the attractor much closer to the saddle point, resulting in the attractor being very unstable. In the following, two numerical methods are applied to verify the above section's theoretical results. One is to obtain the system's stationary probability distribution by numerical simulation first, then convert the stationary probability distribution according to the WKB approximation [47] to derive the global quasi-potential of the system. Another is to examine the MFETs for each attractor, then, according to the Arrheniustype law [31] to determine each attractor's activation energy.

The probability distribution
Traditionally, the solution of the stationary probability distribution function is the main task for the stochastic system. However, for the complexity of the stochastic system with multiple excitations, very few PDFs can be derived analytically. Even for the numerical Monte Carlo simulation, the solution procedure is challenging and time-consuming. In this paper, the stochastic generalized cell mapping (SGCM) method [44] is applied. Rather than simulating the stochastic system for thousands of periods, it only needs to simulate the system for one period of time to derive the system's one-step transition probability matrix. Then, use the one-step transition probability matrix, the probability distribution of the system at any time could be calculated if the initial distribution is given. Figure 7 shows a coup of probability densities of the displacement at different times with the initial distribution given as a uniform distribution. As seen, the probability densities of the displacement derived by the SGCM method agree well with the MC method. Moreover, Fig. 8 gives the stationary joint probability density of velocity and displacement derived by the SGCM method. For both the parameters, the probabilities mainly focus on attractor I, especially for x ¼ 0:57, the stationary probability is almost all absorbed by attractor I.
As is seen, though observing the probability distribution can help to evaluate the attractors' stabilities, it is a vague understanding. Actually, with the probability distribution in hand, a piece of much more intuitive information about the system's stability can be given by a simple transformation of the probability distribution. In the weak noise limit, the probability is entirely defined by the quasi-potential W and is of Boltzmann form [30,42] W x; y ð Þ % ÀD ln p x; y ð Þ ð26Þ Moreover, the quasi-potential difference between the saddle and the attractor is the activation energy used to depict the attractor's stability. Therefore, by converting the stationary probability distribution derived above according to Eq. (26), the system's global quasi-potential can be presented approximately. Figure 9 shows the approximated global quasipotential for both parameters. The higher the difference between the attractor and saddle is, the more stable the attractor is. For x ¼ 0:525, the quasipotential differences for attractor I and attractor N are 0.0000965 and 0.0000911, respectively. And for x ¼ 0:57, the quasi-potential differences for attractor I and attractor N are 0.00109 and 0.000073. It can be seen that they agree with the theoretical results derived in the above section well.

Mean first exit time
The MFET is another critical quantity to measure the attractor's stability, its ability to survive random perturbations. In the weak noise limit D ! 0, the MFET is usually found to obey the Arrhenius law, which Kramer firstly obtained governing the potential system's transition rates in reaction kinetics [31].
DU plays the role of the activation energy or the height of potential well needed to escape from one attractor to another. Several sets of MFET derived from the numerical MC simulation are shown in Fig. 10. By fitting the data according to Eq. (27), the activation energies for each attractor are estimated: for x ¼ 0:525, the value of the activation energy DU I ¼ 0:000114 is found for attractor I and DU N ¼ 0:000107 for attractor N; for x ¼ 0:57, the activation energies of attractor I and N are DU I ¼ 0:001147 and DU N ¼ 0:000082, respectively. As is shown the results also agree well with the theoretical results.
Accompanied by the MFET, there is another essential quantity, the transition rate. It is inversely proportional to the MFET and has been shown of the form W D $ exp ÀDU=D ð Þ [40]. When the stochastic system attains a metastable state, the transitions Fig. 8 The stationary joint probability density of velocity and displacement derived by SGCM method with D ¼ 0:00002. a x ¼ 0:525 and b x ¼ 0:57 between the attractors will reach a balance. For example, for the parameter x ¼ 0:57, the transition rate W IN D from attractor I to attractor N equals approximately to exp À0.001147=D ð Þ , and the transition rate W NI D from attractor N to attractor I equals approximately exp À0.000082=D ð Þ . When the stochastic system reaches a stable state, the transitions between the attractors will reach a balance W IN D P I ¼ W NI D P N . Here, P I is the stationary probability around attractor I and P N around attractor N. Due to W IN D ( W NI D , P I must be much greater than P N to reach the stable state, resulting in the high probability around attractor I in Fig. 8b, almost 100%. At last, for x ¼ 0:525, several escape realizations before transiting into the domain of the attractor N are plotted in Fig. 11. As is shown, when the particles leave the transient chaos, the escape realizations follow along the most probable exit path seemly in a deterministic way. While, in the transient chaos, due to its sensitivity to the noise, only when D ! 0 will the realizations moving along the MPEP be observed.

Conclusion
In this paper, the stochastic stability of a vibro-impact system with multiple external forces is investigated. For the predominant feature of the coexistence of multiple metastable states, it is tried to work out each attractor's activation energy rather than the PDF to Fig. 9 The global quasi-potential derived from the stationary probability densities. a x ¼ 0:525 and b x ¼ 0:57 have a stability analysis. Based on the large deviation theory, the asymptotic analysis for the vibro-impact system in weak noise limit is carried out, and a timevarying Hamilton's equation for the quasi-potential is derived. By applying the theoretical analysis to concrete examples, the activation energies and the MPEPs of the vibro-impact system are obtained, and it shows good agreement with the numerical simulation. It should be indicated that the asymptotic analysis and the derived Hamilton's equation are not specific for the concrete examples, so no matter what kind of the coexistence of multiple metastable states it is always applicable. The successful prediction of the activation energy and MPEP opens a new window for the study of stochastic vibro-impact systems, especially for those with multiple excitations or multiple stable states. Furthermore, for observation of the dramatic difference of the activation energies between the two similar examples (due to the existence of a transient chaos in one of the examples), it urgently demands that the asymptotic analysis should be carried out when studying the stochastic stability of the vibro-impact system.