識(shí)別信號(hào)對(duì)尾流激勵(lì)的葉片氣動(dòng)力Volterra降階模型的影響

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

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

0 引言

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

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

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

1 基于Volterra級(jí)數(shù)的尾流激勵(lì)的葉片氣動(dòng)力降階模型

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

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

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

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

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

1.2 基于Volterra級(jí)數(shù)的尾流激勵(lì)的葉片氣動(dòng)力降階模型

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

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

圖1 尾流激勵(lì)的葉片

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

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

用1階Volterra級(jí)數(shù)展開(kāi)該系統(tǒng)，得到上游尾流激勵(lì)的氣動(dòng)力降階模型

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

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

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

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

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

2 算例驗(yàn)證

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

圖2 尾流激勵(lì)的葉片流場(chǎng)

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

首先用該fluent氣動(dòng)模型辨識(shí)氣動(dòng)力降階模型的核函數(shù)。為簡(jiǎn)化表述，這里只介紹總壓變化引起葉片氣動(dòng)力變化的核函數(shù)其它核函數(shù)可采用相同方法識(shí)別，不再贅述。識(shí)別時(shí)，首先選定1個(gè)穩(wěn)態(tài)進(jìn)口總壓，讓進(jìn)口所有點(diǎn)的總壓都為該值，其它條件不變。穩(wěn)態(tài)進(jìn)口總壓選取3種：111325、118000和 120500 Pa。用fluent求解穩(wěn)態(tài)進(jìn)口總壓下的葉片氣動(dòng)力。3種進(jìn)口下葉片的升力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)進(jìn)口總壓上增加1個(gè)階躍擾動(dòng)，計(jì)算擾動(dòng)引起的氣動(dòng)力變化。采用3種幅值的階躍擾動(dòng)信號(hào)1000、2000和3000 Pa。對(duì)穩(wěn)態(tài)總壓和階躍信號(hào)進(jìn)行組合配對(duì)，分別計(jì)算葉片的升力和核函數(shù)。歸一化的核函數(shù)（單位總壓變化下葉片升力變化）如圖3所示。從圖中可見(jiàn)，算例中不同總壓穩(wěn)態(tài)值和階躍信號(hào)幅值下核函數(shù)基本一致。說(shuō)明總壓穩(wěn)態(tài)值和信號(hào)幅值對(duì)核函數(shù)識(shí)別結(jié)果的影響不顯著。這也進(jìn)一步說(shuō)明在小擾動(dòng)條件下流場(chǎng)系統(tǒng)是弱非線性系統(tǒng)[10]。

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

將圖3的核函數(shù)代入式（5）得到尾流激勵(lì)的葉片氣動(dòng)力降階模型。用Matlab編寫(xiě)氣動(dòng)力降階模型的程序可以預(yù)測(cè)任意總壓變化引起的葉片氣動(dòng)力變化。用該氣動(dòng)力降階模型計(jì)算總壓波動(dòng)（如圖4所示）引起的葉片氣動(dòng)力。該進(jìn)口總壓穩(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。該總壓波動(dòng)的數(shù)據(jù)取自圖2（a）中上下游葉片流場(chǎng)界面處。

圖4 尾流進(jìn)口總壓波形

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

圖5 總壓波動(dòng)下葉片的升力

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

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

圖6 葉片的氣動(dòng)力的誤差

3 總結(jié)

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

[1]趙養(yǎng)正,劉前智,廖明夫.非均勻柵距對(duì)壓氣機(jī)葉片非定常氣動(dòng)力的影響[J].推進(jìn)技術(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]白志剛,劉前智.非均勻柵距對(duì)軸流壓氣機(jī)流動(dòng)影響的數(shù)值研究[J].航空工程進(jìn)展，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]肖大啟，鄭赟，楊慧.軸向間距對(duì)轉(zhuǎn)子葉片氣動(dòng)激勵(lì)的影響[J].航空動(dòng)力學(xué)報(bào)，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]王梅,江和甫,呂文林.在尾流激振情況下葉片振動(dòng)應(yīng)力預(yù)估計(jì)技術(shù)[J].航空動(dòng)力學(xué)報(bào)，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]張陳安,葉正寅,劉鋒，等.進(jìn)口導(dǎo)流葉片對(duì)轉(zhuǎn)子葉片顫振特性的影響[J].推進(jìn)技術(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]趙振華，陳偉，陳敏.對(duì)轉(zhuǎn)渦輪低壓轉(zhuǎn)子葉片的流固耦合數(shù)值分析[J].機(jī)械科學(xué)與技術(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ī)非定常氣動(dòng)力組合建模方法[J].推進(jìn)技術(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.