李航 申永军† 杨绍普 彭孟菲 韩彦军

1) (石家庄铁道大学, 交通工程结构力学行为与系统安全国家重点实验室, 石家庄 050043)

2) (石家庄铁道大学机械工程学院, 石家庄 050043)

非线性动力系统极易发生共振, 在多频激励下可能发生联合共振或组合共振, 目前关于非线性系统的主-超谐联合共振的研究少见报道.本文以Duffing 系统为对象, 研究系统在主-超谐联合共振时的周期运动和通往混沌的道路.应用多尺度法得到系统的近似解析解, 并利用数值方法对解析解进行验证, 结果吻合良好.基于Lyapunov 第一方法得到稳态周期解的稳定性条件, 并分析了非线性刚度对稳态周期解的幅值和稳定性的影响.此外, 由于近似解只能描述周期运动, 不足以描述系统的全局特性, 因而应用Melnikov 方法对系统进行全局分析, 得到系统进入Smale 马蹄意义下混沌的条件, 依据该条件以及主-超谐联合共振的条件选取一组参数进行数值仿真.分岔图和最大Lyapunov 指数显示出两个临界值: 当激励幅值通过第一个临界值时, 异宿轨道破裂, 混沌吸引子突然出现, 系统以激变方式进入混沌; 激励幅值通过第二个临界值时, 系统在混沌态下再次发生激变, 进入另一种混沌态.利用Melnikov 方法考察了第一个临界值在多种频率组合下的变化趋势, 并用数值仿真验证了解析结果的正确性.

1 引 言

Duffing 系统是经典的非线性动力系统, 能够描述转子系统[1−2]、船舶横摇[3−4]和受电弓弓头悬挂系统的垂向振动[5−6]等多种动力学问题, 在微弱信号检测[7−9]和能量采集[10−11]领域也有应用.目前在对Duffing 振子的研究中, 系统的周期振动解及其稳定性受到学者关注.Shen 等[12,13]研究了分数阶Duffing 系统的主共振, 利用平均法得到系统的近似解, 分析了分数阶项对系统稳态响应的影响, 提出等效线性刚度和等效线性阻尼的概念.Agarwal 等[14]研究了高斯白噪声激励下刚度软化Duffing 系统的频响特性, 并通过实验分析了噪声幅值对系统跳跃现象的影响, 发现噪声可以破坏无噪声系统中的迟滞特性.温少芳等[15]研究了一类含分数阶时滞耦合反馈的Duffing 系统, 利用平均法得到系统的一阶近似解, 并分析了反馈系数和分数阶项对系统动力学特性的影响, 发现分数阶时滞耦合反馈同时具有速度时滞反馈和位移时滞反馈的作用.

通往混沌的道路与混沌控制也是Duffing 系统动力学行为研究的重点问题之一.杨晓丽等[16]研究了有界噪声与周期激励下Duffing 系统的混沌运动, 基于Melnikov 方法得到了系统进入混沌运动的必要条件, 讨论了噪声对系统混沌阈值的影响.马少娟等[17]和吴存利等[18]研究了含随机参数的Duffing 系统, 利用多项式逼近方法将随机系统等效为确定性系统, 用确定性系统的方法研究了随机系统的混沌运动.Niu 等[19]研究了一类分数阶Duffing 系统的混沌运动, 利用谐波平衡法得到了原系统的等效整数阶系统, 基于Melnikov 方法得到系统进入混沌的条件.在对Duffing 系统周期振动解和混沌运动的研究中也发展了一些新的研究方法.Liu 等[20]提出推广的广义胞映射法可用于分数阶Duffing 系统的全局分析.Shen 等[21]和Niu等[22]对增量谐波平衡法的推广可用于分数阶系统和强非线性系统的高阶近似分析.

非线性系统极易发生共振, 除主共振外还具有一些特殊的共振现象, 如超谐共振[23]、亚谐共振[24]等, 受多频激励的非线性系统还可能发生组合共振[25]或者联合共振[26,27].尤其在联合共振下,Duffing 系统的稳态周期解可多达7 个分支, 其中包括4 个稳定分支, 系统的响应在稳定分支间的跳跃现象更为复杂.此外, 在多个激励频率不可公约时, Duffing 系统将产生概周期响应[28], 并且随着系统参数改变, 概周期响应容易失稳进入混沌[29].毕勤胜等[30]研究了一类受两个简谐激励且频率不可公约的Duffing 系统, 对系统在非共振情况下通往混沌的道路进行研究, 发现随着激励幅值改变,概周期响应的某一周期分量或某一组合频率成分将发生倍化分岔导致混沌.在实际问题中, 动力系统往往受到多个激励源的作用, 因此对多频激励的非线性系统进行研究是有必要的.鉴于针对非线性系统在主-超谐联合共振情况下的研究极为少见,本文以多频激励Duffing 系统为对象, 考察系统在主-超谐联合共振时的周期运动与混沌运动.通过多尺度法得到系统的近似解析解, 研究周期运动在Lyapunov 意义下的稳定性.基于Melnikov 方法分析系统的全局动力学特性, 研究Duffing 系统在主-超谐联合共振时通往混沌的道路.

2 Duffing 系统的主-超谐联合共振的一次近似解

研究受两个简谐激励的Duffing 系统:

其中 m , c , k , α1, Fi和 ωi( i =1,2 )分 别为系统的质量、线性阻尼系数、线性刚度系数、非线性刚度系数、激励幅值和激励频率.

对系统进行如下坐标变换:

其中 ε 为小参数, 满足 0 <ε ≪1 , ω0为系统的固有频率.

于是(1)式可以等价地写成

为研究系统的主-超谐联合共振, 要求 ω1≈ω0,且 f1为小量, 即

其中 σ1和 σ2为引入的调谐因子, 从而(2)式成为

应用多尺度法[27]研究系统的一次近似解.引入两个时间尺度 T0=t , T1=εt , 并假设(3)式的解具有如下形式:

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

(5a)式的解为

其中 a (T1) 和 β (T1) 分别为慢变振幅和相位.

(6)式可写成复数形式:

将(7)式代入(5b)式得到:

由此得到消除永年项的条件为

从而系统(1)在主-超谐联合共振下的一次近似解可以表示为

其中a 和 β 是 T1的函数, 由(10)式确定.

3 周期运动的稳定性

3.1 Lyapunov 意义下的稳定性条件

相比瞬态响应, 稳态响应的稳定性更值得关注.观察(10)式可以发现, 定常解存在, 即D1a=0和 D1β =0 有解的必要条件是 β −σ1T1和β −3σ2T1是常数, 这时有 D1β =σ1=3σ2.进一步可以得到ω1=3ω2, 即只有当两个激励频率满足此倍数关系时, 才能通过解析方法求得主-超谐联合共振的定常解.

引入调谐因子 σ =σ1=3σ2, 再令 β −σT1=φ ,(10)式可以写成自治的形式:

相应的一次近似解变成:

其中 a 和 φ 是 T1的函数, 由(12)式确定.

为检验一次近似解(13)式的正确性, 利用变步长四阶五级Runge-Kutta 法对系统(1)进行数值仿真并和一次近似解对比, 取系统参数m=1, c=0.02, k =4, α1=0.4, F1=0.1, F2=2,ω2=ω1/3.首先对比稳态响应的幅频曲线, 仿真时长 t =1000 , 取后 2 0% 响应的最大值为稳态幅值,得到解析解与数值解的幅频曲线对比结果如图1所示.然后对比位移时间历程, 取激励频率 ω1=2.2 ,ω2=ω1/3 , 解析解初值取 ( a0,φ0)=(0.1,0) , 对应的数值解初值由(13)式计算, 仿真得到解析解与数值解的位移时间历程对比结果如图2 所示.

图1 幅频曲线对比Fig.1.Comparisons of amplitude-frequency curves.

从图1 和图2 中可见, 解析解和数值解在较宽的频率范围内吻合良好, 即使是瞬态响应也能够较为精确地描述.此外, 从图1 中可以看到系统的稳态响应在一定频率范围内具有多值性, 事实上, 系统的多值性和稳定性都由近似解中的第一部分, 即a cos(ω1t+φ)支配.因此, 在系统周期运动的稳定性分析中只需考察这一部分.

图2 位移时间历程对比 (a) 瞬态响应; (b) 稳态响应Fig.2.Comparisons of displacement time histories: (a) Transient response; (b) steady-state response.

令(12)式中的 D1a=0 , D1φ=0 , 可以得到

定常解的振幅 a¯ 和相位 φ¯ 满足的方程组:

从(14)式导出定常解的幅频响应方程和相频响应方程分别为

下面考察定常解的稳定性, 定义状态向量V =[a,φ]T, 构造向量函数:

其中 P =trJ , Q =det[J].

根据Lyapunov 第一方法[31], 系统在定常解处 渐近稳定的条件是 P <0 且 Q >0.对于阻尼系统, 总有 P <0 , 因此系统的稳态周期运动在Lyapunov 意义下的稳定性条件为

3.2 数值仿真

仍取参数 m =1 , c =0.02 , k =4 , F1=0.1 ,F2=2 , ω2=ω1/3.由(15)式计算定常解的幅频特性和相频特性, 再由(18)式判断对应解支的稳定性, 得到 α1=0.4 和 α1=−0.4 时定常解的幅频曲线和相频曲线分别如图3 和图4 所示.从图3 和图4 可见, 系统最多存在2 个稳定解支和1 个不稳定解支, 且图3(b)的多值区域与图1 吻合.如前所述, 许多工程振动问题都可以由Duffing 系统描述.转子系统在启停过程中可能经过共振区, 这时系统的响应将会在多个稳定解支之间切换, 振幅突然增大或减小, 称为跳跃现象.Wang 等[11]基于Duffing系统设计了一种能量采集装置, 在能量采集实验中观察到这种跳跃现象.

确定稳态响应的多值和不稳定的频率范围是一个重要问题, 利用解析解能够在短时间内计算得到较为精确的结果.利用(15)式和(18)式仿真得到非线性刚度系数 α1对系统的影响如图5 所示.对于刚度软化系统, 随着系统刚度逐渐软化, 共振峰逐渐降低, 多解区域扩大; 对于刚度硬化系统,随着系统逐渐硬化, 共振峰也逐渐降低, 但对多解区域影响不大.

图3 定常解的幅频曲线 (a) α 1=−0.4 , 刚度软化; (b)α1=0.4, 刚度硬化 (圆表示稳定解支, 星号表示不稳定解支)Fig.3.Amplitude-frequency curves of steady-state response:(a) α 1 =−0.4 , stiffness softening; (b) α 1 =0.4 , stiffness hardening (the circles for stable solution and the asterisks for unstable one).

4 混沌运动分析

4.1 系统的混沌阈值

摄动理论得到的近似解具有一定局限性.由于系统在一定条件下可能进入混沌, 此时近似解无法描述系统的响应, 因而需对系统进行全局分析, 确定系统的混沌阈值.Melnikov 方法是检测动力系统混沌阈值的常用解析方法, 通过考察过鞍点的稳定流形和不稳定流形之间的距离, 得到系统的同宿轨道或异宿轨道发生横截相交的数学条件, 从而可以判断系统是否进入Smale 马蹄意义下的混沌.

在此考虑系统(1)刚度软化(即 α1<0 )的情况, 对(1)式进行如下坐标变换

其中 ς 为小参数.从而(1)式可以写成:

将(19)式写成状态方程形式:

图4 定常解的相频曲线 (a) α 1=−0.4 , 刚度软化; (b)α1=0.4, 刚度硬化(圆表示稳定解支, 星号表示不稳定解支)Fig.4.Phase-frequency curves of steady-state response:(a) α 1 = − 0.4 , stiffness softening; (b) α 1 =0.4 , stiffness hardening (the circles for stable solution and the asterisks for unstable one).

当 ς =0 时, 未扰系统为Hamilton 系统:

其Hamilton 量为

由于

图5 非线性系数 α 1 的影响 (a) 刚度软化; (b) 刚度硬化 (圆表示稳定解支, 星号表示不稳定解支)Fig.5.Effects of the nonlinear coefficient α 1 : (a) Stiffness softening; (b) stiffness hardening (the circles for stable solution and the asterisks for unstable one).

设 t =0 时, u1(0)=0 , 由 (25a)式 可以得到从而可以由(25)式解得 u1和 u2满足:

将(19)式写成向量函数形式:

受扰系统(19)的Melnikov 函数为

从而

因此, 系统(1)出现异宿轨道横截相交的必要条件是:

用原系统参数表示为

4.2 数值仿真

由(30)式可知, 对于刚度软化的系统(1), 激励 Fi越强、阻尼 c 越小和非线性刚度 α1越软, 易于使系统的异宿轨道发生横截相交, 从而系统(1)进入Smale 马蹄意义下的混沌.由于激励幅值 F1已限制为小量, 而 F2的选取较为自由, 因而考察 F2对系统动力学特性的影响.取系统参数 m =1 , c =0.3 ,k =1 , α1=−1 , F1=0.2 , ω1=1.1 , ω2=ω1/3 ,代入(30)式得到 F2>0.0623.即在上述的参数下,当 F2>0.0623 时系统(1)可能出现异宿轨道横截相交, 并且 F2越大于此值, 系统越可能进入混沌.但该条件仅为系统进入混沌的必要条件, 还需进一步通过数值仿真才能得到系统进入混沌的临界值.

在不同的激励幅值 F2下利用自由插值的梯形方法[32,33]对系统(1)进行数值仿真, F2的起始值为 0.0600 , 初值为 ( u01,) = ( 0.208010172501537,0.326464098793873), 仿真时长为 t =200π/ω2, 步长为 2 ×10−4π/ω2.仿真得到系统的分岔图如图6所示.奇怪吸引子的相轨迹线在局部上总是以指数速率分离的, Lyapunov 指数通过度量相邻相轨迹线分离或收敛的平均变化率判断系统是否进入混沌, 是一种可靠的定量分析方法.利用Jacobian 法[34]计算得到这组参数下系统的最大Lyapunov 指数随 F2变化的趋势如图7 所示.

图6 系统(1)随 F 2 变化的分岔图 (a) 整体图; (b) 局部放大图Fig.6.Bifurcation diagram of Eq.(1) changing with F 2 :(a) Panoramic view; (b) local enlarged view.

图7 系统(1)的最大Lyapunov 指数随 F 2 变化的趋势 (a) 整体图; (b) 局部放大图Fig.7.The trend of the largest Lyapunov exponent of Eq.(1) changing with F 2 : (a) Panoramic view; (b) Local enlarged view.

图8 F 2 =0.259 时的周期轨道Fig.8.Periodic orbits for F 2 =0.259.

从图6 可见, 当 F2>0.259 时, 系统发生激变进入混沌态, 与图7 中 F2>0.259 时最大Lyapunov指数的骤增吻合.取 F2=0.259 , 初值为 ( u02,) =( 0.110106001335058, 0.248137954024669) , 得到系统的周期轨道如图8.取 F2=0.2595 , 初值为(u03,)= ( 0.091801856843936, 0.382525085869613),仿真得到系统的位移时间历程、相轨迹和Poincare截面如图9 所示.从图8 和图9 可见, 激励幅值F2从 0.259 增加到 0.2595 时异宿轨道破裂, 周期轨道发展成为混沌吸引子, 系统的相轨迹受到左侧鞍点吸引.此外, 从图6 可见系统在 F2>0.268 时再次发生激变, 进入另一种混沌态.取 F2=0.268 , 初值为(u04,= (0.144691373245013, 0.202712823433900),仿真得到系统的位移时间历程和相轨迹如图10 所示, 这时系统的相轨迹被左右两个鞍点吸引, 出现新的混沌吸引子.

图9 F 2 =0.2595 时的混沌状态 (a) Poincare 截面; (b) 位移时间历程; (c) 相轨迹Fig.9.Chaotic state for F 2 =0.2595 : (a) Poincare map;(b) displacement time histories; (c) phase trajectory.

4.3 解析结果与数值仿真结果的关系

在前述的系统参数下, 通过Melnikov 方法解析得到激励幅值 F2的临界值为 0.0623 , 通过数值仿真得到系统首次出现混沌的临界值为 0.259.为进一步研究解析与数值仿真得到的临界值的关系, 仍取系统参数 m =1 , c =0.3 , k =1 , α1=−1 , F1=0.2 ,并且激励频率保持严格的倍数关系 ω1=3ω2, 分别利用数值仿真和Melnikov 方法解析求解不同的频率组合下 F2的临界值, 得到仿真结果与解析结果的对比如图11(a)所示.

图10 F 2 =0.268 时的混沌状态 (a) 位移时间历程; (b) 相轨迹Fig.10.Chaotic state for F 2 =0.268 : (a) Displacement time histories; (b) phase trajectory.

然后考察激励频率为近似倍数关系ω1≈3ω2的情况.其他参数不变, 并且固定 ω1=1.1 , 求解不同的激励频率 ω2下系统出现混沌的临界激励幅值F2, 得到解析与数值仿真结果对比如图11(b)所示.再固定 ω2=0.3667 , 求解不同的激励频率ω1下系统出现混沌的临界激励幅值 F2, 得到结果如图11(c)所示.

从图11 可见, 解析结果能够大致反映 F2的临界值的趋势, 与仿真结果仅在数值上存在差异, 原因是Melnikov 方法解析得到的临界值仅为系统进入混沌的必要条件, 系统的激励幅值 F2需大于该临界值到一定程度才能首次进入混沌.此外, 从图11(a)和图11(c)可见, 系统在激励频率满足严格的倍数关系和近似的倍数关系这两种情况下, 临界值的变化趋势大致相同.

5 结 论

图11 临界激励幅值 F 2 的解析结果与数值仿真结果对比 (a) ω 1 =3ω2 ; (b) ω 1 =1.1 ; (c) ω 2 =0.3667 (星 号代表解析结果, 圆代表数值仿真结果)Fig.11.Comparisons of analytical results and numerical simulation results of critical excitation amplitude F 2 : (a)ω1 =3ω2 ; (b) ω 1 =1.1 ; (c) ω 2 =0.3667 (the asterisks for analytical results and the circles for numerical simulation results).

应用多尺度法得到Duffing 系统的主-超谐联合共振的近似解析解, 并用数值方法进行验证, 解析解的瞬态过程和稳态过程都具有较高精度.基于Lyapunov 第一方法研究了系统稳态周期运动的稳定性, 发现系统最多存在2 个稳定解支和1 个不稳定解支.讨论了非线性刚度对稳态周期运动的幅频响应的影响, 结果表明非线性刚度硬化和软化均能够降低响应的幅值, 并且刚度软化能够进一步影响响应的多值性.基于Melnikov 方法对刚度软化系统进行全局分析, 得到系统的异宿轨道横截相交的必要条件, 发现较大的激励幅值、小阻尼和刚度软化容易使系统进入混沌.选取一组符合主-超谐联合共振条件的参数进行数值仿真, 发现两个临界值.当激励幅值通过第一个临界值时, 系统的异宿轨道破裂, 以激变方式进入混沌态; 激励幅值通过第二个临界值时, 系统再次发生激变, 出现新的混沌吸引子.对多种频率组合的参数分别利用数值仿真和Melnikov 方法解析得到系统首次进入混沌的临界值, 对两种方法所得结果进行了对比, 解析结果在一定程度上能够真实反映系统的临界值变化的趋势, 数值仿真的结果也验证了解析结果的正确性.