基于CFD/CSD耦合的顫振與動(dòng)載荷分析方法

謝 亮，徐 敏，李 杰，蔡天星

(1.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072;2.西北工業(yè)大學(xué) 航天學(xué)院，西安 710072)

顫振是一種極具破壞性的氣動(dòng)彈性現(xiàn)象，可能導(dǎo)致飛行器在短時(shí)間內(nèi)解體性破壞。因此顫振的計(jì)算是氣動(dòng)彈性領(lǐng)域內(nèi)十分重要的一項(xiàng)內(nèi)容。傳統(tǒng)的顫振計(jì)算中采用求解線化速度勢(shì)方程［1]計(jì)算非定常氣動(dòng)力，僅能考慮氣動(dòng)力的線性效應(yīng)，當(dāng)飛行器處于跨聲速、高超聲速階段或者結(jié)構(gòu)出現(xiàn)大變形時(shí)氣動(dòng)力的非線性效應(yīng)比較明顯，則其計(jì)算結(jié)果可信度不足。目前計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)技術(shù)已發(fā)展得較為成熟，它與計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(Computational Structural Dynamics，CSD)的結(jié)合為氣彈分析提供了強(qiáng)有力的工具，故而CFD/CSD耦合方法已經(jīng)成為當(dāng)前氣彈仿真的主流手段［2－4]。然而采用CFD/CSD耦合算法的計(jì)算效率較低，為此建立起CFD/CSD全耦合系統(tǒng)的降階模型(Reduced order model，ROM)有其必要性［5]。氣動(dòng)彈性學(xué)家 Dowell［6]將目前應(yīng)用于氣動(dòng)彈性領(lǐng)域的降階模型分成兩類:一類是基于模態(tài)的方法，另一類是基于信號(hào)的方法。在后一類ROM中，基于Volterra級(jí)數(shù)的降階模型應(yīng)用較為廣泛。在上世紀(jì)90年代，NASA的 Silva［7]首先利用脈沖響應(yīng)辨識(shí) Volterra核，并使用特征實(shí)現(xiàn)算法構(gòu)建非定常氣動(dòng)力模型，取得了不錯(cuò)的效果。其后，Raveh針對(duì)脈沖響應(yīng)存在數(shù)值不穩(wěn)定的問題，改用階躍響應(yīng)辨識(shí)Valterra核，文獻(xiàn)［5] 指出采用階躍辨識(shí)得到的一階核包含了二階核中的非對(duì)角分量，因而能體現(xiàn)更多的非線性。

除計(jì)算顫振臨界速度外，氣動(dòng)彈性分析還要承擔(dān)飛行器在動(dòng)態(tài)氣動(dòng)載荷與自身慣性力耦合作用下的動(dòng)態(tài)載荷計(jì)算的任務(wù)。因?yàn)閯?dòng)載荷的計(jì)算是結(jié)構(gòu)動(dòng)態(tài)設(shè)計(jì)的關(guān)鍵技術(shù)之一，它為結(jié)構(gòu)響應(yīng)計(jì)算、結(jié)構(gòu)的動(dòng)態(tài)設(shè)計(jì)與故障分析提供了可靠的依據(jù)。目前工程上并不作動(dòng)載荷的計(jì)算［8]，僅求解飛行器的靜氣彈響應(yīng)方程，計(jì)算出靜變形情況下的氣動(dòng)載荷，并乘上一個(gè)安全系數(shù)作為飛行器的動(dòng)態(tài)載荷。目前使用CFD/CSD耦合方法進(jìn)行顫振分析已經(jīng)比較成熟，但是用此手段進(jìn)行飛行器動(dòng)載荷的分析還比較少見。僅見有文獻(xiàn)［9] 進(jìn)行了超聲速情況下舵面的動(dòng)載荷分析。

本文采用CFD/CSD耦合算法建立了氣彈仿真系統(tǒng)，詳細(xì)論述了仿真過程中提取動(dòng)態(tài)載荷(包括氣動(dòng)載荷與慣性載荷)的方法?；谙到y(tǒng)辨識(shí)的方法，使用Volterra級(jí)數(shù)建立了降階模型，提供了快速計(jì)算顫振邊界的手段。分別采用CFD/CSD全耦合方法和ROM計(jì)算了AGARD 445.6機(jī)翼的顫振邊界，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比，驗(yàn)證了計(jì)算程序的可靠性。在此基礎(chǔ)上采用ROM分析了一型帶邊條平直翼的顫振邊界，使用CFD/CSD耦合方法計(jì)算了此機(jī)翼在飛行動(dòng)壓下的氣彈響應(yīng)，結(jié)果表明即使是在顫振邊界內(nèi)仍有可能出現(xiàn)極限環(huán)振蕩(Limit cycle oscillations－LCO)。對(duì)于0.8馬赫，迎角13°情況下出現(xiàn)的LCO情況，分析了其響應(yīng)過程中的動(dòng)載荷情況。

1 基于CFD/CSD耦合的氣彈仿真

1.1 計(jì)算結(jié)構(gòu)動(dòng)力學(xué)方程

將結(jié)構(gòu)模型離散化之后通過非保守系統(tǒng)的拉格朗日方程可以得到如下的結(jié)構(gòu)動(dòng)力學(xué)方程:

其中M，C，K分別為結(jié)構(gòu)質(zhì)量矩陣、阻尼矩陣和剛度矩陣，QF為結(jié)構(gòu)外激勵(lì)。本文中采用振型疊加法求解上述方程。通過求解特征方程Kφi=λiMφi，獲得系統(tǒng)的特征值 λ1，…，λN和特征向量 φ1，…，φN，忽略對(duì)系統(tǒng)的實(shí)際響應(yīng)貢獻(xiàn)較小的高階振型，保留低階n個(gè)振型，組成振型矩陣 Φ =［φ1，φ2，…，φn] ，則式(1)可近似為如下的n階方程:

式中，M=ΦTMΦ，C=ΦTCΦ，K=ΦTKΦ分別為廣義質(zhì)量、廣義阻尼和廣義剛度矩陣。ξ為廣義位移，滿足:q為廣義激勵(lì)力。

1.2 非定常氣動(dòng)載荷求解

對(duì)于方程(1)的右端項(xiàng)QF，即非定常氣動(dòng)載荷，采用CFD方法求解。積分形式的Euler方程為:

1.3 CFD/CSD耦合算法

圖1 CFD/CSD耦合方式Fig.1 CFD/CSD coupling method

采用松耦合算法在時(shí)域內(nèi)分別求解流體、結(jié)構(gòu)動(dòng)力學(xué)方程，交錯(cuò)時(shí)間推進(jìn)獲得耦合系統(tǒng)的響應(yīng)，其流程圖如圖1。這種算法的最大好處是能充分利用現(xiàn)有的計(jì)算流體動(dòng)力學(xué)和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)的程序，只要增加流體－結(jié)構(gòu)數(shù)據(jù)交換模塊即可，從而保持程序的模塊化。雖然結(jié)構(gòu)求解器和流體求解器都具有二階及二階以上的精度，但是此種耦合方法由于流體和結(jié)構(gòu)積分的不同步勢(shì)必造成數(shù)值誤差，當(dāng)積分時(shí)間步取得較大時(shí)計(jì)算結(jié)果有可能失真。為解決這個(gè)問題，對(duì)松耦合方法作了如下改進(jìn)［5]:

在氣動(dòng)彈性計(jì)算中非定常氣動(dòng)力和結(jié)構(gòu)位移一般都呈連續(xù)變化。為此，對(duì)圖1中的步驟①做如下修正:通過時(shí)刻的氣動(dòng)力(均已知)做三點(diǎn)拋物插值，擬合得出tn～tn+1時(shí)間內(nèi)氣動(dòng)力關(guān)于時(shí)間的二次函數(shù)。將這個(gè)隨時(shí)間變化的氣動(dòng)力代入結(jié)構(gòu)動(dòng)力學(xué)方程(1)，推進(jìn)求得tn+1時(shí)刻的結(jié)構(gòu)參數(shù)。這個(gè)修正在很大程度上提高了CFD/CSD松耦合計(jì)算的精度。

結(jié)構(gòu)流體數(shù)據(jù)交換采用無限平板樣條(Infinite plate spline，IPS)方法［10]，動(dòng)網(wǎng)格采用無限代數(shù)插值(Transfinite interpolation infinite，TFI)實(shí)現(xiàn)。

1.4 基于Volterra級(jí)數(shù)的非定常氣動(dòng)降階模型。

Volterra級(jí)數(shù)是一個(gè)無窮級(jí)數(shù)，小擾動(dòng)下，Euler方程計(jì)算得到的非定常氣動(dòng)力可以表示為二階Volterra級(jí)數(shù)形式，時(shí)域內(nèi)離散形式為:

采用階躍響應(yīng)辨識(shí)非定常氣動(dòng)力，給CFD/CSD耦合系統(tǒng)如下階躍響應(yīng):可辨識(shí)得到包含二階核對(duì)角分量及部分非對(duì)角分量的近似一階核為:

利用得到的volterra核通過特征實(shí)現(xiàn)算法［5]可以得到以狀態(tài)空間形式表示的非定常氣動(dòng)力的降階模型，此狀態(tài)空間模型以廣義位移為輸入，以廣義力為輸出，通過此模型計(jì)算方程(1)的右端項(xiàng)(即氣動(dòng)力載荷)，結(jié)構(gòu)動(dòng)力學(xué)方程仍然采用模態(tài)疊加法，則可以得到整個(gè)氣動(dòng)彈性系統(tǒng)的低階狀態(tài)空間模型，通過調(diào)整動(dòng)壓分析系統(tǒng)的穩(wěn)定性來求顫振邊界。

1.5 動(dòng)態(tài)載荷的計(jì)算

飛行器動(dòng)態(tài)載荷計(jì)算是結(jié)構(gòu)設(shè)計(jì)和故障診斷的重要環(huán)節(jié)，僅通過靜氣彈結(jié)構(gòu)設(shè)計(jì)有時(shí)并不能滿足動(dòng)強(qiáng)度要求。傳統(tǒng)的氣動(dòng)彈性分析方法通過頻域或時(shí)域識(shí)別方法來確定氣動(dòng)載荷［1]，然而當(dāng)結(jié)構(gòu)變形較大或者飛行器處于跨聲速、高超聲速階段時(shí)，非定常氣動(dòng)力存在明顯的非線性效應(yīng)，只能通過CFD技術(shù)實(shí)現(xiàn)氣動(dòng)載荷提取。本文采用CFD/CSD耦合計(jì)算飛行器的動(dòng)態(tài)響應(yīng)，提取了非定常氣動(dòng)力載荷和慣性載荷，實(shí)現(xiàn)了動(dòng)載荷分析。

對(duì)于式(1)，右端項(xiàng)表示飛行器所受到的氣動(dòng)載荷，除此外，飛行器還受到自身慣性力施加的慣性載荷，即方程(1)左端第一項(xiàng)M·q·。飛行器所受到的真實(shí)載荷即由這兩項(xiàng)組成:

在CFD/CSD耦合計(jì)算的每一個(gè)時(shí)間步中，分別提取出氣動(dòng)載荷QF與慣性載荷，按式(8)求和，得總載荷。

氣動(dòng)載荷的求解依靠CFD系統(tǒng)，在雙時(shí)間推進(jìn)求解流場(chǎng)的每一真實(shí)時(shí)間推進(jìn)過程中，當(dāng)內(nèi)層迭代(即偽時(shí)間步推進(jìn))收斂后，利用物面的壓強(qiáng)沿物面作數(shù)值積分即可求得該時(shí)刻的瞬態(tài)氣動(dòng)載荷。慣性載荷的求解有兩種方案，一種是將求解得到的隨時(shí)間變化的瞬態(tài)氣動(dòng)載荷加載到CSD系統(tǒng)上進(jìn)行動(dòng)態(tài)響應(yīng)分析，從而求解得到整個(gè)響應(yīng)過程中慣性載荷及總載荷，同時(shí)可以獲得此過程中的瞬態(tài)應(yīng)力，文獻(xiàn)［9] 即采用的此種方案。本文采用的是另一種方法:將使用有限元方法求得的CSD系統(tǒng)的總體質(zhì)量陣M導(dǎo)入CFD/CSD耦合系統(tǒng)，利用計(jì)算動(dòng)態(tài)慣性載荷。由于CFD/CSD耦合中CSD系統(tǒng)采用的是模態(tài)迭加法，故而無法直接按求解慣性載荷，但是注意到q=Φξ，則慣性載荷的求解可使用下式:

上式中，下標(biāo)i表示第i階模態(tài)，Φi與ξi分別是該階模態(tài)對(duì)應(yīng)的振型與廣義加速度，n表示截取的模態(tài)階數(shù)。M，Φi由有限元方法進(jìn)行模態(tài)分析得到，而ξi由CFD/CSD耦合方法求解，利用求解得到的廣義速度ξi作數(shù)值差分求解得到。

2 算例分析

2.1 算例1:AGARD 445.6驗(yàn)證算例

分別采用CFD/CSD耦合方法與ROM求解了AGARD 445.6機(jī)翼的顫振邊界。首先采用有限元法對(duì)AGARD 445.6進(jìn)行模態(tài)分析，選擇前四階模態(tài)組成解耦后的結(jié)構(gòu)動(dòng)力學(xué)方程。按前述方法將之與CFD系統(tǒng)耦合建立CFD/CSD時(shí)域仿真系統(tǒng)。建立ROM時(shí)，計(jì)算模態(tài)仍選取前四階模態(tài)，對(duì)CFD/CSD耦合系統(tǒng)施加階躍響應(yīng)，響應(yīng)幅值取1.0E－4，辨識(shí)物理時(shí)間步長(zhǎng)取1.0E－4s，計(jì)算1000個(gè)時(shí)間步長(zhǎng)內(nèi)的非定常氣動(dòng)力，按前述方法建立起系統(tǒng)的低階狀態(tài)空間模型。CFD/CSD耦合算法及ROM計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比見圖2，可見計(jì)算值與實(shí)驗(yàn)值［11]相符較好，且與國內(nèi)外計(jì)算結(jié)果一致［2－5]，由此驗(yàn)證了CFD/CSD耦合算法與程序的準(zhǔn)確性與可靠性。CFD/CSD耦合方法與ROM計(jì)算結(jié)果幾乎重合，表明在顫振計(jì)算中可以使用ROM代替CFD/CSD全耦合系統(tǒng)，實(shí)現(xiàn)顫振邊界的快速求解，對(duì)本算例，CFD/CSD耦合仿真計(jì)算時(shí)間約為降階模型計(jì)算的10～100倍。

圖2 AGARD445.6機(jī)翼顫振邊界Fig.2 Flutter boundary of AGARD 445.6

2.2 算例2:帶邊條的平直翼的顫振及動(dòng)載荷分析

機(jī)翼幾何形狀及有限元模型見圖3，表面氣動(dòng)網(wǎng)格見圖4.采用有限元法其進(jìn)行模態(tài)分析，截取前四階模態(tài)組成解耦后的CSD系統(tǒng)。前四階模態(tài)頻率及類型見表1。各階模態(tài)的變形通過IPS方法插值到氣動(dòng)網(wǎng)格上的變形見圖5，可見變形網(wǎng)格光滑，插值效果優(yōu)良。

圖3 機(jī)翼有限元模型Fig.3 FEM model

圖4 表面氣動(dòng)網(wǎng)格Fig.4 CFD grid of case 2

表1 算例1結(jié)構(gòu)前四階模態(tài)Tab.1 Lower 4 structure modes of case 1

圖5 模態(tài)插值得到機(jī)翼的變形Fig.5 Deformation by interpolation from structure modes

鑒于ROM的準(zhǔn)確性已由AGARD 445.6標(biāo)模算例得到了驗(yàn)證，為求快速獲得結(jié)果，故而此處直接使用ROM進(jìn)行顫振分析，得到此機(jī)翼顫振邊界如圖6所示。

采用CFD/CSD耦合方法計(jì)算了機(jī)翼在飛行動(dòng)壓下的氣彈響應(yīng)情況。由于計(jì)算的顫振臨界動(dòng)壓都遠(yuǎn)遠(yuǎn)大于飛行動(dòng)壓(以海平面的參數(shù)計(jì)算，因此時(shí)的動(dòng)壓是最嚴(yán)荷的)，故而在計(jì)算的數(shù)個(gè)狀態(tài)下，響應(yīng)都是收斂的，然而在馬赫數(shù)0.8，迎角13°的情況下，出現(xiàn)了第一、二、四階模態(tài)(彎曲)收斂，第三階模態(tài)(扭轉(zhuǎn))不收斂的現(xiàn)象，如圖7(圖中也給出了最后幾個(gè)周期內(nèi)二、三、四階廣義位移的局部放大圖)，這與一般意義上的彎扭耦合顫振發(fā)散呈現(xiàn)出不同的特點(diǎn)。

分析位移響應(yīng)曲線，從幅值上判斷，氣彈響應(yīng)過程中開始階段以彎曲響應(yīng)為主，但是彎曲響應(yīng)(一、二、四階模態(tài))是收斂的，而第一階扭轉(zhuǎn)模態(tài)發(fā)散，最后呈現(xiàn)出極限環(huán)振蕩的形式，使響應(yīng)過程變?yōu)橐耘まD(zhuǎn)響應(yīng)為主的階段。為詳細(xì)確定其發(fā)散的原因，給出在一個(gè)第一階扭轉(zhuǎn)模態(tài)周期(0.003 31 s)內(nèi)四個(gè)時(shí)刻的表面壓力系數(shù)Cp的分布，如圖8。發(fā)現(xiàn)在此周期內(nèi)的大部分時(shí)刻，其表面Cp云圖與第三階模態(tài)(圖5(c))基本重合。另一方面，此時(shí)飛行動(dòng)壓為45 352 Pa，遠(yuǎn)遠(yuǎn)小于圖6所示的馬赫數(shù)0.8時(shí)的顫振動(dòng)壓(308 000 Pa)，因此，此種不收斂的響應(yīng)不應(yīng)該歸于顫振發(fā)散。故而可以確定這是在非定常氣動(dòng)力作用下的強(qiáng)迫振動(dòng)。對(duì)于這種情況，由于其在一個(gè)較長(zhǎng)時(shí)間歷程內(nèi)都維持著較小的振幅，因而其對(duì)結(jié)構(gòu)的影響應(yīng)當(dāng)進(jìn)一步通過提取其響應(yīng)過程中的動(dòng)態(tài)載荷來分析。

按1.5節(jié)所述方法提取出響應(yīng)過程中的動(dòng)態(tài)載荷信息，分別計(jì)算慣性載荷、氣動(dòng)載荷，按式(8)求總載荷，并分別求其對(duì)翼根的彎矩，結(jié)果如圖9示。計(jì)算所得氣動(dòng)載荷在以彎曲響應(yīng)為主的階段峰值為275 N·m，進(jìn)入極限環(huán)后峰值為335 N·m，慣性載荷在以彎曲響應(yīng)為主的階段峰值為22 N·m，進(jìn)入極限環(huán)后峰值為11 N·m，總載荷在彎扭耦合階段的峰值為285 N·m，進(jìn)入極限環(huán)后峰值為325 N·m。由圖可見，在響應(yīng)前期的以彎曲響應(yīng)為主的階段，機(jī)翼受到的總載荷大于氣動(dòng)載荷，而當(dāng)響應(yīng)進(jìn)入扭轉(zhuǎn)響應(yīng)為主的階段時(shí)，由于單純的扭轉(zhuǎn)導(dǎo)致的附加氣動(dòng)力起阻振作用，故而總載荷小于氣動(dòng)載荷，是氣動(dòng)載荷與慣性載荷之差。在本算例中，定常氣動(dòng)載荷為230 N·m，動(dòng)態(tài)載荷的峰值是它的1.4倍左右，這與文獻(xiàn)［9] 所給出的超聲速舵面在高超聲速情況下總載荷是定常氣動(dòng)載荷的3倍左右的結(jié)論不同，則可認(rèn)為動(dòng)態(tài)載荷與靜態(tài)載荷之間的關(guān)系與飛行器的構(gòu)型、飛行條件相關(guān)，因此在飛行器設(shè)計(jì)過程中，除了要進(jìn)行靜態(tài)載荷分析外，仍有必要分析動(dòng)態(tài)載荷情況。而使用CFD/CSD耦合的方法提供了提取彈性響應(yīng)過程中動(dòng)態(tài)載荷的有效手段。

圖9 載荷對(duì)機(jī)翼翼根的彎矩Fig.9 Blending moment toward wing root

飛行器也有可能出現(xiàn)極限環(huán)振蕩，對(duì)此應(yīng)當(dāng)予以注意。提取了響應(yīng)過程中的動(dòng)態(tài)載荷信息，結(jié)果表明氣彈響應(yīng)過程中，彎曲響應(yīng)導(dǎo)致的瞬態(tài)總載荷比氣動(dòng)載荷要大，而扭轉(zhuǎn)響應(yīng)過程中的總載荷小于氣動(dòng)載荷，因此在飛行器設(shè)計(jì)過程中，除了進(jìn)行靜態(tài)載荷的計(jì)算外，仍有必要分析動(dòng)態(tài)載荷，而基于CFD/CSD耦合而建立起來的氣動(dòng)彈性時(shí)域仿真系統(tǒng)提供了提取出響應(yīng)過程中動(dòng)態(tài)載荷(包括氣動(dòng)載荷與慣性載荷)的有效手段。

3 結(jié)論

采用CFD/CSD耦合方法計(jì)算了AGARD 445.6機(jī)翼的顫振邊界，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果相符較好，由此驗(yàn)證了本文所采用的CFD/CSD耦合算法的準(zhǔn)確性與可靠性?；赩olterra級(jí)數(shù)建立了CFD/CSD耦合系統(tǒng)的降階模型，實(shí)現(xiàn)了顫振的快速求解，AGARD 445.6算例的結(jié)果表明ROM的計(jì)算結(jié)果與CFD/CSD全耦合方法的結(jié)果基本一致，而效率上提高了一到兩個(gè)數(shù)量級(jí)。

用CFD/CSD耦合方法模擬了機(jī)翼在飛行動(dòng)壓下的氣彈響應(yīng)歷程，結(jié)果表明即使飛行動(dòng)壓遠(yuǎn)小于顫振邊界，由于空氣動(dòng)力的非線性與結(jié)構(gòu)的非線性因素，飛行器也有可能出現(xiàn)極限環(huán)振蕩，對(duì)此應(yīng)當(dāng)予以注意。提取了響應(yīng)過程中的動(dòng)態(tài)載荷信息，結(jié)果表明氣彈響應(yīng)過程中，彎曲響應(yīng)導(dǎo)致的瞬態(tài)總載荷比氣動(dòng)載荷要大，而扭轉(zhuǎn)響應(yīng)過程中的總載荷小于氣動(dòng)載荷，因此在飛行器設(shè)計(jì)過程中，除了進(jìn)行靜態(tài)載荷的計(jì)算外，仍有必要分析動(dòng)態(tài)載荷，而基于CFD/CSD耦合而建立起來的氣動(dòng)彈性時(shí)域仿真系統(tǒng)提供了提取出響應(yīng)過程中動(dòng)態(tài)載荷(包括氣動(dòng)載荷與慣性載荷)的有效手段。

［1] 管 德.非定?？諝鈩?dòng)力計(jì)算［M] .北京:北京航空航天大學(xué)出版社，1991.

［2] Liu F.Calculation of wing flutter by a coupled fluid-structure method［J] .Journal of Aircraft，2001，38(2):334 －342.

［3] Elizabeth M，Lee-Rausch ，Batina H T.Calculation of AGARD wing 445.6 flutter using navier-stokes aerodynamics［R] .AIAA，93－3476，1993.

［4] Bennett R M，Edwards J W.An overview of recent developments in computational aeroelasticity［R] .AIAA，1998，98 －2421.

［5] 曾憲昂.模型降階技術(shù)在氣動(dòng)彈性中的應(yīng)用研究［D] .西安:西北工業(yè)大學(xué)，2008.

［6] Dowell E H，Hall K C.Modeling of fluid-structure interaction［J] .Annual Review Fluid Mechanics，1993，33:445 －490.

［7] Silva W A.Application of nonlinear systems theory to transonic unsteady aerodynamic responses［J] .Journal of Aircraft，1993，30(5):600 －668.

［8] 管 德.飛機(jī)氣動(dòng)彈性力學(xué)手冊(cè)［M] .北京:航空工業(yè)出版社，1994.

［9] 蔡天星，徐 敏，姚偉剛，等.基于CFD/CSD耦合的超聲速舵面動(dòng)載荷分析［J] .工程力學(xué)，2011，28(3):245 －250.

［10] Goura G S L .Time marching analysis of flutter using computational fluid dynamics ［D] . University of Glasgow，2001.

［11] Yates E C.AGARD standard aeroelastic configurations for dynamic response I-Wing 445.6 ［R] .Agard Report No.765，The 61stMeeting of the Structures and Materials Panel at Oberammergau，Germany:September，1985.