超聲速二元機翼隨機混沌運動分析*

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

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

氣動彈性是飛行器設計中備受關注的重要問題，主要氣動彈性現象有顫振、抖振、動力響應等[1-2]。其中，動力響應是指彈性系統受到與自身系統無關的、隨時間任意變化的外界干擾力作用而發生的強迫振動，外界擾動可以是諧和的、周期的、脈沖的或隨機的[3]。在航空動力學領域中,隨機擾動普遍存在，它主要是由飛行器飛行過程中空氣熱力、風力和尾渦等隨機因素相互作用形成的大氣湍流引起的。從強度觀點來看，飛機結構可能在嚴重的湍流中由于超載而遭破壞，中等大小的湍流則是飛機結構疲勞損傷的主要來源[4]。因此，結合氣流擾動的機翼系統動力學特性研究能更好地反映實際工況。近年來，考慮外界擾動的機翼氣動彈性問題正受到越來越多學者的關注。

Poirel和Price[5-8]研究了二元機翼在湍流隨機擾動作用下的顫振問題，考慮了機翼結構線性和立方非線性情況，將湍流近似為高斯隨機過程，從概率密度和最大Lyapunov指數分析了二元機翼的隨機動力學行為，指出了隨機激勵下的顫振點較確定性系統顫振點提前，且與激勵的強度密切相關。文獻[9]采用能量隨機平均法，求解FPK(Fokker-planck-kolmogorov)方程，獲得隨機激勵下的一元機翼分岔點，并分析了系統的隨機Hopf分岔特性。文獻[10]利用攝動法獲得了二元機翼在非高斯色噪聲作用下的Lyapunov指數，分析了二元機翼在隨機激勵下的穩定性。文獻[11]采用隨機平均法得到二元機翼在寬帶噪聲激勵下的Lyapunov指數，探討了隨機噪聲譜密度對機翼穩定性的影響。文獻[12]采用隨機減量法和矩陣束法識別紊流激勵中的模態參數，提出了一種紊流激勵下的顫振邊界預測方法。可以看出，對隨機激勵下的機翼動力學特性研究主要集中在隨機顫振預測、隨機分岔特性研究中。

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

1 二元機翼隨機非線性動力學方程

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

采用第二類拉氏方程獲得受隨機擾動的二元機翼動力學方程為：

(1)

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

(2)

(3)

(4)

A0矩陣中各元素表達式如下：

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

其中

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

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

2 采用中心流形方法對系統進行降維

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

通過Routh-Hurwitz判據得到：Ma<3.985時，系統穩定，響應收斂到平衡點；3.9854.086時，系統發散，因此，Ma*=3.985為系統Hopf分岔點。

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

(5)

根據中心流形定理，得到中心流形函數

(6)

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

(7)

3 二維聯合概率密度函數求解

根據Kapitaniak對隨機混沌特性的定義：二維聯合概率密度分布具有多峰且概率時差圖具有典型的康托集合結構時，系統具有隨機混沌特性。為此首先對系統的二維聯合概率密度函數進行分析。由于方程(7)為非線性方程，故本文聯合利用累積量截斷法、非高斯截斷法對系統的二維聯合概率密度函數進行求解。方程(7)的伊藤隨機微分方程形式如下：

(8)

式中，D1、D1為隨機激勵強度。

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

(10)

(11)

4 二元機翼隨機混沌特性分析

現通過二維聯合概率密度函數曲線、概率時差圖對二元機翼的隨機動力學特性進行研究。由于機翼俯仰自由度的振動幅值比沉浮自由度的振動幅值大，故本文主要研究俯仰角-俯仰角速度的概率密度曲線。現將系統發生顫振時的來流馬赫數Ma=3.985作為分界線，以此將系統劃分為顫振前區、顫振后區，并研究來流馬赫數分別為Ma=3、Ma=5時系統顫振前區和顫振后區的動力學特性。

根據第3節求解二維概率密度函數的方法，獲得系統顫振前區和顫振后區在隨機擾動作用下的概率密度函數。隨機擾動強度根據文獻[17]添加，即弱隨機擾動取值0.01，一般隨機擾動取值0.1，強隨機擾動取值0.5。圖2是在顫振前區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時系統的二維概率密度分布圖；圖3是在顫振后區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可以看出：二維概率密度在顫振前區時，形狀在弱隨機擾動下為分離雙峰，一般隨機擾動下為相鄰雙峰，強隨機擾動下變為單峰；在顫振后區時，形狀在弱隨機擾動下為多峰，一般隨機擾動下變為雙峰，強隨機擾動下變為單峰。并且在顫振前區和顫振后區的二維聯合概率密度形狀在弱隨機擾動下不同，而在強隨機擾動下形狀相似，由此可知系統拓撲結構發生了質變，系統發生了P分岔。

根據概率密度的多峰狀態，為了判斷系統在弱隨機擾動、一般隨機擾動下是否進入了混沌運動狀態，現繪制概率時差圖對系統的動力學特性進行進一步分析，因為二元機翼在一定來流馬赫數下，受一定強度的隨機激勵影響發生分岔、進入混沌運動狀態時間非常短，因此繪制概率時差圖的時延也較小，本文選取0.01 s。圖4為顫振前區，在弱隨機擾動、一般隨機擾動下的概率時差圖；圖5為顫振后區， 在弱隨機擾動、一般隨機擾動下的概率時差圖。

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

(a) D=0.01

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

(a) D=0.01

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

由于強隨機擾動下二維聯合概率密度函數為單峰，根據Kapitaniak定義，此時系統不會發生Kapitaniak定義下的隨機混沌。

綜上可得，不管在顫振前區還是顫振后區，受外部隨機激勵的影響系統動態特性都會發生突變、發生分岔，甚至進入混沌運動狀態。但也會在一定隨機擾動強度下從混沌運動狀態突變回周期運動。因此，在進行二元機翼氣動彈性特性研究時，要考慮不同飛行環境狀態的影響，避免系統發生復雜動力學行為。

5 數值計算驗證

為了驗證分析結果的有效性，現采用數值方法對分析結果進行驗證。首先采用Monte Carlo模擬法對隨機激勵進行模擬，然后利用與前文一致的參數取值，在MATLAB平臺上對狀態方程(4)進行數值求解。不同工況下系統的時域響應曲線、龐加萊截面圖如圖6、圖7所示，系統的最大Lyapunov指數隨來流馬赫數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 不同擾動強度下顫振前區的系統響應曲線及龐加萊截面(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 不同擾動強度下顫振后區的系統響應曲線及龐加萊截面(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指數變化曲線圖Fig.8 The largest lyapunov exponent changes with Ma

從圖6～8可以看出，在弱隨機擾動、一般隨機擾動強度作用下，無論是顫振前區還是顫振后區，系統的俯仰角隨時間均呈無規律變化，同時龐加萊截面由分形結構的大量散點組成，且系統在Ma=3、Ma=5時最大Lyapunov指數均為正值，說明系統進入了混沌運動狀態；而在隨機擾動強度為0.5時，系統在Ma=3、Ma=5時最大Lyapunov指數均為負值，說明系統未進入混沌運動狀態，驗證了前文中Kapitaniak定義下的隨機混沌分析結果，說明了本文近似解析定性分析方法的準確性和有效性。

6 結論

本文通過累積量截斷法、非高斯截斷法及Edgeworth展式等獲得系統的聯合概率密度函數，在Kapitaniak定義下研究了二元機翼的隨機混沌特性，并通過數值方法對計算結果進行驗證。主要研究結果如下：

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

2)采用累積量截斷法和非高斯截斷法獲得系統二維聯合概率密度函數，利用Kapitaniak定義判斷系統的隨機混沌特性。研究發現，當隨機擾動強度為0.01和0.1時，系統在顫振前區和顫振后區均進入了混沌運動狀態；在隨機擾動強度為0.5時，系統沒有進入混沌運動狀態。這說明外部隨機擾動對系統具有很大影響。因此，在二元機翼氣動彈性特性控制時，必須充分考慮外部隨機激勵強度的影響，以增加系統的穩定性和安全性。

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