第一作者趙昱男，博士生，1988年12月生

通信作者王時龍男，博士，教授，1966年8月生

多股簧系統(tǒng)動態(tài)響應(yīng)等效線性化分析

趙昱1,2，王時龍1,2，周杰1,2，孫守利1,2，李川3

(1.重慶大學(xué)機(jī)械傳動國家重點(diǎn)實驗室，重慶400044；2.重慶大學(xué)機(jī)械工程學(xué)院，重慶400044；3. 重慶工商大學(xué)檢測控制集成系統(tǒng)重慶巿工程實驗室，重慶400044)

摘要：分析歸一化Bouc-Wen模型描述的滯遲阻尼能量損耗，建立該模型極限環(huán)在任意變形幅值下能量損耗計算公式及快速計算方法，提出多股簧系統(tǒng)頻率響應(yīng)特性的等效線性化分析方法并用數(shù)值仿真進(jìn)行驗證。結(jié)果表明，等效線性化分析與數(shù)值仿真結(jié)果一致，且計算速度遠(yuǎn)高于數(shù)值仿真。該方法可顯著提高多股簧系統(tǒng)設(shè)計效率，對工程中大量具有近似線性系統(tǒng)響應(yīng)特性的多股簧系統(tǒng)設(shè)計有實際意義。

關(guān)鍵詞：多股簧；Bouc-Wen模型；能量損耗；等效線性化

基金項目：國家自然科學(xué)基金項目(51375508，51375517)；教育部創(chuàng)新團(tuán)隊計劃(IRT1196)；重慶高校創(chuàng)新團(tuán)隊項目(KJTD201313)

收稿日期：2014-07-14修改稿收到日期：2014-09-25

中圖分類號:TH122

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.20.009

Abstract:The energy dissipation of hysteretic damping described by the normalized Bouc-Wen model was studied. An equation for calculating the energy dissipation under arbitrary deformation amplitude was established and a fast evaluation approach was provided as well. On this basis, an equivalent linearization analysis method for the dynamic response of systems with stranded wire helical springs was proposed. Numerical simulations were carried out to verify the proposed method. The results obtained by the proposed method coincide well with those obtained by the numerical method while the proposed method is much more efficient than the numerical method. The equivalent linearization method is able to improve the efficiency of the designing of systems with stranded wire helical springs significantly; therefore, it is of practical value for the designing of many practical systems with approximate linear response characteristics.

Equivalent linearization analysis of the dynamic response of systems with stranded wire helical springs

ZHAOYu1,2,WANGShi-long1,2,ZHOUJie1,2,SUNShou-li1,2,LIChuan3(1. State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing 400044, China；2. College of Mechanical Engineering, Chongqing University, Chongqing 400044, China；3. Engineering Laboratory for Detection, Control and Integrated System, Chongqing Technology and Business University, Chongqing 400044, China)

Key words:stranded wire helical spring; Bouc-Wen model; energy dissipation; equivalent linearization

圖1　多股簧 Fig.1 A stranded wire helical spring

多股簧為由多層多股鋼絲擰制的鋼絲繩卷繞而成的螺旋彈簧[1]，見圖1。由于減振效果優(yōu)良、服役壽命長，常用于自動武器等關(guān)鍵部件[2]。多股簧加工、檢測已實現(xiàn)數(shù)控及自動化[3-4]，并在更多領(lǐng)域受到重視。

多股簧結(jié)構(gòu)復(fù)雜，受載變形時彈簧內(nèi)部鋼絲間相互擠壓并相對滑動，導(dǎo)致呈一定非線性剛度及滯遲阻尼，其動態(tài)響應(yīng)常用一種修正的Bouc-Wen模型描述[2]。Bouc-Wen模型[5-6]用于描述磁流變裝置[7]、鋼絲繩減振器[8]等滯遲阻尼器件，具有顯著的非線性表達(dá)式，使其響應(yīng)分析較困難。Okuizumi等[9]提出的基于多尺度法分析方法只能處理特殊的Bouc-Wen模型且計算過程復(fù)雜；Wong等[10]提出的多諧波分析法涉及大量迭代計算易遇收斂問題且分析耗時較長。工程中設(shè)計多股簧系統(tǒng)時仍采用數(shù)值方法，該方法耗時較長，因而降低設(shè)計分析效率，不利于多股簧推廣應(yīng)用。

等效線性化方法作為非線性系統(tǒng)分析的常用方法分析速度快，且對非線性不很強(qiáng)系統(tǒng)往往能給出滿足工程精度要求結(jié)果。大量觀測表明，常見的多股簧系統(tǒng)中通常不會出現(xiàn)超諧波、次諧波響應(yīng)、分岔等典型非線性響應(yīng)現(xiàn)象，符合等效線性化分析應(yīng)用條件。本文通過分析歸一化Bouc-Wen模型的能量損耗，以單自由度系統(tǒng)為例研究多股簧系統(tǒng)等效線性化分析方法，旨為多股簧系統(tǒng)設(shè)計提供高效分析手段。

1多股簧系統(tǒng)模型

1.1多股簧恢復(fù)力-變形量模型

多股簧軸向受載變形時的恢復(fù)力與變形量之關(guān)系即為多股簧恢復(fù)力-變形量模型。多股簧在周期載荷下的恢復(fù)力-變形量關(guān)系見圖2。

圖2　典型多股簧恢復(fù)力-變形量示意圖 Fig.2 Illustration of the typical restoring force-deformation relation of stranded wire helical spring

多股簧動態(tài)響應(yīng)可用修正的Bouc-Wen模型描述，即

(1)

式中：“·”為對時間t求導(dǎo)；r為總恢復(fù)力；-1<ω<1為由歸一化Bouc-Wen模型描述的滯遲力；Δ為多股簧變形量；t為時間；kEi為剛度系數(shù)；kAi為滯遲阻尼放大系數(shù)；N，M為多項式階數(shù)；σ>1/2，ρ>0，κω>0，n>0均為歸一化Bouc-Wen模型參數(shù)。該諸參數(shù)均可通過文獻(xiàn)[2]參數(shù)識別方法的實驗數(shù)據(jù)識別獲得。

1.2單自由度多股簧系統(tǒng)運(yùn)動方程

圖3　單自由度系統(tǒng)模型 Fig.3 Structure of the single degree of freedom system

由質(zhì)量塊m、單股彈簧ks、線性阻尼c0及多股簧組成的單自由度微型懸掛裝置系統(tǒng)結(jié)構(gòu)見圖3。系統(tǒng)中單股簧用作傳感器，通過測量其載荷即可求得質(zhì)量塊位移。令ks=0，則系統(tǒng)看完簡化為簡單多股簧懸掛系統(tǒng)。

圖3系統(tǒng)中單股彈簧上端固定，多股簧下端與可周期振動基礎(chǔ)相連，質(zhì)量塊與兩彈簧自由端相連。其中yg為基礎(chǔ)激勵位移，y為質(zhì)量塊m與基礎(chǔ)的相對位移。考慮重力及系統(tǒng)靜止時兩彈簧預(yù)壓，設(shè)系統(tǒng)具有零初始條件y(0)=0，ω(0)=0，用式(1)建立系統(tǒng)運(yùn)動方程為

(2)

式中：yg為幅值為Ag頻率為f的簡諧激勵，即

yg(t)=Agsin(2πft)

(3)

為便于分析計算，需對式(2)進(jìn)行無量綱化處理。引入空間尺度L、時間尺度T，定義無量綱量位移為x，無量綱時間為τ，令有、無量綱量關(guān)系為

(4)

用式(2)～式(4)得無量綱運(yùn)動方程為

(5)

式中：αEi為剛度系數(shù)；αAi為滯遲阻尼放大系數(shù)；cd為阻尼系數(shù)；Ω為頻率；F為激勵幅值；即

(6)

αAi=L1-iT2κωkAi/m

(7)

cd=Tc0/m

(8)

Ω=2πfT

(9)

F=T2Ag(ks-Ω2m)/(Lm)

(10)

2歸一化Bouc-Wen模型的能量損耗分析

文獻(xiàn)[11]基于高斯超幾何函數(shù)研究Bouc-Wen模型解析解。用高斯超幾何函數(shù)可建立歸一化Bouc-Wen模型極限環(huán)變形幅值與能量損耗關(guān)系，為等效線性化分析基礎(chǔ)。

2.1歸一化Bouc-Wen模型極限環(huán)

式(5)歸一化Bouc-Wen模型參數(shù)滿足σ>1/2時其極限環(huán)的有界輸入、輸出(BIBO)穩(wěn)定，且輸出z的極限為±1[12]。x取不同幅值時生成的2個極限環(huán)ABCD見圖4。極限環(huán)A、C點(diǎn)為x達(dá)到極值點(diǎn)，B、D點(diǎn)為z=0點(diǎn)。虛線AECF為以A、C為對角頂點(diǎn)的假想矩形。

圖4　歸一化Bouc-Wen模型極限環(huán) Fig.4 Limit cycles of the normalized Bouc-Wen model

式(5)無量綱歸一化Bouc-Wen模型可寫為

(11)

x′、z同號時由式(11)可得AD、BC段解析式，即

(12)

x′、z異號時可得CD、AB段解析式，即

(13)

式中：2F1為高斯超幾何函數(shù)[13]；C1，C2為待定常數(shù)。

(14)

利用Bouc-Wen模型對稱性得D點(diǎn)橫坐標(biāo)xD為

(15)

以曲線AD為分析對象，由式(12)中對應(yīng)曲線AD的常數(shù)C1為

C1=xD

(16)

(17)

將式(15)、(16)代入式(17)并化簡，得

(18)

式中：β為不完全函數(shù)。

2.2歸一化Bouc-Wen模型能量損耗

歸一化Bouc-Wen模型一個周期內(nèi)的能量損耗E即為圖4中極限環(huán)所圍面積，計算式為

E=2(SADx-SABx)

(19)

式中：SADx為曲線AD與x軸所圍面積；SABx為曲線AB與x軸所圍陰影面積，即

(20)

據(jù)式(11)，將式(20)改寫為

(21)

利用超幾何函數(shù)寫出式(21)中積分解析式為

(22)

同理可得SADx為

(23)

將式(22)、(23)代入式(19)求得

(24)

2.3能量損耗快速計算

E=ET+4(X-XT)

(25)

為進(jìn)一步加快等效線性化分析速度，可在[0,XT]區(qū)間內(nèi)均勻產(chǎn)生一系列X值，利用并行算法求出其對應(yīng)的E并保存為數(shù)據(jù)表，分析中需要計算某X值對應(yīng)的E值時即可利用插值法快速得到結(jié)果。

取Bouc-Wen模型參數(shù)為σ=2、ρ=1、n=1/2，利用本方法及數(shù)值仿真所得E-X曲線見圖5。仿真時先通過Runge-Kutta法由式(5)得到極限環(huán)，再用“鞋帶公式”[14]計算極限環(huán)面積。由圖5可知，本文方法分析結(jié)果與數(shù)值仿真結(jié)果吻合較好。

圖5　歸一化Bouc-Wen模型能量損耗 Fig.5 Energy dissipation of the normalized Bouc-Wen model

3系統(tǒng)頻率響應(yīng)特性等效線性化分析

等效線性化分析為用簡單線性系統(tǒng)近似描述復(fù)雜的非線性系統(tǒng)。設(shè)式(5)的等效線性化形式為

x″+(cd+ce)x′+kex+Fsin(Ωτ)=0

(26)

式中：“′”為對無量綱時間τ求導(dǎo)；ke為等效剛度系數(shù)；ce為多股簧滯遲阻尼等效阻尼系數(shù)。

由式(2)知，多股簧剛度特性由其變形量即y的多項式表示。因此，質(zhì)量塊在原點(diǎn)附近小幅振動時可將非線性剛度用其一階近似值代替。系統(tǒng)受零均值簡諧激勵時等效線性系統(tǒng)響應(yīng)亦為零均值。取多股簧響應(yīng)模型中非線性彈性項一次項系數(shù)作為等效剛度，即

ke=αE1

(27)

同理，多股簧恢復(fù)力-變形量模型式(5)中滯遲阻尼非線性放大多項式在小振幅時亦可用一階形式近似

(28)

此時，歸一化Bouc-Wen模型極限環(huán)被一次函數(shù)調(diào)制，形狀見圖6。

圖6　經(jīng)調(diào)制的極限環(huán) Fig.6 Modulated limit cycles

由于歸一化Bouc-Wen模型產(chǎn)生的極限環(huán)對稱，經(jīng)一次函數(shù)調(diào)制后極限環(huán)面積不會改變，因此調(diào)制后極限環(huán)消耗能量Es為

Es=αA0E

(29)

等效阻尼系數(shù)ce的取值應(yīng)使等效線性化系統(tǒng)能量消耗即等效阻尼力做功Es與滯遲阻尼消耗能量Ec相等。等效線性系統(tǒng)在簡諧激勵下的穩(wěn)態(tài)相對位移表達(dá)式為

x=Xsin(Ωτ-φ)

(30)

式中：Ω為激勵圓頻率；φ為響應(yīng)與激勵相位差。

阻尼力Fc為

Fc=ceXΩcos(Ωτ-φ)

(31)

阻尼力做功即系統(tǒng)能量損耗Ec為

(32)

系統(tǒng)等效阻尼系數(shù)ce可表示為

(33)

式(26)等效線性系統(tǒng)穩(wěn)態(tài)響應(yīng)幅頻特性為

(34)

相頻特性為

(35)

聯(lián)立式(18)、(24)(或(25))、(34)可得關(guān)于X、E的代數(shù)方程組。利用Levenberg-Marqardt等非線性方程數(shù)值求解算法[15]可方便求出X、E數(shù)值。

工程中一般更關(guān)心絕對位移幅頻特性，只需將激勵與相對位移響應(yīng)合成即可。據(jù)式(3)、(4)、(30)可寫出穩(wěn)態(tài)響應(yīng)的絕對位移xa為

xa=XΩsin(Ωτ-φ)+(Agsin(Ωτ))/L

(36)

其振幅Xa為

(37)

由分析過程知，等效線性化分析方法將復(fù)雜非線性微分方程求解轉(zhuǎn)化成簡單的代數(shù)方程求解，無需時程分析即可求出系統(tǒng)的頻響特性，因而可提高分析效率；但等效線性化分析中的能量損耗及振幅計算等過程始終針對穩(wěn)態(tài)響應(yīng)進(jìn)行。因此，該方法只適用系統(tǒng)的頻響特性分析，不適用過渡狀態(tài)或時程分析。

4仿真驗證

數(shù)值仿真為檢驗等效線性化分析有效性途徑。用數(shù)值仿真求系統(tǒng)幅頻特性過程為：①用Runge-Kutta法求系統(tǒng)在某一激勵頻率下若干周期的時程響應(yīng)；②求出最后一個周期的響應(yīng)振幅及均值，記錄系統(tǒng)最終狀態(tài)；③以最終狀態(tài)作為初始狀態(tài)再求若干周期的時程響應(yīng)；④計算最后一個周期的響應(yīng)振幅及均值，記錄系統(tǒng)最終狀態(tài)；⑤計算第2、4步記錄的振幅及均值間相對誤差，若誤差小于要求值則改變激勵頻率，重復(fù)①～⑤；若誤差大于要求值則重復(fù)③～⑤。本文數(shù)值仿真及等效線性化分析程序均用Mathematica編寫。

取系統(tǒng)參數(shù)m=15 kg，ks=5 000 N/m，c0=50 Ns/m，kE2=105N/m2，kE1=5 000 N/m2，kA1=40/m，kA0=1/m，κω=7 N，n=0.5，σ=2，ρ=1 000。頻率范圍1～15 Hz，在此范圍內(nèi)共分析40個頻率點(diǎn)的系統(tǒng)響應(yīng)。

分別用等效線性化及數(shù)值仿真法分析系統(tǒng)響應(yīng)，所得系統(tǒng)在激勵幅值A(chǔ)g為2 mm、4 mm、6 mm時相對位移幅頻特性見圖7。由圖7可知，相對位移隨激勵頻率增加先逐漸減小在一定頻率下相對位移接近0后，響應(yīng)振幅逐漸增大到共振，共振后又逐漸減小并趨于與激勵振幅相同。

圖7　相對位移幅頻特性 Fig.7 Magnitude response of the relative displacement

圖8　相對位移相頻特性 Fig.8 Phase response of the relative displacement

利用等效線性化方法直接可獲得系統(tǒng)相頻特性。Ag=4 mm時相對位移相頻特性見圖8。其它激勵幅值下相頻特性與圖8類似。由圖8可知，等效線性系統(tǒng)的相頻特性與一般線性系統(tǒng)類似，響應(yīng)相位差隨頻率增加而增加，共振時相位差達(dá)90°；激勵頻率很高時相位差接近180°，此時絕對位移接近0。

工程中更關(guān)心系統(tǒng)的絕對位移響應(yīng)。在激勵幅值A(chǔ)g為不同值時的絕對位移幅頻特性見圖9。由圖9可知，等效線性化與數(shù)值仿真分析結(jié)果一致。激勵幅值較小時等效線性化分析振幅偏小，隨激勵振幅增加誤差逐漸減?。患钫穹^續(xù)增加時等效線性化振幅較數(shù)值仿真結(jié)果大，誤差逐漸增加。因振幅增大后多股簧恢復(fù)力-變形量模型中非線性因素影響逐漸增大，使等效線性系統(tǒng)對非線性系統(tǒng)近似精度降低。

圖9　絕對位移幅頻特性 Fig.9 Magnitude response of the absolute displacement

等效線性化分析的最大優(yōu)勢為計算耗時少效率高。數(shù)值仿真與等效線性化分析耗時對比見表1。由表1可見，等效線性化耗時遠(yuǎn)小于數(shù)值仿真，由此表明，等效線性化分析為能顯著提高分析效率方法。實際工程中設(shè)計含多股簧系統(tǒng)時可將其作為初步設(shè)計分析方法，設(shè)計系統(tǒng)各參數(shù)大致值，再利用數(shù)值仿真等校核，進(jìn)一步優(yōu)化設(shè)計參數(shù)，提高設(shè)計效率。

表1　分析時間

5結(jié)論

在研究歸一化Bouc-Wen模型能量損耗基礎(chǔ)上提出多股簧系統(tǒng)等效線性化分析方法，結(jié)論如下：

(1)歸一化Bouc-Wen模型極限環(huán)在任意變形幅值下的能量損耗可通過高斯超幾何函數(shù)精確表達(dá)，并可由求解關(guān)于輸出幅值及能量損耗的代數(shù)方程求出。等效線性化分析前可利用并行計算快速建立能量損耗數(shù)據(jù)表，分析過程中利用插值法可快速獲得特定振幅下的能量損耗。

(2)等效線性化分析能避免非線性微分方程求解，而通過代數(shù)方程求解系統(tǒng)的頻率響應(yīng)特性，顯著提高分析效率。

(3)等效線性化分析誤差與激勵幅值有關(guān)，在一定范圍內(nèi)誤差先隨激勵幅值增加而減小，而后隨激勵振幅增加而增大。

(4)等效線性化方法適用分析工程中大量具有近似線性系統(tǒng)響應(yīng)特性的多股簧系統(tǒng)頻率響應(yīng)特性，不適用系統(tǒng)時程分析。

參考文獻(xiàn)

[1]王時龍,雷松,周杰,等. 兩端并圈多股彈簧的沖擊響應(yīng)研究[J]. 振動與沖擊,2011,30(3):64-68.

WANG Shi-long，LEI Song，ZHOU Jie,et al. Impact response of stranded wires helical springs with closed ends[J]. Journal of Vibration and Shock,2011,30(3):64-68.

[2]Zhao Y，Wang S L，Zhou J，et al. Modeling and identification of the dynamic behavior of stranded wire helical springs[J]. Journal of Vibroengineering,2013, 15(1): 326-339.

[3]Peng Y X，Wang S L，Zhou J，et al. Structural design, numerical simulation and control system of a machine tool for stranded wire helical springs[J]. Journal of Manufacturing Systems, 2012, 31(1SI): 34-41.

[4]王時龍,張明明,周杰,等. 振動狀態(tài)下螺旋彈簧運(yùn)動狀態(tài)模型[J]. 機(jī)械工程學(xué)報,2012,48(1):78-83.

WANG Shi-long，ZHANG Ming-ming，ZHOU Jie，et al. Motion model of helical springs under vibrational condition[J]. Journal of Mechanical Engineering,2012,48(1):78-83.

[5]Bouc R. Forced vibration of mechanical systems with hysteresis[C]// Djadkov S. Proceedings of the 4th Conference on Nonlinear Oscillations. Prague：Academia,1967: 315.

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

[7]相恒波,方秦,王瑋,等. 磁流變阻尼器抗爆隔震性能的數(shù)值模擬[J]. 振動與沖擊,2007,26(11):106-111.

XIANG Heng-bo，F(xiàn)ANG Qin，WANG Wei,et al. Numerical simulation on shock isolation of blast-resistant structure with magnetorheological dampers[J]. Journal of Vibration and Shock,2007,26(11):106-111.

[8]Ni Y Q，Ko J M，Wong C W，et al. Modelling and identification of a wire-cable vibration isolator via a cyclic loading test part 1: experiments and model development[J]. Proceedings of the Institution of Mechanical Engineers Part I-Journal of Systems and Control Engineering,1999, 213(I3): 163-171.

[9]Okuizumi N，Kimura K. Multiple time scale analysis of hysteretic systems subjected to harmonic excitation[J]. Journal of Sound and Vibration,2004, 272(3): 675-701.

[10]Wong C W，Ni Y Q，Lau S L. Steady-state oscillation of hysteretic differential model 1. response analysis[J]. Journal of Engineering Mechanics-Asce,1994, 120(11): 2271-2298.

[11]Charalampakis A E，Koumousis V K. On the response and dissipated energy of Bouc-Wen hysteretic model[J]. Journal of Sound and Vibration. 2008, 309(3/4/5): 887-895.

[12]Ikhouane F, Manosa V, Rodellar J. Dynamic properties of the hysteretic Bouc-Wen model[J]. Systems & Control Letters, 2007, 56(3): 197-205.

[13]Olver F W. NIST handbook of mathematical functions[M]. Cambridge: Cambridge University Press, 2010:318.

[14]Braden B. The surveyor’s area formula[J]. The College Mathematics Journal,1986, 17(4): 326-337.

[15]Burden R L，F(xiàn)aires J D. Numerical analysis[M]. Boston: Cengage Learning, 2010.