閆 蓋,方明霞,楊英豪,沈海軍
(1. 上海第二工業(yè)大學(xué) 工學(xué)部, 上海 201209; 2. 同濟大學(xué) 航空航天與力學(xué)學(xué)院, 上海 200092;3. 上海飛機設(shè)計研究院, 上海 201210)
氣動彈性是飛行器設(shè)計中備受關(guān)注的重要問題,主要氣動彈性現(xiàn)象有顫振、抖振、動力響應(yīng)等[1-2]。其中,動力響應(yīng)是指彈性系統(tǒng)受到與自身系統(tǒng)無關(guān)的、隨時間任意變化的外界干擾力作用而發(fā)生的強迫振動,外界擾動可以是諧和的、周期的、脈沖的或隨機的[3]。在航空動力學(xué)領(lǐng)域中,隨機擾動普遍存在,它主要是由飛行器飛行過程中空氣熱力、風(fēng)力和尾渦等隨機因素相互作用形成的大氣湍流引起的。從強度觀點來看,飛機結(jié)構(gòu)可能在嚴(yán)重的湍流中由于超載而遭破壞,中等大小的湍流則是飛機結(jié)構(gòu)疲勞損傷的主要來源[4]。因此,結(jié)合氣流擾動的機翼系統(tǒng)動力學(xué)特性研究能更好地反映實際工況。近年來,考慮外界擾動的機翼氣動彈性問題正受到越來越多學(xué)者的關(guān)注。
Poirel和Price[5-8]研究了二元機翼在湍流隨機擾動作用下的顫振問題,考慮了機翼結(jié)構(gòu)線性和立方非線性情況,將湍流近似為高斯隨機過程,從概率密度和最大Lyapunov指數(shù)分析了二元機翼的隨機動力學(xué)行為,指出了隨機激勵下的顫振點較確定性系統(tǒng)顫振點提前,且與激勵的強度密切相關(guān)。文獻[9]采用能量隨機平均法,求解FPK(Fokker-planck-kolmogorov)方程,獲得隨機激勵下的一元機翼分岔點,并分析了系統(tǒng)的隨機Hopf分岔特性。文獻[10]利用攝動法獲得了二元機翼在非高斯色噪聲作用下的Lyapunov指數(shù),分析了二元機翼在隨機激勵下的穩(wěn)定性。文獻[11]采用隨機平均法得到二元機翼在寬帶噪聲激勵下的Lyapunov指數(shù),探討了隨機噪聲譜密度對機翼穩(wěn)定性的影響。文獻[12]采用隨機減量法和矩陣束法識別紊流激勵中的模態(tài)參數(shù),提出了一種紊流激勵下的顫振邊界預(yù)測方法??梢钥闯?,對隨機激勵下的機翼動力學(xué)特性研究主要集中在隨機顫振預(yù)測、隨機分岔特性研究中。
而對飛行器混沌運動問題的研究,目前主要針對考慮飛行器結(jié)構(gòu)非線性的自治系統(tǒng)。如文獻[13]采用伽遼金法對考慮幾何非線性的長直機翼運動方程進行離散,通過數(shù)值方法分析了機翼的顫振及混沌運動特性。文獻[14]研究分析飛行器操縱面的操縱剛度對混沌運動特性影響較小,而阻尼對系統(tǒng)的混沌特性影響較大。文獻[15]采用數(shù)值模擬方法和預(yù)測程序,研究了不可壓縮流中具有結(jié)構(gòu)二次、三次非線性項的二元機翼系統(tǒng)分岔和混沌特性??梢钥闯觯壳暗难芯慷嗍菑臋C翼非線性結(jié)構(gòu)參數(shù)出發(fā),研究系統(tǒng)的顫振、抖振問題,旨在提高系統(tǒng)的臨界飛行速度、優(yōu)化系統(tǒng)的結(jié)構(gòu)參數(shù)、改善控制方法等。但從外在隨機激勵角度探討二元機翼系統(tǒng)的復(fù)雜動力學(xué)特性較少,且采用解析方法對超聲速二元機翼的隨機混沌特性進行研究更為少見。為此,本文將從定性和定量角度出發(fā),根據(jù)Kapitaniak對隨機混沌的定義(即隨機混沌過程必須滿足兩個條件:①概率密度函數(shù)具有多個最大值;②概率時差圖具有康托集合結(jié)構(gòu)),提出一種半解析方法探討隨機激勵下二元機翼的混沌運動特性。即通過聯(lián)合使用累積量截斷法、非高斯截斷法獲得系統(tǒng)的二維聯(lián)合概率密度函數(shù)及系統(tǒng)的概率時差圖,采用Kapitaniak方法分析系統(tǒng)在不同擾動強度下的隨機混沌特性,并通過數(shù)值方法對分析結(jié)果進行驗證。本文研究不僅可以分析隨機激勵下擾動強度、飛行器結(jié)構(gòu)參數(shù)對系統(tǒng)混沌域的影響,還可以推廣到其他隨機非線性的大型復(fù)雜系統(tǒng)的動力學(xué)研究中,推動高維非線性系統(tǒng)動力學(xué)行為基礎(chǔ)研究的發(fā)展。因此本文研究具有重要的理論和實際意義。
圖1 二元機翼模型Fig.1 Model of a two-dimensional airfoil
采用第二類拉氏方程獲得受隨機擾動的二元機翼動力學(xué)方程為:
(1)
文獻[16-17]研究表明,在飛行馬赫數(shù)為2~5時,活塞理論比較適合機翼氣動力計算。為此本文采用三階活塞理論給出了機翼非線性氣動力和氣動力矩:
(2)
(3)
(4)
A0矩陣中各元素表達式如下:
4χαc∞ρηb2-χαbch)
其中
S1=-Fh+σhξh(t)
S2=-kα2α3-Fα+σαξα(t)
由于方程(4)為4維非線性動力學(xué)方程,直接對其動力學(xué)特性進行研究比較困難,為此首先采用中心流形方法對系統(tǒng)進行降維。由于中心流形定理僅適用于自治系統(tǒng),因此首先要對隨機激勵進行變換,本文采用Monte Carlo法將功率譜密度函數(shù)變換成多項余弦函數(shù)疊加的形式,并通過擴大向量將非自治系統(tǒng)轉(zhuǎn)換為自治系統(tǒng)。
通過Routh-Hurwitz判據(jù)得到:Ma<3.985時,系統(tǒng)穩(wěn)定,響應(yīng)收斂到平衡點;3.985
取μ=Ma-Ma*(Ma*=3.985)為分岔參數(shù),給定系統(tǒng)參數(shù),引入非奇異變換x=py,其中y=(y1,y2,y3,y4)T,p為相對應(yīng)于系統(tǒng)中零平衡點的Jacobi矩陣的第1項實部、第2項實部和虛部以及第4項實部構(gòu)成的方陣,代入方程后系統(tǒng)化成
(5)
根據(jù)中心流形定理,得到中心流形函數(shù)
(6)
o(3)表示三階以上的高階小量,忽略此量,可得到約化方程(7),其中Δ1、Δ2是經(jīng)過變換的隨機激勵。
(7)
根據(jù)Kapitaniak對隨機混沌特性的定義:二維聯(lián)合概率密度分布具有多峰且概率時差圖具有典型的康托集合結(jié)構(gòu)時,系統(tǒng)具有隨機混沌特性。為此首先對系統(tǒng)的二維聯(lián)合概率密度函數(shù)進行分析。由于方程(7)為非線性方程,故本文聯(lián)合利用累積量截斷法、非高斯截斷法對系統(tǒng)的二維聯(lián)合概率密度函數(shù)進行求解。方程(7)的伊藤隨機微分方程形式如下:
(8)
式中,D1、D1為隨機激勵強度。
由式(8)得到相應(yīng)的FPK方程,求得矩方程表達式為:
(10)
(11)
現(xiàn)通過二維聯(lián)合概率密度函數(shù)曲線、概率時差圖對二元機翼的隨機動力學(xué)特性進行研究。由于機翼俯仰自由度的振動幅值比沉浮自由度的振動幅值大,故本文主要研究俯仰角-俯仰角速度的概率密度曲線。現(xiàn)將系統(tǒng)發(fā)生顫振時的來流馬赫數(shù)Ma=3.985作為分界線,以此將系統(tǒng)劃分為顫振前區(qū)、顫振后區(qū),并研究來流馬赫數(shù)分別為Ma=3、Ma=5時系統(tǒng)顫振前區(qū)和顫振后區(qū)的動力學(xué)特性。
根據(jù)第3節(jié)求解二維概率密度函數(shù)的方法,獲得系統(tǒng)顫振前區(qū)和顫振后區(qū)在隨機擾動作用下的概率密度函數(shù)。隨機擾動強度根據(jù)文獻[17]添加,即弱隨機擾動取值0.01,一般隨機擾動取值0.1,強隨機擾動取值0.5。圖2是在顫振前區(qū)Ma=3時,在隨機擾動強度分別為0.01、0.1和
(a) D=0.01
(b) D=0.1
(c) D=0.5圖2 不同隨機擾動強度下的二維概率密度圖(Ma=3)Fig.2 Two-dimensional probability density diagram with various disturbance intensities(Ma=3)
0.5時系統(tǒng)的二維概率密度分布圖;圖3是在顫振后區(qū)Ma=5時不同隨機擾動強度下的二維概率密度分布圖。
(a) D=0.01
(b) D=0.1
(c) D=0.5圖3 不同隨機擾動強度下的二維概率密度圖(Ma=5)Fig.3 Two-dimensional probability density diagram with various disturbance intensities(Ma=5)
從圖2、圖3可以看出:二維概率密度在顫振前區(qū)時,形狀在弱隨機擾動下為分離雙峰,一般隨機擾動下為相鄰雙峰,強隨機擾動下變?yōu)閱畏?;在顫振后區(qū)時,形狀在弱隨機擾動下為多峰,一般隨機擾動下變?yōu)殡p峰,強隨機擾動下變?yōu)閱畏?。并且在顫振前區(qū)和顫振后區(qū)的二維聯(lián)合概率密度形狀在弱隨機擾動下不同,而在強隨機擾動下形狀相似,由此可知系統(tǒng)拓撲結(jié)構(gòu)發(fā)生了質(zhì)變,系統(tǒng)發(fā)生了P分岔。
根據(jù)概率密度的多峰狀態(tài),為了判斷系統(tǒng)在弱隨機擾動、一般隨機擾動下是否進入了混沌運動狀態(tài),現(xiàn)繪制概率時差圖對系統(tǒng)的動力學(xué)特性進行進一步分析,因為二元機翼在一定來流馬赫數(shù)下,受一定強度的隨機激勵影響發(fā)生分岔、進入混沌運動狀態(tài)時間非常短,因此繪制概率時差圖的時延也較小,本文選取0.01 s。圖4為顫振前區(qū),在弱隨機擾動、一般隨機擾動下的概率時差圖;圖5為顫振后區(qū), 在弱隨機擾動、一般隨機擾動下的概率時差圖。
從圖4、圖5可以發(fā)現(xiàn),在顫振前區(qū)和顫振后區(qū),系統(tǒng)受弱隨機擾動強度和一般隨機擾動強度下,概率時差圖均出現(xiàn)了典型的康托集合效應(yīng),即在某一位置附近出現(xiàn)頻率明顯高于其他位置。結(jié)合其二維概率密度呈現(xiàn)多峰形狀,因此系統(tǒng)發(fā)生了Kapitaniak定義下的隨機混沌。
(a) D=0.01
(b) D=0.1圖4 不同擾動強度下顫振前區(qū)的概率時差圖(Ma=3)Fig.4 Probability time difference diagrams with various disturbance intensities(Ma=3)
(a) D=0.01
(b) D=0.1 圖5 不同擾動強度下顫振后區(qū)的概率時差圖(Ma=5)Fig.5 Probability time difference diagrams with various disturbance intensities(Ma=5)
由于強隨機擾動下二維聯(lián)合概率密度函數(shù)為單峰,根據(jù)Kapitaniak定義,此時系統(tǒng)不會發(fā)生Kapitaniak定義下的隨機混沌。
綜上可得,不管在顫振前區(qū)還是顫振后區(qū),受外部隨機激勵的影響系統(tǒng)動態(tài)特性都會發(fā)生突變、發(fā)生分岔,甚至進入混沌運動狀態(tài)。但也會在一定隨機擾動強度下從混沌運動狀態(tài)突變回周期運動。因此,在進行二元機翼氣動彈性特性研究時,要考慮不同飛行環(huán)境狀態(tài)的影響,避免系統(tǒng)發(fā)生復(fù)雜動力學(xué)行為。
為了驗證分析結(jié)果的有效性,現(xiàn)采用數(shù)值方法對分析結(jié)果進行驗證。首先采用Monte Carlo模擬法對隨機激勵進行模擬,然后利用與前文一致的參數(shù)取值,在MATLAB平臺上對狀態(tài)方程(4)進行數(shù)值求解。不同工況下系統(tǒng)的時域響應(yīng)曲線、龐加萊截面圖如圖6、圖7所示,系統(tǒng)的最大Lyapunov指數(shù)隨來流馬赫數(shù)Ma的變化曲線,如圖8所示。
(a) 俯仰角時程曲線(D=0.01)(a) Time history curve of pitch angle(D=0.01) (b) 龐加萊截面(D=0.01)(b) Poincare section(D=0.01)
(c) 俯仰角時程曲線(D=0.1)(c) Time history curve of pitch angle(D=0.1) (d) 龐加萊截面(D=0.1)(d) Poincare section(D=0.1)圖6 不同擾動強度下顫振前區(qū)的系統(tǒng)響應(yīng)曲線及龐加萊截面(Ma=3)Fig.6 Response of the system with various disturbance intensities and Poincare section(Ma=3)
(a) 俯仰角時程曲線(D=0.01)(a) Time history curve of pitch angle(D=0.01) (b) 龐加萊截面(D=0.01)(b) Poincare section(D=0.01)
(c) 俯仰角時程曲線(D=0.1)(c) Time history curve of pitch angle(D=0.1) (d) 龐加萊截面(D=0.1)(d) Poincare section(D=0.1)圖7 不同擾動強度下顫振后區(qū)的系統(tǒng)響應(yīng)曲線及龐加萊截面(Ma=5)Fig.7 Response of the system with various disturbance intensities and Poincare section(Ma=5)
(a) D=0.01
(b) D=0.1
(c) D=0.5圖8 最大lyapunov指數(shù)變化曲線圖Fig.8 The largest lyapunov exponent changes with Ma
從圖6~8可以看出,在弱隨機擾動、一般隨機擾動強度作用下,無論是顫振前區(qū)還是顫振后區(qū),系統(tǒng)的俯仰角隨時間均呈無規(guī)律變化,同時龐加萊截面由分形結(jié)構(gòu)的大量散點組成,且系統(tǒng)在Ma=3、Ma=5時最大Lyapunov指數(shù)均為正值,說明系統(tǒng)進入了混沌運動狀態(tài);而在隨機擾動強度為0.5時,系統(tǒng)在Ma=3、Ma=5時最大Lyapunov指數(shù)均為負值,說明系統(tǒng)未進入混沌運動狀態(tài),驗證了前文中Kapitaniak定義下的隨機混沌分析結(jié)果,說明了本文近似解析定性分析方法的準(zhǔn)確性和有效性。
本文通過累積量截斷法、非高斯截斷法及Edgeworth展式等獲得系統(tǒng)的聯(lián)合概率密度函數(shù),在Kapitaniak定義下研究了二元機翼的隨機混沌特性,并通過數(shù)值方法對計算結(jié)果進行驗證。主要研究結(jié)果如下:
1)采用第二類拉氏方程建立二元機翼隨機非線性動力學(xué)方程,并通過三階活塞理論推導(dǎo)二元機翼的非線性氣動力和氣動力矩。通過中心流形定理對系統(tǒng)進行降維,將2自由度下4維非線性系統(tǒng)的動力學(xué)方程成功降為2維,使系統(tǒng)的隨機非線性動力學(xué)特性研究得到簡化。
2)采用累積量截斷法和非高斯截斷法獲得系統(tǒng)二維聯(lián)合概率密度函數(shù),利用Kapitaniak定義判斷系統(tǒng)的隨機混沌特性。研究發(fā)現(xiàn),當(dāng)隨機擾動強度為0.01和0.1時,系統(tǒng)在顫振前區(qū)和顫振后區(qū)均進入了混沌運動狀態(tài);在隨機擾動強度為0.5時,系統(tǒng)沒有進入混沌運動狀態(tài)。這說明外部隨機擾動對系統(tǒng)具有很大影響。因此,在二元機翼氣動彈性特性控制時,必須充分考慮外部隨機激勵強度的影響,以增加系統(tǒng)的穩(wěn)定性和安全性。
3)通過時域響應(yīng)曲線、龐加萊截面圖和最大Lyapunov指數(shù)等數(shù)值方法對計算結(jié)果進行驗證,兩種具有較好的一致性,說明本文研究具有較高的研究精度。