王子維, 劉東健, 陳 逖

(中國空氣動力研究與發(fā)展中心, 四川綿陽 621000)

引 言

壓氣機(jī)的氣動穩(wěn)定性是影響航空發(fā)動機(jī)性能、 穩(wěn)定性、 安全性和壽命的關(guān)鍵因素之一. 為了避免壓氣機(jī)出現(xiàn)氣動失穩(wěn), 須對其設(shè)計進(jìn)行折衷, 這通常會導(dǎo)致壓氣機(jī)設(shè)計點性能的降低. 通過對壓氣機(jī)的氣動失穩(wěn)現(xiàn)象進(jìn)行研究, 探索其流動機(jī)理, 對于找到可行的抑制失穩(wěn)的流動控制方式, 從而拓寬壓氣機(jī)的失穩(wěn)邊界, 改進(jìn)設(shè)計, 提高壓氣機(jī)設(shè)計點的性能, 有著積極的促進(jìn)作用.

旋轉(zhuǎn)失速是壓氣機(jī)氣動失穩(wěn)的主要表現(xiàn)形式之一. 目前, 旋轉(zhuǎn)失速的流動機(jī)理尚未完全建立起來[1]. 當(dāng)前對于旋轉(zhuǎn)失速的模擬主要分為建模分析和全三維數(shù)值模擬兩種. 通過對壓氣機(jī)的葉排效率和落后角等建模, 可以在一維尺度上完成壓氣機(jī)失速邊界的快速預(yù)測[2]. 然而, 這種方法對于探索壓氣機(jī)失速流動機(jī)理并無作用. 另外, 南航屠寶鋒[3], 郭晉等[4]采用三維非定?？蓧嚎sEuler方程和三維激波盤模型, 發(fā)展了一種可分析風(fēng)扇/壓氣機(jī)三維非定常動態(tài)失速過程的計算模型和計算程序. 這種計算模型可以較準(zhǔn)確地模擬失速前壓氣機(jī)內(nèi)部非定常流動的變化情況, 然而并沒有將黏性的影響考慮在內(nèi). 有較多的研究者對壓氣機(jī)進(jìn)行全三維黏性定常模擬, 通過計算發(fā)散來判斷壓氣機(jī)是否失穩(wěn)[5]. 然而, 這種方法有一定的局限性, 數(shù)值發(fā)散與流動失穩(wěn)并不能完全對應(yīng); 另外, 這種方法也無法模擬壓氣機(jī)失速前后非定常流動的變化過程. 由于壓氣機(jī)旋轉(zhuǎn)失速是全周的、 強(qiáng)非定常性的流動現(xiàn)象, 因此通過全環(huán)非定常模擬方法模擬旋轉(zhuǎn)失速的發(fā)展過程, 是相對更為準(zhǔn)確的方法. 例如， 鞠鵬飛等[6]采用自主研發(fā)的MAP程序?qū)无D(zhuǎn)子NASA Rotor37的失速發(fā)展過程進(jìn)行了模擬, 賀象等[7]采用商業(yè)軟件CFX對北航低速大尺寸單轉(zhuǎn)子的失速發(fā)展過程進(jìn)行了模擬, 李清鵬等[8]采用商業(yè)軟件NUMECA對某亞聲速軸流壓氣機(jī)轉(zhuǎn)子的失速過程進(jìn)行了模擬. 失速過程模擬的難度較大, 極易出現(xiàn)流量持續(xù)下降的計算發(fā)散現(xiàn)象, 由于多排壓氣機(jī)存在葉排間的轉(zhuǎn)靜干涉, 其失速過程的模擬難度更大. 吳艷輝等[9]采用商業(yè)軟件NUMECA對某單級軸流壓氣機(jī)的失速過程進(jìn)行了模擬. 然而, 其僅模擬了流動失穩(wěn)的起始過程, 并沒有模擬出失速由出現(xiàn)先兆到完全發(fā)展的整個過程.

因此, 為了準(zhǔn)確復(fù)現(xiàn)多排壓氣機(jī)中失速發(fā)展的完整過程, 發(fā)展了守恒的動態(tài)重疊網(wǎng)格技術(shù), 以準(zhǔn)確傳遞葉排間的流動信息; 發(fā)展了流量出口邊界條件, 以快速逼近失速邊界; 發(fā)展了節(jié)流閥邊界條件, 以準(zhǔn)確模擬失速發(fā)展的過程.

1 多排壓氣機(jī)旋轉(zhuǎn)失速模擬方法

1.1 守恒的動態(tài)重疊網(wǎng)格技術(shù)

發(fā)展了基于單元相交的守恒重疊網(wǎng)格方法. 以葉排交界面為對稱面, 由于數(shù)值方法是2階的, 因此在葉排交界面處將網(wǎng)格外推兩層虛擬網(wǎng)格. 以相鄰葉排網(wǎng)格單元(即貢獻(xiàn)單元)對虛擬網(wǎng)格單元進(jìn)行切割, 以虛擬網(wǎng)格單元落在各個貢獻(xiàn)單元內(nèi)的體積占虛擬網(wǎng)格單元體積的比例為權(quán)重系數(shù), 對各個守恒單元的守恒變量進(jìn)行加權(quán)平均, 得到虛擬網(wǎng)格單元內(nèi)守恒變量的值.

1.1.1 網(wǎng)格相交法

守恒性重疊網(wǎng)格技術(shù)的核心為網(wǎng)格相交法[10-12]. 對于相交的兩個網(wǎng)格單元Cell 1與Cell 2, 須求出這兩個單元的相交體. 由于結(jié)構(gòu)網(wǎng)格單元與非結(jié)構(gòu)網(wǎng)格單元均可分解為1個或多個四面體單元, 則這兩個單元可以表示為兩組四面體單元{T11,T12,…,T1n}與{T21,T22,…,T2m}. 于是上述問題便轉(zhuǎn)化為求這兩組四面體單元的相交體.

如圖1所示, 四面體P1屬于第1個網(wǎng)格單元, 四面體P2屬于第2個網(wǎng)格單元. 通過Sutherland-Hodgman 算法便可求出兩個四面體單元相交得到的多面體P3.

圖1 四面體相交示意圖Fig. 1 Tetrahedron intersection

將兩組四面體單元兩兩求交, 得到多個相交體{I1,I2…,Imn}, 然后取并集, 便可以得到兩個網(wǎng)格單元的相交體I. 已知{I1,I2…,Imn}中每個體積和體心坐標(biāo), 可以求出I的體積和體心坐標(biāo): {I1,I2…,Imn}的體積之和為I的體積; 以{I1,I2…,Imn}中每個單元的體積除以I的體積作為權(quán)重系數(shù), 對這些單元的體心坐標(biāo)進(jìn)行加權(quán)平均, 可以求出I的體心坐標(biāo).

通過上述步驟可以求出Cell 1與Cell 2的相交體I的體積和體心坐標(biāo). 不失一般性, 假設(shè)Cell 1為貢獻(xiàn)單元, Cell 2為插值單元. Cell 1的格心值為q1, 則加權(quán)后Cell 1的格心值對Cell 2的格心值貢獻(xiàn)為q1VI/VCell 2. 其中,VI為I的體積,VCell 2為Cell 2的體積.

1.1.2 快速的并行動態(tài)重疊網(wǎng)格技術(shù)

與主流做法相同, 基于MPI實現(xiàn)并行的動態(tài)重疊網(wǎng)格計算. 考慮到壓氣機(jī)幾何與運動規(guī)律的特殊性, 故提出兩個技術(shù), 用于提高并行動態(tài)重疊計算的效率.

第1個技術(shù)細(xì)節(jié)如下： 收集緊鄰?fù)粋€葉排交界面的網(wǎng)格塊, 將這些網(wǎng)格塊所屬的進(jìn)程收集起來, 單獨建立通信子； 同時在分區(qū)時, 盡量將不同交界面處的網(wǎng)格塊放在不一樣的進(jìn)程里. 通過這種方式, 在大規(guī)模模擬中, 基本可以做到每一個交界面處的網(wǎng)格塊所屬的進(jìn)程組與其他交界面處網(wǎng)格塊所屬的進(jìn)程組交集為空， 大大簡化了編程工作. 基于分治的思想， 將大的進(jìn)程分組拆分成一個個小的進(jìn)程組， 每個小的進(jìn)程分組處理其對應(yīng)的交界面處的動態(tài)重疊計算， 相比在大的進(jìn)程分組里直接處理所有交界處的動態(tài)重疊計算， 極大地降低了進(jìn)程空等情況出現(xiàn)的頻率, 提高了并行計算效率.

第2個技術(shù)的細(xì)節(jié)如下： 由于壓氣機(jī)固有的時間上的周期性特征, 其插值關(guān)系是周期性出現(xiàn)的, 存儲一個周期的插值關(guān)系便可以避免重復(fù)地進(jìn)行找點和插值關(guān)系計算. 考慮到壓氣機(jī)葉排交界面運動方式為確定的旋轉(zhuǎn)運動, 基于這個特征可以對一個周期內(nèi)插值關(guān)系的建立進(jìn)行進(jìn)一步的優(yōu)化, 在前一步插值關(guān)系的基礎(chǔ)上進(jìn)行較少的計算得到下一步的插值關(guān)系.

1.2 失速邊界逼近技術(shù)

在大流量工況下, 設(shè)置壓力出口邊界條件, 即設(shè)定出口機(jī)匣壁面處靜壓, 并采用徑向平衡方程得到出口靜壓分布. 然后通過逐步提高出口反壓的方式逐步推進(jìn)到壓氣機(jī)近失速點. 以近失速點的流場為初場, 進(jìn)行壓氣機(jī)穩(wěn)定邊界附近的流動模擬. 在壓氣機(jī)穩(wěn)定邊界附近, 考慮到出口反壓較高時, 若依然固定機(jī)匣處為靜壓, 有可能得到不符合物理實際的流場解. 為避免出現(xiàn)這種情況, 設(shè)計流量出口邊界. 通過二分法的形式逐步調(diào)低流量, 以最大限度逼近失速邊界.

以NASA Stage 35為研究對象[13], 其數(shù)模如圖2所示, 其設(shè)計參數(shù)如表1所示.

生成全環(huán)結(jié)構(gòu)化網(wǎng)格, 所有物面第1層網(wǎng)格的厚度保證y+=1, 網(wǎng)格塊數(shù)量為728, 網(wǎng)格單元數(shù)為2.2×107, 轉(zhuǎn)子網(wǎng)格與靜子網(wǎng)格之間通過滑移網(wǎng)格面交換數(shù)據(jù), NASA Stage 35的物面網(wǎng)格見圖3.

圖2 NASA Stage 35數(shù)模Fig. 2 Numerical model of NASA Stage 35

表1 NASA Stage 35的設(shè)計參數(shù)Table 1 Design parameters of NASA Stage 35

圖3 物面網(wǎng)格Fig. 3 Grids of wall surface

基于自主研發(fā)的軸流壓氣機(jī)數(shù)值模擬軟件ASPAC對NASA Stage 35進(jìn)行了堵塞狀態(tài)到失速邊界狀態(tài)的模擬. ASPAC采用有限體積法求解Reynolds時均N-S方程, 對流通量采用Roe格式計算, 黏性通量采用中心差分方法計算, 湍流模型為SA模型. 軟件基于MPI實現(xiàn)了并行計算能力. 對于非定常模擬, 采用Jameson的雙時間步方法. ASPAC在葉輪機(jī)械的數(shù)值模擬上經(jīng)由多個標(biāo)模進(jìn)行了計算驗證, 計算結(jié)果與實驗結(jié)果吻合較好, 計算精度與商業(yè)軟件的計算精度相當(dāng)[14]. 在全環(huán)非定常模擬中, 在壓氣機(jī)轉(zhuǎn)一圈的時間內(nèi)采用2 000個物理時間步, 子迭代步數(shù)為50. 模擬結(jié)果與NASA數(shù)值模擬程序TURBO[15]和實驗結(jié)果的對比如圖4所示. 其中, 實線為3種結(jié)果對應(yīng)的壓比, 綠色虛線與藍(lán)色虛線分別表示ASPAC和TURBO計算結(jié)果與實驗值的誤差, 計算公式如下

其中, 某一流量下對應(yīng)的壓比通過線性插值得到.

圖4 模擬與實驗結(jié)果對比Fig. 4 Comparisons between simulation and experiment

ASPAC和TURBO的計算結(jié)果與實驗結(jié)果的誤差如表2所示, 基于實驗值計算了相對誤差. 在堵點附近, 兩個軟件的模擬結(jié)果與實驗結(jié)果的誤差均較大, 這是由于堵點附近流量-壓比曲線的斜率絕對值極大, 較少的擾動便會導(dǎo)致很大的誤差; 隨著流動向近失速點靠近, 兩個軟件的模擬結(jié)果與實驗結(jié)果的誤差很快縮小.

表2 數(shù)值模擬結(jié)果與實驗結(jié)果誤差Table 2 Differences between simulation and experiment

1.3 節(jié)流閥邊界條件

以近失速狀態(tài)的非定常流場為初場, 進(jìn)行全環(huán)非定常續(xù)算. 通過添加節(jié)流閥邊界條件, 模擬了失速由起始到完全發(fā)展的全過程. 節(jié)流閥邊界條件為

2 模擬結(jié)果

在失速模擬的計算初期, 固定機(jī)匣處靜壓以排除人為擾動的可能, 在流量穩(wěn)定以后, 根據(jù)穩(wěn)定狀態(tài)下的流量和出口壓力, 確定節(jié)流閥開度系數(shù), 并將出口邊界條件改為節(jié)流閥邊界條件. 在失速模擬的整個過程中, 節(jié)流閥開度系數(shù)保持不變.

失速發(fā)展整個過程的流量變化如圖5所示. 在添加節(jié)流閥邊界條件后, 流量在開始幾圈保持穩(wěn)定狀態(tài), 在第10圈的時候, 隨著失速擾動的出現(xiàn), 流量急速下降, 約0.5圈后降到了最低點. 流量在最低點維持到第14圈, 并花費了約3圈的時間進(jìn)入到穩(wěn)定失速狀態(tài). 流量劇烈變化的細(xì)節(jié)方法如圖6所示.

圖5 失速發(fā)展過程的入口流量變化Fig. 5 Variation of inlet mass flow during stall development

圖6 失速發(fā)展過程的入口流量變化(細(xì)節(jié)圖)Fig. 6 Variation of inlet mass flow during stall development(details)

壓氣機(jī)由穩(wěn)定運行到失速完全發(fā)展過程的流量-出口靜壓特性如圖7所示, 流量-總壓比特性如圖8所示. 由圖7左邊一段的曲線可以明顯地看到出口靜壓與流量二次方之間的正比關(guān)系, 這與節(jié)流閥的特性相匹配.

圖7 失速發(fā)展過程的流量-出口靜壓特性圖Fig. 7 Variation of mass flow-outlet static pressure during stall development

圖8 失速發(fā)展過程的流量-總壓比特性圖Fig. 8 Variation of mass flow-total pressure ratio during stall development

然而, 由圖8可以發(fā)現(xiàn), 在失速發(fā)展過程中, 流量與總壓比的關(guān)系要復(fù)雜得多. 由于入口總壓指定為一個大氣壓, 在模擬過程中保持不變, 則總壓比的變化實際上正比于出口總壓的變化. 圖8箭頭所指的近失速狀態(tài)處, 在數(shù)值上減小出口節(jié)流閥開度, 導(dǎo)致壓比快速升高, 流量快速下降, 這個擾動最終導(dǎo)致了旋轉(zhuǎn)失速的產(chǎn)生. 可以看到, 在失速擾動出現(xiàn)以后, 隨著流量的下降, 出口總壓也隨之下降; 在流量由最低點逐步增大并逼近失速完全狀態(tài)的過程中, 出口總壓隨之上升, 但始終低于流量下行過程中對應(yīng)的出口總壓. 最終流動穩(wěn)定在圖8中箭頭所指的旋轉(zhuǎn)失速完全發(fā)展所在的位置, 旋轉(zhuǎn)失速完全發(fā)展后的總壓比和流量均小于近失速點的總壓比和流量.

面向壓氣機(jī)入口, 沿著逆時針方向, 在每個轉(zhuǎn)子葉片通道的機(jī)匣處依次設(shè)置數(shù)值探針, 每個通道的探針與葉片的相對位置均相同. 圖9展示了失速發(fā)展過程中機(jī)匣處探針上的靜壓隨時間的變化. 圖9對靜壓進(jìn)行了歸一化處理, 展示了1, 4, 7, 10, …, 34號轉(zhuǎn)子葉片通道的靜壓隨時間的變化. 可以看到, 失速擾動起始伴隨著機(jī)匣靜壓的躍升, 且不同葉片通道的機(jī)匣處靜壓躍升幾乎是同步的. 經(jīng)過大約7圈的時間, 擾動逐步增大, 并在最后發(fā)展成穩(wěn)定的失速團(tuán). 失速團(tuán)轉(zhuǎn)速與壓氣機(jī)轉(zhuǎn)速方向相同, 約為壓氣機(jī)轉(zhuǎn)速的46.41%.

圖9 不同葉片通道機(jī)匣處的靜壓變化Fig. 9 Variation of static pressure on the shroud in different blade passages

失速先兆階段與完全失速時, 壓氣機(jī)s1流面(95%葉高)的相對Mach數(shù)分布及流線如圖10所示. 如圖10(a)所示, 由于節(jié)流閥開度減小, 流動出現(xiàn)了堵塞, 堵塞的現(xiàn)象發(fā)生在所有葉片通道. 如圖10(b)與(c)所示, 在失速完全發(fā)展階段, 出現(xiàn)了一大一小兩個失速團(tuán). 沿著轉(zhuǎn)子轉(zhuǎn)動方向, 在失速團(tuán)正方向的葉片, 來流沖角減小, 將會解除該葉片周圍的流動分離; 在失速團(tuán)反方向的葉片, 來流沖角增大, 將會導(dǎo)致該葉片出現(xiàn)流動分離; 這兩個效應(yīng)結(jié)合, 導(dǎo)致了失速團(tuán)沿著葉片轉(zhuǎn)動反方向傳播.

(a) Stall inception

由于硬盤存儲空間的限制, 沒有存儲從流動堵塞階段到失速團(tuán)完全發(fā)展階段每一步的流場, 無法詳細(xì)分析從流動堵塞階段到失速完全發(fā)展時壓氣機(jī)內(nèi)部流動的演化細(xì)節(jié). 然而, 從失速發(fā)展過程的流量變化情況, 以及機(jī)匣處采樣點的靜壓變化情況, 可以分析得知, 在通過關(guān)緊節(jié)流閥施加擾動之后, 流量減小, 機(jī)匣處靜壓升高, 隨著流動的部分恢復(fù), 最終形成了兩個典型的失速團(tuán).

壓氣機(jī)由穩(wěn)定運行到失速完全發(fā)展前, 失速發(fā)展中和失速完全發(fā)展時, 軸向Mach數(shù)等值面分別如圖11～13所示.

圖11 轉(zhuǎn)動到第4.5圈時的軸向Mach數(shù)等值面(Ma=-0.1, 靜壓著色)Fig. 11 Iso-surface of axial Mach number at 4.5th revolution(Ma=-0.1, colored by static pressure)

圖12 轉(zhuǎn)動到第12.65圈時的軸向Mach數(shù)等值面(Ma=-0.1, 靜壓著色)Fig. 12 Iso-surface of axial Mach number at 12.65th revolution(Ma=-0.1, colored by static pressure)

圖13 轉(zhuǎn)動到第26.7圈時的軸向Mach數(shù)等值面(Ma=-0.1, 靜壓著色)Fig. 13 Iso-surface of axial Mach number at 26.7th revolution(Ma=-0.1, colored by static pressure)

由于展示的是Ma=-0.1的等值面, 可以認(rèn)為該等值面在一定程度上代表了逆流區(qū)的位置和大小. 可以看到, 逆流區(qū)逐步由轉(zhuǎn)子的吸力面蔓延到壓力面; 同時, 逆流區(qū)逐漸由軸對稱變?yōu)榉禽S對稱. 由圖13可以看到, 當(dāng)旋轉(zhuǎn)失速完全發(fā)展后, 轉(zhuǎn)子區(qū)域出現(xiàn)一大一小兩個失速團(tuán); 同時, 在靜子的壓力面也出現(xiàn)了逆流區(qū), 且與轉(zhuǎn)子區(qū)域的失速團(tuán)相鄰. 最后, 在旋轉(zhuǎn)失速發(fā)展的整個過程中, 逆流區(qū)主要存在于壓氣機(jī)的葉尖區(qū)域.

壓氣機(jī)失速完全發(fā)展后, 95%葉高的熵分布如圖14所示, 高熵區(qū)域為失速團(tuán), 可以清晰地看到一大一小兩個失速團(tuán)， 這兩個失速團(tuán)的結(jié)構(gòu)幾乎不隨時間變化.

圖14 轉(zhuǎn)動到第26.7圈時95%葉高的熵分布Fig. 14 Distribution of entropy at 95% span at 26.7th revolution

3 結(jié)論

ASPAC和TURBO的模擬結(jié)果與實驗結(jié)果吻合較好. 尖峰型擾動幾乎同時發(fā)生在不同的轉(zhuǎn)子葉片通道上. 然后, 擾動逐漸增大, 最終形成的旋轉(zhuǎn)失速團(tuán)轉(zhuǎn)速為壓氣機(jī)轉(zhuǎn)速的46.41%, 同時失速區(qū)的流動結(jié)構(gòu)幾乎不隨時間變化. 在尖峰型擾動后, 逆流從轉(zhuǎn)子吸力面蔓延到轉(zhuǎn)子壓力面, 同時也影響到了靜子. 在整個過程中, 逆流主要存在于葉尖附近.