Stochastic response analysis for nonlinear vibration systems with adjustable stiffness property under random excitation

Nonlinear vibration systems with adjustable stiffness property have attracted considerable attentions for their prominent broadband performances. In the present manuscript, we consider the stochastic dynamical systems with adjustable stiffness and proposed a numerical method for the random responses analysis of the Gaussian white noise excited systems. A multi-dimensional Fokker-Plank-Kolmogorov equation governing the joint probability density of the mechanical states is derived according to the theory of diffusion processes. We solve the multi-dimensional equation using a splitting method and obtained the stationary probability densities and the mean-square responses directly. Two classical nonlinear vibration systems with adjustable stiffness, including the energy harvesting system and the Duffing system with Dahl friction, are presented as examples. Their comparisons with the results from Monte-Carlo simulations illustrate the effectiveness of the proposed procedure for both monostable and bistable cases, even for cases with strong excitation. In addition, the splitting method is efficient for higher-dimensional problem and has advantages of simple implementation, less storage of intermediate values and so on. Hence, in terms of the application scope, the proposed procedure is superior to the current mainstream methods for the random response evaluation of nonlinear vibration systems with adjustable stiffness.


Introduction
Nonlinear systems with noise have attracted considerable attentions over the last few decades. In engineering applications, the vibration systems such as vibration friction systems, energy harvesting systems and vibration suppression systems are all under noise excitation. Meanwhile, the nonlinear vibration systems under random excitation are also utilized in fault diagnosis, parameter identification and structure optimization. Furthermore, those noise-excited systems are generally called stochastic systems [1][2][3][4][5], which are widely studied in response [6][7][8][9][10][11][12][13][14][15][16][17][18][19], stability [20][21][22][23][24], reliability [25][26][27] and so on. Among these studies, the stochastic response evaluation is the most common one. PLOS  In order to master the stochastic response property and to figure out the features of the vibration structure, such as the probability density functions (PDFs) and the mean square responses of nonlinear systems for random excitation, most researches focused on the computation method for the random vibration systems. For example, F. Rudinger [6] utilized the equivalent nonlinearization method to obtain the frequency response function and power spectral density function of the response for nonlinear systems under white noise excitation. Xu et al. [7] used the same approach to investigate the random vibration with inelastic impact subject to Gaussian white noise, so as to analytically obtain the joint probability density of nonlinear system and the statistics of system response. Jin et al. [8] also employed the equivalent nonlinearization technique and imported the generalized harmonic transformation to establish a semi-analytical solution of random response for nonlinear vibration energy harvesters subjected to Gaussian white noise excitation. One of the benefits of this approach was the applicability to the strongly nonlinear, self-excited or parametric systems compared with the equivalent linearization technique. The path-integral method showed its advantages in the computational efficiency, especially for the low-dimensional problems in engineering practice. Hence, Dimentberg et al. [9] applied the path integration method to study the SDOF stochastic vibroimpact problems and acquired the PDFs of the responses. Narayanan and Kumar [10] combined the modified path integration and the finite difference method to solve the Fokker-Planck equation of nonlinear systems subjected to random and harmonic excitations. Zhu and Duan [11] utilized the same method to research the PDFs of the non-linear ship rolling issues. Considering both the richness of the theory and the accuracy of the solution, the stochastic averaging method had prominent position in the stochastic response evaluation. Thus, Xu et al. [12] employed the stochastic averaging method to study the probability density and statistics for random responses of vibro-impact systems with inelastic contact. Gu and Zhu [13] also analyzed the random responses of the vibro-impact systems under Gaussian white noise by using the stochastic averaging method for quasi-Hamiltonian system. Anh and Hieu [14] obtained the response of the Duffing oscillator under combined periodic and random excitations by combining the stochastic averaging method, the equivalent linearization and the finite element method. Jin et al. [15] investigated the PDFs for transient responses of non-linear stochastic systems through the stochastic averaging and Mellin transform. Yang et al. [16] applied the stochastic averaging method for energy envelope to study the PDFs of stationary responses for a vibroimpact Duffing system with bilateral barriers under external and parametric Gaussian white noises. Jiang and Chen [17] utilized the stochastic averaging based on generalized harmonic functions to evaluate the responses of the energy harvesting systems under Gaussian white excitations. In addition, many other effective methods have been proposed. For instance, Rong et al. [18] used the method of multiple scales to investigate the resonant resonance response of a non-linear vibro-impact oscillator subjected to combined deterministic harmonic and random excitations. Yang et al. [19] employed the perturbation method to analyze the stationary responses of Rayleigh vibroimpact oscillator under parametric Poisson white noise.
From the investigations and results in previous studies, the equivalent nonlinearization technique [6][7][8] and the stochastic averaging method [12][13][14][15][16][17] have been widely applied in the random response evaluation due to their accuracy and simplicity. However, both of them have disadvantages that the approximate equivalent exists and is required necessarily, which may influence the applicability of these methods and the accuracy of the results. As an illustration, the stochastic averaging method can be only adopted to predict the random responses of light damping systems subjected to weak excitations. In order to overcome these shortcomings, we apply the splitting method [27][28][29][30] based on the Fokker-Plank-Kolmogorov equation of the original system in this paper. The method is a technique for separating multi-dimensional spatial operators into a sum of one-dimensional operators which is simple to implement, of less intermediate values storage, and flexible with respect to different boundary conditions. Furthermore, it can be efficiently and straightforwardly extended to higher-dimensional, strongly nonlinear, and strong excitation problems.
In the present manuscript, the stochastic responses of the nonlinear vibration systems with adjustable stiffness property for random excitation are investigated using the splitting method. Firstly, we derive a multi-dimensional Fokker-Plank-Kolmogorov equation governing the joint probability densities of the nonlinear systems according to the theory of diffusion processes. Secondly, we utilize a splitting method to solve the multi-dimensional equation and directly obtain the mean-square responses by the joint probability densities. Two classical nonlinear systems with variable stiffness are presented as examples. One is the energy harvesting system, and the other one is the Duffing system with Dahl friction. Finally, we compare the results obtained by splitting method with those from Monte-Carlo simulations (MCS) and evaluate the effectiveness and the applicability of the proposed procedure for both mono-stable and bi-stable cases.

Stochastic responses of adjustable stiffness systems
Consider a nonlinear vibration system with adjustable stiffness described by the following stochastic differential equations, in which, X represents the relative displacement, 2ξ is the viscous damping coefficient, f(X) denotes the nonlinear adjustable stiffness determining whether the system is mono-stable, bistable, or otherwise (subsequent details will be given in section 3), α and β are the non-dimensional parameters, W(t) is Gaussian white noise in the sense of Stratonovich with zero mean and correlation function R(t 2 −t 1 ) = 2Dδ(t 2 −t 1 ), in which δ(•) is the Dirac delta function and 2D is the noise intensity. Eq (1) can be written as Itô stochastic differential equations, where B(t) is an unit Wiener process. Based on the theory of Markov process [31], the transition probability density p(x,y,z,t|x 0 , y 0 ,z 0 ,t 0 ) is governed by the Fokker-Planck-Kolmogorov (FPK) equation corresponding to Itô Eq (2), given by where a 1 ¼ y The initial condition for Eq (3) is and the boundary condition is where s is the boundary surface. The partial differential Eq (3) with the initial and boundary conditions can be solved by using splitting method (see S1 File) which is a technique for separating multi-dimensional spatial operators into a number of one-dimensional operators. According to the method, we obtain the following separating operators and difference equation, in which, p m , p m+1/3 , p m+2/3 , and p m+1 are values of function p(x,y,z,t|x 0 ,y 0 ,z 0 ,t 0 ) at t = mΔt, (m +1/3)Δt, (m+2/3)Δt, and (m+1)Δt, respectively. Eq (7) can be further written in finite difference form with respect to space, and solved by using chasing technique (see S2 File). We acquire the stationary probability density p(x,y,z) while t ! t n , in which t n is a relatively large number. Thus, the stochastic response of the nonlinear system with adjustable stiffness property, including p(x,y),p(x),p(y),E(X 2 ),E(Z 2 ) and so forth, are obtained by in which, p(x,y) is the joint probability density, p(x) and p(y) are the marginal probability densities, E(X 2 ) and E(Z 2 ) are the mean square values of random responses.

Example 1: Vibration energy harvester
The first example considers the model of a generic piezoelectric vibration energy harvester [8,[32][33][34][35] to be simplified as a base-excited one-degree-of-freedom system coupled to a capacitive energy harvesting circuit, as shown in Fig 1. The coupled equation governing the mechanical states and the voltage can be written as where " X represents the relative displacement of mass m and " V is the voltage measured across the equivalent load resistance R. c is the linear viscous damping coefficient, θ is the linear electromechanical coupling coefficient, and C p is the piezoelectric capacitance. € " X b is the base acceleration, and the dot denotes the derivative with respect to time τ. The function " U ð" xÞ represents the potential energy of the mechanical subsystem of the form in which, k 1 > 0 and k 3 > 0 are the linear and nonlinear stiffness coefficients, respectively, and r is an adjustable parameter of linear stiffness. The coupled equations can be further non-dimensionalized by introducing the following non-dimensional parameters, where l c is a length scale, ω 1 represents the fundamental frequency of the degenerated linear mechanical subsystem with k 3 = 0 and r = 0. With these transformations, we express the nondimensional coupled equations as with In this manuscript, the base acceleration excitation € X b is assumed as Gaussian white noise which is defined same as W(t) in Eq (1). The output power " P ¼ " V 2 =R is the most important quantity in energy harvesting. By introducing the transformation P ¼ " the non-dimensional expression of the output power as The parameters of the energy harvesting system (12) are set according to the research of Daqaq [36], the adjustable linear stiffness coefficient r = 0 for mono-stable system and r = 1.1 for bi-stable system, the nonlinear stiffness parameter δ 3 = 0.5, the time constant ratio α = 0.05, the linear non-dimensional electro-mechanical coupling coefficient κ = 0.75, the mechanical damping ratio ξ = 0.01, except as specially provided. Set the time step Δt = 0.02 and the calculating time T = 300 (see S3 File), we obtain the stationary probability densities of  Fig 2. Obviously, as the excitation intensity 2D changed, the responses for mono-stable system differ in magnitude but not in trend. The Monte-Carlo simulations (MCS) are carried out to evaluate the accuracy of the proposed procedure, and the results from MCS are represented by the symbols, which also applies for the figures below. It is obvious that the numerical solutions agree very well with the results from MCS.

mono-stable system displacement and velocity with different excitation intensity and show them in
We compare the stochastic responses of T = 10, T = 30, T = 300 and show them in Fig 3. In which, T = 10 and T = 30 correspond to transient probability densities of system displacement, while T = 300 correspond to stationary probability density of system displacement. It is easily seen that the splitting method is an ideal approach solving the stationary PDFs of stochastic systems but not appropriate for the transient response evaluation. We explain this phenomenon as follows: for the stationary response, terms on the l.h.s. of Eq (7) are equal to zero as there is no time variation in a stationary PDF, while (p (m+1/3)p (m) )/2 = (p (m+2/3)p (m+1/3) )/2 = (p (m+1)p (m+2/3) )/2 = p, being p the stationary PDF; therefore, the three equations in (7) correspond to assume that the space derivatives with respect to X, Y and Z are separately equal to zero, instead of being equal to zero their summation according to the FPK Eq (3). This is a key assumption of the proposed numerical technique. Fig 4 depicts the dependences of the mean-square displacement E[X 2 ] and the mean output power E[P] for energy harvesting on the adjustable parameter of the linear stiffness r. It is easily seen that the mean output power and the mean-square displacement increase with the parameter r, which due to the fact that the larger the adjusting stiffness parameter, the weaker the restoring force and the larger the displacement.
For the situation of r>1 corresponding to the double-well system, the stationary marginal probability densities and the stationary joint probability density of system displacement and velocity are depicted in Figs 5 and 6 for two typical values of excitation intensity, 2D = 10 −3 (relatively small excitation) and 2D = 10 −2 (normal excitation), respectively. It is revealed in  5 that the energy harvesting system possesses bistable potential shape while the motion frequently switches between two potential wells. Fig 6 depicts that the mechanical system possesses mono-stable potential shape, and the stationary probability density of system displacement does not display any bistable potential property. The random responses acquired by the splitting method meet well with those from MCS for the normal excitation (2D = 10 −2 ), Random response of nonlinear vibration systems however, for the situation of relatively small excitation (2D = 10 −3 ), the prediction results on the stationary probability density of system displacement are not very accurate.
In the following section, the influences of the excitation intensity 2D on the mean-square displacement E[X 2 ] and the mean output power E[P] for the bi-stable mechanical system are investigated and shown in Fig 7. With the increase of 2D, the mean-square displacement first decreases and then increases. The results obtained by the splitting method presented low accuracy for the excitation intensity remaining in some certain areas. For arbitrary excitation intensity in the concerned range, however, the splitting method gives high prediction accuracy on the mean output power. Furthermore, the value of the mean output power of the bi-stable system at the adjustable stiffness parameter r = 1.1 is obviously larger than that of the monostable system (0 < r < 1), which means in the random environment, the bi-stable harvesting system with elaborated design parameter can be superior to the mono-stable one.

Example 2: Duffing system with Dahl friction
The second example considers the following Duffing system with Dahl friction subjected to externally random excitation, as shown in Fig 8. The equations of motion of the system are written as where 2ξ is the viscous damping coefficient, k 1 is the adjustable linear stiffness coefficient and k 3 is the nonlinear stiffness coefficient representing the intensity of stiffness nonlinearity, σ 0 Z is the Dahl friction force, f c is the Coulomb friction force, W(t) is Gaussian white noise with intensity 2D. The random responses of the mechanical system (15) are investigated using the splitting method. Numerical results of the mono-stable potential case (k 1 > 0) and bi-stable potential case (k 1 < 0, k 3 > 0) are shown to evaluate the effectiveness and accuracy of the proposed procedure. To illustrate the impact of the Dahl friction on the stochastic responses, the stationary probability densities of the corresponded friction-free systems are provided for comparison.
The system parameters are set according to the study of Wang [37], the excitation intensity 2D = 0.1, the viscous damping coefficient ξ = 0.01, the nonlinear stiffness parameter k 3 = 1, the Dahl friction force parameter σ 0 = 0.06 and Coulomb friction force f c = 0.05, except as specially supplied. We set the time step Δt = 0.05 and the calculating time T = 200, and then acquire the stationary probability densities of system displacement and velocity as shown in Figs 9 and 10 for the mono-stable system and the bi-stable system, respectively. We also illustrate the accuracy of the splitting method via the results from MCS represented by the symbols. It is easily seen that the numerical solutions meet very well with the results from MCS for both the mono-stable and bi-stable systems. Fig 11 depicts the dependence of the mean-square responses on the adjustable linear stiffness parameter k 1 . It is obvious that the acquired results from the splitting method agree well with the MCS solutions for both the mono-stable case and bi-stable case. With the increase of k 1 , the Duffing system with Dahl friction turned from bi-stable case to mono-stable case, and switched at the point k 1 = 0. The mean-square displacement decreases as k 1 increases, however, the mean-square velocity first decreases and then increases.
The relations of the mean-square displacement E[X 2 ] and the mean-square velocity E[V 2 ] to the linear viscous damping coefficient ξ are shown in Fig 12. For both the mono-stable ( Fig  12A) and the bi-stable (Fig 12B) cases, the mean-square responses monotonically decrease with the viscous damping parameter ξ, and the rate of change reduces.

Conclusions
In this article, we investigate the random responses of the nonlinear vibration systems with adjustable stiffness property under Gaussian white noise excitation. Based on the original vibration systems rather than the equivalent systems, we propose a splitting method to solve the FPK equation associated with the Itô stochastic differential equation. We also obtain the stationary probability densities of system displacement and velocity as well as the mean-square displacement and the mean-square velocity.
To verify the accuracy and validate the applicability of the suggested approach, we present two classical nonlinear vibration systems with adjustable stiffness, i.e., the energy harvesting system and the Duffing system with Dahl friction as examples. For both the mono-stable and the bi-stable situation, we compare the results acquired by the splitting method with those from Monte-Carlo simulations. The agreement between the numerical results and MCS results validates the effectiveness of the proposed technique. The influences of system parameters are discussed in detail, including the excitation intensity, the viscous damping coefficient, the nonlinear stiffness coefficient, and especially the linear adjustable stiffness coefficient, on the stationary probability densities and the mean-square responses.
We conclude that the proposed procedure is more advanced than the stochastic averaging technique for the high excitation cases. We are currently working in addressing for engineering interests, such as higher-dimensional problems, strongly nonlinear systems and colored noise excitation situations, and we will report results in future publications.