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

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

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

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

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

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

1.1 計算結(jié)構(gòu)動力學方程

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

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

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

1.2 非定常氣動載荷求解

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

1.3 CFD/CSD耦合算法

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

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

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

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

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

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

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

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

1.5 動態(tài)載荷的計算

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

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

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

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

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

2 算例分析

2.1 算例1:AGARD 445.6驗證算例

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

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

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

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

圖3 機翼有限元模型Fig.3 FEM model

圖4 表面氣動網(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)插值得到機翼的變形Fig.5 Deformation by interpolation from structure modes

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

采用CFD/CSD耦合方法計算了機翼在飛行動壓下的氣彈響應(yīng)情況。由于計算的顫振臨界動壓都遠遠大于飛行動壓(以海平面的參數(shù)計算，因此時的動壓是最嚴荷的)，故而在計算的數(shù)個狀態(tài)下，響應(yīng)都是收斂的，然而在馬赫數(shù)0.8，迎角13°的情況下，出現(xiàn)了第一、二、四階模態(tài)(彎曲)收斂，第三階模態(tài)(扭轉(zhuǎn))不收斂的現(xiàn)象，如圖7(圖中也給出了最后幾個周期內(nèi)二、三、四階廣義位移的局部放大圖)，這與一般意義上的彎扭耦合顫振發(fā)散呈現(xià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)為主的階段。為詳細確定其發(fā)散的原因，給出在一個第一階扭轉(zhuǎn)模態(tài)周期(0.003 31 s)內(nèi)四個時刻的表面壓力系數(shù)Cp的分布，如圖8。發(fā)現(xiàn)在此周期內(nèi)的大部分時刻，其表面Cp云圖與第三階模態(tài)(圖5(c))基本重合。另一方面，此時飛行動壓為45 352 Pa，遠遠小于圖6所示的馬赫數(shù)0.8時的顫振動壓(308 000 Pa)，因此，此種不收斂的響應(yīng)不應(yīng)該歸于顫振發(fā)散。故而可以確定這是在非定常氣動力作用下的強迫振動。對于這種情況，由于其在一個較長時間歷程內(nèi)都維持著較小的振幅，因而其對結(jié)構(gòu)的影響應(yīng)當進一步通過提取其響應(yīng)過程中的動態(tài)載荷來分析。

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

圖9 載荷對機翼翼根的彎矩Fig.9 Blending moment toward wing root

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

3 結(jié)論

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

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

［1] 管 德.非定?？諝鈩恿τ嬎悖跰] .北京:北京航空航天大學出版社，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ù)在氣動彈性中的應(yīng)用研究［D] .西安:西北工業(yè)大學，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] 管 德.飛機氣動彈性力學手冊［M] .北京:航空工業(yè)出版社，1994.

［9] 蔡天星，徐 敏，姚偉剛，等.基于CFD/CSD耦合的超聲速舵面動載荷分析［J] .工程力學，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.