多频激励局部非线性系统响应求解的降维增量谐波平衡法

姚红良，王重阳，王 帆，闻邦椿

（东北大学机械工程与自动化学院，辽宁 沈阳110819）

多频激励局部非线性系统响应求解的降维增量谐波平衡法

姚红良，王重阳，王 帆，闻邦椿

（东北大学机械工程与自动化学院，辽宁 沈阳110819）

针对传统增量谐波平衡法求解多频激励局部非线性系统周期响应耗时太长的问题，提出了降维增量谐波平衡方法。首先通过谐波平衡理论分析了多频激励局部非线性系统响应中各自由度各次谐波的定量对比关系，并且根据该定量对比关系使系统的维数降至与非线性自由度个数相同；其次针对降维后的复数非线性系统推导了多频增量谐波平衡法，以及原系统各自由度各阶响应的还原方法；最后利用双频激励局部非线性悬臂梁系统进行了所提方法的精度和效率验证。结果表明：该方法的精度与传统方法一致，但是在局部非线性自由度较少时其效率远高于传统方法。

局部非线性；周期响应；多频激励；增量谐波平衡法；降维

引 言

工程中具有局部非线性特性的系统或者结构很多，如与外界发生接触的梁或桁架结构、非线性油膜轴承支撑的转子系统等。分析局部非线性系统的稳态响应，对于掌握这些系统的动态特性具有重要的意义。

目前对多自由度系统进行响应分析，所采用的方法主要有两类：时域分析方法如Runge-Kutta方法、Newmark-β法等；频域分析方法如谐波平衡法（Harmonic Balance，HB）［1］、描述函数法（Describing Function，DF）［2］、增量谐波平衡法（Incremental Harmonic Balance，IHB）［3-5］等。在分析相同振动问题时，频域分析方法的速度要比时域分析方法快很多，因为这些方法能够跳过瞬态振动的分析时间而直接进入稳态；同时，频域分析方法还能够直接揭示出系统中频率特性和变化情况。

但是，在分析自由度较多的系统时，频域方法也会遇到计算耗时太多的问题，因此必须考虑采取降维措施。目前有两种应用较多的降维措施：第一种是采用模态降维的方法，应用最普遍的是固定界面模态综合法［6］，即将结构划分为线性和非线性子结构两部分，对其中的线性子结构应用模态截断，仅保留低阶模态，从而降低系统的自由度；第二种方法是利用局部非线性系统的特性——将原系统分为线性子结构和非线性子结构，用非线性子结构中各自由度的运动来描述线性子结构中各自由度的运动，从而将自由度降低。目前第二种方法已经应用于谐波平衡法［7-8］和描述函数法［2］。

目前各种降维措施主要针对单频激励局部非线性系统。实际工程中存在大量的多频激励系统［9-10］，如果采用传统的求解方法进行求解，计算耗时很长，因此必须提出降维的多频激励系统频域分析方法。目前仅有文献［10］提出了一种基于谐波平衡法的多频激励降维分析方法，基于其他频域方法的多频激励降维分析方法尚未见报道。

从算法原理上，增量谐波平衡法比谐波平衡法对强非线性问题、大范围变化参数问题的处理能力更强［11-12］。因此，本文提 出了基于 增量谐 波平衡法的多频激励降维分析方法。首先，介绍了传统多频激励增量谐波平衡方法及其不足；然后，利用局部非线性系统特性进行系统降维；再后，结合降维方法改进了传统的增量谐波平衡法，提出了多频激励降维增量谐波平衡法；最后，以双频激励局部非线性悬臂梁系统为例对该方法的精度和效率进行了验证。

1 局部非线性多自由度系统运动微分方程及增量谐波平衡法求解

1.1 局部非线性多自由度系统的运动微分方程

工程中常见的N自由度非线性系统，受不同频率激励，其动力学方程形式为

式中M，C，K为质量、阻尼、刚度矩阵；fnx，˙（）

x为非线性项；ωi为各激振频率；Fci和Fsi为对应于ωi的激励幅值。

当非线性系统处于稳态时，振动响应x中自由度k的响应可以展开为多重傅里叶级数形式

设ω= ［ω1，…，ωrT］，li=［l1，…，lr］，νi=li·，则式（2）可以简化为

1.2 多频增量谐波平衡法及其不足［9-10］

对于形如式（1）的系数为实数的方程，对x取增量εx，得到

将式（4）代入式（1），并忽略二次以上高次项，得

取τ=［ω1t，…，ωrt］T=［τ1，…，τr］T，对式（5）两边取Galerkin过程，有

即

对式（7）进行迭代，可以求得A0，从而求得x。

如果系统自由度为N个，谐波取到r次，则形成的瞬态刚度矩阵Km的维数为 ［N×（2r+1）］× ［N×（2r+1）］，因此在N比较小时，该方法的效率不成问题，但是当N比较大时，如N=20，r=25，则矩阵维数为1 020×1 020，这导致计算速度下降很快。这是传统增量谐波平衡法的不足。

2 多频激励局部非线性系统的降维

当非线性系统处于稳态时，非线性力fn（x，˙x）也是稳态的，可以展开为

重新排列方程，将与非线性相关的自由度排列到前边，即

将式（11）代入式（10）中的第一行，有

由式（13）得

因此可以得到

由式（12）和（15）可得

式（16）的实部为

式（17）可由增量谐波平衡法求解，求出¯x1后结合式（11）求出2。

3 降维后系统的增量谐波平衡法求解

对形如式（17）的复数系数方程，设

仿照式（6）可以得到：

对式（21）进行迭代，可以求得A0，从而求得x1，再结合式（11）可以求得x2，这样就得到了原系统各自由度的响应。

该方法可以称为降维的增量谐波平衡法，该方法迭代公式中Km的维数比传统多维增量谐波平衡法大为减少。

4 数值试验

4.1 数学模型

如图1所示的悬臂梁系统：设梁截面为圆形，半径为5 mm；梁总长为400 mm，分成20段；梁材料密度为7 850 kg/m3，弹性模量为2.1×1011Pa；α= 0和β=1×10-3为阻尼；频率不同的简谐激振力分别作用于第14节点和17节点，幅值分别为50和30 N；假设在第J1=11节点有限位装置，梁与限位装置作用时受力符合下式

式中e为间隙，kn为接触刚度。

图1 悬臂梁模型图Fig.1 Cantilever beam with single impact

采用有限元方法建立悬臂梁的动力学模型，梁单元采用铁木辛柯模型。节点i的广义坐标为，其中xi为竖直方向位移，θyi为转角。

采用增量谐波平衡法求解系统响应时，首先需要设定响应中频率成分，假设各次谐波取到激励频率的r倍，即

例如当激振频率ω1和ω2分别为100和140rad/s，r=4时，将ωk中的非正值频率和重复频率去掉，可以得到l1和l2的分布如图2所示，共有37个频率分量。

图2 响应中频率选取Fig.2 Thechosenfrequencyintheresponse

4.2 精度分析

4.2.1 单自由度非线性精度分析

取激振频率ω1=250rad/s，ω2=1.4ω1，kn= 2×105N/m，e=0.01m。分别采用Newmark-β方法和本文方法进行响应计算。采用本文方法时，所取得谐波次数分别按r=2，3，4来选择。计算中每个周期分为50段。本文中非线性力属于分段线性非线性力，在采用增量谐波平衡法进行积分时，可以采用对分插值方法，如文献［13-14］所示。

求得悬臂梁中第5点、第11点和第20点响应分别如图3（a），（b）和（c）所示。从图中可以看出，即使所取谐波项次数较少，本文方法也有很高的精度。

式中 δ为本文方法的综合误差，x（tk）和（tk）分别为用Newmark方法和本文方法计算得到的tk时间的响应，m为所取的响应总数。

采用该方法，得到取不同谐波次数时各点精度，列于表1。从表中可以看出，随着所选取谐波项次数的增加，本文方法的精度逐渐提高。

表1 不同谐波项时本文方法精度误差Tab.1 Precisionerrorofthepresentmethodunderdifferent harmonicterms

采用本文方法不仅可以求出非线性发生位置处的响应，还可以很方便地求出其他位置响应中的各次频率成分，如图4所示为第5点、第11点和第20点的频域响应。由图4可以看出，响应中含有多种不可公约的组合频率。

通过该方法可以直接求出激励频率变化时，系统响应中各次谐波分量的变化。例如取ω1=100～400rad/s，ω2=1.4ω1变化时，第5点和第20点的响应三维谱图分别如图5（a）和（b）所示。这是时域方法所不具有的一个优势。

图3 ω1=250rad/s时各节点时域响应Fig.3 Thetime-domainresponseofsomenodeswhenω1=250rad/s

图4 ω1=250 rad/s时各节点频域响应Fig.4 The frequency-domain response of some nodes whenω1=250 rad/s

图5 系统响应的三维谱图Fig.5 The 3D spectrogram of the system responses

4.2.2 多自由度非线性精度分析

保持边界条件、激励条件不变，假设系统有两处发生碰撞，分别是在节点5和节点11，如图6所示。系统接触刚度为kn1=kn2=5×105N/m，间隙分别为e1=0.001 m和e2=0.01 m。

图6 双碰撞位置悬臂梁模型图Fig.6 Cantilever beam with double impact

分别采用传统多频增量谐波平衡法和本文方法进行响应计算，求得悬臂梁中第5点、第11点和第20点响应分别如图7（a），（b）和（c）所示。

比较在两个非线性自由度时本文方法与Newmark方法的精度，得到取不同谐波次数时各点精度对比如表2所示。从表中可以看出，随着所选取谐波项次数的增加，本文方法的精度逐渐提高。

图7 ω1=250 rad/s时各节点时域响应Fig.7 The time-domain response of some nodes whenω1=250 rad/s

表2 不同谐波项时本文方法精度误差Tab.2 Precision error of the present method under different harmonic terms

4.3 效率分析

比较传统的多频增量谐波平衡方法和本文方法的效率。表3是进行单自由度非线性时两种方法效率的比较，表4是进行二自由度非线性时两种方法效率 的比 较。本 文所 采用 的 计 算 机 性 能：CPU 2.83G，内存3.25G，采用编程语言Matlab。

表3 单自由度非线性时运行所需CPU时间对比Tab.3 CPU time for the two methods with 1 nonlinearity location

表4 两个自由度非线性时运行所需CPU时间对比Tab.4 CPU time for the two methods with 2 nonlinearity location

从表中可以看出，在本例中本文方法在非线性个数较少、所取谐波项少时效率很高；当非线性自由度较多，且所取谐波项较多时，方法的优势逐渐减弱。这是因为随着谐波项次数增多，形成的瞬态刚度矩阵的维数也增大，造成效率下降。

但是在整体自由度增多而非线性自由度保持不变时，本文方法的瞬态刚度矩阵维数保持不变，因此效率不变；而传统IHB方法的瞬态刚度矩阵维数增大，效率下降。此时，本文方法的优越性就可以体现出来。

5 结 论

（1）针对多频激励的局部非线性系统响应求解耗时长的问题，提出了降维的增量谐波平衡法，在整体自由度较多而非线性自由度较少时，本文方法具有很高的效率；

（2）该方法由于利用了多自由度系统的局部非线性特性和稳态响应特性，将各自由度的各次谐波响应用局部非线性自由度的响应表示，从而能够使系统维数降低，且降维后系统含有原系统所有动力学特性；

（3）该方法的理论基础为稳态响应各次谐波的谐波平衡理论，因此仅能对稳态响应进行分析，而不能分析混沌运动等非稳态响应。

［1］ Liu J K，Chen F X，Chen Y M.Bifurcation analysis of aero-elastic systems with hysteresis by incremental harmonic balance method［J］.Applied Mathematics and Computation，2012，219（5）：2 398—2 411.

［2］ Wei F，Zheng G T.Multi-harmonic response analysis of systems with local nonlinearities based on describing functions and linear receptance［J］.ASME Journal of Vibration and Acoustics，2010，132（3）：0310041—0310046.

［3］ Shen Y，Yang S，Liu X.Nonlinear dynamics of a spur gear pair with time-varying stiffness and backlash based on incremental harmonic balance method［J］.International Journal of Mechanical Sciences，2006，48 （11）：1 256—1 263.

［4］ 晏致涛，张海峰，李正良.基于增量谐波平衡法的覆冰输电线舞动分析［J］.振动工程学报，2012，25（2）：161—166. Yan Zhitao，Zhang Haifeng，Li Zhengliang.Galloping analysis of iced transmission lines based on incremental harmonic balance method［J］.Journal of Vibration Engineering，2012，25（2）：161—166.

［5］ 王威，宋玉玲，李瑰贤.基于增量谐波平衡法的汽车转向系非线性动力 学 特 性 研究［J］.振动 工 程 学 报，2010，23（3）：355—360. Wang Wei，Song Yuling，Li Guixian.Nonlinear dynamics of automobile steering using incremental harmonic balance method［J］.Journal of Vibration Engineering，2010，23（3）：355—360.

［6］ Kawamura S，Naito T，Zahid H M，et al.Analysis of nonlinear steady state vibration of a multi-degree-offreedom system using component mode synthesis method［J］.Applied Acoustics，2008，69（7）：624—633.

［7］ Bonello P，Minh Hai P.A receptance harmonic balance technique for the computation of the vibration of a whole aero-engine model with nonlinear bearings［J］.Journal of Sound and Vibration，2009，324（1）：221—242.

［8］ Cigeroglu E，Ozguven H N.Nonlinear vibration analysis of bladed disks with dry friction dampers［J］. Journal of Sound and Vibration，2006，295（3-5）：1 028—1 043.

［9］ Guskov M，Sinou J，Thouverez F.Multi-dimensional harmonic balance applied to rotor dynamics［J］.Mechanics Research Communications，2008，35（8）：537—545.

［10］Didier J，Sinou J，Faverjon B.Nonlinear vibrations of a mechanical system with non-regular nonlinearities and uncertainties［J］.Communications in Nonlinear Science and Numerical Simulation，2013，18（11）：3 250—3 270.

［11］Cameron T M，Griffin J H.An alternating frequency/ time domain method for calculating the steady-state response of nonlinear dynamic systems［J］.ASME Journal of Applied Mechanics，1989，56（1）：149—154.

［12］Zhou J X，Zhang L.Incremental harmonic balance method for predicting amplitudes of a multi-dof nonlinear wheel shimmy system with combined Coulomb and quadratic damping［J］.Journal of Sound and Vibration，2005，279（1）：403—416.

［13］Lau S L，Zhang W S.Nonlinear vibrations of piecewise linear systems by incremental harmonic balance method［J］.ASME Journal of Applied Mechanics，1992，59（1）：153—160.

［14］Raghothama A，Narayanan S.Bifurcation and chaos in geared rotor bearing system by incremental harmonic balance method［J］.Journal of Sound and Vibration，1999，226（3）：469—492.

The demension-reductive incremental harmonic balance method for solving the response of local nonlinear dynamic system with multi-frequency excitation

YAO Hong-liang，WANG Chong-yang，WANG Fan，WEN Bang-chun
（School of Mechanical Engineering and Automation，Northeastern University，Shenyang China，110819）

The time-consuming problem of the traditional incremental harmonic balance method can be often encountered when determining the steady state response of local nonlinear system with multi-frequency excitation.The demension-reductive incremental harmonic balance method is proposed to increase the efficiency.For this method，the quantitative comparison relationship of every harmonic in each degree of freedom in the response of the local nonlinear system is analyzed by using the harmonic balance theory，and the dimension of system is reduced to be equal to the dimension of the nonlinear structure by using this relationship.Then the multi-frequency incremental harmonic balance method for the reduced complex system，as well as the response restoring method of every harmonic in each degree of freedom of the original system is deduced.Furthermore，the accuracy and efficiency of the proposed method is verified by using the dual-frequency excitation local nonlinear cantilever beam system.Studies show that the accuracy of the proposed method is in line with the traditional method but the efficiency is much higher with the system of less degree of freedom.

local nonlinearity；steady state response；multi-frequency excitation；incremental harmonic balance method；dimension reduction

O322

：A

1004-4523（2015）05-0741-07

10.16385/j.cnki.issn.1004-4523.2015.05.008

姚红良（1979—），男，博士，副教授。电话：15998389686；E-mail：hlyao@mail.neu.edu.cn

2014-04-15；

2014-11-10

国家重点基础研究发展计划资助项目（2011CB706504）；中央高校基本科研业务费专项资金资助项目（N120403007）；国家自然科学基金资助项目（51005042）