分数阶平方阻尼Mathieu振子的动力学分析

郭建斌, 申永军,2

(1.石家庄铁道大学 机械工程学院，河北 石家庄 050043；2.石家庄铁道大学 交通工程结构力学行为与系统安全国家重点实验室，河北 石家庄 050043)

分数阶微积分被提出以来，众多学者对其开展了详细的探讨，尤其在分数阶微积分定义和计算方法方面取得了显著成果[1-2]。在这个过程中，分数阶微积分也由数学理论研究逐步走向工程应用[3-4]。在动力学领域，利用分数阶微积分描述黏弹性材料的本构关系，提高此类非线性系统振动特性研究的准确性是当下的研究热点之一。如Cao et al[5]通过建立分数阶阻尼模型，研究了分数阶阻尼对系统的影响，展现了分数阶系统特有的动力学特性。申永军等[6-7]研究了分数阶微分项在线性和非线性系统中的作用机理，首次提出了等效线性阻尼和等效线性刚度概念。Xu et al[8]结合摄动法和多尺度法，提出了一种处理随机谐波激励作用下强非线性分数阶系统的新方法。

Mathieu方程作为Hill方程的一种特殊形式，因其复杂的动力学特性，在参激振动研究中得到了广泛应用。例如，温少芳[9]以高速列车弓网系统为对象，建立了含分数阶导数的Mathieu模型，证明了利用分数阶微分项参数表示弹簧参数的可行性。陈予恕等[10]研究了van der Pol-Duffing-Mathieu振子主参数共振的二次近似分岔行为，证明了采用摄动法描述此类系统周期响应和分岔行为的可靠性。此外，在船舶工程领域，非线性阻尼Mathieu振子是用来研究船舶横摇运动的重要模型之一。丁勇等[11]在线性加平方阻尼的基础上，建立了船舶横摇参激振动模型，研究了非线性阻尼对系统主参数共振稳态解的影响。唐友刚等[12]通过分析立方阻尼Mathieu方程，分析了系统的主参数共振，为研究船舶倾覆机理奠定了基础。

综上所述，众多学者已对典型的Mathieu系统作了深入分析，关于分数阶系统振动特性的研究也趋于成熟。在参激系统中引入分数阶微分项，研究其在此类系统中的作用规律，不仅可以完善此类系统模型，同时还可以丰富黏弹性器件的应用场景。因此，建立了含分数阶微分项的平方阻尼Mathieu模型，利用多尺度法研究系统的主共振响应，通过数值仿真分析分数阶微分项对该系统幅频特性的作用效果，所得结果为黏弹性材料在此类系统隔振、减振研究方面的应用提供了理论验证。

1 一次近似解

研究如下含平方阻尼的分数阶Mathieu振子模型

(1)

(2)

式中,Γ(z)为Gamma函数，具有Γ(z+1)=zΓ(z)的特性。

研究强迫激励频率ω≈ω0时的主共振情况，且要求激励幅值F为小量，为方便计算引入ω=ω0+εσ,F=εf,K=εk,σ=O(1) ，f=O(1),k=O(1) 。

式(1)变换为

(3)

式中，ε为小参数，满足0<ε≪1。

采用多尺度法研究系统一次近似解，引入2个时间尺度T0=t、T1=εt，并假设式(3)的解有以下形式

u(t;ε)=u0(T0,T1)+εu1(T0,T1)

(4)

将式(4)代入式(3)，比较ε的同次幂，得到一组偏微分方程

(5)

(6)

式(5)的解为

u0(T0,T1)=a(T1)cos[ω0T0+β(T1)]

(7)

式中，a(T1)、β(T1)为慢变振幅、相位。为方便计算，可将式(7)写成复数形式

u0(T0,T1)=A(T1)eiω0T0+cc

(8)

(9)

将式(8)和式(9)代入式(5)，可得到消除永年项条件

(10)

(11)

(12)

其中

θ=ωT0+β

(13)

再引入

φ=β-σT1

(14)

(15)

(16)

u(t)=acos(ω0T0+β)=acos(ω0T0+σT1+φ)=acos(ωt+φ)

(17)

式中，a和φ由式(15)、式(16)确定。

由式(17)可见，系统主共振近似解的振动频率等于强迫激励频率且是参数激励频率的1/2，此外，相比传统的Mathieu系统，分数阶微分项分量K(p)和C(p)分别作用于周期解的幅值和相位，使其存在定量上的差别，导致了响应幅值的降低和相位的滞后。

为验证近似解的准确性，利用式(1)进行数值仿真，得出的数值解与式(17)计算出的近似解析解作对比。利用文献[2]中介绍的数值方法研究系统(1)，该方法的近似公式为

(18)

(19)

此外，当t=0时，分数阶微分项的初值为

(20)

图1 幅频曲线对比 图2 位移时间历程图对比

2 定常解及其稳定性

(21)

(22)

可由式(21) 、式(22)推导出非零定常解的幅频和相频方程

(23)

(24)

此外，以上求得的非零定常解能否实现还取决于其是否具有渐近稳定性。这里利用Lyapunov第一方法来计算稳态运动的稳定性条件，以此来考察解的稳定性。用慢变振幅a和相位φ定义二维状态向量V=[a,φ]T，构造向量函数

(25)

(26)

λ2-Pλ+Q=0

(27)

式中，P=trJ;Q=det[J]。

(28)

-[C(p)+μ1]-4aμ2ω0<0

(29)

图3 定常解幅频曲线

3 分数阶微分项对系统的影响

3.1 分数阶微分项阶次的影响

首先考察分数阶微分项阶次对系统动力学特性的影响。取一组基本参数ε=0.1,μ1=0.3,μ2=0.3,ω0=1,K=0,F=0.01对系统进行仿真计算。分别选取不同阶次p，利用式(18)计算系统的幅频响应如图4所示。从图4可以看出，随着分数阶阶次p从0.1增加至0.9，系统的稳态响应幅值在逐渐变小，且幅频曲线整体逐渐向低频方向偏移。此外，阶次p还改变了幅频曲线的拓扑结构，系统的多值现象逐渐消失。图5给出F=0.1时阶次p对稳态响应共振峰值的影响，可见此时分数阶微分项对共振幅值有着明显的抑制作用。

图4 分数阶微分阶次p对系统的影响( K=0.05) 图5 阶次p对系统共振峰值的影响(F=0.1)

3.2 分数阶微分项系数的影响

分别考察p→0和p→1 情况下系数K对系统幅频特性的影响，图6～图8中圆圈为稳定解，星号为不稳定解。首先考虑p→1的情况，固定系统参数，依次取系数K为0.01、0.02、0.03和0.05，系统幅频曲线随K的变化情况如图6所示。为方便比较，图7中给出了系统幅频曲线随线性阻尼系数μ1的变化情况。通过对比图6和图7发现，当分数阶阶次p→1时，随着K的逐渐增大，稳态响应的幅值在逐渐缩小，并且K达到一定值时，导致系统幅频曲线形态发生变化，由参数激励和强迫激励共同作用引起的多解现象消失，改变了定常解的稳定性。同样地，随着线性阻尼系数μ1逐渐增大，系统的幅频响应曲线发生了类似变化。以上情况说明，阶次p→1时分数阶微分项系数K对系统的作用几乎等同于线性阻尼系数μ1。

图6 分数阶系数K对幅频曲线的影响(p→1)

图7 线性阻尼系数μ1对幅频曲线的影响(p=0.6,K=0.01)

下面考虑p→0的情况。其他参数不变，取p=0.1，通过改变系数K来观察分数阶微分项对系统幅频特性的影响如图8所示。可以看出，随着分数阶系数K的逐渐增大，幅频曲线逐渐向高频方向偏移，改变了系统的共振频率，但系统的响应幅值并未受到明显影响，说明此时分数阶微分项呈现较强的刚度特性。

图8 分数阶系数对幅频响应的影响(p→0)

4 结论

应用多尺度法研究了强迫激励下分数阶平方阻尼Mathieu振子的主共振，建立了定常解的幅频响应方程。利用Lyapunov理论分析了系统的幅频特性，由于参数激励和强迫激励的共同作用，在共振区域内稳态响应系统至多存在3个解支。此外，通过数值仿真分析了分数阶微分项对系统幅频曲线的影响，发现改变分数阶微分项的阶次或系数可使其对系统幅频特性产生不同程度的影响：p→1时，分数阶微分项呈现出较强的阻尼特性，其对系统的作用几乎等同于线性阻尼，改变系数K主要影响系统响应幅值；p→0时，分数阶微分项呈现出较强的刚度特性，改变系数K主要影响系统的共振频率。以上结果揭示了分数阶微分项(黏弹性器件)在此类系统中的作用规律，验证了其对系统响应特性的作用效果。