Transient Response Analysis of Nonlinear Oscillators With Fractional Derivative Elements Under Gaussian White Noise Using Complex Fractional Moments

Mechanical and structural systems in various engineering applications are exposed to random excitations, such as seismic ground motion, fluctuating pressure in turbulent flow, and hydraulic force of ocean waves. In order to grasp the response characteristics of dynamic systems subjected to random excitation and to assess their reliability, it is important to accurately determine the probability density function (PDF) of the system response. For a nonlinear system under Gaussian white noise, the system response can usually be regarded as a Markov process, and the time evolution of the response PDF is governed by a partial differential equation called Fokker–Planck (FP) equation. Unfortunately, exact solutions of the FP equation are not available except for the stationary response of a very limited class of nonlinear systems [1]. Solutions for transient/non-stationary responses are even more difficult to obtain, owing to their time dependence. Therefore, the development of approximate analytical or numerical techniques for finding system response PDFs is a significant issue in the field of stochastic dynamics.

In the last few decades, fractional calculus has attracted a great deal of attention due to its usefulness as a mathematical modeling tool for non-local or long-memory phenomena. Various topics on fractional calculus in science and engineering can be found in the books [2–7] and the papers [8–10]. One successful application of fractional calculus to engineering is the modeling of viscoelastic materials. Indeed, many experimental results have proven that fractional derivative models with a small number of parameters can describe viscoelastic behavior with a high degree of accuracy [11–18]. Therefore, in mechanical and structural dynamics, the fractional models have become widely used as models of viscoelastic dampers for seismic isolation and vibration control [16,19–23].

Under such circumstances, there has been a growing interest in the response determination of dynamic systems with fractional derivative elements under stochastic excitation. In this regard, Monte Carlo simulation (MCS) comes to mind as a general-purpose method. However, computing the fractional derivative operator requires considerable computational time and resources, since this operator is defined by a convolution integral, i.e., depends on the function’s entire past time history. Hence, direct numerical integration of system equations with fractional-order derivatives is much more computationally demanding than that of systems with only integer-order derivatives. For this reason, in recent years, a lot of efforts have been put into developing analytical techniques for determining the stochastic response of fractional-order systems. Meanwhile, machines and structures including viscoelastic materials have often been modeled as single-degree-of-freedom (SDOF) systems with fractional derivative elements [2,4,8,9,12–14,16]. Therefore, the development of analytical methods for SDOF systems with fractional derivatives is interesting and important from a practical point of view.

For moment analysis, an approach based on convolutional representation of the response was presented to obtain the response variance of SDOF linear systems with fractional derivatives [24,25]. More recently, approximate closed-form solutions of response moments were obtained for fractional linear oscillators to evolutionary stochastic excitation by expressing their impulse response with Prony series [26]. A moment equation approach has also been employed [27–29], where fractional-order systems are approximately converted into coupled integer-order differential equations. In addition, frequency-domain approaches have been developed for stationary moment analysis [30–32].

For the analysis of PDFs, stochastic averaging has been widely utilized. The displacement and velocity responses of fractional systems are no longer a Markov vector due to the history dependence of the fractional terms, even if the excitation is a white noise process. Nevertheless, by applying stochastic averaging, the response amplitude can be approximately treated as a Markov process, resulting in the FP equation governing the amplitude PDF [33]. As an additional benefit, system dimension reduction is also achieved. In this manner, the amplitude PDF of the stationary response of various fractional nonlinear oscillators has been accurately obtained [33–41].

Compared to moment and stationary PDF analyses, the analysis of transient response PDF is a more difficult task. Although the number of studies is still limited at present, efforts are being made to seek an approximate solution of the transient amplitude PDF of fractional nonlinear systems using a path integral technique [42] or a Rayleigh distribution assumption [43]. However, these existing methods have some drawbacks in terms of computational cost and applicability. The path integration generally requires very short time-steps to obtain a highly accurate solution. The Rayleigh distribution assumption for amplitude PDF may not be appropriate for some types of nonlinear systems, such as oscillators possessing limit cycles. Therefore, it remains an important challenge to develop efficient and versatile analytical methods for determining the transient response PDF of nonlinear oscillators with fractional derivative elements.

Meanwhile, a novel statistical moment named complex fractional moment (CFM) was proposed for the characterization of random variables [44,45]. The CFM is a complex-order moment, which is calculated as the Mellin transform of a PDF. Similar to the Fourier transform relation between PDFs and characteristic functions, the inverse Mellin transform of the CFM reconstructs the original PDF. This striking feature of the CFM allows us to estimate PDFs via CFMs. Di Paola presented a CFM-based approach for solving the FP equation corresponding to first-order nonlinear systems excited by Gaussian white noise [46], in which by applying the Mellin transform to the FP equation, it is converted into a set of linear ordinary differential equations with respect to response CFMs. In other words, a partial differential equation is converted into linear ordinary differential equations much easier to solve. Then, solving the CFM equations and taking the inverse Mellin transform on the resulting CFMs, the response PDF is constructed. Through the above method, the PDFs of both transient and stationary responses were accurately determined in the whole domain, including the tail region. This method has also been applied successfully for white noise excitations of other types [47–49]. Jin et al. developed the CFM-based approach for SDOF oscillators with nonlinear damping under Gaussian white noise by combining stochastic averaging [50]. Stochastic averaging leads to the FP equation for response amplitude, which is then solved via CFMs. The CFM method combined with stochastic averaging yielded transient and stationary response PDFs with high accuracy, and the computational time was significantly reduced compared to the corresponding MCS. In the subsequent studies, this method has been extended to a wide class of integer-order nonlinear oscillators, including systems with nonlinearities in both stiffness and damping [51], nonlinear vibro-impact systems [52] and nonlinear systems with time delay [53]. Such great success of the CFM-based method suggests that it can also be a powerful tool in the transient response analysis of fractional-order nonlinear oscillators.

In this paper, a CFM-based analysis is performed to determine the transient response PDF of nonlinear oscillators with fractional derivative elements, whose nonlinearities are included in both stiffness and damping. This analysis is mainly based on an analytical method for integer-order nonlinear oscillators recently developed by some of the authors [51]. First, we aim to derive the FP equation for response amplitude using stochastic averaging. For this purpose, the equivalent natural frequency and equivalent damping coefficient of the oscillator need to be determined. As an important point in the analysis of fractional oscillators, these equivalent coefficients have to include not only the effects of nonlinearities, but also the contributions of the fractional derivative terms. Furthermore, the coefficients must be given in polynomial form of response amplitude to utilize the CFM-based approach. This paper proposes formulas for determining the equivalent coefficients that satisfy the above requirements, based on an equivalent linearization technique. Then, resorting to stochastic averaging, the FP equation is derived. The solution of the FP equation, i.e., the amplitude PDF, is sought through the CFM approach. The joint PDF of displacement and velocity is also obtained from the amplitude PDF. The accuracy and versatility of the present analytical technique are demonstrated through numerical examples of three linear and nonlinear oscillators with fractional derivatives.

The purpose of this paper is to obtain the transient response PDF of the oscillator expressed by Eq. (7). Some of the authors recently developed a CFM-based analytical approach for determining the transient response PDF of an integer-order nonlinear oscillator obtained by removing the fractional derivative term from Eq. (7) [51]. The present analysis for the fractional nonlinear oscillator (7) is also performed based on this analytical approach. An important point that must be addressed to accomplish this analysis is the appropriate determination of the equivalent natural frequency and equivalent damping coefficient for the fractional nonlinear oscillator.

To solve the FP equation derived in the following section via the CFM-based approach, $ωe2(A)$ and $βe(A)$ must be given as polynomials in amplitude $A$ [51]. Unfortunately, this requirement is not met, due to the presence of $ωeα(A)$ and $ωeα−1(A)$ in Eqs. (16) and (17) (⁠$0<α<1$⁠). In order to overcome this problem, $ωeα(A)$ in Eq. (16) and $ωeα−1(A)$ in Eq. (17) are approximated by $ω0α$ and $ω0α−1$⁠, respectively. This analysis assumes $ε≪1$⁠. In Eq. (16), the term to which the approximation is made includes $ε$⁠. On the other hand, Eq. (17) does not include $ε$⁠, but the damping term of the equivalent linear system in Eq. (11) is multiplied by $ε$⁠. Therefore, the effect of the approximation is considered sufficiently small.

Finally, the magnitude of the parameter $ε$ of the oscillator in performing the above analysis is mentioned. In Refs. [61] and [62], the stationary response PDFs and transient response statistics of various integer-order nonlinear systems were analyzed by combining the equivalent linearization and the stochastic averaging. It was shown that these PDFs and statistics can be obtained correctly when $ε$ is less than $0.1$⁠. Thus, the range of $ε$ for which the combination of the equivalent linearization and the stochastic averaging works well is expected to be $ε<0.1$⁠.

In this way, the solution of the FP equation (22) is attained in terms of CFMs by solving Eqs. (34)–(36). The computation time is only a few seconds. When performing a stationary response analysis, linear algebraic equations derived by setting the time derivative term in Eq. (34) to zero are solved in conjunction with Eq. (35).

The displacement PDF $pX(x,t)$ and the velocity PDF $pX˙(x˙,t)$ are determined as the marginal PDFs of $pXX˙(x,x˙,t)$⁠.

In Secs. 4 and 5, the derivation process for the equations used in the present analysis was described in detail. However, when applying the present approach, it is not necessary to perform all of the analyses in Secs. 4 and 5 in sequence. The actual analysis consists of fewer equations and procedures. The practical analytical procedure is summarized as follows:

1. Calculate the equivalent natural frequency $ωe(A)$ and the equivalent damping coefficient $βe(A)$ using Eqs. (18) and (19). $ωe2(A)$ is given in the form of Eqs. (25) and (29), while $βe(A)$ is given in the form of Eq. (30).

2. Substitute $ωe(A)$ and $βe(A)$ into Eq. (34) to get the linear ordinary differential equations for amplitude CFMs $MpA(γs−1,t)(s=−m,…,−1,1,…,m)$⁠.

3. Solve Eq. (34) in combination with Eq. (35) and the initial condition given by Eq. (36) to determine $MpA(γs−1,t)$$(s=−m,…,m)$⁠.

4. Apply the inverse Mellin transform to $MpA(γs−1,t)$⁠, i.e., substitute $MpA(γs−1,t)$ into Eq. (6) to obtain the amplitude PDF $pA(a,t)$⁠.

5. Calculate $a$ and $|J|$ in the joint PDF $pXX˙(x,x˙,t)$ of $X(t)$ and $X˙(t)$ in Eq. (37) via Eqs. (38) and (39).

6. Substitute $pA(a,t)$⁠, $a,$ and $|J|$ into Eq. (37) to obtain $pXX˙(x,x˙,t)$⁠.

The accuracy and versatility of the present analytical approach are assessed through numerical examples of three linear and nonlinear oscillators with fractional derivative elements. To the best of the author’s knowledge, no exact PDF solution has yet been found for oscillators with fractional derivatives, even for stationary responses. For this reason, the analytical results of response PDFs are compared with the results of the pertinent MCS.

The reliability of the analytical method is validated with three different values of the fractional derivative order $α$⁠: $α=0.3,0.5,$ and $0.7$⁠. It is known that $0.3≤α≤0.7$ is frequently adopted as the best fitting value in modeling the dynamic properties of viscoelastic materials [11–15]. For this reason, the above three values are used for the fractional order. Other system and excitation parameters are set as $ω0=1$⁠, $ε=0.01$⁠, $λ=7,$ and $S0=1$⁠. Further, $Δη=0.3$⁠, $ρ=3.7,$ and $m=220$ are adopted in solving Eq. (44).

For comparisons of the analytical results with simulation results, MCS is performed on Eq. (42). The difficulty inherent in simulations of fractional-order systems is that the fractional derivatives depend on the entire past time history of the system response, as can be seen from Eq. (8). As an efficient and accurate numerical technique, this study utilizes a method to approximately convert an equation of motion including fractional-order derivatives into a set of integer-order differential equations [63]. A brief summary of this methodology is provided in Appendix. Once the integer-order differential equations are obtained, standard numerical integration algorithms can be used to generate realizations of the system response. In this study, the fourth-order Runge–Kutta method is employed. The time increment size $Δt$ for the numerical integration is $Δt=10−4$⁠. MCS results are computed from $5×106$ realizations. $Δt$ and the number of realizations were chosen based on a preliminary study in order to obtain as accurate a solution as possible, since the MCS results are used as substitutes for the true response PDFs of the oscillator. The computer system used in the MCS has an Intel Core i9-10,900 K CPU with 32 GB RAM, and the operating system is Microsoft Windows 10. The source code for the MCS was prepared using the c language, which was compiled with the GNU compiler collection (version 12.3). The simulation method described herein is also adopted in the second and third examples shown later.

In Fig. 1, the PDFs $pA(a,t)$⁠, $pX(x,t)$⁠, and $pX˙(x˙,t)$ of amplitude $A(t)$⁠, displacement $X(t)$⁠, and velocity $X˙(t)$ for the case of $α=0.5$ are shown. The solid lines indicate the analytical results, whereas the discrete symbols denote the pertinent MCS results. The three times in the figure were chosen to show the early phase of the response, the stationary phase and a phase in between the two, while also taking into account the ease of viewing the figure. From the standpoint of evaluating system reliability, it is crucial to obtain the tails of response PDFs with high accuracy. Therefore, semi-log plots of $pX(x,t)$ and $pX˙(x˙,t)$ are also presented to confirm the agreement between the analytical and MCS results in the tail region. The results of $pA(a,t)$ demonstrate the validity of the assumptions about $pA(a,t)$ used in deriving Eq. (32) from Eq. (31). Slight deviations can be seen in the peaks of the displacement and velocity PDFs. However, the accuracy of the analytical results is overall quite satisfactory, including the tail region. Although not shown in the figure, good agreement was also achieved for the cases of $α=0.3$ and $α=0.7$ (some results are provided in Fig. 2). Additionally, the authors also carefully checked the validity for various other values of $α$ in the range $0<α<1$⁠, and the analytical method accurately yielded transient response PDF solutions for any $α$⁠, including the cases of $α=0.1$ and $0.9$⁠.

Figure 2 shows a comparison of displacement and velocity PDFs between the three types of fractional order $α$⁠. The following tendencies can be observed in the response PDFs for larger $α$⁠:

• The widths of the PDFs are narrower (i.e., their variance decreases).

• They reach stationary PDFs faster.

Finally, the computation time for the present analytical method was $2.7$ s, whereas the MCS took about $25$ h (the computation times were similar in the second and third examples shown later as well). Thus, the computational time was significantly reduced compared to the corresponding MCS.