計及熱效應(yīng)的小展弦比機(jī)翼氣彈穩(wěn)定性分析＊

趙 衛(wèi)， 韓曉林， 費慶國

（1.東南大學(xué)土木工程學(xué)院 南京，210096） （2.東南大學(xué)江蘇省工程力學(xué)分析重點實驗室 南京，210096）

引 言

隨著高／超聲速飛行器的發(fā)展，高／超聲速氣流引起氣動加熱產(chǎn)生的熱問題越來越受到重視。結(jié)構(gòu)在熱環(huán)境下的動力學(xué)問題是發(fā)展高／超聲速航天飛行器面臨的重要課題。20世紀(jì)50年代末至60年代初，美國對熱氣動彈性問題的研究出現(xiàn)了一個高潮，這些研究為之后航天飛機(jī)的設(shè)計打下了堅實的基礎(chǔ)［1］。Garrick［2］綜述了熱氣動彈性力學(xué)的發(fā)展，在傳統(tǒng)氣彈力三角形的基礎(chǔ)上補(bǔ)充了由熱效應(yīng)引起的作用力，提出了一種新的描述熱氣動彈性力學(xué)各學(xué)科關(guān)系的四面體關(guān)系圖。隨著計算機(jī)運行能力的提高和計算流體動力學(xué)（computational fluid dynamics，簡稱CFD）的發(fā)展，氣動加熱的計算精度和效率得到很大提高，為結(jié)構(gòu)在熱環(huán)境中的動力特性分析提供了準(zhǔn)確的“溫度載荷”。在結(jié)構(gòu)和氣動力耦合方面，Krist等［3］基于RANS方程開發(fā)了CFL3D求解器，為基于CFD的高／超聲速飛行器的結(jié)構(gòu)氣彈穩(wěn)定性分析提供了有效工具。由于CFL3D能夠準(zhǔn)確模擬非定常氣動力，McNamara等［1，4］基于 CFL3D分析了超聲速氣流中氣動彈性穩(wěn)定性，并取得了一定的進(jìn)展。20世紀(jì)90年代初，美國開展了NASP（25倍馬赫數(shù)）、TAV，Hy－Tech，Hyper－X等飛行器研究，為熱氣彈分析提供了大量的實驗數(shù)據(jù)［5-6］。Heeg等［7］針對 NASP 驗證模型提出了熱氣動彈性分析的5個步驟，但其并沒有考慮熱應(yīng)力引起的附加剛度。針對熱效應(yīng)對結(jié)構(gòu)振動特性的影響，史曉明等［8-9］研究了變厚度彈翼的熱動力學(xué)特性，對比有限元分析和實驗結(jié)果，發(fā)現(xiàn)結(jié)構(gòu)前3階固有頻率較常溫均有所下降。在熱氣動彈性數(shù)值計算方面，吳志剛等［10-11］提出了分層求解的思想，采用分開建立方程的松耦合解決方案，即分別進(jìn)行氣動加熱計算、結(jié)構(gòu)溫度場分析、熱模態(tài)分析、氣動彈性分析，在一個時間步上順序求解，從而使問題簡化。基于此分層求解思想，文獻(xiàn)［10］分析了高／超聲速穩(wěn)態(tài)熱環(huán)境下翼面的結(jié)構(gòu)、氣動和熱的耦合動力學(xué)問題，采用單向耦合方法進(jìn)行計算，發(fā)現(xiàn)受熱結(jié)構(gòu)的動力特性和顫振特性均可能發(fā)生變化，尤其是對于根部固支翼面。

筆者以小展弦比機(jī)翼作為研究對象，根部完全固支，探索其在高／超聲速氣流的瞬態(tài)熱環(huán)境下的動力學(xué)特性和氣彈穩(wěn)定性，考慮了熱應(yīng)力引起的剛度效應(yīng)，研究熱應(yīng)力對結(jié)構(gòu)剛度的影響，利用活塞理論計算非定常氣動力，采用p－k法求解氣彈方程［12-15］，分析討論熱效應(yīng)對結(jié)構(gòu)振動特性及顫振邊界的影響，同時還研究了不同邊界約束條件下熱應(yīng)力對結(jié)構(gòu)振動特性和氣彈穩(wěn)定性的影響。

1 有限元模型

小展弦比機(jī)翼根部弦長和尖部弦長分別為1.0m和0.3m，展長為0.4m，前緣后掠角為33.7°，機(jī)翼根部完全固支。有限元模型［16］如圖1所示。

材料為超硬鋁合金7075－Al，該型號鋁合金材主要用于制造飛機(jī)結(jié)構(gòu)及要求強(qiáng)度高、抗腐蝕性強(qiáng)的高應(yīng)力構(gòu)件。假定材料密度不隨溫度變化，ρ＝2 750kg／m3，泊松比μ＝0.3。其他材料結(jié)構(gòu)參數(shù)和熱特性參數(shù)隨溫度變化如表1所示。

表1 不同溫度下7075－Al材料參數(shù)Tab.1 Material properties of 7051－Al at differnet temperature

圖1 小展弦比機(jī)翼有限元模型示意圖Fig.1 Finit element model of low aspectration wing

2 瞬態(tài)溫度場分析

2.1 溫度場分布方程

三維瞬態(tài)溫度場在直角坐標(biāo)系下滿足微分方程

當(dāng)一個方向上（若為z方向）溫度變化為零時，方程就為二維問題的熱傳導(dǎo)方程

若只有z方向有溫度變化，則方程為一維熱傳導(dǎo)方程

其中：kx，ky，kz分別為材料沿x，y，z方向的熱傳導(dǎo)系數(shù)；φ為瞬態(tài)溫度場的場變量φ（x，y，z，t）；ρ為大氣密度；c為比熱；Q為內(nèi)部熱源。

筆者假設(shè)機(jī)翼沿厚度方向的溫度分布變化為0，即z向無溫度梯度，同時機(jī)翼內(nèi)部不產(chǎn)生內(nèi)源，則式（2）可寫成

2.2 瞬態(tài)溫度場分析

利用式（4）計算翼面溫度分布。本研究中“熱源”為熱流密度，假設(shè)機(jī)翼前緣的熱流密度為1.0×105W／m2，考慮空氣對流的影響，對流傳熱系數(shù)h＝10.0W／m2°C，環(huán)境溫度為20°C，選取50，500，1 000，2 000和3 000s時段的溫度視為“溫度載荷”進(jìn)行結(jié)構(gòu)模態(tài)分析和顫振分析，分析不同時刻下的頻率和顫振邊界，并探討其變化規(guī)律。

3 熱模態(tài)分析

熱效應(yīng)對結(jié)構(gòu)振動特性的影響主要是改變結(jié)構(gòu)總體剛度。溫度對結(jié)構(gòu)剛度的影響主要表現(xiàn)在兩個方面：a.高溫下結(jié)構(gòu)材料性能發(fā)生退化，直接降低了結(jié)構(gòu)剛度；b.不均勻的溫度場使結(jié)構(gòu)產(chǎn)生熱變形，熱變形在外界約束條件下會在結(jié)構(gòu)內(nèi)部產(chǎn)生預(yù)拉或預(yù)壓應(yīng)力，或二者同時存在，使結(jié)構(gòu)發(fā)生硬化或軟化。為此，熱環(huán)境下結(jié)構(gòu)剛度矩陣包含原剛度陣受溫度變化影響后的值KT和附加的應(yīng)力剛度Kσ，即

忽略結(jié)構(gòu)阻尼的影響，考慮熱效應(yīng)的結(jié)構(gòu)振動方程為

其中：M為結(jié)構(gòu)的質(zhì)量矩陣；ω為結(jié)構(gòu)的固有頻率；φ為結(jié)構(gòu)的模態(tài)振型。

對兩種工況進(jìn)行分析：工況1為同時考慮熱應(yīng)力引起的附加剛度及材料特性退化引起結(jié)構(gòu)剛度的下降；工況2為只考慮材料特性的退化引起結(jié)構(gòu)剛度的下降。

熱模態(tài)分析分為兩步：a.加載給定加熱時刻下的“溫度載荷”，更新剛度矩陣，選中存儲每一步的剛度矩陣；b.調(diào)用修正后的剛度矩陣進(jìn)行模態(tài)計算。利用SOL 106求解器分析不同時刻溫度載荷下的應(yīng)力，選中“正則化模態(tài)分析”選項進(jìn)行模態(tài)分析，其結(jié)構(gòu)固有頻率如表2所示。對應(yīng)的前4階常溫模態(tài)如圖2所示，常溫模態(tài)頻率與文獻(xiàn)［16］中的數(shù)據(jù)相符。由表2可知，加熱到3ks時，前4階固有頻率較常溫下降了9.4%，17.3%，12.9%和4.5%，其中扭轉(zhuǎn)模態(tài)（第2階和第3階模態(tài)）頻率下降比較大。

表2 不同溫度載荷下機(jī)翼的固有頻率（工況1）Tab.2 Model frequencies at different temperature loads（case 1） Hz

圖2 常溫下機(jī)翼前4階模態(tài)振型圖Fig.2 Mode shapes of the first 4modes

表3 不同溫度載荷下機(jī)翼的固有頻率（工況2）Tab.3 Natural frequencies of the first 4modes at different temperature loads（case 2） Hz

忽略熱應(yīng)力引起的附加剛度效應(yīng)時，固有頻率在不同時刻溫度載荷下的數(shù)據(jù)如表3所示。由表3可知，該工況下翼面加熱到3ks時，前4階固有頻率較常溫時下降了4.4%，4.9%，5.5%和4.3%。

對比兩種工況下的頻率變化曲線見圖3。只考慮材料退化對結(jié)構(gòu)剛度的影響時，結(jié)構(gòu)各階固有頻率隨著氣動加熱的進(jìn)行一直降低。當(dāng)考慮熱應(yīng)力的變化產(chǎn)生附加剛度時，熱效應(yīng)在氣動加熱初期提高了結(jié)構(gòu)固有頻率，表明此時熱應(yīng)力的變化增大了結(jié)構(gòu)的總體剛度；但隨著氣動加熱的進(jìn)行，熱應(yīng)力的引入在很大程度上降低了結(jié)構(gòu)總體剛度，特別是對扭轉(zhuǎn)剛度的影響很大，如圖3中的第2、第3支模態(tài)。

圖3 熱效應(yīng)對結(jié)構(gòu)固有頻率的影響Fig.3 Influences of thermal effect to the natural frequencies

4 顫振計算

4.1 氣彈運動方程

在超聲速、高超聲速非定常氣動力的計算中，活塞理論能很好地滿足工程精度的要求［14］?；钊碚撌且环N無粘非定常氣動理論，同時也是一種簡化的氣動力理論，只適用于在超聲速氣流下對機(jī)翼進(jìn)行顫振分析。實驗證明，馬赫數(shù)在2～5之間，用這種理論對超聲速翼面計算都能得到滿足精度要求的結(jié)果［12-14］。

在等熵條件下由動量方程可得

由于活塞前進(jìn)速度｜v｜?c∞，故有｜v／c∞｜?1，上式展開后則可略去高階微分項。當(dāng)只保留一階項時，稱為一階活塞理論，即為

當(dāng)保留二階項時，可得到二階活塞理論，即

作用在機(jī)翼表面上的壓力分布可得到上、下表面的壓力差為

對于本研究中的均勻厚度機(jī)翼，此時厚度效應(yīng)為零，則上、下翼面壓力差

在熱環(huán)境下，熱氣動彈性顫振方程［10］為

Kj包括熱應(yīng)力引起的附加剛度，其結(jié)果可由式（5）計算可得。針對方程（12），采用p－k法求解，即給定一系列速度，反復(fù)迭代求解顫振數(shù)據(jù)。

4.2 飛行環(huán)境

氣流馬赫數(shù)為3.0，空氣密度ρ＝1.226kg／m3，飛行高度為海平面，空氣密度比為1.0，機(jī)翼攻角為0°。結(jié)構(gòu)單元為實體單元，氣動網(wǎng)格點的位移和力通過無限板樣條方法插值到結(jié)構(gòu)單元上，在NASTRAN中采用p－k法進(jìn)行氣彈分析［15］。

4.3 結(jié)果分析

沿用熱模態(tài)分析中的兩種工況，工況1下，結(jié)構(gòu)在超聲速氣流中的顫振邊界如表4所示。與常溫下的顫振邊界相比，氣動加熱到3ks時，顫振速度和顫振頻率較常溫時下降了33.67%和26.49%。工況2下，顫振邊界如表5所示。相比較而言，忽略熱應(yīng)力的影響時，3ks時段的顫振速度和顫振頻率較常溫時下降了5.11%和6.40%，與表3中固有頻率的下降幅度非常接近。

圖4反應(yīng)出熱應(yīng)力的引入對結(jié)構(gòu)顫振邊界的影響是很大的，盡管加熱初期熱應(yīng)力提高了固有頻率，增大了顫振速度和顫振頻率，但隨著氣動加熱的進(jìn)行，熱應(yīng)力卻在很大程度上降低了顫振邊界，相比只有材料性能退化對顫振邊界的影響下降很多。

表4 不同溫度載荷下的顫振邊界（工況1）Tab.4 Flutter boundaries of the wing at different temperature loads（case 1）

表5 不同溫度載荷下的顫振邊界（工況2）Tab.5 Flutter boundaries of the wing at different temperature loads（case 2）

為了進(jìn)一步探討熱應(yīng)力的變化對顫振邊界的影響，筆者考慮了另外一種邊界約束條件下的顫振邊界，在根部弦處（0.2，0，0）和（0.8，0，0）兩點處固支，其對應(yīng)的顫振邊界如表6所示。

表6 不同溫度載荷下的顫振邊界（兩點約束）Tab.6 Flutter boundaries of the wing under two－node fixed boundary conditions

圖4 熱效應(yīng)對顫振特性的影響Fig.4 Influence of thermal effect to the flutter boundary

分析表6中的顫振數(shù)據(jù)，在邊界約束較少的情況下，其顫振邊界在氣動加熱初期也得到了提高。隨著氣動加熱的進(jìn)行，顫振速度和顫振頻率較常溫時下降了11.72%和12.76%。對比機(jī)翼兩種邊界約束條件下分別在兩種工況下的顫振邊界曲線，如圖5所示。從圖中可以看出，邊界約束越少，熱應(yīng)力的變化對結(jié)構(gòu)顫振邊界的影響也越小。

5 結(jié) 論

1）高溫導(dǎo)致結(jié)構(gòu)材料特性的退化直接減小原始結(jié)構(gòu)剛度，降低結(jié)構(gòu)固有頻率。

2）熱應(yīng)力的引入在氣動加熱初期提高了結(jié)構(gòu)各階固有頻率，然而隨著氣動加熱的進(jìn)行，熱應(yīng)力卻在很大程度上降低了其各階固有頻率，其中扭轉(zhuǎn)模態(tài)頻率下降最為嚴(yán)重。

3）只考慮材料退化對結(jié)構(gòu)剛度的影響時，顫振邊界的下降程度與固有頻率的下降程度接近。

圖5 不同邊界兩種工況顫振邊界曲線Fig.5 Flutter boundaries under two different boundary conditions

4）不同的邊界約束條件下熱應(yīng)力的分布是不同的，邊界約束越小，熱應(yīng)力的變化也會越小，其對結(jié)構(gòu)剛度的影響程度也越??；反之，影響會變大。

5）在翼面溫度分布趨向于穩(wěn)態(tài)的過程中，顫振邊界伴隨著氣動加熱的進(jìn)行不斷發(fā)生變化，不同結(jié)構(gòu)邊界約束條件下，熱應(yīng)力在氣動加熱初期均提高了顫振速度和顫振頻率，但最終都在很大程度上降低了顫振速度和顫振頻率。因此在高／超聲速飛行器設(shè)計過程中，必須要考慮由于氣動加熱引起的熱效應(yīng)對結(jié)構(gòu)振動特性和氣動彈性穩(wěn)定性的影響，確保飛行器的飛行安全。

［1］ McNamara J J，Thuruthimattam B J，F(xiàn)riedmann P P，et al.Hypersonic aerothermoelastic studies for reusable launch vehicles［C］∥45thAIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamic ＆ Materials Conference.Palm Springs，California：AIAA，2004：1－35.

［2］ Garrick I E.A survey of aerothermoelasticity［J］.Aerospace Engineering，1963，22（1）：140－147.

［3］ Krist S L，Biedron R T，Rumsey C L.CFL3Duser′s manual（version 5.0）［R］.［S.l.］：NASA，1998.

［4］ Friedmann P P，McNamara J J，Thuruthimattam B J，et al.Aeroelastic analysis of hypersonic vehicles［J］.Journal of Fluids and Structures，2004，19：681－712.

［5］ Jack J M，Peretz P F，Kenneth G P，et al.Three－dimensional aeroelasticity and aerothermoelastic behavior in hypersonic flow［C］∥Proceedings of 46thAIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics ＆ Materials Conference. Austin，Texas：AIAA，2005：1－47.

［6］ McNamara J J.Aeroelastic and aerothermoelastic behavior of two and three dimensional lifting surface in hypersonic flow［D］.PH.D Thesis：The University of Michigan，2005.

［7］ Heeg J，Zeiler T A，Pototzky A S，et al.Aerothermoelastic analysis of a NASP demonstrator model［R］.［S.l.］：NASA，1993.

［8］ 史曉明，楊炳淵.瞬態(tài)加熱環(huán)境下變厚度板溫度場及熱模態(tài)分析［J］.計算機(jī)輔助工程，2006，15（S）：15－18.

Shi Xiaoming，Yang Bingyuan.Temperature field and mode analysis of flat plate with thermal environment of transient heating［J］.Computer Aided Engineering，2006，15（S）：15－18.（in Chinese）

［9］ 王宏宏，陳懷海，崔旭利，等.熱效應(yīng)對導(dǎo)彈翼面固有振動特性的影響［J］.振動、測試與診斷，2010，30（3）：275－279.

Wang Honghong，Chen Huaihai，Cui Xuli，et al，Thermal effect on dynamic characteristics of a missle wing［J］.Journal of Vibration，Measurement ＆ Vibration，2010，30（3）：275－279.（in Chinese）

［10］吳志剛，惠俊鵬，楊超.高超聲速下翼面的熱顫振工程分析［J］.北京航空航天大學(xué)學(xué)報，2005，31（3）：270－273.

Wu Zhigang，Hui Junpeng，Yang Chao.Hypersonic aerothermoelastic analysis of wings［J］.Journal of Beijing University of Aeronautics and Astronautics，2005，31（3）：270－273.（in Chinese）

［11］張偉偉.超聲速、高超聲速非線性氣動彈性問題研究［D］.西安：西北工業(yè)大學(xué)，2004.

［12］陳桂彬，鄒從青，楊超.氣動彈性設(shè)計基礎(chǔ)［M］.北京：北京航空航天大學(xué)出版社，2004：107－109.

［13］葉正寅，張偉偉，史愛明，等.流固耦合力學(xué)基礎(chǔ)及其應(yīng)用［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2010：107－108.

［14］趙永輝.氣動彈性力學(xué)與控制［M］.北京：科學(xué)出版社，2007：280－282.

［15］Rodden W P，Johnson E H.MSC／NASTRAN version 68，aeroelastic analysis user′s guide［R］.Santa Ana：MSC.Software Corporation，2002.

［16］葉獻(xiàn)輝，楊翊仁.氣動加熱下三角機(jī)翼顫振［J］.西南交通大學(xué)學(xué)報，2008：43（1）：62－66.

Ye Xianhui，Yang Yiren.Flutter of delta wing under aerodynamic heating［J］.Journal of Southwest Jiaotong University，2008，43（1）：62－66.（in Chinese）