机械系统由于设计、加工、安装中的多种因素的共同作用，使得系统中必然存在各种类型的间隙和约束．间隙和约束的存在使得机械系统在激振力作用下产生冲击振动，进而缩减机械系统的使用寿命，甚至引发安全事故，所以研究机械系统中存在的冲击振动对于机械系统的结构设计、噪声抑制以及参数优化有着重要意义[1]．间隙、约束的存在使得此类系统成为了分段光滑的动力学系统，同时表现出了丰富的非线性特性．针对碰撞发生过程的描述是研究冲击振动系统的关键，一种经典的研究碰撞发生过程的方法是假设碰撞发生时间极短，发生碰撞后速度和位移的方向立即改变，用恢复系数来描述碰撞前后速度的变化，即为刚性碰撞．以刚性碰撞描述碰撞过程的大量研究相继展开，同时取得了许多卓有成效的成果．针对含刚性约束的冲击振动系统中存在的擦碰奇异性[2]、颤碰运动[3]、低频黏滞特性以及产生机理[4]的研究为后续对该类系统的研究打下了基础．文献[5]对含双侧刚性限幅约束的单自由度振子做了较为系统的研究；文献[6]对通过刚性碰撞建模的振动落砂机的概周期运动进行了有效控制；文献[7]研究了正弦激励下含间隙约束非光滑振子的分岔特性；文献[8]研究了一类含双侧等间隙刚性约束的两自由度冲击振动系统，分析了对称结构下的非对称分岔行为；文献[9]以含刚性限幅约束的两自由度冲击振子为研究对象，详细分析了迟滞域内周期运动的共存现象，通过胞映射的原理得到了共存周期运动的种类．随着研究的逐渐深入，发现使用恢复系数描述冲击过程的刚性约束力模型很难符合实际的碰撞过程，此时提出了使用弹性约束来描述碰撞过程，这不仅包含了碰撞过程中的受力过程，也包含形变过程．许多研究针对含弹性约束的单自由度冲击振动系统展开，探讨了包括约束刚度对于冲击振动特性的影响[10]，弹性约束的改变对于系统全局动力学特性的影响[11]．文献[12]针对工程领域常见的悬臂结构建立了其动力学模型，分析其从周期运动经擦边分岔向混沌运动转迁的方式．在单自由度冲击振动系统研究的基础上，许多研究针对两自由度含弹性约束的冲击振动系统展开[13-14]，重点分析了周期运动的转迁规律．刚性约束力模型较好的描述了碰撞前后的能量变化过程，但缺少形变过程；弹性约束力模型较好地描述了受力和形变过程，但缺少能量损耗．研究者们希望提出一种能更加贴近描述碰撞过程中形变过程和能量损耗的约束力模型，随着赫兹接触理论的发展，发现通过非线性赫兹力描述的动力学模型和实验中的实测数据有很好的符合度，尤其适用于低速碰撞的情况．文献[15]以自由落体小球为研究对象，比较了不同约束力模型对于小球弹跳的影响．文献[16]在单自由度冲击振动系统中比较了赫兹约束力建模后的动力学系统和实验测量的结果，发现两者具有很好的符合度．文献[17]分析了非线性赫兹约束在悬臂梁动力学系统建模中的应用，发现了分岔中的迟滞现象．文献[18]通过单自由度漂移振子的研究发现了不同约束力模型对于系统动力学特性存在显著影响．文献[19]分析对比了在单自由度齿轮冲击振动系统中不同接触力模型对于齿轮动力学行为的影响．文献[20]在汽车驱动系统中建立了非线性赫兹力模型下的动力学系统，分析了相对接触位移和速度之间的关系，比较了线性接触刚度模型和赫兹接触刚度模型下振动冲击力和时频响应的差异．可见在以往的研究中，研究者们对于冲击振动系统研究中所采用的约束力模型主要为刚性约束和弹性约束，对赫兹约束力模型研究相对较少．少量针对赫兹约束模型的研究中主要以单自由度振子为主，对多自由度赫兹约束模型研究较少．在模型的研究中多借助单参数分岔图、相图等传统的动力学分析方法，而这很难获得系统参数对于系统动态特性的影响．本研究以一类含间隙两自由度冲击振动系统来模拟机械系统中各零部件之间的相互冲击，通过选取两组参数进行双参数协同仿真，对比分析了刚性约束力、弹性约束力、非线性赫兹约束力模型下系统动力学特性的异同，重点研究了非线性赫兹约束力模型下各类周期冲击振动在(ω,δ)-平面上的存在区域及系统参数对其影响．1 力学模型图1为三类不同约束力模型下的两自由度冲击振动系统．在该系统中，质块Mi(i=1，2)与支撑基础之间通过线性弹簧阻尼元件Ki和Ci完成连接固定，此外质块M1和M2通过线性弹簧阻尼元件K3和C3完成连接．质块Mi受到简谐激振力Pisin(ΩT+τ)作用，其中：Pi，Ω，T和τ分别为系统幅值、激励频率、时间和初始相位角．系统中用Xi表示质块Mi的位移，B表示质块M1和M2之间存在的间隙值，那么随着质块运动幅值的增大，当质块M1和M2之间存在的位移差与间隙值相等即X2-X1=B时，发生相互冲击．图1(a)为刚性约束，表示M1和M2于X2-X1=B处发生刚性冲击，冲击过程中的能量损失由速度恢复系数R确定，质块Mi碰撞前后的瞬态速度分别用X˙i-和X˙i+表示；图1(b)为弹性约束，表示M1和M2于X2-X1=B处发生弹性冲击，其约束刚度为K0；图1(c)使用非线性赫兹接触模型描述M1和M2于X2-X1=B处的相互冲击，约束刚度和阻尼分别为Kh和Ch．10.13245/j.hust.239089.F001图1不同约束力下两自由度冲击振动系统的力学模型引入无量纲化和参数位移xi，时间t，阻尼系数ζ，频率ω，固有频率ωn，质量分配比μm，阻尼分配比μc和μc3，刚度分配比μk和μk3，弹性约束刚度分布比μk0，激振力分布比f20，无量纲间隙δ．xi=K1Xi/(P1+P2);  t=Tωn；ζ=C1/(2K1M1);  ω=Ω/ωn;ωn=K1/M1；μm=M2/(M1+M2)；μc=C2/(C1+C2)；μk=K2/(K1+K2)；μk0=K0/(K1+K0)；μc3=C3/(C1+C3)；μk3=K3/(K1+K3)；  f20=P2/(P2+P1)；δ=K1B/(P1+P2)．(1)图1(a)所示的刚性冲击系统力学模型的无量纲运动方程如下：      100μm1-μmx¨1x¨2+2ζ11-μc3-2ζμc31-μc3-2ζμc31-μc32ζ(μc1-μc+μc31-μc3)∙x˙1x˙2+11-μk3-μk31-μk3-μk31-μk3μk1-μk+μk31-μk3x1x2=           1-f20f20sin(ωt+τ)  (x2-x1δ); (2)x˙1+=[(1-μm)-μmR]x˙1-+μm(1+R)x˙2-;    x˙2+=(1+R)(1-μm)x˙1-+[μm-R(1-μm)]x˙2-(x2-x1=δ). (3)建立如图1(b)所示的弹性冲击系统力学模型的无量纲运动方程如下      100μm1-μmx¨1x¨2+2ζ11-μc3-2ζμc31-μc3-2ζμc31-μc32ζ(μc1-μc+μc31-μc3)∙x˙1x˙2+1/(1-μk3)-μk3/(1-μk3)-μk3/(1-μk3)μk/(1-μk)+μk3/(1-μk3)∙x1x2+-f0f0=1-f20f20sin(ωt+τ), (4)式中： f0=[μk0/(1-μk0)](x2-x1-δ)    (x2-x1δ);0    (x2-x1≤δ). (5)建立如图1(c)所示的非线性赫兹冲击系统力学模型的无量纲运动方程如下      100μm1-μmx¨1x¨2+2ζ11-μc3-2ζμc31-μc3-2ζμc31-μc32ζ(μc1-μc+μc31-μc3)x˙1x˙2+11-μk3-μk31-μk3-μk31-μk3μk1-μk+μk31-μk3x1x2+-fhfh=1-f20f20sin(ωt+τ), (6)式中fh为无量纲化的非线性赫兹接触力，fh=      [μkh/(1-μkh)](x2-x1-δ)3/2{1+[2ζμch/(1-μch)](x˙2-x˙1)}    (x2-x1δ);0    (x2-x1≤δ), (7)μkh=K¯h/(K¯h+K1)，K¯h=Kh(P1+P2)/K1，μch=C¯h/(C¯h+C1)，C¯h=Ch(P1+P2)．部分无量纲参数具有固定的取值范围：μm∈(0,1)，μc∈(0,1)，μk∈(0,1)，μk0∈(0，1)，μc3∈(0,1)，μk3∈(0,1)，f20∈(0,1)，为研究系统动力学特性与系统参数的关联关系提供了便利．为便于分析两质块相互冲击的模式类型及其存在区域，使用p/n来表示不同的周期运动，其中：n为系统的振动周期与激振力周期的比值(n=1，2，…)；p为两质块在一个振动周期内发生相互冲击的次数(p=0，1，2，…)．须要特别指出的是：当n=1时，将p/1(p≥1)定义为基本周期运动；而0/1向1/1的转迁规律与p/1向(p+1)/1(p≥1)相类似，从而将0/1也纳入基本周期冲击振动群．为识别双参数平面上各类冲击振动的模式类型及其存在区域，定义以下Poincaré截面：    σp={(x1,x˙1,x2,x˙2,t)∈R4×T|x2-x1=δ,x˙2-x˙10};    σn={(x1,x˙1,x2,x˙2,t)∈R4×T|Δx=x2-x1,mod(t=2π/ω)}. (8)当系统呈现周期振动时，Poincaré截面σn上不动点的个数即为周期数n；Poincaré截面σp上不动点的个数即为碰撞次数p．对应Poincaré截面σp的映射称为冲击映射，可以表示为X(i+1)=f(X(i),μ)，(9)式中：X(i)=(x˙1+(i),x2(i),x˙2+(i),τ(i))T；X(i+1)=(x˙1+(i+1)，x2(i+1),x˙2+(i+1),τ(i+1))T；X∈R4；μ为系统参数，μ∈R8．2 不同约束模型的系统动力学特性使用不同约束力模型来描述冲击振动过程，在形变过程、能量变化上存在诸多不同．为了分析对比刚性约束、弹性约束、非线性赫兹约束在同一两自由度冲击振动系统中对动力学特性的影响，选取基准参数：μm=0.5，μc=0.5，μk=0.5，μk0=0.99，μc3=0.5，μk3=0.5，μch=0.5，μkh=0.99，f20=0，ζ=0.025，R=0.8．选取对系统动力学特性影响较大的两个参数ω和δ作为分岔参数，在(ω,δ)-参数平面上可以得到各类周期运动的分布区域．在(ω,δ)-参数平面内，使用不同颜色表示不同的周期振动类型p/n．须要指出的是：p¯/1表示刚性冲击下带有黏滞特性的完整颤碰振动，p˜/1表示非完整颤碰振动，图中未标识的灰色区域多为混沌运动以及少量的概周期运动和未辨识的周期运动．图2(a)～(c)分别为刚性约束、弹性约束以及非线性赫兹约束下系统在双参数平面上各类周期运动的分布区域，图2(d)为图2(b)的局部放大图．如图2(a)～(c)，在整个(ω,δ)-参数平面对比分析三种约束下周期运动的分布区域，可以发现三种模型在参数平面的左上角和右上角主要为无冲击周期运动0/1；中高频率域周期运动的模式类型和分布区域大致相同，位于左上角的0/1运动和1/1运动分岔线附近存在参数岛，其中包含的主要亚谐运动的模式均为1/2，2/2，2/4；三种模型的差异主要体现在位于参数平面左下角的低频率、小间隙参数域．10.13245/j.hust.239089.F002图2(ω,δ)-参数平面周期运动的模式类型和存在区域在刚性约束下，如图2(a)所示，低频率、小间隙参数域存在基本周期运动群p/1，颤振运动p˜/1，以及黏滞运动p¯/1．随着ω或δ的减小，p/1运动经Real Grazing分岔进入(p+1)/1运动；当冲击振动次数p不断增大且最小冲击速度无限接近但不等于0时，系统进入p˜/1；当系统越过Sliding分岔线进入p¯/1运动时，其表现为p足够大且最小冲击速度为0．图3(a)描述了当间隙值δ=0.1时，系统的单参数分岔图．10.13245/j.hust.239089.F003图3当δ=0.1时系统的单参数分岔图弹性约束条件下，其低频率、小间隙参数域的局部放大如图2(d)所示，相较于刚性约束最大的区别在于无黏滞运动的存在区域；此外，基本周期运动群p/1(p≥2)的带状域随着冲击次数p的增大显著收窄；未识别的灰色区域占据了较大区域，结合单参数分岔图(图3(b))可以发现其中主要为混沌运动．在极低频率域，主要为碰撞次数极高的基本周期运动(p≥10)．非线性赫兹约束条件下，如图2(c)所示，其在低频率、小间隙参数域基本周期运动群p/1(p≥2)的带状域的分布与弹性约束条件下类似，同时也不存在黏滞运动的存在域；在极小频率域范围基本周期运动1/1和2/1占据了较大区域，未识别的灰色区域相对减少．图3(c)所示的单参数分岔图中可以清晰观察到周期运动的转迁过程．综合以上分析，非线性赫兹约束作为一种能较为准确描述冲击振动形变过程和能量损耗的力学模型，其动力学特性与刚性约束、弹性约束在低频率、小间隙参数域存在显著差异，而机械振动系统的低频特性以及低频异常振动的发生一直是冲击振动系统研究的重点和难点，故而有必要对其低频动力学特性进一步深入分析．3 赫兹约束的系统低频动力学特性非线性赫兹约束条件下，系统在低频率、小间隙参数域表现出了其特有的动力学特性，这里对其中周期运动的模式类型、存在区域、转迁机理等进一步深入分析．图4(a)为图2(c)低频率、小间隙参数域的局部细节描述；图4(b)为分别通过激振力频率递减ω⃖和递增ω⃗进行数值仿真后的合成双参数平面图，用符号ω⃡表示，其中深灰色区域为两种相邻基本周期运动共存的区域．从图4(a)可以观察到：基本周期运动1/1，2/1，3/1，4/1主要呈带状分布，同时在其带状存在域中可能嵌入了由周期倍化分岔而产生的2p/2运动，如3/1运动中嵌入的6/2运动、4/1运动中嵌入的8/2运动；之后随着激振力频率的减小，基本周期运动5/1，6/1，7/1的存在域不再呈现为完整的带状分布，各基本周期运动的存在域相互分割．从图5(a)所示的单参数分岔图中可看到：当ω减小时，相邻基本周期运动主要经实擦边分岔进行转迁，其表现为周期数不变，碰撞增加一次．10.13245/j.hust.239089.F004图4非线性赫兹约束低频率、小间隙参数域局部放大图10.13245/j.hust.239089.F005图5当δ=0.15时舌状域的单参数分岔图图6描述了2/1运动经实擦边分岔转迁为3/1运动的相图．相邻基本周期运动p/1和(p+1)/1转迁的过程中存在舌状域和迟滞域两类转迁区域，图4(b)中可以清晰观察到这两类区域．以2/1运动向3/1运动转迁过程中出现的舌状域为例，随着间隙值δ的减小，2/1运动经周期倍化分岔产生了亚谐运动4/2，4/2运动经实擦边分岔产生5/2运动，5/2运动经周期倍化分岔产生10/4运动，该倍化序列使得周期运动进入混沌运动．10.13245/j.hust.239089.F006图6δ=0.15，2/1运动经Grazing分岔转迁为3/1运动相图图7为当δ递减时的单参数分岔，可以清晰观察到该转迁过程，此时舌状域中主要的亚谐运动为4/2和5/2．类似地，3/1运动转迁为4/1运动舌状域中的主要亚谐运动为6/2和7/2．可见相邻基本周期运动转迁过程中的舌状域的上边界一般为p/1运动的周期倍化分岔线，包含的亚谐运动的一般规律为2p/2和(np+1)/n．10.13245/j.hust.239089.F007图7当ω=0.7时系统的单参数分岔图相邻基本周期运动转迁过程中出现的另一类转迁区域为迟滞域，其产生是由于擦边分岔和鞍结分岔转迁过程中的不可逆．p/1运动经实擦边分岔转迁为(p+1)/1运动，而(p+1)/1运动经鞍结分岔转迁为p/1运动．通过ω递减和ω递增得到的共存区域即为迟滞域．(p+1)/1运动的鞍结分岔线为迟滞域的上边界，p/1运动的实擦边分岔线为迟滞域的下边界．以基本周期运动2/1和基本周期运动3/1转迁过程中的迟滞域HR2为例，在图4(b)中可以清晰地观察到该区域．图8为当δ=0.185时，横向穿越迟滞域HR2的合成单参数分岔图，在迟滞域的范围内，2/1运动和3/1运动共存．10.13245/j.hust.239089.F008图8当δ=0.185时横向穿越迟滞域合成单参数分岔图图9为在该单参数合成分岔图中同一点10.13245/j.hust.239089.F009图9当δ=0.185，ω=0.61时，迟滞域内基本周期运动共存(δ,ω)=(0.185，0.61)处采用不同的初值进行数值仿真得到的相图，其中图9(a)为3/1运动，图9(b)为2/1运动．4 赫兹约束模型的参数匹配如图1(c)所示非线性赫兹接触冲击振动系统中基准参数与前文中保持一致，这里分别考虑改变阻尼系数ζ和质量分配比μm对于系统动力学特性的影响．在取值范围内分别选取不同参数，分析在(ω,δ)-参数平面上周期运动的分布区域的变化．阻尼系数ζ是表征系统能量耗散的重要参数．在基准参数中ζ=0.025，这里分别取ζ=0.015，ζ=0.05和ζ=0.1，分析其对周期运动分布区域的影响，如图10所示．当ζ选取较小值0.015时，如图10(a)所示，参数平面左下角基本周期冲击振动p/1(p2)的带状域隆起加剧，分形特征更加明显，极低频区域内呈现的1/1和2/1基本周期冲击振动向低频方向压缩．相邻基本周期冲击振动间夹杂的舌形域面积增大，舌形域中亚谐冲击振动的模式类型更加丰富．可见较小的阻尼系数会增加亚谐运动、概周期运动和混沌运动的存在区域．当ζ取较大值ζ=0.05和ζ=0.1时，如图10(b)和图10(c)所示，在右上角和左上角的区域，为0/1运动，其他区域为p/1(p0)运动．左下角p/1(p2)基本周期冲击振动最大冲击次数减小，极低频区域内呈现的1/1和2/1基本周期冲击振动向高频区域延伸而存在域增大．舌形域存在区域减小，其中包含的亚谐运动、未识别的灰色区域也随之减小．10.13245/j.hust.239089.F010图10不同阻尼系数ζ下(ω,δ)-平面内周期运动的分布区域可见当ζ增大时，基本周期冲击振动的碰撞次数p减小，极低频率域内基本周期运动1/1和2/1向高频区域延伸，其存在区域随之扩大，舌形域减小，亚谐运动、概周期运动及混沌运动的存在区域减小．故较大的阻尼系数可以使得系统在较大参数域处于稳定的基本周期运动机械系统中，各零部件的质量常常不相等，质量分配比μm是表征系统中各零部件质量分配的关键参数，研究其对系统动力学特性的影响具有重要意义．在μm的取值范围μm∈(0，1)内分别选取μm=0.1，0.3，0.7，周期运动在双参数平面内模式类型及其存在区域如图11所示．当μm=0.1时，通过对比图11(a)和图2(c)可以观察到，在右上角的高频区域和左上角的大间隙区域，两质块呈现为0/1无冲击振动，左上角无冲击振动0/1的存在区域向低频方向压缩而缩小，右上角无冲击振动0/1的存在区域向低频方向延伸而扩大；系统在右下角的高频小间隙区域出现了亚谐冲击振动(1/2，2/2，3/2和6/3)；左下角存在的基本周期冲击振动群消失，极低频区域内存在的1/1运动延伸扩大至ω=0.500；整个参数域中未标识的灰色区域急剧减少．10.13245/j.hust.239089.F011图11不同质量分配比μm下(ω,δ)-平面内周期运动的分布区域当μm=0.3时，如图11(b)所示，系统在右上角高频区域中呈现的0/1无冲击振动被分割为两部分，完整性遭到了破坏；中频区域的1/1基本周期冲击振动被0/1无冲击运动代替；中高频区域内，基本周期冲击振动1/1和无冲击振动0/1的转迁边界上嵌入了亚谐冲击振动(主要为1/2)和未标识的灰色区域；基本周期冲击振动群的存在区域向高频大间隙区域延伸而扩大，极低频区域内的1/1运动向高频区域延伸至ω=0.250附近．当μm=0.7时，如图11(c)所示，系统右上角高频较大间隙区域存在的0/1无冲击振动向低频区域延伸而扩大，左上角低频较大间隙区域内存在的0/1无冲击振动向低频区域压缩而缩小；右下角的高频小间隙区域出现了未标识的灰色区域；低频小间隙区域内，基本周期冲击振动群的存在区域被压缩，向低频区域延伸而减小，相邻基本周期冲击振动间嵌入较大的亚谐冲击振动区域和灰色的未标识区域，极低频区域内1/1运动的存在区域被压缩至ω=0.050附近．可以发现：随着质量分配比μm的增大，基本周期振动的存在区域由于大量亚谐运动和未识别灰色区域的嵌入而显著减小，亚谐冲击振动的模式类型更加丰富，极低频区域内存在的1/1基本周期冲击振动的存在区域向低频区域压缩．可见质量分配比的改变对于系统周期运动特性存在显著影响，较小的质量分配比可以使系统在较大参数域处于稳定的基本周期运动．5 结论本研究对比分析了含刚性约束、弹性约束、非线性赫兹约束下两自由度受迫振动系统的冲击振动特性，通过多目标、多参数协同仿真获得了不同约束下系统周期运动的分布区域，重点分析了非线性赫兹约束下系统参数对于系统动力学特性的影响．a. 不同约束力模型描述的碰撞过程会对系统周期运动的模式和分布区域产生影响，特别是在低频率、小间隙参数域存在显著差异．黏滞运动仅存在于刚性约束条件下，弹性约束和非线性赫兹约束下基本周期运动的存在区域大致相同，但其低频转迁规律存在差异．b. 对非线性赫兹约束下系统的低频特性进行了深入分析．在基本周期运动转迁的过程中存在舌状域和迟滞域两类转迁区域．舌状域中主要亚谐运动的模式类型为2p/2和(np+1)/n，其主要经基本周期运动的周期倍化分岔和亚谐运动的擦边分岔产生．迟滞域产生的机理为擦边分岔和鞍结分岔的不可逆性，在迟滞域内相邻基本周期运动共存，初值的不同将导致系统处于不同的稳定基本周期运动．c. 分析了非线性赫兹约束下冲击振动系统中阻尼系数ζ和质量分配比μm对于系统冲击振动特性的影响．ζ主要对系统的低频区域产生影响．较小的ζ使得系统冲击振动加剧，基本周期冲击振动的带状域隆起加剧，亚谐冲击振动、概周期运动和混沌运动增多；当ζ增大时，基本周期运动的模式减少，亚谐运动和混沌运动也同步减少，在整个参数平面内，系统在较大参数域处于低碰撞次数的基本周期运动．当μm较小时，系统在更大区域处于p/1运动；随着μm增大，p/1运动的面积减小，舌形域中亚谐冲击振动的模式类型更加丰富，中高频区域开始出现面积较大的灰色区域，系统的周期运动变得特别复杂而难于预测．因此，较大的阻尼系数和较小的质量分配比可以使得系统在较大参数域处于更加稳定的基本周期运动．