張 珺，李立州，羅 驍，李彩霞，原梅妮

（1.太原學院數(shù)學系，太原030001；2.中北大學機電工程學院，太原030051）

0 引言

航空發(fā)動機中存在氣流導向的靜葉和對外作功的動葉。當動、靜葉相對轉(zhuǎn)動時，上游尾流會使下游葉片氣動力周期性振蕩[1-3]，使其受迫振動、發(fā)生疲勞破壞[4-5]；近年來發(fā)現(xiàn)上游尾流甚至可以改變?nèi)~片的顫振特性[6-8]。因此，研究動、靜葉干涉下葉片的氣動力彈性振動對發(fā)動機設(shè)計有著重要意義。目前，對動、靜葉干涉下葉片的氣動力彈性振動的研究主要采用數(shù)值和試驗方法，效率極低，工程應用不便[9-15]。特別是對于需要反復進行氣動彈性分析的發(fā)動機多學科優(yōu)化、可靠性分析和疲勞設(shè)計優(yōu)化，直接采用數(shù)值和試驗方法的工作量巨大。亟需開發(fā)精確、高效的方法。

氣動力降階模型(aerodynamic ROM)是1種描述擾動和葉片氣動力之間關(guān)系的簡化數(shù)學模型[10]。Silva[10]提出了基于Volterra級數(shù)的非線性氣動力降階模型，用于描述機翼振動引起的氣動力變化；張偉偉[11-12]等基于Volterra級數(shù)模型研究了葉柵的顫振；Liou[13]等用Volterra級數(shù)法研究了跨聲速葉片的顫振；Ekici[14]等用諧波平衡方法研究了葉片的顫振，認為所得結(jié)果與勢流理論結(jié)果符合較好；Ashcroft[15]等用諧波平衡法研究了2維壓氣機葉柵的亞聲速和跨聲速顫振特性，可以準確預測顫振。

但現(xiàn)有的氣動力降階模型的研究集中在葉片和機翼顫振方面，對動、靜葉干涉的葉片氣動彈性振動問題沒有涉及。本文將基于Volterra級數(shù)的氣動力降階模型的方法擴展到動、靜葉干涉問題，建立尾流激勵的葉片氣動力降階模型，并研究穩(wěn)態(tài)條件和識別信號的選取對該氣動力降階模型精度的影響。該氣動力降階模型可以在葉片高低周疲勞設(shè)計、葉片氣動彈性分析和優(yōu)化、葉片壽命可靠性分析中快速提供氣動力荷載。

1 基于Volterra級數(shù)的尾流激勵的葉片氣動力降階模型

1.1 時變系統(tǒng)的Volterra級數(shù)降階模型

任意時變系統(tǒng)Ψ可以用Volterra級數(shù)表示[10]

對于均勻采樣的時變系統(tǒng)，其離散的Volterra級數(shù)模型可以寫為

式中：X[n]=X（t）|t=nΔΤ=X（nΔΤ），為系統(tǒng)的輸入向量；Y[n]=Y（t）|t=nΔΤ=Y（nΔΤ），為系統(tǒng)的輸出向量；n=0,1,…,N，為離散的時間；ΔΤ為時間步長；t為時間；H0為系統(tǒng)穩(wěn)態(tài)響應矩陣；H1、H2、Hl分別為 1、2、l階核函數(shù)矩陣。

由于Volterra級數(shù)的2階以上核函數(shù)的數(shù)量大，核函數(shù)識別的工作量巨大，因此常采用1階核函數(shù)的Volterra級數(shù)降階模型。

1.2 基于Volterra級數(shù)的尾流激勵的葉片氣動力降階模型

基于以上Volterra級數(shù)降階模型，本文建立尾流激勵對葉片氣動力影響的降階模型。以2維葉片為例介紹氣動力降階模型的方法，葉片如圖1所示。由于上、下游葉片相對轉(zhuǎn)動，上游尾流以速度w在進口處移動。上游尾流移動使下游流場和葉片氣動力周期性變化。

該上游尾流激勵的葉片流場可以看作進口條件隨時間變化的系統(tǒng)，進口條件可以用進口總壓p和進氣角β表示，是時間t和坐標y的函數(shù)。因此，尾流激勵的葉片氣動力系統(tǒng)可以表示為

圖1 尾流激勵的葉片

式中：Cl（t）、Cm（t）、Cd（t）分別為葉片的升力、力矩、阻力。

進口總壓和進氣角的變化是由于上游尾流沿進口移動引起的，是沿著進口移動的行波，因此進口各點的總壓和進氣角可以用進口已知點y0的總壓和進氣角表示，即和 β（y,則尾流激勵的葉片氣動力系統(tǒng)可以表示為

用1階Volterra級數(shù)展開該系統(tǒng)，得到上游尾流激勵的氣動力降階模型

式中：H0、U0和V0分別為穩(wěn)態(tài)升力、力矩和阻力；分別為進口總壓變化引起的升力、力矩和阻力變化的1階核函數(shù)；分別為進氣角變化引起的升力、力矩和阻力變化的1階核函數(shù)。

1.3 Volterra級數(shù)氣動力降階模型中核函數(shù)的識別方法

氣動力降階模型的核函數(shù)可以用階躍、脈沖和噪聲等信號識別。用階躍信號辨識出的核函數(shù)比較穩(wěn)定且能反映系統(tǒng)的部分非線性特性，因此在實踐中廣泛采用。階躍信號識別氣動力降階模型的核函數(shù)的方法的步驟是：首先獲得葉片的穩(wěn)態(tài)流場和氣動力；然后在穩(wěn)態(tài)流場上增加單位階躍的擾動，求在該階躍擾動下葉片氣動力的響應；最后用如下公式計算降階模型的核函數(shù)

式中：y[n]為由單位階躍擾動引起的葉片氣動力響應；

穩(wěn)態(tài)流場的選擇應當使氣動力的穩(wěn)態(tài)值接近尾流激勵下氣動力波動的平均值，而且階躍信號的選擇應當符合小擾動假設(shè)。氣動力響應可以是氣動力試驗的結(jié)果，也可以是CFD仿真的結(jié)果。將識別好的核函數(shù)代入氣動力降階模型（式（5）），就可預測任意尾流下該葉片的氣動力。

2 算例驗證

為驗證本文方法，以圖2（a）中下游葉片為研究對象，建立上游尾流對該葉片氣動力影響的降階模型。取流場進口的總壓和進氣角為輸入，取下游中間葉片的氣動力為輸出。尾流的總壓和進氣角的波形數(shù)據(jù)取自如圖2所示流場上下游葉片流場交界面。上游葉片以10 m/s的速度沿著進口向上移動。

圖2 尾流激勵的葉片流場

采用fluent求解，空氣，理想氣體，Spallart-Allmaras黏性模型，無滑移壁面，葉片進、出口如圖2所示。出口壓力為101325 Pa，溫度為300 K。求解時，先給定1個穩(wěn)態(tài)進口總壓和穩(wěn)態(tài)進氣角，讓進口所有點的總壓和進氣角都為該值，求得穩(wěn)態(tài)流場。在穩(wěn)態(tài)流場求解的基礎(chǔ)上，用fluent的udf程序調(diào)整不同時刻流場進口的壓力和進氣角，以此實現(xiàn)尾流波形在進口的移動。

首先用該fluent氣動模型辨識氣動力降階模型的核函數(shù)。為簡化表述，這里只介紹總壓變化引起葉片氣動力變化的核函數(shù)其它核函數(shù)可采用相同方法識別，不再贅述。識別時，首先選定1個穩(wěn)態(tài)進口總壓，讓進口所有點的總壓都為該值，其它條件不變。穩(wěn)態(tài)進口總壓選取3種：111325、118000和 120500 Pa。用fluent求解穩(wěn)態(tài)進口總壓下的葉片氣動力。3種進口下葉片的升力H0分別為 -435、-702、-797 N；阻力 V0分別為 145.96014、246.88905、285.71491 N；力矩 U0分別為 -6.4133183、-10.636482、-12.217384 N·m。

在各穩(wěn)態(tài)進口總壓上增加1個階躍擾動，計算擾動引起的氣動力變化。采用3種幅值的階躍擾動信號1000、2000和3000 Pa。對穩(wěn)態(tài)總壓和階躍信號進行組合配對，分別計算葉片的升力和核函數(shù)。歸一化的核函數(shù)（單位總壓變化下葉片升力變化）如圖3所示。從圖中可見，算例中不同總壓穩(wěn)態(tài)值和階躍信號幅值下核函數(shù)基本一致。說明總壓穩(wěn)態(tài)值和信號幅值對核函數(shù)識別結(jié)果的影響不顯著。這也進一步說明在小擾動條件下流場系統(tǒng)是弱非線性系統(tǒng)[10]。

圖3 核函數(shù)識別結(jié)果

將圖3的核函數(shù)代入式（5）得到尾流激勵的葉片氣動力降階模型。用Matlab編寫氣動力降階模型的程序可以預測任意總壓變化引起的葉片氣動力變化。用該氣動力降階模型計算總壓波動（如圖4所示）引起的葉片氣動力。該進口總壓穩(wěn)態(tài)值為120500 Pa，穩(wěn)態(tài)升力為-797 N，穩(wěn)態(tài)阻力為285.71491 N，穩(wěn)態(tài)力矩為-12.217384 N·m。該總壓波動的數(shù)據(jù)取自圖2（a）中上下游葉片流場界面處。

圖4 尾流進口總壓波形

用降階模型計算得到的總壓波動引起的葉片氣動力如圖5所示。圖中減去了穩(wěn)態(tài)氣動力。圖中第1個數(shù)字代表辨識該核函數(shù)時的穩(wěn)態(tài)總壓，第2個數(shù)字代表階躍信號的幅值。為方便比較，圖中也給出了用fluent計算出的葉片氣動力。從圖中可見，所建立的降階模型對尾流激勵的葉片氣動力描述較好，得到的葉片氣動力的波形和振幅與用CFD模型得到的結(jié)果基本一致；用各核函數(shù)得到氣動力結(jié)果基本相同，說明穩(wěn)態(tài)總壓和階躍信號對核函數(shù)的精度影響不大。

圖5 總壓波動下葉片的升力

用降階模型計算出的葉片氣動力的誤差如圖6所示。該誤差是降階模型計算出的氣動力與CFD模型計算出的氣動力的差的絕對值。從圖中可見，在初始沖擊階段誤差較小，在周期振蕩階段誤差較大。需要進一步提高周期振蕩階段的降階模型的精度。

另外需要說明的是：識別本例中氣動力降階模型的核函數(shù)所用CFD計算時間為6 h，降階模型計算時間小于1 s，而尾流激勵下葉片氣動力CFD計算時間為32 h。

圖6 葉片的氣動力的誤差

3 總結(jié)

本文基于Volterra級數(shù)，建立了尾流激勵的葉片氣動力降階模型，并研究了穩(wěn)態(tài)條件和階躍信號幅值選取對尾流激勵的氣動力降階模型辨識精度的影響。結(jié)果表明：所建立的Volterra級數(shù)的氣動力降階模型對尾流激勵的葉片氣動力描述較好，得到的葉片氣動力的波形和振幅與用CFD模型得到的結(jié)果基本一致，所建立的氣動力降階模型能夠正確表征尾流對葉片的激勵作用；不同穩(wěn)態(tài)條件和階躍信號幅值下核函數(shù)和葉片氣動力響應基本相同，穩(wěn)態(tài)條件和階躍信號幅值對核函數(shù)的影響不顯著；該氣動力降階模型可以快速估計尾流激勵下葉片的氣動力，非常適合在葉片高低周疲勞設(shè)計、葉片氣動彈性分析和優(yōu)化、葉片壽命可靠性分析中快速提供氣動力荷載，而不需要反復進行CFD計算。

[1]趙養(yǎng)正,劉前智,廖明夫.非均勻柵距對壓氣機葉片非定常氣動力的影響[J].推進技術(shù)，2007，28（2）：167-170.ZHAO Yangzheng,LIU Qianzhi,LIAO Mingfu.Effects of uneven blade spacing on unsteady blade forcing in an axial compressor[J].Journal of Propulsion Technology,2007，28（2）：167-170.(in Chinese)

[2]白志剛,劉前智.非均勻柵距對軸流壓氣機流動影響的數(shù)值研究[J].航空工程進展，2010，1（1）：90-94.BAI Zhigang,LIU Qianzhi.Numerical study of the effects of uneven blade spacing on the flow in an axial compressor[J].Advances in Aeronautical Science and Engineering,2010，1（1）：90-94.(in Chinese)

[3]肖大啟，鄭赟，楊慧.軸向間距對轉(zhuǎn)子葉片氣動激勵的影響[J].航空動力學報，2012，27（10）:2307-2313.XIAO Daqi,ZHENG Yun,YANG Hui.Influence of axial spacing on aerodynamic excitation of rotor blade[J].Journal of Aerospace Power,2012，27（10）:2307-2313.(in Chinese)

[4]王梅,江和甫,呂文林.在尾流激振情況下葉片振動應力預估計技術(shù)[J].航空動力學報，2007，22（4）：608-613.WANG Mei,JIANG Hefu,LYU Wenlin.Method to predict the blade vibration stress induced by wake flow[J].Journal of Aerospace Power,2007，22（4）：608-613.(in Chinese)

[5]Lau Y L,Leung R C K,So R M C.Vortex-induced vibration effect on fatigue life estimate of turbine blades[J].Journal of Sound and Vibration,2007,307:698-719.

[6]張陳安,葉正寅,劉鋒，等.進口導流葉片對轉(zhuǎn)子葉片顫振特性的影響[J].推進技術(shù)，2010，31（3）：335-339.ZHANG Chen'an,YE Zhengyin,LIU Feng,et al.Investigations on flutter characteristics of rotor blade with IGV/fan interactions[J].Journal of Propulsion Technology,2010，31（3）：335-339.(in Chinese)

[7]Liu Z Y.Numerical researches on aeroelastic problem of a rotor due to IGV/Fan Interaction[R].AIAA-2009-865.

[8]趙振華，陳偉，陳敏.對轉(zhuǎn)渦輪低壓轉(zhuǎn)子葉片的流固耦合數(shù)值分析[J].機械科學與技術(shù)，2012，31（9）：1424-1428.ZHAO Zhenhua,CHEN Wei,CHEN Min.Numerical analysis of low pressure rotorblade using fluid-structure interaction method in counter-rotating turbine[J].Mechanical Science and Technology for Aerospace Engineering,2012，31（9）：1424-1428.(in Chinese)

[9]Clark W S,Hall K C.A Time-linearised Navier-Stokes analysis of stall flutter[J].ASME Journal of Turbomachinery,2000,122（3）:467-476.

[10]Walter A Silva.Reduced order models based on linear and nonlinear aerodynamic impulse responses[J].AIAA Journal，1999，233：1-11.

[11]張偉偉，蘇丹，張陳安，等.一種基于CFD的葉輪機非定常氣動力組合建模方法[J].推進技術(shù)，2012，33（1）：37-41.ZHANG Weiwei,SU Dan,ZHANG Chen'an,et al.A CFD-based compositional methodology of unsteady aerodynamic modeling for turbomachinery[J].Journal of Propulsion Technology,2012，33（1）：37-41.(in Chinese)

[12]Dan Su,Weiwei Zhang,Mingsheng Ma,et al.An efficient coupled method of cascade flutter investigation based on reduced order model[R].AIAA-2013-3618.

[13]MengSing Liou,Weigang Yao.Flutter analysis for turbomachinery using volterra series[R].ASME 2014-GT-25474.

[14]Kivanc Ekici,Robert E Kielb,Kenneth C Hall.The effect of aerodynamic asymmetries on turbomachinery flutter[J].Journal of Fluids and Structures，2013（36）：1-17.

[15]Graham Ashcroft,Christian Frey,Hans-Peter Kersken.On the development of a harmonic balance method for aeroelastic analy-sis[C]//the 11th World Congress on Computational Mechanics(WCCM XI),Barcelona,Spain,2014.