基于時滯慣性流形的二維平面壁板非線性氣動彈性分析

梅冠華，張家忠，席 光

（西安交通大學 能源與動力工程學院，西安 710049）

氣動彈性領(lǐng)域主要研究氣動力和彈性結(jié)構(gòu)變形之間的相互作用。氣動力產(chǎn)生結(jié)構(gòu)變形，而結(jié)構(gòu)變形又影響著氣動力，這種氣動力和變形間的耦合作用導致了流體和固體復雜的動力學響應。如顫振、混沌等。這些響應的深入分析對于揭示流固耦合問題的本質(zhì)，對于該類非線性振動的避免、控制和利用都具有十分重要的意義［1］。

氣動彈性問題中的非線性來源于氣動力和結(jié)構(gòu)體兩方面：流場中的激波、分離會導致結(jié)構(gòu)變形與氣動力的非線性關(guān)系；幾何大變形、材料非線性會引起結(jié)構(gòu)體受力與其變形的非線性關(guān)系［2-4］。理論上為準確描述這種非線性氣動彈性問題，需綜合考慮氣動力和結(jié)構(gòu)體的非線性效應。對于目前的氣動彈性分析，一部分考慮了氣動力的非線性，結(jié)構(gòu)體采用線性假設(shè)；一部分考慮了結(jié)構(gòu)體的非線性，氣動力采用線性假設(shè)；綜合考慮兩方面非線性特性的分析較為少見。

作為典型的非線性氣動彈性問題，對于二維平面壁板顫振模型，平板變形采用Von Karman幾何大變形理論，氣動力采用修正一階活塞理論近似，推導出系統(tǒng)的非線性偏微分控制方程。Dowell［5］運用傳統(tǒng)Galerkin方法（TGM）研究了系統(tǒng)不同控制參數(shù)對于顫振邊界和振動幅值的影響，并與實驗結(jié)果進行對比。然而模態(tài)數(shù)目的選取是TGM中的關(guān)鍵問題，若選取過少則不能獲得精確解，若選取過多則耗費大量計算時間。慣性流形的提出簡化了系統(tǒng)離散方程的維數(shù)，在保證計算精度的同時減少了計算量。慣性流形（IM）、近似慣性流形（AIM）和時滯慣性流形（IMD）的概念反映了非線性發(fā)展方程的高階模態(tài)與低階模態(tài)間的相互作用規(guī)律［6-8］。IMD提出非線性發(fā)展方程的高階模態(tài)不僅依賴于當前低階模態(tài)分量值，還依賴于高階分量某個時間間隔以前的值［9-10］。IMD較AIM更加符合實際情況，不僅存在條件大為寬松，且具有良好的收斂性和穩(wěn)定性［11］。基于IMD思想的非線性Galerkin方法被引入到二維平面壁板顫振模型的求解中，與TGM結(jié)果對比發(fā)現(xiàn)精度不變且節(jié)省計算時間和資源。

對計算結(jié)果，分別以無量綱動壓和無量綱壓縮內(nèi)力為分岔參數(shù)，無量綱振幅為響應作出分岔圖形，發(fā)現(xiàn)陣發(fā)性通向混沌的途徑，混沌區(qū)域表現(xiàn)出周期窗口和自相似特性。通過對系統(tǒng)的相圖、位移的FFT頻譜以及Lyapunov指數(shù)的分析，驗證系統(tǒng)存在穩(wěn)定、屈曲、諧調(diào)和非諧調(diào)運動四種典型類型。對非諧調(diào)運動的詳細研究發(fā)現(xiàn)系統(tǒng)存在倍周期運動、準周期運動和混沌運動等豐富的非線性響應。研究結(jié)果為識別和進一步控制此類非線性現(xiàn)象提供了依據(jù)。

1 力學模型

二維平面壁板顫振模型如圖1所示：平板長度為a，為方便理論分析，假設(shè)寬度b→∞，則問題簡化成二維平板，兩端初始拉伸內(nèi)力N0，平板厚度為h，單位長度質(zhì)量為ρ，彈性模量E，泊松比μ，上表面沿x方向來流速度為U∞，馬赫數(shù)M∞?1，空氣密度ρ∞，邊界條件為兩端簡支。

在Kirchhoff-Love假設(shè)下，僅考慮板的橫向振動w（x，t），應變位移關(guān)系采用Von Karman幾何非線性大變形理論表示，給出平板運動方程［5］：

圖1 二維平面壁板顫振模型Fig.1 2-Dimensional panel flutter model

引入無量綱參數(shù)：

系統(tǒng)控制方程（1）無量綱化為：

邊界條件（4）無量綱化為：

2 數(shù)值求解方法

根據(jù)傳統(tǒng)Galerkin方法，構(gòu)造滿足邊界條件（8）的基函數(shù)列：sin（iπ），i=1，2，…。平板的位移可以設(shè)為此函數(shù)列的疊加形式：

得到k時刻后N個高階模態(tài)方程組為：

此代數(shù)方程組實現(xiàn)了后N個高階模態(tài)用前N個低階模態(tài)表示，且加入了時間滯后。這樣后N個高階模態(tài)不用通過繁瑣的數(shù)值積分，而是直接求解代數(shù)方程組得到，從而節(jié)省了計算時間。

3 數(shù)值分析

為了獲得系統(tǒng)的精確解至少需要4到6階模態(tài)，對于4階以上模態(tài)系統(tǒng)的理論分析較為困難，多以數(shù)值模擬為主。系統(tǒng)的三個控制參數(shù)/M∞、無量綱拉伸內(nèi)力Rx和無量綱動壓λ中∞影響微小，Rx和λ影響顯著。由于壓縮內(nèi)力是系統(tǒng)趨于不穩(wěn)定的因素，因為方便分析，定義無量綱壓縮內(nèi)力R=-Rx/π2，注意R與Rx的方向相反。取定∞=0.01，分別選取λ和R作為分岔參數(shù)，采用6階模態(tài)進行數(shù)值計算，以=0.75為參考點，計算時間充分長以消除瞬態(tài)響應，系統(tǒng)達到穩(wěn)態(tài)響應后，若為平面穩(wěn)定或屈曲則取其無量綱變形，若為振動則取其無量綱速度=0時刻的無量綱振幅，即取相圖上=0直線與系統(tǒng)軌道交點處的值，以此為響應繪制系統(tǒng)分岔圖。對于特定參數(shù)下系統(tǒng)響應，采用圖相圖、穩(wěn)態(tài)的FFT頻譜圖以及Lyapunov指數(shù)綜合分析。

頻譜圖中橫坐標為無量綱頻率，縱坐標為幅值，其反映系統(tǒng)所含不同頻率正弦振動分量所占比重，直觀上若頻譜圖中僅有單個脈沖的峰值，則僅含單獨頻率分量，為諧調(diào)運動；若頻譜圖中存在有限個脈沖峰值，則為有限個頻率下正弦振動的疊加，若頻率間呈倍數(shù)關(guān)系，即存在公約數(shù)，則為倍周期運動，若頻率間不存在倍數(shù)關(guān)系，即頻率成分不可互約，則為準周期運動；若存在無窮個峰值密布于頻譜圖上，則此狀態(tài)為混沌，在相圖上表現(xiàn)為有限區(qū)域內(nèi)雜亂的非封閉曲線。

混沌的更精確代數(shù)判定可用Lyapunov指數(shù)［14］，其反應系統(tǒng)初始微小偏差隨時間趨于無窮的指數(shù)膨脹或收縮狀況，正值表示膨脹，零值表示周期運動，負值表示收縮，其絕對值大小表示收縮或膨脹程度的強弱。一般而言系統(tǒng)最大Lyapunov指數(shù)為正值則可認為達到混沌狀態(tài)。采用經(jīng)典單位球演化伴隨GSR正交化置換新矢量的算法計算六階模態(tài)下系統(tǒng)全部Lyapunov指數(shù)，若出現(xiàn)正值則認為達到混沌狀態(tài)，不同參數(shù)下混沌系統(tǒng)的前三個Lyapunov指數(shù)如表1所示。

表1 不同參數(shù)下混沌系統(tǒng)的前三個Lyapunov指數(shù)Tab.1 First three Lyapunov Exponents of chaos systems at different parameters

3.1 TGM與IMD計算結(jié)果的比較

R=0即平板兩端不受初始內(nèi)力的狀況下，當無量綱動壓λ逐漸增大跨過臨界值λcr=345后，平板發(fā)生Hopf分岔，由靜態(tài)穩(wěn)定轉(zhuǎn)化為極限環(huán)振動，且振幅隨λ的增大而增大，如圖2（a）所示。FFT頻譜分析表明系統(tǒng)僅含單獨頻率，且該頻率隨著λ的增大而增大，如圖2（b）所示。與TGM結(jié)果對比發(fā)現(xiàn)IMD結(jié)果保留原有精度。

圖 2 R=0，振幅和頻率隨λ變化圖Fig.2 Amplitude and frequencyvs.λ at R=0

λ=0即無氣動載荷狀況下，隨R的變化如圖3所示。當R逐漸增大跨過臨界值Rcr時，平板由靜態(tài)穩(wěn)定轉(zhuǎn)化為屈曲失穩(wěn)，到達兩個變形位置中的一個。且變形量隨R增大而增大。對比發(fā)現(xiàn)IMD結(jié)果與TGM結(jié)果吻合較好。

3.2 λ變化對系統(tǒng)響應的影響

固定無量綱壓縮內(nèi)力R=4，考慮無量綱動壓λ變化對系統(tǒng)響應的影響。以λ為分岔參數(shù)，為響應的分岔圖如圖4所示。

當0＜λ＜114時，系統(tǒng)表現(xiàn)為屈曲，到達兩個變形位置中的一個，以λ=100為例如圖5（a）所示；

當114＜λ＜159時，系統(tǒng)演變?yōu)榛煦纾驭?120為例如圖5（b）所示，前兩個Lyapunov指數(shù)為正值如表1所示；

當159＜λ＜173時，系統(tǒng)表現(xiàn)為倍周期運動。以λ=160為例如圖5（c）所示，首頻率1.59 Hz占據(jù)主要地位，各頻率間呈現(xiàn) 1、2、3、4、5 倍頻關(guān)系，2 倍和 5 倍頻貢獻十分微小，如表2所示；

當173＜λ＜227時，系統(tǒng)演化為混沌狀態(tài)，以λ=210為例如圖5（d）所示，其前兩個Lyapunov指數(shù)為正值如表1所示；

當227＜λ＜249時，系統(tǒng)發(fā)生倒分岔最終演化為單周期振動，以λ=260為例如圖5（e）所示，首頻率占據(jù)了絕對統(tǒng)治地位如表2所示。

圖3 λ =0，屈曲變形隨R變化圖Fig.3 Buckling deformation vs.R at λ =0

表2 R=4，不同λ下系統(tǒng)主要頻率分量及其幅值Tab.2 Main frequencies and their amplitudes of different λ at R=4

圖4 R=4，隨λ的分岔圖Fig.4 Bifurcation diagram of with λ at R=4

3.3 R變化對系統(tǒng)響應的影響

固定無量綱動壓λ=200，考慮無量綱壓縮內(nèi)力R變化對系統(tǒng)響應的影響。以R為分岔參數(shù)，w為響應繪制分岔圖如圖6所示，其中圖6（a）為整體圖，圖6（b）為局部放大圖。

圖5 R=4，不同λ下系統(tǒng)相圖及頻譜分析Fig.5 System phase portraits and spectra analysis of different λ at R=4

當0≤R≤1.875時，系統(tǒng)表現(xiàn)為平面穩(wěn)定，以R=1.0為例如圖7（a）所示，對于任何初始擾動，系統(tǒng)最終回復到靜態(tài)平衡位置；

當1.876≤R≤3.325 時，系統(tǒng)發(fā)生 Hopf分岔，表現(xiàn)為諧調(diào)振動，以R=2.5為例如圖7（b）所示，振動僅含單獨頻率3.29 Hz；

當3.326≤R≤3.623時，系統(tǒng)發(fā)生一系列倍周期-準周期間歇型的分岔，即為陣發(fā)性途徑最終通向混沌：以R=3.5，R=3.55為代表的倍周期運動如圖7（c）～圖7（d）所示，其各頻率成分間呈現(xiàn)出倍數(shù)關(guān)系如表3 所示；以R=3.43，R=3.54 和R=3.60 為代表的準周期運動如圖7（e）～圖7（g）所示，其各頻率成份間并不存在倍數(shù)關(guān)系，如表3所示。這是由于系統(tǒng)的非線性特性，不同頻率的模態(tài)間存在耦合，導致在特定狀態(tài)下單獨模態(tài)或多個模態(tài)占據(jù)主導地位，系統(tǒng)則表現(xiàn)為單頻率或是多個頻率組合的振動形式，故而出現(xiàn)如上復雜的振動行為。

當3.624≤R≤5時，系統(tǒng)通過上述間歇型分岔最終進混沌狀態(tài)，混沌圖形中的自相似結(jié)構(gòu)以及周期亮窗等特性呈現(xiàn)在分岔局部放大圖6（b）中。以R=4.5為例的混沌運動如圖7（h）所示，其前兩個Lyapunov指數(shù)均大于零如表1所示。

表3 λ=200不同R下系統(tǒng)主要頻率分量及其幅值Tab.3 Main frequencies and their amplitudes of different R at λ=200

圖6 λ=200，隨R分岔圖Fig.6 Bifurcation diagram of with R at λ=200

4 結(jié)論

（1）將基于IMD的非線性Galerkin方法引入到二維平面壁板顫振模型的求解中，結(jié)果表明IMD在保留了TGM精度的同時，可以縮減求解維數(shù)，大量節(jié)省計算時間，為IMD求解類似問題提供了借鑒。

（2）系統(tǒng)的兩個主要參數(shù)，即λ和R將響應分為多個區(qū)域，即：二維平面壁板氣動彈性系統(tǒng)在不同參數(shù)下表現(xiàn)出平面穩(wěn)定、屈曲變形失穩(wěn)、諧調(diào)振動、倍周期振動、準周期振動以及混沌等豐富的非線性現(xiàn)象。

（3）分別以λ和R為分岔參數(shù)，為響應得出系統(tǒng)分岔圖形，通過對其分析發(fā)現(xiàn)存在陣發(fā)性通向混沌的途徑，觀察到混沌中的自相似結(jié)構(gòu)以及周期窗口現(xiàn)象。

（4）Lyapunov指數(shù)是精確判斷混沌的代數(shù)判據(jù)，出現(xiàn)正值時即可認定系統(tǒng)處于混沌狀態(tài)，F(xiàn)FT所得頻譜圖反映了系統(tǒng)所含頻率分量及其所占比重，可以定性地識別系統(tǒng)響應。

（5）本文所采用的方法可以用于指導該類非線性氣動彈性問題機理的進一步深入研究，所得出的結(jié)論可對該類現(xiàn)象的識別和控制提供依據(jù)。進而合理利用自激振動現(xiàn)象實現(xiàn)延遲激波產(chǎn)生、延遲分離點和增升減阻等效果。

［1］Kang W，Zhang J Z，F(xiàn)eng P H.Numerical simulation and Aeroelastic Analysis of a Local Flexible Airfoil at Low Reynolds Numbers［A］.In：The 8th Asian Computational Fluid Dynamics Conference［C］.Hong Kong，China，2010：10-14.

［2］Dowell E H. Some recent advances in nonlinear aeroelasticity：fluid-structure interaction in the 21st century［A］.In：51st AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Conference，18th AIAA/ASME/AHS Adaptive Structures Conference［C］//Orlando，F(xiàn)lorida，USA，2010：12-15.

［3］Dowell E H，Edwards J，Strganac T.Nonlinear aeroelasticity［J］.Journal of Aircraft，2003，40（5）：857-874.

［4］Dowell E H，Tang D.Nonlinear aeroelasticity and unsteady aerodynamics［J］.AIAA Journal，2002，40（9）：1697-1707.

［5］ Dowell E H.Nonlinear oscillations of a fluttering plate I［J］.AIAA Journal，1966，4（7）：1267-1275.

［6］ Zhang J Z，Liu Y，Lei P F，et al.Dynamic snap-through buckling analysis of shallow arches under impact load based on approximate inertialmanifolds［J］. Dynamics of Continuous， Discreteand ImpulsiveSystems， SeriesB（DCDIS-B），2007，14（S5）：287-291.

［7］Foias C， SellG R. Inertialmanifolds fornonlinear evolutionary equations［J］.Journal of Differential Equations，1988，73（2）：309-353.

［8］ Foias C，Manley O，Temam R.On the interaction of small eddies in Ttwo-dimensional turbulence flows［J］.Math Modeling and Numerical Analysis，1988，M 2AN，22（1）：93-114.

［9］ Debussche A，Temam R.Inertial manifolds with delay［J］.Appl.Math.Lett，1995，8（2）：21-24.

［10］戴正德，郭柏靈.慣性流形與近似慣性流形［M］.北京：科學出版社，2000.

［11］張家忠，陳麗鶯，梅冠華，等.基于時滯慣性流形的淺拱動力屈曲研究［J］.振動與沖擊，2009，28（6）：100-103.

［12］Ashley H，Zartarian G.Piston theory-a new aerodynamic tool for the aaeroelastician［J］.Journal of the Aeronautical Science，1956，23（12）：1109-1118.

［13］ Zhang W W，Ye Z Y，Zhang C A，et al.Supersonic flutter analysis based on a local piston theory［J］.AIAA Journal，2009，47（10）：2321-2328.

［14］ Wolf A，Swift J B，Swinney H L，et al.Determining lyapunov exponents from a time series［J］.Physica D，1985，16（3）：285-317.