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

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

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

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

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

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

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

0 引 言

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

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

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

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

1 計算方法

1.1 控制方程

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

寫成半離散化的格式

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

1.2 諧波平衡方法

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

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

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

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

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

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

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

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

由此，通過時域控制方程式（1）可得頻域控制方程

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

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

其中D為頻域算子。

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

1.3 數(shù)值計算方法

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

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

1.4 邊界條件

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

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

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

2.1 NACA0012葉柵

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

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

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

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

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

2.2 平板葉柵

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

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

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

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

利用本文的諧波平衡方法對某壓氣機(jī)典型截面葉柵氣彈穩(wěn)定性進(jìn)行了分析。壓氣機(jī)葉型如圖5所示，其弦長c=0.0777m，安裝角γ=46.29°，稠度τ= 2.55。其計算條件為：進(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)振動，扭轉(zhuǎn)中心位于葉片前緣，振動扭轉(zhuǎn)角度幅度為α=3°，相鄰葉片振動相位差σ= 0°，折合頻率k=0.222。計算采用二階諧波格式，一個周期內(nèi)等時間距選擇五個時間層。

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

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

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

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

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

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

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

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

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

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

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

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

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

3 諧波階次對計算效率的影響

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

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

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

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

表1 不同階次諧波平衡方法計算耗時對比Table 1 Comparisonof computational effiency for different number of harmonics

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

4 結(jié) 論

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

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

（JC201245）

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

施永強(qiáng)，楊青真，黃秀全，等.一種快速模擬振蕩葉柵非定常流的數(shù)值方法［J］.空氣動力學(xué)學(xué)報，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.