張恒銘, 唐弋棣, 衛(wèi)海洋, 林彬彬

(1.中國民航飛行學(xué)院航空工程學(xué)院，廣漢 618300；2.L.A.Turbine渦輪機(jī)械公司，瓦倫西亞 91355；3.北京動力機(jī)械研究所，北京 100074)

葉輪機(jī)中的流動現(xiàn)象具有強(qiáng)烈的非定常性，而引起這種非定常性的因素很多。其中，由葉排相對運(yùn)動和流動參數(shù)周向不均勻性所引起的非定常排間干擾是研究的熱點(diǎn)。排間干擾主要包含了周向靜壓不均勻所引起的勢擾動和尾跡擾動。這些擾動對壓氣機(jī)的性能有顯著的影響。中外的學(xué)者在非定常流動方面開展了大量的研究[1-5]。他們的研究表明排間干擾對葉輪機(jī)性能產(chǎn)生了不可忽略的影響，例如尾跡作用于下游流場時不但會與下游的邊界層、尾緣渦和端壁二次流等發(fā)生相互作用，且自身還會在下游的輸運(yùn)過程中發(fā)生拉伸、變形，形成“耗散”和“恢復(fù)效應(yīng)”。排間擾動還是發(fā)動機(jī)噪聲的最主要來源之一，研究表明葉片排本身的湍流現(xiàn)象(上游尾跡、角渦和葉尖泄漏渦等)作用于相鄰葉排時，會誘發(fā)湍流干涉噪聲[6-8]。此外，排間擾動還會影響葉輪機(jī)葉片的氣彈特性，誘發(fā)顫振，引起受迫響應(yīng)問題[9-10]。因此，葉輪機(jī)的排間干擾效應(yīng)十分重要，要進(jìn)行相關(guān)的數(shù)值研究，必須進(jìn)行非定常的數(shù)值仿真。

然而，對于葉輪機(jī)來說，常規(guī)的非定常計算通常采用時域的直接推進(jìn)方法，往往需要進(jìn)行多通道，甚至全通道的數(shù)值模擬[11-12]，計算量巨大，對計算資源的消耗往往比定常模擬高兩個數(shù)量級以上。因此，時域方法雖然比較精確，可以捕捉擾動與主流之間的非線性作用。但時域方法的耗時問題，限制了其在工程中的應(yīng)用。

另一類可用于葉輪機(jī)非定常求解的方法是頻域方法，頻域方法通常是指將非定常變量表達(dá)為周期性函數(shù)的形式，再通過推導(dǎo)和變換最終將方程組表達(dá)為頻域形式的方程組，從而可以對非定常流場進(jìn)行求解，其最大的優(yōu)點(diǎn)在于求解速度快。目前，頻域方法中，使用較為廣泛的有He等[13]提出的非線性諧波方法(nonlinear harmonic method，NHM)，Hall等[14]提出的諧波平衡方法(harmonic balance method，HB)以及Gopinath等[15]提出的時間譜方法(time spectral method，TSM)。鑒于頻域方法的求解效率較高，其已經(jīng)被一些研究人員用來考察葉輪機(jī)中的非定常現(xiàn)象[16-17]。但是頻域方法的精確性不如時域方法，精確度與采用諧波的數(shù)量關(guān)系較大。同時，采用過多的諧波數(shù)會帶來計算資源的浪費(fèi)。因此，諧波數(shù)量的合理選取是一個值得研究的問題，目前，這方面的研究較少。

將諧波平衡方法引入到葉輪機(jī)的非定常數(shù)值模擬中，并以一對轉(zhuǎn)風(fēng)扇為研究對象，通過對比常規(guī)的時域非定常結(jié)果，考察了諧波數(shù)量對諧波平衡法精度和可靠性的影響。

1 物理模型和計算方法

葉輪機(jī)中的非定常數(shù)值仿真方法可分為時域非定常和頻域非定常兩種類型。時域非定常算法計算精度相對較高，但計算量巨大，往往難以被工程應(yīng)用所接受。而頻域非定常算法雖然會損失部分精度，但其可以大大減少計算資源的消耗。頻域非定常計算方法主要有時間線化方法、非線性諧波方法和諧波平衡法。其中，諧波平衡法能夠捕捉主流和擾動之間的強(qiáng)烈非線性擾動，且前提假設(shè)較少，方程簡單易于編程，進(jìn)而得到廣泛應(yīng)用。

1.1 諧波平衡法

諧波平衡方法是一種用于求解周期性非定常流動的高效方法。該方法的核心思想是將變量分解成總體平均項和時間擾動項，然后再將擾動項以傅里葉級數(shù)進(jìn)行近似。非定常流動變量可表示為

(1)

式(1)中：W為非定常守恒變量;W0為總體平均量;An和Bn為第n個諧波的傅里葉系數(shù);ωn是第n個諧波的角頻率。令W1，W2,…，W2N+1為2N+1個時刻的變量值，則有式(2)成立:

(2)

式(2)中:

(3)

W*=[W1W2…W2N+1]T

(4)

(5)

對于變量W*，其滿足N-S方程，寫成半離散形式為

(6)

式(6)中：V為網(wǎng)格單元體積；把式(2)代入式(6)有：

(7)

(8)

式(8)即為需要進(jìn)行求解的諧波平衡控制方程，τ為虛擬時間。容易知道，相對于定常控制方程，其多了一個由諧波產(chǎn)生的源項VDW*，且其控制方程數(shù)目是定常的2N+1倍。因此，理論上講，其計算量略高于定常的2N+1倍，相比時域非定常計算，其效率明顯高得多。

1.2 雙時間步方法

時域非定常計算，采用的是雙時間步法，對于半離散形式的方程，在時間上應(yīng)用二階的向后差分有

(9)

為對式(9)進(jìn)行求解，引入了虛擬時間步，將物理時間層的替換為虛擬時間層的,并添加虛擬時間推進(jìn)導(dǎo)數(shù)項，得：

(10)

式(10)即為用于求解時域非定常的控制方程，虛擬時間步的迭代稱為內(nèi)迭代，當(dāng)計算收斂時，Wn+1=Wm。

1.3 計算方法

在空間離散上，采用了中心差分的JST格式，為二階精度，其計算效率較高，穩(wěn)定性好，在工程中已經(jīng)被廣泛應(yīng)用和驗證。時間推進(jìn)上采用了隱式方法，采用了一階L-F分裂格式進(jìn)行了對流通量的線化處理，黏性項線化進(jìn)行了簡化，保留了對角項。此外，對于諧波控制方程，還需對諧波源項進(jìn)行線化處理，這里參考Sicot的處理方式[18]，完成這些線化以后就得到了求解的線性代數(shù)方程組，對整個線性代數(shù)方程組的求解采用了上下對稱的逐次超松弛迭代方法(LUSSOR)。對于隱式離散的方程組，系數(shù)矩陣進(jìn)行分裂后為

(L+D+U)ΔW*=-R*-VDW*-VDΔW*

(11)

式(11)中：L為下三角矩陣;U為上三角矩陣;D為對角陣。LUSSOR迭代進(jìn)行求解時的每一步包含了Smax次的向前掃掠和向后掃掠。即:

(12)

式(12)中:s為最內(nèi)層迭代;l為次內(nèi)層迭代;m為外部迭代。其具體的算法過程如圖1所示。

圖1 諧波Lussor算法流程

對于雙時間步法的非定常計算，方程的空間離散、數(shù)值方法跟諧波平衡法保持一致。而對應(yīng)的離散線性方程組，也采用LUSSOR方法進(jìn)行迭代。不同的是，雙時間步方法中由于沒有諧波源項，不需要對諧波源項進(jìn)行特殊處理，因此，其算法流程中沒有圖1中的次內(nèi)迭代部分(深色背景部分)，即l層的迭代。

2 研究對象及計算設(shè)置

研究對象為一對轉(zhuǎn)風(fēng)扇，該風(fēng)扇的設(shè)計流量為230 kg/s，風(fēng)扇進(jìn)口處的外徑為0.635 m，輪轂比為0.247。前轉(zhuǎn)子轉(zhuǎn)速為5 200 r/min，葉尖切線速度是345.8 m/s，后轉(zhuǎn)子轉(zhuǎn)速為-4 100 r/min，葉尖切線速度為-272.6 m/s。前轉(zhuǎn)子的設(shè)計壓比為1.53，后轉(zhuǎn)子的設(shè)計壓比為1.45，總設(shè)計壓比為2.218 5。如圖2所示為對轉(zhuǎn)風(fēng)扇的三維示意圖。

圖2 對轉(zhuǎn)風(fēng)扇三維圖

時域計算采用了雙時間步法，物理時間步取為后轉(zhuǎn)子所受前轉(zhuǎn)子擾動周期的1/20，其也是前轉(zhuǎn)子所受擾動周期的1/30。而對于諧波平衡方法，按照采用諧波數(shù)目的不同共進(jìn)行了7次數(shù)值模擬，其上下游擾動諧波數(shù)均保持相同，諧波數(shù)從1取到7，這樣以分析隨諧波數(shù)變化時，諧波方法精度的變化。

由于在周期性邊界上使用了相位滯后條件，因此只需要使用一個葉片通道網(wǎng)格。而在固壁上應(yīng)用了壁面函數(shù)來處理邊界層的流動，其對壁面第一層網(wǎng)格大小的依賴性較低，則H形網(wǎng)格即可以滿足數(shù)值計算的需要，對轉(zhuǎn)風(fēng)扇的三維網(wǎng)格如圖3所示。

圖3 對轉(zhuǎn)風(fēng)扇三維網(wǎng)格

網(wǎng)格在周向和徑向的網(wǎng)格數(shù)分別為57、65，兩排網(wǎng)格數(shù)總共約80×104。葉尖間隙部分采用了削尖處理，間隙在徑向上的網(wǎng)格數(shù)為4。

諧波平衡法的計算中采用了并行策略，鑒于諧波方法的良好并行特性，在計算資源足夠的情況下，其對計算時間的消耗隨諧波數(shù)的增加幾乎不變，完成對轉(zhuǎn)風(fēng)扇的非定常數(shù)值仿真耗時約為2.5 h。而同樣采用并行計算的時域非定常數(shù)值模擬則消耗了約28 h。即時域非定常的耗時約為諧波平衡方法的11倍，這證明了諧波平衡方法的高效性。

3 計算結(jié)果和分析

排間干擾中非常重要的一種擾動是上游尾跡的擾動，熵等值線圖可以用來顯示尾跡的傳播。圖4和圖5給出了20%和80%葉高的熵等值線對比。從圖4、圖5中可以看出，隨著諧波數(shù)的增加，上下游轉(zhuǎn)子通道交界面處熵等值線的連續(xù)性變得越來越好，對于20%截面來說，當(dāng)諧波數(shù)為1時，尾跡在交界面處的連續(xù)性最差，后轉(zhuǎn)子計算域中交界面附近處的尾跡形狀幾乎不可辨識，隨著諧波數(shù)提高到2，交界面處的尾跡逐漸清晰可辨，但連續(xù)性還是稍差。諧波數(shù)到3時，其連續(xù)性也已經(jīng)變得很好，之后，再隨著諧波數(shù)增加，其連續(xù)性會繼續(xù)變好，但整體的變化已經(jīng)不大。而對于80%葉高處的熵圖來說，當(dāng)諧波數(shù)為1時，雖然交界面附近傳播到下游的尾跡比較清楚，但是尾跡的寬度和強(qiáng)度都跟上游計算域有明顯差別，跟20%截面的情況類似，當(dāng)諧波數(shù)到3時，其連續(xù)性也得到了很大的改善，之后隨著諧波數(shù)的增加，其連續(xù)性的提升也并不大。

此外，通過對比還可以發(fā)現(xiàn)，諧波計算和時域計算的尾跡在下游葉排通道中的形態(tài)基本保持一致。對于諧波數(shù)為1的情況，由于諧波數(shù)較小，對尾跡的形狀刻畫得不夠細(xì)致，因此使得尾跡寬度模擬得較大，但對于諧波數(shù)3甚至更多的情況，其吻合得很好，這就證明了諧波方法的有效性，同時，這也跟Hall等[14]及Sicot等[19]的驗證分析相吻合。

為了進(jìn)一步分析諧波方法的準(zhǔn)確性，圖6、圖7給出了前、后轉(zhuǎn)子葉尖間隙和50%葉高尾跡附近處的靜壓變化對比。影響靜壓變化的擾動因素很多，勢擾動、尾跡擾動都很大程度影響著其變化規(guī)律。

URANS表示時域非定常，HB N=1，2,…，7表示諧波數(shù)

圖5 80%葉高瞬時熵等值線圖

對于前轉(zhuǎn)子來說，其僅受勢擾動的影響，因此其變化規(guī)律相對來說比較單一。這也能從圖6中反映出來。無論對于葉尖間隙處還是尾跡附近處的靜壓變化，其都是一個類正弦曲線。對于間隙內(nèi)的監(jiān)測點(diǎn)來說，其變化曲線的形狀隨著諧波數(shù)的提升并沒有發(fā)生大的變化，只是其平均值逐漸上升，變得越來越接近時域的變化曲線，分析認(rèn)為其主要原因在于隨著諧波數(shù)的提高，第二排計算的精確程度越來越高，由此影響了第一排靜壓的平均值，但其隨時間的變化規(guī)律基本沒有發(fā)生改變。而對于前轉(zhuǎn)子葉中尾跡附近的監(jiān)測點(diǎn)來說，隨著諧波數(shù)提高，其變化的規(guī)律有一些細(xì)微的改變，而平均值則沒有太大的波動，而當(dāng)諧波數(shù)達(dá)到3以上時，靜壓變化規(guī)律基本不再隨諧波數(shù)的變化而改變?？傮w來說，前轉(zhuǎn)子所受排間擾動的影響較為簡單，用較少的諧波數(shù)就可以達(dá)到比較滿意的結(jié)果，尤其對于葉中區(qū)域來說，2個諧波數(shù)已經(jīng)基本可以滿足工程應(yīng)用要求。

后轉(zhuǎn)子所受的擾動相對來說較為復(fù)雜，勢擾動、熵擾動和渦擾動同時存在于后轉(zhuǎn)子通道內(nèi)，對比圖6、圖7也能發(fā)現(xiàn)其變化規(guī)律更為復(fù)雜。其次，后轉(zhuǎn)子的靜壓變化幅值也比前轉(zhuǎn)子大很多，尤其對于間隙部分的流動來說，幾乎大一個數(shù)量級。

對于間隙內(nèi)的靜壓變化來說，諧波數(shù)為1時的結(jié)果跟時域的結(jié)果依然相差較大。而當(dāng)諧波數(shù)達(dá)到2時，計算結(jié)果準(zhǔn)確度大大提升，而當(dāng)諧波數(shù)達(dá)到4以后整個靜壓變化曲線形態(tài)已跟時域非常接近。總體來說，除了諧波數(shù)為1的情況，其靜壓變化相對于時域的誤差并不大，而最大誤差出現(xiàn)在0.4周期處，該處的變化是一個周期內(nèi)變化最復(fù)雜的地方。這表明對于比較復(fù)雜的變化，諧波平衡方法更容易產(chǎn)生誤差，但其變化的趨勢還是跟時域結(jié)果保持了一致。

圖6 前轉(zhuǎn)子監(jiān)測點(diǎn)的靜壓變化

圖7 后轉(zhuǎn)子監(jiān)測點(diǎn)的靜壓變化

對于后轉(zhuǎn)子尾跡附近的監(jiān)測點(diǎn)，從圖6(b)可以看出，其變化規(guī)律在所有監(jiān)測點(diǎn)中是最為復(fù)雜的，它的變化曲線跟正弦形態(tài)差別很大。因此，諧波的計算誤差會較大，對比時域的變化就可以發(fā)現(xiàn)這種情況，即使對于諧波數(shù)為7的情況，其仍有一定的誤差，但這種誤差相對于平均靜壓值來說并不算大(不到5%)；另一方面，諧波模擬的靜壓變化趨勢基本能夠與時域保持一致，且當(dāng)諧波數(shù)達(dá)到2或3時，其變化規(guī)律就跟時域的靜壓變化規(guī)律達(dá)到了較好的相似。Gopinath等也發(fā)現(xiàn)了類似的現(xiàn)象[20]，并認(rèn)為用諧波平衡方法來研究排間的干涉效應(yīng)是可行的。

4 結(jié)論

將諧波平衡方法引入到了葉輪機(jī)械的非定常數(shù)值模擬中，通過諧波源項的隱式線化處理和對應(yīng)隱式迭代方法的改進(jìn)，實(shí)現(xiàn)了葉輪機(jī)的高效非定常數(shù)值模擬。

然后將諧波平衡方法應(yīng)用到了某對轉(zhuǎn)風(fēng)扇的非定常數(shù)值模擬中，結(jié)果發(fā)現(xiàn)諧波平衡方法對計算時間的消耗僅為時域非定常計算的9%。進(jìn)一步對比時域非定常的計算結(jié)果，分析了諧波平衡方法的有效性和可靠性。結(jié)果發(fā)現(xiàn)，諧波平衡方法能夠較好捕捉排間的非定常擾動信息，能夠清晰刻畫勢擾動和尾跡的傳播過程。當(dāng)諧波數(shù)達(dá)到3或者以上時，其計算精度就基本可以滿足工程應(yīng)用的需要。