趙紅亮,張文強,邱佳慧,張敏,杜娟,*,聶超群
1.華北電力大學 能源動力與機械工程學院,北京 102206
2.北京理工大學 機電學院,北京 100081
3.中國科學院 工程熱物理研究所 數字孿生研究中心,北京 100190
4.中國科學院 先進能源動力重點實驗室,北京 100190
5.中國科學院 輕型動力創(chuàng)新研究院,北京 100190
6.中國科學院大學 工程科學學院,北京 100049
航空發(fā)動機實際運行(如戰(zhàn)斗機發(fā)射導彈、飛機編隊飛行、艦載機在航母上垂直/短距起降、客機飛過閃電區(qū)等)中都會在發(fā)動機進口產生溫度畸變。吸入高溫燃氣或者蒸汽后會降低發(fā)動機穩(wěn)定裕度,易于引發(fā)失穩(wěn)[1-3]。研究總溫畸變對發(fā)動機(壓氣機)的穩(wěn)定性和性能影響具有重要意義。為此,國內外的研究機構和研究人員們利用試驗、理論分析和數值仿真手段對發(fā)動機(壓氣機)進氣總溫畸變展開了大量研究工作并積累了豐富成果。
試驗方面,美國航空航天局劉易斯研究中心的Rudey 和Antl[4]利用氫氣燃燒器產生總溫畸變,對某雙軸渦扇發(fā)動機進行了總溫畸變試驗;試驗結果表明高壓壓氣機是對溫度畸變最敏感的部件,總溫畸變強度和高溫區(qū)范圍均會造成穩(wěn)定邊界下移。中國燃氣渦輪研究院的葉培梁和劉大響[5]對某小涵道比渦扇發(fā)動機進行了周向總溫畸變下穩(wěn)定性評定試驗;試驗確定了對該發(fā)動機的穩(wěn)定性影響最嚴重的一組溫度畸變參數,發(fā)現(xiàn)失速首先發(fā)生在高壓壓氣機中。
理論分析方法是基于壓氣機的數學模型把進氣溫度畸變作為輸入參數,利用模型的輸出參數實現(xiàn)快速分析壓氣機性能變化。目前研究進氣畸變對壓氣機穩(wěn)定性影響的模型有平行壓氣機模型[6]、激盤模型[7]、體積力模型[8-9]等。美國阿諾德工程發(fā)展中心的Hale 等[10-12]基于體積力思想,將通過流線曲率法獲得的葉輪機械源項耦合到三維歐拉求解器中,發(fā)展了三維計算程序TEACC(Turbine Engine Analysis Compressor Code),利用該程序定量評估了進氣畸變對某三級軍用風扇性能和穩(wěn)定性的影響。郭晉等[13]通過結合三維體積力和經典二維平行壓氣機模型提出了一種混合保真度計算模型;利用該模型以較低計算成本定量評估了復雜進氣畸變對某大涵道比發(fā)動機的性能影響,并捕捉到內部流場的三維特征。李秋實等[14]利用體積力模型將葉片排和進氣道-風扇相互作用等作為源項耦合到COMSOL-CFD 程序中,建立了進氣道-風扇耦合模擬方法;該模型較好地模擬出了進氣道-風扇系統(tǒng)的特性性能及流場性能。葉巍等[15-16]通過前人提出的畸變傳遞模型[17]研究了進氣溫度畸變對某發(fā)動機穩(wěn)定性和性能的影響;結果表明穩(wěn)定裕度的損失更多取決于畸變強度,而不是畸變的周向范圍;較高轉速下發(fā)動機抗溫度畸變能力明顯強于中低轉速;在發(fā)射的溫度畸變強度下發(fā)動機已失去氣動穩(wěn)定性,必須使用防喘系統(tǒng)。張百靈等[18]以劉易斯研究中心設計的5 級壓氣機的前3 級為研究對象建立了一個壓縮系統(tǒng)的模型,結合體積力方法,通過求解非穩(wěn)態(tài)三維歐拉方程較為準確地計算了周向總溫畸變的傳遞特性;分析發(fā)現(xiàn)周向總溫畸變經過壓氣機后強度衰減了19%,高溫區(qū)的周向位置偏移約80°,且總溫畸變強度在轉子內發(fā)生局部增大現(xiàn)象,第1 級轉子失速的可能性最大。
數值仿真方面,Weston 等[19]對進氣總壓畸變條件下的某3 級風扇展開全環(huán)數值仿真,研究表明經過風扇后總壓畸變引發(fā)了總溫畸變,越接近失速工況點,總溫畸變強度越大,這與轉子葉片周圍輸入的功有關,運行在近失速工況點時轉過全環(huán)過程中第1 級轉子的單個葉片做功量變化范圍可達25%。李志平等[20]利用數值模擬手段對某單級軸流壓氣機進行了不同畸變強度的壓力-溫度組合畸變對穩(wěn)定性和性能影響的研究;通過近失速工況點轉子葉頂區(qū)二維流線的分析發(fā)現(xiàn)組合畸變下壓氣機失穩(wěn)機理仍是泄漏流從轉子葉片前緣溢出。尤延鋮等[21]總結概述了國內外在航空發(fā)動機進氣總溫畸變方面的研究工作,包括進氣總溫畸變來源、溫度畸變發(fā)生裝置及各類總溫畸變評定標準,最后列舉分析了針對總溫畸變的發(fā)動機防喘控制技術優(yōu)缺點,指出對進氣總溫畸變產生機理、傳遞和影響等研究急需加強。
通過上述研究發(fā)現(xiàn)在進氣總溫畸變試驗方面主要是針對整機展開研究,但由于整機試驗成本昂貴、周期漫長等問題,近年來公開的文獻較少;采用理論分析可實現(xiàn)進氣畸變對發(fā)動機(壓氣機)穩(wěn)定性和性能影響的快速預測,但是對內部流場細節(jié)研究有所不足;此外,目前總溫畸變下的高精度非定常計算工作較少,主要聚焦于壓氣機進口總壓畸變引發(fā)的總溫畸變或總壓總溫復合畸變條件下,畸變區(qū)在壓氣機內部的傳播規(guī)律及對壓氣機性能的影響。而針對進氣總溫畸變條件下壓氣機性能變化的背后物理機理和失速現(xiàn)象有待進一步深入研究。相關試驗數據[22]表明徑向溫度畸變影響較小,通常忽略不計。因此本文對Darmstadt 跨聲單級壓氣機展開了周向總溫畸變條件下全環(huán)非定常數值模擬研究,分析失速過程中壓氣機內部流場變化,以期揭示總溫畸變下該壓氣機發(fā)生失速的動力學機理。
研究對象為德國Darmstadt 跨聲單級壓氣機試驗臺(如圖1 所示)。該壓氣機由德國MTU 航空發(fā)動機公司和德國航空太空中心設計并優(yōu)化[23-25],轉子葉片在輪轂附近彎曲度較高且在葉頂附近很薄,整個工作范圍內靜子葉片沒有流動分離跡象且出口角分布均勻,是渦扇高壓壓氣機前級的典型代表。壓氣機采用800 kW 直流驅動,通過變速器改變軸速度,最大轉速為20 500 r/min,最大扭矩為350 N·m。壓氣機測量截面有級進口面(ME15)、轉子進口面(ME21)、級出口面(ME30),壓氣機性能特性是根據級進口面和級出口面的面平均參數計算得到的,已公布的試驗數據包含100%和65%設計轉速結果。表1總結了試驗臺的關鍵參數,更詳細的數據見文獻[26-27]。
表1 Darmstadt 跨聲壓氣機的關鍵參數Table 1 Key parameters of Darmstadt transonic compressor
圖1 Darmstadt 跨聲壓氣機試驗臺Fig.1 Test rig of Darmstadt transonic compressor
通過Ansys CFX 軟件求解雷諾平均Navier-Stokes 方程。計算設置采用了標準k-ω湍流模型[28],壁面為絕熱無滑移,流體介質為理想氣體。對于均勻來流,采用具有周期性邊界條件的單通道定常計算,進口邊界條件設置為公布的試驗測量值[27],即進口總壓為99 kPa、進口總溫為295 K、湍流強度為4%、湍流長度尺度為0.09 m。轉靜交界面類型為Stage(Mixing-Plane)。對進氣總溫畸變采用全環(huán)非定常計算,進口總溫設置為180°周向范圍的500 K 高總溫區(qū)和295 K 的低總溫區(qū)(如圖2 所示,圖中θ為周向位置),轉靜交界面類型為Transient Rotor Stator。出口邊界條件采用固定背壓為30 kPa 的噴嘴,通過調整噴嘴喉部大小獲取不同工況點的特性,大量研究[29-31]已證實此出口邊界條件的設置可在壓氣機較小流量工況點使計算穩(wěn)定,適合壓氣機失穩(wěn)計算。求解格式均為高階求解模式,非定常計算采用隱式雙時間步,轉子葉片轉過一個葉片通道所需的物理時間步長為80 步,每個物理時間步長的虛擬時間步長為10 步。
圖2 全環(huán)計算域示意圖Fig.2 Schematic of annulus computation domain
葉片網格是通過NUMECA/AutoGrid5 生成,采用O4H 拓撲,根據文獻[32-33],網格最大壁面距離y+<3。針對單通道計算域進行了粗糙、中等、精細3 種不同網格方案的模擬測試。3 組網格方案的網格節(jié)點在流向、周向、徑向具體分布情況見表2。將100%轉速下計算特性和試驗數據進行對比,如圖3 所示,可見粗糙網格雖可求解出更小流量點,但整體特性上級總壓比和等熵效率低于試驗數據;隨網格量增大,特性線逐漸向右上角移動,折合流量、級總壓比和等熵效率的數值更大,與試驗值吻合較好;同時中等網格和精細網格求解的失速裕度及特性差距很小,沒有明顯變化。從求解精度和計算資源角度綜合考慮,最終選擇了中等網格量,即單通道計算域的網格量約為247 萬,全環(huán)計算域的網格量約為5 900 萬。
表2 網格流線、周向、徑向分布Table 2 Mesh distributions at streamline, circumferential and radial directions
圖3 不同網格下壓氣機特性與試驗數據[27]對比Fig.3 Comparison of compressor characteristics for different meshes with test data[27]
驗證過程是將數值模擬結果和試驗數據進行對比,量化數值仿真的精度,也是進行研究的第一步。65%和100%轉速壓氣機級總壓比和等熵效率與試驗數據的對比如圖4 所示,可看出數值模擬得到的特性和試驗值吻合較好,同時能準確捕捉不同轉速下壓氣機的近失速點和堵塞點。
圖4 數值結果與試驗數據[27]對比Fig.4 Comparison of numerical results with test data[27]
為進一步驗證數值模擬方法的準確性,圖5 給出了100%轉速下最高效率(Peak Efficiency,PE)點和近失速工況(Near Stall,NS)點的級出口面參數徑向分布對比,圖中mC為折合流量。不同工況點級出口面的總壓比和總溫比徑向分布對比如圖5(a)和圖5(b)所示,模擬值的總壓比在葉根和葉尖處相較于試驗數據略大,而總溫比整體吻合較好。從PE 點到NS 點級總壓比徑向分布變得不均勻,NS 點靠近葉根和葉頂區(qū)域,均有較大的總壓分布;而級總溫比在60%葉高以上溫梯度大且變化范圍較大。盡管數值模擬結果與試驗數據在部分葉高有一定偏差,但結果與其他研究人員采用不同求解器的計算結果[27,33]趨勢基本一致。
圖5 級出口面(ME30)參數與試驗數據[27]對比Fig.5 Comparison of parameters at stage exit (ME30)with test data[27]
試驗臺的靜子輪轂空腔會存在泄漏流,根據文獻[34],考慮0.5%的泄漏流量后數值模擬結果和試驗數據在徑向分布上的對比會更好,工作中未考慮泄漏流的影響。
最后從內部流場角度對3 個不同的典型工況(近失速點、最高效率點和近堵塞點)的靜壓參數進行對比。數值模擬結果提取的99.9%葉高靜壓云圖與試驗轉子葉頂區(qū)的壁面靜壓云圖[27]對比如圖6 所示,圖中LE 為前緣(Leading Edge),TE 為尾緣(Trailing Edge),可見數值模擬能精確捕捉脫體激波、葉頂泄漏流的軌跡和通道激波的位置。圖3~圖6 的對比充分驗證了數值仿真設置的正確性。
圖6 數值模擬所得內部流場與試驗數據[27]的對比Fig.6 Comparison of internal flow fields obtained from numerical simulation with test data[27]
文獻[35]的數值研究表明誘發(fā)該壓氣機失速的先兆波為突尖型,計算模擬捕捉到前緣溢出和尾緣回流的流動特征。在失速模擬中從失速前最后一個穩(wěn)定解開始節(jié)流。失速特征如下:在約1 rev 后壓比上升,質量流量開始下降;在3 rev 左右達最小值后質量流量和壓比又開始增加,在約6 rev 后達最終失速狀態(tài)。通過失速模擬中攻角變化的研究確定了該壓氣機失速擾動發(fā)生在轉子葉頂區(qū)域,且在失速先兆階段動葉周圍的流動發(fā)生前緣溢出和強烈的尾緣回流現(xiàn)象,為突尖型失速先兆。模擬過程中失速團的周向傳播速度約為1/3 轉子轉速。文獻[36]的試驗研究表明該壓氣機以突尖型先兆啟動失速,形成一個單獨的失速團并繞環(huán)旋轉。試驗中,100%轉速下從穩(wěn)定極限點關小閥門。在約5 rev 時失速團尺寸最大,對應局部流量最小,傳播速度約為0.44 倍轉子轉速;隨后失速團尺寸開始減小,約15 rev 時失速團尺寸達最小,對應局部流量最大,傳播速度約為0.52 倍轉子轉速。試驗觀測到隨流量減小失速團尺寸增大、傳播速度變慢的特性,失速團的尺寸和傳播速度隨壓氣機接近穩(wěn)定的失速點而震蕩。
為分析進氣總溫畸變影響壓氣機性能的物理機理,表3 展示了PE 點(折合流量為14.58 kg/s)壓氣機進口面高溫畸變區(qū)和低溫非畸變區(qū)關鍵參數的對比,可見高溫畸變區(qū)總溫為500 K,低溫非畸變區(qū)總溫為295 K;總壓均為99 kPa;高溫區(qū)氣流密度為0.656 8 kg·m-3,低溫區(qū)氣流密度1.109 4 kg·m-3;高溫區(qū)軸向速度為141 m·s-1,而低溫區(qū)軸向速度為109 m·s-1;高溫區(qū)密流為93 kg·m-2·s-1,低溫區(qū)密流為121 kg·m-2·s-1??梢娍倻鼗儗簹鈾C的主要影響為高溫區(qū)氣流密度大幅減小,軸向速度增大。同時高溫畸變區(qū)工作在低折合轉速線上,因此總溫畸變下壓氣機工作特性點往左下角移動。
表3 高、低溫區(qū)關鍵參數對比Table 3 Comparison of key parameters in high and low temperature regions
為總結總溫畸變下壓氣機內部流場隨不同工況點的變化規(guī)律,對不同周向位置處時均總壓徑向分布情況進行分析。在總溫畸變非定常計算的PE 點(折合流量為14.58 kg/s)和NS 點(折合流量為13.46 kg/s)的轉子出口面,每隔45°周向角度提取總壓徑向數據進行對比。周向角度位置定義如圖7 所示,背景為NS 點轉子出口面時均總壓云圖。
圖7 轉子出口面周向角度位置Fig.7 Position of circumferential angle at rotor exit surface
圖8 (a)、圖8(b)分別展示了總溫畸變下PE工況點0°~180°和180°~360°不同周向位置的總壓徑向分布情況,可見在0°、45°、90°、135°和180°的高溫區(qū)周向位置,總壓徑向分布逐漸變窄。不同周向位置的總壓徑向分布形式比較類似,葉中區(qū)域分布均勻,而葉頂葉根區(qū)域由于端壁邊界層、角區(qū)等導致流動損失進而使總壓存在較大梯度;在180°、225°、270°、315°和360°的周向位置,總壓徑向分布逐漸變寬。
圖8 PE 點轉子出口面不同周向位置處總壓徑向分布Fig.8 Radial distributions of total pressure at different circumferential positions of rotor exit for PE condition
圖9(a)、圖9(b)分別展示了總溫畸變下NS工況點0°~180°和180°~360°不同周向位置的總壓徑向分布情況,可見在0°、45°、90°、135°和180°的高溫區(qū)周向位置,總壓徑向分布呈現(xiàn)逐漸變窄的趨勢。在180°、225°、270°、315°和360°的低溫區(qū)周向位置,總壓徑向分布呈現(xiàn)逐漸變寬的趨勢且分布較為不均勻。轉子出口面總壓云圖如圖7所示,在225°周向位置處,葉根區(qū)域的總壓較小,而葉頂區(qū)域存在一個高總壓區(qū),總壓徑向分布不均勻。這是由于畸變引起氣流角沿周向分布不同,且沿徑向也發(fā)生遷移,從而加功量不同,導致總壓在不同周向位置的徑向分布差異。
圖9 NS 點轉子出口面不同周向位置處總壓徑向分布Fig.9 Radial distribution of total pressure at different circumferential positions of rotor exit for NS condition
圖10 為不同葉高處和具有相同噴嘴喉部面積的、均勻來流的、絕對周向氣流角的周向分布對比。定義逆轉子旋轉方向為絕對周向氣流角正方向。結合圖11 中的速度三角形,可見在0°~180°范圍絕對周向氣流角α由正值減小到負值,相對氣流角β逐漸減小,因此攻角減小,葉片載荷減小,總壓減小。在180°~360°范圍,絕對周向氣流角α由負值增大到正值,相對氣流角β逐漸增大,因此攻角增大,葉片載荷增大,總壓增大。并且對于該壓氣機而言,總溫畸變下周向氣流角變化范圍在葉根區(qū)域最大,從10%葉高到90%葉高逐步減小。以上研究結果表明,在轉子由低溫非畸變區(qū)轉入高溫畸變區(qū)時攻角最大,因此猜測此處可能最先發(fā)生失穩(wěn)。為驗證這一猜想,接下來進一步對失速過程進行了分析。
圖10 轉子進口絕對周向氣流角Fig.10 Absolute circumferential flow angle at rotor inlet
圖11 速度三角形中絕對周向氣流角Fig.11 Absolute circumferential flow angle in velocity triangle
進氣總溫畸變條件下,壓氣機的級總壓比特性曲線如圖12 所示。從圖12 中可見由一個穩(wěn)定解過渡到失速狀態(tài),總壓比和折合流量先略微升高,然后迅速下降,再回升,以此循環(huán)往復呈現(xiàn)周期性變化特征。
圖12 總溫畸變下壓氣機特性Fig.12 Compressor performance characteristics with total temperature distortion
進氣總溫畸變條件下失速過程中壓氣機進口物理流量隨轉數的變化情況如圖13 所示,物理流量和折合流量的換算關系為
圖13 總溫畸變下失速過程壓氣機進口物理流量Fig.13 Physical flow rate at compressor inlet during stall process under total temperature distortion
式中:mC為折合流量;m為物理流量為級進口總壓;為級進口總溫。
從10 rev 開始由一個穩(wěn)定的工況解輕微減小噴嘴喉部面積。由圖13 可見從11 rev 開始物理流量迅速從10.8 kg/s 下降到7.6 kg/s,隨后升高到11.3 kg/s,往復變化,且物理流量的震蕩幅度有輕微減小,變化趨勢與圖12 中壓氣機特性變化相符合。
圖14 展示了不同轉數t下轉子出口面的軸向速度分布云圖,可見在11 rev 軸向速度沿周向分布較均勻,未出現(xiàn)回流;在13 rev 出現(xiàn)2 個較大失速團和較小失速團(黑虛框),且低軸向速度區(qū)域近乎占據了整個轉子出口面一半周向范圍,此時壓氣機已經處于失速狀態(tài);轉到16 rev 和20 rev時失速團合并為1 個,且占據的周向范圍大幅縮??;轉到23 rev 和27 rev 時失速團變?yōu)?~2 個較大失速團和若干小失速團,大約占據整個轉子出口面1/3 周向范圍。結合圖13 物理流量隨轉數時間的震蕩變化情況看,當物理流量較?。ㄈ?3 rev 和27 rev)時失速團的數量較多、占據范圍較大,堵塞程度較嚴重;當物理流量較大(如16 rev 和20 rev)時失速團數量較少、占據范圍較小,堵塞程度較輕;由此計算出流量振蕩頻率約為75.0 Hz。
圖14 失速過程的軸向速度云圖Fig.14 Axial velocity contours during stall process
為進一步分析失速過程中失速團的體積大小和空間分布,圖15 展示了不同轉數瞬態(tài)解的轉子域中軸向速度為0 的三維等值面分布云圖。在全環(huán)域瞬態(tài)解中生成軸向速度為0 的等值面,這個等值面可看作失速團的邊界,等值面圍成的空間大小可看作失速團的體積,這種方法是由Zhang 和Vahdati[31]提出的。從圖15 中可清晰看到進氣總溫畸變條件下壓氣機失速過程中失速團的空間分布集中在葉頂區(qū)域。隨著轉子旋轉時間推進,失速團的空間位置在葉頂沿周向移動,失速團的體積也在不斷變化。
圖15 失速過程軸向速度為0 的等值面云圖Fig.15 Iso-surface of zero axial velocity contours during stall process
為進一步確定失速先兆最先出現(xiàn)的周向位置和先兆波類型,在數值計算設置中沿周向均勻布置了12 支絕對靜壓探針監(jiān)測轉子前緣的靜壓值變化。探針距離葉片前緣為10%轉子軸向弦長,探針的徑向位置為99%葉高,探針布局方式和序號如圖16 所示。
圖16 靜壓探針的周向布局Fig.16 Circumferential arrangement of static pressure probes
圖17 展示了12 支絕對靜壓探針的壓力信號隨時間的波動情況,可看出CH10 探針(即壓氣機轉子從低總溫非畸變區(qū)轉入高總溫畸變區(qū)的周向位置處)最先監(jiān)測到明顯壓力波動,隨后順轉子旋轉方向的探針壓力波動幅值增大,表明轉子葉片前緣的非定常波動增強,壓氣機由穩(wěn)定的近失速工況點轉入失速狀態(tài)。失速先兆的周向傳播速度約為88.9%轉子轉速,失速初期失速團的傳播速度約為66.0%轉子轉速。圖18(a)和圖18(b)展示了11.0 rev 和11.5 rev 時98%葉高的速度矢量和軸向速度云圖,可見在11.0 rev 時只觀測到葉頂泄漏流從轉子前緣輕微溢出;在11.5 rev 時可看到葉頂泄漏流在轉子前緣明顯溢出和尾緣回流現(xiàn)象,此時葉片通道發(fā)生堵塞,最終形成失速擾動。這符合Vo 等[37]通過單/多通道壓氣機定常/非定常數值計算提出的突尖型失速先兆的準則:葉頂泄漏流從轉子葉片前緣溢出和尾緣跨通道回流。因此總溫畸變下誘發(fā)旋轉失速的先兆波為突尖型失速先兆。
圖17 周向靜壓探針的壓力信號Fig.17 Pressure signals of circumferential static pressure probes
圖18 不同時刻速度矢量和軸向速度云圖Fig.18 Velocity vector and axial velocity contour at different times
對德國Darmstadt 跨聲壓氣機進行了均勻進氣條件下單通道定常數值模擬,數值結果與試驗數據吻合較好,驗證了計算設置的準確性;進行了總溫畸變條件下的壓氣機全環(huán)非定常數值模擬,對失速過程進行了研究,得到以下結論。
1) 在周向范圍180°、高溫畸變區(qū)500 K 的周向總溫畸變條件下壓氣機總壓比-流量特性大幅下降,在失速過程中特性曲線呈周期性震蕩現(xiàn)象,當物理流量較小時對應壓比較小,失速團數量較多、占據范圍較大,堵塞程度較重;當物理流量較大時對應壓比較大,失速團數量較少、占據范圍較小,堵塞程度較輕。
2) 總溫畸變下不同周向位置處的總壓徑向分布不同。順轉子葉片旋轉方向總壓在高溫區(qū)逐漸減小,在低溫區(qū)逐漸增大,且在葉根區(qū)域變化最為明顯。
3) 總溫畸變下最先發(fā)生壓力擾動的周向位置是轉子轉入高溫畸變區(qū)處,誘發(fā)旋轉失速的先兆波為突尖型失速先兆。
4) 總溫畸變下失速先兆的周向傳播速度約為88.9%轉子轉速,失速初期失速團的傳播速度約為66.0%轉子轉速,高于均勻進氣下失速團的傳播速度。