IHB法在多自由度Bouc-Wen滯回非線性系統響應特性研究中的應用

趙 倩，劉子良，姚紅良，聞邦椿

(1.上海理工大學 機械工程學院，上海 200093；2.東北大學 機械工程與自動化學院，沈陽 110819)

滯回非線性常用于描述位移或速度與恢復力之間、應變與應力之間的滯后關系，此類系統廣泛存在于振動利用、土木工程等領域[1]。滯回非線性一般來自工程材料特別是復合材料的非線性特性、接觸面摩擦特性和結合面接觸變形等。由于滯回振動系統是一種多值性的非解析系統，其動力學理論分析比較困難[2]。特別是近年來，隨著微/納米技術的迅速發展，基于智能材料(壓電陶瓷、磁致伸縮材料、形狀記憶合金)的執行器在航空航天技術、機械系統隔振、地震工程、土木工程等領域得到廣泛應用，然而智能材料的滯回特性作為一種非常規的非線性現象，其多映射、非平滑和率相關性使相關問題的研究變得更加復雜[3]。

以往關于滯回非線性系統的研究可歸結為描述滯回系統的力學模型、滯回模型的參數辨識、滯回系統的求解方法及動特性分析[4]三個方面。到目前為止，已發展了若干力學模型描述滯回系統，而Bouc-Wen模型[5-7]作為一種光滑曲線型滯回模型，其既包含非線性阻尼，又包含非線性剛度，對各種光滑的滯回曲線都能較好地近似描述，同時也能很好地反應智能材料阻尼器[8-11]、壓電致動器[12]、超磁致伸縮致動器[13]等的動態性能，因此，廣泛應用于工程實際當中。

對于滯回系統的動力學分析方法，最先以等價線性化法[14]為主，然而其分析結果難以真實反應系統的非線性本質，且僅適用于分析低自由度系統。隨著計算機技術的發展，數值積分方法(如Newmark-β法，Runge-Kutta法)變成一種有效的分析手段，但對于多自由度非線性系統來說往往耗時太長。然而對由微分方程控制的動力系統，求解其穩態解的數值方法主要還有多尺度法[15]、KBM法[16]、IHB法[17]等。IHB法作為一種重要的頻域分析方法，其優點在于可以跳過時域方法中的瞬態過程，直接進入穩態過程，因此計算效率高，耗時少。它是一種適于求解微分方程動力系統穩態周期解的半解析半數值方法，可以直接給出頻域分布結果，與弧長法[18]結合還可以求得不穩定解。

關于IHB法在滯回非線性系統分析方面的應用，1985年Pierre等[19]通過引入符號函數的技巧，將IHB法成功地應用于具有Den Hartog理想干摩擦環節的非線性振動系統；曲敬鐳將文獻[19]的結果推廣至具有雙線性遲滯特性的非線性系統；白鴻柏等[20]進一步將IHB法推廣應用于三次非線性黏性阻尼雙線性遲滯振動系統；熊蕊等[21]將IHB法用于求解二階遲滯非線性控制系統簡諧激勵下的響應。

然而以往對滯回非線性系統的研究主要針對單自由度或兩自由度系統。考慮到工程實際中多為多自由度系統，本文力圖將IHB法推廣至含Bouc-Wen滯回力的多自由度非線性系統的研究。故針對懸臂梁系統，引入Bouc-Wen滯回非線性模型，提出將滯回力作為一個增加的自由度的方式，重新建立梁結構與滯回力相互作用的運動微分方程；然后將IHB法推廣至求解該類系統的穩態運動周期解，同時引入弧長法解決由滯回非線性引起的跳躍現象；通過與Newmark-β法進行比較，體現了該法在效率上的優越性，并分析了一些滯回特性及頻域響應特性。

1 Bouc-Wen滯回力-懸臂梁系統數學模型

1.1 Bouc-Wen模型及其數學描述

Bouc于1967年介紹了一種至今廣為利用的微分形式的數學模型來模擬滯回非線性系統；1976年Wen等改進了這個模型，并使其系統化地可以表示各種滯回線的滯回非線性系統。其表達式為

(1)

以往對于具有此模型的振動系統常用以下單自由度振動微分方程方程描述

(2)

1.2 Bouc-Wen滯回力-懸臂梁系統模型

懸臂梁結構廣泛應用于橋梁、建筑等工程領域，如懸臂橋、塔吊等。在實際工程分析中，大部分受力部件都可簡化為懸臂梁。故針對如圖1所示的懸臂梁系統進行研究。將Bouc-Wen模型引入到系統中，從而建立Bouc-Wen滯回力-懸臂梁系統的動力學模型。

懸臂梁采用彈性軸段單元建立，忽略軸向變形，每節點包括兩個移動方向和兩個轉動方向的自由度。梁全長300 mm，橫截面直徑10 mm，均勻劃分為6段，外激勵作用在節點6處y方向。系統的微分方程為

(3)

圖1 Bouc-Wen滯回力-懸臂梁模型Fig.1 Model of cantilever with Bouc-Wen hysteresis

假設Bouc-Wen滯回力模型作用于節點5處，即將非線性阻尼力施加于自由度i，則方程(3)可化為

(4)

將z視作變量，并取n=1，將式(4)轉化為如下形式

(5)

其中，

取

(6)

則式(5)可化為

(7)

在方程(7)兩邊同除以系數k，可以簡化得如下形式

(8)

為便于推導，令

(9)

將方程式(8)簡化為

(10)

2 增量諧波平衡法-弧長法簡介

對于式(10)描述的系統，當處于穩態運動時，根據諧波平衡理論，節點i處的周期響應xi可以展開成如下傅里葉級數形式

(11)

式中：Ai={ai0ai1bi1…airbir}T,CS={1 cosωtsinωt…cos(rωt) sin(rωt)}。

x=SA

(12)

令τ=ωt，式(10)化為

(13)

設Δx、Δω為響應x和ω的增量，即可采用增量形式表示為

x=x0+Δx,ω=ω0+Δω

(14)

將式(14)代入式(13)，忽略二次以上高次項，可以得到

(15)

設δx=SδA為x的變分，對式(15)兩邊取Galerkin過程，有

(16)

(17)

即

KmcΔA+RmcAΔω=Rm1A+Rm2

(18)

其中，

(19)

弧長法是Newton-Raphson迭代方法的一種改進迭代方法，與IHB法相結合能夠求得系統的不穩定穩態解。

針對式(18)得到的平衡方程，應用弧長法，取

Ψ(A,ω)=-Rm1A-Rm2

(20)

將式(20)化為

(21)

式中：Kt=Km，q=RmcA。

增加控制方程

ΔATΔA+Δλφ2ω2=Δl

(22)

式中：φ為調整ω與ΔX的一個參數。

通過聯立式(21)和(22)求解ΔA和Δλ，從而可以迭代求得傅里葉系數矩陣A，再由式(12)得到響應x。

3 數值仿真與分析

采用IHB法求解系統的穩態運動周期解(求解過程均取4次諧波項)，得到不同控制參數、位置處的滯回曲線、幅頻響應及三維譜圖，研究系統的滯回特性及頻域響應特性。仿真參數：外激勵角頻率600 rad/s，幅值5 N，系數αf=1×104，A=1，k=1×103。阻尼矩陣采用瑞利阻尼C=aM+bK的形式，其中a=0，b=1×10-4。

3.1 方法的精度和效率比較

為了驗證方法的精確性和時效性，與Newmark-β數值積分方法進行比較，取α=0.3,β=-0.7，時域響應及滯回響應結果見圖2，兩者的結果基本一致，驗證了算法的精度；從計算效率上來看，求解時域響應時，Newmark-β法用時7.223 s，而IHB法僅用時0.834 s，計算效率較數值積分方法提高了8.66倍，體現了該法保證精度的前提下在計算效率方面的優勢。

3.2 滯回特性研究

滯回曲線作為一種重要的滯回特性分析工具，其所圍成的面積為結構在一個周期里所消耗的能量，因此可直觀地反映結構系統的阻尼特性和耗能能力。分別選取①α>0,β>0且α≥β；②α>0,β>0且α<β；③α>0,β<0且|α|≥|β|；④α>0,β<0且|α|<|β|四種不同的情況進行分析，得到不同位置處的滯回曲線，情況①～③分別見于圖3～圖5，第四種情況不再是光滑的滯回曲線，系統將產生不穩定解。

圖2 兩種方法所得結果比較Fig.2 Comparison of results obtained by the two methods

圖3 β>0且α≥β時的滯回曲線Fig.3 Hysteretic curves at β>0 and α≥β

圖4 β>0且α<β時的滯回曲線Fig.4 Hysteretic curves at β>0 and α<β

圖5 β<0且|α|≥|β|時的滯回曲線Fig.5 Hysteretic curves at β<0 and |α|≥|β|

由圖3～圖4可以看出，當β>0時，系統呈現軟特性，即系統恢復力隨位移絕對值的增加而減小；且當α≥β(圖3)時，滯回曲線形狀飽滿，耗能能力強，系統位移相對較小；當α<β(圖4)時，滯回曲線形狀狹窄，耗能能力減弱，系統位移相對較大。由圖5可以看出，當α<β且|α|≥|β|時，系統呈現硬特性，即系統恢復力隨位移絕對值的增加而增加。對比不同位置處的滯回曲線，激勵位置處形狀更為飽滿，耗能能力更強，其它位置(節點3)處形狀與非線性滯回力作用處類似，位移相對較小。

3.3 頻域響應特性研究

采用IHB法對系統進行頻域響應特性研究。分別取β>0和β<0兩種情況，得到不同參數組合下的幅頻響應曲線如圖6所示。同時給出了三維譜圖，見圖7。

(a) β>0

(b) β<0圖6 不同控制參數的幅頻響應曲線Fig.6 Amplitude-frequency curves with different parameters

圖6(a)中，當β>0時，系統呈現軟特性，隨α/β的比值減小，軟特性越明顯，與圖3～圖4的滯回特性分析結果相一致；且α/β的比值越大時，響應的幅值越小，證明此時耗能減振能力更強，對共振幅值的抑制效果增強。

圖6(b)中，當β<0時，系統呈現硬特性，隨α/|β|的比值減小，硬特性越明顯，與圖4的滯回特性分析結果一致；且當α/|β|<1時，系統出現了跳躍現象，對應了滯回特性分析的第四種情況；與圖6(a)(β>0)相同的是，α/|β|的比值越大，響應的幅值越小，系統的耗能減振能力以及對共振幅值的抑制效果越強。

在三維譜圖中，由于滯回力的作用系統中出現了3×倍頻。在圖7(a)～圖(b)兩種情況下，3×倍頻較弱；但圖7(c)～圖(d)中，即當β<0且α/|β|<1時系統出現了明顯的跳躍和幅值放大現象，特別是圖(d)中，產生了更加豐富的2×、4×倍頻成分，當外激勵角速度高于600 rad/s時，3×倍頻超過工頻的振動能量，需予以重視。

(a) β>0(α=0.75, β=0.25)

(b) β<0 (α=0.75, β=-0.25)

(c) β<0 (α=0.25, β=-0.75)

(d) β<0 (α= 0.15, β= -0.85)圖7 不同控制參數的三維譜圖Fig.7 Three-dimensional spectrums with different control parameters

4 結 論

(1) 提出將滯回力作為一個增加的自由度的方式，重新建立梁結構-滯回力相互作用的振動微分方程，將Bouc-Wen模型引入原系統，建立了Bouc-Wen滯回力-懸臂梁系統的數學模型。

(2) 將IHB法推廣至含Bouc-Wen模型的多自由度滯回非線性系統響應的求解，引入弧長法可以解決由遲滯非線性引起的跳躍現象。通過數值仿真表明，該方法能有效地用于此類系統響應特性的研究，并且在保證精度的前提下，效率得到了大幅提高。

(3) 通過本文方法分析了一些滯回特性，當β<0且α/|β|<1時，系統中出現跳躍現象，為不穩定情況，工程中需引起重視并盡量避免；當α/|β|>1時，系統的耗能減振能力以及對共振幅值的抑制效果好，工程中應力圖采用。

參 考 文 獻

[1] 聞邦椿, 李以農, 張義民, 等.振動利用工程[M].北京: 科學出版社, 2005.

[2] 李鴻光, 何旭, 孟光.Bouc-Wen滯回系統動力學特性的仿真研究[J].系統仿真學報, 2004, 16(9): 2009-2011.

LI Hongguang, HE Xu, MENG Guang.Numerical simulation for dynamic characteristics of Bouc-Wen hysteretic system[J].Journal of System Simulation, 2004, 16(9): 2009-2011.

[3] 趙新龍, 汪佳麗.結合誤差變換的Bouc-Wen 遲滯非線性系統反步控制器設計[J].控制理論與應用, 2014, 31(8): 1094-1098.

ZHAO Xinlong, WANG Jiali.Backstepping control with error transformation for Bouc-Wen hysteresis nonlinear system[J].Control Theory and Applications, 2014, 31(8): 1094-1098.

[4] LI H G, MENG G.Nonlinear dynamics of a SDOF oscillator with Bouc-Wen hysteresis[J].Chaos, Solitons and Fractals, 2007, 34: 337-343.

[5] BOUC R.Forced vibration of mechanical system with hysteresis[C].Proceedings of the Fourth Conference on Nonlinear Oscillations, Prague, Czechoslovakia, 1967.

[6] WEN Y K.Method for random vibration of hysteretic systems[J].ACSE of the Engineering Mechanics Division, 1976, 102(2): 249-263.

[7] ISMAIL M, IKHOUANE F, RODELLAR J.The hysteresis Bouc-Wen model, a survey[J].Archives of Computational Methods in Engineering, 2009, 16(2): 161-188.

[8] SPENCER B F, DYKE S J, SAIN M K.Phenomenological model magneto rheological dampers[J].Journal of Engineering Mechanics, 1997, 123(3): 230-238.

[9] BHATTACHARJEE A, CHATTERJEE A.Dissipation in Bouc-Wen model: Small amplitude, large amplitude and two-frequency forcing[J].Journal of Sound and Vibration, 2013, 332: 1807-1819.

[10] KWOK N M, HA Q P, NGUYEN M T, et al.Bouc-Wen model parameter identification for a MR fluid damper using computationally efficient GA[J].ISA Transactions, 2007, 46: 167-179.

[11] TALATAHARI S, KAVEH A, MOHAJER RAHBARI N.Parameter identification of Bouc-Wen model for MR fluid dampers using adaptive charged system search optimization[J].Journal of Mechanical Science and Technology, 2012, 26(8): 2523-2534.

[12] ZHU W, RUI X T.Hysteresis modeling and displacement control of piezoelectric actuators with the frequency- dependent behavior using a generalized Bouc-Wen model[J].Precision Engineering, 2016, 44: 299-307.

[13] GUO Y, MAO J.Modeling of hysteresis nonlinearity based on generalized bouc-wen model for GMA[J].Springer Berlin Heidelberg, 2010, 6424(6): 147-158.

[14] 白鴻柏, 張培林, 黃協清.黏性阻尼滯遲振子簡諧激勵響應的等效線性化計算方法研究[J].振動與沖擊, 2000, 19(4): 44-47.

BAI Hongbai, ZHANG Peilin.HUANG Xieqing.Equivalent linearization approximate method for response computation of a bilinear hysteretic vibration with viscous damping under sinusoidal excitation[J].Journal of Vibration and Shock, 2000, 19(4): 44-47.

[15] 曹樹謙, 付耀東, 賀微波.基于三次遲滯模型的超聲電機圓環定子主共振響應[J].振動與沖擊, 2009, 28(1): 36-41.

CAO Shuqian, FU Yaodong, HE Weibo.Primary resonance of a ring-type stator of ultrasonic motors based on cubic polynomial hysteresis model[J].Journal of Vibration and Shock, 2009, 28(1): 36-41.

[16] 侯東曉, 劉彬, 時培明.一類滯后相對轉動動力學方程的分岔特性及其解析近似解[J].物理學報, 2009, 58(9): 5942-5949.

HOU Dongxiao, LIU Bin, SHI Peiming.The bifurcation of a kind of relative rotational dynamic equation with hysteresis and its approximate solution[J].Acta Physica Sinca, 2009, 58(9): 5942-5949.

[17] ZHAO Qian, YAO Hongliang, XU Qi et al.Prediction method for steady-state response of local rubbing blade-rotor systems[J].Journal of Mechanical Science and Technology, 2015, 29(4): 1537-1545.

[18] FERREIRA J V, SERPA A L.Application of the arc-length method in nonlinear frequency response[J].Journal of Sound and Vibration, 2005, 284：133-149.

[19] PIEREE C, FERRI A A, DOWELL E H.Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method[J].ASME Journal of Applied Mechanics, 1985, 52(4): 959-964.

[20] 白鴻柏, 黃協清.三次非線性粘性阻尼雙線性滯遲振動系統IHB分析方法[J].西安交通大學學報, 1998, 32(10): 35-38.

BAI Hongbai, HUANG Xieqing.Response of forced vibration of bilinear hysteretic oscillator[J].Journal of Xi’an Jiaotong University, 1998, 32(10): 35-38.

[21] 熊蕊, 劉向東, 耿潔, 等.一類二階遲滯非線性控制系統簡諧激勵響應的IHB分析方法[J].振動與沖擊, 2013, 32(13): 99-105.

XIONG Rui, LIU Xiangdong, GENG Jie, et al.IHB method for response computation of a case of second-order hysteretic nonlinear control systems under a sinusoidal excitation[J].Journal of Vibration and Shock, 2013, 32(13): 99-105.