机械系统中发生的碰撞或振动不仅会产生大量噪音，还将严重影响系统的安全性与可靠性[1-3]；因此，作为机械设备中重要的隔振减振元件，橡胶等黏弹性材料因其同时具备超弹性和黏弹性而被广泛应用于各种振动系统中，如机车车辆二系悬挂橡胶堆、空气弹簧橡胶囊、钢轨扣件弹性垫板．因Kelvin模型(弹簧与阻尼并联)和麦克斯韦模型(弹簧与阻尼串联)均无法同时反映橡胶材料的松弛特性与蠕变特性[4]，所以目前常用二者的耦合模型表征橡胶等黏弹性材料的力学特性，如Zener模型[5]、Berg模型[6]、Dzierek模型[7]及分数导数模型[8-9]．相比其他复杂模型，Zener模型具有较少的系统参数和较低的动力学计算难度，是开展橡胶隔振系统非线性动力学特性研究的较佳模型．文献[10]采用谐波平衡法求解了Zener非线性黏弹材料的幅频特性，并进行了相关试验验证．文献[11]分析了非线性Zener隔振系统在多种分岔诱导下周期运动的共存与转迁机理．虽然目前已有不少文献涉及黏弹性模型动力学特性的分析[4-12]，但针对黏弹性系统稳定性与分岔集的研究仍比较少见．该现状主要由于黏弹性模型中半自由度节点的耦合作用，增大了此类系统的求解与判稳难度．文献[13]采用分段线性麦克斯韦模型表征抗蛇形减震器，并对比分析了多种抗蛇行减振器节点刚度对高速列车运动稳定性的影响规律．文献[14]提出了一种改进Mathieu方程不稳定边界的获取方法．文献[15]针对一类Duffing系统，以切线铅直为鞍结分岔点的判定依据获取了系统的鞍结分岔点，并分析了常数激励与简谐激励对系统鞍结分岔集的影响规律．考虑到非线性幅频响应中出现的分岔类型基本都是鞍结分岔，而鞍结分岔对系统的影响不仅表现在幅频响应曲线鞍结分岔点处的突变上，还表现在系统因鞍结分岔产生的多不变集共存上．欲降低工程实际中此类运动状态突变对系统寿命、稳定性和安全性的不利影响，应尽量避免机械系统工作在上述多不变集共存区间内．因此，有必要对这种广泛存在于幅频响应且具有重要影响的分岔点进行捕捉与分析．本研究采用非线性Zener模型表征橡胶等黏弹性隔振系统的动力学特性，将系统幅频响应稳定性与鞍结分岔集的获取联系起来．首先建立了系统的运动微分方程并无量纲化，然后利用谐波平衡法求解了系统主振动下的幅频响应方程及骨架线方程，最后借助Floquet理论与Mathieu方程相关理论获得了系统周期解稳定的边界条件，并以此为基础得到了系统在不同双参平面内的鞍结分岔集与多不变集共存区．1 系统的力学模型和运动微分方程设K为质量块与地基之间非线性弹簧的刚度系数，T为外激励时间，Ω为质量块所受外激励频率．图1所示为非线性Zener隔振系统的力学模型．质量块M通过麦克斯韦模型和弹性力Fr=K(X+εX3)的非线性弹簧并联在基座上，并因谐激励Fsin(ΩT)的作用而垂向振动．节点是位于Zener模型左侧麦克斯韦模型中刚度系数为K0的线性弹簧与阻尼系数为C0的阻尼中间的无质量质点．以系统静平衡位置为坐标原点，X 和Y分别表示质量块和节点的位移．10.13245/j.hust.238700.F001图1非线性Zener隔振系统力学模型图1所示系统的运动微分方程为MX¨+K0(X-Y)+K(X+εX3)=Fsin(ΩT);K0(Y-X)+C0Y˙=0. (1)经无量纲化后，系统运动微分方程变形为x¨+x-y+μkx+μkk1x3=psin(ωt);y-x+2ξy˙=0, (2)式中：x=K0X/Fs，y=K0Y/Fs，Fs为外激励幅值F的对照参数，用以体现橡胶隔振系统的幅值相关性；系统固有频率Ω0=K0/M；无量纲时间t=Ω0T；无量纲激励频率ω=Ω/Ω0；刚度比μk= K/K0；阻尼系数ξ=C0/(2MΩ0)；非线性系数k1=εf2，f=Fs/K0；激励幅值p=F/Fs．2 系统的幅频特性与骨架曲线设图1所示系统的主振动为x(t)=Asin(ωt+φ);y(t)=Bsin(ωt+θ), (3)式中：A与φ分别为质量块的振动幅值和相角；B与θ分别为节点的振动幅值和相角；ω为无量纲化后的外激励频率．将式(3)代入式(2)，可得关于质量块振动幅值、相角及节点振动幅值、相角的耦合方程组，解耦所得方程组便可计算出质量块M在主振动下的幅频响应方程T2A2+4ξ2ω2(1+T)2A2=(1+4ξ2ω2)p2, (4)式中：T=μk-ω2+(3μkk1A2)/4；p为激励幅值；ξ为阻尼系数．令式(4)中激励幅值与阻尼系数均为0，再结合质量块振动幅值A并不常为0，可得能反映非线性Zener隔振系统等效固有频率ωe随质量块振动幅值A变化的骨架曲线方程ωe =μk+(3μkk1A2)/4. (5)基于式(4)与式(5)，利用幅频响应曲线和骨架曲线刻画非线性Zener隔振系统的基本动力学特性．考虑到橡胶等黏弹性材料的幅值敏感性[4]，侧重讨论p对非线性Zener隔振系统动力学特性的影响．图2为当系统参数取μk=1,ξ=0.15,k1=10.13245/j.hust.238700.F002图2质量块M的骨架曲线与幅频响应曲线0.3，及p=0.2，0.5，1.0，2.0，3.0时质量块M的幅频响应曲线，图中：虚线为系统骨架曲线；点线表示系统幅频响应的不稳定周期解．由系统骨架线方程可知：系统骨架曲线实质上是双曲函数ω2/μk-(3k1A2)/4=1位于第一象限的部分，并在图2中呈现为一族幅频响应曲线峰值点的轨迹．该骨架线的起点由μk决定(μk取决于K，K0及M)，而其弯曲程度还受到非线性系数k1的影响．由图2可以看出：随着p的增大，质量块幅频响应曲线整体上移、共振峰值逐渐增大、共振滞后区逐步加宽．此外，随着骨架曲线向右弯曲，质量块幅频响应曲线从近似线性(图2中p=0.2的情况)转变为出现硬弹簧特性的共振滞后区(图2中p=0.5的情况)．在系统共振滞后区的边界处(两个鞍结分岔点处)，系统幅频响应曲线发生跳跃现象．以p=3时的幅频响应为例，当ω由小增大时，质量块振动幅值在下跃点CJ处到达峰值并发生向下跳跃，最终呈现B→CJ→D→E的幅频路径，而当ω由大减小时，质量块振动幅值将以E→FJ→G→B的路径变化，在上跳点FJ处发生向上跳跃．在同一参数正反扫频的两个跳跃点CJ与FJ对应的频率范围[ωCJ,ωFJ]内，系统将呈现2个稳定周期解与1个不稳定周期解的共存，而随着p的增大，系统共振滞后区迅速扩张，使系统在较宽频率范围内呈现出稳定不变集与不稳定不变集的共存．3 系统周期解的稳定性判别系统稳定性的方法有很多，本研究采用Floquet理论[15]对非线性Zener隔振系统的稳定性进行判定．首先，在系统主共振响应(式(3))中引入一个非常小的扰动量δ^x(t)与δ^y(t)，则系统扰动解可表示为x^(t)=Asin(ωt+φ)+δ^x(t);y^(t)=Bsin(ωt+θ)+δ^y(t). (6)将式(6)代入式(2)，略去其中关于扰动量的高阶项并化简，可得基于解x的线性变分方程为δ^x″(t)+δ^x(t)-δ^x(t)/1+(2ξω)2+μkδ^x(t)+3μkk1A2[1-cos2(ωt+φ)]δ^x(t)/2=0. (7)令变量τ=ωt+φ，δ(τ)=δ^x(ωt+φ)，通过变量代换有ωδ'(τ)=δ^x'(t)，ω2δ'(τ)=δ^x″(t)，得到如下形式的Mathieu方程[14]，即    ω2δ'(τ)+[1+μk-1/1+(2ξω)2]δ(τ)+3μkk1A2(1-cos2τ)δ(τ)/2=0. (8)设式(8)特解为δ(τ)=acos(τ)+bsin(τ), (9)式中a和b为质量块扰动响应的两个幅值分量．将式(9)代入式(8)，略去其中的高阶项，并使方程两侧谐波项系数相等，有1+μk-11+(2ξω)2-ω2a+34μkk1A2a=0;1+μk-11+(2ξω)2-ω2b+94μkk1A2b=0.进而可得决定非线性Zener隔振系统周期解稳定性的判别式Δ=Ju00Jd=JuJd, (10)式中：Ju=1+μk-1/1+(2ξω)2-ω2+3μkk1A2/4；Jd=1+μk-1/1+(2ξω)2-ω2+9μkk1A2/4．结合Floquet理论与Mathieu方程的相关理论，系统周期解稳定边界对应于判别式Δ=0，当Δ0时，系统周期解是稳定的，而当Δ0时，系统周期解处于不稳定状态．通过判别式Δ=0还可将系统周期解稳定边界划分为上稳定边界Su(Ju=0)(式(10)左上角元素)与下稳定边界Sd(Jd=0)(式(10)右下角元素)两部分．系统上稳定边界Su与下稳定边界Sd所形成的区域即系统不稳定振动幅值的分布范围，系统在该区域内为不稳定状态．图3为参数μk=1,ξ=0.15,k1=0.3，系统的稳定边界及不稳定区域．图中：白色区域表示系统当前周期解是稳定的；蓝绿色区域为系统当前参数下不稳定周期解的存在范围；位于该区域上下两侧的包络线即系统的上稳定边界Su与下稳定边界Sd．由图3可以看出：系统幅频响应将在满足上稳定边界判别式Ju=0的点(ω,A)处向下跳跃，即系统上稳定边界Su和下稳定边界Sd实质是不同外激励幅值下幅频响应曲线鞍结分岔下跃点CJ与鞍结分岔上跳点FJ所构成的集合，且上稳定边界Su与下稳定边界Sd交于点pud=(ωud,0)，其中ωud={ωud|1+μk-1/1+(2ξωud)2-ωud2=0,ωud≥0}，其大小主要取决于μk，但ξ的选取会使其发生较小的偏移，如系统在μk=1,ξ=0.15时的ωud=1.022．10.13245/j.hust.238700.F003图3系统的稳定边界及不稳定区域受系统鞍结分岔诱导，系统幅频响应曲线在其共振滞后区内呈现出2个稳定不变集与1个不稳定不变集的共存，而随着激励幅值p的不断增加，系统不同激励幅值下的不稳定不变集(图3点线)形成了上述不稳定区域．在该不稳定区域内，系统对外界条件比较敏感，将根据不同的初始状态向同激励频率下共存的上下两个稳定幅值转迁．由图3还可以看出：随着p逐渐增大，系统不稳定区域自点pud处出现后便不断向高频方向扩展．因此，对于非线性Zener隔振系统，只要激励幅值p≠0，便一定存在与之对应的、或宽或窄的多不变集共存区间．再结合式(10)，系统稳定边界虽然不受p的影响，但系统共振滞后区的分布还依赖于p的选取，因此有必要研究当参数变化时，系统多不变集共存区间的分布规律．4 系统双参平面内的鞍结分岔集作为非线性动力系统幅频响应中最普遍的分岔类型，鞍结分岔不仅可以通过切线铅直这一典型图像特征捕捉，还可由前文所得稳定边界与幅频响应曲线的交点获取：联立系统幅频响应方程(式(4))与上稳定边界条件Ju=0(ω0,A≥0)，即可计算出系统当前参数下鞍结分岔下跃点CJ对应的坐标(ωCJ,ACJ)；同理，系统幅频响应曲线与下稳定边界在第一象限内的交点坐标也即系统当前参数下鞍结分岔上跳点FJ对应的坐标(ωFJ,AFJ)．一方面，该方法将系统幅频响应稳定性与鞍结分岔集关联起来，避免因采用切线铅直带来的高阶隐函数求导问题，降低了非线性系统幅频响应中鞍结分岔点的获取难度，拓展了稳定性判定结果的应用价值；另一方面，该方法基于稳定边界的划分将系统不同双参平面内的鞍结分岔集划分为上跳鞍结分岔集与下跃鞍结分岔集(其实质为不同系统参数下，系统上跳频率ωFJ与下跃频率ωCJ所构成的集合)．图4为当参数μk=1,ξ=0.15,k1=0.3时，p对系统幅频响应曲线的影响规律．图中：红色五角星连线表示不同p下质量块主共振峰的变化轨迹，蓝色三角形为当ω=0.1时质量块的振动幅值，用其在不同p下点的轨迹表征系统在低频范围内的动态特性．参数对系统幅频响应曲线影响规律的表述方式与此同理．由图4可以看出：p从0.2逐渐增大至3的过程中，系统低频范围内的振动幅值有所增加、主共振峰显著增高、非线性摇头现象愈发明显、不稳定周期解对应的频率范围逐渐加宽．10.13245/j.hust.238700.F004图4p对系统幅频响应曲线的影响图5为当参数μk=1,ξ=0.15,k1=0.3时，系统在(p,ω)双参平面内的鞍结分岔集与多不变集共存区．图中：浅灰色区域表示系统当前参数下有且仅有1个稳定周期解，两条黑色曲线分别为不同激励幅值时系统的上跳鞍结分岔集Nu与下跃鞍结分岔集Nd，图中Nu曲线和Nd曲线包络的区域即系统在(p,ω)双参平面内的多不变集共存区．由图5可以看出：随着p的增加，系统上跳鞍结分岔集Nu与下跃鞍结分岔集Nd从交点pud=(1.021 7，0)处开始向高频方向不断延伸，致使系统多不变集共存区间同步向高频方向扩张．由此可见：虽然系统此时有较大的非线性系数，但当p0.6时，系统为近似线性系统，ωFJ与ωCJ非常接近，系统不稳定周期解及多不变集共存区间占据极小的频率范围．10.13245/j.hust.238700.F005图5系统在(p,ω)平面内的鞍结分岔集与多不变集共存区图6为当参数μk=1,ξ=0.15,p=2时，非线性系数k1对系统幅频响应曲线的影响规律(k1∈[0.05,2])．由图6可以看出：当系统非线性由弱变强时，系统低频范围内的振动幅值逐渐减小、主共振峰值先减后增、系统不稳定周期解及多不变集共存区间所占频率范围越来越宽．由此可见：适当增大非线性系数对系统振动幅值有一定的削减作用，但过大的k1会使系统呈现出过多的硬特性，同时也会使多不变集共存现象充斥于系统的中高频区域，将严重影响系统的隔振效果．10.13245/j.hust.238700.F006图6非线性系数k1对系统幅频响应曲线的影响图7为当参数μk=1,ξ=0.15,p=2时，系统在(k1,ω)双参平面内的鞍结分岔集与多不变集共存区．由图7可以看出：在激励幅值p=2的情况下，足够小的非线性系数也会使系统幅频响应呈现出多不变集共存．因为系统Nu曲线与Nd曲线的交点坐标与k1无关，所以系统上跳鞍结分岔集Nu与下跃鞍结分岔集Nd仍从点pud=(1.021 7,0)开始向高频方向扩张．随着k1的增加，系统Nd曲线迅速向高频方向延伸，而Nu曲线向高频方向的延伸在k1达到一定程度后逐渐放缓．10.13245/j.hust.238700.F007图7系统在(k1,ω)平面内的鞍结分岔集与多不变集共存区图8为当参数ξ=0.15,k1=0.3,p=2时，μk对系统幅频响应曲线的影响规律(μk∈[0.02,5])．由图8可以看出：当μk很小时，系统主共振峰也出现在低频区附近，使系统低频区的振幅显著增大．随着μk逐渐增大，系统等效固有频率ωe逐渐远离低频区，主共振峰值与低频范围内系统的振动幅值明显减小．结合式(2)可知：μk与k1同时作用于质量块位移的三次方，所以μk增大也会加剧系统的非线性摇头现象、拓宽系统不稳定周期解的存在范围．因此，增大μk可使弱非线性Zener隔振系统呈现出较强的非线性特性，而减小μk也可削弱强非线性Zener隔振系统骨架曲线的弯曲程度．10.13245/j.hust.238700.F008图8刚度比对系统幅频响应曲线的影响图9为当参数ξ=0.15,k1=0.3,p=2时，系统在(μk,ω)双参平面内的鞍结分岔集与多不变集共存区．结合图3与式(10)可知：μk是影响系统上跳鞍结分岔集与下跃鞍结分岔集交点的主要因素．当μk=0时，交点频率ωud=0，即μk从0增大到5的过程中，系统上跳鞍结分岔集Nu与下跃鞍结分岔集Nd是从原点开始向高频方向延伸的，系统上跳频率ωFJ与下跃频率ωCJ均迅速增大，多不变集共存区间也随之向高频方向逐步偏移．因此，μk的选取从很大程度上决定了系统上下稳定边界的交点和系统Nu曲线与Nd曲线的交点．10.13245/j.hust.238700.F009图9系统在(μk,ω)平面内的鞍结分岔集与多不变集共存区图10为当参数μk=1,k1=0.3,p=2时，阻尼系数ξ对系统幅频响应曲线的影响规律(ξ∈[0.02,0.50])．由图10可以看出：ξ对系统低频范围内振动幅值的影响较小，但对系统幅频响应主共振峰与多不变集共存区间有较大的影响．当ξ从0.02增加至0.50的过程中，系统主共振峰值与多不变集共存区间长度均先减后增，因而可基于此寻找到一个最佳阻尼系数ξopt，使系统振动幅值及多不变集共存区间达到最小．10.13245/j.hust.238700.F010图10阻尼系数对系统幅频响应曲线的影响图11为当其他参数保持不变，ξ从0.06逐渐增大至0.40时，系统的幅频响应曲线簇．基于幅频响应方程，通过公式∂A(ω,ξ)/∂ω=0，可得不同阻尼系数下质量块幅频响应主共振峰随ω变化的曲线(见图11所示绿色曲线)，再由方程∂A(ω,ξ)/10.13245/j.hust.238700.F011图11系统幅频响应曲线簇(ξ∈[0.06,0.40])(∂ω∂ξ)=0，可得绿色曲线最低点的坐标．通过观察图11发现：系统任意阻尼系数下的质量块幅频响应曲线都会通过一个公共点Q=(2.258,4)，且该公共点恰好为绿色曲线最低点．因此，总可以找到一条幅频响应曲线，使其峰值点恰好为公共点Q，此时系统幅频响应主共振峰达到全域最低，称该阻尼系数为系统当前参数下的最佳阻尼系数ξopt1，称其对应的幅频响应曲线为系统当前参数下的最佳幅频响应曲线．图11所示黑色曲线即当ξ取ξopt1=0.223时，系统当前参数下的最佳幅频响应曲线，此时系统不变集共存区间也处于较窄状态．图12为当参数μk=1,k1=0.3,p=2时，系统在(ξ,ω)双参平面内的鞍结分岔集与多不变集共存区．由图12可以看出：系统上跳鞍结分岔集对ξ的改变并不敏感，但当ξ逐渐增加时，系统下跃点先快速向低频方向移动，后又迅速向高频方向转迁，并在ξopt2=0.205处形成最小的系统多不变集共存区间，此时质量块振动幅值也处于较低状态．对比发现：当μk=1,k1=0.3,p=2时，使系统主共振峰最低的最佳阻尼系数ξopt1与使系统有最小多不变集共存区间的最佳阻尼系数ξopt2相差很小，工程实际中可根据具体目标在二者之间取舍．10.13245/j.hust.238700.F012图12系统在(ξ,ω)平面内的鞍结分岔集与多不变集共存区5 结论a．系统骨架曲线实质上是关于参数μk,k1,ω和A的双曲函数在第一象限的部分；随着激励幅值的不断增大，质量块幅频响应曲线从近似线性转变为出现多不变集共存的共振滞后区．b．在系统多不变集共存区内，系统对外界条件比较敏感，其不稳定周期解将依初始状态的设定而稳定到与其共存的上下两个稳定周期解上；虽然系统稳定边界不受激励幅值p影响，但系统共振滞后区的分布还依赖于p的选取．c．利用系统上下稳定边界与幅频响应曲线的交点即系统鞍结分岔点这一特征，在系统稳定性判定的基础上顺势进行鞍结分岔集的获取．此举拓展了系统周期解稳定性判定的应用价值，同时也为非线性系统幅频响应中鞍结分岔集与多不变集共存区的获取提供了一种便捷方法．d．系统刚度比μk的选取不仅直接决定了骨架曲线的起点，还从一定程度上决定了系统上下稳定边界的交点及系统上跳鞍结分岔集Nu与下跃鞍结分岔集Nd的交点．阻尼系数ξ的改变对系统上跳鞍结分岔集Nu的影响不大，但通过不同ξ下系统幅频响应曲线簇与下跃鞍结分岔集Nd的分布，可以得到系统当前参数下的最佳阻尼系数ξopt1与ξopt2，使系统幅频响应曲线具有最低的主共振峰或最窄的多不变集共存区间．