閆 蓋,方明霞,楊英豪,沈海軍

(1. 上海第二工業(yè)大學(xué) 工學(xué)部， 上海 201209； 2. 同濟(jì)大學(xué) 航空航天與力學(xué)學(xué)院， 上海 200092；3. 上海飛機(jī)設(shè)計(jì)研究院， 上海 201210)

氣動彈性是飛行器設(shè)計(jì)中備受關(guān)注的重要問題，主要?dú)鈩訌椥袁F(xiàn)象有顫振、抖振、動力響應(yīng)等[1-2]。其中，動力響應(yīng)是指彈性系統(tǒng)受到與自身系統(tǒng)無關(guān)的、隨時間任意變化的外界干擾力作用而發(fā)生的強(qiáng)迫振動，外界擾動可以是諧和的、周期的、脈沖的或隨機(jī)的[3]。在航空動力學(xué)領(lǐng)域中,隨機(jī)擾動普遍存在，它主要是由飛行器飛行過程中空氣熱力、風(fēng)力和尾渦等隨機(jī)因素相互作用形成的大氣湍流引起的。從強(qiáng)度觀點(diǎn)來看，飛機(jī)結(jié)構(gòu)可能在嚴(yán)重的湍流中由于超載而遭破壞，中等大小的湍流則是飛機(jī)結(jié)構(gòu)疲勞損傷的主要來源[4]。因此，結(jié)合氣流擾動的機(jī)翼系統(tǒng)動力學(xué)特性研究能更好地反映實(shí)際工況。近年來，考慮外界擾動的機(jī)翼氣動彈性問題正受到越來越多學(xué)者的關(guān)注。

Poirel和Price[5-8]研究了二元機(jī)翼在湍流隨機(jī)擾動作用下的顫振問題，考慮了機(jī)翼結(jié)構(gòu)線性和立方非線性情況，將湍流近似為高斯隨機(jī)過程，從概率密度和最大Lyapunov指數(shù)分析了二元機(jī)翼的隨機(jī)動力學(xué)行為，指出了隨機(jī)激勵下的顫振點(diǎn)較確定性系統(tǒng)顫振點(diǎn)提前，且與激勵的強(qiáng)度密切相關(guān)。文獻(xiàn)[9]采用能量隨機(jī)平均法，求解FPK(Fokker-planck-kolmogorov)方程，獲得隨機(jī)激勵下的一元機(jī)翼分岔點(diǎn)，并分析了系統(tǒng)的隨機(jī)Hopf分岔特性。文獻(xiàn)[10]利用攝動法獲得了二元機(jī)翼在非高斯色噪聲作用下的Lyapunov指數(shù)，分析了二元機(jī)翼在隨機(jī)激勵下的穩(wěn)定性。文獻(xiàn)[11]采用隨機(jī)平均法得到二元機(jī)翼在寬帶噪聲激勵下的Lyapunov指數(shù)，探討了隨機(jī)噪聲譜密度對機(jī)翼穩(wěn)定性的影響。文獻(xiàn)[12]采用隨機(jī)減量法和矩陣束法識別紊流激勵中的模態(tài)參數(shù)，提出了一種紊流激勵下的顫振邊界預(yù)測方法?？梢钥闯觯瑢﹄S機(jī)激勵下的機(jī)翼動力學(xué)特性研究主要集中在隨機(jī)顫振預(yù)測、隨機(jī)分岔特性研究中。

而對飛行器混沌運(yùn)動問題的研究，目前主要針對考慮飛行器結(jié)構(gòu)非線性的自治系統(tǒng)。如文獻(xiàn)[13]采用伽遼金法對考慮幾何非線性的長直機(jī)翼運(yùn)動方程進(jìn)行離散，通過數(shù)值方法分析了機(jī)翼的顫振及混沌運(yùn)動特性。文獻(xiàn)[14]研究分析飛行器操縱面的操縱剛度對混沌運(yùn)動特性影響較小，而阻尼對系統(tǒng)的混沌特性影響較大。文獻(xiàn)[15]采用數(shù)值模擬方法和預(yù)測程序，研究了不可壓縮流中具有結(jié)構(gòu)二次、三次非線性項(xiàng)的二元機(jī)翼系統(tǒng)分岔和混沌特性?？梢钥闯?，目前的研究多是從機(jī)翼非線性結(jié)構(gòu)參數(shù)出發(fā)，研究系統(tǒng)的顫振、抖振問題，旨在提高系統(tǒng)的臨界飛行速度、優(yōu)化系統(tǒng)的結(jié)構(gòu)參數(shù)、改善控制方法等。但從外在隨機(jī)激勵角度探討二元機(jī)翼系統(tǒng)的復(fù)雜動力學(xué)特性較少，且采用解析方法對超聲速二元機(jī)翼的隨機(jī)混沌特性進(jìn)行研究更為少見。為此，本文將從定性和定量角度出發(fā)，根據(jù)Kapitaniak對隨機(jī)混沌的定義(即隨機(jī)混沌過程必須滿足兩個條件：①概率密度函數(shù)具有多個最大值；②概率時差圖具有康托集合結(jié)構(gòu))，提出一種半解析方法探討隨機(jī)激勵下二元機(jī)翼的混沌運(yùn)動特性。即通過聯(lián)合使用累積量截?cái)喾?、非高斯截?cái)喾ǐ@得系統(tǒng)的二維聯(lián)合概率密度函數(shù)及系統(tǒng)的概率時差圖，采用Kapitaniak方法分析系統(tǒng)在不同擾動強(qiáng)度下的隨機(jī)混沌特性，并通過數(shù)值方法對分析結(jié)果進(jìn)行驗(yàn)證。本文研究不僅可以分析隨機(jī)激勵下擾動強(qiáng)度、飛行器結(jié)構(gòu)參數(shù)對系統(tǒng)混沌域的影響，還可以推廣到其他隨機(jī)非線性的大型復(fù)雜系統(tǒng)的動力學(xué)研究中，推動高維非線性系統(tǒng)動力學(xué)行為基礎(chǔ)研究的發(fā)展。因此本文研究具有重要的理論和實(shí)際意義。

1 二元機(jī)翼隨機(jī)非線性動力學(xué)方程

圖1 二元機(jī)翼模型Fig.1 Model of a two-dimensional airfoil

采用第二類拉氏方程獲得受隨機(jī)擾動的二元機(jī)翼動力學(xué)方程為：

(1)

文獻(xiàn)[16-17]研究表明，在飛行馬赫數(shù)為2～5時，活塞理論比較適合機(jī)翼氣動力計(jì)算。為此本文采用三階活塞理論給出了機(jī)翼非線性氣動力和氣動力矩：

(2)

(3)

(4)

A0矩陣中各元素表達(dá)式如下：

4χαc∞ρηb2-χαbch)

其中

S1=-Fh+σhξh(t)

S2=-kα2α3-Fα+σαξα(t)

2 采用中心流形方法對系統(tǒng)進(jìn)行降維

由于方程(4)為4維非線性動力學(xué)方程，直接對其動力學(xué)特性進(jìn)行研究比較困難，為此首先采用中心流形方法對系統(tǒng)進(jìn)行降維。由于中心流形定理僅適用于自治系統(tǒng)，因此首先要對隨機(jī)激勵進(jìn)行變換，本文采用Monte Carlo法將功率譜密度函數(shù)變換成多項(xiàng)余弦函數(shù)疊加的形式，并通過擴(kuò)大向量將非自治系統(tǒng)轉(zhuǎn)換為自治系統(tǒng)。

通過Routh-Hurwitz判據(jù)得到：Ma<3.985時，系統(tǒng)穩(wěn)定，響應(yīng)收斂到平衡點(diǎn)；3.9854.086時，系統(tǒng)發(fā)散，因此，Ma*=3.985為系統(tǒng)Hopf分岔點(diǎn)。

取μ=Ma-Ma*(Ma*=3.985)為分岔參數(shù)，給定系統(tǒng)參數(shù)，引入非奇異變換x=py，其中y=(y1,y2,y3,y4)T，p為相對應(yīng)于系統(tǒng)中零平衡點(diǎn)的Jacobi矩陣的第1項(xiàng)實(shí)部、第2項(xiàng)實(shí)部和虛部以及第4項(xiàng)實(shí)部構(gòu)成的方陣，代入方程后系統(tǒng)化成

(5)

根據(jù)中心流形定理，得到中心流形函數(shù)

(6)

o(3)表示三階以上的高階小量，忽略此量，可得到約化方程(7)，其中Δ1、Δ2是經(jīng)過變換的隨機(jī)激勵。

(7)

3 二維聯(lián)合概率密度函數(shù)求解

根據(jù)Kapitaniak對隨機(jī)混沌特性的定義：二維聯(lián)合概率密度分布具有多峰且概率時差圖具有典型的康托集合結(jié)構(gòu)時，系統(tǒng)具有隨機(jī)混沌特性。為此首先對系統(tǒng)的二維聯(lián)合概率密度函數(shù)進(jìn)行分析。由于方程(7)為非線性方程，故本文聯(lián)合利用累積量截?cái)喾?、非高斯截?cái)喾▽ο到y(tǒng)的二維聯(lián)合概率密度函數(shù)進(jìn)行求解。方程(7)的伊藤隨機(jī)微分方程形式如下：

(8)

式中，D1、D1為隨機(jī)激勵強(qiáng)度。

由式(8)得到相應(yīng)的FPK方程，求得矩方程表達(dá)式為：

(10)

(11)

4 二元機(jī)翼隨機(jī)混沌特性分析

現(xiàn)通過二維聯(lián)合概率密度函數(shù)曲線、概率時差圖對二元機(jī)翼的隨機(jī)動力學(xué)特性進(jìn)行研究。由于機(jī)翼俯仰自由度的振動幅值比沉浮自由度的振動幅值大，故本文主要研究俯仰角-俯仰角速度的概率密度曲線。現(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ū)在隨機(jī)擾動作用下的概率密度函數(shù)。隨機(jī)擾動強(qiáng)度根據(jù)文獻(xiàn)[17]添加，即弱隨機(jī)擾動取值0.01，一般隨機(jī)擾動取值0.1，強(qiáng)隨機(jī)擾動取值0.5。圖2是在顫振前區(qū)Ma=3時，在隨機(jī)擾動強(qiáng)度分別為0.01、0.1和

(a) D=0.01

(b) D=0.1

(c) D=0.5圖2 不同隨機(jī)擾動強(qiáng)度下的二維概率密度圖(Ma=3)Fig.2 Two-dimensional probability density diagram with various disturbance intensities(Ma=3)

0.5時系統(tǒng)的二維概率密度分布圖；圖3是在顫振后區(qū)Ma=5時不同隨機(jī)擾動強(qiáng)度下的二維概率密度分布圖。

(a) D=0.01

(b) D=0.1

(c) D=0.5圖3 不同隨機(jī)擾動強(qiáng)度下的二維概率密度圖(Ma=5)Fig.3 Two-dimensional probability density diagram with various disturbance intensities(Ma=5)

從圖2、圖3可以看出：二維概率密度在顫振前區(qū)時，形狀在弱隨機(jī)擾動下為分離雙峰，一般隨機(jī)擾動下為相鄰雙峰，強(qiáng)隨機(jī)擾動下變?yōu)閱畏?；在顫振后區(qū)時，形狀在弱隨機(jī)擾動下為多峰，一般隨機(jī)擾動下變?yōu)殡p峰，強(qiáng)隨機(jī)擾動下變?yōu)閱畏?。并且在顫振前區(qū)和顫振后區(qū)的二維聯(lián)合概率密度形狀在弱隨機(jī)擾動下不同，而在強(qiáng)隨機(jī)擾動下形狀相似，由此可知系統(tǒng)拓?fù)浣Y(jié)構(gòu)發(fā)生了質(zhì)變，系統(tǒng)發(fā)生了P分岔。

根據(jù)概率密度的多峰狀態(tài)，為了判斷系統(tǒng)在弱隨機(jī)擾動、一般隨機(jī)擾動下是否進(jìn)入了混沌運(yùn)動狀態(tài)，現(xiàn)繪制概率時差圖對系統(tǒng)的動力學(xué)特性進(jìn)行進(jìn)一步分析，因?yàn)槎獧C(jī)翼在一定來流馬赫數(shù)下，受一定強(qiáng)度的隨機(jī)激勵影響發(fā)生分岔、進(jìn)入混沌運(yùn)動狀態(tài)時間非常短，因此繪制概率時差圖的時延也較小，本文選取0.01 s。圖4為顫振前區(qū)，在弱隨機(jī)擾動、一般隨機(jī)擾動下的概率時差圖；圖5為顫振后區(qū)， 在弱隨機(jī)擾動、一般隨機(jī)擾動下的概率時差圖。

從圖4、圖5可以發(fā)現(xiàn)，在顫振前區(qū)和顫振后區(qū)，系統(tǒng)受弱隨機(jī)擾動強(qiáng)度和一般隨機(jī)擾動強(qiáng)度下，概率時差圖均出現(xiàn)了典型的康托集合效應(yīng)，即在某一位置附近出現(xiàn)頻率明顯高于其他位置。結(jié)合其二維概率密度呈現(xiàn)多峰形狀，因此系統(tǒng)發(fā)生了Kapitaniak定義下的隨機(jī)混沌。

(a) D=0.01

(b) D=0.1圖4 不同擾動強(qiáng)度下顫振前區(qū)的概率時差圖(Ma=3)Fig.4 Probability time difference diagrams with various disturbance intensities(Ma=3)

(a) D=0.01

(b) D=0.1 圖5 不同擾動強(qiáng)度下顫振后區(qū)的概率時差圖(Ma=5)Fig.5 Probability time difference diagrams with various disturbance intensities(Ma=5)

由于強(qiáng)隨機(jī)擾動下二維聯(lián)合概率密度函數(shù)為單峰，根據(jù)Kapitaniak定義，此時系統(tǒng)不會發(fā)生Kapitaniak定義下的隨機(jī)混沌。

綜上可得，不管在顫振前區(qū)還是顫振后區(qū)，受外部隨機(jī)激勵的影響系統(tǒng)動態(tài)特性都會發(fā)生突變、發(fā)生分岔，甚至進(jìn)入混沌運(yùn)動狀態(tài)。但也會在一定隨機(jī)擾動強(qiáng)度下從混沌運(yùn)動狀態(tài)突變回周期運(yùn)動。因此，在進(jìn)行二元機(jī)翼氣動彈性特性研究時，要考慮不同飛行環(huán)境狀態(tài)的影響，避免系統(tǒng)發(fā)生復(fù)雜動力學(xué)行為。

5 數(shù)值計(jì)算驗(yàn)證

為了驗(yàn)證分析結(jié)果的有效性，現(xiàn)采用數(shù)值方法對分析結(jié)果進(jìn)行驗(yàn)證。首先采用Monte Carlo模擬法對隨機(jī)激勵進(jìn)行模擬，然后利用與前文一致的參數(shù)取值，在MATLAB平臺上對狀態(tài)方程(4)進(jìn)行數(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 不同擾動強(qiáng)度下顫振前區(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 不同擾動強(qiáng)度下顫振后區(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可以看出，在弱隨機(jī)擾動、一般隨機(jī)擾動強(qiáng)度作用下，無論是顫振前區(qū)還是顫振后區(qū)，系統(tǒng)的俯仰角隨時間均呈無規(guī)律變化，同時龐加萊截面由分形結(jié)構(gòu)的大量散點(diǎn)組成，且系統(tǒng)在Ma=3、Ma=5時最大Lyapunov指數(shù)均為正值，說明系統(tǒng)進(jìn)入了混沌運(yùn)動狀態(tài)；而在隨機(jī)擾動強(qiáng)度為0.5時，系統(tǒng)在Ma=3、Ma=5時最大Lyapunov指數(shù)均為負(fù)值，說明系統(tǒng)未進(jìn)入混沌運(yùn)動狀態(tài)，驗(yàn)證了前文中Kapitaniak定義下的隨機(jī)混沌分析結(jié)果，說明了本文近似解析定性分析方法的準(zhǔn)確性和有效性。

6 結(jié)論

本文通過累積量截?cái)喾?、非高斯截?cái)喾癊dgeworth展式等獲得系統(tǒng)的聯(lián)合概率密度函數(shù)，在Kapitaniak定義下研究了二元機(jī)翼的隨機(jī)混沌特性，并通過數(shù)值方法對計(jì)算結(jié)果進(jìn)行驗(yàn)證。主要研究結(jié)果如下：

1)采用第二類拉氏方程建立二元機(jī)翼隨機(jī)非線性動力學(xué)方程，并通過三階活塞理論推導(dǎo)二元機(jī)翼的非線性氣動力和氣動力矩。通過中心流形定理對系統(tǒng)進(jìn)行降維，將2自由度下4維非線性系統(tǒng)的動力學(xué)方程成功降為2維，使系統(tǒng)的隨機(jī)非線性動力學(xué)特性研究得到簡化。

2)采用累積量截?cái)喾ê头歉咚菇財(cái)喾ǐ@得系統(tǒng)二維聯(lián)合概率密度函數(shù)，利用Kapitaniak定義判斷系統(tǒng)的隨機(jī)混沌特性。研究發(fā)現(xiàn)，當(dāng)隨機(jī)擾動強(qiáng)度為0.01和0.1時，系統(tǒng)在顫振前區(qū)和顫振后區(qū)均進(jìn)入了混沌運(yùn)動狀態(tài)；在隨機(jī)擾動強(qiáng)度為0.5時，系統(tǒng)沒有進(jìn)入混沌運(yùn)動狀態(tài)。這說明外部隨機(jī)擾動對系統(tǒng)具有很大影響。因此，在二元機(jī)翼氣動彈性特性控制時，必須充分考慮外部隨機(jī)激勵強(qiáng)度的影響，以增加系統(tǒng)的穩(wěn)定性和安全性。

3)通過時域響應(yīng)曲線、龐加萊截面圖和最大Lyapunov指數(shù)等數(shù)值方法對計(jì)算結(jié)果進(jìn)行驗(yàn)證，兩種具有較好的一致性，說明本文研究具有較高的研究精度。