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

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

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

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

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

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

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

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

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

1.1 非定場氣動力計算方法

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

(1)

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

1.2 機翼動力學方程

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

(2)

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

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

(3)

(4)

式中

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

(5)

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

1.3 算例驗證

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

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

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

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

2 結果與分析

2.1 主要振動模態

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

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

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

2.2 顫振速度和頻率

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

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

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

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

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

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

2.3 顫振誘發機理研究

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

將第二～四階振動模態分別固定,重新仿真得到顫振速度時機翼其他模態廣義位移的時間響應曲線,如圖5所示。圖6給出了第四階模態固定時機翼翼尖位置的垂向彎曲位移響應及其FFT頻譜。可以看出,當第二階或第三階模態固定時,廣義位移均振動衰減,意味著顫振消失。當第四階模態固定時,廣義位移和翼尖位移均呈現近似等幅響應,機翼

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

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

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

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

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

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

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

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

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

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

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

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

3 結 論

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

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

2) 通過對具有不同材料剛度和密度的機翼進行數值模擬,對比其顫振速度和頻率,發現適當地增大材料剛度或者減小材料密度,可以提高大展弦比機翼的顫振速度和頻率。