Dynamic analysis of the response of Duffing-type oscillators subject to interacting parametric and external excitations

Abstract

Prediction of the response of nonlinear dynamical systems under interacting parametric and external excitations is important in designing systems such as sensors, amplifiers or energy harvesters, to achieve the desired performance. This paper concerns the nonlinear forced Mathieu equation, with linear damping and a 2:1 ratio between the parametric and external excitation frequencies. The Method of Varying Amplitudes (MVA) is employed to derive approximate analytical expressions for the response of the system. Both single-term and double-term solutions are developed: it is seen that, employing the double-term approximation, the MVA can accurately predict the response of the system over a wide range of frequencies and system parameters, showing a maximum of 0.2% deviation from numerical results obtained by direct integration of the equation of motion. This is in contrast with most of the available theoretical approaches such as the conventional Method of Multiple Scales, which can predict the response accurately only for a narrow range of system parameters and excitation frequencies. Furthermore, it is seen that the response is bounded, and analytical expressions for the frequency and amplitude of the upper bound are developed: this is unlike other methods which predict unbounded response, unless nonlinear damping is considered. Analytical expressions for the response are developed, and results are verified with numerical results obtained from direct integration of the equation of motion. Numerical examples are presented, showing good agreement with results obtained by the MVA.

Appendices

Appendix A: second approximation of MMS assuming Duffing nonlinearity term $\eta$ to be $O(\varepsilon )$

In Appendix, results from the second-order MMS are briefly reviewed [33]. To obtain the frequency response equation of a dynamical system under pure parametric excitation with no external excitation ($d = 0$) using the MMS, assuming the system parameters including the Duffing nonlinearity term are all small and of order $\varepsilon$, we express Eq. (1) as

$$ \ddot{u} + \varepsilon \beta_{\varepsilon } \dot{u} + \omega_{0}^{2} \left( {1 + \varepsilon P_{\varepsilon } \cos \left( {\varOmega t} \right)} \right)u + \varepsilon \eta_{\varepsilon } u^{3} = 0, $$

(A1)

where

$$ \varepsilon \beta_{\varepsilon } = \beta ,\;\varepsilon P_{\varepsilon } = P,\;\varepsilon \eta_{\varepsilon } = \eta . $$

(A2)

To the second approximation of the MMS, the response of the system can be written as [33]

$$ u\left( {T_{0} ,T_{1} ,T_{2} } \right) = u_{0} \left( {T_{0} ,T_{1} ,T_{2} } \right) + \varepsilon u_{1} \left( {T_{0} ,T_{1} ,T_{2} } \right) + \varepsilon^{2} u_{2} \left( {T_{0} ,T_{1} ,T_{2} } \right) + O\left( {\varepsilon^{3} } \right), $$

(A3)

where $T_{0} = t$, $T_{1} = \varepsilon t$ and $T_{2} = \varepsilon^{2} t$ are time scales. Defining the operator

$$ D_{n} = \frac{\partial }{{\partial T_{n} }},\;\;n = 0,1,2,..., $$

(A4)

the time derivatives can be expressed as

$$ \frac{{\text{d}}}{{{\text{d}}t}} = D_{0} + \varepsilon D_{1} + \varepsilon^{2} D_{2} + O(\varepsilon^{3} ), $$

(A5)

$$ \frac{{{\text{d}}^{2} }}{{{\text{d}}t^{2} }} = D_{0}^{2} + 2\varepsilon D_{0} D_{1} + \varepsilon^{2} \left( {2D_{0} D_{1} + D_{1}^{2} } \right) + O\left( {\varepsilon^{3} } \right), $$

(A6)

Substituting Eq. (A3) into Eq. (A1), considering Eqs. (A4)–(A6) and equating coefficients of like powers of $\varepsilon$, we obtain

$$ O(1):\;\;D_{0}^{2} u_{0} + \omega_{0}^{2} u_{0} = 0, $$

(A7)

$$ O(\varepsilon ):D_{0}^{2} u_{1} + \omega_{0}^{2} u_{1} = - 2D_{0} D_{1} u_{0} - \beta_{\varepsilon } D_{0} u_{0} - \eta_{\varepsilon } u_{0}^{3} - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \left( {e^{{i\varOmega T_{0} }} + e^{{ - i\varOmega T_{0} }} } \right)u_{0} , $$

(A8)

$$ \begin{aligned} O(\varepsilon^{2} ):D_{0}^{2} u_{2} + \omega_{0}^{2} u_{2} = & - D_{1}^{2} u_{0} - 2D_{0} D_{2} u_{0} - 2D_{0} D_{1} u_{1} - 3\eta_{\varepsilon } u_{0}^{2} u_{1} \\ & \;\; - \beta_{\varepsilon } \left( {D_{1} u_{0} + D_{0} u_{1} } \right) - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \left( {e^{{i\varOmega T_{0} }} + e^{{ - i\varOmega T_{0} }} } \right)u_{1} . \\ \end{aligned} $$

(A9)

The solution of Eq. (A7) is assumed to be of the form

$$ u_{0} = A\left( {T_{1} ,T_{2} } \right)e^{{i\omega_{0} T_{0} }} + \overline{A}\left( {T_{1} ,T_{2} } \right)e^{{ - i\omega_{0} T_{0} }} , $$

(A10)

where $\overline{A}\left( {T_{1} ,T_{2} } \right)$ is the complex conjugate of $A\left( {T_{1} ,T_{2} } \right)$. Substituting Eq. (A10) into Eq. (A8), we obtain

$$ \begin{aligned} D_{0}^{2} u_{1} + \omega_{0}^{2} u_{1} = & \left( { - 2i\omega_{0} D_{1} A - i\omega_{0} \beta_{\varepsilon } A - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \overline{A}e^{{i\left( {\Omega - 2\omega_{0} } \right)T_{0} }} - 3\eta_{\varepsilon } A^{2} \overline{A}} \right)e^{{i\omega_{0} T_{0} }} \\ & \;\; - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } Ae^{{i\left( {\varOmega + \omega_{0} } \right)T_{0} }} - \eta_{\varepsilon } A^{3} e^{{3i\omega_{0} T_{0} }} + {\text{CC}}, \\ \end{aligned} $$

(A11)

where ${\text{CC}}$ is the complex conjugate of the preceding terms. Considering $\varOmega \approx 2\omega_{0}$, the coefficients of $\exp \left( {i\omega_{0} T_{0} } \right)$ in Eq. (A11) represent the secular terms which can be eliminated when

$$ D_{1} A = - \frac{1}{2}\beta_{\varepsilon } A + \frac{i}{4}\omega_{0} P_{\varepsilon } \overline{A}e^{{i\left( {\varOmega - 2\omega_{0} } \right)T_{0} }} + \frac{3i}{{2\omega_{0} }}\eta_{\varepsilon } A^{2} \overline{A}. $$

(A12)

Consequently, the particular solution of $u_{1}$ in Eq. (A11) is obtained as

$$ u_{1} = \frac{{\omega_{0}^{2} P_{\varepsilon } A}}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}e^{{i\left( {\varOmega + \omega_{0} } \right)T_{0} }} + \frac{1}{{8\omega_{0}^{2} }}\eta_{\varepsilon } A^{3} e^{{3i\omega_{0} T_{0} }} + CC. $$

(A13)

Substituting $u_{0}$ and $u_{1}$ from Eqs. (A10) and (A13) into Eq. (A9), we obtain

$$ \begin{aligned} D_{0}^{2} u_{2} + \omega_{0}^{2} u_{2} = & \left( { - D_{1}^{2} A - \beta_{\varepsilon } D_{1} A - 2i\omega_{0} D_{2} A - \frac{{\omega_{0}^{4} P_{\varepsilon }^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{3}{{8\omega_{0}^{2} }}\eta_{\varepsilon }^{2} \overline{A}^{2} A^{3} - \frac{{3\omega_{0}^{2} \eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}e^{{i\left( {\varOmega - 2\omega_{0} } \right)T_{0} }} - \frac{1}{16}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\left( {\varOmega - 2\omega_{0} } \right)T_{0} }} } \right)e^{{i\omega_{0} T_{0} }} \\ & \;\; + {\text{NST}} + {\text{CC}} \\ \end{aligned} $$

(A14)

where ${\text{NST}}$ denotes the non-secular terms. Defining a frequency detuning parameter $\sigma$ such that

$$ \varOmega - 2\omega_{0} = \varepsilon \sigma , $$

(A15)

eliminating the secular terms in Eq. (A14) give

$$ D_{1}^{2} A = - \beta_{\varepsilon } D_{1} A - 2i\omega_{0} D_{2} A - \frac{{\omega_{0}^{4} P_{\varepsilon }^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{3}{{8\omega_{0}^{2} }}\eta_{\varepsilon }^{2} \overline{A}^{2} A^{3} - \frac{{3\omega_{0}^{2} \eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}e^{{i\sigma T_{1} }} - \frac{1}{16}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\sigma T_{1} }} . $$

(A16)

Furthermore, Eq. (A12) can be rewritten as

$$ \begin{aligned} D_{1}^{2} A = & \frac{1}{16}A\left( {\omega_{0}^{2} P_{\varepsilon }^{2} + 4\beta_{\varepsilon }^{2} } \right) - \frac{3i}{{\omega_{0} }}\beta_{\varepsilon } \eta_{\varepsilon } A^{2} \overline{A} - \frac{9}{{4\omega_{0}^{2} }}\eta_{\varepsilon }^{2} A^{3} \overline{A}^{2} \\ & \;\; - \left( {\frac{i}{4}\omega_{0} \beta_{\varepsilon } P_{\varepsilon } \overline{A} + \frac{1}{4}\omega_{0} P_{\varepsilon } \sigma \overline{A} + \frac{3}{8}\eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} } \right)e^{{i\sigma T_{1} }} \\ & \;\; + \frac{3}{8}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\sigma T_{1} }} . \\ \end{aligned} $$

(A17)

Also, from Eq. (A16) and considering Eq. (12) we obtain

$$ \begin{aligned} D_{1}^{2} A = & \frac{1}{2}\beta_{\varepsilon }^{2} A - \frac{3i}{{2\omega_{0} }}\beta_{\varepsilon } \eta_{\varepsilon } A^{2} \overline{A} - 2i\omega_{0} D_{2} A - \frac{{\omega_{0}^{4} P_{\varepsilon }^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} \\ & \;\; - \frac{3}{{8\omega_{0}^{2} }}\eta_{\varepsilon }^{2} A^{3} \overline{A}^{2} - \left( {\frac{i}{4}\omega_{0} \beta_{\varepsilon } P_{\varepsilon } \overline{A} + \frac{{3\omega_{0}^{2} \eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}} \right)e^{{i\sigma T_{1} }} \\ & \;\; - \frac{1}{16}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\sigma T_{1} }} . \\ \end{aligned} $$

(A18)

Combining Eqs. (A17) and (A18) and cancelling $D_{1}^{2} A$ give

$$ \begin{aligned} 2i\omega_{0} D_{2} A = & \frac{1}{4}\beta_{\varepsilon }^{2} A - \frac{1}{16}\omega_{0}^{2} P_{\varepsilon }^{2} A - \frac{{\omega_{0}^{4} P_{\varepsilon }^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} + \frac{3i}{{2\omega_{0} }}\beta_{\varepsilon } \eta_{\varepsilon } A^{2} \overline{A} \\ & \;\; + \frac{15}{{8\omega_{0}^{2} }}\eta_{\varepsilon }^{2} A^{3} \overline{A}^{2} + \left( {\frac{1}{4}\omega_{0} P_{\varepsilon } \sigma \overline{A} + \frac{3}{8}\eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} - \frac{{3\omega_{0}^{2} \eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}} \right)e^{{i\sigma T_{1} }} \\ & \;\; - \frac{7}{16}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\sigma T_{1} }} . \\ \end{aligned} $$

(A19)

Multiplying Eq. (A12) by $2i\omega_{0} \varepsilon$ and Eq. (A19) by $\varepsilon^{2}$ and combining the resulting equations yield

$$ \begin{aligned} 2i\omega_{0} \left( {\varepsilon D_{1} A + \varepsilon^{2} D_{2} A} \right) + & \varepsilon \left( {i\omega_{0} \beta_{\varepsilon } A + \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \overline{A}e^{{i\sigma T_{1} }} + 3\eta_{\varepsilon } A^{2} \overline{A}} \right) \\ + & \varepsilon^{2} \left( {\frac{1}{16}\omega_{0}^{2} P_{\varepsilon }^{2} A + \frac{{\omega_{0}^{4} P_{\varepsilon }^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{1}{4}\beta_{\varepsilon }^{2} A - \frac{3i}{{2\omega_{0} }}\beta_{\varepsilon } \eta_{\varepsilon } A^{2} \overline{A} - \frac{15}{{8\omega_{0}^{2} }}\eta_{\varepsilon }^{2} A^{3} \overline{A}^{2} - \frac{1}{4}\omega_{0} P_{\varepsilon } \sigma \overline{A}e^{{i\sigma T_{1} }} - \frac{3}{8}\eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} e^{{i\sigma T_{1} }} + \frac{{3\omega_{0}^{2} \eta_{\varepsilon } P_{\varepsilon } A\overline{A}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}e^{{i\sigma T_{1} }} + \frac{7}{16}\eta_{\varepsilon } P_{\varepsilon } A^{3} e^{{ - i\sigma T_{1} }} } \right) = 0. \\ \end{aligned} $$

(A20)

Taking Eq. (A2) into account, Eq. (A20) can be rewritten as

$$ \begin{aligned} 2i\omega_{0} \frac{{{\text{d}}A}}{{{\text{d}}t}} - & \frac{1}{4}\omega_{0} P\overline{A}\left( {\varOmega - 4\omega_{0} } \right)e^{{i\left( {\varOmega - 2\omega_{0} } \right)t}} - \frac{1}{4}\beta^{2} A + i\omega_{0} \beta A + 3\eta A^{2} \overline{A} \\ + & \frac{1}{16}\omega_{0}^{2} p^{2} A + \frac{{\omega_{0}^{4} p^{2} A}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{3i}{{2\omega_{0} }}\beta \eta A^{2} \overline{A} - \frac{15}{{8\omega_{0}^{2} }}\eta^{2} A^{3} \overline{A}^{2} \\ + & \left( {\frac{{3\omega_{0}^{2} }}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{3}{8}} \right)\eta PA\overline{A}^{2} e^{{i\left( {\varOmega - 2\omega_{0} } \right)t}} + \frac{7}{16}\eta PA^{3} e^{{ - i\left( {\varOmega - 2\omega_{0} } \right)t}} = 0, \\ \end{aligned} $$

(A21)

where $A$ is a function of time and can be expressed in a polar form

$$ A\left( t \right) = \frac{1}{2}a\left( t \right)e^{i\lambda \left( t \right)} , $$

(A22)

where $a\left( t \right)$ and $\lambda \left( t \right)$ are real. substituting Eq. (A22) into Eq. (A21), we obtain

$$ \begin{aligned} - 256\omega_{0} a\dot{\lambda } + & \left( {\omega_{0}^{2} p^{2} + \frac{{32\omega_{0}^{4} p^{2} }}{{\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - 32\beta^{2} } \right)a + 96\eta a^{3} \\ - & \frac{15}{{\omega_{0}^{2} }}\eta^{2} a^{5} + i\left( {256\omega_{0} \dot{a} + 128\omega_{0} \beta a - \frac{48}{{\omega_{0} }}\beta \eta a^{3} } \right) \\ & + \left( {\left( {\frac{{48\omega_{0}^{2} }}{{\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - 12} \right)\eta Pa^{3} - 32\omega_{0} P\left( {\varOmega - 4\omega_{0} } \right)a} \right)e^{{i\left( {\left( {\varOmega - 2\omega_{0} } \right)t - 2\lambda } \right)}} \\ + & 14\eta Pa^{3} e^{{ - i\left( {\left( {\varOmega - 2\omega_{0} } \right)t - 2\lambda } \right)}} = 0. \\ \end{aligned} $$

(A23)

Applying the transformation

$$ \tau \left( t \right) = \left( {\varOmega - 2\omega_{0} } \right)t - 2\lambda \left( t \right) $$

(A24)

into Eq. (A23) and separating the resultant real and imaginary parts, the system of equations

$$ \begin{aligned} a\dot{\tau } = & \left( {\varOmega - 2\omega_{0} - \frac{1}{16}\omega_{0} p^{2} - \frac{{\omega_{0}^{3} p^{2} }}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} + \frac{1}{{4\omega_{0} }}\beta^{2} } \right)a + \frac{15}{{128\omega_{0}^{3} }}\eta^{2} a^{2} \\ & \;\; - \frac{3}{{4\omega_{0} }}\eta a^{3} - \left( {\left( {\frac{{3\omega_{0} }}{{8\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} + \frac{1}{{64\omega_{0} }}} \right)\eta Pa^{3} - \frac{1}{4}P\left( {\varOmega - 4\omega_{0} } \right)a} \right)\cos \left( \tau \right) = 0, \\ \end{aligned} $$

(A25)

$$ \begin{aligned} \dot{a} = & - \frac{1}{2}\beta a + \frac{3}{{16\omega_{0}^{2} }}\beta \eta a^{3} \\ & \;\; - \left( {\left( {\frac{{3\omega_{0} }}{{16\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - \frac{13}{{128\omega_{0} }}} \right)\eta Pa^{3} - \frac{1}{8}P\left( {\varOmega - 4\omega_{0} } \right)a} \right)\sin \left( \tau \right), \\ \end{aligned} $$

(A26)

is obtained. Consequently, solving Eqs. (A25) and (A26) for the steady state ($\dot{a} = \dot{\tau } = 0$), the frequency response Eq. (44) is obtained (with symbol $a$ replaced by $A$ for consistency with the results obtained by the MVA).

Appendix B: second approximation of MMS assuming Duffing nonlinearity $\eta$ to be $O(\varepsilon^{2} )$

To the second approximation of the MMS, assuming $\eta$ is of order $\varepsilon^{2}$, the equation of motion is scaled as

$$ \ddot{u} + \varepsilon \beta_{\varepsilon } \dot{u} + \omega_{0}^{2} \left( {1 + \varepsilon P_{\varepsilon } \cos \left( {\varOmega t} \right)} \right)u + \varepsilon^{2} \eta_{\varepsilon } u^{3} = 0, $$

(B1)

where $\varepsilon^{2} \eta_{\varepsilon } = \eta$. Considering the system parameters to be of different orders of ε has been used commonly in the literature for various problems [25, 32, 33]. Here, $\eta_{\varepsilon }$ is considered to be of order $\varepsilon^{2}$, which only implies the degree of smallness for Duffing nonlinearity. Considering Eqs. (A3)–(A6) and following a similar approach, Eq. (A7) still holds true for coefficients of order $\varepsilon^{0}$, while the equations representing coefficients of order $\varepsilon$ and $\varepsilon^{2}$ change to

$$ D_{0}^{2} u_{1} + \omega_{0}^{2} u_{1} = - 2D_{0} D_{1} u_{0} - \beta_{\varepsilon } D_{0} u_{0} - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \left( {e^{{i\varOmega T_{0} }} + e^{{ - i\varOmega T_{0} }} } \right)u_{0} , $$

(B2)

$$ \begin{aligned} D_{0}^{2} u_{2} + \omega_{0}^{2} u_{2} = & - D_{1}^{2} u_{0} - 2D_{0} D_{2} u_{0} - 2D_{0} D_{1} u_{1} \\ & - \beta_{\varepsilon } \left( {D_{1} u_{0} + D_{0} u_{1} } \right) - \frac{1}{2}\omega_{0}^{2} P_{\varepsilon } \left( {e^{{i\varOmega T_{0} }} + e^{{ - i\varOmega T_{0} }} } \right)u_{1} - \eta_{\varepsilon } u_{0}^{3} , \\ \end{aligned} $$

(B3)

respectively. The particular solution for $u_{1}$ in Eq. (B2) is

$$ u_{1} = \frac{{\omega_{0}^{2} P_{\varepsilon } A}}{{2\varOmega \left( {\varOmega + 2\omega_{0} } \right)}}e^{{i\left( {\varOmega + \omega_{0} } \right)T_{0} }} + CC. $$

(B4)

Following an approach similar to that of Eqs. (A14)–(A22) in Appendix A, the modulation equation is obtained as

$$ \begin{aligned} 12\eta a^{3} - & 4\omega_{0} Pa\left( {\varOmega - 4\omega_{0} } \right)e^{{i\left( {\left( {\varOmega - 2\omega_{0} } \right)t - 2\lambda } \right)}} + 32i\omega_{0} \dot{a} + 16i\omega_{0} \beta a \\ & + \omega_{0}^{2} p^{2} a + \frac{{4\omega_{0}^{4} p^{2} a}}{{\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} - 4\beta^{2} a - 32a\omega_{0} \dot{\lambda } = 0. \\ \end{aligned} $$

(B5)

Consequently, applying the transformation (A24), the system of equations

$$ \begin{aligned} a\dot{\tau } = & a\left( {\varOmega - 2\omega_{0} } \right) - \frac{{3\eta a^{3} }}{{4\omega_{0} }} - \frac{1}{16}\omega_{0} p^{2} a - \frac{{\omega_{0}^{3} p^{2} a}}{{4\varOmega \left( {\varOmega + 2\omega_{0} } \right)}} \\ & + \frac{1}{4}\beta^{2} a + \frac{1}{2}Pa\left( {\varOmega - 4\omega_{0} } \right)\cos \left( \tau \right), \\ \end{aligned} $$

(B6)

$$ \dot{a} = - \frac{1}{2}\beta a + \frac{1}{8}Pa\left( {\varOmega - 4\omega_{0} } \right)\sin \left( \tau \right), $$

(B7)

is obtained. Solving Eqs. (B6) and (B7) for the steady state ($\dot{a} = \dot{\tau } = 0$), the frequency response Eq. (46) is obtained (with symbol $a$ replaced by $A$ for consistency with the results obtained by the MVA).