謝 亮,徐 敏,李 杰,蔡天星
(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)過程中的動載荷情況。
將結(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)的右端項QF,即非定常氣動載荷,采用CFD方法求解。積分形式的Euler方程為:
圖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)。
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)定性來求顫振邊界。
飛行器動態(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ù)值差分求解得到。
分別采用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
機翼幾何形狀及有限元模型見圖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)載荷(包括氣動載荷與慣性載荷)的有效手段。
采用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.