張博涵,王 強,謝業(yè)平
(1.北京航空航天大學能源與動力工程學院,北京100083;2.中國航發(fā)沈陽發(fā)動機研究所,沈陽110015)
發(fā)動機和進氣道是飛機推進系統(tǒng)的2大主要部件,發(fā)動機要吸入空氣,為飛機提供不同飛行狀態(tài)所需的動力,而進氣道則需要保證在不同飛行條件下流入發(fā)動機的空氣流量和空氣品質(zhì)等,從而保證發(fā)動機正常工作。因此,進氣道與發(fā)動機之間需要有良好的適應性關系,否則將對發(fā)動機乃至整個推進系統(tǒng)的穩(wěn)定性和綜合性能產(chǎn)生嚴重影響[1]。
現(xiàn)代戰(zhàn)斗機機身和推進系統(tǒng)的整合涵蓋了從飛機穩(wěn)定性到控制,再到進氣道/發(fā)動機相容性的多種問題[2]。為此,在試驗和數(shù)值模擬過程中需要使用大量試驗和計算資源,并應用各種分析和輔助計算工具。其中進氣道/發(fā)動機一體化的試驗,尤其是關于進氣道與發(fā)動機相容性的試驗,通常需要在風洞和高空試驗設施中進行,試驗周期長,在前期準備和試驗過程中耗資巨大[3]。相比而言,利用數(shù)值模擬方法可以在試驗前對進氣道和發(fā)動機的相容性進行預先研究,且可以在一定程度上對試驗無法模擬的飛行條件進行計算分析。因此,針對試驗出現(xiàn)的問題,可以轉(zhuǎn)而尋找準確、高效的數(shù)值模擬方法,來研究進氣道和發(fā)動機的匹配問題。
對于發(fā)動機單獨的數(shù)值模擬已經(jīng)有較為成熟的計算方法,也發(fā)展了一系列快速而精確的計算程序。但由于壓氣機與進氣道氣流流動的尺度相差較大,當采用精確的計算方法進行聯(lián)合計算時,耗時長、效率低、占用的計算資源巨大。由此發(fā)展了對壓氣機部分進行?;姆椒?,現(xiàn)階段對風扇的模化方法主要采用激盤模型(Actuator Disk,AD)和徹體力模型,且經(jīng)過了多年的發(fā)展,2種模型已經(jīng)可以滿足進氣道和發(fā)動機匹配計算的精度要求。激盤模型主要是將風扇的葉片?;?個無厚度的盤面,原本受風扇作用的氣流在盤面上完成葉片對氣流的作用,此時在激盤上下游都為正常的進氣道內(nèi)流,使壓氣機與進氣道的尺度相近,提高計算效率[4-6]。而徹體力模型由于在?;瘯r將風扇部分的控制方程變?yōu)闊o黏條件,因此進氣道與?;蟮娘L扇只能通過進氣道出口的參數(shù)相互傳遞邊界條件,而無法進行聯(lián)合數(shù)值計算[7-9]。
在Joo W G[10]建立激盤模型過程中由于忽略了葉柵的厚度,即忽略了葉柵內(nèi)存儲質(zhì)量和能量的能力以及葉片的徑向力,且激盤軸向的位置選擇都會對計算結(jié)果產(chǎn)生影響。Joo討論了激盤模型的解與激盤位置的相關性,同時也研究了對出口氣流角的修正問題。本文基于Joo W G的激盤模型建立了1種分裂激盤模型(Divided Actuator Disk,DAD),應用分裂激盤模型對進氣畸變條件下的3級軸流壓氣機性能進行分析。
用進氣道和轉(zhuǎn)子流場對1個激盤模型進行簡單描述,如圖1所示。風扇被?;?個平面激盤,位于某一合適的軸向位置,例如在葉片的后緣或者在弦平面上。原本真實存在的葉柵變?yōu)檗D(zhuǎn)子上、下游的流場,葉片上下游的流場參數(shù)靠強加在激盤兩側(cè)的邊界條件來耦合,從而代替真實風扇對氣流的作用。即假想風扇葉片對氣流的作用只存在于?;钠矫婕けP,激盤的上下游均為正常的進氣道內(nèi)流。
圖1 激盤模型
將這樣的激盤的上下游聯(lián)系起來需要激盤的邊界條件滿足以下5個條件[10]:質(zhì)量守恒;徑向動量守恒;滯止焓守恒;指定的相對出口氣流角;指定的熵增。第1、3個條件忽略了葉片排存儲質(zhì)量和能量的能力,第2個條件則忽略了葉片作用在流體上的徑向力,第4個條件是將風扇對氣流的加功作用通過指定的相對出口氣流角確定。
由于常規(guī)激盤模型的假設,在使用激盤模型時需要根據(jù)具體情況對模型進行修正,因此存在一定的局限性。激盤模型的定義是將整個風扇模化成1個無厚度的平面激盤,風扇對氣流的作用無法準確、完全地由激盤模擬,因此需要對激盤模型進行相應的修正。當平面激盤選擇在葉柵的不同位置時,激盤的徑向長度發(fā)生變化,與在葉柵后緣和真實風扇徑向距離不一致,造成與真實加功量不一致。最終影響激盤模型計算的精確性。
為此,本文在Joo W G建立的激盤模型基礎上,將1個平面激盤分開,分別貼合到葉柵的前后緣,發(fā)展了分裂激盤模型,可以有效地消除因徑向力假設而產(chǎn)生的誤差。
分裂激盤模型的網(wǎng)格邊界分別與葉片前、后緣的掃掠形狀一致(如圖2所示),可避免激盤盤面形狀以及位置引起的誤差,同時也與真實風扇前后氣流的流動情況有一定相關性。根據(jù)特定假設,可將葉片前、后緣網(wǎng)格單元面上的參數(shù)建立起相關控制方程,參與流場的數(shù)值迭代。
假設葉柵內(nèi)氣流流線滿足以下條件,且激盤后的網(wǎng)格也是依此關系由激盤前網(wǎng)格自動生成的。
圖2 分裂激盤模型總體
式中:R為激盤盤面單元格處流線的半徑位置;Rmax、Rmin分別為機匣和輪轂的半徑。
這樣前、后激盤盤面位置的參數(shù)可沿同一條流線建立方程。
由滯止熵守恒條件可得
式中:I為相對坐標系下的滯止熵;ω為風扇角速度;下標1、2分別指葉片前、后緣位置。
由流量守恒條件可得
式中:A為葉片進、出口網(wǎng)格單元面面積;Wˇ為旋轉(zhuǎn)坐標系下A所在網(wǎng)格體的速度向量;dˇ為葉片進、出口網(wǎng)格單元面的法向單位向量,如圖3所示。
圖3 分裂激盤參數(shù)
式中:α為盤面位置流線與軸線的夾角;β為相對坐標系下的葉片出口氣流角;dz為dˇ軸向分量;Wz為相對坐標系下速度的軸向分量,與絕對速度軸向分量Vz相同。
相對坐標系下的滯止熵又可表示為
由此得到關于密度ρ的方程為
式中:m為通過網(wǎng)格單元面的單位流量。
由式(7)求解密度,然后得到流場中其他參數(shù)的值。
同時由于轉(zhuǎn)動的葉片與氣流的相對速度可能會大于聲速,這種情況下葉柵會發(fā)生喘振,前面描述的分裂激盤模型中不包括葉片的這種流動堵塞作用,因此無法預測葉片通道的堵塞。本文使用1種基于堵塞葉片中2維流動特性的簡單模型來確定堵塞流量。2維葉片通道中堵塞流動如圖4所示。假設整個葉片通道流動條件是均勻的,則可簡單計算出喉道的有效面積。
圖4 在葉片通道內(nèi)的堵塞流動
忽略聲速線上游的損失,可得到通過葉片通道的最大質(zhì)量流量,流量守恒為
式中:A1和A*為激盤入口和喉道橫截面積處流管軸向的橫截面積;上標rel指相對于葉片;F*為合適的可壓流函數(shù)。
由于前面假設了葉片通道的流線滿足式(1),則該流量守恒公式可直接對流線兩側(cè)單元格參數(shù)進行控制。
分裂激盤模型的邊界條件都通過這個堵塞模型式(8)規(guī)定的最大流量進行判斷和控制。如果流量比堵塞流量小,那么分裂激盤的進口邊界條件由原模型給出;否則進口單元格的流量就固定為堵塞的值,這時通過使用堵塞的流量計算得到出口單元格的流動變量,直接指定出口單元格的流量并繼續(xù)計算。
本文以NASA Rotor 67風扇轉(zhuǎn)子作為研究對象[11]。Rotor 67轉(zhuǎn)子設計用于軸向流動,工作時作為孤立轉(zhuǎn)子,沒有進口導葉和靜葉。由于葉型數(shù)據(jù)公開且有詳細的試驗結(jié)果,因此有大量針對該轉(zhuǎn)子的理論和數(shù)值研究工作[12-13]。其基本參數(shù)見表1,具體的幾何參數(shù)與試驗數(shù)據(jù)參考文獻[14]。
表1 Rotor 67轉(zhuǎn)子基本參數(shù)
根據(jù)文獻[14]中的代碼自動生成葉片與流道的幾何參數(shù),如圖5所示。
圖5 Rotor 67和流道的幾何模型
由于激盤模型原則上是基于3維數(shù)值模擬的簡化模擬方法,因此本文中驗證部分分裂激盤模型的參數(shù)依據(jù)A.Arnone的計算結(jié)果[15]。通過使用可變系數(shù),采用隱式殘差平滑和完全多重網(wǎng)格的方法,提高了4階Runge-Kutta法的計算效率。使用多重網(wǎng)格方法對NASA Rotor 67轉(zhuǎn)子計算,結(jié)果如圖6所示。從圖中可見,模擬結(jié)果與NASA公開的試驗結(jié)果一致性較好。
根據(jù)Arnone方法的計算結(jié)果,選取風扇在峰值效率點工作時的出口相對氣流角徑向分布和絕熱效率,建立分裂激盤模型。其中葉片前、后緣和相對出口氣流角均使用二次曲線進行擬合,沿軸向截面的網(wǎng)格如圖7所示。
圖6 NASA rotor 67轉(zhuǎn)子在設計轉(zhuǎn)速時的性能
圖7 分裂激盤模型沿流向截面的網(wǎng)格
邊界條件設置與3維模擬一致,壓力入口總壓為101325 Pa,總溫為288.15 K,壁面速度無滑移。設置不同出口反壓,模擬結(jié)果如圖8所示。從圖中可見,模型計算的壓比整體略高于全3維模擬結(jié)果,從接近失速到接近堵塞范圍內(nèi)與3維模擬的誤差最大約3%,在接近最高效率點的一致性很好。在接近失速以及失速時由于模型的局限性,模型的計算結(jié)果與全3維模擬結(jié)果偏差較大,這是不可避免的。在風扇正常工作范圍內(nèi),模型的壓比稍大于3維模擬的,究其原因主要是模型采用的絕熱效率為3維模擬中的峰值效率,所以模型計算的壓比整體偏高[16]。但由于峰值絕熱效率是評價風扇性能的重要參數(shù),因此在考慮了模型針對個體的精確性和普適性后,堅持選擇風扇的峰值效率作為模型的輸入?yún)?shù)。
圖8 壓比特性對比
在最高效率點的模型方法與3維模擬計算的壓比沿相對葉高的分布如圖9所示。從圖中可見,在葉根和葉尖處模型的結(jié)果與3維模擬結(jié)果偏差較大,主要是由于激盤模型對邊界層流動的處理方法有局限性,且葉尖部分在3維模擬中考慮了轉(zhuǎn)子的葉尖間隙,而在模型計算中沒考慮,以上多種原因共同導致壓比偏差。
本文模型與Joo W G的激盤模型對Rotor 67轉(zhuǎn)子的壓比特性模擬結(jié)果對比如圖10所示,其中激盤模型的激盤位置位于風扇出口。從圖中可見,分裂激盤模型改進了激盤模型的結(jié)果。風扇在正常工作范圍內(nèi),分裂激盤模擬的壓比特性更接近于3維模擬結(jié)果,主要由于分裂激盤模型的邊界基本擬合了風扇進、出口的幾何形狀,而激盤模型的邊界為環(huán)形平面,引起相對出口氣流角的誤差,進一步產(chǎn)生風扇加功量不同的問題。而且分裂激盤模型不需要選擇激盤盤面的位置和后續(xù)的修正,相對于激盤模型能夠更加準確地模擬風扇的特性。
圖9 最高效率點下總壓比徑向分布
圖10 分裂激盤模型與激盤模型的壓比特性對比
由于在建立分裂激盤模型時,控制方程的時間復雜度并沒有成比例增加,監(jiān)控計算過程的殘差如圖11所示。從圖中可見,2種模型殘差的收斂速度量級相當。由于本文模型對3維定常流場的N-S方程采用3階精度有限體積法離散,隱式LU-SGS(虛)時間推進法迭代求解。對流通量采用Roe格式計算,擴散通量采用中心差分格式計算。湍流模型為CG K-epsilon線性湍流模型,模型參數(shù)輸運方程采用混合隱式迭代加解析方法求解,因此可以通過K和epsilon的殘差量級判斷計算是否收斂。圖中K與epsilon的殘差量最終都在10-4量級以下,可認為計算達到收斂。
圖11 Rotor67轉(zhuǎn)子計算殘差對比
本文選用的多級壓氣機模型為NASA Compressor 74A。該壓氣機為5級核心壓氣機,由于文獻[17]中研究了前3級的總體性能,因此本文也僅對壓氣機的前3級進行模擬和分析。壓氣機基本結(jié)構(gòu)如圖12所示,基本設計參數(shù)見表2。
圖12 NASA Compressor 74A壓氣機前3級葉型流道結(jié)構(gòu)
表2 壓氣機前3級基本參數(shù)
分別對無畸變和畸變條件下的壓氣機氣動穩(wěn)定性進行分析,以檢測分裂激盤模型對于穩(wěn)態(tài)總壓畸變的模擬能力。其中畸變條件的進口總壓給定如圖13所示。為模擬進氣道出口畸變情況,設置為既存在徑向畸變又存在周向畸變。無畸變區(qū)總壓為101325 Pa,低壓畸變區(qū)總壓為無畸變區(qū)的90%。
圖13 進口總壓分布
壓氣機前3級葉柵進、出口的總壓分布如圖14所示。從圖中可見,進口的總壓畸變對葉柵區(qū)域的總壓分布產(chǎn)生影響,畸變區(qū)氣流的徑向分布引起壓氣機內(nèi)氣流徑向重新分布;葉尖周向總壓的分布也顯示出經(jīng)過葉柵后,氣流發(fā)生摻混;經(jīng)過前3級葉柵后,總壓畸變有一定改善。
圖14 壓氣機前3級進、出口的總壓分布
在設計轉(zhuǎn)速下進口總壓畸變與無畸變的壓氣機流量-壓升特性曲線對比如圖15所示。從圖中可見,總壓畸變導致壓氣機總壓比減小,且畸變時壓氣機的特性曲線整體向失速邊界靠近,表明該模型能模擬總壓畸變對壓氣機性能的影響。
圖15 壓氣機壓比特性對比
由于已經(jīng)驗證了分裂激盤模型可以很好地模擬風扇正常工作時的特性,分別對進氣道和風扇的工作特性進行分析,得到共同匹配點后,用分裂激盤模型代替3維風扇模型進行進氣道-風扇聯(lián)合流場的數(shù)值模擬研究。
本文所研究的進氣道幾何形狀已經(jīng)確定,且斜板可根據(jù)飛行條件調(diào)整角度,如圖16所示。第1道斜板前和第2、3道斜板間的角度是固定的,只有第1、2道斜板間可根據(jù)來流條件不同與發(fā)動機工作狀態(tài)不同繞接縫轉(zhuǎn)動[18]。
對不同來流馬赫數(shù)和斜板角度的進氣道工作特性進行分析,來流靜壓均為22700 Pa。選擇的典型進氣道斜板角度與來流馬赫數(shù)條件分別為:(1)來流馬赫數(shù)為1.33,斜板角度為-2.7°;(2)來流馬赫數(shù)為1.33,斜板角度為 2°;(3)來流馬赫數(shù)為 2.25,斜板角度為13.14°。得到接近臨界狀態(tài)進氣道各項參數(shù)見表3。
圖16 斜板角度為-2°的進氣道幾何截面
表3 接近臨界狀態(tài)進氣道各項參數(shù)
與該變幾何進氣道匹配的是1.5級風扇,即1排導向葉片、1排動葉和1排靜葉。單級風扇定常模擬所用的計算網(wǎng)格如圖17所示。在不同轉(zhuǎn)速的計算過程中,導向葉片與軸向的夾角始終為 0°。
圖17 1.5級風扇定常計算所用網(wǎng)格
根據(jù)圖17葉片和風扇流道的幾何模型,建立計算域,在進口導葉前和靜葉后延伸2個弦的長度作為邊界條件的緩沖。使用RANS方法進行數(shù)值模擬,由于計算方法已廣泛用于風扇3維計算,因此無需對計算結(jié)果與計算方法進行驗證。得到流場的詳細參數(shù)與風扇的絕熱效率、等熵效率等特性參數(shù)。綜合各轉(zhuǎn)速下的風扇定常特性,得到風扇壓比特性曲線如圖18所示。其中設計最大轉(zhuǎn)速為7975 r/min。
圖18 風扇總壓比特性曲線
由前2節(jié)所得進氣道與風扇單獨的工作特性,根據(jù)換算流量的定義
對于上節(jié)的風扇工作特性,由于設置的進口條件與標準大氣一致,可以直接將結(jié)果作為換算流量和換算轉(zhuǎn)速使用。換算流量匹配如圖19所示,其中3條與x軸平行的線分別代表3個進氣道換算流量:128.44、110.94、76.87 kg/s。
由此選擇斜板角度為2°,Ma=1.33的進氣道與設計最大轉(zhuǎn)速為7975 r/min的風扇,理論上流量匹配。應用分裂激盤模型對該狀態(tài)下的風扇葉柵進行模化,進氣道-風扇聯(lián)合流場的計算域如圖20所示,真實的網(wǎng)格劃分至動葉的前緣之前,動葉和靜葉之間間隙和靜葉后的網(wǎng)格均由程序自動生成,如圖21所示。
圖19 換算流量匹配
圖20 聯(lián)合流場的拓撲結(jié)構(gòu)
圖21 風扇部分自動生成的網(wǎng)格結(jié)構(gòu)
對進氣道-壓氣機聯(lián)合流場進行數(shù)值模擬,可得聯(lián)合內(nèi)流特性,如圖22所示。從圖中可見,聯(lián)合流場的最大流量要略小于單獨進氣道的內(nèi)流最大流量。分析可知主要原因是進氣道出口的總壓畸變影響了風扇的性能,導致風扇能夠通過的最大流量受到限制。
對相同邊界條件下的單獨進氣道與聯(lián)合流場數(shù)值模擬結(jié)果對比,內(nèi)流馬赫數(shù)與進氣道出口的總壓分別如圖23、24所示。
圖22 發(fā)動機流量與風扇出口反壓的關系
圖23 馬赫數(shù)對比
圖24 進氣道出口總壓對比
由單獨進氣道與聯(lián)合流場的結(jié)果對比可知,加入分裂激盤模型進行聯(lián)合計算對進氣道內(nèi)流特性的影響可以忽略。由于理論上已經(jīng)證明進氣道與風扇達到匹配狀態(tài),聯(lián)合流場與單獨的數(shù)值模擬結(jié)果的一致性證明了分裂激盤模型可對進氣道/風扇的聯(lián)合流場進行一體化分析。
本文基于激盤模型建立了1個相對更準確的分裂激盤模型,減小了激盤模型假設產(chǎn)生的誤差。通過對Rotor 67轉(zhuǎn)子的數(shù)值模擬驗證,分析多級壓氣機在穩(wěn)態(tài)總壓畸變的性能影響,得到以下結(jié)論:
(1)分裂激盤模型對風扇的?;椒ㄊ蛊洳恍枰鄳男拚瘸R?guī)激盤模型更精確,同時計算收斂速度沒有明顯減慢。且不同于常規(guī)激盤模型激盤前、后流場參數(shù)無法使用的問題,該模型由于沒有對風扇幾何進行?;?,因此激盤前、后的流場參數(shù)有一定的參考價值。
(2)分別計算進氣道和風扇的工作特性,通過折合流量這一參數(shù),得到進氣道與風扇流量匹配的理論狀態(tài)點。通過分裂激盤模型對風扇?;?,并與進氣道進行聯(lián)合計算,得到的流場參數(shù)驗證了使用分裂激盤模型對進發(fā)一體化問題的可行性。
綜上所述,該模型在原激盤模型的基礎上提高了計算精度,更準確地描述了葉柵級之間的流動規(guī)律,同時有效地對針對進氣道-風扇的聯(lián)合流場進行性能分析。相對于其他常見的簡化模型如黏性徹體力模型,分裂激盤模型可對進氣道-壓氣機流場進行聯(lián)合數(shù)值計算,使其成為分析進發(fā)一體化問題的有力工具。