Stationary Response of a Kind of Nonlinear Stochastic Systems with Variable Mass and Fractional Derivative Damping

: Viscoelasticity and variable mass are common phenomena in Micro-Electro-Mechanical Systems (MEMS), and could be described by a fractional derivative damping and a stochastic process, respectively. To study the dynamic inﬂuence cased by the viscoelasticity and variable mass, stationary response of a kind of nonlinear stochastic systems with stochastic variable-mass and fractional derivative, damping is investigated in this paper. Firstly, an approximately equivalent system of the studied nonlinear stochastic system is presented according to the Taylor expansion technique. Then, based on stochastic averaging of energy envelope, the corresponding Fokker–Plank–Kolmogorov (FPK) equation is deduced, which gives an approximated analytical solution of stationary response. Finally, a nonlinear oscillator with variable mass and fractional derivative damping is proposed in numerical simulations. The approximated analytical solution is compared with Monte Carlo numerical solution, which could verify the effectiveness of the obtained results.


Introduction
As the rapid development of Micro-Electro-Mechanical Systems (MEMS) in recent years, the size of micro-devices is declining to micron, even nanometer, scale. In nanometer scale, some particular microscopic phenomena always occurs, such as viscoelastic damping and variable mass. For example, in materials like plastics and nano-wires, the most adequate kind of damping is viscoelastic, which is usually described by the fractional derivative damping [1][2][3][4]. Thus, the dynamics of the micro-system with fractional derivative damping have attracted amount of attention, such as the stability, stochastic response, bifurcation analysis, and so on [5][6][7]. However, except for fractional derivative damping, variable mass is another common phenomenon in MEMS of nanometer scale, which will be considered altogether. According to the MEMS application, the stationary response of nonlinear stochastic systems with a variable mass and a fractional derivative damping is studied in this paper.
Variable mass phenomenon commonly occurs in physics and engineering fields, such as comet, jet planes, and rockets, which are all missing mass in their motions [8]. Within these applications, the variable mass was always modeled by continuous functions in deterministic systems, whose dynamics researches were developed in the last decades [9]. All the above variable mass phenomena are described by deterministic functions according to their inherent regularity. However, in the micro world, variable mass phenomenon is always characterized as a stochastic one. When the size of micro-devices is declining to micron, even nanometer, scale, the random adsorption and desorption of the particles around the micro-devices could significantly increase or decrease its quality, resulting in the variable mass phenomenon [10,11]. It is able to change the dynamical characteristics of the micro-devices and cause the damage or even failure of the device performance.
For example, the variable mass of atomic force microscope probe can affect its detection accuracy [12]. Thus, the stochastic dynamics of nonlinear systems with variable mass have also attracted the attention of scholars.
The first step to study nonlinear stochastic systems with variable mass is to construct a model of variable mass. Since the size and mass of MEMS itself are very tiny, the random adsorption and desorption of surrounding particles would significantly change the system's mass. Thus, this variable mass phenomenon should be described as a stochastic process. The initial models of random variable mass oscillator were always simplified by dichotomous or trichotomous noise [13][14][15][16][17]. However, the variable mass modeling method based on dichotomous or trichotomous noise limits the mass change within two or three states, which can not fully show the randomness of the variable mass of micro-devices. In order to better simulate the randomness of the variable mass, many scholars used Gaussian white noise to model the mass disturbance. Wang et al. [18] utilized Gaussian white noise to describe the random mass disturbance and proposed the stochastic averaging technique of the variable mass system based on Hamilton theory, which provided a method to analyze the random responses of stochastic systems with Gaussian-white-noise-based variable mass. Then, Qiao et al. [19,20] presented a group of stochastic stability conditions of the variable mass system under the Gaussian white noise excitation by using the stochastic averaging method and the maximum Lyapunov exponent method. Li [21] analyzed the stochastic response of a vibro-impact system with variable mass under the Gaussian white noise excitation, as well as its stochastic P-bifurcation characteristics.
Motivated by the above discussion, fractional derivative damping and variable mass are both considered in this paper. In this case, the stationary response of a kind of nonlinear stochastic system with fractional derivative damping and variable mass is studied. Based on the Taylor expansion technique, an approximately equivalent equation of the nonlinear stochastic system is given. Then, the corresponding stochastic averaging method is provided by following the diffusion approach to energy envelope. By constructing the corresponding Fokker-Plank-Kolmogorov (FPK) equation, an approach to calculate its approximate analytical solution of a stationary response is proposed. Besides, the effectiveness of the obtained results is verified by simulation of the Monte Carlo numerical solution.
The rest of this paper is organized as follows: in Section 2, some preliminaries are introduced including the fractional derivative and problem statement; then, the approximately equivalent system is proposed in Section 3; a stochastic averaging method is proposed in Section 4; an example is listed to verify the obtained results in Section 5; and, finally, the conclusions are given in Section 6.

Definition 1 ([31]
). The R-L fractional derivative of a continuous function x : (0, + ∞) −→ R with order q > 0 is defined as: where k is a positive integer and k − 1 ≤ q ≤ k, and Γ(·) is the gamma function, i.e., The G-L definition is a discrete form of fractional derivative, which is equivalent to the discretized R-L definition and is introduced as follows.
In addition, the Caputo definition of a fractional derivative is introduced in the following Definition 3.

Definition 3 ([33]
). The Caputo fractional derivative of a continuous function x : (0, + ∞) −→ R with order q > 0 is defined as: where k is an positive integer and k − 1 ≤ q ≤ k.
The Caputo definition owns the practical physical meaning in the initial values and has more advantages than R-L or G-L definition. Thus, Caputo fractional derivative is employed to describe the fractional derivative damping (viscoelastic damping) in this paper.
Moreover, the Laplace transformation of the Caputo fractional derivative is given as, where L{ ·} is a Laplace transform operator with satisfying X(s) = L{x(t)}, and k is an positive integer with k − 1 ≤ q ≤ k. Especially, with null initial conditions, the Laplace transformation of the Caputo fractional derivative is obtained as

Problem Statement
Consider the following single-degree-of-freedom nonlinear stochastic system [19] with a variable mass and a fractional derivative damping, where differentiable g(x) is the nonlinear restoring force with g(0) = 0 and g (0 , and D f is the intensity of the noise ξ f (t). Besides, m represents a variable mass described by wherem is a positive parameter, ξ m (t) denotes Gaussian white noise with satisfying E[ξ m (t)ξ m (t + τ)] = 2D m δ(τ) and D m is the intensity of the noise ξ m (t). In this paper, it is assumed that both stochastic processes ξ f (t) and ξ m (t) are independent with each other. In addition,ẍ represents the second derivative, and ζ C 0 D q t x denotes the fractional derivative damping force with 0 < q < 2 and ζ > 0, which satisfies the Caputo definition in Definition 3.
To study the stationary response of System (6), the corresponding stochastic averaging method is studied. Due to the variable mass and fractional derivative damping, the classical stochastic averaging method could not be applied directly. Thus, the studied System (6) should be transformed into an approximately equivalent system, such that the equivalent form satisfies an improved stochastic averaging method. In the next section, the approximately equivalent form of System (6) is investigated.

Approximately Equivalent Equation of the Fractional Derivative Damping
According to the above section, the conventional stochastic averaging method is not available for System (6). Thus, the approximately equivalent form of System (6) is studied in this section. First, an approximately equivalent method of the fractional derivative damping C 0 D q t x is proposed as follows. (6) can be approximately rewritten by where r * is the characteristic root with the largest real part of the following characteristic equation, Proof of Theorem 1. The restoring force g(x) is differentiable and satisfies g(0) = 0, g (0) > 0. Thus, the linearized homogeneous equation of System (6) in the neighborhood of the origin becomesmẍ Then, according to the Laplace transformation in Equation (5), the characteristic equation of System (10) can be obtained as Equation (9).
Due to m, ζ, g (0) > 0, the real parts of characteristic roots in characteristic Equation (9) must be negative. Suppose that Ω is the set of the characteristic roots in characteristic Equation (9), and r * is the characteristic root with the largest real part, i.e., r * = argmax r∈Ω Re(r).
In dynamic analysis, r * is a critical index especially for stability analysis. Thus, the approximately equivalent equation of the fractional derivative damping s q can be obtained based on the Taylor expansion technique around the characteristic root r * , i.e., Omitting the higher term of (s − r * ) 2 in System (11), it is approximately equivalent to Denoter * as the conjugate of r * . Then,r * is also a characteristic root and satisfies Equation (12). Because r * andr * both satisfy Equation (12), it gives the following equation, i.e., s q ≈ r q * +r q * 2 + q r q−1 * +r According to the inverse Laplace transformation of System (13), it obtains Equation (8) as follows, This completes the proof.

Remark 1.
The approximately equivalent equation of the fractional derivative damping C 0 D q t x in Theorem 1 is accurate in the neighborhood of the characteristic root r * . Thus, the approximately equivalent Equation (8) is much more effective when the difference of real parts of characteristic roots in characteristic Equation (9) is relatively small.

Approximately Equivalent Equation of the Nonlinear Stochastic System
According to the approximately equivalent Equation (8) of the fractional derivative damping, the approximately equivalent equation of the nonlinear stochastic System (6) could be obtained in this subsection.
The approximately equivalent Equation (18) is close to the traditional form of secondorder stochastic systems. Thus, the corresponding stochastic averaging method is available, which is presented in the next section.

Stochastic Averaging of a Nonlinear Stochastic System with Variable Mass and Fractional Derivative Damping
In this section, a stochastic averaging of a nonlinear stochastic system with variable mass and fractional derivative damping is proposed based on the diffusion approach to the energy envelope. Then, according to the obtained stochastic averaging, the stationary response of such a nonlinear stochastic system is discussed.
According to the approximately equivalent Equation (18), the corresponding Stratonovich stochastic differential equation can be obtained [18], • dB m (t), where B m (t) and B f (t) are independent standard Wiener processes. Because y 1 , y 2 in the above Stratonovich stochastic differential equation are independent of B m (t), B f (t), no Wong--Zakai correction term exists. Then, the corresponding Itô stochastic differential equations of System (18) is given, Define a quasi Hamiltonian of System (20) as, Then, System (20) can be rewritten as the following form, dB m (t), According to the Itô differential rule and Systems (21) and (22), the Itô equation of the quasi Hamiltonian H can be obtained, It is noted that H is a slowly varying process and y 1 is a rapidly varying process. Based on the Khasminskii's theorem [34], the slowly varying process H weakly converges to a Markov process. Then, the Itô equation of the approximate Markov process (20) can be obtained by the time averaging of Equation (23), i.e., In Equation (24), m H and σ H 1 , σ H 2 are respectively drift and diffusion coefficients, which are introduced in [18] and described by and where the region of integration Θ is denoted by and The stationary probability density solution p(H) of Equation (24) can be obtained by solving the FPK equation as follows, Substituting in ∂ ∂t p(H) = 0, the solution p(H) is given by where C is a normalization constant defined by Then, the joint probability density of generalized displacement and momentum is p(y 1 , y 2 ) = p(H) T(H) H=H(y 1 ,y 2 ) .
In addition, the marginal probability densities of the generalized displacement and momentum can be calculated by the following integrations,

Numerical Examples
In this section, an example is given to show the effectiveness of the proposed results. Consider the following nonlinear oscillator with variable mass and fractional derivative damping, The corresponding parameters are initially determined as According to the characteristic Equation (9), the characteristic roots are approximately calculated as −0.035 ± 1.035. The highest differential order of System (34) is 2, which is equal to the number of the corresponding characteristic roots. The characteristic roots of System (34) are approximately −0.035 ± 1.035, who own the same real parts. In this case, the proposed approximate method is accurate. Then, based on Equation (8) and Theorem 1, it has the approximately equivalent equation rewritten by (1.0088 + 0.1ξ m (t))ẍ + 0.0713ẋ + 1.0826x + 2x 3 = 0.2ξ f (t).
which means m e = 1.0088, ζ e = 0.0713, g e (x) = 1.0826x + 2x 3 referring to System (14). Due to Equation (21) According to Equation (30), the stationary probability density solution p(H) of quasi Hamiltonian (36) can be obtained, which is shown in Figure 1. In Figure 1, the solid curve is the analytical solution of p(H) in quasi Hamiltonian (36) based on the proposed stochastic averaging method. The ' * ' is the numerical solution of the stationary probability density solution p(H) according to the Monte Carlo simulations. In addition, Figures 2 and 3 show the stationary probability densities of y 1 and y 2 in Equation (33), respectively. Similarly, the solid curve and the ' * ' represent the corresponding analytical solution and numerical solution, respectively. In Section 3, the approximately equivalent equations of the fractional derivative damping and the mass disturbance are given to deduce the corresponding stochastic averaging method. Multiple approximation methods may have an influence in the accuracy of the obtained method. According to Figures 1-3, the numerical solutions agree well with the analytical solutions, which could verify the effectiveness of the proposed approximate analytical method of the stochastic averaging.

Conclusions
In this paper, a class of nonlinear stochastic systems with variable mass and fractional derivative damping is studied in order to reveal the dynamic properties related to viscoelasticity and variable mass. Especially, the stationary response of the studied model is analyzed based on a group of approximately equivalent systems and stochastic averaging. Firstly, an approximately equivalent equation of the fractional derivative damping is proposed based on the Taylor expansion technique, such that the fractional derivative damping could be transformed into an approximate second-order system. Besides, the approximately equivalent equation of the nonlinear stochastic system is provided with the first-order Taylor approximation of the mass disturbance. Then, according to the obtained approximately equivalent equation, the corresponding stochastic averaging method is given based on the diffusion approach to energy envelope. In this case, the analytic method of the stationary probability density solution is proposed by solving the obtained FPK equation. At last, an example of nonlinear oscillator with variable mass and fractional derivative damping is given to illustrate the effectiveness of the obtained results.