史曉鳴，楊炳淵，李海東，唐國安

（1.復(fù)旦大學(xué) 力學(xué)與工程科學(xué)系，上海 200433；2.上海航天技術(shù)研究院 第八設(shè)計部，上海 200233）

0 引 言

隨著大機(jī)動、大過載的要求以及新的發(fā)射方式的應(yīng)用，飛行攻角顯著增加。通常在大攻角飛行時顫振邊界降低，顫振危險性比零攻角時更為突出，掌握飛行器大攻角顫振分析方法對于工程設(shè)計具有重要的現(xiàn)實意義。大攻角顫振分析的關(guān)鍵在于非定常氣動力的計算，大攻角顫振氣動力是一個高度復(fù)雜的非線性非定常流動問題，勢流理論將不再適用［1－2］。非定常CFD技術(shù)能真實地模擬大攻角非定常流場的特征，能考慮組合體部件間的干擾和大攻角帶來的非線性，但其計算規(guī)模大、耗時多、效率低下，在工程設(shè)計中全面推廣仍存在很多困難［3］。目前工程設(shè)計中常用的非定常氣動力計算方法有核函數(shù)法、面元法、超聲速格網(wǎng)法、細(xì)長體理論、活塞理論和干擾因子法等［4－5］，這些氣動力模型都是以勢流理論為基礎(chǔ)的線化模型，不能考慮大攻角情況，部件間干擾的處理也比較粗糙，不能完全滿足工程設(shè)計的需求。

飛行器工程設(shè)計迫切需要一種既能適合于較大攻角復(fù)雜外形飛行器顫振非定常氣動力計算，又能快速得到計算結(jié)果的方法。而上世紀(jì)90年代至今發(fā)展起來的當(dāng)?shù)亓骰钊碚摚?－7］及其與CFD技術(shù)的結(jié)合［8－9］，較好地解決了大攻角翼面的非定常氣動力的計算問題。本文基于當(dāng)?shù)亓骰钊碚撏茖?dǎo)了旋成體機(jī)身和任意外形三維機(jī)翼非定常氣動力計算公式，利用CFD數(shù)值求解定常流場當(dāng)?shù)亓鲄?shù)，使當(dāng)?shù)亓骰钊碚撏茝V用于大攻角三維機(jī)翼和旋成體機(jī)身，并可考慮部件間的干擾，從而可進(jìn)行全機(jī)組合體大攻角顫振計算。以翼－身組合體外形飛行器為例，完成了Ma＝2～7，攻角α＝0°～20°的顫振計算，計算結(jié)果經(jīng)與非定常氣動力采用CFD數(shù)值方法計算的CFD／CSD直接耦合的顫振計算結(jié)果對比，精度滿足工程設(shè)計需要，其計算速度卻提高了百倍以上。

1 翼－身組合體飛行器的動力學(xué)模型

考慮由機(jī)身、機(jī)翼、垂尾組成的面對稱布局翼－身組合體。忽略結(jié)構(gòu)阻尼，由分枝模態(tài)法可導(dǎo)出結(jié)構(gòu)的運(yùn)動微分方程［10－11］：

式中，M和K分別為結(jié)構(gòu)的模態(tài)質(zhì)量和模態(tài)剛度矩陣，QA為廣義氣動力矢量，由機(jī)身、機(jī)翼、垂尾各部件的廣義氣動力組成：

式中Δp為飛行器表面非定常壓力分布，Φ為各個分枝模態(tài)振型函數(shù)組成的列矢量。

2 非定常氣動力

本節(jié)討論式（2）中機(jī)身、機(jī)翼、垂尾各個部件的非定常氣動力計算方法。考慮圖1所示飛行器表面任意曲面，過曲面上任意點o作切面P，假設(shè)當(dāng)?shù)亓魉俣萔L沿切面的任一方向，以o為原點作局部坐標(biāo)系oξηζ，其中ξ軸方向與VL方向重合，ζ軸與過o點的外法線方向重合。由當(dāng)?shù)亓骰钊碚摚嫔先我恻c的非定常氣動力為：

式中ρL、cL為曲面任意點的當(dāng)?shù)貧饬髅芏群吐曀伲琖m為物面擾動引起的下洗。假設(shè)垂直于物面的振動表示為ζ＝ζ（ξ，η），則物面振動速度和轉(zhuǎn)角引起的下洗為

圖1 任意曲面上的局部坐標(biāo)系Fig.1 Local coordinate system

對于常見的軸對稱旋成體機(jī)身，采用圖2所示的柱坐標(biāo)系，表面方程可用母線方程r＝R（x）表示，機(jī)身垂直于機(jī)身軸線的振動為Z＝Z（x，t），其外法線對z軸的方向余弦為

式中μb＝［1＋（?R／?x）］－1／2為機(jī)身母線切線與x軸夾角的余弦，因此垂直于機(jī)身表面的振動為

圖2 機(jī)身坐標(biāo)系Fig.2 Coordinate system of fuselage

對于如圖3所示翼面，則采用機(jī)翼局部坐標(biāo)系，設(shè)對稱翼型的翼面以中性面為參考平面的型面方程為H＝H，），以離開中性面為正，翼面垂直于中性面振動為＝（，t）。則垂直于翼型面上、下表面的振動：

式中上翼面取正號，下翼面取負(fù)號，μw＝［1＋（?H／?）2＋（?H／?）2］－1／2為該點翼型面外法線對軸的方向余弦。

圖3 翼面坐標(biāo)系Fig.3 Coordinate system of wing

式（6）和（7）代入式（4），并利用方向?qū)?shù)公式，可得垂直于機(jī)身表面的下洗

以及垂直于機(jī)翼上、下表面的下洗

式（8）、（9）中，式中VxL為機(jī)身表面當(dāng)?shù)貧饬魉俣仍趚向的分量；VL、VL分別為機(jī)翼表面當(dāng)?shù)貧饬魉俣鹊?、向分量。同理可得垂尾與式（9）相同的下洗表達(dá)式，只是式中坐標(biāo)替換為x替換為z，位移替換為Y。將機(jī)身、機(jī)翼和垂尾的下洗代入式（3）得到非定常壓力分布，再代入式（2），并做模態(tài)坐標(biāo)變換，可得廣義氣動力Q、Q和Q。

通過CFD數(shù)值仿真求解不同馬赫數(shù)不同攻角下飛行器定常流的方法即可獲得式（3）、（8）、（9）中的當(dāng)?shù)亓鲄?shù)。因此，各種復(fù)雜外形飛行器部件之間的互相干擾、流場三維效應(yīng)以及大攻角的非線性等都可在當(dāng)?shù)亓鲄?shù)中得到反映，從而使得當(dāng)?shù)亓骰钊碚撃軌蜻m應(yīng)復(fù)雜外形飛行器和部件的大攻角非定常氣動力計算和顫振計算，計算精度得以提高而且時間成本大大下降。

3 算例與分析

以文獻(xiàn)［10］已完成模態(tài)計算的面對稱布局翼－身組合體飛行器為本文算例，各分枝模態(tài)如下選?。簩ΨQ分枝取剛體沉浮、俯仰和前三階彈性模態(tài)；反對稱分枝取剛體滾轉(zhuǎn)和前四階彈性模態(tài)；垂尾分枝取前五階彈性模態(tài)。表1給出各分枝模態(tài)頻率，振型可參見文獻(xiàn)［10］。采用頻域復(fù)特征值方法搜索顫振臨界參數(shù)。

用本文結(jié)合CFD的當(dāng)?shù)亓骰钊碚摲椒ㄍ瓿闪撕Ｆ矫娓叨认?，Ma＝2～7，攻角α＝0°～20°的全機(jī)組合體顫振計算。圖4為Ma＝4來流下不同攻角的顫振邊界，圖中同時給出非定常氣動力采用全數(shù)值仿真計算的CFD／CSD直接耦合方法的計算結(jié)果，本文結(jié)果與其相比，在臨界動壓數(shù)值上最大誤差為28.46%，換算為臨界速度的最大誤差僅13.14%，計算精度滿足工程設(shè)計要求。本文方法一個馬赫數(shù)點的計算時間在1h左右，而CFD／CSD全數(shù)值仿真方法的計算時間超過100h，在滿足精度要求的情況下，計算效率提高100倍以上，更適合于在工程設(shè)計中應(yīng)用。圖5為不同來流馬赫數(shù)和攻角下的顫振邊界，從中可見攻角對顫振邊界的影響不可忽略，尤其是來流為高馬赫數(shù)時，攻角增大使得顫振邊界大幅下降。

表1 各個分枝結(jié)構(gòu)固有頻率（Hz）Table 1 Frequencies of each branch（Hz）

圖4 Ma＝4來流狀態(tài)下不同攻角顫振邊界Fig.4 Flutter critical dynamic pressure at different angles of attack when Ma＝4

圖5 不同馬赫數(shù)下顫振邊界隨攻角變化Fig.5 Flutter critical dynamic pressure at different angles of attack and Mach numbers

4 結(jié) 論

以當(dāng)?shù)亓骰钊碚摓榛A(chǔ)，假設(shè)當(dāng)?shù)亓餮仫w行器表面任意切線方向，推導(dǎo)了旋成體機(jī)身和任意外形三維機(jī)翼非定常壓力分布計算公式，結(jié)合CFD數(shù)值求解定常當(dāng)?shù)亓鲄?shù)，形成了可用于全機(jī)組合體大攻角超聲速顫振計算的方法。算例結(jié)果表明本文方法計算精度接近全CFD／CSD數(shù)值仿真方法，計算效率則大幅度提高，適合工程設(shè)計計算要求。

算例結(jié)果顯示攻角對顫振邊界有很大影響，飛行器設(shè)計中需充分予以考慮。尤其是高馬赫數(shù)時，攻角增大使得顫振邊界大幅下降。

［1］楊炳淵，史曉鳴，梁強(qiáng).高速有翼導(dǎo)彈多場耦合動力學(xué)的研究和進(jìn)展（上）［J］.強(qiáng)度與環(huán)境，2008，35（5）：55－63.（YANG Bing－yuan，SHI Xiao－ming，LIANG Qiang.Investigation and development of the multi－physics coupling dynamics on the hypersonic winged missiles （I）［J］.Structure＆EnvironmentEngineering，2008，35（5）：55－63.）

［2］葉正寅，謝飛.中等及大迎角下的翼型氣動彈性性質(zhì)研究［J］.風(fēng)機(jī)技術(shù)，2009，（2）：9－13.（YE Zheng－yin，XIE Fei.Research on the aeroelastic characteristics of an airfoil in moderate and high incidences［J］.Compressor，Blower＆FanTechnology，2009，（2）：9－13.）

［3］SCHUSTER D M，LIU D D.Computational aeroelasticity：success，progress，challenge［J］.Journal of Aircraft，2003，40（5）：843－856.

［4］McNAMARA J J，F(xiàn)RIEDMANN P P.Aeroelastic and aerothermoelastic analysis of hypersonic vehicles：current status and future trends［R］.AIAA 2007－2013.

［5］CHAE S，HODGES D H.Dynamics and aeroelastic analysis of missile［R］.AIAA 2003－1968.

［6］楊炳淵，宋偉力.應(yīng)用當(dāng)?shù)亓骰钊碚摰拇蠊ソ巧γ骖澱駳鈩恿Ρ磉_(dá)式［J］.力學(xué)季刊，1999，20（3）：223－228.（YANG Bing－yuan，SONG Wei－li.Expressions about aerodynamic forces of flutter for wing with high angle of attack by local flow piston theory［J］.ChineseQuarterly ofMechanics，1999，20（3）：223－228.）

［7］楊炳淵，宋偉力.前緣激波脫體的大迎角翼面顫振工程計算方法［J］.振動與沖擊，2002，22（4）：100－103.（YANG Bing－yuan，SONG Wei－li.Engineering approximation of high attack angle flutter for wing with bow wave on leading edge［J］.JournalofVibrationand Shock，2002，22（4）：100－103.）

［8］張偉偉，葉正寅.基于當(dāng)?shù)亓骰钊碚摰臍鈩訌椥杂嬎惴椒ㄑ芯浚跩］.力學(xué)學(xué)報，2005，37（5）：632－639.（ZHANG Wei－wei，YE Zheng－yin.Numerical method of aeroelasticity based on local piston theory［J］.Acta MechanicaSinica，2005，37（5）：632－639.）

［9］ZHANG Wei－wei，YE Zheng－yin，ZHANG Chen－an.Supersonic flutter analysis based on a local piston theory［J］.AIAAJournal，2009，47（10）：2321－2328.

［10］史曉鳴，許泉，唐國安.基于分枝模態(tài)法的面對稱布局飛行器結(jié)構(gòu)動力學(xué)建模［J］.上海航天，2011，28（2）：27－31.（SHI Xiao－ming，XU Quan，TANG Guo－an.The dynamics modeling of plane symmetrical vehicle structural based on branch mode method［J］.AerospaceShanghai，2011，28（2）：27－31.）

［11］楊炳淵，史曉鳴，梁強(qiáng).高速有翼導(dǎo)彈多場耦合動力學(xué)的研究和進(jìn)展（下）［J］.強(qiáng)度與環(huán)境，2008，35（6）：56－62.（YANG Bing－yuan，SHI Xiao－ming，LIANG Qiang.Investigation and development of the multi－physics coupling dynamics on the hypersonic winged missiles （II）［J］.Structure＆EnvironmentEngineering，2008，35（6）：56－62.）