多頻激勵局部非線性系統響應求解的降維增量諧波平衡法

姚紅良，王重陽，王 帆，聞邦椿

（東北大學機械工程與自動化學院，遼寧 沈陽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）