含平方阻尼项的Mathieu-Duffing系统混沌与分岔研究

解加全, 王海军, 师 玮, 张佳乐, 霍逸婷, 曹佳琳, 高 蔷

(1. 太原师范学院 数学与统计学院, 山西 晋中 030619;2. 太原师范学院 智能优化计算与区块链技术山西省重点实验室, 山西 晋中 030619;3. 太原理工大学 机械与运载工程学院, 太原 030024)

Mathieu-Duffing方程是Mathieu方程和Duffing方程的结合和合并,Duffing方程是描述共振现象、调和振动、次调和振动、拟周期振动、概周期振动、奇异吸引子和混沌现象的简单数学模型,Duffing 方程在船舶横摇、微弱信号检测、振动能量采集等问题中被广泛应用[1-5]。Mathieu方程作为 Hill 方程一种特殊形式,是一种含参数激励的线性振动模型,经常被用来处理一些参数共振现象,在参激振动研究中得到了广泛的应用[6-9]。工程实际中的许多非线性振动问题的数学模型都可以转化为Mathieu-Duffing方程来研究,如医学、脉冲爆炸现象、信号检测,船的横摇运动、结构振动、化学键的破坏等[10-14]。侯东晓等[15]利用多尺度法求解一类具有Mathieu-Duffing振子的两质量相对转动系统发生主共振-基本参数共振的分岔响应方程。牛江川等[16]研究了含分数阶微分项的Mathieu-Duffing振子在简谐激励下的亚谐共振。邢真慈等[17]运用多尺度方法研究了参数激励下带有时滞反馈的随机Mathieu- Duffing方程的主参数共振响应问题。沈建和等[18]研究了Mathieu-Duffing振子混沌轨道至任意目标周期轨道的控制问题。

自 1963年Lorenz发现第一个混沌吸引子以来,混沌理论得到了迅猛的发展,18世纪以来,科学家对天体力学、流体力学和非线性振动中的一些失稳现象的研究发现了分岔现象,分岔和混沌也是非线性科学领域重要的研究内容,在图像处理、保密通信、计算机网络、金融系统、信号处理、交通、生命科学等领域取得了迅猛发展和广泛应用[19-24]。冯进钤等[25]证明了随机Duffing单边约束系统同样存在丰富的倍周期分岔现象。楼京俊等[26]研究了多频激励下的软弹簧型 Duffing系统的混沌动力学。朱为国等[27]研究在温度场, 机械场与电磁弹性耦合作用下四边固定支撑矩形薄板的分岔与混沌运动问题。马少娟等[28]应用Laguerre正交多项式逼近法研究了含有随机参数的双势阱Duffing系统的分岔和混沌行为。振动共振现象由Landa和McClintock首次提出,目前已经有部分学者对振动共振现象进行了研究。例如,Yang等[29]研究了具有分数阶位势非线性的过阻尼系统的振动共振。彭荣荣[30]运用多尺度法研究了一类含有外激力和五次非线性恢复力的Duffing系统。李欣业等[31]运用平均法研究了Duffing-Van der Pol振子的主参数共振响应及其时滞反馈控制问题。在现代航海中,轮船受到风浪的影响,会左右横摇,严重时会使轮船倾覆。郭建斌等[32]在对船舶横摇运动研究中,使用平方阻尼Mathieu振子作为模型并加入分数阶项,提供减震思路。

综上所述,本文利用多尺度法研究了平方阻尼的Mathieu-Duffing 振子在简谐激励下的主共振响应,建立主共振幅频响应方程的解析表达式,对定常解的稳定性进行分析,应用 Melnikov 方法得到系统进入Smale马蹄意义下混沌的条件,通过数值仿真分析了系统不同参数对混沌运动的影响,所得结果为抑制混沌运动的研究方向提供了理论验证。

1 近似解析解

在船舶横摇中研究如下含有平方阻尼的Mathieu-Duffing振子模型

(1)

对系统作如下变量替换

ε=β/m,2ε=γ/m,f1=F/m

式中:ε为小参数;ω0为系统的无阻尼固有频率。则式(1)可以转换为

(2)

为研究系统的主共振情况,要求ω≈ω0,且f1为小量,即

f1=εf,ω=ω0+εσ,f=O(1),σ=O(1)

其中,σ为引入的调谐因子,则式(2)可以转换为

(3)

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

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

(4)

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

(5)

式(5)的解为

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

(6)

式中,a(T1)和η(T1)分别为慢变振幅和相位。

式(7)也可以写成复数形式

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

(7)

(8)

(9)

(10)

其中

(11)

引入

(12)

(13)

(14)

系统的一次近似解为

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

(15)

式中,a和φ由式(14)和式(15)确定。

2 定常解及其稳定性

(16)

(17)

进一步可以从式(17)和式(18)推导出稳态运动的幅频响应方程和相频响应方程分别为

(18)

(19)

下面利用Lyapunov第一方法考察定常解的稳定性,定义状态向量V=[a,φ]T,利用式(14)和式(15)构造向量函数

(20)

(21)

λ2-Pλ+Q=0

(22)

式中,P=trJ,Q=det[J]。

处渐近稳定的条件是P<0且Q>0,因此系统的稳态响应在Lyapunov意义下的稳定性条件为

(23)

3 混沌运动的Melnikov条件

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

ευ=γ/m,δε=F/m,ε=β1/m

其中ε为小参数,从而式(1)可以写成

(24)

将式(25)写成状态方程形式

(25)

(26)

当ε=0时,未扰系统为Hamilton系统

(27)

(28)

其Hamilton量为

(29)

势函数为

(30)

由于

(31)

(32)

(33)

(34)

(35)

将式(25)写成向量函数的形式

(36)

则受扰系统(25)的Melnikov函数为

(37)

整理得

(38)

从而根据Melnikov定理可知,系统出现异宿轨道横截相交的必要条件是

(39)

4 数值仿真

为了验证近似解析解的正确性,选取一组参数,根据式(14)和式(15)得到系统幅频响应曲线。图1显示数值解与解析解在振动形式上是一致的,其幅值均出现了跳跃现象,表明解析解与数值解有很高的拟合度。

图1 幅频曲线对比,其他计算参数为m=1;k=4;c=0.2;α1=0.7; γ=0.7;F=0.3;β=0.3;μ1=0.4;μ2=0.3Fig.1 Amplitude-frequency curve comparison, other calculation parameters are m=1;k=4;c=0.2;α1=0.7;γ=0.7;F=0.3; β=0.3;μ1=0.4;μ2=0.3

由式(19)计算定常解的幅频特性,再由式(24)判断对应解支的稳定性,得到定常解的幅频特性曲线如图2所示。由于系统同时受到强迫激励和参数激励的联合作用,导致系统至多存在3个周期解支,其中1支为稳定解,2支为不稳定解。

图2 定常解的幅频响应,其他计算参数为m=1;k=0.4;c=0.01; α1=0.05;γ=0.2;F=0.01;β=0.1;μ1=0.01;μ2=0.01Fig.2 The amplitude-frequency response of the steady solution, other calculation parameters are m=1;k=0.4;c=0.01;α1=0.05;γ=0.2;F=0.01;β=0.1;μ1=0.01;μ2=0.01

图3给出非线性阻尼取不同值的幅频曲线,可以清晰的看到随着非线性阻尼的增大,系统总响应的幅值减小,由此说明本系统中非线性阻尼项的存在对系统振动起抑制作用,取值愈大抑制愈强。

图3 不同非线性阻尼的幅频曲线,其他参数为m=1;k=4;c=0.2; α1=0.7;γ=0.7;F=0.3;β=0.3Fig.3 Amplitude-frequency curves with different nonlinear damping, and other calculated parameters are m=1;k=4; c=0.2;α1=0.7;γ=0.7;F=0.3;β=0.3

选取激励频率ω作为分岔参数,选取一组参数,根据数值解得到ω∈[0,1]系统的分岔图,如图4所示。随着激励频率ω变化,系统发生了倍周期分岔,且系统中间交替发生多次的倍周期分岔。

图4 激励频率ω的全局分岔图,其他计算参数为m=0.5;k=0.02; c=0.02;γ=1.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05; F=0.05Fig.4 The global bifurcation diagram of excitation frequency ωand other calculation parameters are m=0.5;k=0.02;c=0.02; γ=1.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05;F=0.05

图5给出系统在ω=0.8时的Poincare截面图,由图知当ω=0.8时,Poincare 截面有一个吸引点, 因此系统是一个典型的单周期运动。

图5 ω=0.8的Poincare截面,其他计算参数为 m=0.5;k=0.02; c=0.02;γ=1.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05;F=0.05Fig.5 The Poincare cross section when ω=0.8, and other calculated parameters are m=0.5;k=0.02;c=0.02;γ=1.2; α1=0.2;β=0.3;μ1=0.05;μ2=0.05;F=0.05

图6给出系统在ω=0.09时的Poincare截面图,相轨迹在有界区域内反复缠绕而不封闭,Poincare截面非有限点集,非封闭曲线,具有明显的分形结构,以上特征均说明此时系统处于混沌状态。

图6 ω=0.09的Poincare截面,其他计算参数为m=0.5;k=0.02; c=0.02;γ=1.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05;F=0.05Fig.6 The Poincare cross section when ω=0.09, and other calculated parameters are m=0.5;k=0.02;c=0.02;γ=1.2; α1=0.2;β=0.3;μ1=0.05;μ2=0.05;F=0.05

选取激励幅值F作为分岔参数,选取一组参数,由数值解得到F∈[0.04,1]系统分岔密度分布图,如图7所示。随着激励幅值F变化,系统发生倍周期分岔,系统从 1-1 单周期发生倍周期分岔转迁为混沌运动,且系统中间交替发生多次的倍周期分岔,最终回归到周期运动。

图7 激励幅值F的分岔密度分布图,其他计算参数为m=1;k=0.02; c=0.04;γ=0.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05; ω=0.09Fig.7 Bifurcation density distribution diagram of excitation amplitude F, and other calculation parameters are m=1;k=0.02;c=0.04;γ=0.2;α1=0.2;β=0.3;μ1=0.05; μ2=0.05;ω=0.09

图8给出系统在F=0.4时的Poincare截面图,由图知当F=0.4时, Poincare 截面有一个吸引点,系统是一个典型的单周期运动。

图8 F=0.4的Poincare截面,其他计算参数为 m=1;k=0.02; c=0.04;γ=0.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05;ω=0.09Fig.8 The Poincare cross section when F=0.4, and other calculated parameters are m=1;k=0.02;c=0.04;γ=0.2; α1=0.2;β=0.3;μ1=0.05;μ2=0.05;ω=0.09

图9给出系统在F=0.093时的Poincare截面图。相轨迹在有界区域内反复缠绕回旋而不封闭,Poincare截面非有限点集,非封闭曲线,具有明显的分形结构,以上特征均表明此时系统处于混沌状态。

图9 F=0.093的Poincare截面,其他计算参数为m=1;k=0.02; c=0.04;γ=0.2;α1=0.2;β=0.3;μ1=0.05;μ2=0.05;ω=0.09Fig.9 The Poincare cross section when F=0.093, and other calculated parameters are m=1;k=0.02;c=0.04;γ=0.2; α1=0.2;β=0.3;μ1=0.05;μ2=0.05;ω=0.09

图4表明激励频率变化引起系统倍周期分岔,并呈现不同频率下倍周期分岔的变化的趋势。图5与图6取不同频率利用Poincare截面展现倍周期分岔,验证图4分岔的正确性。

图7显示外激励幅值变化引起系统倍周期分岔,并展现出不同幅值下倍周期分岔的演化方向。图8与图9取不同激励幅值利用Poincare截面呈现倍周期分岔,验证图7分岔的正确性。

5 结 论

本文研究一类含平方阻尼项的Mathieu-Duffing系统在简谐激励下的动力学特性。利用 Lyapunov理论分析系统的幅频特性,由于系统同时受到强迫激励和参数激励的联合作用,导致系统至多存在3个周期解支,其中1支为稳定解、2支为不稳定解。激励频率和激励幅值的变化均可导致系统经倍周期分岔进入混沌运动,且激励幅值越大、激励频率越小越容易发生混沌运动。平方项阻尼项的存在对系统振动有抑制作用,取值愈大对系统振动抑制愈强。因此,适当降低激励幅值、增大激励频率有利于抑制混沌运动的发生。