含單面限位局部非線性結構的主頻響應分析

孫 櫻， 萬雨婷， 陳力奮， 郭其威

(1. 復旦大學 航空航天系，上海 200433；2. 上海市空間飛行器機構重點實驗室，上海 201109)

間隙和限位是實際工程結構中常見的、典型的非線性形式，常見的有鉸間隙、齒輪間隙等，如大型網(wǎng)架式可展空間結構中有許多含間隙的運動副，可能導致空間結構展開時的非線性振動、磨損以及結構各部分展開不同步等問題[1]；限位則往往是考慮安全或功能需求而對物體運動過程中的位置進行了約束，如為了防止橋梁上部結構在地震中發(fā)生落梁破壞而在伸縮縫處安裝的限位裝置[2]。大量研究表明，系統(tǒng)中的間隙或限位對結構動力學特性影響顯著[3-5]。

含間隙結構動力學分析主要涉及間隙內碰撞過程的處理和整體結構動力學方程的求解。目前，根據(jù)結構部件相對運動關系的假設不同，可以將間隙內碰撞過程歸結為四類模型[6]：連續(xù)接觸模型、有限元模型、經(jīng)典碰撞模型和接觸變形模型。Hertz基于接觸變形模型，對于能量耗散只考慮變形位移和變形速度的綜合作用，針對碰撞接觸面較小的靜態(tài)接觸問題，提出了將間隙內碰撞過程的動力學效應近似地等效為分段線性的非線性模型?？紤]到限位約束的結構特征，國內外學者一般也將其簡化為分段線性系統(tǒng)，如厲行軍等[7]將隔振裝置中的限位器簡化為帶限位的單自由度單層隔沖系統(tǒng)。

對于系統(tǒng)中含有分段線性等非線性元件的結構動力學方程，線性理論中的模態(tài)分析法或譜分析法不能適用，一般只能采用半解析方法或數(shù)值方法。常見的半解析方法有諧波平衡法、描述函數(shù)法、Volterra級數(shù)法[8]等。Maezawa等[9]首先將諧波平衡法應用于分析單自由度分段線性問題；Shaw等[10]將帶限位器的小車運動簡化為分段線性單自由度系統(tǒng)，運用映射手段獲得了系統(tǒng)的周期解析；Moussi等[11]對帶雙邊限位器的倒立擺系統(tǒng)用諧波平衡法和漸近法得到正弦激勵下的相圖，研究了第一、第二階共振頻率同激勵幅值的關系；吳志強等[12]將彈性碰撞振動系統(tǒng)簡化為對稱的分段線性模型，運用平均法求得幅頻特性與相頻特性，基于約束分岔理論對系統(tǒng)進行穩(wěn)定性分析。數(shù)值方法則能夠精確地提供時域和頻域的信息。盧緒祥等[13]建立了含對稱間隙結構的碰撞振動動力學模型，采用Runge-Kutta法數(shù)值求解，碰撞振動系統(tǒng)隨著頻率比、激振力幅值和間隙的改變，出現(xiàn)了周期、倍周期和混沌運動；Onesmus等[14]研究了不同接觸力模型下的含兩間隙的平面剛體結構的動力學特性；林道福等[15]給出了帶有限位器的浮筏隔振系統(tǒng)在局部坐標系下的剛度矩陣表達，用紐馬克積分法計算結構在基礎沖擊激勵下的運動響應。半解析法的特點在于計算效率快，但如果結構過于復雜，一般難以得到其理論解形式；而數(shù)值解法對各種結構運動均可計算，但耗時較長，對相關參數(shù)的選取較為困難。

綜合國內外學者的研究成果，對系統(tǒng)中含有分段線性非線性元件的結構，目前的研究主要還是局限在單自由度或者少自由度系統(tǒng)，針對一般的有限元模型的分析方法尚不成熟。本文主要對于含單面限位分段線性模型的局部非線性結構的主頻響應分析方法展開研究，采用描述函數(shù)表征單面限位非線性元件的頻域特性，基于迭代的逆陣更新方法快速求解大型多自由度結構的主頻響應，通過與實驗結果的比較，驗證本文所提出分析方法的有效性。

1 基于描述函數(shù)的逆陣更新方法

受周期性諧波激勵的含有限個非線性元件的局部非線性結構振動微分方程為

(1)

式中：N(t)為結構中的非線性內力向量。由于非線性內力的存在，使得式(1)所示的非線性系統(tǒng)頻響特性和激勵方式、激勵幅值有關；系統(tǒng)響應可能包含激勵頻率以外的諧波成分，且振動能量在這些頻率上的分布情況和激勵幅值密切相關。本文采用描述函數(shù)對非線性內力進行頻域特性表征，針對局部非線性結構利用迭代的逆陣更新方法進行非線性頻響的快速計算。

1.1 非線性內力的描述函數(shù)表示

非線性內力向量N(t)中的第k個分量Nk表示節(jié)點k與節(jié)點j之間的非線性內力在第k個自由度上的分量nkj之和，即

(2)

式中：k=1,2,…，N；j=k為一端接地情形。如式(1)所示的非線性系統(tǒng)受周期性諧波Fsin(ωt)激勵時，設其試探解為

x(t)=Xsin(ωt+θ)

(3)

節(jié)點k與j之間的頻域相對位移Ykj為

(4)

式中：Xk為節(jié)點k的頻域位移響應。

將式(3)所示的試探解代入式(1)，可知各非線性內力形如nkj=nkj(sinωt,cosωt)，將其展開成簡諧形式

(5)

含有限個非線性元件的局部非線性結構一般關注其主頻響應，其對應的非線性內力也僅取其基頻項，即

nkj?Akj(X,ω)sin[ωt+φ1(X,ω)]=Akj(X,ω)sinψ(6)

式中：ψ=ωt+φ1，將Akj(X,ω)用傅里葉積分表示

(7)

令

Dkj(X,ω)

則有

(9a)

(9b)

式(8)所定義的Dkj(X,ω)即為文獻[16]中所示的非線性控制中的描述函數(shù)，若已知節(jié)點k與節(jié)點j之間的非線性內力，則可由式(9a)和式(9b)計算描述函數(shù)Dkj(X,ω)。由式(8)不難發(fā)現(xiàn)，描述函數(shù)Dkj(X,ω)是節(jié)點k與節(jié)點j之間，一個周期內的平均非線性內力與相對位移之比，相當于等效線性剛度

為了評價等效線性系統(tǒng)式(6)和原非線性系統(tǒng)式(5)的貼近程度，引入誤差函數(shù)

(10)

當系統(tǒng)式(1)受到Fsin(ωt)={Fk(t)}激勵時，由文獻[17]可知，當且僅當

φfkj nkj(τ)=φfkj nkj,lin(τ)

(11)

成立時的等效線性系統(tǒng)式(6)才能使式(10)所示的誤差E取值最小。式(11)中φ為相對輸入fkj(t)=Fk(t)-Fj(t)與非線性內力nkj(t)的互相關函數(shù)

(12)

不難驗證，系統(tǒng)式(1)受到Fsin(ωt)激勵時，式(5)所示的原非線性內力與式(6)所示的等效非線性內力滿足式(11)，由此表明，式(6)為最優(yōu)準線性化系統(tǒng)。

對于結構中第k個自由度存在如圖1所示的單面限位非線性約束的系統(tǒng)，忽略摩擦及碰撞接觸時的動量和能量損失，非線性內力nkk由質量塊與限位彈簧接觸前后剛度系數(shù)發(fā)生突變的非光滑因素而導致，可簡化為圖2所示的雙線性模型，其中δ為間隙寬度，k1為接觸前剛度，接觸后剛度則為k1+k2，滿足

(13)

圖1 單面限位示意圖Fig.1 Diagram of single-sided limitation system

圖2 單面限位非線性模型Fig.2 Characteristic of nonlinear internal force

由式(9a)和式(9b)，可得單面限位非線性內力nkk(t)的描述函數(shù)為

此即等效非線性內力nkk(t)滿足式(8)的頻域表征。

1.2 逆陣更新方法(Inverse Matrix Update Method, IMUM)

將式(3)代入式(1)，考慮到式(2)和式(6)，則結構的振動微分方程式(1)可以化為

[-ω2M+iωC+K+Δ(X)]·X(ω)=F

(15)

其中非線性內力向量N(t)的頻域幅值向量

G(X,ω)=Δ(X)·X(ω)

(16)

當結構僅含有限個非線性元件時，如在節(jié)點k處存在接地非線性元件，在節(jié)點i和j以及節(jié)點m與n之間存在非線性元件，則Δ(X)矩陣具有如下稀疏特性

上述Δ(X)可以進一步寫成表示非線性元件位置的向量δ以及非線性內力描述函數(shù)的低階矩陣D的乘積

Δ(X)=δ·D·δT

(18)

其中用描述函數(shù)表示的非線性程度矩陣

(19)

非線性位置向量

(20)

第i列 第j列 第k列 第m列 第n列

式中：M為系統(tǒng)中非線性元件的個數(shù)。記α=(-ω2M+iωC+K)-1為結構中線性部分的動柔度矩陣，令

θ(X)=[α-1+Δ(X)]-1

(21)

稱θ(X)為非線性結構的擬動柔度矩陣。利用Sherman-Morrison法則[18]，由式(18)可得

(22)

由此即把局部非線性結構擬動柔度矩陣θ(X)的N階矩陣求逆問題轉化為M階矩陣的求逆。一般來說，局部非線性結構中，非線性元件個數(shù)M遠小于結構的自由度數(shù)N，即M<

X(ω)=θ(X)·F

(23)

利用式(22)能極大地提高式(23)求解的計算效率，稱該方法為逆陣更新方法[19]。

1.3 非線性頻響的迭代求解

式(23)所示為關于非線性頻響X(ω)的非線性代數(shù)方程組，將通過構造迭代格式進行求解。為了改善迭代收斂的穩(wěn)定性，引入加權參數(shù)λ，構造Mann迭代[20]格式

X(j+1)(ω)=λ(j)X(j+1)*(ω)+(1-λ(j))X(j)(ω),λ(j)∈(0,1)

(24)

式中：X(j+1)*(ω)=θ(X(j)(ω))·F具體算法流程見圖3。

圖3 逆陣更新的算法流程Fig.3 Flow chart for the computational algorithm IMUM

實際計算結果表明，λ取值越小迭代格式式(24)越容易收斂，但收斂速度越慢；且對于不同的頻率點ωi，迭代收斂的差異性較大，如當激勵頻率接近系統(tǒng)固有頻率時較難收斂。為此，采用一種自適應的λ取值算法[21]：如果在頻率ωi時迭代不收斂，取λ(j)(ωi)=λ(j)(ωi)·r(其中0

2 單面限位單自由度系統(tǒng)的仿真研究

研究如圖4所示含單面限位的T字型梁結構，系統(tǒng)由試件、實驗平臺、激振器、數(shù)據(jù)控制采集系統(tǒng)、加速度和力傳感器組成。試件為一矩形截面懸臂梁，一端用壓塊固定，靠近另一自由端處有片梁限位；片梁用壓塊在兩端固定；片梁與懸臂梁之間的空隙，構成一個一端接地的單面限位非線性約束。實驗中通過調節(jié)片梁高度控制限位間隙δ的大小。

圖4 含單間隙懸臂梁結構Fig.4 Cantilever beam with single-sided limitation

本節(jié)先通過對懸臂梁進行單自由度建模，采用逆陣更新方法(IMUM)求解懸臂梁自由端的位移頻響，并與龍格—庫塔(RK45)法的計算結果進行比較。

將圖4(c)所示懸臂梁簡化為單自由度系統(tǒng)，接地端的間隙限位非線性約束模型如圖1所示。結構參數(shù)為：質量m=0.05 kg，由靜力學實驗測得分段線性剛度分別為k1=300 N/m，k2=1 200 N/m，取間隙寬度δ=1.4 mm，設阻尼c=0.4 N·s/m；激勵為單頻簡諧激勵P=Fsin(2πft)，F(xiàn)=0.08 N，f為掃頻頻率，掃頻范圍9～16 Hz，掃頻間隔0.02 Hz。分別用IMUM和RK45進行由低到高和由高到低兩個方向的掃頻分析。

采用Runge-Kutta法進行數(shù)值計算時采用以下措施：

(1)模擬掃頻時，在初始頻率點處，初始條件為零位移、零初速度；前一個頻率點(已達到穩(wěn)定)最后時刻的運動狀態(tài)(位移和速度)作為后續(xù)頻率點的初始條件。

(2)在每個頻率點處，截取系統(tǒng)相鄰的幾段時域響應，比較其幅值的大小，以此判斷系統(tǒng)是否已經(jīng)達到穩(wěn)定狀態(tài)。一旦判定系統(tǒng)穩(wěn)定，就進入下一個頻率點。

圖5所示為兩種方法的計算結果。不難發(fā)現(xiàn)，兩種方法得到的正向和反向掃頻位移頻響曲線均表現(xiàn)出明顯的跳躍現(xiàn)象：由低到高跳躍點在14 Hz附近，由高到低跳躍點在13 Hz附近；且都有明顯的非光滑轉折點。

圖6所示為兩種方法計算結果的比較，結果表明了IMUM方法計算單面限位非線性結構頻響特性的有效性。圖中峰值附近略有差異，初步分析原因應在于本文的IMUM方法僅考慮了響應中的主頻成分。

圖5 單自由度單面限位系統(tǒng)位移頻響曲線Fig.5 Displacement frequency response of SDOF structure with single-sided limitation

圖6 IMUM和Runge-Kutta法系統(tǒng)位移頻響比較Fig.6 Displacement frequency response comparison between IMUM and RK4

比較兩種方法的計算時間可以發(fā)現(xiàn)，在相同計算條件下，RK45耗時2 045.35 s，IMUM耗時0.78 s。顯然，在計算精度保證的前提下，IMUM方法極大地提高了非線性頻響分析的計算效率。

進一步分析定頻激勵下的時域響應。選取跳躍點所在區(qū)域的頻率點13.5 Hz進行定頻激勵，計算時間長度370.37 s內的位移響應，期間進行4次沖擊響應的模擬分析，如圖7(b)所示。分別取一段低幅值穩(wěn)態(tài)和高幅值穩(wěn)態(tài)時的時域響應信號(圖7(c)和圖7(d))，觀察到兩個顯著不同的穩(wěn)態(tài)響應幅值1.19 mm和1.80 mm。由于結構位移幅值超過間隙寬度，運動過程受到單面限位的限制，高幅值響應曲線具有上下不對稱性；而結構在反向掃頻的該頻率點處，位移幅值尚未因接近間隙寬度而發(fā)生幅值點的跳躍，表現(xiàn)為低幅值響應曲線上下對稱。

圖7 雙穩(wěn)態(tài)現(xiàn)象Fig.7 Bistable phenomenon

3 含單面限位的懸臂梁多自由度結構頻響分析

實驗采用的是如圖4所示的單面限位懸臂梁結構，相關參數(shù)見表1。懸臂梁與片梁的接觸剛度由重力-位移的靜力學試驗測得，取1 200 N/m。激勵點、響應點和間隙位置距離懸臂梁自由端分別為262 mm，6 mm，30 mm，如圖8(a)所示。

表1 懸臂梁相關參數(shù)

在ANSYS中建立單面限位懸臂梁模型，采用beam188單元，沿長度方向每10 mm劃分一單元，共62個單元；力傳感器、加速度傳感器、片梁位置的單元號

圖8 實驗結構多自由度有限元模型Fig.8 Finite element model of experiment structure

分別為36，62和59，如圖8(b)所示。加速度傳感器質量為10 g，相當于懸臂梁上的一個集中質量，因此結構整體質量和剛度矩陣維數(shù)均為372×372。

通過改變單面限位間隙δ、激勵力幅值F的大小和不同的掃頻方式，研究間隙寬度和激勵力對單面限位非線性約束梁結構頻響特征的影響。

仿真分析的工況參數(shù)與掃頻實驗一致：①激勵力幅值F=0.2 N和0.4 N時的懸臂梁(未接觸到片梁)和間隙寬度δ=7 mm時約束梁的頻響特征；②間隙寬度δ=7 mm，F(xiàn)依次為1.0 N,1.4 N,1.8 N和2.0 N，對應實驗中墊放3片鋼尺；③激勵力幅值F=0.8 N，間隙寬度δ依次為4 mm,7 mm,10 mm,12 mm和20 mm，對應實驗中墊放鋼尺片數(shù)分別為4，3，2，1和0。

每種工況下，進行由低到高和由高到低兩種掃頻模式。

3.1 非線性阻尼識別

考慮單面限位間隙約束懸臂梁中的非線性阻尼，假設線性與非線性阻尼均為跟剛度矩陣有關的比例阻尼：C=αK線性剛度+βM+υK非線性剛度。線性剛度部分對應阻尼不隨工況的改變而改變，非線性剛度部分對應阻尼隨工況的改變而變化。

由工況①中F=0.2 N和0.4 N時懸臂梁實驗測得的頻響曲線，可以大致識別出線性比例阻尼系數(shù)α=2×10-4和β=0.2。圖9為IMUM仿真計算與實驗曲線的對比，其中共振頻率和幅值：F=0.2 N時為15.2 Hz和3.4 mm；F=0.4 N時為15.2 Hz和7.5 mm，體現(xiàn)出較好的線性性。

工況①間隙寬度為7 mm，對應于實驗中墊放3片鋼尺，施加外激勵力幅值為F=0.8 N時，由約束梁實驗測得的頻響曲線，可以大致識別出非線性比例阻尼系數(shù)υ=4×10-7。圖10為IMUM仿真計算與實驗曲線的對比，圖中的頻響曲線特征(共振峰值和幅值的轉折點)顯示出用IMUM法仿真分析的有效性。

(a) (b)圖9 懸臂梁線性比例阻尼識別結果Fig.9 Identification results of the linear proportional damping for the cantilever beam

圖10 F=0.8 N，間隙寬度為7 mm的約束梁非線性阻尼系數(shù)識別結果Fig.10 Identification results of the nonlinear proportional damping for the cantilever beam with F=0.8 N and δ=7 mm

3.2 仿真計算與實驗結果的對比

利用工況①識別出的比例阻尼系數(shù)，針對工況②給出的4種激勵力幅值情況，計算梁結構的位移頻響曲線，IMUM仿真計算與實驗結果對比如圖11所示。針對工況③給出的5種間隙寬度情況，計算梁結構的位移頻響曲線，仿真計算與實驗結果對比如圖12所示。

從圖11可知，隨著激勵力幅值的增大，結構的頻響特征有下述特點：

1)激勵力幅值F>0.6 N時，無論是仿真結果還是實驗結果，結構的一階固有頻率值和共振峰幅值隨著激勵力幅值的增加而增大，結構表現(xiàn)出明顯的硬彈簧特征；

2)懸臂梁與片梁碰撞后，位移幅值曲線在此處形成拐點，顯示出間隙非線性的非光滑特征。

從圖12可知，隨著墊放鋼尺片數(shù)的增加，即間隙寬度的減小，結構的頻響特征表現(xiàn)為：

1)間隙寬度大于等于10 mm，此時片梁不起作用，結構為懸臂梁；

2)間隙寬度小于7.8 mm，即墊放鋼尺片數(shù)大于等于3片時，隨著間隙寬度的減小(鋼尺片數(shù)的增加)，結構的一階固有頻率值增大，共振峰幅值減小，結構表現(xiàn)出硬彈簧特征；

3)懸臂梁與片梁碰撞后，位移幅值曲線也相應在此處形成拐點。

圖11 約束梁從低到高掃頻位移幅值曲線，間隙寬度7 mmFig.11 Displacement frequency responses for the T-shaped beam by forward sweep with δ=7 mm

圖12 約束梁從低到高掃頻位移幅值曲線，F(xiàn)=0.8 NFig.12 Displacement frequency responses for the T-shaped beam by forward sweep with F=0.8 N

選取工況②條件下，即間隙位置固定為7 mm(對應實驗中墊放3片鋼尺)，F(xiàn)依次為0.2 N，0.4 N,0.6 N,1.0 N,1.6 N,2.0 N進行低高和高低兩種方式仿真計算，得到不同激勵力幅值下的位移幅值曲線見圖13。

圖13 激勵力幅值不同時，從低到高和從高到低掃頻的位移幅值曲線Fig.13 Displacement frequency responses for the T-shaped beam by forward and backward sweep with F=0.2 N，0.4 N,0.6 N,1.0 N,1.6 N,2.0 N

由圖13可知：

(1)激勵力幅值F達到0.6 N時，位移曲線出現(xiàn)“跳躍”現(xiàn)象，從低到高掃頻與從高到低掃頻兩條曲線不重合，出現(xiàn)跳躍區(qū)間，并且隨著激勵力幅值的增大，“跳躍”現(xiàn)象越明顯，跳躍區(qū)間越大；

(2)激勵力幅值F=0.2 N和0.4 N時，從低到高掃頻與從高到低掃頻兩條曲線基本重合，結構表現(xiàn)為線性特征。

綜合圖11～圖13不難發(fā)現(xiàn)，在一定范圍內，結構的間隙寬度越小，激勵力幅值越大，主梁結構出現(xiàn)頻率漂移現(xiàn)象，表現(xiàn)出的硬剛度特性越明顯。單面限位的存在，使結構的頻響曲線出現(xiàn)明顯的“跳躍”現(xiàn)象，呈現(xiàn)非線性振動的雙穩(wěn)態(tài)特性。

4 結 論

本文基于描述函數(shù)和逆陣更新方法對含單面限位局部非線性結構進行基頻響應分析。單自由度模型計算中，通過與RK45數(shù)值方法的對比，顯示出IMUM迭代算法的高效性和準確性；在多自由度有限元模型基頻響應分析中，通過非線性阻尼系數(shù)的識別，單面限位約束梁的IMUM仿真分析與實驗結果呈現(xiàn)了較好的一致性。分析結果表明，外激勵幅值的增加和間隙寬度的減小對結構動態(tài)特性的影響有著相似的規(guī)律，都會使約束梁結構出現(xiàn)頻率漂移現(xiàn)象，表現(xiàn)出越來越硬的剛度特性；間隙限位的存在，使結構的頻響曲線出現(xiàn)明顯的“跳躍”現(xiàn)象，呈現(xiàn)非線性振動的雙穩(wěn)態(tài)特性。