王子維, 劉東健, 陳 逖
(中國(guó)空氣動(dòng)力研究與發(fā)展中心, 四川綿陽(yáng) 621000)
壓氣機(jī)的氣動(dòng)穩(wěn)定性是影響航空發(fā)動(dòng)機(jī)性能、 穩(wěn)定性、 安全性和壽命的關(guān)鍵因素之一. 為了避免壓氣機(jī)出現(xiàn)氣動(dòng)失穩(wěn), 須對(duì)其設(shè)計(jì)進(jìn)行折衷, 這通常會(huì)導(dǎo)致壓氣機(jī)設(shè)計(jì)點(diǎn)性能的降低. 通過(guò)對(duì)壓氣機(jī)的氣動(dòng)失穩(wěn)現(xiàn)象進(jìn)行研究, 探索其流動(dòng)機(jī)理, 對(duì)于找到可行的抑制失穩(wěn)的流動(dòng)控制方式, 從而拓寬壓氣機(jī)的失穩(wěn)邊界, 改進(jìn)設(shè)計(jì), 提高壓氣機(jī)設(shè)計(jì)點(diǎn)的性能, 有著積極的促進(jìn)作用.
旋轉(zhuǎn)失速是壓氣機(jī)氣動(dòng)失穩(wěn)的主要表現(xiàn)形式之一. 目前, 旋轉(zhuǎn)失速的流動(dòng)機(jī)理尚未完全建立起來(lái)[1]. 當(dāng)前對(duì)于旋轉(zhuǎn)失速的模擬主要分為建模分析和全三維數(shù)值模擬兩種. 通過(guò)對(duì)壓氣機(jī)的葉排效率和落后角等建模, 可以在一維尺度上完成壓氣機(jī)失速邊界的快速預(yù)測(cè)[2]. 然而, 這種方法對(duì)于探索壓氣機(jī)失速流動(dòng)機(jī)理并無(wú)作用. 另外, 南航屠寶鋒[3], 郭晉等[4]采用三維非定常可壓縮Euler方程和三維激波盤模型, 發(fā)展了一種可分析風(fēng)扇/壓氣機(jī)三維非定常動(dòng)態(tài)失速過(guò)程的計(jì)算模型和計(jì)算程序. 這種計(jì)算模型可以較準(zhǔn)確地模擬失速前壓氣機(jī)內(nèi)部非定常流動(dòng)的變化情況, 然而并沒(méi)有將黏性的影響考慮在內(nèi). 有較多的研究者對(duì)壓氣機(jī)進(jìn)行全三維黏性定常模擬, 通過(guò)計(jì)算發(fā)散來(lái)判斷壓氣機(jī)是否失穩(wěn)[5]. 然而, 這種方法有一定的局限性, 數(shù)值發(fā)散與流動(dòng)失穩(wěn)并不能完全對(duì)應(yīng); 另外, 這種方法也無(wú)法模擬壓氣機(jī)失速前后非定常流動(dòng)的變化過(guò)程. 由于壓氣機(jī)旋轉(zhuǎn)失速是全周的、 強(qiáng)非定常性的流動(dòng)現(xiàn)象, 因此通過(guò)全環(huán)非定常模擬方法模擬旋轉(zhuǎn)失速的發(fā)展過(guò)程, 是相對(duì)更為準(zhǔn)確的方法. 例如, 鞠鵬飛等[6]采用自主研發(fā)的MAP程序?qū)无D(zhuǎn)子NASA Rotor37的失速發(fā)展過(guò)程進(jìn)行了模擬, 賀象等[7]采用商業(yè)軟件CFX對(duì)北航低速大尺寸單轉(zhuǎn)子的失速發(fā)展過(guò)程進(jìn)行了模擬, 李清鵬等[8]采用商業(yè)軟件NUMECA對(duì)某亞聲速軸流壓氣機(jī)轉(zhuǎn)子的失速過(guò)程進(jìn)行了模擬. 失速過(guò)程模擬的難度較大, 極易出現(xiàn)流量持續(xù)下降的計(jì)算發(fā)散現(xiàn)象, 由于多排壓氣機(jī)存在葉排間的轉(zhuǎn)靜干涉, 其失速過(guò)程的模擬難度更大. 吳艷輝等[9]采用商業(yè)軟件NUMECA對(duì)某單級(jí)軸流壓氣機(jī)的失速過(guò)程進(jìn)行了模擬. 然而, 其僅模擬了流動(dòng)失穩(wěn)的起始過(guò)程, 并沒(méi)有模擬出失速由出現(xiàn)先兆到完全發(fā)展的整個(gè)過(guò)程.
因此, 為了準(zhǔn)確復(fù)現(xiàn)多排壓氣機(jī)中失速發(fā)展的完整過(guò)程, 發(fā)展了守恒的動(dòng)態(tài)重疊網(wǎng)格技術(shù), 以準(zhǔn)確傳遞葉排間的流動(dòng)信息; 發(fā)展了流量出口邊界條件, 以快速逼近失速邊界; 發(fā)展了節(jié)流閥邊界條件, 以準(zhǔn)確模擬失速發(fā)展的過(guò)程.
發(fā)展了基于單元相交的守恒重疊網(wǎng)格方法. 以葉排交界面為對(duì)稱面, 由于數(shù)值方法是2階的, 因此在葉排交界面處將網(wǎng)格外推兩層虛擬網(wǎng)格. 以相鄰葉排網(wǎng)格單元(即貢獻(xiàn)單元)對(duì)虛擬網(wǎng)格單元進(jìn)行切割, 以虛擬網(wǎng)格單元落在各個(gè)貢獻(xiàn)單元內(nèi)的體積占虛擬網(wǎng)格單元體積的比例為權(quán)重系數(shù), 對(duì)各個(gè)守恒單元的守恒變量進(jìn)行加權(quán)平均, 得到虛擬網(wǎng)格單元內(nèi)守恒變量的值.
1.1.1 網(wǎng)格相交法
守恒性重疊網(wǎng)格技術(shù)的核心為網(wǎng)格相交法[10-12]. 對(duì)于相交的兩個(gè)網(wǎng)格單元Cell 1與Cell 2, 須求出這兩個(gè)單元的相交體. 由于結(jié)構(gòu)網(wǎng)格單元與非結(jié)構(gòu)網(wǎng)格單元均可分解為1個(gè)或多個(gè)四面體單元, 則這兩個(gè)單元可以表示為兩組四面體單元{T11,T12,…,T1n}與{T21,T22,…,T2m}. 于是上述問(wèn)題便轉(zhuǎn)化為求這兩組四面體單元的相交體.
如圖1所示, 四面體P1屬于第1個(gè)網(wǎng)格單元, 四面體P2屬于第2個(gè)網(wǎng)格單元. 通過(guò)Sutherland-Hodgman 算法便可求出兩個(gè)四面體單元相交得到的多面體P3.
圖1 四面體相交示意圖Fig. 1 Tetrahedron intersection
將兩組四面體單元兩兩求交, 得到多個(gè)相交體{I1,I2…,Imn}, 然后取并集, 便可以得到兩個(gè)網(wǎng)格單元的相交體I. 已知{I1,I2…,Imn}中每個(gè)體積和體心坐標(biāo), 可以求出I的體積和體心坐標(biāo): {I1,I2…,Imn}的體積之和為I的體積; 以{I1,I2…,Imn}中每個(gè)單元的體積除以I的體積作為權(quán)重系數(shù), 對(duì)這些單元的體心坐標(biāo)進(jìn)行加權(quán)平均, 可以求出I的體心坐標(biāo).
通過(guò)上述步驟可以求出Cell 1與Cell 2的相交體I的體積和體心坐標(biāo). 不失一般性, 假設(shè)Cell 1為貢獻(xiàn)單元, Cell 2為插值單元. Cell 1的格心值為q1, 則加權(quán)后Cell 1的格心值對(duì)Cell 2的格心值貢獻(xiàn)為q1VI/VCell 2. 其中,VI為I的體積,VCell 2為Cell 2的體積.
1.1.2 快速的并行動(dòng)態(tài)重疊網(wǎng)格技術(shù)
與主流做法相同, 基于MPI實(shí)現(xiàn)并行的動(dòng)態(tài)重疊網(wǎng)格計(jì)算. 考慮到壓氣機(jī)幾何與運(yùn)動(dòng)規(guī)律的特殊性, 故提出兩個(gè)技術(shù), 用于提高并行動(dòng)態(tài)重疊計(jì)算的效率.
第1個(gè)技術(shù)細(xì)節(jié)如下: 收集緊鄰?fù)粋€(gè)葉排交界面的網(wǎng)格塊, 將這些網(wǎng)格塊所屬的進(jìn)程收集起來(lái), 單獨(dú)建立通信子; 同時(shí)在分區(qū)時(shí), 盡量將不同交界面處的網(wǎng)格塊放在不一樣的進(jìn)程里. 通過(guò)這種方式, 在大規(guī)模模擬中, 基本可以做到每一個(gè)交界面處的網(wǎng)格塊所屬的進(jìn)程組與其他交界面處網(wǎng)格塊所屬的進(jìn)程組交集為空, 大大簡(jiǎn)化了編程工作. 基于分治的思想, 將大的進(jìn)程分組拆分成一個(gè)個(gè)小的進(jìn)程組, 每個(gè)小的進(jìn)程分組處理其對(duì)應(yīng)的交界面處的動(dòng)態(tài)重疊計(jì)算, 相比在大的進(jìn)程分組里直接處理所有交界處的動(dòng)態(tài)重疊計(jì)算, 極大地降低了進(jìn)程空等情況出現(xiàn)的頻率, 提高了并行計(jì)算效率.
第2個(gè)技術(shù)的細(xì)節(jié)如下: 由于壓氣機(jī)固有的時(shí)間上的周期性特征, 其插值關(guān)系是周期性出現(xiàn)的, 存儲(chǔ)一個(gè)周期的插值關(guān)系便可以避免重復(fù)地進(jìn)行找點(diǎn)和插值關(guān)系計(jì)算. 考慮到壓氣機(jī)葉排交界面運(yùn)動(dòng)方式為確定的旋轉(zhuǎn)運(yùn)動(dòng), 基于這個(gè)特征可以對(duì)一個(gè)周期內(nèi)插值關(guān)系的建立進(jìn)行進(jìn)一步的優(yōu)化, 在前一步插值關(guān)系的基礎(chǔ)上進(jìn)行較少的計(jì)算得到下一步的插值關(guān)系.
在大流量工況下, 設(shè)置壓力出口邊界條件, 即設(shè)定出口機(jī)匣壁面處?kù)o壓, 并采用徑向平衡方程得到出口靜壓分布. 然后通過(guò)逐步提高出口反壓的方式逐步推進(jìn)到壓氣機(jī)近失速點(diǎn). 以近失速點(diǎn)的流場(chǎng)為初場(chǎng), 進(jìn)行壓氣機(jī)穩(wěn)定邊界附近的流動(dòng)模擬. 在壓氣機(jī)穩(wěn)定邊界附近, 考慮到出口反壓較高時(shí), 若依然固定機(jī)匣處為靜壓, 有可能得到不符合物理實(shí)際的流場(chǎng)解. 為避免出現(xiàn)這種情況, 設(shè)計(jì)流量出口邊界. 通過(guò)二分法的形式逐步調(diào)低流量, 以最大限度逼近失速邊界.
以NASA Stage 35為研究對(duì)象[13], 其數(shù)模如圖2所示, 其設(shè)計(jì)參數(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)格之間通過(guò)滑移網(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è)計(jì)參數(shù)Table 1 Design parameters of NASA Stage 35
圖3 物面網(wǎng)格Fig. 3 Grids of wall surface
基于自主研發(fā)的軸流壓氣機(jī)數(shù)值模擬軟件ASPAC對(duì)NASA Stage 35進(jìn)行了堵塞狀態(tài)到失速邊界狀態(tài)的模擬. ASPAC采用有限體積法求解Reynolds時(shí)均N-S方程, 對(duì)流通量采用Roe格式計(jì)算, 黏性通量采用中心差分方法計(jì)算, 湍流模型為SA模型. 軟件基于MPI實(shí)現(xiàn)了并行計(jì)算能力. 對(duì)于非定常模擬, 采用Jameson的雙時(shí)間步方法. ASPAC在葉輪機(jī)械的數(shù)值模擬上經(jīng)由多個(gè)標(biāo)模進(jìn)行了計(jì)算驗(yàn)證, 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好, 計(jì)算精度與商業(yè)軟件的計(jì)算精度相當(dāng)[14]. 在全環(huán)非定常模擬中, 在壓氣機(jī)轉(zhuǎn)一圈的時(shí)間內(nèi)采用2 000個(gè)物理時(shí)間步, 子迭代步數(shù)為50. 模擬結(jié)果與NASA數(shù)值模擬程序TURBO[15]和實(shí)驗(yàn)結(jié)果的對(duì)比如圖4所示. 其中, 實(shí)線為3種結(jié)果對(duì)應(yīng)的壓比, 綠色虛線與藍(lán)色虛線分別表示ASPAC和TURBO計(jì)算結(jié)果與實(shí)驗(yàn)值的誤差, 計(jì)算公式如下
其中, 某一流量下對(duì)應(yīng)的壓比通過(guò)線性插值得到.
圖4 模擬與實(shí)驗(yàn)結(jié)果對(duì)比Fig. 4 Comparisons between simulation and experiment
ASPAC和TURBO的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的誤差如表2所示, 基于實(shí)驗(yàn)值計(jì)算了相對(duì)誤差. 在堵點(diǎn)附近, 兩個(gè)軟件的模擬結(jié)果與實(shí)驗(yàn)結(jié)果的誤差均較大, 這是由于堵點(diǎn)附近流量-壓比曲線的斜率絕對(duì)值極大, 較少的擾動(dòng)便會(huì)導(dǎo)致很大的誤差; 隨著流動(dòng)向近失速點(diǎn)靠近, 兩個(gè)軟件的模擬結(jié)果與實(shí)驗(yàn)結(jié)果的誤差很快縮小.
表2 數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果誤差Table 2 Differences between simulation and experiment
以近失速狀態(tài)的非定常流場(chǎng)為初場(chǎng), 進(jìn)行全環(huán)非定常續(xù)算. 通過(guò)添加節(jié)流閥邊界條件, 模擬了失速由起始到完全發(fā)展的全過(guò)程. 節(jié)流閥邊界條件為
在失速模擬的計(jì)算初期, 固定機(jī)匣處?kù)o壓以排除人為擾動(dòng)的可能, 在流量穩(wěn)定以后, 根據(jù)穩(wěn)定狀態(tài)下的流量和出口壓力, 確定節(jié)流閥開度系數(shù), 并將出口邊界條件改為節(jié)流閥邊界條件. 在失速模擬的整個(gè)過(guò)程中, 節(jié)流閥開度系數(shù)保持不變.
失速發(fā)展整個(gè)過(guò)程的流量變化如圖5所示. 在添加節(jié)流閥邊界條件后, 流量在開始幾圈保持穩(wěn)定狀態(tài), 在第10圈的時(shí)候, 隨著失速擾動(dòng)的出現(xiàn), 流量急速下降, 約0.5圈后降到了最低點(diǎn). 流量在最低點(diǎn)維持到第14圈, 并花費(fèi)了約3圈的時(shí)間進(jìn)入到穩(wěn)定失速狀態(tài). 流量劇烈變化的細(xì)節(jié)方法如圖6所示.
圖5 失速發(fā)展過(guò)程的入口流量變化Fig. 5 Variation of inlet mass flow during stall development
圖6 失速發(fā)展過(guò)程的入口流量變化(細(xì)節(jié)圖)Fig. 6 Variation of inlet mass flow during stall development(details)
壓氣機(jī)由穩(wěn)定運(yùn)行到失速完全發(fā)展過(guò)程的流量-出口靜壓特性如圖7所示, 流量-總壓比特性如圖8所示. 由圖7左邊一段的曲線可以明顯地看到出口靜壓與流量二次方之間的正比關(guān)系, 這與節(jié)流閥的特性相匹配.
圖7 失速發(fā)展過(guò)程的流量-出口靜壓特性圖Fig. 7 Variation of mass flow-outlet static pressure during stall development
圖8 失速發(fā)展過(guò)程的流量-總壓比特性圖Fig. 8 Variation of mass flow-total pressure ratio during stall development
然而, 由圖8可以發(fā)現(xiàn), 在失速發(fā)展過(guò)程中, 流量與總壓比的關(guān)系要復(fù)雜得多. 由于入口總壓指定為一個(gè)大氣壓, 在模擬過(guò)程中保持不變, 則總壓比的變化實(shí)際上正比于出口總壓的變化. 圖8箭頭所指的近失速狀態(tài)處, 在數(shù)值上減小出口節(jié)流閥開度, 導(dǎo)致壓比快速升高, 流量快速下降, 這個(gè)擾動(dòng)最終導(dǎo)致了旋轉(zhuǎn)失速的產(chǎn)生. 可以看到, 在失速擾動(dòng)出現(xiàn)以后, 隨著流量的下降, 出口總壓也隨之下降; 在流量由最低點(diǎn)逐步增大并逼近失速完全狀態(tài)的過(guò)程中, 出口總壓隨之上升, 但始終低于流量下行過(guò)程中對(duì)應(yīng)的出口總壓. 最終流動(dòng)穩(wěn)定在圖8中箭頭所指的旋轉(zhuǎn)失速完全發(fā)展所在的位置, 旋轉(zhuǎn)失速完全發(fā)展后的總壓比和流量均小于近失速點(diǎn)的總壓比和流量.
面向壓氣機(jī)入口, 沿著逆時(shí)針?lè)较? 在每個(gè)轉(zhuǎn)子葉片通道的機(jī)匣處依次設(shè)置數(shù)值探針, 每個(gè)通道的探針與葉片的相對(duì)位置均相同. 圖9展示了失速發(fā)展過(guò)程中機(jī)匣處探針上的靜壓隨時(shí)間的變化. 圖9對(duì)靜壓進(jìn)行了歸一化處理, 展示了1, 4, 7, 10, …, 34號(hào)轉(zhuǎn)子葉片通道的靜壓隨時(shí)間的變化. 可以看到, 失速擾動(dòng)起始伴隨著機(jī)匣靜壓的躍升, 且不同葉片通道的機(jī)匣處?kù)o壓躍升幾乎是同步的. 經(jīng)過(guò)大約7圈的時(shí)間, 擾動(dòng)逐步增大, 并在最后發(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
失速先兆階段與完全失速時(shí), 壓氣機(jī)s1流面(95%葉高)的相對(duì)Mach數(shù)分布及流線如圖10所示. 如圖10(a)所示, 由于節(jié)流閥開度減小, 流動(dòng)出現(xiàn)了堵塞, 堵塞的現(xiàn)象發(fā)生在所有葉片通道. 如圖10(b)與(c)所示, 在失速完全發(fā)展階段, 出現(xiàn)了一大一小兩個(gè)失速團(tuán). 沿著轉(zhuǎn)子轉(zhuǎn)動(dòng)方向, 在失速團(tuán)正方向的葉片, 來(lái)流沖角減小, 將會(huì)解除該葉片周圍的流動(dòng)分離; 在失速團(tuán)反方向的葉片, 來(lái)流沖角增大, 將會(huì)導(dǎo)致該葉片出現(xiàn)流動(dòng)分離; 這兩個(gè)效應(yīng)結(jié)合, 導(dǎo)致了失速團(tuán)沿著葉片轉(zhuǎn)動(dòng)反方向傳播.
(a) Stall inception
由于硬盤存儲(chǔ)空間的限制, 沒(méi)有存儲(chǔ)從流動(dòng)堵塞階段到失速團(tuán)完全發(fā)展階段每一步的流場(chǎng), 無(wú)法詳細(xì)分析從流動(dòng)堵塞階段到失速完全發(fā)展時(shí)壓氣機(jī)內(nèi)部流動(dòng)的演化細(xì)節(jié). 然而, 從失速發(fā)展過(guò)程的流量變化情況, 以及機(jī)匣處采樣點(diǎn)的靜壓變化情況, 可以分析得知, 在通過(guò)關(guān)緊節(jié)流閥施加擾動(dòng)之后, 流量減小, 機(jī)匣處?kù)o壓升高, 隨著流動(dòng)的部分恢復(fù), 最終形成了兩個(gè)典型的失速團(tuán).
壓氣機(jī)由穩(wěn)定運(yùn)行到失速完全發(fā)展前, 失速發(fā)展中和失速完全發(fā)展時(shí), 軸向Mach數(shù)等值面分別如圖11~13所示.
圖11 轉(zhuǎn)動(dòng)到第4.5圈時(shí)的軸向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)動(dòng)到第12.65圈時(shí)的軸向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)動(dòng)到第26.7圈時(shí)的軸向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)子的吸力面蔓延到壓力面; 同時(shí), 逆流區(qū)逐漸由軸對(duì)稱變?yōu)榉禽S對(duì)稱. 由圖13可以看到, 當(dāng)旋轉(zhuǎn)失速完全發(fā)展后, 轉(zhuǎn)子區(qū)域出現(xiàn)一大一小兩個(gè)失速團(tuán); 同時(shí), 在靜子的壓力面也出現(xiàn)了逆流區(qū), 且與轉(zhuǎn)子區(qū)域的失速團(tuán)相鄰. 最后, 在旋轉(zhuǎn)失速發(fā)展的整個(gè)過(guò)程中, 逆流區(qū)主要存在于壓氣機(jī)的葉尖區(qū)域.
壓氣機(jī)失速完全發(fā)展后, 95%葉高的熵分布如圖14所示, 高熵區(qū)域?yàn)槭賵F(tuán), 可以清晰地看到一大一小兩個(gè)失速團(tuán), 這兩個(gè)失速團(tuán)的結(jié)構(gòu)幾乎不隨時(shí)間變化.
圖14 轉(zhuǎn)動(dòng)到第26.7圈時(shí)95%葉高的熵分布Fig. 14 Distribution of entropy at 95% span at 26.7th revolution
ASPAC和TURBO的模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好. 尖峰型擾動(dòng)幾乎同時(shí)發(fā)生在不同的轉(zhuǎn)子葉片通道上. 然后, 擾動(dòng)逐漸增大, 最終形成的旋轉(zhuǎn)失速團(tuán)轉(zhuǎn)速為壓氣機(jī)轉(zhuǎn)速的46.41%, 同時(shí)失速區(qū)的流動(dòng)結(jié)構(gòu)幾乎不隨時(shí)間變化. 在尖峰型擾動(dòng)后, 逆流從轉(zhuǎn)子吸力面蔓延到轉(zhuǎn)子壓力面, 同時(shí)也影響到了靜子. 在整個(gè)過(guò)程中, 逆流主要存在于葉尖附近.