陳恒 王揚渝 金江明

(浙江工業(yè)大學機械工程學院, 杭州 310014)

引言

機翼顫振是飛行中常見可能帶來災難性后果的氣動彈性現(xiàn)象,屬于自激振動.機翼的顫振主要是由機翼的結(jié)構非線性和氣動非線性相耦合引起.由于結(jié)構非線性和氣動非線性的存在,振動幅值限制在一定范圍內(nèi),其主要形式為極限環(huán)振動.極限環(huán)振動是氣動彈性線性的重要問題,造成飛行器結(jié)構的危害,如F-16和F/A18戰(zhàn)機在超音速飛行時均會產(chǎn)生極限環(huán)振動[1].對于極限環(huán)顫振,人們在理論與實驗研究方面已經(jīng)做了不少研究.Gilliat等[2]研究了結(jié)構非線性和失速狀態(tài)下氣動非線性的運動,得到極限環(huán)顫振.Jones和Roberts[3]識別了帶雙線性的非線性氣動彈性結(jié)構的極限環(huán)振動.Tang和Dowell[4]研究了帶掛裝的機翼產(chǎn)生間隙非線性結(jié)構的振動響應.趙永輝等[5]進行了大展弦比夾芯翼大攻角顫振分析.崔鵬等研究了新型運輸機的顫振特征,發(fā)現(xiàn)翼梢裝置對機翼顫振的不利影響[6].周秋萍和邱志平研究了外掛連接具有初偏間隙非線性的機翼顫振問題,并對其顫振進行區(qū)間分析得到顫振邊界曲線[7].

機翼結(jié)構往往帶有間隙非線性、雙線性和碰摩等強非線性特征,存在分岔和混沌等動力特征,減振難度顯著增大.目前其減振方法包括主動控制和被動控制形式.主動控制需測量振動狀態(tài),通過控制策略反饋到控制器,實現(xiàn)對振動的抑制[8].對于非光滑振動系統(tǒng),線性剛度吸振器不再適用[9],Vakakis等[10]提出非線性靶能量轉(zhuǎn)移的概念,耦合非線性能量阱(nonlinear energy sink,以下簡稱NES),實現(xiàn)非線性被動控制.McFarland等[11]實驗驗證NES在強非線性系統(tǒng)中的減振效果.Lee等[12]使用NES成功抑制了范德波爾方程的極限環(huán)振動,并應用NES技術對機翼顫振進行抑制,研究了一個二自由度機翼結(jié)構通過NES的減振效果,得到三種不同的抑制機理[13].

另外,研究氣動彈性系統(tǒng)在不同流速下的整體振動特征,常用的方法是計算其分岔圖.對于預測和分析機翼結(jié)構分岔的方法主要有范式法,中心流方法[14],胞映射法[15],諧波平衡法[16]等.目前,數(shù)值法在氣動彈性系統(tǒng)分析和分岔計算中等到廣泛應用和發(fā)展,如G. Dimitriadis[17]用數(shù)值法計算了帶間隙非線性的氣動彈性結(jié)構的分岔,沈華勛[18]用數(shù)值法研究柔性機翼的穩(wěn)定性問題.而對于氣動彈性系統(tǒng)的強非線性,不同系統(tǒng)具有很強的獨特性,難以用一種通用的方法實現(xiàn)分岔的分析.本文帶控制界面的機翼存在采用雙線這種非光滑的剛度,通過耦合NES實現(xiàn)其顫振抑制,通過基于弧長數(shù)值連續(xù)的方法構造分岔圖,驗證其減振效果.

1 二元機翼基于非線性能量阱的減振模型

考慮不可壓縮流中的具有浮沉、俯仰和控制截面轉(zhuǎn)動三個自由度的二元機翼顫振系統(tǒng)模型(如圖1(a)所示,其中NES質(zhì)量m=0).圖l中e是空氣動力焦點到彈性軸的距離;b是半弦長;a是從半弦點到彈性軸的距離;c是從半弦點到控制截面的距離.h,α和β分別為浮沉位移,俯仰位移和控制截面的轉(zhuǎn)動位移.Kh和Kα表示h和α方向的彈性系數(shù);c1和c2表示h和α方向的非線性剛度系數(shù);Kβ為控制截面轉(zhuǎn)動的剛度,用雙線性表示(如圖1(b)所示).U為氣流速度.運動方程為:

(1)

其中M為機翼質(zhì)量,Sα和Sβ分別為機翼主體和控制截面的質(zhì)量不平衡,Iα和Iβ為機翼主體和控制截面的轉(zhuǎn)動慣量.

圖1 帶控制截面二元機翼模型(a)機翼截面; (b)機翼控制截面雙線性扭轉(zhuǎn)剛度Fig. 1 Two-dimensional wing with control surface model(a) Schematic of the typical wing sections with control-surface; (b) Bilinear torsional stiffness of wing control-surface

(2)

當轉(zhuǎn)角位小角度時力矩Mα≈eLα,Mβ≈eLβ.

當振動結(jié)構耦合本征非線性剛度的結(jié)構會形成能量單向不可逆轉(zhuǎn)移,達到減振效果[10],非線性能量阱結(jié)構本質(zhì)特征是要求耦合結(jié)構的剛度為本征非線性,如三次方非線性剛度的NES結(jié)構可以形成非線性能量阱結(jié)構.在二元機翼模型上增加非線性能量阱結(jié)構(單自由度的三次方非線性剛度的NES結(jié)構),實現(xiàn)對機翼顫振的抑制(如圖1(a)所示,其中NES質(zhì)量m≠0).在原來三自由度的基礎上增加一個自由度,形成四自由度系統(tǒng).其運動方程如下:

(3)

方程(3)的無量綱形式如下:

ελ(y′-δα′-v′)+C(y-δα-v)3=0

δελ(δα′+v′-y′)+C(y-δα-v)3=0

εv″+ελ(v′+δα′-y′)+C(v+δα-y)3=0

(4)

其中δ=0.3.

2 減振前后數(shù)值響應的整體特征

為了描述機翼結(jié)構減振前后的整體響應特征,在不同風速下分別數(shù)值模擬系統(tǒng)響應,記錄響應穩(wěn)定后的極值形成峰值-峰值圖.構造峰值-峰值圖的具體步驟為,給定初始條件,風速從[0.8 1.0]范圍內(nèi)200等分,在每個不同風速下數(shù)值模擬系統(tǒng)的響應,在位移風速坐標下記錄τ∈[8000 8500]內(nèi)響應的極值點(通過數(shù)值響應可認定若系統(tǒng)為穩(wěn)定,則當時間τ>8000時響應已趨于穩(wěn)定).考慮初始大擾動和初始小擾動兩種初始值情況構造峰值-峰值圖,小擾動初值為在浮沉自由度y=0.001,y′=0,其他自由度α=β=v=0,α′=β′=v′=0,其數(shù)值響應的峰值-峰值圖如圖2所示(黑色點為減振前響應極值,紅色點表示減振后響應極值).大擾動初值為在浮沉自由度y=0.1,y′=0,其他自由度α=β=v=0,α′=β′=v′=0,其數(shù)值響應的峰值-峰值圖如圖3所示(黑色點為減振前響應極值,藍色點表示減振后響應極值).

減振前系統(tǒng),對于大初值和小初值始擾動,其峰值-峰值圖保持基本相同(圖2和圖3黑色點).初步反應了系統(tǒng)在耦合NES之前的分岔結(jié)構,可以發(fā)現(xiàn)當u*=0.88時系統(tǒng)發(fā)生Hopf分岔.由峰值-峰值圖可以發(fā)現(xiàn)當風速小于0.903時,控制截面振動幅值小于0.3系統(tǒng)響應為光滑振動,如果風速大于0.903系統(tǒng)導致非光滑振動.當風速介于0.903和0.918之間時響應為周期運動形成極限環(huán).當風速介于0.918與0.922之間以及0.923與0.941之間將出現(xiàn)混沌現(xiàn)象.當風速大于0.903除上述區(qū)域外系統(tǒng)都將發(fā)生準周期運動.

圖2 小初始值下減振前后峰值-峰值圖對比Fig. 2 Peak-peak diagram of the sytem with/without NESat small initial conditions

圖3 大初始值下減振前后峰值-峰值圖對比Fig. 3 Peak-peak diagram of the sytem with/without NES atlarge initial conditions

圖2表示小擾動初值情況下減振前后系統(tǒng)的峰值-峰值圖對比.從圖2可以看出當初值為小擾動時,系統(tǒng)減振效果非常顯著,減振后Hopf分岔點為0.892明顯大于減振前的0.88,這表明如果風速小于0.892,增加NES后系統(tǒng)可以完全抑制顫振.如果風速介于0.892和0.908之間,則減振后系統(tǒng)為極限環(huán)運動,相對于減振前其極限環(huán)的幅值大大減少.當風速值大于0.908則系統(tǒng)減振后為準周期振動,但是其響應的幅值遠小于減振前.

圖3表示大擾動初值情況下減振前后系統(tǒng)的峰值-峰值圖對比.從圖3可以看出當初值為大擾動時,Hopf分岔點同樣為0.892大于減振前的0.88,這表明如果風速小于0.892,增加NES后系統(tǒng)可以完成抑制顫振.風速介于0.892和0.919之間,則減振后系統(tǒng)為極限環(huán)運動,相對于減振前其極限環(huán)的幅值降低.當風速值大于0.919則系統(tǒng)減振前后振幅無明顯變化,減振效果不佳.

3 減振前后系統(tǒng)分岔圖的計算和對比

數(shù)值連續(xù)法用于求解帶一個或多個參數(shù)的非線性問題[19],在求解氣動彈性力學問題中已有相應的軟件包如AUTO[20]和MatCont[21],Roberts等[22]利用其求解雙線性問題的解,Dimitriadis求解了間隙非線性機翼振動的分岔問題[18].數(shù)值連續(xù)法主要用于求解如下形式的非線性代數(shù)方程的解:

f(x,λ)=0

(5)

其中f表示非線性函數(shù)向量,x表示未知數(shù)向量,λ表示變量.給定方程在λ0時已知解x0,利用數(shù)值連續(xù)法尋找在變量λ0+Δλ時的解x0+Δx,其中Δλ為微小變化.

在本文中需要求解的是常微分方程,其形式為:

(6)

其中x=[x1(t),…,xn(t)]T,f=[f1,…,fn],λ為變量參數(shù).使用數(shù)值連續(xù)法需要把常微分方程轉(zhuǎn)化為代數(shù)方程.基于系統(tǒng)的周期解分析分岔問題,把方程(5)轉(zhuǎn)化為一系列的非線性代數(shù)方程:

F(x(t),λ)=0

(7)

其中F表示非線性函數(shù).在特定的參數(shù)λ下,給定一個假設的初始值x0,基于牛頓迭代可以得到其精確解x.

(8)

x1=x0+Δx

(9)

其中Δx表示對x0的改進,經(jīng)過多次迭代可以得到精確值.

把常微分方程轉(zhuǎn)化為代數(shù)方程常用的方法為打靶法,這實際上是一種邊界值方法.若方程的解x(t)是周期解,那么數(shù)值連續(xù)法就轉(zhuǎn)化為求解如下的邊界值問題.

x(0)=x(T)

(10)

常微分方程轉(zhuǎn)化為:

F(x,T,λ)=x(0)-x(T)=0

(11)

給定初始值x(0),通過數(shù)值積分得到x(T),計算函數(shù)F的值.對于求解微分方程的響應可以使用龍格庫塔等數(shù)值方法.假設系統(tǒng)參數(shù)為λ0時,響應為x0(t)周期為T0,那么當參數(shù)為λ=λ0+Δλ時,

=-F(x0(0),T0,λ)

(12)

(13)

方程(11)表示n個方程有n+1個未知數(shù).在這個方程中,使用相位固定法,最簡單的就是假設初值中一個元素始終為0,比如x1(0)=0,這樣未知數(shù)就降為n個.

(14)

其中δx是微小量,δxi表示為n×n矩陣δ×I的第i列.當x(0)的第一個元素均為0時,計算雅可比矩陣從i=2開始,結(jié)果矩陣為n×(n-1).

上述方法為基本的數(shù)值連續(xù)法,稱為自然參數(shù)連續(xù).對于非線性問題,基于自然參數(shù)的連續(xù)法求解中其解的分支可能出現(xiàn)折返,交叉甚至消失.出現(xiàn)這類情況會產(chǎn)生奇異點,自然參數(shù)的連續(xù)法將失效.更加有效的方法是采用弧長連續(xù)(arclength continuation).弧長數(shù)值連續(xù)選用參數(shù)s不再采用系統(tǒng)參數(shù)λ,對應的步長采用Δs.當數(shù)值連續(xù)進行到分岔附近則需要減少步長,以免丟失真實的解.當遠離分岔時若采用很小步長,則影響方法的計算效率.故需采用變步長的策略[15].弧長為x(0),T和λ對應s的改變速率.

(15)

把Δλ看成另外一個未知數(shù),可以得到:

(16)

(17)

此時系統(tǒng)為n+1各方程,可以求解n+1個未知數(shù).

(18)

(19)

對于周期解,另外一個重要性質(zhì)是其穩(wěn)定性.對其初始條件施加一個微小擾動,運動偏離原來周期運動則為不穩(wěn)定,反之為穩(wěn)定.對于一個周期為T的運動軌跡x(t,x(0)),其在t=T時的雅可比矩陣為:

(20)

對初始條件施加微擾動Δx(0),對擾動后的解x(T,x(0)+Δx(0))泰勒展開:

Δx(T)=ΦT(x(0))Δx(0)+o(‖Δx(0)‖2)

(21)

其中Δx(T)=x(T,x(0)+Δx(0))-x(T,x(0)),表示經(jīng)過一個周期后運動的偏離量,如果經(jīng)過m個周期.

Δx(mT)=[ΦT(x(0))]mΔx(0)+o(‖Δx(0)‖2)

(22)

周期運動的穩(wěn)定性取決于矩陣ΦT的2n個特征值,即Floquet乘子.若Floquet乘子至少有一個值大于1,則周期運動為不穩(wěn)定,反之則系統(tǒng)穩(wěn)定.

給定初始的一個周期解,基于上述的弧長數(shù)值連續(xù),可以得到系統(tǒng)在不同參數(shù)下周期解,及其穩(wěn)定性分析,構造系統(tǒng)減振前后的分岔圖,如圖4所示.圖6中紅線和綠線表示減振前系統(tǒng)的分岔圖,其中紅線為穩(wěn)定解,綠線為不穩(wěn)定解.藍線及灰色虛線為減振后的分岔圖,藍線為穩(wěn)定解,虛線表示不穩(wěn)定解.圖5為減振后分岔圖和峰值-峰值圖的對比,圖5左邊一欄為下分支與數(shù)值響應的峰值對比,圖5右側(cè)一欄為上分支和數(shù)值響應的對比.由圖5可得當分岔圖為穩(wěn)定分支則其幅值和數(shù)值響應的峰值非常吻合,其不穩(wěn)定分支可出現(xiàn)兩種情形即準周期運動(小擾動初值風速大于0.908)和混沌現(xiàn)象(大擾動初值風速大于0.957).由圖4可以得到增加NES后系統(tǒng)Hopf分岔點出現(xiàn)在μ*=0.892,這和峰值-峰值圖的結(jié)論一致.在Hopf分岔點出現(xiàn)以前表示系統(tǒng)可以完全抑制顫振.當0.892<μ*<0.908時,分岔圖中存在一個穩(wěn)定下分支,當系統(tǒng)響應為小擾動時運動會趨于下分支.另外0.894<μ*<0.90時存在一個穩(wěn)定的上分支,當系統(tǒng)響應為小擾動時響應會趨于上分支.此時減振前后均為極限環(huán)運動,減振措施可以使得極限環(huán)的幅值更小.由于系統(tǒng)雙線性存在,在不穩(wěn)定分支對于初值為小擾動才生準周期運動,對于大擾動初值則出現(xiàn)混沌現(xiàn)象.

圖4 減振前后分岔圖對比Fig. 4 Comparison of Bifurcation diagram for the systems with/without NES

通過峰分岔圖計算對比,結(jié)合峰值-峰值圖,發(fā)現(xiàn)添加NES作為減振措施當系統(tǒng)為小擾動時可以實現(xiàn)很好的顫振抑制效果.對于大擾動的情景μ*<0.892時可以完全抑制顫振,0.892<μ*<0.918時可以降低極限環(huán)振動幅值.而當μ*>0.918則無明顯減振效果抑制失敗,從分岔圖及峰值-峰值圖可以發(fā)現(xiàn)響應幅值與減振前基本保持相同,不會導致減振失敗的負面效應.機翼顫振抑制大致可分為三種情況:(a)風速小于Hopf分岔點出現(xiàn)之前顫振完全抑制,圖6所示為抑制前后的時間序列對比;(b)0.892<μ*<0.908時,減振前后分岔圖均為穩(wěn)定的極限環(huán)運動,振動幅值大幅降低,其時間序列對比如圖7所示.(c)μ*>0.908初始小擾動,減振前后分岔圖均不穩(wěn)定分支.如圖8所示,減振前響應為非光滑的準周期振動,而減振后為光滑小幅值準周期運動.

圖5 分岔圖與峰值-峰值圖對比Fig. 5 Comparison between Peak-Peak diagram and bifurcation diagram

圖6 μ*=0.888,初始小擾動減振前后時間序列對比Fig. 6 Comparison of time series for the systems with/without NES when μ*=0.888 with small initial conditions

圖7 μ*=0.895,初始大擾動減振前后時間序列對比Fig. 7 Comparison of time series for the systems with/without NES when μ*=0.895 with large initial conditions

圖8 μ*=0.935,初始小擾動減振前后時間序列對比Fig. 8 Comparison of time series for the systems with/without NES when μ*=0.935 with small initial conditions

4 結(jié)論

帶控制截面的機翼結(jié)構簡化為三自由度二元剛性氣動彈性模型,機翼結(jié)構耦合非線性能量阱,實現(xiàn)對機翼顫振的抑制.建立機翼減振前后的運動方程,通過數(shù)值模擬構造減振前后峰值-峰值圖,初步反應其整體的減振效果.減振后Hopf分岔點為0.892明顯大于減振前的0.88,這表明如果風速小于0.892,增加NES后系統(tǒng)可以完成抑制顫振.小擾動情況或者大擾動初值在風速低于0.919時,系統(tǒng)減振效果顯著.

通過弧長數(shù)值連續(xù)法和結(jié)合Floquet算子構造減振前后的分岔圖,并研究了其穩(wěn)定性.通過Hopf分岔點的推遲及對比減振前后穩(wěn)定和不穩(wěn)定分支的存在,驗證了三類不同形式的顫振抑制:顫振完全抑制(風速小于Hopf分岔點出現(xiàn)之前),極限環(huán)部分抑制(減振前后分極限環(huán)運動幅值大幅降低),準周期運動抑制(減振前響應為非光滑的準周期振動,而減振后為光滑的小幅值準周期運動).

1Bunton R, Denegri Jr C. Limit cycle oscillation characteristics of fighter aircraft.JournalofAircraft,EngineeringNote, 2000,37:916～918

2Jones D P, Roberts I, Gaitonde A L. Identification of limit cycles for piecewise nonlinear aeroelastic systems.JournalofFluidsandStructures, 2007,23:1012～1028

3Tang D, Dowell E H. Experimental and theoretical study of gust response for a wing-store model with freeplay.JournalofSoundandVibration, 2006,295:659～684

4趙永輝,胡海巖. 大展弦比夾芯翼大攻角顫振分析. 振動工程學報, 2004,17(1):25～30 (Zhao Y H, Hu H Y. Flutter Analysis of a high-aspect-ratio sandwich wing under large angle of attack.JournalofVibrationEngineering, 2004,17(1):25～30 (in Chinese))

5崔鵬,韓景龍. 新型運輸機機翼的顫振特性分析. 振動工程學報, 2011,24(2):192～199 (Cui P, Han J L. Flutter analysis of new transport—type wings.JournalofVibrationEngineering, 2011,24(2):192～199 (in Chinese))

6周秋萍,邱志平. 機翼帶外掛系統(tǒng)極限環(huán)顫振的區(qū)間分析. 航空學報, 2010,31(3):514～518 (Zhou Q P, Qiu Z P. Interval analysis for limit cycle flutter of a wing with an external store.ActaAeronauticaetAstronauticaSinica, 2010,31(3):514～518 (in Chinese))

7Tanaka N. Impact vibration control using a semi-active damper.JournalofSoundandVibration, 1992,158(2):277～292

8Tanaka N. Experiment of the dynamic damper with a preview action.JSMEInternationalJournal, 1987,53(487):650～655

9Lee Y S, Kerschen G, Vakakis A F, et al. Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment.PhysicaD, 2005,204(1):41～69

10 Kerschen G, Vakakis A F, Panagopoulos P N, et al. Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment.PhysicaD, 2005,204(1):41～69

11 Lee Y S, Vakakis A F, Bergman L A, et al. Suppression of limit cycle oscillations in the van der pol oscillator by means of passive nonlinear energy sinks (NESs).StructuralControlandHealthMonitoring, 2006,13:41～75

12 Lee Y S, Vakakis A F, Bergman L A, et al. Suppressing aeroelastic instability using broadband passive targeted energy transfers, Part 1: Theory.AIAAJournal, 2007,45(3):693～711

13 Vio G A, Cooper J E. Limit cycle oscillation prediction for aeroelastic systems with discrete bilinear stiffness.InternationalJournalofAppliedMathematicsandMechanics, 2005,1:100～119

14 Liu L, Wong Y S, Lee B H K. Application of the centre manifold theory in non-linear aeroelasticity.JournalofSoundandVibration, 2000,234(4):641～659

15 Levitas J, Weller T, Singer J. Poincare-like simple cell mapping for non-linear dynamical systems.JournalofSoundandVibration, 1994,176:641～662

16 Raghothama A, Narayanan S. Non-linear dynamics of a two-dimensional airfoil by incremental harmonicbalance method.JournalofSoundandVibration, 1999,226:493～517

17 Dimitriadis G. Bifurcation analysis of aircraft with structural nonlinearity and freeplay using numerical continuation.JournalofAircraft, 2008,45(3):893～905

18 沈華勛. 飛翼式柔性飛機縱向動力學建模與穩(wěn)定性分析. 動力學與控制學報, 2016,14(3):241～246 (Sheng H X. Longitudinal danamic modeling and stability analysis of very flexible flying wings.JounalofDynamicsandControl, 2016,14(3):241～246 (in Chinese))

19 Kubicek M. Dependence of solution of nonlinear systems on a parameter.ACMTransactionsonMathematicalSoftware2, 1976:98～107

20 Doedel E J, Champneys A R, Fairgrieve T F, et al. AUTO 97-AUTO 2000: Continuation and bifurcation software for ordinary differential equations (with HomCont): User′s Guide, Technology Report, Concordia University, Montreal, Canada, 2000

21 Govaerts W, Dhooge A, Kuznetsov Y A. MATCONT: A matlab package for numerical bifurcation of ODEs.ACMTransactionsonMathematicalSoftware, 2003,29(2):141～16

22 Roberts I, Jones D P, Lieven N A J, et al. Analysis of piecewise linear aeroe-lastic systems using numerical continuation.ProceedingsoftheIMechEPartGJournalofAerospaceEngineering, 2002,216(1):1～11