论文部分内容阅读
摘要： 针对传统增量谐波平衡法求解多频激励局部非线性系统周期响应耗时太长的问题，提出了降维增量谐波平衡方法。首先通过谐波平衡理论分析了多频激励局部非线性系统响应中各自由度各次谐波的定量对比关系，并且根据该定量对比关系使系统的维数降至与非线性自由度个数相同；其次针对降维后的复数非线性系统推导了多频增量谐波平衡法，以及原系统各自由度各阶响应的还原方法；最后利用双频激励局部非线性悬臂梁系统进行了所提方法的精度和效率验证。结果表明：该方法的精度与传统方法一致，但是在局部非线性自由度较少时其效率远高于传统方法。
关键词： 局部非线性； 周期响应； 多频激励； 增量谐波平衡法； 降维
中图分类号： O322 文献标志码： A 文章编号： 1004-4523（2015）05-0741-07
引 言
工程中具有局部非线性特性的系统或者结构很多，如与外界发生接触的梁或桁架结构、非线性油膜轴承支撑的转子系统等。分析局部非线性系统的稳态响应，对于掌握这些系统的动态特性具有重要的意义。
目前对多自由度系统进行响应分析，所采用的方法主要有两类：时域分析方法如RungeKutta方法、Newmarkβ法等；频域分析方法如谐波平衡法（Harmonic Balance， HB）[1]、描述函数法（Describing Function， DF）[2]、增量谐波平衡法（Incremental Harmonic Balance， IHB）[35]等。在分析相同振动问题时，频域分析方法的速度要比时域分析方法快很多，因为这些方法能够跳过瞬态振动的分析时间而直接进入稳态；同时，频域分析方法还能够直接揭示出系统中频率特性和变化情况。
但是，在分析自由度较多的系统时，频域方法也会遇到计算耗时太多的问题，因此必须考虑采取降维措施。目前有两种应用较多的降维措施：第一种是采用模态降维的方法，应用最普遍的是固定界面模态综合法[6]，即将结构划分为线性和非线性子结构两部分，对其中的线性子结构应用模态截断，仅保留低阶模态，从而降低系统的自由度；第二种方法是利用局部非线性系统的特性——将原系统分为线性子结构和非线性子结构，用非线性子结构中各自由度的运动来描述线性子结构中各自由度的运动，从而将自由度降低。目前第二种方法已经应用于谐波平衡法[78]和描述函数法[2]。
目前各种降维措施主要针对单频激励局部非线性系统。实际工程中存在大量的多频激励系统[910]，如果采用传统的求解方法进行求解，计算耗时很长，因此必须提出降维的多频激励系统频域分析方法。目前仅有文献[10]提出了一种基于谐波平衡法的多频激励降维分析方法，基于其他频域方法的多频激励降维分析方法尚未见报道。
从算法原理上，增量谐波平衡法比谐波平衡法对强非线性问题、大范围变化参数问题的处理能力更强[1112]。因此，本文提出了基于增量谐波平衡法的多频激励降维分析方法。首先，介绍了传统多频激励增量谐波平衡方法及其不足；然后，利用局部非线性系统特性进行系统降维；再后，结合降维方法改进了传统的增量谐波平衡法，提出了多频激励降维增量谐波平衡法；最后，以双频激励局部非线性悬臂梁系统为例对该方法的精度和效率进行了验证。
采用有限元方法建立悬臂梁的动力学模型，梁单元采用铁木辛柯模型。节点i的广义坐标为xiθyi，其中xi为竖直方向位移，θyi为转角。
采用增量谐波平衡法求解系统响应时，首先需要设定响应中频率成分，假设各次谐波取到激励频率的r倍，即ωk=l1ω1+l2ω2， l1，l2≤r
（23） 例如当激振频率ω1和ω2分别为100和140 rad/s，r=4时，将ωk中的非正值频率和重复频率去掉，可以得到l1和l2的分布如图2所示，共有37个频率分量。
采用本文方法不仅可以求出非线性发生位置处的响应，还可以很方便地求出其他位置响应中的各次频率成分，如图4所示为第5点、第11点和第20点的频域响应。由图4可以看出，响应中含有多种不可公约的组合频率。
通过该方法可以直接求出激励频率变化时，系统响应中各次谐波分量的变化。例如取ω1=100～400 rad/s，ω2=1.4ω1变化时，第5点和第20点的响应三维谱图分别如图5（a）和（b）所示。这是时域方法所不具有的一个优势。
保持边界条件、激励条件不变，假设系统有两处发生碰撞，分别是在节点5和节点11，如图6所示。系统接触刚度为kn1=kn2=5×105 N/m，间隙分别为e1=0.001 m和e2=0.01 m。
分别采用传统多频增量谐波平衡法和本文方法进行响应计算，求得悬臂梁中第5点、第11点和第20点响应分别如图7（a），（b）和（c）所示。
比较在两个非线性自由度时本文方法与Newmark方法的精度，得到取不同谐波次数时各点精度对比如表2所示。从表中可以看出，随着所选取谐波项次数的增加，本文方法的精度逐渐提高。
4.3 效率分析
比较传统的多频增量谐波平衡方法和本文方法的效率。表3是进行单自由度非线性时两种方法效率的比较，表4是进行二自由度非线性时两种方法效率的比较。本文所采用的计算机性能：CPU2.83G，内存3.25G，采用编程语言Matlab。
从表中可以看出，在本例中本文方法在非线性个数较少、所取谐波项少时效率很高；当非线性自由度较多，且所取谐波项较多时，方法的优势逐渐减弱。这是因为随着谐波项次数增多，形成的瞬态刚度矩阵的维数也增大，造成效率下降。
但是在整体自由度增多而非线性自由度保持不变时，本文方法的瞬态刚度矩阵维数保持不变，因此效率不变；而传统IHB方法的瞬态刚度矩阵维数增大，效率下降。此时，本文方法的优越性就可以体现出来。
5 结 论 （1） 针对多频激励的局部非线性系统响应求解耗时长的问题，提出了降维的增量谐波平衡法，在整体自由度较多而非线性自由度较少时，本文方法具有很高的效率；
（2） 该方法由于利用了多自由度系统的局部非线性特性和稳态响应特性，将各自由度的各次谐波响应用局部非线性自由度的响应表示，从而能够使系统维数降低，且降维后系统含有原系统所有动力学特性；
（3） 该方法的理论基础为稳态响应各次谐波的谐波平衡理论，因此仅能对稳态响应进行分析，而不能分析混沌运动等非稳态响应。
参考文献：
[1] Liu J K， Chen F X， Chen Y M. Bifurcation analysis of aeroelastic 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. Multiharmonic 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 timevarying 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 multidegreeoffreedom 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 aeroengine 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（35）：1 028—1 043.
[9] Guskov M， Sinou J， Thouverez F. Multidimensional 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 nonregular 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 steadystate 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 multidof 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.
Abstract： The timeconsuming problem of the traditional incremental harmonic balance method can be often encountered when determining the steady state response of local nonlinear system with multifrequency excitation. The demensionreductive 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 multifrequency 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 dualfrequency 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.
Key words： local nonlinearity； steady state response； multifrequency excitation； incremental harmonic balance method； dimension reduction
作者简介： 姚红良（1979—），男，博士，副教授。电话： 15998389686； Email： hlyao@mail.neu.edu.cn
关键词： 局部非线性； 周期响应； 多频激励； 增量谐波平衡法； 降维
中图分类号： O322 文献标志码： A 文章编号： 1004-4523（2015）05-0741-07
引 言
工程中具有局部非线性特性的系统或者结构很多，如与外界发生接触的梁或桁架结构、非线性油膜轴承支撑的转子系统等。分析局部非线性系统的稳态响应，对于掌握这些系统的动态特性具有重要的意义。
目前对多自由度系统进行响应分析，所采用的方法主要有两类：时域分析方法如RungeKutta方法、Newmarkβ法等；频域分析方法如谐波平衡法（Harmonic Balance， HB）[1]、描述函数法（Describing Function， DF）[2]、增量谐波平衡法（Incremental Harmonic Balance， IHB）[35]等。在分析相同振动问题时，频域分析方法的速度要比时域分析方法快很多，因为这些方法能够跳过瞬态振动的分析时间而直接进入稳态；同时，频域分析方法还能够直接揭示出系统中频率特性和变化情况。
但是，在分析自由度较多的系统时，频域方法也会遇到计算耗时太多的问题，因此必须考虑采取降维措施。目前有两种应用较多的降维措施：第一种是采用模态降维的方法，应用最普遍的是固定界面模态综合法[6]，即将结构划分为线性和非线性子结构两部分，对其中的线性子结构应用模态截断，仅保留低阶模态，从而降低系统的自由度；第二种方法是利用局部非线性系统的特性——将原系统分为线性子结构和非线性子结构，用非线性子结构中各自由度的运动来描述线性子结构中各自由度的运动，从而将自由度降低。目前第二种方法已经应用于谐波平衡法[78]和描述函数法[2]。
目前各种降维措施主要针对单频激励局部非线性系统。实际工程中存在大量的多频激励系统[910]，如果采用传统的求解方法进行求解，计算耗时很长，因此必须提出降维的多频激励系统频域分析方法。目前仅有文献[10]提出了一种基于谐波平衡法的多频激励降维分析方法，基于其他频域方法的多频激励降维分析方法尚未见报道。
从算法原理上，增量谐波平衡法比谐波平衡法对强非线性问题、大范围变化参数问题的处理能力更强[1112]。因此，本文提出了基于增量谐波平衡法的多频激励降维分析方法。首先，介绍了传统多频激励增量谐波平衡方法及其不足；然后，利用局部非线性系统特性进行系统降维；再后，结合降维方法改进了传统的增量谐波平衡法，提出了多频激励降维增量谐波平衡法；最后，以双频激励局部非线性悬臂梁系统为例对该方法的精度和效率进行了验证。
采用有限元方法建立悬臂梁的动力学模型，梁单元采用铁木辛柯模型。节点i的广义坐标为xiθyi，其中xi为竖直方向位移，θyi为转角。
采用增量谐波平衡法求解系统响应时，首先需要设定响应中频率成分，假设各次谐波取到激励频率的r倍，即ωk=l1ω1+l2ω2， l1，l2≤r
（23） 例如当激振频率ω1和ω2分别为100和140 rad/s，r=4时，将ωk中的非正值频率和重复频率去掉，可以得到l1和l2的分布如图2所示，共有37个频率分量。
采用本文方法不仅可以求出非线性发生位置处的响应，还可以很方便地求出其他位置响应中的各次频率成分，如图4所示为第5点、第11点和第20点的频域响应。由图4可以看出，响应中含有多种不可公约的组合频率。
通过该方法可以直接求出激励频率变化时，系统响应中各次谐波分量的变化。例如取ω1=100～400 rad/s，ω2=1.4ω1变化时，第5点和第20点的响应三维谱图分别如图5（a）和（b）所示。这是时域方法所不具有的一个优势。
保持边界条件、激励条件不变，假设系统有两处发生碰撞，分别是在节点5和节点11，如图6所示。系统接触刚度为kn1=kn2=5×105 N/m，间隙分别为e1=0.001 m和e2=0.01 m。
分别采用传统多频增量谐波平衡法和本文方法进行响应计算，求得悬臂梁中第5点、第11点和第20点响应分别如图7（a），（b）和（c）所示。
比较在两个非线性自由度时本文方法与Newmark方法的精度，得到取不同谐波次数时各点精度对比如表2所示。从表中可以看出，随着所选取谐波项次数的增加，本文方法的精度逐渐提高。
4.3 效率分析
比较传统的多频增量谐波平衡方法和本文方法的效率。表3是进行单自由度非线性时两种方法效率的比较，表4是进行二自由度非线性时两种方法效率的比较。本文所采用的计算机性能：CPU2.83G，内存3.25G，采用编程语言Matlab。
从表中可以看出，在本例中本文方法在非线性个数较少、所取谐波项少时效率很高；当非线性自由度较多，且所取谐波项较多时，方法的优势逐渐减弱。这是因为随着谐波项次数增多，形成的瞬态刚度矩阵的维数也增大，造成效率下降。
但是在整体自由度增多而非线性自由度保持不变时，本文方法的瞬态刚度矩阵维数保持不变，因此效率不变；而传统IHB方法的瞬态刚度矩阵维数增大，效率下降。此时，本文方法的优越性就可以体现出来。
5 结 论 （1） 针对多频激励的局部非线性系统响应求解耗时长的问题，提出了降维的增量谐波平衡法，在整体自由度较多而非线性自由度较少时，本文方法具有很高的效率；
（2） 该方法由于利用了多自由度系统的局部非线性特性和稳态响应特性，将各自由度的各次谐波响应用局部非线性自由度的响应表示，从而能够使系统维数降低，且降维后系统含有原系统所有动力学特性；
（3） 该方法的理论基础为稳态响应各次谐波的谐波平衡理论，因此仅能对稳态响应进行分析，而不能分析混沌运动等非稳态响应。
参考文献：
[1] Liu J K， Chen F X， Chen Y M. Bifurcation analysis of aeroelastic 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. Multiharmonic 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 timevarying 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 multidegreeoffreedom 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 aeroengine 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（35）：1 028—1 043.
[9] Guskov M， Sinou J， Thouverez F. Multidimensional 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 nonregular 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 steadystate 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 multidof 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.
Abstract： The timeconsuming problem of the traditional incremental harmonic balance method can be often encountered when determining the steady state response of local nonlinear system with multifrequency excitation. The demensionreductive 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 multifrequency 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 dualfrequency 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.
Key words： local nonlinearity； steady state response； multifrequency excitation； incremental harmonic balance method； dimension reduction
作者简介： 姚红良（1979—），男，博士，副教授。电话： 15998389686； Email： hlyao@mail.neu.edu.cn