SUPER-HARMONIC AND SUB-HARMONIC SIMULTANEOUS RESONANCES OF FRACTIONAL-ORDER DUFFING OSCILLATOR
-
摘要: 研究了含分数阶微分项的达芬（Duffing）振子的超谐与亚谐联合共振.采用平均法得到了系统的一阶近似解析解，提出了超、亚谐联合共振时等效线性阻尼和等效线性刚度的概念.建立了联合共振定常解幅频曲线的解析表达式，并对联合共振幅频响应的近似解析解和数值解进行了比较，二者吻合良好，证明了求解过程及近似解析解的正确性.然后，将等效线性阻尼和等效线性刚度的概念与传统整数阶系统进行比较，证明分数阶微分项不仅起着阻尼的作用同时还起着刚度的作用.最后，通过数值仿真研究了不同的分数阶微分项系数和阶次对联合共振幅频曲线多值性和跳跃现象的影响，并与单一频率下超谐共振或亚谐共振进行了对比.研究发现，分数阶微分项系数与阶次不仅影响着系统的响应幅值、共振频率，同时还对系统的周期解个数、发生区域面积、发生先后等有重要影响.并且，在不同的基本参数下该系统分别表现出单独超谐共振、单独亚谐共振以及超谐共振和亚谐共振同时存在的现象.这些结果对系统动力学特性的研究具有重要意义.Abstract: The super-harmonic and sub-harmonic simultaneous resonance of Duffing oscillator with fractional-order derivative is studied in this paper. The first-order approximate analytical solution is obtained by averaging method. The definitions of equivalent linear damping coefficient and equivalent linear stiffness efficient for super-harmonic and subharmonic simultaneous resonance are presented. The analytical amplitude-frequency equation for steady-state solution of simultaneous resonance is established. A comparison of the analytical solution with the numerical results is made, and their satisfactory agreement verifies the correctness and higher-order precision of the approximately analytical results. Then, a further comparison between the fractional-order and traditional integer-order Duffing oscillator is fulfilled through the definitions of equivalent linear damping coefficient and equivalent linear stiffness coefficient, and the results prove that the fractional-order parameters has the effects of both damping and stiffness, which is similar in other fractionalorder system. At last the numerical simulation is used to analyze the effects of different fractional-order parameters on multi-value characteristics and jumping phenomena of the amplitude-frequency curve under simultaneous resonance, and the differences between the super-harmonic and sub-harmonic resonances under single-frequency excitation are analyzed in detail. It could be found that the fractional-order parameters not only affect the response amplitude and resonant frequency of the system, but also has significant influence on the number, existing area and the occurrence order of periodic solutions. Moreover, single super-harmonic resonance, single sub-harmonic resonance and both existing of these two resonances could be respectively found under different basic parameters, which is important to study the dynamic characteristics of the similar system.
-
引言
20世纪70年代，研究人员发现分形几何、幂律现象和记忆过程等相关现象能够与分数阶微积分建立起密切的联系，而且分数阶微积分还可以是一种很好的描述与刻画手段，由此将分数阶微积分的研究推向高潮[1-4].迄今为止，国内外的学者对于非线性动力学的研究仍然没有停息[5-8]，有关分数阶微积分的应用研究得到迅速发展[9-12]，越来越多的学者投入到了这个新兴的领域中.
在工程问题中，分数阶微积分的应用领域主要分为两类, 一类是在控制系统中引入分数阶微积分项来改善控制效果.如无人机滑模姿态的控制[13]、永磁同步电动机的双闭环控制[14]、速度反馈的PID控制[15]、集成神经网络的控制[16]等领域.另一类是用来模拟含记忆特性的工程材料的真实本构关系.如隔振器件的建模(油气悬架的建模[17]、空气悬架的建模[18]、磁流变阻尼器的建模[19])、单相逆变器的建模[20]、非定常蠕变本构模型的研究[21]、黏弹性材料变形研究[22]等领域.另外达芬方程作为一类典型的动力学系统，能够描述物理与工程领域中丰富的非线性动力学行为而称为经典模型.在实际工程中，例如船的横摇运动、结构振动、转子的故障检测等诸多非线性振动的数学模型都可以转化为达芬方程.近几年，由于分数阶微分模型能够在较大频率范围内描述材料的力学性能而成为倍受重视的本构模型，因此本文将分数阶微分项引入到达芬振子中可以更加真实地反应系统特性.例如：陈林聪等[23]研究含分数阶导数的达芬振子在谐和与宽带噪声联合激励下平稳响应的P分岔现象时，发现分数阶微分项的变化可以导致系统发生随机P分岔.申永军等[24-28]通过平均法对含分数阶微分的线性单自由度振子进行了动力学研究，提出了等效线性阻尼和等效线性刚度的概念，并对分数阶达芬振子的主共振和超谐共振进行了研究，得到了不同的分数阶微分项对系统响应特性的影响规律.温少芳等[29]在分数阶时滞反馈对达芬振子动力学特性的影响研究中发现，在分数阶和时滞耦合作用下，反馈系数和分数阶阶次等参数对系统动力学特性有着重要的影响，同时发现分数阶时滞反馈能同时起到位移时滞反馈和速度时滞反馈的作用.韦鹏等[30]对分数阶达芬振子的亚谐共振同样做出了相关的研究.除此之外，分数阶微积分在理论上的研究也有极大的应用价值，如Atanackovic等[31]研究了分数阶微积分欧拉-拉格朗日方程；王学彬等[32]以R-L型和Caputo型两种常见的分数阶导数定义为例，研究了如何利用拉普拉斯变换求解分数阶微分方程.陈明杰[33]提出了一种崭新的信号分析工具即分数阶傅里叶变换，并用经典的傅里叶变换的观点对分数阶傅里叶变换进行了充分的解释等等.
由于能够给出系统响应特性与参数之间的直接对应关系，所以解析研究在分数阶微积分的研究中占据着重要的地位.目前，在大量关于分数阶微积分近似解析解的研究文献中，大多数是将分数阶微分项当做系统的特殊阻尼来处理，这是不准确的.前述关于分数阶达芬振子的研究大多以单一谐波激励为主，较少涉及到多频激励.而在实际工程中，多频激励的发生并不少见.而且，多频激励下非线性系统具有更加丰富的现象，尤其是联合共振较为突出.Nayfeh的专著[34]以整数阶达芬振子的联合共振为例，简单介绍了其中的复杂动力学现象.本文以多频激励下含分数阶微分项的达芬振子为对象，研究了分数阶微分项对达芬振子3次超谐与1/3次亚谐联合共振时动力学特性的影响.采用平均法求解该系统的一阶近似解析解并通过数值解验证了结果的正确性.通过等效线性阻尼与等效线性刚度的概念，研究了分数阶微分项在不同的系数和阶次下对系统动力学行为尤其是多值性和跳跃现象的影响.
1. 分数阶达芬振子的3次超谐与1/3次亚谐联合共振的一阶近似解析解
研究如下含有分数阶微分项的达芬振子模型
$$\ddot {x}(t) + kx(t) + c\dot {x}(t) + \alpha _1 x^3(t) + K_1 D^p\left[x(t) \right] = \\ \qquad F_1 \cos \left( \omega _1 t + \beta _1 \right) + F_2 \cos \left( \omega _2 t + \beta _2 \right)$$ (1) 式中， $m, k, c, \alpha _1$ 分别为系统的质量、线性刚度、线性阻尼、非线性刚度系数； $F_1, \omega _1, \beta _1 , F_2, \omega _2, \beta _2 $ 分别为两项简谐激振力的幅值、频率及相位； $K_1$ $(K_1 > 0)$ 和 $p \ (0 \leqslant p \leqslant 1)$ 分别是分数阶微分项的系数和阶次.分数阶微分的定义有多种，这里采用Caputo形式
$$D^p\left[{x(t)} \right] = \dfrac{1}{\Gamma \left( {1 - p} \right)} \int_0^t \dfrac{\dot {x}(u)}{\left( {t - u} \right)^p}{\rm{d}}u$$ (2) 其中 $\Gamma (z)$ 为伽马函数，具有迭代性质 $\Gamma ( z+1 )=z\Gamma(z)$ .
对系统进行如下变量代换
$$\omega _0 = \sqrt {\dfrac{k}{m}} , \ \ 2\varepsilon \mu = \dfrac{c}{m} , \ \ \varepsilon \alpha = \dfrac{\alpha _1 }{m} , \ \ \varepsilon k_1 = \dfrac{K_1 }{m} \\ F_a = \dfrac{F_1 }{m} , \ \ F_b = \dfrac{F_2 }{m}$$ 式 (1) 可变为
$$\ddot x(t)+\omega^2_0x(t)+2\varepsilon\mu\dot x(t)+\varepsilon \alpha x^3(t)+\varepsilon k_1D^p [x(t)]=\\ \qquad F_a\cos (\omega_1t+\beta_1)+F_b \cos (\omega_2 t+\beta_2)$$ (3) 式(3) 从形式上满足了利用平均法求解系统近似解析解的要求.
研究3次超谐与1/3次亚谐联合共振，也就是系统的两个外部激励频率分别接近固有频率的1/3倍和3倍时的情况，即 $\omega _1 \approx \omega _0/3$ , $\omega _2 \approx 3\omega _0 $ .
为了定量地表示这两个近似等式的接近程度，引入调谐参数
$$ 9\omega _1^2 = \omega _0^2 + \varepsilon \sigma _1 $$ (4a) $$ \dfrac{1}{9}\omega _2^2 = \omega _0^2 + \varepsilon \sigma _2 $$ (4b) 因此
$$ \omega _1^2 = \dfrac{1}{81}\omega _2^2 + \varepsilon \left( {\sigma _1 - \sigma _2 } \right) $$ (4c) 将式 (4a) 代入式 (3) 可得
$$\ddot {x}(t) + 9\omega _1^2 x(t) = \\ \qquad \varepsilon \sigma _1 x(t) - 2\varepsilon \mu \dot {x}(t) - \varepsilon \alpha x^3(t) - \varepsilon k_1 D^p\left[{x(t)} \right] +\\ \qquad F_a \cos \left( {\omega _1 t + \beta _1 } \right) + F_b \cos \left( {\omega _2 t + \beta _2 } \right)$$ (5) 式(5) 可改写为
$$\ddot {x}(t) + 9\omega _1^2 x(t) = \varepsilon \left[{P_1 \left( {x, \dot {x}} \right) + P_2 \left( {x, \dot {x}} \right)} \right] + \\ \qquad F_a \cos \left( {\omega _1 t + \beta _1 } \right) + F_b \cos \left( {\omega _2 t + \beta _2 } \right)$$ (6) 其中
$$P_1 \left( {x, \dot {x}} \right) = \sigma _1 x(t) - 2\mu \dot {x}(t) - \alpha x^3(t) $$ (7a) $$P_2 \left( {x, \dot {x}} \right) = - k_1 D^p\left[{x(t)} \right] $$ (7b) 假设系统的解具有如下形式
$$x = a\cos \varphi + B_1 \cos \left( {\omega _1 t + \beta _1 } \right) + \\ \qquad B_2 \cos \left( {9\omega _1 t + \beta _2 } \right)$$ (8a) $$\dot {x} = - 3\omega _1 a\sin \varphi - B_1 \omega _1 \sin \left( {\omega _1 t + \beta _1 } \right) - \\ \qquad 9\omega _1 B_2 \sin \left( {9\omega _1 + \beta _2 } \right)$$ (8a) 其中， $\varphi = 3\omega _1 t + \theta$ , $B_1 = \dfrac{F_a }{ \omega _0^2 - \omega _1^2 }$ , $B_2 = \dfrac{F_b }{ \omega _0^2 - \omega _2^2 }$ .
将式(8) 中的解代入原系统，根据平均法可得
$$\dot {a} = - \dfrac{1}{3\omega _1 }\left[{P_1 \left( {a, \theta } \right) + P_2 \left( {a, \theta } \right)} \right]\sin \varphi $$ (9a) $$a\dot {\theta } = - \dfrac{1}{3\omega _1 }\left[{P_1 \left( {a, \theta } \right) + P_2 \left( {a, \theta } \right)} \right]\cos \varphi $$ (9b) 其中
$$P_1 ( a, \theta ) = \varepsilon \Bigg \{ \sigma _1 \Bigg [a\cos \varphi + B_1 \cos \left( {\dfrac{\varphi- \theta }{3} + \beta _1 } \right) + \\ \qquad B_2 \cos \left( {3\varphi -3\theta + \beta _2 } \right) \Bigg] - 2\mu \Bigg [- 3\omega _1 a\sin \varphi- \\ \qquad B_1 \omega _1 \sin \left( {\dfrac{\varphi - \theta }{3} + \beta _1 } \right) - 9B_2 \omega _1 \sin ( 3\varphi - \\ \qquad 3\theta + \beta _2 )\Bigg] -a \Bigg [a\cos \varphi + B_1 \cos \left( {\dfrac{\varphi- \theta }{3}} \right) + \\ \qquad B_2 \cos \left( {3\varphi -3\theta + \beta _2 } \right) \Bigg]^3 \Bigg \}$$ (10a) $$P_2 \left( {a, \theta } \right) = - \varepsilon k_1 D^p\Bigg[a\cos \varphi + B_1 \cos \left( {\dfrac{\varphi- \theta }{3} + \beta _1 } \right) + \\ \qquad B_2 \cos \left( {3\varphi -3\theta + \beta _2 } \right) \Bigg]$$ (10b) 采用平均法对式(9) 在区间 $\left[{0 , T} \right]$ 上进行积分平均，可以得到有关振幅 $a$ 和相位θ的平均方程
$$\dot {a} = - \dfrac{1}{3\omega _1 T}\int_0^T {\left[{P_1 \left( {a, \theta } \right) + P_2 \left( {a, \theta } \right)} \right]} \sin \varphi {\rm{d}}\varphi $$ (11a) $$a\dot {\theta } = - \dfrac{1}{3\omega _1 T}\int_0^T {\left[{P_1 \left( {a, \theta } \right) + P_2 \left( {a, \theta } \right)} \right]} \cos \varphi {\rm{d}}\varphi $$ (11b) 其中积分上限 $T$ 可以取 $ 2 \pi $ (若 $P_i \left( {a, \theta } \right) \ \left( {i = 1, 2} \right)$ 是周期函数)，或者 $T = \infty $ (若 $P_i \left( {a, \theta } \right) \ \left( {i = 1, 2} \right)$ 是非周期函数).
因此，式(11) 中第1部分的积分为
$$\dot {a}_1 = - \dfrac{1}{6\omega _1 \pi }\int_0^{2 \pi } P_1 \left( {a, \theta } \right)\sin \varphi {\rm{d}}\varphi =\\ \qquad - \varepsilon \mu a + \dfrac{\varepsilon aB_1^3 \sin \left( {\theta - 3\beta _1 } \right)}{24\omega _1 } + \dfrac{\varepsilon aB_2 a^2\sin \left( {3\theta - \beta _2 } \right)}{8\omega _1 }$$ (12a) $$a\dot {\theta }_1 = - \dfrac{1}{6\omega _1 \pi }\int_0^{2\pi } {P_1 \left( {a, \theta } \right)\cos \varphi {\rm{d}}\varphi } =\\ \qquad \dfrac{\varepsilon \alpha a^3}{8\omega _1 } + \dfrac{\varepsilon \alpha a\left( {B_1^2 + B_2^2 } \right)}{4\omega _1 } + \dfrac{\varepsilon \alpha B_1^3 \cos \left( {\theta - 3\beta _1 } \right)}{24\omega _1 } +\\ \qquad \dfrac{\varepsilon \alpha B_2 a^2\cos \left( {3\theta - \beta _2 } \right)}{8\omega _1 } - \dfrac{\varepsilon a\sigma _1 }{6\omega _1 }$$ (12b) 计算式(11) 中第2部分的积分
$$\dot {a}_2 = \mathop {\lim }\limits_{T \to \infty } - \dfrac{1}{3\omega _1 T}\int_0^T P_2 \left( {a, \theta } \right)\sin \varphi {\rm{d}}\varphi =\\ \qquad \mathop {\lim }\limits_{T \to \infty } \dfrac{\varepsilon k_1 }{3\omega _1 T}\int_0^T D^p\Bigg [a\cos \varphi + B_1 \cos \left( {\dfrac{\varphi - \theta }{3} + \beta _1 } \right) +\\ \qquad B_2 \cos \left( {3\varphi -3\theta + \beta _2 } \right) \Bigg]\sin \varphi {\rm{d}}\varphi$$ (13a) $$a\dot {\theta }_2 = \mathop {\lim }\limits_{T \to \infty } - \dfrac{1}{3\omega _1 T}\int_0^T P_2 \left( {a, \theta } \right)\cos \varphi {\rm{d}}\varphi =\\ \qquad \mathop {\lim }\limits_{T \to \infty } \dfrac{\varepsilon k_1 }{3\omega _1 T}\int_0^T D^p\Bigg [a\cos \varphi + B_1 \cos \left( {\dfrac{\varphi- \theta }{3} + \beta _1 } \right) +\\ \qquad B_2 \cos \left( {3\varphi -3\theta + \beta _2 } \right) \Bigg]\cos \varphi {\rm{d}}\varphi$$ (13b) 引入两个基本公式
$$ B_3 = \mathop {\lim }\limits_{T \to \infty } \int_0^T {\dfrac{\sin \left( {\omega t} \right)}{t^p}}{\rm{d}}t = \omega ^{p - 1}\Gamma \left( {1 - p} \right)\cos \left( {\dfrac{p\pi }{2}} \right) $$ (14a) $$ B_4 = \mathop {\lim }\limits_{T \to \infty } \int_0^T {\dfrac{\cos \left( {\omega t} \right)}{t^p}}{\rm{d}}t = \omega ^{p - 1}\Gamma \left( {1 - p} \right)\sin \left( {\dfrac{p\pi }{2}} \right) $$ (14b) 将 $\dot {a}_2 $ 分为3部分，根据式(2) 的Caputo定义形式，其中第1项为
$$\dot {a}_{21} = \mathop {\lim }\limits_{T \to \infty } \dfrac{\varepsilon k_1 }{3\omega _1 T}\int_0^T {D^p\left[{a\cos \varphi } \right]} \sin \varphi {\rm{d}}\varphi =\\ \qquad - \dfrac{\varepsilon ak_1 }{\Gamma \left( {1 - p} \right)T}\mathop {\lim }\limits_{T \to \infty } \int_0^T {\left[{\int_0^t {\dfrac{\sin \left( {3\omega _1 u + \theta } \right)}{\left( {t-u} \right)^p}{\rm{d}}u} } \right]} \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t$$ (15) 利用坐标变换 $s = t - u$ , ${\rm{d}} s = -{\rm{d}}u$ 得到
$$\dot {a}_{21} = - \dfrac{\varepsilon ak_1 }{\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T \cdot \\ \qquad \left[{\int_0^t {\dfrac{\sin \left( {3\omega _1 t + \theta-3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right] \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t =\\ \qquad - \dfrac{\varepsilon ak_1 }{\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\cdot \\ \qquad \left\{ {\int_0^T {\left[{\int_0^t {\dfrac{\sin \left( {3\omega _1 t + \theta } \right)\cos \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} } \right. \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t - \\ \qquad \int_0^T {\left[{\int_0^t {\dfrac{\cos \left( {3\omega _1 t + \theta } \right)\sin \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t \Bigg\}$$ (16) 再将 $\dot {a}_{21} $ 分为两部分，前半部分为 $\dot {a}_{211} $
$$\dot {a}_{211} = - \dfrac{\varepsilon ak_1 }{\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\left[ {\int_0^t {\dfrac{\cos \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad \sin ^2\left( {3\omega _1 t + \theta } \right){\rm{d}}t =\\ \qquad - \dfrac{\varepsilon ak_1 }{2\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\left[ {\int_0^t {\dfrac{\cos \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad \left[{1-\cos \left( {6\omega _1 t + 2\theta } \right)} \right]{\rm{d}}t =\\ \qquad - \dfrac{\varepsilon ak_1 }{12\omega _1 \Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\left[ {\int_0^t {\dfrac{\cos \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad {\rm{d}}\left[{6\omega _1 t-\sin \left( {6\omega _1 t + 2\theta } \right)} \right]$$ (17) 对上式进行分部积分可以得到
$$\dot {a}_{211} = - \dfrac{\varepsilon ak_1 }{12\omega _1 \Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^t {\dfrac{\cos \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} \cdot \\ \qquad \left[{6\omega _1 t-\sin \left( {6\omega _1 t + 2\theta } \right)} \right]\left| {\begin{array}{l} T \\ 0 \end{array}} \right. + \dfrac{\varepsilon ak_1 }{12\omega _1 \Gamma \left( {1 - p} \right)}\cdot \\ \qquad \mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\dfrac{\cos \left( {3\omega _1 t} \right) \left[{6\omega _1 t-\sin \left( {6\omega _1 t + 2\theta } \right)} \right]}{t^p}}{\rm{d}}t$$ (18) 利用极限的性质， $\dot {a}_{211} $ 式中后半部为0.再将式(14b)代入上式，得到
$$\dot {a}_{211} = - \dfrac{\varepsilon ak_1 }{2\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \left[{\int_0^t {\dfrac{\cos \left( {3\omega _1 s} \right)}{s^p}ds} } \right] =\\ \qquad - \dfrac{\varepsilon ak_1 \left( {3\omega _1 } \right)^{p - 1}}{2}\sin \left( {\dfrac{p\pi }{2}} \right)$$ (19) 计算 $\dot {a}_{21} $ 的后半部分 $\dot {a}_{212} $
$$ \dot {a}_{212} = \dfrac{\varepsilon ak_1 }{\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T} \cdot \\ \qquad \int_0^T {\left[{\int_0^t {\dfrac{\cos \left( {3\omega _1 t + \theta } \right)\sin \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t =\\ \qquad \dfrac{\varepsilon ak_1 }{2\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\left[ {\int_0^t {\dfrac{\sin \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad \sin \left( {6\omega _1 t + 2\theta } \right){\rm{d}}t = \\ \qquad - \dfrac{\varepsilon ak_1 }{12\omega _1 \Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\int_0^T {\left[ {\int_0^t {\dfrac{\sin \left( {3\omega _1 s} \right)}{s^p}{\rm{d}}s} } \right]} \cdot \\ \qquad {\rm{d}}\cos \left[{6\omega _1 t + 2\theta } \right]$$ (20) 对式(20) 进行分部积分可以得到
$$\dot {a}_{212} = - \dfrac{\varepsilon ak_1 }{12\omega _1 \Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\Bigg[\int_0^t {\dfrac{\sin \left( {3\omega _1 s} \right)}{s^p}}{\rm{d}}s \cdot \\ \qquad \cos ( 6 \omega _1 t + 2\theta )\left| {\begin{array}{l} T \\ 0 \end{array}} \right. - \int_0^T (\sin \left( {3\omega _1 t} \right)\cdot \\ \qquad \cos \left( {6\omega _1 t + 2\theta } \right) ) \Big / t^p{\rm{d}}t \Bigg]$$ (21) 根据极限的性质，当 $T \to \infty $ 时 $\dot {a}_{212} = 0$ .
下面计算 $\dot {a}_2 $ 中的第2项，根据上述分析过程可得
$$\dot {a}_{22} = \mathop {\lim }\limits_{T \to \infty } \dfrac{\varepsilon k_1 }{3\omega _1 T}\int_0^T {D^p\left[{B_1 \cos \left( {\dfrac{\varphi-\theta }{3} + \beta _1 } \right)} \right]} \sin \varphi {\rm{d}}\varphi = \\ \qquad \mathop {\lim }\limits_{T \to \infty } - \dfrac{\varepsilon k_1 B_1 }{\Gamma \left( {1 - p} \right)T}\int_0^T {\left[{\int_0^t {\dfrac{\sin \left( {\omega _1 u + \beta _1 } \right)}{\left( {t-u} \right)^p}{\rm{d}}u} } \right]} \cdot \\ \qquad \sin \left( {3\omega _1 t + \theta } \right){\rm{d}}t$$ (22) 对上式进行坐标变换及分部积分可以得到
$$\dot {a}_{22} = - \dfrac{\varepsilon k_1 B_1 }{\Gamma \left( {1 - p} \right)}\mathop {\lim }\limits_{T \to \infty } \dfrac{1}{T}\Bigg\{ \int_0^t {\dfrac{\cos \left( {\omega _1 s} \right)}{s^p}}{\rm{d}}s\cdot \\ \qquad \dfrac{ \cos \left( {\beta _1 - \theta - \omega _1 t} \right) - \cos \left( {\omega _1 t} \right)\cos \left( {\beta _1 + \theta + 2\omega _1 t} \right) } {2\omega _1 } \cdot \\ \qquad \sin \left( {\omega _1 t} \right)\left| {\begin{array}{l} T \\ 0 \end{array}} \right. - \int_0^T \Big\{ \left[{\cos \left( {\beta _1- \theta- \omega _1 t} \right)- } \right. \\ \qquad \left. {\cos ^2\left( {\omega _1 t} \right) \cos \left( {\beta _1 + \theta + 2\omega _1 t} \right)} \right]\sin \left( {\omega _1 t} \right) \Big \} \Big / \sqrt{2\omega _1 t^p}{\rm{d}}t\Bigg\}$$ (23) 根据极限的性质，当 $T \to \infty $ 时 $\dot {a}_{22} = 0$ .
同理可知
$$ \dot {a}_{23} = \mathop {\lim }\limits_{T \to \infty } \dfrac{\varepsilon k_1 }{3\omega _1 T}\int_0^T {D^p\left[{B_2 \cos \left( {3\varphi-3\theta + \beta _2 } \right)} \right]}\cdot \\ \qquad \sin \varphi {\rm{d}}\varphi = 0$$ (24) 因此
$${a}_2 = - \dfrac{\varepsilon ak_1 \left( {3\omega _1 } \right)^{p - 1}}{2}\sin \left( {\dfrac{p\pi }{2}} \right)$$ (25a) 同理，采用类似的推导过程可以得到
$$a\dot {\theta }_2 = \dfrac{\varepsilon ak_1 \left( {3\omega _1 } \right)^{p - 1}}{2}\cos \left( {\dfrac{p\pi }{2}} \right)$$ (25b) 联立式(25) 和式(12)，可以得到
$$ \dot {a} = - \varepsilon \mu a - \dfrac{\varepsilon ak_1 \left( {3\omega _1 } \right)^{p - 1}}{2}\sin \left( {\dfrac{p\pi }{2}} \right) + \dfrac{\varepsilon \alpha B_1^3 \sin \left( {\theta - 3\beta _1 } \right)}{24\omega _1 } + \\ \qquad \dfrac{\varepsilon \alpha B_2 a^2\sin \left( {3\theta - \beta _2 } \right)}{8\omega _1 }$$ (26a) $$ a\dot {\theta } = \dfrac{\varepsilon ak_1 \left( {3\omega _1 } \right)^{p - 1}}{2}\cos \left( {\dfrac{p\pi }{2}} \right) + \dfrac{\varepsilon \alpha a^3}{8\omega _1 } + \\ \qquad \dfrac{\varepsilon \alpha a\left( {B_1^2 + B_2^2 } \right)}{4\omega _1 } + \dfrac{\varepsilon \alpha B_1^3 \cos \left( {\theta - 3\beta _1 } \right)}{24\omega _1 } + \\ \qquad \dfrac{\varepsilon \alpha B_2 a^2\cos \left( {3\theta - \beta _2 } \right)}{8\omega _1 } - \dfrac{\varepsilon a\sigma _1 }{6\omega _1 }$$ (26b) 将原参数代入上式，则有
$$ \dot {a} = - \dfrac{a}{2m}C\left( p \right) + \dfrac{\alpha _1 B_1^3 \sin \left( {\theta - 3\beta _1 } \right)}{24\omega _1 m} + \\ \qquad \dfrac{\alpha _1 B_2 a^2\sin \left( {3\theta - \beta _2 } \right)}{8\omega _1 m}$$ (27a) $$ a\dot {\theta } = - \dfrac{3a\omega _1 }{2} + \dfrac{a}{6m\omega _1 }K\left( p \right) + \dfrac{\alpha _1 }{8\omega _1 m}\left[{a^3 + } \right. \\ \qquad \left. {2a\left( {B_1^2 + B_2^2 } \right)} \right] + \dfrac{\alpha _1 B_1^3 \cos \left( {\theta - 3\beta _1 } \right)}{24\omega _1 m} + \\ \qquad \dfrac{\alpha _1 B_2 a^2\cos \left( {3\theta - \beta _2 } \right)}{8\omega _1 m}$$ (27b) 式中
$$ C\left( p \right) = c + K_1 \left( {3\omega _1 } \right)^{p - 1}\sin \left( {\dfrac{p\pi }{2}} \right) $$ (28a) $$ K\left( p \right) = k + K_1 \left( {3\omega _1 } \right)^p\cos \left( {\dfrac{p\pi }{2}} \right) $$ (28b) 上式 $C( p )$ 和 $K( p )$ 分别定义为3次超谐与1/3次亚谐联合共振时的等效线性阻尼及等效线性刚度.
分析式(28) 可知，分数阶项的系数与阶次均出现在等效线性阻尼与等效线性刚度中.这与文献[20]所提出的分数阶微分项不仅起着阻尼的作用，同时还起着线性刚度的作用这一结论相一致.它的系数 $K_1 $ 和阶次 $p$ 对于上式的等效线性阻尼和等效线性刚度有着重要的影响.当 $K_1 \to 0$ 时，系统将会由分数阶趋近于传统的整数阶形式，并且与等效线性阻尼和等效线性刚度呈线性关系，影响着系统的共振频率和振幅的大小.当 $K_1 $ 增大时，等效线性阻尼与等效线性刚度都会随之增大，因此系统的响应振幅将会减小，共振频率也会随之改变；反之亦然.同时，分数阶微分项的阶次 $p$ 也对该系统有重要的影响.当 $p \to 1$ 时，分数阶微分项几乎等同于线性阻尼达到极大值 $c + K_1 $ ，系统的响应幅值将降到最低；当 $p \to 0$ 时，分数阶微分项几乎等同于线性刚度达到极大值 $k + K_1 $ ，系统的共振频率增大，响应振幅相应地减小.不难看出等效线性阻尼与等效线性刚度与激励的频率也存在着一定的关系.
2. 联合共振定常解
非线性系统的稳态振动在工程实际情况中具有较大的应用价值，因此这里研究3次超谐与1/3次亚谐联合共振情况下的定常解(即稳态解).令式(27) 中的 $\dot {a} = 0$ 和 $a\dot {\theta } = 0$ 则有
$$\dfrac{\alpha _1 B_1^3 \sin \left( {\theta - 3\beta _1 } \right)}{24\omega _1 m} + \dfrac{\alpha _1 B_2 a^2\sin \left( {3\theta - \beta _2 } \right)}{8\omega _1 m} = \dfrac{a}{2m}C(p)$$ (29a) $$ \dfrac{\alpha _1 B_1^3 \cos \left( {\theta - 3\beta _1 } \right)}{24\omega _1 m} + \dfrac{\alpha _1 B_2 a^2\cos \left( {3\theta - \beta _2 } \right)}{8\omega _1 m} = \dfrac{3a\omega _1 }{2} - \\ \qquad \dfrac{a}{6m\omega _1 }K\left( p \right) - \dfrac{\alpha _1 }{8\omega _1 m}\left[{a^3 + 2a\left( {B_1^2 + B_2^2 } \right)} \right]$$ (29b) 为了简化计算过程，且不影响对结果的分析，现假设 $\beta _1 = \beta _2 = 0$ .
令
$$H_1 = \dfrac{\alpha _1 B_1^3 }{24\omega _1 m} , \ \ H_2 = \dfrac{\alpha _1 B_2 }{8\omega _1 m} , \ \ H_3 = \dfrac{C\left( P \right)}{2m} \\ H_4 = \dfrac{\alpha _1 }{8\omega _1 m} , \ \ H_5 = \dfrac{3\omega _1 }{2} - \dfrac{K\left( P \right)}{6m\omega _1 } - \dfrac{\alpha _1 \left( {B_1^2 + B_2^2 } \right)}{4\omega _1 m}$$ 式(29) 可简化成
$$H_1 \sin \theta + H_2 a^2\sin 3\theta = H_3 a $$ (30a) $$H_1 \cos \theta + H_2 a^2\cos 3\theta = - H_4 a^3 + H_5 a $$ (30b) 现将式(30a)与式(30b)平方相加再利用2倍角公式可得
$$\sin \theta = \Big [\Big (- H_3^2 a^2- H_5^2 a^2 + 2H_4 H_5 a^4- H_4^2 a^6 + \\ \qquad H_1^2 + H_2^2 a^4 + 2H_1 H_2 a^2\Big ) \Big / \Big (4H_1 H_2 a^2 \Big ) \Big]^{1/2}$$ (31) 利用三倍角公式，式(30a)可写成
$$\left( {H_1 + 3H_2 a^2} \right)\sin \theta - 4H_2 a^2\sin ^3 \theta = H_3 a$$ (32) 将式(31) 代入上式平方后可消去θ，经整理得到定常解的幅频曲线方程
$$\left( {H_1 + 3H_2 a^2} \right)^2 \left( {\dfrac{ - H_3^2 a^2 - H_5^2 a^2 + 2H_4 H_5 a^4 - H_4^2 a^6}{4H_1 H_2 a^2} + } \right. \\ \qquad \left. {\dfrac{H_1^2 + H_2^2 a^4 + 2H_1 H_2 a^2}{4H_1 H_2 a^2}} \right) - 8H_2 a^2 \left( {H_1 + 3H_2 a^2} \right)\cdot \\ \qquad \left( {\dfrac{ - H_3^2 a^2 - H_5^2 a^2 + 2H_4 H_5 a^4 - H_4^2 a^6 + H_1^2 + H_2^2 a^4}{4H_1 H_2 a^2} + } \right. \\ \qquad \left. {\dfrac{2H_1 H_2 a^2}{4H_1 H_2 a^2}} \right)^2 + \left( {4H_2 a^2} \right)^2 \left( {\dfrac{ - H_3^2 a^2 - H_5^2 a^2 + 2H_4 H_5 a^4}{4H_1 H_2 a^2}} \right. - \\ \qquad \left. {\dfrac{H_4^2 a^6 + H_1^2 + H_2^2 a^4 + 2H_1 H_2 a^2}{4H_1 H_2 a^2}} \right)^3 = \left( {H_3 a} \right)^2$$ (33) 及相频曲线方程
$$\theta = \arcsin \Big \{ \Big [\Big (-H_3^2 a^2-H_5^2 a^2 + 2H_4 H_5 a^4-H_4^2 a^6 + \\ \qquad H_1^2 + H_2^2 a^4 + 2H_1 H_2 a^2\Big ) \Big / \Big (4H_1 H_2 a^2 \Big ) \Big] ^{1/2}\Big \}$$ (34) 令 $a^2 = \bar {a}$ ，将式(33) 改写成关于 $\bar {a}$ 的多项式
$$A_1 \bar {a} ^7 + B_1 \bar {a} ^6 + C_1 \bar {a} ^5 + D_1 \bar {a} ^4 + E_1 \bar {a} ^3 + \\ \qquad F_1 \bar {a} ^2 + G_1 \bar {a} + K_1 = 0$$ (35) 至此，得到了关于 $\bar {a}$ 的多项式方程.考虑到多项式系数的表达式非常复杂(见附录)，所有推导过程可用Mathematica软件方便完成.
3. 数值仿真
现对式(1) 所表示的系统进行数值求解，选取基本系统参数 $m = 5$ , $k = 10$ , $c = 0.2$ , $\alpha _1 = 3$ , $F_1 = 10$ , $F_2 = 1$ , $K_1 = 0.2$ , $p = 0.4$ , 根据式(8) 和式(35) 得到该系统的幅频曲线，如图 1中实线所示，图中横轴为激励频率 $\omega $ ，纵轴为响应振幅 $a$ .为了验证解析结果的正确性，参照文献[2]中介绍的分数阶幂级数法，对该分数阶达芬系统进行数值分析.首先，引入数值解法的近似公式
$$D^p\left[{x\left( {t_l } \right)} \right] \approx h^{ - p}\sum\limits_{j = 0}^l {C_j^p} x\left( {t_{l - j} } \right)$$ (36) ![]()
图 1 系统近似解与数值解幅频曲线的比较Figure 1. The comparison of the amplitude-frequency curve by the approximately analytical solution with that by numerical integration其中 $t_l = lh$ 为时间采样点， $h$ 为时间步长， $C_j^p $ 为分数阶二项式系数，具有如下迭代关系
$$C_0^p = 1 , \quad C_j^p = \left( {1 - \dfrac{1 + p}{j}} \right)C_{j - 1}^p$$ (37) 采用Matlab软件进行数值计算，选取步长 $h = 0.000 5$ ，计算时间为400 s，将前320 s响应值略去，取后80 s的响应最大幅值为稳定幅值，所得数值计算结果如图 1中圆圈所示.可见，两者符合较好，说明本文得到的结果具有较高的精度.
选取一组基本系统参数 $m = 5$ , $k = 45$ , $c = 0.2$ , $\alpha _1 = 45$ , $F_1 = 40$ , $F_2 = 2$ , $K_1 = 1$ , $p = 0.5$ , 根据式(8) 和式(35) 得到该系统的幅频曲线，如图 2所示.图 2中整数阶曲线为 $K_1 = 0$ 时的情况.显然，在该系统发生3次超谐与1/3次亚谐联合共振的情况下，当第一激励 $F_1 $ 大于第二激励 $F_2 $ 数倍时，该系统出现了与分数阶达芬振子超谐共振相似的情况，表现出3次超谐共振[35]的特性.即分数阶微分项引起了该系统等效线性阻尼和等效线性刚度的增大，使该分数阶达芬振子幅频曲线中的振幅减小，幅频曲线的弯曲程度增加(即改变了系统的频率特性).
![]()
图 2 传统整数阶与分数阶振子的幅频曲线比较Figure 2. Comparison of amplitude-frequency curves between the traditional integer-order and the fractional-order oscillators当分数阶微分项系数 $K_1 $ 选取不同值时，得到的幅频曲线如图 3所示.随着 $K_1 $ 的逐渐增大，该系统的幅值呈现出减小趋势.同时，该系统的等效线性刚度在逐渐增大，导致系统的共振频率增大.因此该系统的分数阶微分项系数 $K_1 $ 对系统响应的影响体现在两个方面，即共振振幅和共振频率.
![]()
图 3 不同分数阶系数下幅频曲线的比较Figure 3. Comparison of amplitude-frequency curves in different fractional coefficients分析不同分数阶阶次 $p$ 的影响，选取 $K_1 = 1$ , $p$ 分别取0.1, 0.25, 0.5, 0.75, 0.9时的幅频曲线如图 4所示.如图可知，随着 $p$ 的逐渐增大，导致等效线性阻尼也随之增大，系统的共振幅值随之减弱，固有频率也逐渐减小；当 $p$ 达到一定程度时系统的多解和跳跃现象将会消失.因此，分数阶微分项的阶次 $p$ 对系统的共振振幅和共振频率同样有重要影响.
![]()
图 4 不同分数阶阶次下幅频曲线的比较Figure 4. Comparison of amplitude-frequency curves in different fractional orders选取另一组基本参数 $m = 5$ , $k = 45$ , $c = 0.2$ , $\alpha _1 = 15$ , $F_1 = 2$ , $F_2 = 200$ , $K_1 = 1$ , $p = 0.5$ , 根据式(8) 和式(35) 便可以得到该系统的幅频曲线如图 5所示，图 5中整数阶曲线为 $K_1 = 0$ 时的情况.由图 5可见，当第一激励 $F_1 $ 小于第二激励 $F_2 $ 数倍时，该系统出现了与分数阶达芬振子亚谐共振相似的情况，表现出1/3次亚谐共振[35]的特性.即分数阶微分项引起了该系统等效线性阻尼和等效线性刚度的增大，使系统幅频曲线中的振幅减小及共振频率增大，幅频曲线向右偏移.
![]()
图 5 传统整数阶与分数阶振子的幅频曲线比较Figure 5. Comparison of amplitude-frequency curves of the traditional integer-order and the fractional-order oscillators选取 $p = 0.5$ , $K_1 $ 分别取0.5, 1, 1.5时的幅频曲线如图 6(a)所示； $K_1 = 1 $ , 并且 $p$ 分别取0, 0.5, 1时的幅频曲线如图 6(b)所示.分析图 6(a)可知，随着分数阶微分项系数 $K_1 $ 逐渐增大，导致系统等效线性阻尼及等效线性刚度也随之增大，因此系统的响应幅值逐渐减小、系统的共振频率逐渐增大，同时系统的共振区域逐渐减小且幅频曲线向右偏移.分析图 6(b)可知，随着分数阶微分项阶次 $p$ 逐渐增大，系统等效线性阻尼逐渐增大、等效线性刚度逐渐减小，因此系统的响应幅值及共振频率逐渐减小，同时系统的共振区域显著减小，并且在不同参数下幅频曲线出现了相交现象.通过上述分析可知，分数阶微分项系数与阶次对系统的动力学特性有着重要影响.
![]()
图 6 不同分数阶参数下幅频曲线的比较Figure 6. Comparison of amplitude-frequency curves under different fractional-order parameters改变系统的基本参数为 $m = 2$ , $k = 15$ , $c = 0.2$ , $\alpha _1 = 80$ , $F_1 = 3 000$ , $F_2 = 15 000$ , $K_1 = 0$ , 该系统为传统整数阶达芬振子的3次超谐与1/3次亚谐联合共振，图 7(a)和图 7(b)分别给出了联合共振下的幅频曲线.由图可知，这种联合共振最多可同时出现7个解，这与单个解存在时的情况有明显的不同.分析图中的多解情况可得：(1) 1个非平凡的稳定解；(2) 3个非平凡解，其中一个是不稳定的；(3) 5个非平凡解，其中两个是不稳定的；(4) 7个非平凡解(P1~P7)，其中3个是不稳定的.将这种联合共振与单独存在的3次超谐共振和1/3次亚谐共振进行比较后可以看出，图 7中B和C两枝类似于亚谐共振，而A和D两枝则类似于超谐共振[36].
图 8为3次超谐与1/3次亚谐联合共振在 $p = 0.5$ , $K_1 $ 分别取0, 10, 20情况下的幅频曲线.不难看出，这种出现多解现象的幅频曲线在改变分数阶微分项系数 $K_1 $ 时出现了和上述相似的结果：随着 $K_1 $ 的增大，导致系统的等效线性阻尼及等效线性刚度也随之增大，因此幅频曲线中的幅值减小，系统的共振频率增大.不仅如此，该系统的多值性也发生了一系列变化：随着 $K_1 $ 的逐渐增大，导致联合共振出现3个解的区域逐渐减小，5个解的发生区域将提前出现并且逐渐扩大，同时7个解的发生区域将向后推移.
![]()
图 8 不同分数阶系数 $K_1 $ 下的幅频曲线比较Figure 8. Comparison of amplitude-frequency curves in different fractional-order parameters图 9为3次超谐与1/3次亚谐联合共振在 $K_1 = 10$ ， $p$ 分别取0.25, 0.5, 1情况下的幅频曲线.如图分析可得，当分数阶微分项阶次 $p$ 逐渐增大到1时，系统的等效线性阻尼也随之增大，不仅降低了系统的幅值，同时对系统的多值性也有很大的影响.随着 $p$ 的逐渐增大，此联合共振出现3个解的区域不断减小，出现5个解的区域扩大并且7个解的发生区域后移.当 $p$ 增大到一定程度时系统的多解性消失.
![]()
图 9 不同分数阶阶次 $p$ 下的幅频曲线比较Figure 9. Comparison of amplitude-frequency curves in different fractional-order parameters4. 结论
本文利用平均法对分数阶达芬振子的3次超谐与1/3次亚谐联合共振的响应进行了研究，通过等效线性阻尼及等效线性刚度的概念对该系统在不同分数阶微分项系数 $K_1 $ 及阶次 $p$ 的情况下分别做出了分析，并与传统整数阶达芬振子进行了比较.结果发现，分数阶微分项阶次 $p$ 和系数 $K_1 $ 通过等效线性阻尼及等效线性刚度对系统起着刚度和阻尼的作用，不仅仅影响着系统的响应幅值、共振频率，同时还对系统的多值性个数、发生区域面积、发生先后等有重要影响，对系统动力学特性的研究具有重要意义.
$$A_1 = - \dfrac{H_4^6 }{4H_1^3 H_2 }\\ B_1 = \dfrac{3H_2 H_4^4 }{4H_1^3 } + \dfrac{3H_4^5 H_5 }{2H_1^3 H_2 },\\ C_1 = - \dfrac{3H_2^3 H_4^2 }{4H_1^3 } - \dfrac{3H_3^2 H_4^4 }{4H_1^3 H_2 } - \dfrac{3H_2 H_4^3 H_5 }{H_1^3 } - \dfrac{15H_4^4 H_5^2 }{4H_1^3 H_2 }\\ D_1 = \dfrac{H_2^5 }{4H_1^3 } + \dfrac{3H_2 H_3^2 H_4^2 }{2H_1^3 } + \dfrac{H_4^4 }{4H_1 H_2 } + \dfrac{3H_2^3 H_4 H_5 }{2H_1^3 } + \\ \qquad \dfrac{3H_3^2 H_4^3 H_5 }{H_1^3 H_2 } + \dfrac{9H_2 H_4^2 H_5^2 }{2H_1^3 } + \dfrac{5H_4^3 H_5^3 }{H_1^3 H_2 } \\ E_1 = - \dfrac{3H_2^3 H_3^2 }{4H_1^3 } + \dfrac{H_2 H_4^2 }{4H_1 } - \dfrac{3H_3^4 H_4^2 }{4H_1^3 H_2 } - \dfrac{3H_2 H_3^2 H_4 H_5 }{H_1^3 }-\\ \qquad \dfrac{H_4^3 H_5 }{H_1 H_2 } - \dfrac{3H_2^3 H_5^2 }{4H_1^3 } - \dfrac{9H_3^2 H_4^2 H_5^2 }{2H_1^3 H_2 } - \dfrac{3H_2 H_4 H_5^3 }{H_1^3 } -\\ \qquad \dfrac{15H_4^2 H_5^4 }{4H_1^3 H_2 }$$ $$F_1 = - \dfrac{H_2^3 }{2H_1 } + \dfrac{3H_2 H_3^4 }{4H_1^3 } + \dfrac{H_4^2 }{2} + \dfrac{H_3^2 H_4^2 }{2H_1 H_2 } - \dfrac{H_2 H_4 H_5 }{2H_1 } + \\ \qquad \dfrac{3H_3^4 H_4 H_5 }{2H_1^3 H_2 } + \dfrac{3H_2 H_3^2 H_5^2 }{2H_1^3 } + \dfrac{3H_4^2 H_5^2 }{2H_1 H_2 } + \dfrac{3H_3^2 H_4 H_5^3 }{H_1^3 H_2 } + \\ \qquad \dfrac{3H_2 H_5^4 }{4H_1^3 } + \dfrac{3H_4 H_5^5 }{2H_1^3 H_2 } \\ G_1 = \dfrac{H_2 H_3^2 }{4H_1 } - \dfrac{H_3^6 }{4H_1^3 H_2 } - H_4 H_5 - \dfrac{H_3^2 H_4 H_5 }{H_1 H_2 } + \dfrac{H_2 H_5^2 }{4H_1 } - \\ \qquad \dfrac{3H_3^4 H_5^2 }{4H_1^3 H_2 } - \dfrac{H_4 H_5^3 }{H_1 H_2 } - \dfrac{3H_3^2 H_5^4 }{4H_1^3 H_2 } - \dfrac{H_5^6 }{4H_1^3 H_2 } \\ K_1 = \dfrac{H_1 H_2 }{4} - \dfrac{H_3^2 }{2} + \dfrac{H_3^4 }{4H_1 H_2 } + \dfrac{H_5^2 }{2} + \dfrac{H_3^2 H_5^2 }{2H_1 H_2 } + \dfrac{H_5^4 }{4H_1 H_2 }$$ -
![]()
图 1 系统近似解与数值解幅频曲线的比较
Figure 1. The comparison of the amplitude-frequency curve by the approximately analytical solution with that by numerical integration
![]()
图 2 传统整数阶与分数阶振子的幅频曲线比较
Figure 2. Comparison of amplitude-frequency curves between the traditional integer-order and the fractional-order oscillators
![]()
图 3 不同分数阶系数下幅频曲线的比较
Figure 3. Comparison of amplitude-frequency curves in different fractional coefficients
![]()
图 4 不同分数阶阶次下幅频曲线的比较
Figure 4. Comparison of amplitude-frequency curves in different fractional orders
![]()
图 5 传统整数阶与分数阶振子的幅频曲线比较
Figure 5. Comparison of amplitude-frequency curves of the traditional integer-order and the fractional-order oscillators
![]()
图 6 不同分数阶参数下幅频曲线的比较
Figure 6. Comparison of amplitude-frequency curves under different fractional-order parameters
![]()
图 8 不同分数阶系数 $K_1 $ 下的幅频曲线比较
Figure 8. Comparison of amplitude-frequency curves in different fractional-order parameters
-
[1] Kilbas AA, Srivastava HM, Trujillo JJ. Theory and Applications of Fractional Differential Equations. Amsterdam:Elsevier, 2006
[2] Petras I. Fractional-order Nonlinear Systems:Modeling, Analysis and Simulation. Beijing:Higher Education Press, 2011:18-19
[3] Rossikhin YA, Shitikova MV. Application of fractional calculus for dynamic problems of solid mechanics:novel trends and recent results. Applied Mechanics Reviews, 2010, 63:010801-1-52 doi: 10.1115/1.4000563
[4] Machado JAT, Kiryakova V, Mainardi F. Recent history of fractional calculus. Communications in Nonlinear Science and Numerical Simulation, 2011, 16:1140-1153 doi: 10.1016/j.cnsns.2010.05.027
[5] 陈章耀, 王亚茗, 张春.双状态切换下BVP振子的复杂行为分析.力学学报, 2016, 48(4):953-962 http://lxxb.cstam.org.cn/CN/abstract/abstract145830.shtml Chen Zhangyao, Wang Yaming, Zhang Chun. Complicated behaviors as well as the mechanism in BVP oscillator with switches related to two states. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(4):953-962 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract145830.shtml
[6] 苏二龙, 罗建军.高超声速飞行器横侧向失稳非线性分岔分析.力学学报, 2016, 48(5):1192-1201 http://lxxb.cstam.org.cn/CN/abstract/abstract146016.shtml Su Erlong, Luo Jianjun. Nonlinear bifurcation analysis of lateral loss of stability for hypersonic vehicle. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5):1192-1201 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract146016.shtml
[7] 高雪, 陈前, 刘先斌.一类分段光滑隔振系统的非线性动力学设计方法.力学学报, 2016, 48(1):192-200 doi: 10.6052/0459-1879-15-099 Gao Xue, Chen Qian, Liu Xianbin. Nonlinear dynamics design for piecewise smooth vibration isolation system. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1):192-200 (in Chinese) doi: 10.6052/0459-1879-15-099
[8] 王浩宇, 吴勇军. 1:1内共振对随机振动系统可靠性的影响.力学学报, 2015, 47(5):807-813 doi: 10.6052/0459-1879-15-058 Wang Haoyu, Wu Yongjun. The influence of one-to-one internal resonance on reliability of random vibration system. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(5):807-813 (in Chinese) doi: 10.6052/0459-1879-15-058
[9] 李占龙, 孙大刚, 韩斌慧.基于分数阶导数的黏弹性减振系统时频特性.应用基础与工程科学学报, 2017, 25(1):187-198 http://www.cnki.com.cn/Article/CJFDTOTAL-YJGX201701017.htm Li Zhanlong, Sun Dagang, Han Binhui. Time and frequency features of viscoelastic vibration damping system based on fractional derivative. Journal of Basic Science and Engineering, 2017, 25(1):187-198 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YJGX201701017.htm
[10] 蔡伟, 陈文.复杂介质中任意阶频率依赖耗散声波的分数阶导数模型.力学学报, 2016, 48(6):1265-1280 http://lxxb.cstam.org.cn/CN/abstract/abstract146097.shtml Cai Wei, Chen Wen. Fractional derivative modeling of frequency-dependent dissipative mechanism for wave propagation in complex media. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(6):1265-1280 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract146097.shtml
[11] 韦鹏, 申永军, 杨绍普.分数阶van der Pol振子的超谐共振.物理学报, 2014, 63(1):47-58 http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201401006.htm Wei Peng, Shen Yongjun, Yang Shaopu. Super-harmonic resonance of fractional-order van der Pol oscillator. Acta Phys. Sin, 2014, 63(1):47-58 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201401006.htm
[12] 林世敏, 许传炬.分数阶微分方程的理论和数值方法研究.计算数学, 2016, 38(1):1-24 http://www.cnki.com.cn/Article/CJFDTOTAL-JSSX201601001.htm Lin Shimin, Xu Chuanju. Theoretical and numerical investigation of fractional differential equations. Mathematica Numerica Sinica, 2016, 38(1):1-24 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-JSSX201601001.htm
[13] 谭健, 周洲, 祝小平.飞翼布局无人机分数阶积分滑模姿态控制.控制理论与应用, 2015, 32(5):607-614 http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201505005.htm Tan Jian, Zhou Zhou, Zhu Xiaoping. Attitude control for flying wing unmanned aerial vehicles based on fractional order integral sliding-mode. Control Theory and Applications, 2015, 32(5):607-614 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201505005.htm
[14] 文家燕, 高远, 刘传国.永磁同步电动机的双闭环分数阶控制研究.微特电机, 2016, 44(1):34-38 http://www.cnki.com.cn/Article/CJFDTOTAL-WTDJ201601009.htm Wen Jiayan, Gao Yuan, Liu Chuanguo. Control of PMSM via a double closed-loop fractionalorder control strategy. Small and Special Electrical Machines, 2016, 44(1):34-38 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-WTDJ201601009.htm
[15] 牛江川, 申永军, 杨绍普.基于速度反馈分数阶PID控制的达芬振子的主共振.力学学报, 2016, 48(2):422-429 doi: 10.6052/0459-1879-15-332 Niu Jiangchuang, Shen Yongjun, Yang Shaopu. Primary resonance of Duffing oscillator with fractional-order PID controller based on velocity feedback. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2):422-429 (in Chinese) doi: 10.6052/0459-1879-15-332
[16] 张碧陶, 高福荣, 姚科.集成神经网络与自适应算法的分数阶滑模控制.控制理论与应用, 2016, 33(10):1373-1377 http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201610014.htm Zhang Bitao, Gao Furong, Yao Ke. Neural network and adaptive algorithm-based fractional order sliding mode controller. Control Theory and Applications, 2016, 33(10):1373-1377 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201610014.htm
[17] 孙会来, 金纯, 张文明.基于分数阶微积分的油气悬架建模与试验分.振动与冲击, 2014, 34(17):167-172 http://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201417031.htm Sun Huilai, Jin Chun, Zhang Wenming. Modeling and tests for a hydro-peneumatic suspension based on fractional calculus. Journal of Vibration and Shocck, 2014, 34(17):167-172 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201417031.htm
[18] 吴光强, 黄焕军, 叶光湖.基于分数阶微积分的汽车空气悬架半主动控制.农业机械学报, 2014, 45(7):19-25 doi: 10.6041/j.issn.1000-1298.2014.07.004 Wu Guangqiang, Huang Huanjun, Ye Guanghu. Semi-active control of automotive air suspension based on fractional calculus. Transactions of The Chinese Society of Agricultural Machinery, 2014, 45(7):19-25 (in Chinese) doi: 10.6041/j.issn.1000-1298.2014.07.004
[19] 陈丙三, 曾寿金, 江吉彬.磁流变阻尼器的分数阶模型及其减振系统分析.机械设计与制造, 2012 (7):219-221 http://www.cnki.com.cn/Article/CJFDTOTAL-JSYZ201207080.htm Chen Bingsan, Zeng Shoujin, Jiang Jibin. Fractional calculus modeling of the megnetorheological damper and the analysis of its damping system. Machinery Design and Manufacture, 2012 (7):219-221 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-JSYZ201207080.htm
[20] 马冬冬, 王志强, 王进君.单相逆变器分数阶建模及分析.电测与仪表, 2017, 54(6):106-112 http://www.cnki.com.cn/Article/CJFDTOTAL-DCYQ201706018.htm Ma Dongdong, Wang Zhiqiang, Wang Jinjun. Fractional order modeling and analysis of single phase inverter. Electrical Measurement & Instrumentation, 2017, 54(6):106-112 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-DCYQ201706018.htm
[21] 何志磊, 朱珍德, 朱明礼.基于分数阶导数的非定常蠕变本构模型研究.岩土力学, 2016, 37(3):737-744 http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201603017.htm He Zhilei, Zhu Zhende, Zhu Mingli. An unsteady creep constitutive model based on fractional order derivatives. Rock and Soil Mechanics, 2016, 37(3):737-744 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201603017.htm
[22] 段晓梦, 殷德顺, 安丽媛.基于分数阶微积分的黏弹性材料变形研究.中国科学:物理学·力学·天文学, 2013, 43(8):971-977 http://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201308012.htm Duan Xiaomeng, Yin Deshun, An Liyuan. The deformation study in viscoelastic materials based on fractional order calculus. Scientia Sinica Physica, Mechanica & Astronomica, 2013, 43(8):971-977 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201308012.htm
[23] 陈林聪, 朱位秋.谐和与宽带噪声联合激励下含分数导数型阻尼的达芬振子的平稳响应.应用力学学报, 2010, 27(1):517-521 http://www.cnki.com.cn/Article/CJFDTOTAL-YYLX201003018.htm Chen Lincong, Zhu Weiqiu. Stationary response of Duffing oscillator with fractionalderivative damping under combined harmonic and wide band noise excitations. Journal of Applied Mechanics, 2010, 27(1):517-521 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YYLX201003018.htm
[24] 申永军, 杨绍普, 邢海军.含分数阶微分的线性单自由度振子的动力学分析.物理学报, 2012, 61(11):158-163 http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201211022.htm Shen Yongjun, Yang Shaopu, Xing Haijun. Dynamical analysis of linear single degreeof-freedom oscillator with fractional-order derivative. Acta Phys Sin, 2012, 61(11):158-163 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201211022.htm
[25] 申永军, 杨绍普, 邢海军.含分数阶微分的线性单自由度振子的动力学分析(Ⅱ).物理学报, 2012, 61(15):55-63 http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201211022.htm Shen Yongjun, Yang Shaopu, Xing Haijun. Dynamical analysis of linear SDOF oscillator with fractional-order derivative(Ⅱ). Acta Phys Sin, 2012, 61(15):55-63 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201211022.htm
[26] Shen YJ, Yang SP, Xing HJ, et al. Primary resonance of Duffing oscillator with fractional-order derivative. Communications in Nonlinear Science and Numerical Simulation, 2012, 17(7):3092-3100 doi: 10.1016/j.cnsns.2011.11.024
[27] Shen Yongjun, Yang Shaopu, Xing Haijun, et al. Primary resonance of Duffing oscillator with two kinds of fractional-order derivatives. International Journal of Non-Linear Mechanics, 2012, 47(9):975-983 doi: 10.1016/j.ijnonlinmec.2012.06.012
[28] 申永军, 杨绍普, 邢海军.分数阶Duffng振子的超谐共振.力学学报, 2012, 44(4):762-768 http://lxxb.cstam.org.cn/CN/abstract/abstract143486.shtml Shen Yongjun, Yang Shaopu, Xing Haijun. Super-harmonic resonance of fractional-order Duffing oscillator. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(4):762-768 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract143486.shtml
[29] 温少芳, 申永军, 杨绍普.分数阶时滞反馈对达芬振子动力学特性的影响.物理学报, 2016, 65(9):166-175 http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201609022.htm Wen Shaofang, Shen Yongjun, Yang Shaopu. Dynamical analysis of Duffing oscillator with fractional-order feedback with time delay. Acta Phys. Sin, 2016, 65(9):166-175 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-WLXB201609022.htm
[30] 韦鹏, 申永军, 杨绍普.分数阶达芬振子的亚谐共振.振动工程学报, 2014, 28(6):811-818 Wei Peng, Shen Yongjun, Yang Shaopu, Sub-harmonic resonance of Duffing oscillator with fractional-order derivative. Journal of Vibration Engineering, 2014, 28(6):811-818 (in Chinese)
[31] Atanackovic TM, Konjik S, Pilipovic S. Variational problems with fractional derivatives:Euler-Lagrange equations. Journal of Physics A Mathematical & Theoretical, 2011, 41(9):1937-1940 http://www.academia.edu/13044261/Variational_problems_with_fractional_derivatives_Euler_Lagrange_equations
[32] 王学彬.拉普拉斯变换方法解分数阶微分方程.西南师范大学学报(自然科学版), 2016, 41(7):7-12 http://www.cnki.com.cn/Article/CJFDTOTAL-XNZK201607002.htm Wang Xuebin. On laplace transform method for solving fractional differential equations. Journal of Southwest China Normal University(Natural Science Edition), 2016, 41(7):7-12 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-XNZK201607002.htm
[33] 陈明杰.分数阶傅里叶变换的数值实现.重庆大学学报(自然科学版), 2003, 26(5):129-132 http://www.cnki.com.cn/Article/CJFDTOTAL-FIVE200305031.htm Chen Mingjie. A numerical algorithms of fractional fourier transform. Journal of Chongqing University (Natural Science Edition), 2003, 26(5):129-132 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-FIVE200305031.htm
[34] 奈弗AH, 穆克DT.非线性振动.北京:高等教育出版社, 1980 Nayfeh AH, Mook DT. Nonlinear oscillations. Beijing:High Education Press, 1980 (in Chinese)
[35] 胡海岩. 应用非线性动力学. 北京: 航空工业出版社, 2000 Hu Haiyan. Applied Nonlinear Dynamics. Beijing:Aviation Industry Press, 2000 (in Chinese)
[36] 褚亦清, 李翠英. 非线性振动分析. 北京: 北京理工大学出版社,1996 Chu Yiqing, Li Cuiying. Nonlinear Vibration Analysis. Beijing: Beijing Institute of Technology Press, 1996 (in Chinese)
-
期刊类型引用(12)
1. 王翠艳，石少轩，王明昊，陈恩利，曲兴业. 含黏弹性阻尼的架桥机制动系统非线性动力学特性分析. 江苏大学学报(自然科学版). 2025(02): 221-230 . 
百度学术
2. 唐建花，李向红，王敏，申永军，李壮壮. 广义分数阶van der Pol-Duffing振子的动力学响应与隔振效果研究. 振动与冲击. 2022(01): 10-18 . 百度学术
3. 陈恩利，王明昊，王美琪，常宇健. 含有分数阶阻尼的SD振子系统非线性动力学响应特征研究. 振动工程学报. 2022(05): 1068-1075 . 百度学术
4. 秦浩，温少芳，申永军，邢海军，王军. 基于Melnikov方法的分数阶Duffing振子混沌阈值解析研究. 振动与冲击. 2021(06): 33-40+134 . 百度学术
5. 薛程，夏兆旺，卢志伟，鞠福瑜，茅凯杰. 分数阶半主动颗粒阻尼隔振系统动力学特性分析. 振动与冲击. 2021(21): 194-200 . 百度学术
6. 郭建斌，申永军，李航. 分数阶拟周期Mathieu方程的动力学分析. 力学学报. 2021(12): 3366-3375 . 本站查看
7. 李航，申永军，李向红，韩彦军，彭孟菲. Duffing系统的主-亚谐联合共振. 力学学报. 2020(02): 514-521 . 本站查看
8. 武薇，申永军，杨绍普. 基于排列熵理论的非线性系统特征提取研究. 振动与冲击. 2020(07): 67-73 . 百度学术
9. 薛程，卢志伟，苏战发，许祥曦，夏兆旺. 分数阶双层隔振系统动力学特性分析. 船舶力学. 2020(11): 1487-1494 . 百度学术
10. 赵珧冰，林恒辉，黄超辉，陈林聪. 温度场中悬索受多频激励组合联合共振响应研究. 振动与冲击. 2019(03): 207-213 . 百度学术
11. 姜源，申永军，温少芳. 分数阶van der Pol振子的超谐与亚谐联合共振. 振动工程学报. 2019(05): 863-873 . 百度学术
12. 黄超辉，赵珧冰. 多频激励下悬索非线性共振特性受温度影响分析. 应用力学学报. 2019(06): 1435-1441+1526 . 百度学术
其他类型引用(22)
下载:
计量
- 文章访问数: 1637
- HTML全文浏览量: 279
- PDF下载量: 905
- 被引次数: 34
