基于Peters-ONERA模型的失速顫振特性研究

梁佳驊, 白俊強, 李國俊

(西北工業(yè)大學 航空學院, 陜西 西安 710072)

失速顫振是一種氣動彈性失穩(wěn)現(xiàn)象,常見于高空長航時(HALE)飛行器的大柔性機翼、直升機的槳葉以及渦輪機的葉片上[1]。失速顫振現(xiàn)象不僅會降低飛行器的氣動效率,而且會導致結構的疲勞甚至破壞。由于失速顫振常常伴隨著以流動分離為特征的動態(tài)失速現(xiàn)象,因此具有很強的氣動非線性[2]?？紤]到失速顫振問題存在著復雜的非線性因素以及研究人員對失速顫振機理缺乏深入的理解,因此有必要對失速顫振問題開展更加深入的研究工作,以便為失速顫振主動控制提供理論基礎和技術支撐。

在當今的氣動彈性研究領域中,將計算流體力學(computational fluid dynamics,CFD)與計算結構力學(computational structure dynamics,CSD)耦合是進行非線性顫振分析的重要手段之一[3]。然而由于實際工程問題的復雜性,所構建的物理模型往往十分復雜,CFD/CSD方法存在耗時長、效率低等劣勢[4],而其中主要的計算代價來源于非線性氣動力的求解。這使得CFD/CSD方法難以廣泛適用于失速顫振的研究工作中。

為了簡化氣動力,Peters等人[5]基于Glauert分解方法發(fā)展了非定常氣動力的有限狀態(tài)(finite-state)理論和入流(induced-flow)理論。相比其他氣動力模型,該方法計算精度較高,工程實際應用性較強。此外,該方法具有時域表達形式,可以加入動態(tài)失速修正項來模擬非線性氣動力。當前針對二維翼型已經(jīng)提出了一些較為成熟的半經(jīng)驗失速模型,例如ONERA失速模型和Beddoes-Leishman失速模型等。文獻[2]總結了當前工程和研究中常用的失速模型,其中ONERA模型形式較為簡單,應用較為廣泛。ONERA模型通過適當?shù)木€性有理近似,將動失速氣動力表達為一個二階非線性常微分方程的形式。此外,ONERA模型以線化氣動力系數(shù)和靜態(tài)氣動力系數(shù)之差作為輸入量,可以描述氣動力的延遲效應和超調效應。

Peters等人[6]采用ONERA失速模型對直升機葉片的動態(tài)失速問題展開了研究。Mcalister等人[7]基于ONERA失速模型對二維翼型動態(tài)失速中的非定常氣動力響應進行預測,并與實驗值進行對比。大量研究結果表明ONERA失速模型可以較好地模擬動態(tài)失速中的非定常氣動力。Tang等人[8]通過引入小角度假設,將ONERA模型轉化到拉氏域中獲得線性非定常氣動力,并從頻域角度對失速顫振問題展開研究。研究結果表明氣動彈性系統(tǒng)在大攻角下的失穩(wěn)模態(tài)發(fā)生了轉換,但研究人員并未對此開展更加詳細的機理解釋。Laxman等人[9]基于ONERA失速模型對二元機翼的混沌響應特性進行研究。Beedy等人[10]基于諧波分解方法和Levenberg-Marquardt非線性最小二乘法,發(fā)展了一種擬合ONERA非線性參數(shù)的新方法,并基于此進行失速顫振邊界的預測。

在國內,劉湘寧等人[11]采用ONERA失速模型,結合諧波平衡法,研究了大展弦比復合材料機翼的顫振穩(wěn)定性。劉廷瑞[12]采用ONERA失速模型,基于H∞魯棒控制進行風輪機葉片的失速顫振抑制研究工作。任勇生等人[13]采用ONERA模型進行了復合材料薄壁梁的非線性氣彈穩(wěn)定性研究工作。孫智偉等人[14]首次綜合考慮了Peters氣動力模型和ONERA失速模型,提出了一種新的適用于失速顫振研究的氣動力模型,并結合魯棒控制開展了顫振抑制方面的研究工作。

上述大多是基于ONERA失速模型開展的失速顫振邊界預測和主動顫振抑制研究工作,鮮有研究工作者基于Peters-ONERA模型,對大攻角下的失速顫振臨界特性和不同初始攻角下氣動彈性系統(tǒng)的分岔現(xiàn)象展開研究。此外,基于Peters-ONERA模型對失速顫振中靜氣彈求解穩(wěn)定性的研究工作也很少見。本文基于Peters-ONERA氣動力模型,建立了系統(tǒng)的狀態(tài)空間方程。首先通過動態(tài)失速算例驗證了氣動力模型,并研究了亞松弛迭代對靜氣彈求解穩(wěn)定性的影響;其次通過頻域方法對大攻角下的失速顫振臨界特性展開研究;最后通過時域方法研究了失速顫振中的極限環(huán)振蕩和分岔現(xiàn)象,并研究了初始擾動對系統(tǒng)響應的影響。

1 氣動力建模

本節(jié)基于Peters氣動力模型和ONERA失速模型分別對線性氣動力和非線性氣動力進行建模。由于上述氣動力模型都具有時域表達形式,因此易于和結構運動方程耦合,形成完備的狀態(tài)空間形式。

1.1 Peters線性模型

圖1展示了一個具有兩自由度的典型二維翼型。其中Q為氣動中心,P為剛心,G為重心。結構坐標系參考原點取在剛心P處,浮沉位移h向下為正,俯仰位移θ抬頭為正。剛心P在翼型弦長中點后ab處,其中a為無量綱長度,b為半弦長。

圖1 典型兩自由度二維翼型

氣動彈性系統(tǒng)關于彈性軸的升力和力矩的線性部分如下[15]:

(1)

在方程(1)中,ρ代表大氣密度,U為無窮遠來流速度,α0代表翼型的初始迎角。在線性升力L0的表達式中,第一項為無環(huán)量升力,其依賴于翼型的速度與加速度;第二項為有環(huán)量升力,其通過λ0考慮了翼型后緣尾渦的影響。其中λ0表示平均入流誘導速度,Peters將其近似用有限個入流狀態(tài)量來表示[5]:

λ0=0.5bTλ

(2)

式中，λ滿足如下微分方程[15]:

(3)

1.2 ONERA失速模型

在ONERA模型中,Γn(n代表升力L或者力矩M)表示由于失速帶來的附加環(huán)量,可以通過二階常微分方程表示。需要指出的是,方程(4)的力矩參考點是翼型的四分之一弦長位置。

(4)

在方程(4)中,氣動力的非線性特性可以通過輸入值ΔCn(n代表升力L或者力矩M)得到。而ΔCn是關于攻角α的函數(shù),一般通過實驗方法或者CFD方法得到。半經(jīng)驗參數(shù)ξn,ωn和ηn(n代表升力L或者力矩M)分別代表了氣動力阻尼、頻率,以及氣動力延遲。參數(shù)可以表示成如下形式(n代表升力L或者力矩M):

ξn=ξ0+ξ2(ΔCn)2

ωn=ω0+ω2(ΔCn)2

ηn=η0+η2(ΔCn)2

(5)

為了得到方程(5)中各個系數(shù)的典型值,必須基于某個特定翼型的實驗結果或CFD結果進行參數(shù)辨識。Beedy等人[10]認為由非定常升力響應數(shù)據(jù)擬合得到的非線性參數(shù)可以用于非定常力矩的預測。如表1所示,給出了部分ONERA非線性參數(shù)的典型值[16],其中Ma代表無窮遠來流馬赫數(shù)。

表1 ONERA非線性參數(shù)的典型值

綜合以上討論,可以得到氣動彈性系統(tǒng)關于彈性軸處的總升力和力矩:

(6)

式中,Cl0和Cm0表示在零攻角下翼型的升力和力矩系數(shù),L0和M0為氣彈系統(tǒng)關于彈性軸升力和力矩的線性部分,具體形式可參照(1)式。

孫智偉[17]認為,忽略阻力不會對動力學特性的定性分析產(chǎn)生較大的影響;張健等人[18]指出,如果不需要精確地預測顫振邊界和系統(tǒng)極限環(huán)響應,可以忽略阻力因素。本文中沒有考慮阻力因素,可以在一定程度上簡化氣動力模型的復雜程度。

2 氣動彈性系統(tǒng)求解方法

對于圖 1所示的氣動彈性系統(tǒng),基于Lagrange方程可以得到結構運動方程的矩陣形式表達如下[19]:

(7)

式中,M代表質量矩陣,G代表阻尼矩陣,K代表剛度矩陣,δ代表廣義位移,f代表作用在系統(tǒng)上的廣義力。

為了得到氣動力控制方程與結構運動方程的全耦合形式,便于失速顫振問題的分析,將廣義氣動力f寫成如下的矩陣表達形式[14]:

(8)

式中，δ代表廣義位移。

(9)

在方程(9)中,矩陣H0和H1的階數(shù)和λ的維數(shù)有很大關系。為了保證計算的精度,λ的維數(shù)可以取4～8[5],但λ的維數(shù)取得過大可能會導致數(shù)值不穩(wěn)定[20]。文獻[14]建議取維數(shù)n=6。

方程(9)即為氣動彈性系統(tǒng)狀態(tài)空間方程的完整形式。從時域分析的角度來說,此問題轉化為一階非線性常微分方程組的求解問題,可以選用多種求解技術,本文采用歐拉預估-校正方法進行時域推進求解;從頻域分析的角度來說,此問題轉化為矩陣特征值的求解問題,本文通過求解矩陣H1的特征根得到了系統(tǒng)的特征根軌跡。值得注意的是,由于矩陣H1是一個非線性矩陣,因此為了進行根軌跡分析,需要對矩陣H1進行線性化處理。本文采用雅克比矩陣法對H1進行線性化處理。

3 算例分析

3.1 動態(tài)失速算例驗證

Mcalister等人[7]指出,在減縮頻率k小于0.15的情況下,ONERA失速模型對非線性氣動力的模擬具有一定的可信度。對于更高的減縮頻率,工程上一般不給予考慮,因為實際機翼也很難達到如此之高的振動頻率[17]。

圖2 升力系數(shù)遲滯曲線對比圖(α=10°+15°sinωt)

本文選用減縮頻率k分別為 0.1和0.05的2個算例進行研究,采用歐拉預估-校正方法對氣動力進行時域推進求解。圖2給出了數(shù)值解和實驗值對比圖。由圖2可知,當減縮頻率為0.1和0.05時,數(shù)值模擬結果在上、下行與實驗值吻合較好。此外,該模型對失速攻角位置的預測能力具有一定的可信度。這與Mcalister等人得到的結論是一致的。然而,該模型未能很好地捕捉到實驗結果中的尖點區(qū)域。這是因為基于ONERA失速模型得到的結果在數(shù)學上是嚴格的二階導數(shù)連續(xù)解,其無法模擬突跳的尖點等不可導區(qū)域。但從總體上來說,當引入ONERA失速模型后,失速后氣動力的變化趨勢和變化現(xiàn)象都能夠被比較清晰地捕捉到。

3.2 考慮亞松弛迭代的靜氣彈問題研究

對于非線性顫振問題,無論是在頻域內進行特征根軌跡分析,還是在時域內進行推進求解,首先應得到系統(tǒng)在某一條件下的初始穩(wěn)態(tài)解,即靜氣彈求解問題。本小節(jié)通過一個算例,對比采用不同迭代方法得到的靜氣彈結果,以說明亞松弛迭代在靜氣彈求解中的意義,并采用亞松弛迭代法,給出不同初始攻角下,系統(tǒng)的靜氣彈解隨來流速度的變化曲線。

對于圖1的氣動彈性系統(tǒng)式(9)給出了系統(tǒng)的狀態(tài)空間形式。為了得到系統(tǒng)的靜氣動彈性控制方程,只需要將所有關于時間的導數(shù)項全部置零,最終得到系統(tǒng)的靜氣動彈性控制方程如下[14]:

(K-A3)·δ=F2Γ+C0

D2Γ+D0=0

(10)

式中,K為剛度陣,C0為零攻角氣動力系數(shù)矩陣,δ為廣義位移向量,Γ為失速氣動力向量。需要注意的是,矩陣D0和D2是關于廣義位移向量δ的復雜分段函數(shù)。如果在(10)式的第二式中直接將Γ解出,并帶入第一式中去求解δ,算式將非常復雜。故本文采用迭代方法對兩式進行依次間接求解。

以文獻[8]中的算例為例,如圖3所示展示出靜態(tài)氣動力系數(shù)隨攻角變化曲線。

圖3 靜態(tài)氣動力系數(shù)曲線

由圖3可知,翼型的靜失速攻角為13°。靜態(tài)升力和力矩系數(shù)曲線的斜率在此處都有突變。圖 4展示出采用簡單迭代法和亞松弛迭代法得到的結果,其中亞松弛因子取0.4,翼型的初始攻角為18°,圖中虛線代表了失速攻角的位置。

電力電纜設備故障的類型比較多，在當前技術應用中需要做好具體技術分析工作，結合系統(tǒng)本身流程要求和探測要求等，在故障分析的過程中采用循序漸進的方式進行處理。最重要的是強化系統(tǒng)日常維護處理，保證電力電纜設備發(fā)揮穩(wěn)定性，為行業(yè)進步奠定基礎。

圖4 靜氣彈收斂歷程對比圖

當采用簡單迭代法,靜氣動彈性解在來流速度為14.5 m/s和15.0 m/s時,均有不同程度的振蕩。這是由于當來流速度為14.5 m/s和15.0 m/s時,解在靜失速攻角附近反復迭代。而一旦涉及非線性迭代,差值過大將導致迭代過程的不穩(wěn)定。當采用亞松弛迭代法時,系統(tǒng)解的穩(wěn)定性提高,收斂曲線變得光滑。在來流速度為14.5 m/s和15.0 m/s時,解收斂到一個穩(wěn)定值。然而,當來流速度為14.0 m/s時可以看出,當采用簡單迭代法時,解在第三步已經(jīng)收斂;當采用亞松弛迭代法時,解在第十步才基本收斂。采用亞松弛迭代法固然可以提高系統(tǒng)解的穩(wěn)定性,但同時也會降低系統(tǒng)解的收斂速度,需要權衡考慮。

采用亞松弛迭代技術,通過數(shù)值模擬預測在不同初始攻角下,系統(tǒng)靜氣彈解隨來流速度的變化情況,亞松弛因子取0.4。如圖5所示展示出計算結果,其中點劃線代表了失速攻角的位置:

圖5 系統(tǒng)靜氣彈解變化曲線

由圖5可知,當初始攻角一定時,來流速度越大,氣動載荷越大,系統(tǒng)的靜氣彈解數(shù)值的絕對值也就越大;在來流速度一定時,初始攻角越大,曲線的斜率越大,表明靜氣彈解的發(fā)散速度越快。值得注意的是:每條曲線在13°處的斜率都是不連續(xù)的,斜率有明顯下降。這是因為:當俯仰位移小于靜失速攻角時,系統(tǒng)中的氣動力只有線性氣動力的參與;當俯仰位移大于靜失速攻角時,非線性氣動力將參與迭代過程,這將導致曲線斜率的突變。由于非線性氣動力的參與,當攻角增加時,氣動力增加趨勢有所減緩,從而減緩了系統(tǒng)靜氣彈解的發(fā)散趨勢。

3.3 失速顫振的顫振臨界特性研究

圖6給出了系統(tǒng)在初始攻角為0°時的頻域響應結果。其中a分支代表俯仰模態(tài);b分支代表沉浮模態(tài);c分支代表線性氣動力項;d分支代表尾渦脫落項;e分支代表了非線性氣動力項。

由圖6a)中的根軌跡圖可知,隨著來流速度的增加,系統(tǒng)的俯仰模態(tài)分支a穿過虛軸,而A點對應的來流速度為16.9 m/s。此時系統(tǒng)由于俯仰模態(tài)失穩(wěn)而發(fā)生顫振。從圖6b)可以看出,在顫振臨界點附近區(qū)域,系統(tǒng)的沉浮模態(tài)和俯仰模態(tài)頻率在不斷靠近,具有典型的經(jīng)典顫振特性。此外,線性氣動力c在B點的失穩(wěn)代表著系統(tǒng)的靜氣動彈性發(fā)散,對應的來流速度為22.9 m/s,即系統(tǒng)靜氣彈發(fā)散速度為22.9 m/s。

下面研究大攻角下的系統(tǒng)響應。如圖7～9所示,展示了初始攻角為25°,32°和35°時的系統(tǒng)頻域響應結果。從圖7～9可以看出,在大攻角情況下,系統(tǒng)顫振特性和經(jīng)典顫振特性是不相同的。

圖6 初始攻角為0°系統(tǒng)頻域響應結果 圖7 初始攻角為25°系統(tǒng)頻域響應結果

圖8 初始攻角為32°系統(tǒng)頻域響應結果 圖9 初始攻角為35°系統(tǒng)頻域響應結果

當翼型的初始攻角為25°時,由圖7a)的根軌跡圖可知,在顫振邊界附近,系統(tǒng)的非線性氣動力模態(tài)分支e有失穩(wěn)的趨勢,系統(tǒng)的俯仰模態(tài)分支a穿過虛軸失穩(wěn)。從圖7b)可以看出,在顫振臨界點附近,系統(tǒng)的俯仰模態(tài)分支a和非線性氣動力模態(tài)分支e頻率在不斷靠近,表明系統(tǒng)的俯仰模態(tài)和非線性氣動力模態(tài)有耦合的趨勢。

當翼型的初始攻角為32°時,由圖8a)的根軌跡圖可知,隨著氣動力非線性的增強,系統(tǒng)的俯仰模態(tài)分支a和沉浮模態(tài)分支b相互吸引,根軌跡有相互靠近的趨勢,但此時系統(tǒng)依然發(fā)生的是由于俯仰模態(tài)失穩(wěn)主導的失速顫振。

當翼型的初始攻角為35°時,由圖9a)的根軌跡圖可知,系統(tǒng)的俯仰模態(tài)分支a和沉浮模態(tài)分支b繼續(xù)相互吸引靠近,最終發(fā)生模態(tài)分支交換。此時系統(tǒng)從俯仰模態(tài)失穩(wěn)轉變?yōu)槌粮∧B(tài)失穩(wěn)。由圖9b)可以看出, 在顫振臨界點附近,系統(tǒng)的沉浮模態(tài)分支b和非線性氣動力模態(tài)分支e頻率在不斷靠近,這表明系統(tǒng)的沉浮模態(tài)和非線性氣動力模態(tài)有耦合的趨勢。

綜上所述,在翼型大攻角的失速顫振問題中,非線性氣動力模態(tài)與結構模態(tài)的耦合作用可能導致結構模態(tài)失穩(wěn),從而誘發(fā)系統(tǒng)的單自由度顫振。

3.4 失速顫振的極限環(huán)振蕩和分岔現(xiàn)象研究

本小節(jié)依然采用3.3小節(jié)的結構參數(shù),進行失速顫振的極限環(huán)振蕩和分岔現(xiàn)象的研究。圖10給出了系統(tǒng)在初始攻角為14°,來流速度為10.1 m/s和10.7 m/s時的時域響應結果。

圖10 不同來流速度下的系統(tǒng)響應

對于經(jīng)典顫振問題,當來流速度超過顫振臨界速度時,系統(tǒng)響應幅值為無限大值,此時系統(tǒng)發(fā)生了動氣彈發(fā)散現(xiàn)象。對于失速顫振問題,當來流速度超過顫振臨界速度時,系統(tǒng)響應幅值為有限值,此時系統(tǒng)處于極限環(huán)振蕩狀態(tài),如圖10b)所示。極限環(huán)振蕩是非線性系統(tǒng)的典型特性。當飛行器結構發(fā)生極限環(huán)振蕩時,雖然不會導致飛機的直接解體和破壞,但會帶來較為嚴重的飛機結構疲勞問題。

圖11 不同初始攻角下的系統(tǒng)分岔曲線

圖11展示了不同初始攻角下,系統(tǒng)俯仰模態(tài)幅值隨來流速度變化的分岔曲線,其中取初始靜氣彈平衡位置為零位置。每條分岔曲線都以小振幅的顫振臨界狀態(tài)為起始端,動氣彈發(fā)散狀態(tài)為結束端。當翼型的初始攻角為6°時,俯仰模態(tài)振幅在來流速度超過顫振臨界速度后急劇增大,出現(xiàn)超臨界霍夫分岔,幅值范圍較小;當初始攻角增加到8°和14°時,俯仰模態(tài)振幅在來流速度超過顫振臨界速度后增速減緩,同時模態(tài)幅值范圍變大;當初始攻角增加到18°時,俯仰模態(tài)振幅在來流速度超過顫振臨界速度后有一個小的階躍,此后模態(tài)幅值增速略有增加。當初始攻角增加到28°和30°時,俯仰模態(tài)振幅的變化趨勢和小攻角時十分相似。

綜上所述,初始攻角的改變會顯著影響系統(tǒng)俯仰模態(tài)的振幅分岔曲線特性。

3.5 初始擾動對系統(tǒng)響應的影響研究

當考慮失速顫振時,氣動彈性系統(tǒng)本身作為非線性系統(tǒng),其響應會受到初始擾動的影響。為此,有必要研究初始擾動對氣動彈性系統(tǒng)響應的影響。在本算例中,依然采用3.3小節(jié)的結構參數(shù),通過改變初始擾動值得到系統(tǒng)的響應情況。在數(shù)值模擬時,取初始靜氣彈平衡位置為零位置,初始攻角為25°。如圖12所示給出了系統(tǒng)在不同初始擾動下俯仰模態(tài)幅值隨來流速度的變化曲線。其中圖12a)為全局圖,圖12b)為圖12a)中圓圈區(qū)域的局部放大圖。

圖12 不同初始速度擾動下的系統(tǒng)響應

在初始擾動為2×10-3時,系統(tǒng)受到的外界擾動較小,當來流速度在6.3 m/s到6.8 m/s范圍內時,系統(tǒng)做小幅值的極限環(huán)振蕩,而當來流速度在6.8 m/s到7.3 m/s范圍內,系統(tǒng)響應是收斂的,模態(tài)幅值為0;當來流速度超過7.3 m/s時,模態(tài)幅值又迅速大幅度增加。當初始擾動為0.1,0.05和0.02時,系統(tǒng)受到的外界擾動有所增大,俯仰模態(tài)幅值曲線幾乎是完全重合的。在來流速度為6.3～6.8 m/s范圍內時,系統(tǒng)做小幅值的極限環(huán)振蕩。而當來流速度在6.75 m/s附近時,模態(tài)幅值迅速大幅度增加,此后沿著光滑的曲線逐漸上升,系統(tǒng)響應繼續(xù)持續(xù)著極限環(huán)振蕩狀態(tài),并未出現(xiàn)收斂現(xiàn)象。

綜上所述,當初始擾動在0.02到0.1范圍時,系統(tǒng)對擾動變化的敏感度較弱;而當初始擾動在2×10-3到0.02之間時,系統(tǒng)對擾動的變化非常敏感。當系統(tǒng)的初始擾動變大時,系統(tǒng)原先穩(wěn)定的狀態(tài)可能被激發(fā)為極限環(huán)振蕩狀態(tài)。

4 結 論

本文基于Peters-ONERA模型進行了氣動彈性系統(tǒng)的氣動力建模,并耦合結構運動方程,建立了氣動彈性系統(tǒng)的狀態(tài)空間方程,采用時域和頻域方法進行失速顫振特性的研究工作。研究結果表明:

1) 本文所采用的氣動力模型可以準確捕捉動態(tài)失速氣動力的主要特征。

2) 采用亞松弛迭代法可以有效地抑制靜氣彈解在迭代過程中的振蕩現(xiàn)象,增強靜氣彈求解的穩(wěn)定性。

3) 在大攻角的失速顫振問題中,非線性氣動力模態(tài)與結構模態(tài)的耦合作用可能導致此結構模態(tài)的失穩(wěn),從而誘發(fā)系統(tǒng)的單自由度顫振。

4) 在失速顫振中,初始攻角的改變會顯著影響系統(tǒng)的分岔特性。

5) 在不同的初始擾動范圍內,氣動彈性系統(tǒng)對擾動的敏感度是不同的。當系統(tǒng)的初始擾動變大時,系統(tǒng)原先穩(wěn)定的狀態(tài)可能被激發(fā)為極限環(huán)振蕩狀態(tài)。