多頻激勵(lì)局部非線性系統(tǒng)響應(yīng)求解的降維增量諧波平衡法

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

（東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽110819）

多頻激勵(lì)局部非線性系統(tǒng)響應(yīng)求解的降維增量諧波平衡法

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

（東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽110819）

針對(duì)傳統(tǒng)增量諧波平衡法求解多頻激勵(lì)局部非線性系統(tǒng)周期響應(yīng)耗時(shí)太長的問題，提出了降維增量諧波平衡方法。首先通過諧波平衡理論分析了多頻激勵(lì)局部非線性系統(tǒng)響應(yīng)中各自由度各次諧波的定量對(duì)比關(guān)系，并且根據(jù)該定量對(duì)比關(guān)系使系統(tǒng)的維數(shù)降至與非線性自由度個(gè)數(shù)相同；其次針對(duì)降維后的復(fù)數(shù)非線性系統(tǒng)推導(dǎo)了多頻增量諧波平衡法，以及原系統(tǒng)各自由度各階響應(yīng)的還原方法；最后利用雙頻激勵(lì)局部非線性懸臂梁系統(tǒng)進(jìn)行了所提方法的精度和效率驗(yàn)證。結(jié)果表明：該方法的精度與傳統(tǒng)方法一致，但是在局部非線性自由度較少時(shí)其效率遠(yuǎn)高于傳統(tǒng)方法。

局部非線性；周期響應(yīng)；多頻激勵(lì)；增量諧波平衡法；降維

引 言

工程中具有局部非線性特性的系統(tǒng)或者結(jié)構(gòu)很多，如與外界發(fā)生接觸的梁或桁架結(jié)構(gòu)、非線性油膜軸承支撐的轉(zhuǎn)子系統(tǒng)等。分析局部非線性系統(tǒng)的穩(wěn)態(tài)響應(yīng)，對(duì)于掌握這些系統(tǒng)的動(dòng)態(tài)特性具有重要的意義。

目前對(duì)多自由度系統(tǒng)進(jìn)行響應(yīng)分析，所采用的方法主要有兩類：時(shí)域分析方法如Runge-Kutta方法、Newmark-β法等；頻域分析方法如諧波平衡法（Harmonic Balance，HB）［1］、描述函數(shù)法（Describing Function，DF）［2］、增量諧波平衡法（Incremental Harmonic Balance，IHB）［3-5］等。在分析相同振動(dòng)問題時(shí)，頻域分析方法的速度要比時(shí)域分析方法快很多，因?yàn)檫@些方法能夠跳過瞬態(tài)振動(dòng)的分析時(shí)間而直接進(jìn)入穩(wěn)態(tài)；同時(shí)，頻域分析方法還能夠直接揭示出系統(tǒng)中頻率特性和變化情況。

但是，在分析自由度較多的系統(tǒng)時(shí)，頻域方法也會(huì)遇到計(jì)算耗時(shí)太多的問題，因此必須考慮采取降維措施。目前有兩種應(yīng)用較多的降維措施：第一種是采用模態(tài)降維的方法，應(yīng)用最普遍的是固定界面模態(tài)綜合法［6］，即將結(jié)構(gòu)劃分為線性和非線性子結(jié)構(gòu)兩部分，對(duì)其中的線性子結(jié)構(gòu)應(yīng)用模態(tài)截?cái)?，僅保留低階模態(tài)，從而降低系統(tǒng)的自由度；第二種方法是利用局部非線性系統(tǒng)的特性——將原系統(tǒng)分為線性子結(jié)構(gòu)和非線性子結(jié)構(gòu)，用非線性子結(jié)構(gòu)中各自由度的運(yùn)動(dòng)來描述線性子結(jié)構(gòu)中各自由度的運(yùn)動(dòng)，從而將自由度降低。目前第二種方法已經(jīng)應(yīng)用于諧波平衡法［7-8］和描述函數(shù)法［2］。

目前各種降維措施主要針對(duì)單頻激勵(lì)局部非線性系統(tǒng)。實(shí)際工程中存在大量的多頻激勵(lì)系統(tǒng)［9-10］，如果采用傳統(tǒng)的求解方法進(jìn)行求解，計(jì)算耗時(shí)很長，因此必須提出降維的多頻激勵(lì)系統(tǒng)頻域分析方法。目前僅有文獻(xiàn)［10］提出了一種基于諧波平衡法的多頻激勵(lì)降維分析方法，基于其他頻域方法的多頻激勵(lì)降維分析方法尚未見報(bào)道。

從算法原理上，增量諧波平衡法比諧波平衡法對(duì)強(qiáng)非線性問題、大范圍變化參數(shù)問題的處理能力更強(qiáng)［11-12］。因此，本文提 出了基于 增量諧 波平衡法的多頻激勵(lì)降維分析方法。首先，介紹了傳統(tǒng)多頻激勵(lì)增量諧波平衡方法及其不足；然后，利用局部非線性系統(tǒng)特性進(jìn)行系統(tǒng)降維；再后，結(jié)合降維方法改進(jìn)了傳統(tǒng)的增量諧波平衡法，提出了多頻激勵(lì)降維增量諧波平衡法；最后，以雙頻激勵(lì)局部非線性懸臂梁系統(tǒng)為例對(duì)該方法的精度和效率進(jìn)行了驗(yàn)證。

1 局部非線性多自由度系統(tǒng)運(yùn)動(dòng)微分方程及增量諧波平衡法求解

1.1 局部非線性多自由度系統(tǒng)的運(yùn)動(dòng)微分方程

工程中常見的N自由度非線性系統(tǒng)，受不同頻率激勵(lì)，其動(dòng)力學(xué)方程形式為

式中M，C，K為質(zhì)量、阻尼、剛度矩陣；fnx，˙（）

x為非線性項(xiàng)；ωi為各激振頻率；Fci和Fsi為對(duì)應(yīng)于ωi的激勵(lì)幅值。

當(dāng)非線性系統(tǒng)處于穩(wěn)態(tài)時(shí)，振動(dòng)響應(yīng)x中自由度k的響應(yīng)可以展開為多重傅里葉級(jí)數(shù)形式

設(shè)ω= ［ω1，…，ωrT］，li=［l1，…，lr］，νi=li·，則式（2）可以簡化為

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

對(duì)于形如式（1）的系數(shù)為實(shí)數(shù)的方程，對(duì)x取增量εx，得到

將式（4）代入式（1），并忽略二次以上高次項(xiàng)，得

取τ=［ω1t，…，ωrt］T=［τ1，…，τr］T，對(duì)式（5）兩邊取Galerkin過程，有

即

對(duì)式（7）進(jìn)行迭代，可以求得A0，從而求得x。

如果系統(tǒng)自由度為N個(gè)，諧波取到r次，則形成的瞬態(tài)剛度矩陣Km的維數(shù)為 ［N×（2r+1）］× ［N×（2r+1）］，因此在N比較小時(shí)，該方法的效率不成問題，但是當(dāng)N比較大時(shí)，如N=20，r=25，則矩陣維數(shù)為1 020×1 020，這導(dǎo)致計(jì)算速度下降很快。這是傳統(tǒng)增量諧波平衡法的不足。

2 多頻激勵(lì)局部非線性系統(tǒng)的降維

當(dāng)非線性系統(tǒng)處于穩(wěn)態(tài)時(shí)，非線性力fn（x，˙x）也是穩(wěn)態(tài)的，可以展開為

重新排列方程，將與非線性相關(guān)的自由度排列到前邊，即

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

由式（13）得

因此可以得到

由式（12）和（15）可得

式（16）的實(shí)部為

式（17）可由增量諧波平衡法求解，求出ˉx1后結(jié)合式（11）求出2。

3 降維后系統(tǒng)的增量諧波平衡法求解

對(duì)形如式（17）的復(fù)數(shù)系數(shù)方程，設(shè)

仿照式（6）可以得到：

對(duì)式（21）進(jìn)行迭代，可以求得A0，從而求得x1，再結(jié)合式（11）可以求得x2，這樣就得到了原系統(tǒng)各自由度的響應(yīng)。

該方法可以稱為降維的增量諧波平衡法，該方法迭代公式中Km的維數(shù)比傳統(tǒng)多維增量諧波平衡法大為減少。

4 數(shù)值試驗(yàn)

4.1 數(shù)學(xué)模型

如圖1所示的懸臂梁系統(tǒng)：設(shè)梁截面為圓形，半徑為5 mm；梁總長為400 mm，分成20段；梁材料密度為7 850 kg/m3，彈性模量為2.1×1011Pa；α= 0和β=1×10-3為阻尼；頻率不同的簡諧激振力分別作用于第14節(jié)點(diǎn)和17節(jié)點(diǎn)，幅值分別為50和30 N；假設(shè)在第J1=11節(jié)點(diǎn)有限位裝置，梁與限位裝置作用時(shí)受力符合下式

式中e為間隙，kn為接觸剛度。

圖1 懸臂梁模型圖Fig.1 Cantilever beam with single impact

采用有限元方法建立懸臂梁的動(dòng)力學(xué)模型，梁單元采用鐵木辛柯模型。節(jié)點(diǎn)i的廣義坐標(biāo)為，其中xi為豎直方向位移，θyi為轉(zhuǎn)角。

采用增量諧波平衡法求解系統(tǒng)響應(yīng)時(shí)，首先需要設(shè)定響應(yīng)中頻率成分，假設(shè)各次諧波取到激勵(lì)頻率的r倍，即

例如當(dāng)激振頻率ω1和ω2分別為100和140rad/s，r=4時(shí)，將ωk中的非正值頻率和重復(fù)頻率去掉，可以得到l1和l2的分布如圖2所示，共有37個(gè)頻率分量。

圖2 響應(yīng)中頻率選取Fig.2 Thechosenfrequencyintheresponse

4.2 精度分析

4.2.1 單自由度非線性精度分析

取激振頻率ω1=250rad/s，ω2=1.4ω1，kn= 2×105N/m，e=0.01m。分別采用Newmark-β方法和本文方法進(jìn)行響應(yīng)計(jì)算。采用本文方法時(shí)，所取得諧波次數(shù)分別按r=2，3，4來選擇。計(jì)算中每個(gè)周期分為50段。本文中非線性力屬于分段線性非線性力，在采用增量諧波平衡法進(jìn)行積分時(shí)，可以采用對(duì)分插值方法，如文獻(xiàn)［13-14］所示。

求得懸臂梁中第5點(diǎn)、第11點(diǎn)和第20點(diǎn)響應(yīng)分別如圖3（a），（b）和（c）所示。從圖中可以看出，即使所取諧波項(xiàng)次數(shù)較少，本文方法也有很高的精度。

式中 δ為本文方法的綜合誤差，x（tk）和（tk）分別為用Newmark方法和本文方法計(jì)算得到的tk時(shí)間的響應(yīng)，m為所取的響應(yīng)總數(shù)。

采用該方法，得到取不同諧波次數(shù)時(shí)各點(diǎn)精度，列于表1。從表中可以看出，隨著所選取諧波項(xiàng)次數(shù)的增加，本文方法的精度逐漸提高。

表1 不同諧波項(xiàng)時(shí)本文方法精度誤差Tab.1 Precisionerrorofthepresentmethodunderdifferent harmonicterms

采用本文方法不僅可以求出非線性發(fā)生位置處的響應(yīng)，還可以很方便地求出其他位置響應(yīng)中的各次頻率成分，如圖4所示為第5點(diǎn)、第11點(diǎn)和第20點(diǎn)的頻域響應(yīng)。由圖4可以看出，響應(yīng)中含有多種不可公約的組合頻率。

通過該方法可以直接求出激勵(lì)頻率變化時(shí)，系統(tǒng)響應(yīng)中各次諧波分量的變化。例如取ω1=100～400rad/s，ω2=1.4ω1變化時(shí)，第5點(diǎn)和第20點(diǎn)的響應(yīng)三維譜圖分別如圖5（a）和（b）所示。這是時(shí)域方法所不具有的一個(gè)優(yōu)勢。

圖3 ω1=250rad/s時(shí)各節(jié)點(diǎn)時(shí)域響應(yīng)Fig.3 Thetime-domainresponseofsomenodeswhenω1=250rad/s

圖4 ω1=250 rad/s時(shí)各節(jié)點(diǎn)頻域響應(yīng)Fig.4 The frequency-domain response of some nodes whenω1=250 rad/s

圖5 系統(tǒng)響應(yīng)的三維譜圖Fig.5 The 3D spectrogram of the system responses

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

保持邊界條件、激勵(lì)條件不變，假設(shè)系統(tǒng)有兩處發(fā)生碰撞，分別是在節(jié)點(diǎn)5和節(jié)點(diǎn)11，如圖6所示。系統(tǒng)接觸剛度為kn1=kn2=5×105N/m，間隙分別為e1=0.001 m和e2=0.01 m。

圖6 雙碰撞位置懸臂梁模型圖Fig.6 Cantilever beam with double impact

分別采用傳統(tǒng)多頻增量諧波平衡法和本文方法進(jìn)行響應(yīng)計(jì)算，求得懸臂梁中第5點(diǎn)、第11點(diǎn)和第20點(diǎn)響應(yīng)分別如圖7（a），（b）和（c）所示。

比較在兩個(gè)非線性自由度時(shí)本文方法與Newmark方法的精度，得到取不同諧波次數(shù)時(shí)各點(diǎn)精度對(duì)比如表2所示。從表中可以看出，隨著所選取諧波項(xiàng)次數(shù)的增加，本文方法的精度逐漸提高。

圖7 ω1=250 rad/s時(shí)各節(jié)點(diǎn)時(shí)域響應(yīng)Fig.7 The time-domain response of some nodes whenω1=250 rad/s

表2 不同諧波項(xiàng)時(shí)本文方法精度誤差Tab.2 Precision error of the present method under different harmonic terms

4.3 效率分析

比較傳統(tǒng)的多頻增量諧波平衡方法和本文方法的效率。表3是進(jìn)行單自由度非線性時(shí)兩種方法效率的比較，表4是進(jìn)行二自由度非線性時(shí)兩種方法效率 的比 較。本 文所 采用 的 計(jì) 算 機(jī) 性 能：CPU 2.83G，內(nèi)存3.25G，采用編程語言Matlab。

表3 單自由度非線性時(shí)運(yùn)行所需CPU時(shí)間對(duì)比Tab.3 CPU time for the two methods with 1 nonlinearity location

表4 兩個(gè)自由度非線性時(shí)運(yùn)行所需CPU時(shí)間對(duì)比Tab.4 CPU time for the two methods with 2 nonlinearity location

從表中可以看出，在本例中本文方法在非線性個(gè)數(shù)較少、所取諧波項(xiàng)少時(shí)效率很高；當(dāng)非線性自由度較多，且所取諧波項(xiàng)較多時(shí)，方法的優(yōu)勢逐漸減弱。這是因?yàn)殡S著諧波項(xiàng)次數(shù)增多，形成的瞬態(tài)剛度矩陣的維數(shù)也增大，造成效率下降。

但是在整體自由度增多而非線性自由度保持不變時(shí)，本文方法的瞬態(tài)剛度矩陣維數(shù)保持不變，因此效率不變；而傳統(tǒng)IHB方法的瞬態(tài)剛度矩陣維數(shù)增大，效率下降。此時(shí)，本文方法的優(yōu)越性就可以體現(xiàn)出來。

5 結(jié) 論

（1）針對(duì)多頻激勵(lì)的局部非線性系統(tǒng)響應(yīng)求解耗時(shí)長的問題，提出了降維的增量諧波平衡法，在整體自由度較多而非線性自由度較少時(shí)，本文方法具有很高的效率；

（2）該方法由于利用了多自由度系統(tǒng)的局部非線性特性和穩(wěn)態(tài)響應(yīng)特性，將各自由度的各次諧波響應(yīng)用局部非線性自由度的響應(yīng)表示，從而能夠使系統(tǒng)維數(shù)降低，且降維后系統(tǒng)含有原系統(tǒng)所有動(dòng)力學(xué)特性；

（3）該方法的理論基礎(chǔ)為穩(wěn)態(tài)響應(yīng)各次諧波的諧波平衡理論，因此僅能對(duì)穩(wěn)態(tài)響應(yīng)進(jìn)行分析，而不能分析混沌運(yùn)動(dòng)等非穩(wěn)態(tài)響應(yīng)。

［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］ 晏致濤，張海峰，李正良.基于增量諧波平衡法的覆冰輸電線舞動(dòng)分析［J］.振動(dòng)工程學(xué)報(bào)，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］ 王威，宋玉玲，李瑰賢.基于增量諧波平衡法的汽車轉(zhuǎn)向系非線性動(dòng)力 學(xué) 特 性 研究［J］.振動(dòng) 工 程 學(xué) 報(bào)，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，F(xiàn)averjon 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

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃資助項(xiàng)目（2011CB706504）；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目（N120403007）；國家自然科學(xué)基金資助項(xiàng)目（51005042）