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

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

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

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

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

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

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

1 基于描述函數的逆陣更新方法

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

(1)

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

1.1 非線性內力的描述函數表示

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

(2)

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

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

(3)

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

(4)

式中：Xk為節點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]中所示的非線性控制中的描述函數，若已知節點k與節點j之間的非線性內力，則可由式(9a)和式(9b)計算描述函數Dkj(X,ω)。由式(8)不難發現，描述函數Dkj(X,ω)是節點k與節點j之間，一個周期內的平均非線性內力與相對位移之比，相當于等效線性剛度

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

(10)

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

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

(11)

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

(12)

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

對于結構中第k個自由度存在如圖1所示的單面限位非線性約束的系統，忽略摩擦及碰撞接觸時的動量和能量損失，非線性內力nkk由質量塊與限位彈簧接觸前后剛度系數發生突變的非光滑因素而導致，可簡化為圖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)的描述函數為

此即等效非線性內力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)

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

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

Δ(X)=δ·D·δT

(18)

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

(19)

非線性位置向量

(20)

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

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

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

(21)

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

(22)

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

X(ω)=θ(X)·F

(23)

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

1.3 非線性頻響的迭代求解

式(23)所示為關于非線性頻響X(ω)的非線性代數方程組，將通過構造迭代格式進行求解。為了改善迭代收斂的穩定性，引入加權參數λ，構造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，迭代收斂的差異性較大，如當激勵頻率接近系統固有頻率時較難收斂。為此，采用一種自適應的λ取值算法[21]：如果在頻率ωi時迭代不收斂，取λ(j)(ωi)=λ(j)(ωi)·r(其中0

2 單面限位單自由度系統的仿真研究

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

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

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

將圖4(c)所示懸臂梁簡化為單自由度系統，接地端的間隙限位非線性約束模型如圖1所示。結構參數為：質量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=0.08 N，f為掃頻頻率，掃頻范圍9～16 Hz，掃頻間隔0.02 Hz。分別用IMUM和RK45進行由低到高和由高到低兩個方向的掃頻分析。

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

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

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

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

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

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

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

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

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

圖7 雙穩態現象Fig.7 Bistable phenomenon

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

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

表1 懸臂梁相關參數

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

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

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

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

仿真分析的工況參數與掃頻實驗一致：①激勵力幅值F=0.2 N和0.4 N時的懸臂梁(未接觸到片梁)和間隙寬度δ=7 mm時約束梁的頻響特征；②間隙寬度δ=7 mm，F依次為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，對應實驗中墊放鋼尺片數分別為4，3，2，1和0。

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

3.1 非線性阻尼識別

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

由工況①中F=0.2 N和0.4 N時懸臂梁實驗測得的頻響曲線，可以大致識別出線性比例阻尼系數α=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，體現出較好的線性性。

工況①間隙寬度為7 mm，對應于實驗中墊放3片鋼尺，施加外激勵力幅值為F=0.8 N時，由約束梁實驗測得的頻響曲線，可以大致識別出非線性比例阻尼系數υ=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的約束梁非線性阻尼系數識別結果Fig.10 Identification results of the nonlinear proportional damping for the cantilever beam with F=0.8 N and δ=7 mm

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

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

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

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

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

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

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

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

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

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

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

選取工況②條件下，即間隙位置固定為7 mm(對應實驗中墊放3片鋼尺)，F依次為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時，位移曲線出現“跳躍”現象，從低到高掃頻與從高到低掃頻兩條曲線不重合，出現跳躍區間，并且隨著激勵力幅值的增大，“跳躍”現象越明顯，跳躍區間越大；

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

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

4 結 論

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