国忠金 夏丽莉 张伟,3
(1.北京工业大学 机电学院, 北京 100124) (2.泰山学院 数学与统计学院, 泰安 271000)(3.机械结构非线性振动与强度北京市重点实验室, 北京 100124)
非线性动力学与振动分析对机械、结构等动力学问题研究是非常重要的,它能够全面了解和准确预测系统运动特性.近年来, 双自由度非线性振动系统的振动频率与周期响应已被广泛研究[1-5].通常,寻求此类复杂方程的精确解析解是比较困难的.因此,多自由度振动系统的频率、响应分析主要集中在近似分析方法方面,诸如:摄动法[3,4]、MLP法[5]、同伦摄动[6]、同伦分析法[7,8]等.
谐波平衡法是不受小参数约束应用最广泛的定量分析方法.诸如:钱长照[9]运用谐波平衡法研究了含有单向离合器、两滑轮及附件的轮-带驱动系统稳定稳态周期响应. 杨荣刚等[10]基于谐波平衡法研究了摆线钢球行星传动系统的基频稳态响应及动态特性. 谐波平衡法对于一阶近似解求解很方便, 但精度不高. 因此, 许多研究者将谐波平衡法进行了一些推广, 发展了一些诸如增量谐波平衡[11]、摄动-增量[12]、谐波-能量平衡[13]、线性谐波平衡[14]、牛顿谐波平衡[15,16]、余量谐波平衡[17]、多层余量谐波平衡[18]等方法. 同伦[19]是拓扑理论的一个基本概念,用于描述两个对象间的连续变化.因此,借鉴同伦思想,通过引入影射参数,可将原始非线性问题扩展为一簇易于求解的问题.余量谐波平衡[17,18]通过引入阶层参数, 融合同伦思想到谐波平衡方法中, 继而余量延拓,易获得高阶近似解.
本文针对两类双自由度质点动力学运动系统构建了其振动频率、稳态响应求解的余量谐波平衡解程序,得到系统的高阶余量谐波平衡近似解频率及响应表达式,并与已有结果进行了比较分析.
考虑如图1所示, 运动在光滑平面上的两质点系统:
假定非线性恢复力函数为:
k1(x2-x1)+k2(x2-x1)3
振动系统的动能与势能分别为:
(1)
其中,k1和k2分别为线性、非线性弹簧刚度系数,x1和x2分别为两等质量质点的位移函数.
图1 线性、非线性刚性连接的两质点动力学模型Fig.1 Two-mass dynamic model connected by linear and nonlinear stiffnesses
继而运用欧拉-拉格朗日方程,可获得运动在光滑平面上的两质点运动方程为[2,4]:
(2)
考虑如图2所示, 运动在光滑平面上且与线性和非线性弹簧连接的两质点系统:
图2 具有固定边界的两质点动力学模型Fig.2 Two-mass dynamic system connected with the fixed boundaries
类似于模型系统(2),此振动系统的动能和势能分别为:
(3)
其中,k为线性弹簧刚度系数.
同样,运用欧拉-拉格朗日方程,可获得具有固定边界的两质点动力学系统为[2,4]:
(4)
假定ω是模型系统(4)的未知振动频率,引入变量τ=ωt,得:
(5)
(6)
根据文献[3],引入如下变量:
u=x1,v=x2-x1
(7)
方程(5),(6)变为:
mω2u″+ku-k1v-k2v3=0
(8)
mω2(u″+v″)+k(u+v)+k1v+k2v3=0
(9)
且:
u(0)=x10,u′(0)=0,v(0)=x20-x10,v′(0)=0
(10)
将上述三式整理可得:
mω2v″+(k+2k1)v+2k2v3=0
(11)
v(0)=x20-x10A,v′(0)=0
(12)
基于方程(11)的对称性,其周期解具有如下基本解级数形式:
{cos[(2k+1)τ|k=0,1,2,…}
(13)
为方便计算,引入阶层参数p,并将系统稳态下解响应及振动频率设为:
v(τ)=v0(τ)+pv1(τ)+p2v2(τ)+…
(14)
其中ωi(i=0,1,2,…)为未知频率.
根据方程(11)与初始条件(12),初始谐波解设为如下形式:
v0(τ)=Acos(τ),τ=ω0t
(15)
将(15)代入(11)式,得初始余量项:
(16)
根据伽辽金法,消除久期项,易获得初始谐波近似频率及周期响应为:
v0(τ)=Acos(ω0t)
(17)
上述近似公式与振幅-频率公式[4]结果一致.
将初始近似(17)代入余量项(16)时,cos(3τ)系数非零.因此,我们将(14)代入(11)合并阶层参数p的一次系数,得:
(18)
式(18)为关于未知ω1和v1(τ)的线性方程,根据周期解级数形式(13)及初始条件,假定:
v1(τ)=a11[cos(τ)-cos(3τ)]
(19)
将(19)代入(18)式,并消除初始余量项,得:
R1(τ) =Γ1(τ)+R0(τ)
(20)
将初始余量项引入(20)式,提高解的精确性.
根据伽辽金法,消除久期项,解cos(τ)和cos(3τ)的系数方程得1-阶余量谐波近似为:
v(1)(τ) =v(0)(τ)+v1(τ)
=(A+a11)cos(τ)-a11cos(3τ),
τ=ω(1)t
(21)
将1-阶余量谐波近似(21)代入余量项(20)时,cos(5τ)系数非零. 因此,我们将(14)代入(11)合并阶层参数p的2次系数,得:
(22)
式子(22)为关于未知ω2和v2(τ)的线性方程,根据周期解级数形式(3)及初始条件,假定:
v2(τ)=a21[cos(τ)-cos(3τ)]+
a22[cos(τ)-cos(5τ)]
(23)
将(23)代入(22)式,并消除1-阶余量项,得:
R2(τ) =Γ2(τ)+R1(τ)
(24)
根据伽辽金法,消除久期项,我们解cos(τ),cos(3τ)和cos(5τ)的系数方程得2-阶余量谐波近似为:
v(2)(τ) =v(0)(τ)+v1(τ)+v2(τ)
=(A+a11+a21+a22)cos(τ)-
(a11+a21)cos(3τ)-a22cos(5τ),
τ=ω(1)t
(25)
类似于上述求解过程.一般地,k-阶余量谐波近似:
v(k)(τ)=v(k-1)(τ)+vk(τ),
v(k-1)(τ)=v(k-2)(τ)+vk-1(τ),
v(0)=Acos(τ),ω(0)=ω0,k=2,3,4,….
(26)
一方面,将方程(5),(6)相加得:
(27)
令ρ(t)=x1(t)+x2(t),方程(27)变为一个二阶线性微分方程:
(28)
初值条件为:
(29)
其解为:
ρ
(30)
另一方面,根据上述余量谐波平衡解程序,已获得变量v(t)的各阶近似响应.如:2-阶余量谐波平衡解为:
x1-x2=(x10-x20)cos(ω(2)t)+a11[cos(ω(2)t)-
cos(3ω(2)t)]+a21[cos(ω(2)t)-
cos(3ω(2)t)]+a22[cos(ω(2)t)-
cos(5ω(2)t)]
(31)
由(30),(31)得模型(4)的2-阶近似响应为:
(32)
(33)
同样地,模型系统(2)及(4)的高阶近似响应易由(26),(30)~(33)获得.
本部分结合实例分析,图解了模型系统(2)、(4)的各阶余量谐波平衡解频率以及与已有文献的结果比较.系统(4)的精确解频率[4,15]为:
(34)
针对系统(2)、(4),表1~2给出了各阶近似余量谐波平衡解频率及本文2-阶余量谐波近似结果与已有文献结果的比较. 其中相对误差定义为
(35)
从表1a和2a可以看出,模型系统(2)和(4)在稳态下的近似解与精确频率的误差随着余量谐波平衡解阶数的增大而大大减小,这表明余量谐波平衡解程序的收敛性.本文给出的2阶近似振动频率结果与已有文献相比:本文的结果比改进的振幅-频率公式[4],初始同伦分析解[8],Pade同伦近似[8],3-阶线性谐波平衡解[14],3-阶牛顿谐波平衡解[15]等在各类参数下结果更加精确,与精确解的相对误差已降低.诸如,对于系统(2),在参数m=10,k1=50,k2=-0.01,x10=-20,x20=40下,上述文献方法的解相对误差依次为3.14%,2.93%,0.13%, 0.019%, 0.018%,本文的2阶-余量谐波平衡解频率与精确解间相对误差为:0.14%.对于系统(4),在参数m=10,k=50,k1=70,k2=90,x10=20,x20=-40下,上述文献方法的解相对误差依次为2.22%, 2.22%, 0.068%, 0.008%, 0.0071%, 本文的2阶-余量谐波平衡解频率与精确解的相对误差为:-0.0057%.
表1a 模型系统(2)的各阶余量谐波平衡解频率Table 1a Various-orders residue harmonic balance solution frequencies of model system (2)
表1b 模型系统(2)的近似解频率结果比较Table 1b Comparisons of approximate frequency for model system (2)
表2a 模型系统(4)的各阶余量谐波平衡解频率Table 2a Various-orders residue harmonic balance solution frequencies of model system (4)
表2b 模型系统(4)的近似解频率结果比较Table 2b Comparisons of approximate frequency for model system (4)
为进一步图解本文获得响应的有效性,图3~5分别显示了系统(2)和(4)在不同参数下系统的时域振幅曲线比较.
图3在系统(2)参数m=5,k1=10,k2=10,x10=20,x20=30下的近似解析响应为:
x1(t)=25-4.7796cos(ω(2)t)-0.2117cos(3ω(2)t)-
0.0087cos(5ω(2)t),ω(2)=17.0663,
x2(t)=25+4.7796cos(ω(2)t)+0.2117cos(3ω(2)t)+
0.0087cos(5ω(2)t),ω(2)=17.0663.
图4在系统(4)参数m=1,k=1,k1=1,k2=1,x10=5,x20=1下的近似解析响应为:
x1(t)=3cos(t)+1.9129cos(ω(2)t)+
0.0837cos(3ω(2)t)+0.0034cos(5ω(2)t),ω(2)=10.8684,
x2(t)=3cos(t)-1.9129cos(ω(2)t)-
0.0837cos(3ω(2)t)-0.0034cos(5ω(2)t),ω(2)=10.8684.
图5在系统(4)参数m=50,k=1,k1=1,k2=1,x10=10,x20=5下的近似解析响应为:
图3a 解析近似响应x1(t)与数值解比较, 当m=5,k1=10,k2=10,x10=20,x20=30Fig.3a Comparison of analytical solution x1(t) with the numerical one for case m=5,k1=10,k2=10,x10=20,x20=30
图3b 解析近似响应x2(t)与数值解比较, 当m=5,k1=10,k2=10,x10=20,x20=30.Fig.3b Comparison of analytical approximate solution x2(t) with the numerical one for the case m=5,k1=10,k2=10,x10=20,x20=30
图4a 解析近似响应x1(t)与数值解比较, 当m=1,k=1,k1=1,k2=1,x10=5,x20=1Fig.4a Comparison of analytical approximate solution x1(t) with the numerical one for the case m=1,k=1,k1=1,k2=1,x10=5,x20=1
图4b 解析近似响应x2(t)与数值解比较, 当m=1,k=1,k1=1,k2=1,x10=5,x20=1Fig.4b Comparison of analytical approximate solution x2(t) with the numerical one for the case m=1,k=1,k1=1,k2=1,x10=5,x20=1
图5a 解析近似响应x1(t)与数值解比较, 当m=50,k=1,k1=1,k2=1,x10=10,x20=5Fig.5a Comparison of analytical approximated solution x1(t) with the numerical one for the case m=50,k=1,k1=1,k2=1,x10=10,x20=5
图5b 解析近似响应x2(t)与数值解比较, 当m=50,k=1,k1=1,k2=1,x10=10,x20=5Fig.5b Comparison of analytical approximate solution x2(t) with the numerical one for the case m=50,k=1,k1=1,k2=1,x10=10,x20=5
0.0990cos(3ω(2)t)+0.0038cos(5ω(2)t),ω(2)=0.883305,
0.0990cos(3ω(2)t)-0.0038cos(5ω(2)t),ω(2)=0.883305.
从图3~5不难看出,本文给出的近似响应与数值解在系统(2)、(4)不同参数下都吻合得相当好.
基于谐波平衡,通过融合同伦思想优势,构建了不含小参数,适用于求解双自由度非线性两质点振动系统的高阶近似余量谐波平衡解程序.不同于线性谐波平衡与牛顿谐波平衡方法,本文解程序在每一阶近似中均消除了上一阶的谐波余量,高阶振动频率近似表达仅需初始谐波近似,不需根据前一阶近似进行调整.理论上,任何精度的高阶近似均能依次获得.相应地,系统的高阶近似响应表达易通过求解、整合获得.结果显示,本文获得2-阶余量谐波平衡近似振动频率比已有的在改进振幅-频率公式法、初始同伦分析法、3-阶线性谐波平衡法、3-阶牛顿谐波平衡法等方面更加精确,与精确解的相对误差在不同参数下均降低了.相应地,本文获得的2-阶解析近似响应与数值解吻合得相当好.