一种改进的KBM法求解非线性振动方程

2018-08-01王磊佳张鹄志祝明桥

振动与冲击 2018年13期

王磊佳, 张鹄志, 胡辉, 祝明桥

(1.湖南科技大学土木工程学院，湖南湘潭 411201;2.湖南科技大学结构抗风与振动控制湖南省重点实验室，湖南湘潭 411201)

非线性振动理论的研究对象主要包括不同参数、初始条件对系统解的影响，不同振动系统的振动规律，分析周期解的形式，研究解的稳定性等方面的问题[1-4]。近似解析法作为研究非线性振动解的主要方法之一[5-6]，主要通过应用数值计算方法的思想构建近似解，且研究对象仅限于弱非线性振动系统，其中具有代表性的方法有：Lindstedt-Poincare(L-P)法[7]、多尺度法[8]、平均法[9]、KBM法[10]。

关于非线性振动系统的求解，当前国内外学者进行了不少研究，如：Cheung等[11]提出了一种改进的L-P方法求解强非线性系统；杨志安等[12]将理论与实践结合总结出了一种对多尺度法的改进理论；陈立群等[13]将平均法应用到了求解多自由度的非线性振动系统；Lim等[14]通过改进Mickens迭代方法求解非线性振动方程拓宽了εA2的取值区间；He[15]拓展了摄动理论将其应用在求解强非线性振动系统领域；Nayfeh[16]求解出了Dffing方程固有频率的精确解。

在已有的方法中，多尺度法和L-P法都是建立在摄动法[17]的基础上，摄动法不仅可以计算周期振动，而且适用于耗散系统的衰减振动；不仅可以计算稳态响应，而且可以求出瞬态过程。平均法作为一种新的计算方法，它与摄动法相比，可以避免摄动法的许多繁琐的中间过程，快速获得结果，但因为计算中精度仅能达到与ε同阶的一次近似解，从而难以满足高精度定量计算的要求。KBM法是在平均法与摄动法的基础上提出来的，吸取它们各自的优点，并在求解含有小参数ε的弱非线性振动系统方面得到了广泛的应用，但随着ε或振幅a的增大，KBM法也难以满足精度要求，因而需要展开进一步的研究。

1 改进的KBM法

弱非线性自治方程为

(1)

式中:ε为一个小参数。当ε≠0时，由于式(1)右边摄动项的存在，使得非线性方程的解中除频率为ω0的主谐波外，还含有微小的高次谐波，且振幅与频率均与ε相关，并缓慢变化，因此，运用常数变易法可以构造出该非线性自治方程的一个级数解

(2)

式中:xi(a,ψ)(i=1,2,…)为ψ的以2π为周期的函数，且a和ψ均是随时间t缓慢变化的函数，可由以下微分方程确定

(3)

(4)

式(3)和式(4)为经典KBM法的基本方程，当ε无限小时，由于该方法的二次近似解精度较高，而被广泛的应于求解非线性自治方程。随着ε的增大，该方法解的精度也会随之下降，工程中通过增加求解ε的幂次可以提高精度，但是当ε增大到一定程度时仍会出现求解精度过低甚至无法求解的问题。

为了解决这个问题，本文将式(4)的等式两边同时平方，得到：

(5)

式中频率分量ωi(a)(i=1,2,3…)未知。

经典KBM法在预设解的形式时存在某种任意性，再通过解的周期性要求以消除该任意性。这种方法在求解非线性振动问题中是行之有效的，因而可求得周期解或非周期解。经典KBM法假定频率ω与ε相关，但忽视两者的具体关系，式(4)最终在摄动过程中逐步得以确定。据此可知，本文中对式(5)的假定同样能满足这一系列要求。

为了消除长期项并保证Ai和ωi的唯一性，则函数xi(a,ψ)必定不是关于sinψ和cosψ的谐调方程，由式(3)和式(5)得

(6)

(7)

式(7)两边再次对时间t求导，整理得：

(8)

将式(2)和式(8)代入式(1)，左边整理得：

(9)

(10)

式(9)和式(10)恒等，则两式ε的同次幂系数必相等。因为xi是关于时间t的周期函数，所以可以确定频率分量ωi(i=1,2,…)。为消除长期项，令式(9)和式(10)中sinψ和cosψ项的系数为零。在实际计算中只需取级数解的前几项为近似解就可以达到较高的精度，本文算例中只取到ε2项，整理得：

引入边界条件，代入式(11)和式(12)即可求解方程的近似解和频率解。

由以上的推导过程不难发现，改进后的KBM法不仅求解过程简单，利于验算，而且适用于电算。

2 算例

Duffing方程系统是一种典型的非线性系统，工程实际中的许多非线性振动问题的数学模型都可以转化为该方程来研究，Duffing方程一般表示为

(13)

边界条件为

(14)

由初始条件可知

x=x0+εx1+ε2x2+…+εixi,

(15)

ω为非线性振动的固有频率，xi,ωi待确定。

(16)

将式(15)和式(16)代入式(13)化简得

(ω2-εω1-ε2ω2-…)(x0+εx1+ε2x2+…)+

ε(x0+εx1+ε2x2+…)3=0

(17)

应用和的立方公式，将式(17)整理成关于ε的幂函数

(18)

式(18)中ε≠0，而方程要衡等于零,则必有ε的系数分别等于零。整理得

(19)

(20)

(21)

由式(14)的初始条件可分别求得

x0=acosψ

(22)

(23)

(24)

由式(22)、式(23)和式(24)求得式(15)前两项的近似解为

(25)

(26)

其中ψ=ωt+ψ0,由Meirovitch[18]的讨论方法可以确定出a和ψ0与初始条件间的关系，令ψ0=0代入式(15)、式(25)得

(27)

a关于ε的幂级数为

a=a0+εa1+ε2a2+…

(28)

将式(28)代入式(27)得

(29)

将式(29)表示成ε的幂次函数，并忽略高于ε2的项。由初始条件ε≠0,可以推得ε前的系数等于0，由此理论可以得到

(30)

(31)

求得

(32)

将上式代入式(28)整理得

(33)

将式(33)代入式(26)和式(25)，并将其整理成关于ε的幂函数得

(34)

(35)

式(35)整理得

(36)

经典KBM法关于Dffing方程的二次近似解为

(37)

(38)

3 结果与讨论

3.1 频率ω的近似解分析

本节以求解Duffing方程为例，将本文改进后的KBM法与以下几种常见方法进行对比。

(39)

(40)

(41)

(42)

Duffing方程的精确解为

(43)

表1 几种不同方法求得的ω近似解的比较

3.2 解析解x(t)的对比分析

为对比实参数ε和振幅a分别取不同值时经典KBM法、He法、Lim法与本文改进后的KBM法求解Duffing方程二次近似解的精度，本节将四阶龙格-库塔法[19]求得的函数在一个周期内的数值解与上述四种方法求解的结果进行对比，如图1所示。其中ω0=1。

(a) a=0.1,ε=1

(b) a=0.1,ε=10

(d) a=10,ε=1

(e) a=10,ε=10

(f) a=10,ε=100

(g) a=100,ε=1

(h) a=100,ε=100

(i) a=100,ε=10 000

图1中(a)、(b)和(c)对比发现，当振幅a=0.1时，四种方法求得的近似解在ε=1时都非常接近数值解；当ε=10时，He法求得的近似解因误差较大而不能在实际工程中应用；当ε=100时，改进后的KBM法和Lim法的二次近似解依然非常接近数值解，但经典KBM法的二次近似解开始偏离数值解。由此表明，当a较小且ε较大时，改进后的KBM法和Lim法的近似解都具有较高的精度。

图1中(d)、(e)和(f)对比可以看出，当振幅a=10，ε=1时，经典KBM法的近似解已经完全偏离数值解从而不能应用于求解非线性振动方程；同时随着ε取值的增大，Lim法的近似解也开始偏离数值解；由(g)、(h)和(i)对比可知，虽然Lim法的近似解已经偏离数值解，但偏离的范围没有随着ε的增加，而明显增大，也正是因为如此，Lim法自提出以来被广泛的应用于求解此类非线性振动方程的周期解，但Lim法的迭代过程繁琐，且需要根据解的精度要求调整迭代次数；改进后KBM法与Lim法相比，不仅求解简单，且解的精度高于Lim法。

图1中(a)、(d)和(g)，(b)、(e)和(h)，(c)、(f)和(i)对比发现，当ε一定时，随着a值的增大，Lim法的近似与改进后的KBM法的近似解相比，Lim法在数值解上下摆动受振幅取值影响较大，改进后的KBM法的近似解一直非常接近数值解。由此证明，在求解Duffing方程中，改进后的KBM法不受振幅a取值的影响。