分数阶单自由度线性振子双侧对称碰撞振动分析

刘汝逾，牛江川,2，申永军,2，杨绍普,2

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

碰撞振动系统是一种非光滑系统，存在许多光滑系统不具备的复杂现象。由于碰撞振动的复杂性，吸引了许多学者对其进行研究。碰撞振动广泛存在于机械设备中，例如刨煤机刨头接触煤炭时会发生碰撞[1]和轮船航行中会与楔形物发生碰撞[2]。吕小红等[3]通过数值仿真研究了一类具有双侧刚性约束的单自由度碰撞振动系统，分析了系统由周期运动向混沌转迁的分岔过程。张继业等[4]研究了单自由度自治系统的碰撞振动，得到了碰撞振动周期运动存在的条件及周期运动稳定的充分性判据。罗冠炜等[5]研究了一类单自由度碰撞振动系统的动力学行为，描述了系统周期运动的规律，分析了映射奇异性对系统局部分岔和全局分岔的影响。冯进钤等[6-7]以典型的Duffing单边碰撞系统为例，研究了系统的颤振分岔，并对系统中的混沌鞍进行了细致的分析。李得洋等[8]利用数值仿真方法分析了单自由度含对称弹性约束碰撞的振动系统，得到了在两参数平面内周期运动的分布及转迁规律。乐源等[9]讨论了一类单自由度双面碰撞振子的对称型周期n-2运动以及非对称型周期n-2运动，并通过分析对称系统的Poincaré映射的对称性证明了对称型周期运动只能发生音叉分岔。Andreaus等[10]建立了含有双侧间隙和保险杠的单自由度双侧碰撞模型，分析了谐波激励强度、间隙宽度、保险杠刚度和伪共振频率等最重要的参数对单自由度系统响应的影响。Kumar等[11]研究了具有刚性势垒的Duffing-van der Pol振荡器在白噪声激励下的分岔问题，利用非光滑变量变换将非光滑碰撞振动系统映射到连续相平面上，并考虑了两侧刚性约束不对称布置的情况。Reboucas等[12]使用逐点映射和标准平均与非平滑变换相结合，分析了具有恢复系数的单自由度模型的振动冲击响应。Xie等[13]通过随机平均和Mellin变换的方法研究了非线性振动系统在高斯白噪声激励下的近似瞬态响应。Li[14]针对冲击振动系统的不连续性，引入非光滑变换克服不连续性，然后将随机平均技术应用于变换后的随机变质量系统，通过求解Fokker-Planck-Kolmogorov方程，可以得到近似的解析响应，以Van der Pol变质量振动冲击系统为例，说明了该方法的有效性。Rigaud等[15]提出了由齿轮啮合的振动冲击引起齿轮嘎嘎声并进行了试验研究，设计了一个特定的试验装置分析正齿轮在输入速度波动下的非线性动力学行为。Reboucas等[16]采用谐波线性化、平均化和数值模拟的方法，对单自由度振动冲击振荡器的频率响应进行了分析，并考虑了3种不同的冲击力模型。Chen等[17]针对高斯白噪声激励下单自由度振动冲击系统响应的封闭形式平稳概率密度函数，提出了一种新的求解方法,采用Zhuravlev-Ivanov变换将单侧势垒的振动冲击振荡器转换为无势垒的振荡器。Nguyen等[18]采用数学建模和动力学分析的方法，研究了惯性质量和激励频率对自成移动系统的Duffing冲击振动装置的耦合影响。Yan等[19]研究了单侧和双侧软约束振动冲击胶囊系统在质量比、刚度比、接触间隙、外激励幅值和频率等控制参数变化下的动力学特性，得到单侧约束胶囊的动力学行为主要是周期性的，当左侧约束刚度增大时，具有两侧约束的其他两个胶囊的动力学响应变得复杂的结论。Xu等[20]研究了具有单边非零偏移屏障的非线性振动系统的多值响应和动力学特性，发现在某些情况下，冲击系统可能具有两个或4个稳态解。然后讨论了稳态解的分类，结果表明非平凡的稳态解可能会因鞍结分岔和Hopf分岔而失去稳定性。Xu等[21]给出了一般非线性振动冲击振子的Melnikov方法，以Duffing振动冲击振荡器为例说明了该方法的应用，并通过相图、Poincaré截面和分岔图对所得结果进行了验证。

由于分数阶微积分可以提供一个强有力的工具对工程中的黏弹性等行为进行建模，近年来关于分数阶微积分的研究受到越来越多的关注。廖少锴等[22]研究了非线性分数微积分求解的单步数值积分算法，讨论了振子自由振动及强迫振动下参数变化对振子非线性特性的影响。杨建华等[23]利用谐波平衡法得到了含分数阶导数阻尼的一类线性系统在简谐激励下系统响应的近似解。黄灿等[24]利用分数阶导数算子研究了线性分数阶振动系统在谐波激励下的稳态响应，并讨论了分数阶导数项对刚度和阻尼的影响。银花等[25]基于精细积分方法，提出了一种新的具有分数阶导数型本构关系的黏弹性结构动力响应的数值计算方法。王学彬[26]给出了两种常见分数阶导数即Riemann-Liouille分数阶导数和Caputo分数阶导数的拉普拉斯变换公式，并给出了具体实例。闫启方等[27]运用黏弹性理论、分数阶导数的性质，研究了分数阶Kelvin模型描述的黏弹性材料的松弛特性、蠕变特性和动态力学行为。刘艳芹[28]将分数阶微分算子引入到黏弹性介质中的阻尼振动中建立了分数阶非线性振动方程，并利用Adomian分解方法借助Mathematica软件求解了该类分数阶阻尼振动方程的近似解，研究了振子运动与方程中分数阶导数的关系。牛江川等[29]研究了基于速度反馈分数阶PID(proportional integral derivative)控制的单自由度线性振子的自由振动，利用平均法得到了系统的近似解析解。Yang等[30]研究了高斯白噪声激励下两类分数阶导数冲击振荡器的模型，利用非光滑变换和随机平均方法得到解析近似解，并对其随机分岔进行了讨论。

关于分数阶黏弹性材料的应用和碰撞振动系统的研究成果已经非常丰富，但是关于附加黏弹性材料碰撞振动系统的研究却很少。将分数阶微积分引入到碰撞振动系统的研究，可以更深入地理解附加黏弹性材料的碰撞振动行为。本文研究了系统关键参数对系统周期运动的影响，可以为附加黏弹性材料的传动齿轮、铸造振动筛和压力成型机等机械设备存在的非线性问题提供借鉴，避免这些系统在混沌状态下工作，提高机械设备的动态性能，实现对复杂设备的振动和噪声的有效控制，为碰撞振动系统的优化设计和改造提供一定的依据。

本文利用分数阶微分描述黏弹性材料的黏弹性，然后基于分数阶单自由度线性双侧刚性碰撞模型，研究双侧对称碰撞振动系统在简谐激励下的稳定性和分岔行为。对于碰撞振动的稳态解，利用平均法得到分数阶线性系统的等效刚度和等效阻尼；对于碰撞振动的瞬态解，利用迭代法得到更精确的瞬态固有频率。在此基础上，分析了双侧对称碰撞振动系统的近似解析解。根据近似解析解，分析了对称n-1-1周期运动的存在条件，并利用Poincaré映射研究了n-1-1周期运动的稳定性。进一步分析了当外界激励频率、分数阶阶次和间隙变化时该系统的分岔行为。

1 近似解析解

分数阶单自由度双侧碰撞模型，如图1所示。其中包含刚度系数为k的线性弹簧、阻尼系数为c的线性阻尼和系数为K1(K1>0)的黏弹性材料。振子m受到周期性的激励力Fcos(ωt+θ)作用，其中振幅为F、激励频率为ω、初始相位为θ。

图1 含间隙的单自由度双侧碰撞物理模型Fig.1 Single-degree-of-freedom bilateral vibro-impact physical model with gaps

当没有碰撞时的数学模型可以表示为

(1)

式中，Dp[x(t)]为x(t)关于t的p阶导数(0≤p≤1)。这里我们采用Caputo[31]定义

(2)

式中，Γ(z)为Gamma函数，满足Γ(z+1)=zΓ(z)。

当碰撞发生时满足以下关系

(3)

由于存在阻尼和分数阶导数，因此在谐波激励下系统解需要分别考虑稳态解和瞬态解。首先假设系统运动无碰撞，然后研究其稳态解。应用平均法得到两个参数Ceq(p)和Keq(p)，它们表示新的阻尼系数和新的刚度系数[32-33]。

Ceq(p)=c+K1ωp-1sin(pπ/2)

(4)

Keq(p)=k+K1ωpcos(pπ/2)

(5)

然后得到式(1)稳态解的一阶近似等价整数阶微分方程

(6)

此外，研究式(1)的瞬态解则无需考虑外部激励，将运动微分方程重写为

(7)

类似地，系统自由振动的新阻尼系数和新刚度系数可以表示为

(8)

(9)

从而得到式(1)瞬态解的一阶近似等价整数阶方程

(10)

(11)

(12)

(13)

(14)

(15)

(16)

系统的速度可以由式(15)和式(16)求导得到。

2 数值解验证

图2 近似解析解和数值解拟合Fig.2 Approximate analytical solution and numerical solution fitting

在不考虑碰撞的情况下，将数值解与解析解进行比较。结果表明，近似等效整数阶系统的稳态解具有较高的精度，而瞬态解的精度较低。另外，分数阶系数K1越大，瞬态解的拟合误差越大。从碰撞振动系统的特征可以看出，碰撞振动系统的每个部分都是瞬态运动。但是，当K1相对较小时，误差要小得多。考虑到振动的影响，我们只需改变上述参数K1即可分析系统的拟合误差。当K1=0.5，K1=1.0和K1=1.5时，系统数值解和近似解析解的局部放大拟合图，如图3～图5所示。随着K1的增大，瞬态解的拟合误差将会增大。

图3 K1=0.5的瞬态解的拟合误差Fig.3 Fitting error with K1=0.5

图4 K1=1.0的瞬态解的拟合误差Fig.4 Fitting error with K1=1.0

图5 K1=1.5的瞬态解的拟合误差Fig.5 Fitting error with K1=1.5

3 n-1-1周期运动的存在性分析

在合适的参数下，系统可以表现为n-p-q周期运动。其中：n为外界激励周期数；p和q分别为振子与刚性约束W1和W2碰撞的次数。撞击面被视为Poincaré截面，并表示为

(17)

假设振子与刚性约束W1之间的第一次碰撞发生在时刻t=t1，则振子移动到刚性约束W2处并在时刻t=t2与W2发生碰撞，具体的流程如图6所示。根据图6所示的周期性n-1-1运动的特性，系统必须满足

图6 对称周期运动相图Fig.6 Phase portrait of symmetric periodic motion

(18)

图的对称周期运动相图Fig.7 Phase portrait of symmetric periodic motion with

4 n-1-1周期运动的Poincaré映射及其稳定性

本章利用Poincaré映射研究n-1-1周期运动的稳定性。根据式(1)，等效平面系统可以表示为

(19)

把式(19)的二维系统改成三维系统，可以得到

(20)

同时将Poincaré横截面定义为

(21)

根据双侧对称碰撞的特征，碰撞振子的运动轨迹在一个周期内分为4个阶段：①碰撞振子离开刚性约束W1并移动到刚性约束W2；②碰撞振子与刚性约束W2发生碰撞；③碰撞振子离开刚性约束W2并移动到刚性约束W1；④碰撞振子与刚性约束W1发生碰撞。

因此，可以将这4个映射表示为

(22)

根据式(13)和式(14)并令τ(0)=θ-φ，可以得到Poincaré映射为

(23)

假设DP1，DP2，DP3和DP4分别对应4个Poincaré映射的连续阶段的线性化矩阵，则可以有

(24)

当碰撞振子与刚性约束碰撞时，只改变速度而不改变位移，且速度恢复系数为r，可以得到

(25)

则与P1光滑矢量场相关的局部流可表示为

(26)

(27)

映射关系可以建立为

(28)

(29)

为了进行以下计算，可以将隐式函数式(29)表示为

(30)

(31)

因此，可以得到DP1为

DP1=

(32)

DP3=

(33)

然后可以得到碰撞振动系统的线性化矩阵为

DP=DP4·DP3·DP2·DP1

(34)

应用DP的特征值可以判断碰撞振动系统周期运动的稳定性。当两个特征值都位于单位圆内时，系统的周期运动是稳定的。只要一个特征值不处于单位元内，周期运动在不动点就不稳定。根据第3章中的参数，可得两个特征值为λ1=0.313 020和λ2=0.424 067，均位于单位元内。可以得出结论：其n-1-1周期运动在不动点处是稳定的，也可以从图7中得到验证。

5 分岔分析与讨论

激励频率ω对分数阶碰撞振动系统的动力学行为有很大的影响。首先取第2章的系统参数，仅更改K1=1.5。首先选取激励频率ω作为分岔变量，随着激励频率的变化，根据数值解可以得到系统在碰撞面上的分岔图，如图8所示。可以看出，当激励频率ω变化时，系统发生不连续擦边分岔。当ω为3.000～1.000时，系统表现为1-1-1周期运动、1-2-2周期运动、1-3-3周期运动、1-4-4周期运动，然后是黏滑振动。

图8 外部激励频率ω的分岔图Fig.8 Bifurcation diagram with ω

当ω=2.682和ω=2.683时的相图，如图9所示。结果表明，系统从1-1-1周期运动经过不连续擦边分岔为1-2-2周期运动，而且在分岔前后都是稳定的周期运动。ω=2.682和ω=2.683的位移时间历程，如图10和图11所示。当ω=1.780和ω=1.420时同样也发生不连续擦边分岔。

图9 不连续擦边分岔的相图Fig.9 Phase portrait of discontinuous boundary bifurcation

图10 ω=2.682的位移时间历程Fig.10 Displacement time history of ω=2.682

图11 ω=2.683的位移时间历程Fig.11 Displacement time history of ω=2.683

图12 分数阶阶次p的分岔图Fig.12 Bifurcation diagram of fractional order p

当p介于0.390～1.000时，系统表现为对称的1-1-1周期运动。例如：当p=0.800时，系统表现为对称的1-1-1周期运动，其相图如图14所示。当p=0.300时，系统表现为非对称的1-1-1周期运动，其相图如图13所示。当系统参数改变时，系统会发生擦边分岔。在上述参数给定的情况下，分岔的分数阶临界值为p=0.276，如图15所示。

图13 p=0.800的对称周期运动相图Fig.13 Phase portrait of symmetric periodic motion with p=0.800

图14 p=0.300的非对称周期运动相图Fig.14 Phase portrait of asymmetric periodic motion with p=0.300

图15 p=0.276的擦边分岔相图Fig.15 Phase portrait of grazing bifurcation with p=0.276

当p<0.276时，系统从非对称周期1-1-1运动转变为拟周期运动。随着p减小，系统表现为混沌运动。当p=0.100时，可以得到系统的混沌吸引子，如图16所示。

图16 p=0.100的混沌运动Fig.16 Chaotic motion with p=0.100

图17 关于碰撞间隙Δ的分岔图Fig.17 Bifurcation diagram of collision gap Δ

图18 Δ=0.850的对称周期运动相图Fig.18 Phase portrait of symmetric periodic motion with Δ=0.850

图19 θ=0,Δ=0.850的非对称周期运动相图Fig.19 Phase portrait of asymmetric periodic motion with θ=0,Δ=0.850

图20 θ=π,Δ=0.850的非对称周期运动相图Fig.20 Phase portrait of asymmetric periodic motion with θ=π,Δ=0.850

图21 Δ=0.630的拟周期运动相图Fig.21 Phase portrait of quasiperiodic motion with Δ=0.630

6 结 论

本文研究了具有两个刚性约束的分数阶线性碰撞振动系统在正弦激励下的稳定性和分岔行为。通过引入分数阶系统的近似等价整数阶系统，得到了系统的近似解析解。基于近似解析解得到了双侧对称n-1-1周期运动的存在条件，并用Poincaré映射分析了n-1-1周期运动的稳定性。通过改变外激励频率、分数阶阶次和间隙，利用数值解对双侧对称碰撞振动系统的分岔行为进行了深入的分析。结果表明，当激励频率变化时，系统会出现不连续擦边分岔，且系统在分岔前后始终是稳定的周期运动。此外，在另一组参数下，随着分数阶阶次的减小，出现了音叉分岔和擦边分岔。另外，当以间隙为控制参数时，随着间隙的减小，对称不动点变得不稳定，而且会产生一对反对称倍周期序列，将导致拟周期运动和混沌。该模型的计算过程和结论可为分数阶碰撞振动系统的动力学行为研究提供参考。