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

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

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

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

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

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

1.1 計算結構動力學方程

將結構模型離散化之后通過非保守系統的拉格朗日方程可以得到如下的結構動力學方程:

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

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

1.2 非定常氣動載荷求解

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

1.3 CFD/CSD耦合算法

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

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

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

結構流體數據交換采用無限平板樣條(Infinite plate spline，IPS)方法［10]，動網格采用無限代數插值(Transfinite interpolation infinite，TFI)實現。

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

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

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

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

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

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

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

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

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

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

2 算例分析

2.1 算例1:AGARD 445.6驗證算例

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

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

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

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

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

圖4 表面氣動網格Fig.4 CFD grid of case 2

表1 算例1結構前四階模態(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標模算例得到了驗證，為求快速獲得結果，故而此處直接使用ROM進行顫振分析，得到此機翼顫振邊界如圖6所示。

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

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

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

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

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

3 結論

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

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

［1] 管 德.非定常空氣動力計算［M] .北京:北京航空航天大學出版社，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] 曾憲昂.模型降階技術在氣動彈性中的應用研究［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.