Dynamical analysis of a fractional-order nonlinear two-degree-of-freedom vehicle system by incremental harmonic balance method

Abstract

At present, there are few studies considering both nonlinear and fractional characteristics of suspension in vehicle systems. In this paper, a fractional nonlinear model of a quarter vehicle with two-degree-of-freedom (2-DOF) is innovatively proposed to describe the suspension system containing the viscoelastic material metal rubber. Given the lack of a general calculation scheme for the multi-degree-of-freedom fractional-order incremental harmonic balance method (IHBM), a general calculation scheme for the 2-DOF incremental harmonic balance method for nonlinear systems with fractional order is derived. The nonlinear dynamical properties of the presented model are acquired using this method. The accuracy of the proposed method is verified through a comparison with the power series expansion method. Afterward, the effects of the various parameters on the dynamic performance are analyzed. The vibration peak value of the fractional-order model is significantly higher than that of the integer-order model (IOM). Therefore, the suspension parameters should be designed with a margin when using IOM.

Introduction

Suspension is an important part of a vehicle, affecting its ride comfort, handling stability, safety, and other performances. The accuracy of the suspension model directly determines the dynamic response characteristics of the vehicle. Viscoelastic damping materials play important roles in damping and vibration control of automobile suspension. When modeling viscoelastic such as hydropneumatic suspension, air suspension, and magnetorheological suspension, the problem of the dissipation force equation of viscoelastic damping material in a dynamic system needs to be solved. The dissipative force is not only related to the historical velocity of viscoelastic material but also transformed into a nonlinear eigenvalue problem in the frequency domain. ¹ Both fractal model and fractional-order model can be used to describe viscoelastic materials. The fractal derivative is a local operator. After transforming the fractal derivative, it can be seen that the fractal derivative reflects the memory of the operator through the nonlinear term. However, the fractional-order derivative is a global operator and describes the historical dependence of the system evolution through the convolution integral. ² The two-scale fractal theory is now widely used to model various phenomena arising in porous media or un-smooth boundaries, but it is to be noted that the exact solutions are generally unreachable. ³ The fractional-order model (FOM) has exceptional memory function and can accurately describe the constitutive relationship of viscoelastic materials in the time domain and frequency domain. ⁴ Fractional calculus has promoted the development of nonlinear dynamics. For instance, professor EL-DIB^5–7 used enhanced perturbation analysis and other algorithms to analyze problems in fractional dynamics. Many scholars established the dynamic system model of viscoelastic suspension by using fractional calculus theory. Professor Oustaloup ⁸ proposed for the first time that the damper in the traditional passive suspension should be replaced by a mechanical system with fractional damping, thus establishing the CRONE suspension. Based on the characteristics of multiphase medium mechanics of hydropneumatic suspension, the fractional Bagley–Torvik equation was established in Ref. 9, and the feasibility of the FOM in hydropneumatic suspension system modeling was verified by simulation and experiment. Based on the theory of fluid mechanics and the ideal gas equation of state, and the advantages of the fractional viscoelastic material model, the multi-branch hydropneumatic suspension system was modeled fractionally.^10,11 The feasibility and effectiveness of the hydropneumatic suspension model are further verified based on the fractional order after comparing the fractional-order and integer-order simulation curves with the actual experimental curves.

The FOM is more accurate than the integer-order model (IOM) in studying the suspension dynamic behavior, which is worthy of further study.

Nonlinearity is a crucial feature of an automobile vertical vibration system. Many scholars have analyzed the dynamic characteristics of a two-degree-of-freedom (2-DOF) 1/4 nonlinear vehicle model.^12–16 Results with high accuracy can be obtained when analyzing nonlinear motion with multi-scale techniques.^17–20 In terms of quantitative research, due to the strong nonlinearity and strong coupling of automobile suspension models, methods such as the perturbation method, average method, and multi-scale method that treat the system as having weak nonlinearity are no longer applicable. Giver this problem, a large number of scholars have done related research to put forward some methods such as the homotopy perturbation method, ²¹ incremental harmonic balance method (IHBM), He’s frequency formulation, ²² and so on. The incremental harmonic balance method can study weak nonlinearity and strong nonlinearity at the same time, and the convergence accuracy can be flexibly controlled. Thus, this method has become an important one for nonlinear system research. Sheng and Wu ²³ used IHBM to analyze the dynamic behavior of the 2-DOF vehicle model and obtained the steady-state periodic solution of a nonlinear vehicle system. Zhou et al. ²⁴ used IHBM to study the quadratic and cubic stiffness, and damping is contained in the suspension system and tire meantime. Zhou et al.^12,25 used IHBM and the Newmark-β method to compare the dynamic response of a nonlinear quarter car, and IHBM was used to study the piecewise leaf spring for the rear suspension of a truck. But the classical IHB method does not consider the fractional-order nonlinear terms to investigate the 2-DOF system. This paper tackles this problem by deriving the incremental term of the fractional-order nonlinear term of the 2-DOF fractional nonlinear system. It is still rare to find the dynamics of 2-DOF systems in nonlinear automotive systems with fractional differential terms. Therefore, IHBM is adopted in this paper to study the 2-DOF nonlinear vehicle suspension model with fractional order to obtain its frequency response characteristics and analyze its dynamic characteristics and the influence of various parameters on the dynamic performance.

In this paper, the FOM of a vehicle system is established, and the common format of 2-DOF IHBM of a fractional system is derived. Then, the amplitude–frequency response characteristics of the FOM 2-DOF vehicle system are analyzed.

Model of the fractional-order nonlinear two-degree-of-freedom vehicle system

According to the literature, ²⁶ metal rubber is a viscoelastic material, lying between the ideal elastomer and the ideal cohesive body, and its material mechanical properties are time-dependent. The elastic restoring force, the damping force with memory, and the complex damping force are represented by the equivalent viscoelastic damping element, and the fractional-order differential term is introduced into the constitutive relation of metal rubber. A simple nonlinear dynamic model of the metal rubber system is obtained and verified by experiments. Experimental results show that this model can accurately describe the mechanical properties of metal rubber with fewer parameters and is suitable for a wide frequency range. Its constitutive model is expressed as equation (1).

$g\left(x\left(t\right),\dot{x}\left(t\right),t\right)=k_{1}x\left(t\right)+k_3{x}^3\left(t\right)+c\dot{x}\left(t\right)+h\left[D^{p}x\left(t\right)\right]$

(1)

where $x\left(t\right)$ is the deformation of metal rubber spring, $k_1$ is the constant partial stiffness coefficient of linear elastic restoring force, $k_3$ is the stiffness coefficient of cubic nonlinear elastic restoring force, $c$ is the first-order viscous damping coefficient, h is the fractional viscoelastic damping coefficient, and $p$ is the fractional derivative order.

Many different definitions of fractional calculus have been given. The definitions of Caputo and Riemann–Liouville fractional differential are shown in equations (2) and (3), respectively. Caputo-type fractional calculus is used in this paper.

${\hspace{0.17em}}_{t_0}^C{D}_t^{p}x\left(t\right)=\frac{1}{\Gamma \left(1-p\right)}{\int }_{t_0}^t \frac{x^{\prime }\left(\tau \right)}{{\left(t-\tau \right)}^p}d\tau$

(2)

${\hspace{0.17em}}_{t_0}^{RL}{D}_t^{p}f\left(t\right)=\frac{1}{\Gamma \left(n-p\right)}\frac{d^n}{d{t}^n}{\int }_{t_0}^t \frac{f\left(\tau \right)}{{\left(t-\tau \right)}^{1+p-n}}d\tau$

(3)

where $\Gamma \left(x\right)$ is a gamma function, and $\Gamma \left(x+1\right)=x\Gamma \left(x\right)$ .

In this paper, metal rubber is used as the suspension material to study the 1/4 automobile suspension model. Figure 1 shows the 1/4 automobile suspension system with the metal rubber substituting the oil, and a fractional nonlinear 2-DOF 1/4 automobile model is constructed according to the system, as shown in Figure 2.

Figure 1.

1/4 physical model of automobile suspension including metal rubber.

Figure 2.

Fractional-order nonlinear 2-DOF suspension model.

In this figure, $m_1$ and $m_2$ are the mass of the body and tire, respectively; $z_0=F\hspace{0.17em}\cos \left(\omega t^{\prime }\right)$ is the road excitation input; and $\mathrm{\omega }$ is the excitation frequency. $z_1$ and $z_2$ are the tire displacement and body displacement, respectively, and k _t is the tire stiffness. The suspension system material is a viscoelastic material, so the system includes nonlinear elastic restoring force and linear damping restoring force and the viscous restoring force expressed by fractional-order terms as follows: ${hD}_t^p\left(z_1-z_2\right)$ , $h$ is the coefficient of fractional differential term and $D_t^p\left(z_1-z_2\right)$ is the $p\left(0\le p\le 1\right)$ derivative of $\left(z_1-z_2\right)$ with respect to $t$ . Caputo fractional-order definition is adopted in this paper.

With the higher-order nonlinear elastic restoring force and the higher-order nonlinear damping term being ignored, the expression of the nonlinear elastic restoring force is $F_k\left(z_1-z_2\right)=k_1\left(z_1-z_2\right)+k_3{\left(z_1-z_2\right)}^3$ . k₁ and k₃ are the nonlinear spring force coefficients of suspension. The linear damping force is $c\left({\dot{z}}_1-{\dot{z}}_2\right)$ , and $c$ is the damping coefficient of suspension.

According to Newton’s second law and the constitutive model of metal rubber, the motion differential equation of the suspension system can be obtained, as shown in equation (4).

$\{\begin{array}{l}{m}_1{\ddot{z}}_1+k_1\left(z_1-z_2\right)+k_3{\left(z_1-z_2\right)}^3+c\left({\dot{z}}_1-{\dot{z}}_2\right)+h{D}_t^p\left(z_1-z_2\right)=0\\ m_2{\ddot{z}}_2+k_t\left(z_2-z_0\right)-k_1\left(z_1-z_2\right)-k_3{\left(z_1-z_2\right)}^3-c\left({\dot{z}}_1-{\dot{z}}_2\right)-h{D}_t^p\left(z_1-z_2\right)=0\end{array}$

(4)

Let $x_1=z_1,y_1=z_1-z_2$ . Then, formula (4) can be written as formula (5).

$\{\begin{array}{l}{m}_1{\ddot{x}}_1+k_1{y}_1+k_3{y}_1^3+c{\dot{y}}_1+h{D}_t^p{y}_1=0\\ m_2\left({\ddot{x}}_1-{\ddot{y}}_1\right)+k_t{x}_1-k_t{y}_1-k_1{y}_1-k_3{y}_1^3-c{\dot{y}}_1-h{D}_t^p{y}_1=k_{t}F\hspace{0.17em}\cos \left(\mathrm{\omega }{t}^{\prime }\right)\end{array}$

(5)

Equation (5) is normalized, and let

$x=\frac{x_1}{k_{t}F},y=\frac{y_1}{k_{t}F},{\omega }_1=\sqrt{\frac{k_1}{m_1}},\nu =\frac{h}{m_1}{,\omega }_2=\sqrt{\frac{k_t}{m_2}},\mu =\frac{m_1}{m_2},\delta =\frac{1}{m_2},\gamma =\frac{k_3{\left(k_{t}F\right)}^2}{m_1},\xi =\frac{c}{2m_1{\omega }_2},$

$\mathrm{\Omega }=\frac{\omega }{{\omega }_1},t={\omega }_1{t}^{\prime }$

Then, (6) is deduced as

$\{\begin{array}{l}\ddot{x}+{\omega }_1^{2}y+\gamma y^3+2\xi {\omega }_2\dot{y}+\nu D_t^{p}y=0\\ \ddot{x}-\ddot{y}+{\omega }_2^{2}x-{\omega }_2^{2}y-{\omega }_1^2\mu y-\gamma \mu y^3-2\xi {\omega }_2\mu \dot{y}-\nu \mu D_t^{p}y=\delta \hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\end{array}$

(6)

in the form of a matrix

$\overline{\mathit{M}}\ddot{\mathit{X}}+\overline{\mathit{C}}\dot{\mathit{X}}+\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3\right)\mathit{X}+\overline{\mathit{H}}{D}_t^p\mathit{X}=\overline{\mathit{P}}$

(7)

where

$\mathit{X}={\left[\begin{array}{ll}x& y\end{array}\right]}^{\mathrm{T}};\overline{\mathit{M}}=\left[\begin{array}{ll}1& 0\\ 1& -1\end{array}\right];\overline{\mathit{C}}=\left[\begin{array}{ll}0& 2\xi {\omega }_2\\ 0& -2\xi {\omega }_2\mu \end{array}\right];{\overline{\mathit{K}}}_1=\left[\begin{array}{ll}0& {\omega }_1^2\\ {\omega }_2^2& -{\omega }_2^2-{\omega }_1^2\mu \end{array}\right];{\overline{\mathit{K}}}_3=\left[\begin{array}{ll}0& \gamma y^2\\ 0& -\gamma \mu y^2\end{array}\right];$

$\overline{\mathit{H}}=\left[\begin{array}{ll}0& \nu \\ 0& -\nu \mu \end{array}\right];\overline{\mathit{P}}=\left[\begin{array}{l}0\\ \delta \hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\end{array}\right]$

Let

$\tau =\mathrm{\Omega }t,{\mathit{X}}^{\prime }=\frac{d\mathit{X}}{d\tau },{\mathit{X}}^{″}=\frac{d^2\mathit{X}}{d{\tau }^2}$

We have

$\{\begin{array}{l}\dot{\mathit{X}}=\frac{d\mathit{X}}{dt}=\mathrm{\Omega }{\mathit{X}}^{\prime }\\ \ddot{\mathit{X}}=\frac{d^2\mathit{X}}{d{t}^2}={\mathrm{\Omega }}^2{\mathit{X}}^{″}\\ D_t^p\mathit{X}={\mathrm{\Omega }}^p{D}_{\tau }^p\mathit{X}\end{array}$

(8)

Formula (9) can be obtained by substituting formula (8) into formula (7).

${\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{X}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{X}}^{\prime }+\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3\right)\mathit{X}+{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathit{X}=\overline{\mathit{P}}$

(9)

Incremental harmonic balance method computation scheme for fractional-order nonlinear two-degree-of-freedom

Incremental process

Suppose X ₀ is the solution of the system (9) at a certain moment under the road surface excitation of fixed frequency ω, and its neighboring state can be expressed as an increment.

$\mathit{X}={\mathit{X}}_0+\mathrm{\Delta }\mathit{X}$

(10)

where

${\mathit{X}}_0={\left[\begin{array}{ll}{x}_0& y_0\end{array}\right]}^{\mathrm{T}},\mathrm{\Delta }\mathit{X}={\left[\begin{array}{ll}\mathrm{\Delta }x& \mathrm{\Delta }y\end{array}\right]}^{\mathrm{T}}$

Formula (11) can be obtained by substituting formula (10) into formula (9) through Taylor expansion and by omitting the higher-order small quantities.

${\mathrm{\Omega }}^2\overline{\mathit{M}}\mathrm{\Delta }{\mathit{X}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}\mathrm{\Delta }{\mathit{X}}^{\prime }+\left({\overline{\mathit{K}}}_1+3{\overline{\mathit{K}}}_3^{\left(1\right)}\right)\mathrm{\Delta }\mathit{X}+{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathrm{\Delta }\mathit{X}=\overline{\mathit{P}}-\left\{{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{X}}_0^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{X}}_0^{\prime }+\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3^{\left(1\right)}\right){\mathit{X}}_0+{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p{\mathit{X}}_0\right\}$

(11)

where

${\overline{\mathit{K}}}_3^{\left(1\right)}=\left[\begin{array}{ll}0& y_0^2\\ 0& -k_3{y}_0^2\end{array}\right]$

Harmonic balance process

The exact solution ${\mathit{X}}_0$ and the increment $\mathrm{\Delta }\mathit{X}$ Fourier expansion are shown in (12).

$\{\begin{array}{l}{x}_{j0}=a_{j0}+{\displaystyle \sum }_{k=1}^N\left(a_{jk}\cos \hspace{0.17em}k\tau +b_{jk}\sin \hspace{0.17em}k\tau \right)\\ \mathrm{\Delta }{x}_j=\mathrm{\Delta }{a}_{j0}+{\displaystyle \sum }_{k=1}^N\left(\mathrm{\Delta }{a}_{jk}\cos \hspace{0.17em}k\tau +\mathrm{\Delta }{b}_{jk}\sin \hspace{0.17em}k\tau \right)\\ j=\mathrm{1,2}\hspace{0.17em}k=1,...,N\end{array}$

(12)

Definition:

$\begin{array}{l}{\mathit{A}}_j={\left[\begin{array}{ll}\begin{array}{lll}{a}_{j0},& a_{j1},& a_{j2}\end{array},\begin{array}{ll}...,& a_{jN}\end{array},& b_{j1},\begin{array}{lll}{b}_{j2},& ...,& b_{j\mathrm{N}}\end{array}\end{array}\right]}^{\mathrm{T}}\\ \mathrm{\Delta }{\mathit{A}}_j={\left[\begin{array}{ll}\begin{array}{lll}\mathrm{\Delta }{a}_{j0},& \mathrm{\Delta }{a}_{j1},& \mathrm{\Delta }{a}_{j2},\end{array}\begin{array}{ll}...,& \mathrm{\Delta }{a}_{jN},\end{array}& \mathrm{\Delta }{b}_{j1},\mathrm{\Delta }{b}_{j2},\end{array}...,\mathrm{\Delta }{b}_{jN}\right]}^{\mathrm{T}}\\ {\mathit{C}}_S={\left[\begin{array}{ll}\begin{array}{lll}1,& \cos \hspace{0.17em}\tau ,& \cos \hspace{0.17em}2\tau ,\end{array}\begin{array}{ll}...,& \cos \hspace{0.17em}N\tau \end{array},& \sin \hspace{0.17em}\tau ,\end{array}\sin 2\tau ,...,\sin N\tau \right]}^{\mathrm{T}}\end{array}$

Let

$\mathit{S}=\left[\begin{array}{ll}{\mathit{C}}_S& 0\\ 0& {\mathit{C}}_S\end{array}\right];\mathit{A}=\left[\begin{array}{l}{\mathit{A}}_1\\ {\mathit{A}}_2\end{array}\right];\mathrm{\Delta }\mathit{A}=\left[\begin{array}{l}\mathrm{\Delta }{\mathit{A}}_1\\ \mathrm{\Delta }{\mathit{A}}_2\end{array}\right]$

Then, the exact solution and increment can be expressed as formula (13).

$\{\begin{array}{l}{\mathit{X}}_0=\mathit{S}\mathit{A}\\ \mathrm{\Delta }\mathit{X}=\mathit{S}\mathrm{\Delta }\mathit{A}\end{array}$

(13)

When the Galerkin average process is applied to integrate formula (11), the fractional differential term is an aperiodic function. Thus, the fractional differential terms should be integrated and averaged within the period T (T = ∞), and the period of other functions is $2\pi$ . Therefore, through integration in $2\pi$ , equation (14) can be obtained.

$\begin{array}{l}\frac{1}{2\pi }{\int }_0^{2\pi }{\left(\delta \mathrm{\Delta }\mathit{X}\right)}^{\mathrm{T}}\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}\mathrm{\Delta }{\mathit{X}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}\mathrm{\Delta }{\mathit{X}}^{\prime }+\left({\overline{\mathit{K}}}_1+3{\overline{\mathit{K}}}_3^{\left(1\right)}\right)\mathrm{\Delta }\mathit{X}\right]d\tau +\frac{1}{T}{\int }_0^T{\left(\delta \mathrm{\Delta }\mathit{X}\right)}^{\mathrm{T}}{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathrm{\Delta }\mathit{X}d\tau \\ =\frac{1}{2\pi }{\int }_0^{2\pi }{\left(\delta \mathrm{\Delta }\mathit{X}\right)}^{\mathrm{T}}\left\{\overline{\mathit{P}}-\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{X}}_0^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{X}}_0^{\prime }+\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3^{\left(1\right)}\right){\mathit{X}}_0\right]\right\}d\tau -\frac{1}{T}{\int }_0^T{\left(\delta \mathrm{\Delta }\mathit{X}\right)}^{\mathrm{T}}{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p{\mathit{X}}_{0}d\tau \end{array}$

(14)

Equation (13) is substituted into integral equation (14) to obtain equation (15).

$\begin{array}{l}{\left(\delta \mathrm{\Delta }\mathit{A}\right)}^{\mathrm{T}}\left\{\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{S}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{S}}^{\prime }+({\overline{\mathit{K}}}_1+3{\overline{\mathit{K}}}_3^{\left(1\right)}\right]d\tau \right\}\mathrm{\Delta }\mathit{A}+{\left(\delta \mathrm{\Delta }\mathit{A}\right)}^{\mathrm{T}}\left(\frac{1}{T}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau \right)\mathrm{\Delta }\mathit{A}\\ ={\left(\delta \mathrm{\Delta }\mathit{A}\right)}^{\mathrm{T}}\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left\{\overline{\mathit{P}}-\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{S}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{S}}^{\prime }+\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3^{\left(1\right)}\right)\mathit{S}\right]\mathit{A}\right\}d\tau -{\left(\delta \mathrm{\Delta }\mathit{A}\right)}^{\mathrm{T}}\left(\frac{1}{T}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau \right)\mathit{A}\end{array}$

(15)

Furthermore, formula (16) is obtained.

$\begin{array}{l}\left\{\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{S}}^{″}+\omega \overline{\mathit{C}}{\mathit{S}}^{\prime }+({\overline{\mathit{K}}}_1+3{\overline{\mathit{K}}}_3^{\left(1\right)}\right]d\tau +\frac{1}{T}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}{\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau \right\}\mathrm{\Delta }\mathit{A}\\ =\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left\{\overline{\mathit{P}}-\left[{\mathrm{\Omega }}^2\overline{\mathit{M}}{\mathit{S}}^{″}+\mathrm{\Omega }\overline{\mathit{C}}{\mathit{S}}^{\prime }+({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3^{\left(1\right)}\right]\mathit{A}\right\}d\tau -\left[\frac{1}{T}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}\left({\mathrm{\Omega }}^p\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}\right)\mathit{A}d\tau \right]\end{array}$

(16)

A linear group of equations of $∆\mathit{A}$ is obtained

$\mathit{N}\mathrm{\Delta }\mathit{A}=\mathit{R}$

(17)

where

$\mathit{N}={\mathit{N}}_1+{\mathit{N}}_2^p$

${\mathit{N}}_1=\frac{1}{2\pi }\left\{{\mathrm{\Omega }}^2{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\overline{\mathit{M}}{\mathit{S}}^{″}d\tau +\mathrm{\Omega }{\int }_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\overline{\mathit{C}}{\mathit{S}}^{\prime }d\tau +{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left({\overline{\mathit{K}}}_1+3{\overline{\mathit{K}}}_3^{\left(1\right)}\right)d\tau \right\}$

${\mathit{N}}_2^{\mathit{p}}=\frac{{\mathrm{\Omega }}^p}{\mathrm{T}}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau$

$\mathit{R}={\mathit{R}}_1+{\mathit{R}}_2^p$

${\mathit{R}}_1=\frac{1}{2\pi }\left\{{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\overline{\mathit{P}}d\tau -\left[{\mathrm{\Omega }}^2{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\overline{\mathit{M}}{\mathit{S}}^{″}d\tau +\mathrm{\Omega }{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\overline{\mathit{C}}{\mathit{S}}^{\prime }d\tau +{{\displaystyle \int }}_0^{2\pi }{\mathit{S}}^{\mathrm{T}}\left({\overline{\mathit{K}}}_1+{\overline{\mathit{K}}}_3^{\left(1\right)}\right)\mathit{S}d\tau \right]\times \mathit{A}\right\}$

${\mathit{R}}_2^p=-\frac{{\mathrm{\Omega }}^p}{T}{{\displaystyle \int }}_0^T{\mathit{S}}^{\mathrm{T}}\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau \times \mathit{A}$

${\mathit{N}}_1$ and ${\mathit{R}}_1$ are iteration formats of integer-order differential terms in IHBM, which can be calculated by directly substituting parameters into the calculation, while ${\mathit{N}}_2^p$ and ${\mathit{R}}_2^p$ contain fractional-order differential terms and their corresponding parts in IHBM.

First, ${\mathit{N}}_2^p$ is simplified as formula (18).

${\mathit{N}}_2^P=\frac{{\mathrm{\Omega }}^p}{T}{\int }_0^T{\mathit{S}}^{\mathrm{T}}\overline{\mathit{H}}{D}_{\tau }^p\mathit{S}d\tau$

(18)

$=\frac{{\mathrm{\Omega }}^p}{T}{\int }_0^T\left[\begin{array}{ll}{\mathit{C}}_S^{\mathrm{T}}& 0\\ 0& {\mathit{C}}_S^{\mathrm{T}}\end{array}\right]\left[\begin{array}{ll}0& \nu \\ 0& -\nu \mu \end{array}\right]D_{\tau }^p\left[\begin{array}{ll}{\mathit{C}}_S& 0\\ 0& {\mathit{C}}_S\end{array}\right]d\tau$

$=\frac{{\mathrm{\Omega }}^p}{T}{\int }_0^T\left[\begin{array}{ll}0& \nu {\mathit{C}}_S^{\mathrm{T}}\\ 0& -\nu \mu {\mathit{C}}_S^{\mathrm{T}}\end{array}\right]\left[\begin{array}{ll}{D}_{\tau }^p{\mathit{C}}_S& 0\\ 0& D_{\tau }^p{\mathit{C}}_S\end{array}\right]d\tau$

$=\frac{{\mathrm{\Omega }}^p}{T}{\int }_0^T\left[\begin{array}{ll}0& \nu {\mathit{C}}_S^{\mathrm{T}}{D}_{\tau }^p{\mathit{C}}_S\\ 0& -\nu \mu {\mathit{C}}_S^{\mathrm{T}}{D}_{\tau }^p{\mathit{C}}_S\end{array}\right]d\tau$

$=\left[\begin{array}{ll}0& {\mathit{M}}_2^p\\ 0& -\mu {\mathit{M}}_2^p\end{array}\right]$

${\mathit{M}}_2^p$ can be written as (19).

${\mathit{M}}_2^P=\frac{{\mathrm{\Omega }}^p}{T}{\int }_0^T\nu {\mathit{C}}_S^{\mathrm{T}}{D}_{\tau }^p{\mathit{C}}_{\mathrm{S}}d\tau =\left[\begin{array}{ll}{\left[{\mathit{M}}_{11}\right]}^P& {\left[{\mathit{M}}_{12}\right]}^P\\ {\left[{\mathit{M}}_{21}\right]}^P& {\left[{\mathit{M}}_{22}\right]}^P\end{array}\right]$

(19)

In formula

${\left[{\mathit{M}}_{11}\right]}_{ij}^p=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \hspace{0.17em}i\tau D_{\tau }^p\cos \hspace{0.17em}j\tau d\tau \hspace{0.17em}i=\mathrm{0,1},...,N\hspace{1em}j=\mathrm{0,1},...,N$

(20a)

${\left[{\mathit{M}}_{12}\right]}_{ij}^p=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \hspace{0.17em}i\tau D_{\tau }^p\sin \hspace{0.17em}j\tau d\tau \hspace{0.17em}i=\mathrm{0,1},...,N\hspace{1em}j=1,...,N$

(20b)

${\left[{\mathit{M}}_{21}\right]}_{ij}^p=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\nu {\mathrm{\Omega }}^p\hspace{0.17em}\sin \hspace{0.17em}i\tau D_{\tau }^p\cos \hspace{0.17em}j\tau d\tau \hspace{0.17em}i=1,...,N\hspace{1em}j=\mathrm{0,1},...,N$

(20c)

${\left[{\mathit{M}}_{22}\right]}_{ij}^p=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\nu {\mathrm{\Omega }}^p\hspace{0.17em}\sin \hspace{0.17em}i\tau D_{\tau }^p\sin \hspace{0.17em}j\tau d\tau \hspace{0.17em}i=1,...,N\hspace{1em}j=1,...,N$

(20d)

According to the definition of Caputo-type calculus and the power series expansion method, formula (21) can be calculated according to Ref. 27.

${\left[{\mathit{M}}_{11}\right]}_{ij}^p=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \hspace{0.17em}i\tau D_{\tau }^p\cos \hspace{0.17em}j\tau \mathrm{d}\tau =\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left\{\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \hspace{0.17em}i\tau \times \left[\frac{1}{\Gamma \left(1-p\right)}{\int }_0^{\tau }\frac{-j\hspace{0.17em}\sin \hspace{0.17em}ju}{{\left(\tau -u\right)}^p}du\right]\right\}d\tau$

(21)

Formula (21) becomes (22) by introducing $s=\tau -u$ and $du=-ds$ .

${\left[{\mathit{M}}_{11}\right]}_{ij}^p=\frac{-j\nu {\mathrm{\Omega }}^p}{\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left\{\cos \hspace{0.17em}i\tau \times \left[{\int }_0^{\tau }\frac{\sin \hspace{0.17em}\left(j\tau -js\right)}{s^p}ds\right]\right\}d\tau$

(22)

$=\frac{-j\nu {\mathrm{\Omega }}^p}{\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left(\cos \hspace{0.17em}i\tau \hspace{0.17em}\sin \hspace{0.17em}j\tau {\int }_0^{\tau }\frac{\cos \hspace{0.17em}js}{s^p}ds\right)d\tau +\frac{-j\nu {\mathrm{\Omega }}^p}{\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left(\cos \hspace{0.17em}i\tau \hspace{0.17em}\cos \hspace{0.17em}j\tau {\int }_0^{\tau }\frac{\sin \hspace{0.17em}js}{s^p}ds\right)d\tau$

For definitions (22), the first part is A₁ and the second is A₂. Formulas (23) and (24) can be obtained by integration in parts.

$A_1=\frac{-j\nu {\mathrm{\Omega }}^p}{\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}\left\{{\int }_0^T\left\{\frac{1}{2}\left[\sin \left(i+j\right)\tau -\sin \left(i-j\right)\tau \right]{\int }_0^{\tau }\frac{\cos \hspace{0.17em}js}{s^p}ds\right\}d\tau \right\}$

(23)

$=\frac{-j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}\left[\frac{\cos \left(i+j\right)\tau }{\left(i+j\right)}-\frac{\cos \left(i-j\right)\tau }{\left(i-j\right)}\right]{{\displaystyle \int }}_0^{\tau }\frac{\cos \hspace{0.17em}js}{s^p}ds{|}_0^T-\frac{j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{{\displaystyle \int }}_0^{\tau }\left[\frac{\cos \left(i+j\right)\tau }{\left(i+j\right)}-\frac{\cos \left(i-j\right)\tau }{\left(i-j\right)}\right]\frac{\cos \hspace{0.17em}j\tau }{{\tau }^p}d\tau$

$A_2=\frac{j\nu {\mathrm{\Omega }}^p}{\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}\left\{{\int }_0^T\left\{\frac{1}{2}\left[\cos \left(i+j\right)\tau +\cos \left(i-j\right)\tau \right]{\int }_0^{\tau }\frac{\sin \hspace{0.17em}js}{s^p}ds\right\}d\tau \right\}$

(24)

$=\frac{j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left[\frac{\sin \left(i+j\right)\tau }{\left(i+j\right)}+\frac{\sin \left(i-j\right)\tau }{\left(i-j\right)}\right]{\int }_0^{\tau }\frac{\sin \hspace{0.17em}js}{s^p}ds{|}_0^T-\frac{j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left[\frac{\sin \left(i+j\right)\tau }{\left(i+j\right)}+\frac{\sin \left(i-j\right)\tau }{\left(i-j\right)}\right]\frac{\sin \hspace{0.17em}j\tau }{{\tau }^p}d\tau$

Let the first part of A₁ be B_A11, the second part be B_A12, the first part of A₂ be B_A21, and the second part be B_A22.

According to the two basic formulas (25)^28–30

$\underset{T\to \infty }{\lim }{\int }_0^T\frac{\sin \hspace{0.17em}\left(js\right)}{s^p}ds=j^{p-1}\Gamma \left(1-p\right)\cos \left(\frac{p\pi }{2}\right)$

(25a)

$\underset{T\to \infty }{\lim }{\int }_0^T\frac{\cos \hspace{0.17em}\left(js\right)}{s^p}dt=j^{p-1}\Gamma \left(1-p\right)\sin \left(\frac{p\pi }{2}\right)$

(25b)

Formula (25a) and (25b) is substituted into (23), (24), and then (26a), (26b), (27a), and (27b) are obtained.

$B_{A11}=\frac{\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2}\underset{T\to \infty }{\lim }\frac{1}{T}\times \left[\frac{\cos \left(i+j\right)T}{\left(i+j\right)}-\frac{\cos \left(i-j\right)T}{\left(i-j\right)}\right]=0$

(26a)

$B_{A12}=-\frac{j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left[\frac{\cos \left(i+j\right)\tau \hspace{0.17em}\cos \hspace{0.17em}j\tau }{\left(i-j\right){\tau }^p}-\frac{\cos \left(i-j\right)\tau \hspace{0.17em}\cos \hspace{0.17em}j\tau }{\left(i-j\right){\tau }^p}\right]d\tau =0$

(26b)

$B_{A21}=\frac{\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)}{2}\underset{T\to \infty }{\lim }\frac{1}{T}\times \left[\frac{\sin \left(i+j\right)T}{\left(i+j\right)}+\frac{\sin \left(i-j\right)T}{\left(i-j\right)}\right]=\{\begin{array}{l}0,i\ne j\\ \frac{\nu \mathrm{\Omega }{j}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2},i=j\end{array}$

(27a)

$B_{A22}=-\frac{j\nu {\mathrm{\Omega }}^p}{2\Gamma \left(1-p\right)}\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\left[\frac{\sin \left(i+j\right)\tau \hspace{0.17em}\sin \hspace{0.17em}j\tau }{\left(i+j\right){\tau }^p}+\frac{\sin \left(i-j\right)\tau \hspace{0.17em}\sin \hspace{0.17em}j\tau }{\left(i-j\right){\tau }^p}\right]d\tau =0$

(27b)

The simultaneous equations (22) to (27) can be used to calculate M ₁₁. Similarly, M ₁₂, M ₂₁, and M ₂₂ can be obtained, as shown in formula (28).

${\left[{\mathit{M}}_{11}\right]}_{ij}^p=\{\begin{array}{ll}0,& i\ne j\\ \frac{\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)}{2},& i=j\end{array}$

(28a)

${\left[{\mathit{M}}_{12}\right]}_{ij}^p=\{\begin{array}{ll}0,& i\ne j\\ \frac{\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2},& i=j\end{array}$

(28b)

${\left[{\mathit{M}}_{21}\right]}_{ij}^p=\{\begin{array}{ll}0,& i\ne j\\ \frac{-\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2},& i=j\end{array}$

(28c)

${\left[{\mathit{M}}_{22}\right]}_{ij}^p=\{\begin{array}{ll}0,& i\ne j\\ \frac{\nu {\mathrm{\Omega }}^p{j}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)}{2},& i=j\end{array}$

(28d)

Expansion (28a)–(28d) gives the form (29a)–(29d).

${\left[{\mathit{M}}_{11}\right]}^p=\frac{\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)}{2}{\left[\begin{array}{llll}{0}^p& 0& \cdots & 0\\ 0& 1^p& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & N^p\end{array}\right]}_{\left(N+1)\times (N+1\right)}$

(29a)

${\left[{\mathit{M}}_{12}\right]}^p=\frac{\nu {\mathrm{\Omega }}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2}{\left[\begin{array}{llll}0& 0& \cdots & 0\\ 1^p& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & N^p\end{array}\right]}_{\left(N+1\right)\times N}$

(29b)

${\left[{\mathit{M}}_{21}\right]}^p=\frac{\nu {\mathrm{\Omega }}^p\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)}{2}{\left[\begin{array}{llll}0& -1^p& \cdots & 0\\ 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& -N^p\end{array}\right]}_{N\times \left(N+1\right)}$

(29c)

${\left[{\mathit{M}}_{22}\right]}^p=\frac{\nu {\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)}{2}{\left[\begin{array}{llll}{1}^p& 0& \cdots & 0\\ 0& 2^p& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& N^p\end{array}\right]}_{N\times N}$

(29d)

The iteration formula of ${\mathit{N}}_2^p$ is obtained from the above equation. Similarly, the iteration formula of ${\mathit{R}}_2^p$ is obtained.

Thus, according to equation (16), when the initial value ${\mathit{A}}_0$ is given, loop iteration is performed until the precision requirements are met. Therefore, the periodic approximate solution of arbitrary precision of the system can be obtained (Figure 3 and Table 1).

Figure 3.

Incremental harmonic balance calculation process.

Table 1.

Simulation parameters.

Suspension parameters	Numerical value
m1/kg	240
m2/kg	30
k1/(N·m−1)	16,000
kt/(N·m−1)	160,000
k3/(N·m−1)	10,000
c/(N·s·m−1)	1000
h	1000
p	0.5

Stability analysis

Most of the discussion about the stability of periodic solutions of IHBM is based on Floquet theory. ³¹ But the fractional-order part in this paper has difficulties when using this theory. A general method is applied to approximate the fractional-order part with sufficient precision in this paper.

In order to facilitate analysis, the fractional terms in equation (5) are represented by equivalent damping coefficient $c_{eq}$ and equivalent stiffness coefficient $k_{eq}$ as follows:

${hD}_t^{p}y=c_{eq}\dot{y}+k_{eq}y$

(30)

where

$\{\begin{array}{l}{c}_{eq}={h\mathrm{\Omega }}^{p-1}\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)\\ k_{eq}={h\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)\end{array}$

(31)

The derivation of equation (31) is described in Ref. 28 and Ref. 32.

This means that the Duffing oscillator with fractional-order derivative which is described by Caputo’s definition can be equivalent by the equivalent damping coefficient and the equivalent stiffness coefficient in equation (31). It can be observed from equation (31) that the equivalent damping coefficient $c_{eq}$ and the equivalent stiffness coefficient $k_{eq}$ are all affected by the fractional coefficient h and the fractional order p. Furthermore it could be concluded that, if the fractional order $p\to 1$ the fractional-order derivative would be almost changed into viscous damping, on the other hand, if the fractional order $p\to 0$ the fractional-order derivative would be almost changed into linear stiffness. This result is the same as the intuitionistic consideration. ²⁸

Substitute equation (31) into (30):

${hD}_t^{p}y={h\mathrm{\Omega }}^{p-1}\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)\dot{y}+{h\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)y$

(32)

In order to obtain the steady-state solutions of $x_{10}$ and $x_{10}$ by IHBM, the small perturbations $x_1$ and $x_2$ are introduced as follows:

$x_1=x_{10}+∆x_1,x_2=x_{20}+∆x_2$

(33)

Substituting $x_1$ and $x_2$ into equation (5), and then we can sort out an equation in matrix form:

$\dot{\mathit{Z}}=\mathit{Q}\left(\tau \right)$

(34)

where

$\mathit{Z}={\left[∆x_1∆{\dot{x}}_1∆x_2∆{\dot{x}}_2\right]}^{\mathrm{T}}$

$\dot{\mathit{Z}}={\left[∆{\dot{x}}_1∆{\ddot{x}}_1∆{\dot{x}}_2∆{\ddot{x}}_2\right]}^{\mathrm{T}}$

$\mathit{Q}\left(\tau \right)=\left[\begin{array}{l}\begin{array}{ll}0& 1\\ 0& 0\end{array}\hspace{0.17em}\begin{array}{ll}0& 0\\ q_{23}& q_{24}\end{array}\\ \begin{array}{ll}0& 0\\ q_{41}& 0\end{array}\hspace{0.17em}\begin{array}{ll}0& 1\\ q_{43}& q_{44}\end{array}\end{array}\right]$

$q_{23}=-\left(k_1+{3k}_3{y}_1^2+{h\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)\right)\bullet \frac{1}{{\mathrm{\Omega }}^2{m}_1}$

$q_{24}=-\left(c+h{\mathrm{\Omega }}^{p-1}\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)\right)\bullet \frac{1}{\mathrm{\Omega }{m}_1}$

$q_{41}=\frac{k_t}{{\mathrm{\Omega }}^2{m}_2}$

$q_{43}=-\left[k_t+k_1+{3k}_3{y}_1^2+{h\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)\right]\bullet \frac{1}{{\mathrm{\Omega }}^2{m}_2}-\left(k_1+{3k}_3{y}_1^2+{h\mathrm{\Omega }}^p\hspace{0.17em}\cos \left(\frac{p\pi }{2}\right)\right)\bullet \frac{1}{{\mathrm{\Omega }}^2{m}_1}$

$q_{44}=-\left[c+h{\mathrm{\Omega }}^{p-1}\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)\right]\bullet \frac{1}{\mathrm{\Omega }{m}_2}-\left(c+h{\mathrm{\Omega }}^{p-1}\hspace{0.17em}\sin \left(\frac{p\pi }{2}\right)\right)\bullet \frac{1}{\mathrm{\Omega }{m}_1}$

Since ${\mathit{X}}_0$ is a periodic function of $\tau$ , the period is T = $2\pi \mathrm{\Omega }$ , each element of Q is also a periodic function with the same period as ${\mathit{X}}_0$ . For equation (34), the following fundamental solution matrix exists:

$\mathit{Y}=\left[\begin{array}{l}\begin{array}{lll}{y}_{11}& y_{12}& y_{13}\\ y_{21}& y_{22}& y_{23}\\ y_{31}& y_{32}& y_{33}\end{array}\hspace{0.17em}\begin{array}{l}{y}_{14}\\ y_{24}\\ y_{34}\end{array}\\ \begin{array}{ll}{y}_{41}& y_{42}\end{array}\hspace{0.17em}\begin{array}{ll}{y}_{43}& y_{44}\end{array}\end{array}\right]$

(35)

According to the Floquet theory, there is a transition matrix so that the periodic solution stability can be decided by the eigenvalues of the matrix ${\mathit{T}}_M$ .

$\mathit{Y}\left(\tau +T\right)={\mathit{T}}_M\mathit{Y}\left(\tau \right)$

(36)

Each period T is equally divided into N parts, where ${∆}_k={\tau }_k-{\tau }_{k-1}$ , in each part time interval, the period coefficient matrix $\mathit{Q}\left(\tau \right)$ is replaced by the constant coefficient matrix ${\mathit{Q}}_k$ .

${\mathit{Q}}_k=\frac{1}{{∆}_k}{\int }_{{\tau }_{k-1}}^{{\tau }_k}\mathit{Q}\left(\tau \right)d\tau$

(37)

Therefore, the transformation matrix ${\mathit{T}}_M$ can be expressed as

${\mathit{T}}_M={\displaystyle \prod }_{i=1}^N\mathit{exp}\left({\mathit{Q}}_i{∆}_i\right)$

(38)

The terms $\mathit{exp}\left(Q_k{∆}_k\right)$ are expanded with the Taylor formula for programming:

${\mathit{T}}_M={\displaystyle \prod }_{i=1}^N\left(I+{\displaystyle \sum }_{j=1}^M\frac{{\left[{∆}_i{\mathit{Q}}_i\right]}^j}{j!}\right)$

(39)

where M is the number of terms in the approximation of the constant matrix exponential. ³³ It is proved that M has sufficient precision when it is greater than or equal to 4. In this paper, M is set to 4. The stability of the periodic solution can be decided by the modulus of the eigenvalues of the matrix ${\mathit{T}}_M$ . When none of the moduli is greater than one, the periodic solution is considered stable, otherwise it is considered unstable.

Simulation results and discussions

Harmonics number selection and result verification

The number of harmonics directly affects the accuracy of steady-state response. A high number N corresponds not only to a high calculation accuracy but also to a long calculation time. A compromise needs to be reached between the number of harmonics. In this paper, N = 4, 5, and 6 are selected for calculation, respectively. The simulation parameters are shown in Table 1. The comparisons of the calculated results are shown in Figure 4. There is a pronounced deviation of the amplitude–frequency response curve when N = 4 and N = 5. However, no significant difference is observed when N = 5 and N = 6. The magnitude of the deviation after enlarging the image is only 10⁻⁶ times, indicating that the accuracy of N = 5 is very high. Therefore, this paper chooses the result of N = 5 for analysis and research, which can guarantee the accuracy requirement with faster calculation speed.

Figure 4.

Amplitude–frequency response curves with different harmonic numbers. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

The power series expansion method is used for numerical simulation, in which the step size is $\frac{\pi }{100\omega }$ , 300 cycles are calculated, and the maximum amplitude of the following 150 cycles is taken as the stable amplitude.

The IHBM results are compared with the numerical simulation power series expansion method (PSEM) results, as shown in Figure 5. The result calculated through the IHMB is almost identical with the numerical method at the amplitude of spring–mass vibration.

Figure 5.

Comparisons of amplitude–frequency response curves of PSEM and IHBM. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Although a slight difference exists between IHBM and the numerical method in the mid and high frequency of suspension dynamic deflection, it remains within a narrow margin. Therefore, the 2-DOF IHBM algorithm with fractional order is feasible and highly accurate.

Influence of fractional-order parameters on the dynamic performance of the vehicle suspension system

When the system contains a fractional-order differential term, two specific parameters exist: the order of the differential term and the coefficient of the differential term. Their effects on the dynamic performance of the system are studied.

Effect of the fractional order on the dynamic behavior of the vehicle system

The system takes the fractional order of 0, 0.2, 0.4, 0.6, 0.8, and 1.0 to study. The body vibration amplitude–frequency response curve after changing the system parameters is shown in Figure 6. As can be seen from the figure, a low-fractional order corresponds to a high resonance peak value, which is consistent with the physical meaning that a low-fractional order corresponds to a high proportion of the stiffness part and a low proportion of the damping part.

Figure 6.

Comparisons of amplitude–frequency responses of different fractional orders. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

An obvious contrast exists between the amplitude–frequency response curves of the fractional order and those of the integer order, thereby indicating that the fractional-order differential term is more accurate for describing the system and closer to the actual system.

Effect of fractional coefficient on the dynamic behavior of the vehicle system

The fractional differential coefficients of 500, 1000, 1500, and 2000 were used in the system for research. The amplitude–frequency response curve of body vibration after changing the parameters is shown in Figure 7. The diagram suggests that a low fractional-order coefficient corresponds to a high amplitude of the system vibration, a high resonance peak, and a more uniform increase.

Figure 7.

Comparisons of amplitude–frequency responses for different fractional-order coefficients. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Analysis of the influence of parameters on nonlinear characteristics

Effect of nonlinear stiffness on suspension nonlinear dynamic performance

The system exhibits nonlinear characteristics because of the cubic stiffness term in the suspension. IHBM is suitable not only for weakly nonlinear systems but also for calculations of strongly nonlinear systems. With a change in the value of cubic stiffness factor k₃, when k₃ = 1000, 2000, 3000, 4000, and 5000 kN/m, the amplitude–frequency response curve is calculated as shown in Figure 8. This figure also shows that when the excitation frequency is low, the deflection of resonance peaks to the right becomes more significant with the increase in k₃. The nonlinear characteristics are not distinct at 1000 kN/m. IHBM can also obtain continuous amplitudes. Starting from 2000 kN/m, the amplitude–frequency response curve shows typical nonlinear characteristics at the low-frequency resonance peak. That exhibits multi-solution and discontinuous jumping phenomena and bends to the high frequency. The vehicle suspension system is strongly nonlinear under low-frequency excitation, and the resonance peak area should be avoided in the design to reduce the nonlinear effect.

Figure 8.

Comparisons of amplitude–frequency responses for different nonlinear stiffness coefficients. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Effect of external excitation amplitude on the dynamic behavior of the vehicle system

With k₃ = 3000 kN/m and other parameters unchanged, the calculated amplitude–frequency characteristic curve is shown in Figure 9 when the excitation amplitudes are changed to F = 0.006 m, F = 0.008 m, and F = 0.01 m. The calculation not only shows that the excitation amplitude affects the resonance peak value of the suspension system but also suggests that a high excitation amplitude corresponds to more obvious nonlinear multi-value characteristics of the system. As shown in the figure, when F = 0.01 m, the coordinate points of the discontinuous region in the forward and reverse scanning can be known to have a double solutions region and are no longer notable when F is reduced in the double solution region. However, the excitation amplitude in the high-frequency section will affect the resonance peak of the system. A high excitation amplitude corresponds to a high resonance peak value, but the nonlinear characteristics of the system are not notable.

Figure 9.

Comparisons of amplitude–frequency responses of different excitation amplitudes. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Effect of mass ratio on the dynamic behavior of the vehicle system

The calculated amplitude–frequency response curve under the constant condition of k₃ = 3000 kN/m when the body mass is changed to m₁ = 200, 220, 240, 260, or 280 kg is shown in Figure 10. The figure shows that the body mass change has almost no effect on the high-frequency resonance peak but a significant impact on the low-frequency resonance peak. As the mass increases, the resonance peak shifts to the left, and the nonlinear multi-value area tends to be more significant. For example, when m₁ = 280 and 240 kg, the discontinuous position coordinates of forward sweeping and reverse sweeping appear in the figure. A large body mass corresponds to a multi-value region.

Figure 10.

Comparisons of amplitude–frequency response of different body mass. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Effect of tire quality on suspension nonlinear performance

When the tire quality is modified, and m₂ = 30, 40, 50, 60, or 70 kg is taken, the amplitude–frequency response curve is calculated as shown in Figure 11. As the figure presents, tire quality does not affect low-frequency resonance peaks, and the amplitude–frequency characteristics of the low-frequency band are identical. Different quality tires in the low-frequency band influence the high-frequency resonance peak profoundly. If the tire quality is changed, then the resonance peak frequencies vary. However, the nonlinear characteristics of the system are not affected.

Figure 11.

Comparisons of amplitude frequency response of different tire masses. (a) Amplitude–frequency response curve of body displacement. (b) Amplitude–frequency response of suspension dynamic deflection.

Stability analysis

When k3 is set to 5,000,000, there is a strong nonlinearity in the low-frequency band. This part is used as an example to discuss the periodic solution stability.

As shown in Figure 12, the curve of vehicle body vibration change with excitation frequency $\mathrm{\omega }$ in low-frequency band is shown, in which the green part is the stable solution and the black dotted line is the unstable part. Similarly, other curves in the paper also have similar characteristics.

Figure 12.

Stability of the periodic solution.

Discussion

The amplitude–frequency curve is a description of the transmission characteristics of the input road excitation in the frequency domain, and it can reflect both the nonlinear characteristics and steady-state of the vehicle suspension system.

(1) The amplitude–frequency response characteristics are significantly affected by the order and the coefficient of the fractional-order system. A low-order corresponds to a high amplitude of system vibration and a high resonance peak value. A high coefficient corresponds to a low resonance peak of the frequency response, which differs from the integer-order comparison. The results show that fractional-order parameters profoundly affect the resonance peak value of the system but have no effect on the resonance frequency.

(2) With the increase in the nonlinear term $k_3$ , the nonlinear characteristics of the system become more pronounced. The system exhibits nonlinear characteristics, and unstable solution regions appear. Multi-value solutions are obtained when the frequency is swept in both directions. The resonance peaks tend to deviate to the right. The nonlinear characteristic of the low-frequency band is more significant than the high-frequency band.

Varying the excitation amplitude as the control parameter, a high excitation amplitude corresponds to the great resonance peak value of the system, and a high excitation amplitude corresponds to more significant nonlinear characteristics in the low-frequency band.

The body mass $m_1$ significantly influences the low-frequency resonance peak and frequency but does not influence the high-frequency resonance peak. Meantime, the tire mass $m_2$ influences the high-frequency resonance peak and frequency profoundly but does not influence the low-frequency resonance peak.

Depending on different system parameters, the nonlinear stiffness, the excitation amplitude, and mass ratio have different degrees of effect on the nonlinear dynamic response. The interaction makes the dynamic behaviors of the vehicle suspension system more complicated.

Compared with the previous research work on 2-DOF vehicle systems, the IHBM extended to fractional multi-DOF can not only reflect the nonlinear characteristics of the vehicle system but also analyze the impact of fractional parameters on the performance of the system.

Conclusion

This paper innovatively proposed the nonlinear FOM of a 2-DOF vehicle system, and the general calculation format of the 2-DOF incremental harmonic balance method containing the fractional-order nonlinear system is derived. Then, the amplitude–frequency response characteristics of a 2-DOF nonlinear vehicle suspension system with fractional order under external harmonic excitation are studied by using this scheme. The effects of fractional parameters and other variables on the nonlinear characteristics of the system are also studied. The conclusions are as follows: (1) Because the vehicle suspension system model is FOM, a margin must be set when designing a suspension using IOM. (2) The key parameters affect the dynamic behaviors of the FOM suspension system.

Incremental harmonic balance method is a semi-analytical and semi-numerical method. The results have intuitive expressions, which not only clarify the physical meaning of the parameters but also make it easy to analyze the influence of parameter changes on the vibration characteristics of the system. For complex, strongly nonlinear vibration problems, the successful promotion and application of the fractional-order incremental harmonic balance method have high theoretical value and broad application prospects in the fractional-order vibration engineering design.