一種快速模擬振蕩葉柵非定常流的數(shù)值方法

施永強(qiáng)，楊青真，黃秀全，郭 霄

（西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院，陜西西安 710072）

一種快速模擬振蕩葉柵非定常流的數(shù)值方法

施永強(qiáng)，楊青真，黃秀全，郭 霄

（西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院，陜西西安 710072）

為了提高葉輪機(jī)內(nèi)周期性非定常流數(shù)值計(jì)算效率，發(fā)展了一種計(jì)算振蕩葉柵非線性非定常流的諧波平衡方法，將周期性非定常流動(dòng)參數(shù)采用傅里葉級(jí)數(shù)表示，從而將傳統(tǒng)的定常流動(dòng)數(shù)值計(jì)算方法應(yīng)用于非定常流問題求解中，大大提高了非定常流計(jì)算效率。計(jì)算結(jié)果表明，發(fā)展的諧波平衡方法與線性方法、非線性時(shí)域法計(jì)算結(jié)果吻合較好，且其計(jì)算效率比時(shí)域法高1個(gè)數(shù)量級(jí)。另外，系統(tǒng)研究了諧波階次對(duì)計(jì)算效率與計(jì)算精度的影響，認(rèn)為在氣流未出現(xiàn)嚴(yán)重分離時(shí)，應(yīng)采用低階諧波平衡方法，以減少計(jì)算耗時(shí)。此外，從壓氣機(jī)振蕩葉柵非定常流場(chǎng)計(jì)算結(jié)果可以看出，葉片吸力面誘導(dǎo)渦的生成頻率與振蕩頻率并不一致，渦生成頻率要高于振蕩頻率。

振蕩葉柵；諧波平衡方法；非定常流；雙時(shí)間法；形狀修正法

0 引 言

葉輪機(jī)內(nèi)流場(chǎng)具有高度復(fù)雜的非定常性。傳統(tǒng)的二維軸對(duì)稱氣動(dòng)設(shè)計(jì)理論和全三維定常流場(chǎng)計(jì)算方法，都無(wú)法反映這一真實(shí)的流動(dòng)特性。用這種定常流體系處理具有強(qiáng)烈非定常特征的葉輪機(jī)設(shè)計(jì)問題時(shí)，從一般意義上講，會(huì)導(dǎo)致對(duì)非定常流“潛力”的某種忽略而使性能有所損失。葉輪機(jī)械中的葉片是彈性體，在氣動(dòng)力、結(jié)構(gòu)等因素影響下，會(huì)發(fā)生振動(dòng)。強(qiáng)的振動(dòng)現(xiàn)象會(huì)導(dǎo)致葉片在短時(shí)間內(nèi)斷裂，造成嚴(yán)重的后果。為了深入研究葉柵非定常流動(dòng)，探索葉輪機(jī)械強(qiáng)迫振動(dòng)、顫振的發(fā)作機(jī)理，就必須對(duì)振蕩葉柵非定常流場(chǎng)以及葉片所產(chǎn)生的周期性非定常氣動(dòng)力的性質(zhì)給予準(zhǔn)確、高效地模擬與預(yù)測(cè)。

目前，利用計(jì)算流體力學(xué)方法對(duì)振蕩葉柵非定常繞流及葉片所受氣動(dòng)力進(jìn)行數(shù)值模擬，進(jìn)而分析振蕩葉片與流場(chǎng)間的相互作用，已成為研究葉片強(qiáng)迫振動(dòng)與顫振的重要手段。在振蕩葉柵非定常流動(dòng)數(shù)值模擬研究方面，傳統(tǒng)的分析方法是時(shí)間推進(jìn)方法。但是，由于時(shí)間推進(jìn)方法計(jì)算過(guò)程中要采用統(tǒng)一的時(shí)間步長(zhǎng)推進(jìn)，因此計(jì)算工作量特別大，計(jì)算相當(dāng)耗時(shí)。1991年，Jameson提出了雙時(shí)間方法［1］，提高了非定常流動(dòng)計(jì)算效率。與時(shí)間推進(jìn)方法相比，這種方法計(jì)算速度是傳統(tǒng)方法的7～9倍［2］，因此被廣泛應(yīng)用于非定常流動(dòng)數(shù)值模擬。

雖然雙時(shí)間方法能夠加速非定常流動(dòng)數(shù)值計(jì)算，提高計(jì)算效率，但對(duì)于復(fù)雜的葉輪機(jī)械非定常流動(dòng)來(lái)說(shuō)，其計(jì)算工作量還是相當(dāng)巨大的。為此，國(guó)外學(xué)者在這方面做了大量的研究工作，發(fā)展了高效的頻域計(jì)算方法求解非定常流動(dòng)。Verdon［3］，Whitehead［4］和Hall［5］等人分別提出了模擬二維葉柵非定常流動(dòng)的時(shí)間線化方法；Hall和Lorence發(fā)展了三維線化非定常流動(dòng)計(jì)算方法［6］；2005年，Ekici和Hall采用時(shí)間線化方法分析了多級(jí)葉輪機(jī)械顫振問題［7］。時(shí)間線化方法在求解葉輪機(jī)械非定常流動(dòng)問題時(shí)，僅需要模擬單個(gè)通道的流動(dòng)狀況，且其計(jì)算時(shí)間與定常流動(dòng)計(jì)算時(shí)間相當(dāng)，因此，被廣泛應(yīng)用于葉輪機(jī)非定常流動(dòng)數(shù)值模擬中。然而，時(shí)間線化方法也存在著不可忽視的缺點(diǎn)——忽略了非定常流動(dòng)的非線性現(xiàn)象。為此，Hall提出了諧波平衡方法用于模擬葉柵中非線性非定常流動(dòng)，計(jì)算結(jié)果表明這種方法計(jì)算非定常流比傳統(tǒng)的時(shí)域方法快1到2個(gè)數(shù)量級(jí)，而且這種方法可以精確地計(jì)算非線性很強(qiáng)的非定常流場(chǎng)［8］。此后，Hall將這種方法應(yīng)用于各種葉輪機(jī)非定常流動(dòng)問題分析當(dāng)中，取 得了 一 系列 的 成果［9-12］。

本文分別采用諧波平衡方法和雙時(shí)間方法求解二維Navier-Stokes方程，對(duì)振蕩葉柵非定常粘性流動(dòng)進(jìn)行了數(shù)值模擬，對(duì)比分析了諧波平衡方法的計(jì)算效率與計(jì)算精度。同時(shí)，將偽邊界形狀修正方法引入非定常流動(dòng)計(jì)算程序中，大大縮短了振蕩葉柵存在相差時(shí)的非定常流計(jì)算時(shí)間。此外，采用諧波平衡方法分析了某壓氣機(jī)葉柵在非定常氣動(dòng)力作用下的氣動(dòng)彈性穩(wěn)定性，探討了振蕩葉柵相位差對(duì)氣彈穩(wěn)定性的影響。

1 計(jì)算方法

1.1 控制方程

在直角坐標(biāo)系下，對(duì)于運(yùn)動(dòng)網(wǎng)格，積分形式的二維雷諾時(shí)均Navier-Stokes方程可寫為：

寫成半離散化的格式

其中U=［ρ，ρu，ρv，e］T，E、F是無(wú)粘通量，Eν、Fν是粘性通量，u和v 分別是x 和y 方向的速度分量，e=分別為 兩個(gè)方向 的網(wǎng)格移動(dòng)速度。

1.2 諧波平衡方法

諧波平衡方法最早由Hall［8］提出，其主要思想是將周期性函數(shù)采用傅里葉級(jí)數(shù)表示。振蕩葉柵中的流體流動(dòng)是周期性函數(shù)，因此流場(chǎng)中的守恒變量可以展開成傅里葉級(jí)數(shù)的形式：

其中，Rn（x，y）、Un（x，y）、Vn（x，y）和En（x，y）分別是四個(gè)守恒變量的n階傅里葉系數(shù)，ω是周期函數(shù)的角頻率。

為了更清楚地表示式（3）中的變量，我們以U 表示任意變量，并以三角形式表示傅里葉級(jí)數(shù)，在忽略N 以上的高階項(xiàng)之后，即有

其中，A0、An和Bn是變量U 的傅里葉系數(shù)。

在一個(gè)周期內(nèi)，選取2N+1個(gè)時(shí)間點(diǎn)，式（4）可寫成矩陣形式

其中U*表示在2N+1個(gè)時(shí)間點(diǎn)上的守恒變量向量，E表示離散傅里葉逆變換算子，~U 表示守恒變量的傅里葉系數(shù)向量。由此式（5）可寫為

相反，守恒變量的傅立葉系數(shù)向量可通過(guò)傅立葉變換得到

其中，E-1表示離散傅立葉變換矩陣。

由此，通過(guò)時(shí)域控制方程式（1）可得頻域控制方程

這里，E*、F*分別為x、y方向的通量。

式（8）表示了2N+1個(gè)控制方程組，不同時(shí)刻的控制方程通過(guò)時(shí)間導(dǎo)數(shù)項(xiàng)耦合。時(shí)間導(dǎo)數(shù)項(xiàng)可表示為：

其中D為頻域算子。

由此，把時(shí)間導(dǎo)數(shù)項(xiàng)作為源項(xiàng)可得諧波平衡控制方程：

1.3 數(shù)值計(jì)算方法

方程（10）為一個(gè)周期內(nèi)2N+1個(gè)時(shí)間點(diǎn)上的控制方程，其方程個(gè)數(shù)為4×（2N+1）。為了求解諧波平衡方程，可以引入一個(gè)偽時(shí)間項(xiàng)，這樣方程就可以采用傳統(tǒng)的時(shí)間推進(jìn)方法進(jìn)行迭代求解。由此可得

方程（11）采用四步Runge-Kutta法在偽時(shí)間上推進(jìn)，并加入二階和四階人工黏性系數(shù)。采用當(dāng)?shù)貢r(shí)間步長(zhǎng)法、雙重網(wǎng)格技術(shù)和隱式殘差光順方法加速計(jì)算收斂。

1.4 邊界條件

邊界條件包括進(jìn)口邊界條件，出口邊界條件，壁面邊界條件和振蕩葉柵非定常流動(dòng)中的相差周期邊界條件。上游邊界給出總溫、總壓、氣流角，當(dāng)來(lái)流為超聲速時(shí)，用超聲葉柵唯一迎角原理，根據(jù)來(lái)流馬赫數(shù)確定進(jìn)口氣流角。即指定進(jìn)口馬赫數(shù)來(lái)代替指定進(jìn)口氣流角。下游邊界給定出口反壓，在軸向分速度為亞聲速時(shí)，出口壓力為反壓，其他參數(shù)外插得到；如果軸向分速為超聲速，氣動(dòng)參數(shù)外插得到。物面邊界條件：采用無(wú)滑移條件，對(duì)于振蕩葉柵，物面邊界條件即為u=umg，v=vmg。相差周期邊界采用偽 邊界形狀修正法［13］對(duì)周期邊界進(jìn)行處理，從而使得 在相鄰葉片存在相位差時(shí)，也能采用單通道進(jìn)行計(jì)算。

2 數(shù)值計(jì)算及結(jié)果分析

本文應(yīng)用諧波平衡方法分別對(duì) NACA0012葉柵、平板葉柵和某壓氣機(jī)典型截面葉柵非定常流場(chǎng)進(jìn)行了數(shù)值模擬。其中，NACA0012葉柵和壓氣機(jī)葉柵做純扭轉(zhuǎn)振動(dòng)，平板葉柵做純彎曲振動(dòng)。

2.1 NACA0012葉柵

由 NACA0012翼型組成的葉柵，其計(jì)算條件如下：進(jìn)口馬赫數(shù)M1=0.593，雷諾數(shù)Re=3.21×106，進(jìn)口氣流角αm=0°，柵距／弦長(zhǎng)=1.0。葉柵作扭轉(zhuǎn)振動(dòng)，扭轉(zhuǎn)幅度α=2°，振動(dòng)中心為葉片前緣，相鄰葉片振動(dòng)相位差σ=0°，折合頻率k=ωc／（2u∞）=0.2，式中ω為振動(dòng)角速度，c為葉片弦長(zhǎng)，u∞為上游遠(yuǎn)場(chǎng)邊界速度。

分別采用二階諧波平衡方法和雙時(shí)間法數(shù)值模擬振蕩葉柵非定常流動(dòng)，其計(jì)算結(jié)果如圖1～圖3所示。其中一階非定常壓力系數(shù)定義為：

式中，p1為壓力的一階諧波分量，ρ∞、u∞分別為上游遠(yuǎn)場(chǎng)氣體密度和速度，Am為扭轉(zhuǎn)振動(dòng)幅度。

圖1 NACA0012翼型上表面一階非定常壓力系數(shù)Fig.1 Unsteady pressure coefficient distribution on upper surface of airfoil NACA0012

圖2 NACA0012翼型升力系數(shù)隨時(shí)間變化曲線Fig.2 Lift coefficient variation of airfoil NACA0012 in a period

圖3 NACA0012翼型升力系數(shù)隨扭轉(zhuǎn)角度變化曲線Fig.3 Lift coefficient variation of airfoil NACA0012 with the torsion angle

由圖1可以看出，諧波平衡方法的計(jì)算結(jié)果與雙時(shí)間方法計(jì)算結(jié)果幾乎完全一致。另外與參考文獻(xiàn)［14］的計(jì)算結(jié)果吻合程度也較好，僅在前緣處略有不同。圖2為一個(gè)周期內(nèi)翼型升力系數(shù)的變化曲線。圖中也分別展示了諧波平衡方法、雙時(shí)間方法以及參考文獻(xiàn)［14］的線性方法計(jì)算結(jié)果。其中雙時(shí)間方法是截取了已經(jīng)計(jì)算穩(wěn)定的最后一個(gè)周期的計(jì)算結(jié)果。由圖可以看出，諧波平衡方法和雙時(shí)間方法的升力系數(shù)變化也幾乎完全一致，其升力系數(shù)的波動(dòng)范圍和相位都相互吻合。與參考文獻(xiàn)［14］相比，諧波平衡方法升力系數(shù)峰值為0.081，雙時(shí)間方法為0.0805，文獻(xiàn)［14］中為0.0763，稍有不同，但其相位一致。圖3為升力系數(shù)隨扭轉(zhuǎn)角度的變化規(guī)律，其諧波平衡方法和雙時(shí)間方法的計(jì)算結(jié)果也非常吻合。

2.2 平板葉柵

為了驗(yàn)證葉片間存在相差時(shí)的計(jì)算結(jié)果，本文以參考文獻(xiàn)［15］的平板葉柵為研究對(duì)象，驗(yàn)證本文的計(jì)算方法。平板葉柵的弦長(zhǎng)c=0.076m，柵距／弦長(zhǎng)= 1.3，進(jìn)口馬赫數(shù)M1=0.65，進(jìn)口氣流角αm=0°，迎角i=0°。葉柵做彎曲振動(dòng)，即垂直于弦線方向做平移振動(dòng)，其彎曲振動(dòng)幅度為Ba=1%c，振動(dòng)頻率f=250Hz（折合平率k=0.57），相鄰葉片振動(dòng)相位差σ=90°。

圖4 平板葉柵非定常壓差系數(shù)Fig.4 Unsteady pressure jump coefficient distribution for a flat plate cascade

本文采用形狀修正方法處理周期相差邊界，因此計(jì)算采用了一個(gè)通道。計(jì)算采用二階諧波平衡方法，每個(gè)周期設(shè)置五個(gè)時(shí)間層。圖4為計(jì)算所得的平板葉柵表面非定常壓差系數(shù)。非定常壓差系數(shù)定義為為翼型表面壓差一階諧波分量。圖中將諧波平衡方法的計(jì)算結(jié)果與文獻(xiàn)［15］中的利用線性理論計(jì)算出的結(jié)果進(jìn)行了對(duì)比，其吻合程度較好。

2.3 某壓氣機(jī)葉柵氣動(dòng)彈性分析

利用本文的諧波平衡方法對(duì)某壓氣機(jī)典型截面葉柵氣彈穩(wěn)定性進(jìn)行了分析。壓氣機(jī)葉型如圖5所示，其弦長(zhǎng)c=0.0777m，安裝角γ=46.29°，稠度τ= 2.55。其計(jì)算條件為：進(jìn)口馬赫數(shù)M1=0.75，進(jìn)氣角（與軸向夾角）αm=52°，進(jìn)口總壓p*=7042.13Pa，進(jìn)口總溫T*=232.26K，出口背壓pb=4998.9Pa。

圖5 壓氣機(jī)葉型Fig.5 The compressor arifoil

葉片發(fā)生扭轉(zhuǎn)振動(dòng)，扭轉(zhuǎn)中心位于葉片前緣，振動(dòng)扭轉(zhuǎn)角度幅度為α=3°，相鄰葉片振動(dòng)相位差σ= 0°，折合頻率k=0.222。計(jì)算采用二階諧波格式，一個(gè)周期內(nèi)等時(shí)間距選擇五個(gè)時(shí)間層。

為了研究葉片的氣彈穩(wěn)定性，定義非定常力矩系數(shù)：

圖6為計(jì)算得到的非定常壓差系數(shù)和非定常力矩系數(shù)虛部沿弦向的分布曲線。由圖6可知，此葉柵在相位差為0°時(shí)，氣動(dòng)力在葉表大部分區(qū)域做負(fù)功，僅在前緣和尾緣附近做正功。

圖6 壓氣機(jī)葉型非定常力矩系數(shù)Fig.6 Unsteady moment coefficient distribution for compressor airfoil

圖7為定常計(jì)算得到的葉柵通道馬赫數(shù)云圖和流線圖，圖8～圖11分別為0T，1／4T，1／2T 和3／4T時(shí)，葉柵通道馬赫數(shù)云圖和流線圖。由圖可見，在定常狀態(tài)時(shí)，葉型吸力面并沒有發(fā)生分離；當(dāng)葉片發(fā)生振動(dòng)時(shí)，在任何時(shí)刻，葉型吸力面總是存在分離渦。這說(shuō)明，葉片一旦發(fā)生振動(dòng)，很容易在其吸力面產(chǎn)生振動(dòng)誘導(dǎo)渦，從而致使壓氣機(jī)的氣動(dòng)性能下降。從非定常流動(dòng)瞬時(shí)流線圖可以看出，在葉型吸力面，一般

圖7 定常計(jì)算葉柵通道馬赫數(shù)云圖及流線圖Fig.7 Steady flow field around the cascades

圖8 0T時(shí)刻馬赫數(shù)分布云圖及流線圖Fig.8 Unteady flow field around the cascades（0T）

圖9 1／4T時(shí)刻馬赫數(shù)分布云圖及流線圖Fig.9 Unteady flow field around the cascades（1／4T）

圖10 1／2T時(shí)刻馬赫數(shù)分布云圖及流線圖Fig.10 Unteady flow field around the cascades（1／2T）

圖11 3／4T時(shí)刻馬赫數(shù)分布云圖及流線圖Fig.11 Unteady flow field around the cascades（3／4T）

存在兩個(gè)分離渦，在1／2T 時(shí)刻，存在三個(gè)分離渦。這說(shuō)明，振動(dòng)誘導(dǎo)的分離渦產(chǎn)生的頻率大于振動(dòng)的頻率，在一個(gè)振動(dòng)周期，可產(chǎn)生若干個(gè)分離渦。對(duì)不同時(shí)刻的流線圖進(jìn)行對(duì)比可以看出，在瞬時(shí)迎角較小的1／4T時(shí)刻，吸力面依然存在兩個(gè)分離渦，而在0T 時(shí)刻，吸力面上一個(gè)分離渦已經(jīng)運(yùn)動(dòng)到尾緣處，即將脫落，其流場(chǎng)品質(zhì)比1／4T 時(shí)刻稍好。因此，在振動(dòng)情況下，瞬時(shí)迎角較小，并不代表其流場(chǎng)品質(zhì)較好，流場(chǎng)結(jié)構(gòu)與分離渦產(chǎn)生的頻率有很大關(guān)系。

為了對(duì)壓氣機(jī)葉柵的氣彈特性進(jìn)行全面分析，本文計(jì)算了相鄰葉片振動(dòng)相位差σ=30°，60°，90°，120°，150°，180°，210°，240°，270°，300°和330°時(shí)的非定常流場(chǎng)，計(jì)算得到了不同振動(dòng)相位差時(shí)葉柵非定常力矩系數(shù)CM，如圖12所示。

由圖12可以看出，本葉柵氣彈失穩(wěn)可能性最大的點(diǎn)并不在σ=0°處。因此，在葉輪機(jī)械氣彈穩(wěn)定性分析時(shí)，需對(duì)其不同振動(dòng)相差情況進(jìn)行分析，以便準(zhǔn)確把握其氣動(dòng)彈性特性。此外，從圖中可知，本葉柵在振動(dòng)相位差位于18°～166°之間時(shí)，有可能發(fā)生氣彈失穩(wěn)。

圖12 壓氣機(jī)葉柵非定常力矩系數(shù)隨相鄰葉片振動(dòng)相位差變化情況Fig.12 Unsteady moment coefficient of compressor airfoil at different interblade phase angles

3 諧波階次對(duì)計(jì)算效率的影響

在計(jì)算非定常流動(dòng)時(shí)，本文所采用的計(jì)算方法一個(gè)周期內(nèi)最少可以只設(shè)置三個(gè)時(shí)間層來(lái)計(jì)算，即一階諧波格式。為了研究諧波平衡方法階次對(duì)計(jì)算效率的影響，本文分別針對(duì)2.1節(jié)的NACA0012組成的葉柵算例采用不同階次的諧波平衡方法進(jìn)行計(jì)算。

圖13是分別采用一階、二階、三階、四階諧波平衡方法計(jì)算得到的葉柵表面非定常壓差系數(shù)。由圖可見，采用不同階次的方法計(jì)算結(jié)果幾近完全相同。

圖13 不同階次諧波平衡方法非定常壓差系數(shù)比較Fig.13 Comparison of unsteady pressure jump coefficient for different number of harmonics

表1是不同階次諧波平衡方法的計(jì)算耗時(shí)對(duì)比。在這里假設(shè)定常計(jì)算所用時(shí)間為1，不同階次諧波平衡方法計(jì)算耗時(shí)與定常計(jì)算耗時(shí)的比值即為無(wú)量綱計(jì)算耗時(shí)。定常計(jì)算與非定常計(jì)算取同樣的計(jì)算步數(shù)，均推進(jìn)4000步。

表1 不同階次諧波平衡方法計(jì)算耗時(shí)對(duì)比Table 1 Comparisonof computational effiency for different number of harmonics

由表1可見，前三階格式的諧波平衡方法計(jì)算耗時(shí)與定常計(jì)算處于同一量級(jí)。而其他學(xué)者認(rèn)為，采用時(shí)域方法的計(jì)算耗時(shí)要比諧波平衡方法的計(jì)算耗時(shí)高一到兩個(gè)量級(jí)。由于計(jì)算所采用的收斂評(píng)判角度不同，因此無(wú)法做出準(zhǔn)確對(duì)比。但從本文計(jì)算過(guò)程來(lái)看，單通道計(jì)算時(shí)，諧波平衡方法的計(jì)算效率確實(shí)相當(dāng)高，時(shí)域方法計(jì)算耗時(shí)約為二階諧波平衡方法計(jì)算耗時(shí)的20～30倍左右，可見二階諧波平衡方法的計(jì)算效率要比時(shí)域方法的計(jì)算效率高一個(gè)量級(jí)。如果存在相位差時(shí)，采用多通道的時(shí)域計(jì)算方法將因?yàn)橛?jì)算區(qū)域成倍增加而使得計(jì)算耗時(shí)也幾乎成倍增加。因此，在存在相位差時(shí)，采用諧波平衡方法的計(jì)算效率比時(shí)域法的計(jì)算效率可以高達(dá)兩個(gè)量級(jí)。另外，從不同階次的諧波平衡方法耗時(shí)對(duì)比來(lái)看，階次每高一階，其耗時(shí)平均增加2.39。因此，在非定常流計(jì)算時(shí)，應(yīng)根據(jù)計(jì)算對(duì)象的不同，選擇適當(dāng)?shù)闹C波平衡方法，既保證精度要求，又能提高計(jì)算效率。

4 結(jié) 論

本文發(fā)展了基于諧波平衡方法計(jì)算振蕩葉柵非定常流動(dòng)的程序，通過(guò)算例對(duì)這種方法進(jìn)行了驗(yàn)證，得出以下結(jié)論：（1）通過(guò)諧波平衡方法與雙時(shí)間方法的計(jì)算結(jié)果對(duì)比表明，兩種計(jì)算方法得到的翼型表面一階非定常壓力系數(shù)分布和升力系數(shù)吻合的非常好，從而驗(yàn)證了程序的可靠性和諧波平衡方法的準(zhǔn)確性。（2）通過(guò)計(jì)算對(duì)比證明諧波平衡方法是一種非常高效的非定常流計(jì)算方法，本文二階諧波平衡方法的計(jì)算效率比雙時(shí)間方法高一個(gè)量級(jí)，因此可以廣泛應(yīng)用于周期性非定常流動(dòng)數(shù)值模擬當(dāng)中。（3）對(duì)于周期性較強(qiáng)的非定常流動(dòng)，可以采用較少的階次就可以滿足計(jì)算精度，可大大提高計(jì)算效率。（4）從壓氣機(jī)葉柵氣彈穩(wěn)定性分析結(jié)果來(lái)看，氣彈失穩(wěn)點(diǎn)不一定在0°相位差處。因此，在氣彈穩(wěn)定性分析時(shí)，要對(duì)不同相差情況進(jìn)行分析。

［1］ JAMESON A.Time dependent calculations using multigrid，with applications to unsteady flows past airfoils and wings［A］.AIAA 10th Computational Fluid Dynamics Conference［C］.Honolulu，Hawaii，1991.

［2］ HU Y C，ZHOU X H.Application of dual-time stepping method to calculating unsteady viscous flow around oscillating cascade［J］.Journal of Northwestern Polytechnical University，2003，21（3）：269-272.（In Chinese）胡運(yùn)聰，周新海.振蕩葉柵粘性非定常繞流計(jì)算的雙時(shí)間法［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2003，21（3）：269-272.

［3］ VERDON J M，CAPAR J R.Development of a linear unsteady aerodynamic analysis for finite-deflection subsonic cascades［J］.AIAA Journal，1982，20（9）：1259-1267.

［4］ WHITEHEAD D S.A finite-element solution of unsteady twodimensional flow in cascades［J］.International Journal for Numerical Methods in Fluids，1990，10（1）：13-34.

［5］ HALL K C，CLARK W S.Linearized euler predictions of unsteady aerodynamic loads in cascades［J］.AIAA Journal，1993，31（3）：540-550.

［6］ HALL K C，LORENCE C B.Calculation of three-demension unsteady flows in turbomachinery using the linearized harmonic euler equations［J］.Journal of Turbomachinery-Transactions of ASME，1993，15（4）：800-809.

［7］ EKICI K，VOYTOVYCH D M，HALL K C.Time-linearized navier-stokes analysis of flutter in multistage turbomachines ［A］.43rd AIAA Aerospace Sciences Meeting and Exhibit［C］.Reno，NV，2005.

［8］ HALL K C，THOMAS J P，Cl ARK W S.Computation of unsteady nonlinear flows in cascades using a harmonic balance technique［J］.AIAA Journal，2002，40（5）：879-886.

［9］ EKICI K，HALL K C.Nonlinear analysis of unsteady flows in multistage turbomachines using the harmonic balance technique ［A］.44th AIAA Aerospace Sciences Meeting and Exhibit［C］.Reno，NV，2006.

［10］EKICI K，HALL K C.Nonlinear analysis of unsteady flows in multistage turbomachines using harmonic balance［J］.AIAA Journal，2007，45（5）：1047-1057.

［11］EKICI K，KIELB R E，HALL K C.The effect of aerodynamic asymmetries on turbomachinery flutter［A］.47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition［C］.2009，Orlando，F(xiàn)lorida.

［12］THOMASJ P，CUSTER C H，DOWELL E H，et al.F-16 fighter aeroelastic computations using a harmonic balance implementation of the OVERFLOW 2 Flow Solver［A］.51st AIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics，and Materials Conference［C］.2010，Orlando，F(xiàn)lorida.

［13］YANG Q Z，SHI Y Q，HUANG X Q，et al.Numerical analysis of three-dimensional unsteady flow around oscillating cascades ［J］.ACTA Aerodynamica Sinica，2006，24（1）：62-65.（In Chinese）楊青真，施永強(qiáng)，黃秀全，等.振蕩葉柵三維非定常流動(dòng)數(shù)值分析［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2006，24（1）：62-65.

［14］HUFF D L.Numerical simulations of unsteady，viscous，transonic flow over isolated and cascaded airfoils using a deforming grid［A］.AIAA 19th Fluid Dynamics，Plasma Dynamics and Lasers Conference［C］.1987／Honolulu，Hawaii.

［15］HE L.An Euler solution for unsteady flows around oscillating blades ［J］.Transactions of American Society of Mechianical Engineers：Journal of Turbomachinery，1990，112（4）：714-722.

A numerical approach for fast simulation of unsteady flow around oscillating cascades

SHI Yongqiang，YANG Qingzhen，HUANG Xiuquan，GUO Xiao

（School of Power&Energy，Northwestern Polytechnical University，Xi’an 710072，China）

A harmonic balance method for the analysis of nonlinear unsteady flow around oscillating cascades is developed to improve computational efficiency for the simulation of periodic flow in turbomachinery.The periodic unsteady flow parameters can be represented by a Fourier series.In this way，the unsteady problem can be solved using conventional computational fluid dynamic method for steady problem.The analysis exploits the fact that，calculated results by harmonic balance are in good agreement with those from dual-time stepping method and linear method，and present method is highly efficient；it is about one order of magnitude faster than conventional dual-time stepping method in the present computation.The effects of the number of harmonics are also studied.It shows that lower order harmonic balance methods are sufficient to simulate periodic unsteady flow accurately except the flow fields including serious separated domain.So，lower order harmonic balance methods could be adopted to analysis unsteady flow in turbomachinery to reduce the computational cost.In addition，the characteristic of unsteady flow around a compressor cascade is investigated，and the mechanism of generation of separated vortex and the vortex shedding frequency are also analyzed.

oscillating cascade；harmonic balance method；unsteady flow；dual-time stepping method；shape correction method

V221.3

Adoi：10.7638／kqdlxxb-2012.0146

0258-1825（2014）04-0481-07

2012-09-06；

2013-02-25

國(guó)家自然科學(xué)基金（50376053，50706039，11302179）；陜西省自然科學(xué)基礎(chǔ)研究基金（2013JQ1009）；西北工業(yè)大學(xué)基礎(chǔ)研究基金

（JC201245）

施永強(qiáng)（1980-），男，甘肅平?jīng)鋈?，副教授，博士，研究方向：葉輪機(jī)氣動(dòng)熱力學(xué).E-mail：yqshi＠nwpu.edu.cn

施永強(qiáng)，楊青真，黃秀全，等.一種快速模擬振蕩葉柵非定常流的數(shù)值方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（4）：481-487.

10.7638／kqdlxxb-2012.0146. SHI Y Q，YANG Q Z，HUANG X Q，et al.A numerical approach for fast simulation of unsteady flow around oscillating cascades［J］.ACTA Aerodynamica Sinica，2014，32（4）：481-487.