大展弦比NACA0012彈性平直機翼顫振特性研究

余秋陽, 李興華, 李新濤, 徐勝金

(1.清華大學(xué) 航天航空學(xué)院, 北京 100083; 2.哈爾濱理工大學(xué) 理學(xué)院, 黑龍江 哈爾濱 150080)

高空長航時飛行器因可執(zhí)行大氣監(jiān)測、通信中繼、高空偵察、反導(dǎo)防御等多種軍用或民用任務(wù),受到廣泛關(guān)注。為了提高氣動效率,在追求結(jié)構(gòu)輕量化的同時,這類飛行器還普遍采用大展弦比的機翼設(shè)計[1]。展弦比的增大會導(dǎo)致機翼結(jié)構(gòu)剛度的降低,在氣動載荷的作用下會出現(xiàn)較大變形和振動,變形和振動反過來又影響機翼表面流動,產(chǎn)生復(fù)雜的氣動彈性現(xiàn)象,如靜發(fā)散和顫振[2-4]。

對于大展弦比機翼顫振問題,現(xiàn)有的時域分析方法普遍采用分區(qū)思想,將機翼的非定常氣動力計算和結(jié)構(gòu)動力學(xué)分析在時間層面上交替進行,同時在機翼表面上進行數(shù)據(jù)傳遞。在每一個時間步內(nèi),首先求解作用在機翼上的氣動力載荷,然后求解機翼的變形和振動響應(yīng)并更新機翼的形狀,為下一時刻的氣動力計算做準備。

大展弦比機翼的非定常氣動力計算方法可以分為片條理論、面元法和CFD數(shù)值模擬三大類。片條理論的主要思想是將大展弦比機翼沿展向劃分為若干片條單元,采用Theodoren模型[5]、ONERA模型[6]以及有限狀態(tài)理論[7]等二維翼型氣動理論求解每個片條的氣動力,根據(jù)展向片條之間的氣動干擾進行三維修正,由此獲得整塊機翼受到的氣動力載荷[8]。渦格法[9]是一種具有代表性的面元法,該方法由勢流方程導(dǎo)出,將機翼離散成沿弦向和展向分布的附著渦,并在機翼后緣布置自由渦,利用翼面無穿透條件求解渦強,積分得到機翼的氣動特性[10]。隨著計算機技術(shù)的發(fā)展,機翼顫振時域分析逐漸采用以Euler方程或N-S方程為基礎(chǔ)的CFD數(shù)值模擬方法來計算非定常氣動力[11-12],這類方法直接求解流動基本方程,除了氣動力以外還能獲得瞬時流場信息。由此可得,目前常用的大展弦比機翼的氣動力求解方法大多基于勢流方程或Euler方程,忽略了氣體黏性對氣動力和流場的影響,應(yīng)用到顫振時域分析時可能無法得到合理的計算結(jié)果。

在大展弦比機翼顫振時域分析過程中,研究人員常用梁模型來求解機翼的結(jié)構(gòu)動力學(xué)響應(yīng),該方法根據(jù)機翼的結(jié)構(gòu)外形特征,將機翼抽象為非線性細長梁,再利用有限元方法進行數(shù)值求解,如本征梁模型[13]和梁單元模型[14]等。Hodges-Dowell梁[13]是一種典型的本征梁模型,該模型直接在變形后的梁坐標系下建立動力學(xué)變量特征,極大地減少了坐標的轉(zhuǎn)換。但是,除了梁的速度以外,該模型還需要預(yù)先求解梁的內(nèi)力,未知量個數(shù)大幅增加。共軛旋轉(zhuǎn)(CR)法[14]是最常用的梁單元模型之一,該方法的核心思想是將梁的變形分解為剛體的平動和轉(zhuǎn)動,以及單元坐標系內(nèi)的線彈性變形,動力學(xué)平衡方程在單元坐標系內(nèi)建立,通過坐標轉(zhuǎn)換矩陣構(gòu)造總體切線剛度矩陣[15]。CR法將梁變形本身的位移作為變量,物理意義明確,缺點是計算過程中需要不斷更新切線剛度矩陣,計算效率有所降低。除了梁模型以外,另一類常用的結(jié)構(gòu)動力學(xué)分析方法是模態(tài)法[16],該方法以結(jié)構(gòu)的主要固有振型為廣義坐標系,基于模態(tài)坐標建立動力學(xué)方程。相比于梁模型,模態(tài)法的優(yōu)勢在于建立的方程形式更為簡單,并且廣義質(zhì)量和廣義剛度矩陣都是對角矩陣,計算效率也更高。

綜上所述,大展弦比機翼的顫振時域分析主要關(guān)注顫振速度和頻率[17-18],部分研究也有涉及根部攻角[19]、展弦比[20]、扭轉(zhuǎn)剛度[12]等參數(shù)對機翼顫振特性的影響。機翼的變形和振動響應(yīng)有助于認識和揭示顫振誘發(fā)機理,目前仍然缺乏詳細的研究。本文以展弦比為16的NACA0012彈性平直機翼為研究對象,結(jié)合三維雷諾平均N-S方程、S-A湍流模型[21]和機翼動力學(xué)方程,建立了機翼顫振時域分析方法。利用該方法對機翼的顫振特性、顫振過程中的變形和振動響應(yīng),以及瞬時流場進行數(shù)值模擬和分析,探究機翼發(fā)生顫振的原因,研究材料剛度和密度對機翼顫振速度和頻率的影響。

1 大展弦比機翼顫振時域分析方法

1.1 非定場氣動力計算方法

本文采用自研CFD軟件求解展弦比為16的NACA0012彈性平直機翼的非定常氣動力,其控制方程為積分形式的守恒型三維非定常雷諾平均N-S方程,直角坐標系下的形式為

(1)

式中,守恒向量項Q為 [ρρuρvρwe]T,ρ為氣體密度,(uvw)為直角坐標系下速度分量,e為單位質(zhì)量氣體總能量。F,G和H為3個方向的對流通量矢量項,Fv,Gv和Hv為3個方向的黏性通量矢量項。Ω表示控制體。對流通量項采用迎風(fēng)型的Roe格式進行離散,黏性通量項采用標準的二階中心格式進行離散,時間推進采用雙時間步法,為了封閉方程引入S-A湍流模型[21]。

1.2 機翼動力學(xué)方程

本文以展弦比為16的NACA0012彈性平直機翼為研究對象,通過模態(tài)法[16]求解機翼的變形和振動響應(yīng)

(2)

式中：Δr為機翼的變形矢量;N為所選取的固有振型階數(shù);Φi(x,y,z)為第i階模態(tài)的固有振型矢量;ξi(t)為第i階模態(tài)的廣義位移。

機翼動力學(xué)方程可用矩陣形式表示為

(3)

(4)

式中

采用基于預(yù)估-校正技術(shù)的四階雜交線性多步法[22]對(4)式進行時域推進求解,將(4)式展開為(5)式,Δt為計算時間步長。該方法通過氣動力插值技術(shù)將氣動力部分由隱式格式轉(zhuǎn)變?yōu)轱@式格式,實現(xiàn)了每個物理時間步內(nèi)僅需求解一次非定常流場,兼具隱式穩(wěn)定性好和顯式效率高的優(yōu)點。

(5)

通過有限元模態(tài)分析得到大展弦比機翼的結(jié)構(gòu)參數(shù)矩陣(M,G和K)和固有振型,通過CFD數(shù)值模擬得到氣動力(F)。由于CFD翼面網(wǎng)格點和模態(tài)分析控制點一般是不重合的,在每一個時間步內(nèi),通過結(jié)構(gòu)固有振型在翼面處進行數(shù)據(jù)傳遞,即壓力分布通過振型轉(zhuǎn)化為廣義氣動力,進而求解機翼動力學(xué)方程獲得廣義位移,廣義位移通過振型再轉(zhuǎn)化為翼面物理位移,作為CFD下一時刻求解的邊界條件。模態(tài)分析得到的結(jié)構(gòu)固有振型通過徑向基函數(shù)(RBF)插值法[23]轉(zhuǎn)換為CFD翼面網(wǎng)格點描述下的振型分布,并采用基于徑向基函數(shù)插值的網(wǎng)格變形技術(shù)[23]實現(xiàn)機翼變形后CFD網(wǎng)格的運動。

1.3 算例驗證

AGARD445.6機翼顫振特性風(fēng)洞實驗是由NASA蘭利研究中心在其跨聲速動態(tài)風(fēng)洞中完成的[24]。目前,該機翼模型已經(jīng)成為國際上顫振計算程序考核的標準算例[25]。

該機翼模型沿流向翼型為NACA65A004,展弦比為1.644,根梢比為0.659 2,1/4弦線后掠角為45°,來流攻角為0°。機翼模型采用桃花心木制成,顫振分析取前四階振動模態(tài),其中,第一階和第三階模態(tài)為彎曲模態(tài),第二階和第四階模態(tài)為扭轉(zhuǎn)模態(tài)。

圖1 AGARD445.6機翼的顫振速度和頻率

圖1給出了本文仿真得到的AGARD445.6機翼顫振速度和頻率隨馬赫數(shù)Ma的變化關(guān)系,并與實驗結(jié)果[24]和文獻結(jié)果[25]比較。可以看出,顫振速度和頻率均出現(xiàn)了明顯的凹坑現(xiàn)象。本文結(jié)果在亞聲速狀態(tài)(Ma=0.499,0.678,0.901和0.96)時和實驗結(jié)果吻合很好,在超聲速區(qū)(Ma=1.072和1.141)略高于實驗結(jié)果,但是仍然明顯優(yōu)于文獻[25]的結(jié)果,表明本文建立的機翼顫振時域分析方法準確可靠。

2 結(jié)果與分析

2.1 主要振動模態(tài)

展弦比為16的NACA0012彈性平直機翼顫振是本文研究的重點。機翼模型的幾何參數(shù)、材料參數(shù)和飛行環(huán)境如表1所示。顫振分析取前五階振動模態(tài),表2給出了模態(tài)名稱和固有頻率,機翼的第一、二和五階模態(tài)為垂向彎曲模態(tài),第三階為扭轉(zhuǎn)模態(tài),第四階為弦向彎曲模態(tài)。

表1 大展弦比NACA0012彈性平直機翼模型參數(shù)

表2 大展弦比NACA0012彈性平直機翼前五階振動模態(tài)的固有頻率

2.2 顫振速度和頻率

不同來流速度U∞時大展弦比NACA0012彈性平直機翼的廣義位移時間響應(yīng)曲線如圖2所示?？梢钥闯?第一和二階振動模態(tài)廣義位移(ξ1和ξ2)在U∞為20 m/s和28.5 m/s時出現(xiàn)了一定幅度的靜漂移,說明由于展弦比較大,且材質(zhì)輕柔,機翼發(fā)生了靜彎曲變形。當U∞為20 m/s時,各階模態(tài)廣義位移響應(yīng)均表現(xiàn)為振動衰減。當U∞提高至28.5 m/s時,廣義位移均表現(xiàn)為近似等幅響應(yīng)。進一步提高U∞至40 m/s時,廣義位移表現(xiàn)出明顯的振動發(fā)散趨勢。綜合以上計算結(jié)果,可以得出,機翼在U∞為28.5 m/s時處于臨界穩(wěn)定狀態(tài),即顫振速度Uf約為28.5 m/s。

利用結(jié)構(gòu)固有振型將機翼的廣義位移轉(zhuǎn)化為翼面物理位移。圖3a)給出了顫振速度為28.5 m/s)時大展弦比NACA0012彈性平直機翼翼尖位置的垂向彎曲位移響應(yīng),圖3b)則給出了對應(yīng)的FFT頻譜?？梢钥闯?位移響應(yīng)主頻率為22 rad/s。因此,可以得出,機翼的顫振頻率ωf為22 rad/s。

圖4給出了來流速度U∞=40 m/s時機翼的整體變形特征以及z=15 m展向位置的壓力云圖。可以看出,當t=0.00 s時,機翼未發(fā)生變形,壓力關(guān)于機翼弦向平面對稱,從前緣到尾緣,壓力先減小然后增大。隨著時間t的推進,機翼的振動幅值逐漸增大,變形愈發(fā)明顯,其局部剖面的形狀和扭轉(zhuǎn)角都發(fā)生變化,機翼背風(fēng)面附近的負壓區(qū)也逐漸增大。

圖2 不同來流速度U∞時大展弦比NACA0012彈性平直機翼的廣義位移時間響應(yīng)曲線

圖3 U∞=28.5 m/s時大展弦比NACA0012彈性平直機翼翼尖位置的垂向彎曲位移響應(yīng)及其FFT頻譜

圖4 U∞=40 m/s時大展弦比NACA0012彈性平直機翼的整體變形特征以及z=15 m展向位置的壓力云圖

2.3 顫振誘發(fā)機理研究

頻率重合理論的顫振機理為:機翼2個振動模態(tài)的頻率隨著來流速度的增大而相互接近,這時這兩階模態(tài)的耦合作用增強,最終引發(fā)顫振[16,26]。通過比較和分析大展弦比NACA0012彈性平直機翼的顫振頻率(22 rad/s)和各階振動模態(tài)固有頻率,可以看出,第一階和第五階模態(tài)頻率均與顫振頻率相差較大,顫振頻率基本介于第二階和第三階模態(tài)頻率之間,因此,很有可能是這兩階模態(tài)耦合導(dǎo)致了機翼顫振,這也符合經(jīng)典顫振的誘發(fā)機理,即彎扭耦合作用[26]。然而,第三階和第四階模態(tài)頻率十分接近,第二～四階模態(tài)對機翼顫振的影響仍不明確,需要作進一步詳細分析。

將第二～四階振動模態(tài)分別固定,重新仿真得到顫振速度時機翼其他模態(tài)廣義位移的時間響應(yīng)曲線,如圖5所示。圖6給出了第四階模態(tài)固定時機翼翼尖位置的垂向彎曲位移響應(yīng)及其FFT頻譜?？梢钥闯?當?shù)诙A或第三階模態(tài)固定時,廣義位移均振動衰減,意味著顫振消失。當?shù)谒碾A模態(tài)固定時,廣義位移和翼尖位移均呈現(xiàn)近似等幅響應(yīng),機翼

圖5 不同振動模態(tài)固定時U∞=28. 5 m/s時大展弦比NACA0012彈性平直機翼其他模態(tài)廣義位移時間響應(yīng)曲線

圖6 第四階模態(tài)固定時U∞=28.5 m/s時大展弦比NACA0012彈性平直機翼翼尖位置的垂向彎曲位移響應(yīng)及其FFT頻譜

顫振速度和頻率沒有變化。以上結(jié)果說明,第二階和第三階模態(tài)的耦合作用是機翼發(fā)生顫振的根本原因。

2.4 材料剛度對顫振速度和頻率的影響

材料剛度和密度是機翼設(shè)計的重要參數(shù)。本節(jié)研究材料剛度對大展弦比NACA0012彈性平直機翼顫振速度和頻率的影響。將機翼的垂向彎曲剛度、弦向彎曲剛度和扭轉(zhuǎn)剛度等比例同時變化,選取了6組剛度值,如表3所示,其中RS為剛度比,表征機翼材料剛度的變化?？梢钥闯?最小的一組為原剛度的25%,即RS為0.25,最大的一組剛度的RS為1.5。相應(yīng)的,圖7給出了這6組剛度下機翼前五階振動模態(tài)的固有頻率,隨著剛度的增大,模態(tài)固有頻率也在逐漸升高。

圖8給出了不同剛度比下大展弦比NACA0012彈性平直機翼的顫振速度和頻率?？梢钥闯?對于最小的一組剛度,其顫振速度比原剛度時減小了約50%,顫振頻率減小了約49%。對于最大的一組剛度,其顫振速度和頻率分別提高了約23%和20%。這說明了隨著材料剛度的增大,機翼的顫振速度和頻率都會提高。

表3 大展弦比NACA0012彈性平直機翼的不同材料剛度

圖7 不同剛度比下大展弦比NACA0012彈性圖8 不同剛度比下大展弦比NACA0012彈性 平直機翼前五階振動模態(tài)的固有頻率 平直機翼顫振速度和頻率

2.5 材料密度對顫振速度和頻率的影響

本節(jié)研究材料密度對大展弦比NACA0012彈性平直機翼顫振速度和頻率的影響。定義ρs/ρ為密度比,用于表征機翼材料密度ρs的變化。圖9給出了不同密度比下機翼前五階振動模態(tài)的固有頻率,可以看出,最小的ρs為原材料密度的25%,最大的為原來的2.5倍。隨著ρs/ρ的增大,模態(tài)的固有頻率逐漸降低。

圖9 不同密度比下大展弦比NACA0012彈性圖10 不同密度比下大展弦比NACA0012彈性 平直機翼前五階振動模態(tài)的固有頻率 平直機翼顫振速度和頻率

圖10給出了不同密度比下大展弦比NACA0012彈性平直機翼的顫振速度和頻率?？梢钥闯?對于最小的材料密度,原ρs時的顫振速度和頻率分別比其降低了約16%和51%。對于最大的ρs,其顫振速度和頻率分別比原ρs下降低了約2%和39%。說明了隨著材料密度的增大,顫振速度會先明顯降低然后趨于不變,顫振頻率會一直降低。

3 結(jié) 論

本文采用以三維雷諾平均N-S方程和S-A湍流模型為基礎(chǔ)的CFD數(shù)值模擬方法來求解展弦比為16的NACA0012彈性平直機翼的非定常氣動力,結(jié)合機翼動力學(xué)方程,建立了機翼顫振時域分析方法,進而對機翼的顫振特性、振動響應(yīng)、整體變形特征和瞬時流場特征進行了數(shù)值模擬研究,主要得到以下結(jié)論:

1) 由于展弦比較大,且材質(zhì)輕柔,機翼會發(fā)生一定幅度的靜彎曲變形。機翼顫振由第二階(垂向彎曲)和第三階(扭轉(zhuǎn))振動模態(tài)的耦合作用所致,顫振頻率也介于這兩階模態(tài)的固有頻率之間。

2) 通過對具有不同材料剛度和密度的機翼進行數(shù)值模擬,對比其顫振速度和頻率,發(fā)現(xiàn)適當?shù)卦龃蟛牧蟿偠然蛘邷p小材料密度,可以提高大展弦比機翼的顫振速度和頻率。