解亞東，朱 賢，丁建國

（中航工業(yè)商用航空發(fā)動機有限責(zé)任公司，上海201108）

大涵道比風(fēng)扇轉(zhuǎn)子的非定常數(shù)值模擬

解亞東，朱 賢，丁建國

（中航工業(yè)商用航空發(fā)動機有限責(zé)任公司，上海201108）

為了更好地設(shè)計大涵道比風(fēng)扇轉(zhuǎn)子及外涵導(dǎo)葉(O G V)，采用非線性諧波法對某大涵道比風(fēng)扇轉(zhuǎn)子及O G V進行了非定常和定常數(shù)值計算，2種計算方法獲得的風(fēng)扇外涵級的流量、增壓比、絕熱效率有一定差異，并分析了性能變化的原因。通過對其典型截面的非定常流場進行快速傅里葉分解，在頻域內(nèi)分析了風(fēng)扇轉(zhuǎn)子和O G V之間的轉(zhuǎn)靜干擾效應(yīng)，以及交界面典型脈動量的周向波形，從頻域的角度分析了影響性能的原因；同時觀察到一些時域中難呈現(xiàn)的現(xiàn)象。結(jié)果表明：通過非定常數(shù)值計算，大涵道比風(fēng)扇外涵級絕熱效率曲線隨時間變化呈現(xiàn)正弦曲線的形狀，幅值約為0.59個百分點；風(fēng)扇外涵級絕熱效率最高、最低點分別出現(xiàn)在風(fēng)扇尾跡通過靜子通道的約3/4和1/4處。

大涵道比風(fēng)扇轉(zhuǎn)子；非線性諧波法；非定常計算；時域分析；頻域分析

0 引言

國內(nèi)外風(fēng)扇數(shù)值模擬普遍采用3維定常和非定常模擬技術(shù)。定常模擬比非定常模擬的計算量小，廣泛應(yīng)用于工程設(shè)計，但是在處理多葉片排交界面時，采用摻混面平均的方式而忽略了上下游的非定常相互作用。非定常模擬相對完整地保留了其技術(shù)細節(jié)，但是計算量龐大，難以用于工程設(shè)計。1985年，美國NASA劉易斯研究中心的Adamczyk[1]引入了3個平均算子，建立了描述葉輪機內(nèi)確定性非定常流動的通道平均方程(Average Passage Equations)，目的是通過引入確定性應(yīng)力來反映轉(zhuǎn)靜干涉對時均流動的影響。正如雷諾平均帶來的“雷諾應(yīng)力”模化難題一樣，Adamczyk的通道平均方法引入了“確定應(yīng)力”的建模難題。為此，20多年來，國內(nèi)外學(xué)者進行了許多創(chuàng)新性研究[2-6]。1992年，Giles提出了漸進求解的思路[7]，其后，在其啟發(fā)下進行了創(chuàng)造性發(fā)展[8-9]，開發(fā)了非線性諧波法(frequency based nonlinear harmonic method，以下簡稱NHM)。該方法將非定常流場分成時均和脈動流場2部分進行耦合求解，分別通過求解含有決定性應(yīng)力的時均流動方程和在頻域內(nèi)N階諧波得到。在確定應(yīng)力建模方法方面較嚴(yán)謹、精度較高，不需要經(jīng)驗參數(shù)，應(yīng)用范圍較廣。在工程應(yīng)用中[10-12]，一般只需要前2～3階諧波[13-14]就可以較好地模擬確定性應(yīng)力。

本文以某大涵道比風(fēng)扇轉(zhuǎn)子及外涵出口導(dǎo)葉（OGV）為第1級，進行了定常和非定常計算，并在時域?qū)用嫔戏治隽思夐g性能的差異。通過對典型截面非定常流場的FFT分解得到各物理量的頻譜圖域，分析了轉(zhuǎn)靜干涉現(xiàn)象和典型截面脈動量周向波形。

1 數(shù)值方法

1.1 定常和非定常計算

采用商用軟件NUMECA對某風(fēng)扇轉(zhuǎn)子及OGV進行了定常和非定常計算。其中風(fēng)扇和OGV的葉片數(shù)目分別為18、55。

為了減小網(wǎng)格對定常和非定常計算結(jié)果的影響，均采用相同的網(wǎng)格數(shù)目及拓撲結(jié)構(gòu)，湍流模型一致，且內(nèi)外涵道均采用靜壓邊界條件，如圖1所示。

圖1 3維網(wǎng)格

1.2 快速傅里葉分解

在非定常計算結(jié)果中，空間每點物理量非定常時變曲線包含64個點，延拓為8個周期后進行快速傅里葉變換（FFT）計算，來提高采樣精度，并保存前3階諧波的實部和虛部。

OGV進口處某一點軸向速度的時域和頻域曲線分別如圖2、3所示。從圖3中可見，在0 Hz處的幅值（軸向速度的均值）最高。OGV軸向速度頻域各基頻和倍頻曲線如圖4所示，是前3階頻率下軸向速度脈動值的放大顯示。其中基頻及各倍頻以外的頻率對應(yīng)的幅值很小，可見時域曲線的周期性良好。

圖2 OGV軸向速度的時域曲線

圖3 OGV軸向速度的頻域曲線

圖4 OGV軸向速度頻域各基頻和倍頻曲線

2 時域分析

2.1 非定常流場時域流場的分析

風(fēng)扇轉(zhuǎn)子和OGV組合為外涵的第1級，其非定常計算的級絕熱效率分布如圖5所示，OGV非定常計算總壓恢復(fù)系數(shù)分布如圖6所示。從圖中可見，級絕熱效率分布類似于正弦曲線，幅值約為0.59個百分點，與試驗結(jié)果[15]吻合。級最高、最低效率點分別出現(xiàn)在轉(zhuǎn)子轉(zhuǎn)過0.765倍、1/4的OGV通道周向弧長處。風(fēng)扇轉(zhuǎn)子非定常計算絕熱效率分布如圖7所示。從圖中可見，非定常計算對風(fēng)扇轉(zhuǎn)子的效率影響較小，由于風(fēng)扇轉(zhuǎn)子和OGV的軸向長度是1個風(fēng)扇轉(zhuǎn)子中徑的弦長，尾跡的影響僅僅往下游傳播，勢流影響在向上游傳播時經(jīng)過較長距離，耗散較多，對風(fēng)扇轉(zhuǎn)子效率的影響較小。這也可從后面的頻域分析得知。

圖5 風(fēng)扇轉(zhuǎn)子和OGV非定常級絕熱效率分布

圖6 OGV非定常計算總壓恢復(fù)系數(shù)分布

圖7 風(fēng)扇轉(zhuǎn)子非定常計算絕熱效率分布

由于OGV受風(fēng)扇轉(zhuǎn)子尾跡和勢流影響較大，會出現(xiàn)總壓恢復(fù)系數(shù)類似于級效率分布曲線。OGV吸力面相對靜壓分布如圖8所示。從圖中可見，Ps、Paver作為無量綱物理量來顯示細微差別，分別指靜壓和轉(zhuǎn)靜交界面處的平均靜壓。在吸力面葉中位置，最大效率比最小效率處的最大靜壓降略微靠近前緣。OGV壓力面相對靜壓分布如圖9所示。從圖中可見，在70%葉高位置，最大比最小效率處的高靜壓分布范圍更廣。說明了尾跡在掃掠OGV通道的不同位置會帶來不同的靜壓分布，這也是導(dǎo)致效率高低的原因。

圖8 OGV吸力面相對靜壓分布

圖9 OGV壓力面相對靜壓分布

在設(shè)計點處且級壓比不變的情況下，非定常計算外涵流量和風(fēng)扇轉(zhuǎn)子絕熱效率比定常計算的分別增加0.00401%、0.07%。非定常計算OGV總壓恢復(fù)系數(shù)比定常計算的減小0.83%，級絕熱效率提高0.05%。2種計算結(jié)果非常近似。這主要是由于風(fēng)扇轉(zhuǎn)子和OGV軸向距離太遠，尾跡和勢流對性能的影響較小。

3 頻域分析

3.1 頻域流場的分析

葉中截面軸向速度第1階諧波實部如圖10所示。從圖中可見，OGV圖像呈條帶形，走向，與尾跡方向相同。尾跡在掃掠過的區(qū)域的特點是以時間為橫坐標(biāo)的流場參數(shù)的時變曲線初始相位角近似，N階諧波的實部/虛部是幅值和該角度的正弦函數(shù)/余弦函數(shù)的乘積，所以在該角度下進行FFT分解會保證實部（或虛部）在尾跡掃掠過區(qū)域的取值近似相等。條文呈正負交替也是由于正、余弦函數(shù)正、負的特點所致。

圖10 葉中截面軸向速度1階諧波實部

通過非定常的參數(shù)計算得出，OGV計算域內(nèi)包含3～4條條帶。OGV通道長度為0.5048 m，軸線速度約為180 m/s，掃掠周期為8.6557×10-4s，通過時間為2.802×10-3s，OGV區(qū)域內(nèi)尾跡條數(shù)為后者與前者的比值3.239，與頻域內(nèi)的條代數(shù)基本吻合。

近輪轂截面和機匣截面周向速度的第1階諧波實部分布分別如圖11、12所示。從圖中可見，條帶分布規(guī)律大體與葉中的類似。由于上游轉(zhuǎn)子出口的氣流角和軸向速度周向分布的不同導(dǎo)致條帶走向和幅值存在差異。其中在葉根處軸向速度的幅值要高于葉尖處的，中間截面最大，即在葉尖處的尾跡影響比葉中和葉根都偏小，所以需要改善風(fēng)扇轉(zhuǎn)子的葉中和葉根。以下為葉中典型界面的分析。

圖11 5%葉展截面周向速度第1階諧波實部

圖12 95%葉展截面周向速度第1階諧波實部

軸向速度第1階諧波的幅值分布如圖13所示。由于交界面是OGV中非定常擾動的源頭，非定常尾跡干擾通過此處向下游傳播，所以從圖中可見在交界面處幅值最大。同時幅值延軸向衰減，如圖14所示。

圖13 軸向速度第1階諧波幅值

圖14 軸向速度第1階幅值沿軸向衰減

軸向速度第 2階諧波實部如圖15所示。從圖中可見，頻率增加，對應(yīng)的波長縮短，導(dǎo)致條紋比第1階諧波密。軸向速度第3階諧波實部如圖16所示。從圖中可見，條紋變得更密。

圖15 軸向速度第2階諧波實部

圖16 軸向速度第3階諧波實部

綜上所述，諧波頻率越高，實部幅值越小且軸向衰減速度越快。

第1～3階諧波實部進口處的周向波形如圖17所示。從圖中可見，隨著諧波階數(shù)的增加，所包含的周期性也增加。

圖17 OGV進口周向速度前3階諧波實部周向波形

切向、徑向速度和密度等因素受尾跡影響，其頻譜圖譜與軸向速度類似，這里暫不討論。靜壓的分部則明顯不同。葉中截面靜壓第1階諧波實部如圖18所示。從圖中可見，由于靜壓受尾跡影響小，更多受葉片勢流干擾。

圖18 葉中截面靜壓第1階諧波實部

上文針對OGV通道中的分部進行分析，下文進行風(fēng)扇轉(zhuǎn)子通道中軸向速度分析。

葉中截面軸向速度第1階諧波實部如圖19所示。從圖中可見，風(fēng)扇轉(zhuǎn)子通道中的第1階諧波也呈條帶狀。軸向速度的第1階諧波幅值也是在交界面處最大，但是幅值很小，表明下游OGV對上游風(fēng)扇轉(zhuǎn)子的勢干擾較弱，幅值在風(fēng)扇轉(zhuǎn)子通道中衰減更快，如圖20所示。

圖19 葉中截面軸向速度第1階諧波實部

圖20 葉中截面軸向速度第1階諧波幅值

4 結(jié)論

（1）在設(shè)計點處，采用3階非線性諧波法進行了非定常計算；與定常計算相比，大涵道比風(fēng)扇轉(zhuǎn)子和OGV的流量、壓比、絕熱效率有一定差異，但由于其軸向間距較大，其差異不大。

（2）通過非定常計算可以得出大涵道比風(fēng)扇轉(zhuǎn)子和OGV級絕熱效率曲線隨時間變化呈現(xiàn)正弦曲線，幅值約為0.59個百分點。

（3）通過非定常計算可以得出風(fēng)扇轉(zhuǎn)子和OGV級絕熱效率最高、最低點分別出現(xiàn)在風(fēng)扇尾跡通過靜子通道約3/4和1/4處。

綜上所述，通過對頻域圖譜的研究，從新的視角分析了風(fēng)扇外涵道級非定常流場。

[1]Adamczyk J J.Model equation for simulation flows in multi stage turbomachinery[R].ASME 85-GT-226.

[2]Adamczyk J J.A model for closing the inviscid form of the average passage equation system[J].Journal of Turbinemachinery，1986，108（2）：180-186.

[3]Rhie C M，Gleixner A J，Spear D A，et al.Development and application of a multistage Navier-Stokes solver-PartⅠ：multistage modeling using body-forces and deterministic stresses[R].ASME 95-GT-342.

[4]Sondak D L，Dorney D J，Davis R L.Modeling turbomachinery unsteadiness with lumped deterministic stresses[R].AIAA-1996-2570.

[5]AGVD Wall，JR Kadambi，JJ Adamczyk.A transport model for the deterministic stresses associated with turbomachinery blade row interaction[J].Journal of Turbomachinery，2000，122（4）：593-603.

[6]季路成.軸流式葉輪機轉(zhuǎn)子/靜子干擾非定常流動探索研究[D].北京：北京航空航天大學(xué)，1998. JI Lucheng.Research on rotor/stator interaction in axial turbomachinery [D].Beijing：Beihang University，1998.（in Chinese）

[7]Giles M.An approach for multistage calculations incorporating unseadiness[R].ASME 92-GT-282.

[8]He L，Ning W.Efficient approach for analysis of unsteady viscous flows in turbomachines[J].AIAA Journal，1998，36（11）：2005-2012.

[9]Chen T，Vasanthakumar P，He L.Analysis of unsteady blade row interaction using nonlinear harmonic approach[J].Journal of Propulsion and Power，2001，17（3）：651-658.

[10]Martens S，Shin H，Gliebe P.Rotor wake unsteady flow field hot-wire measurements in the universal propulsion simulator(UPS)at NASA[R].GE Aircraft Engines TM 97-67.

[11]GanzU W.Experimental investigation of the unsteady flow characteristics in the Boeing 18 inch fan rig[R].AIAA-99-1886.

[12]Verhoff V G.Three dimensional laserwindow formation[R]. NASA-RP-1992-1280.

[13]趙軍，柳陽威，劉寶杰.ASA67級非定常流場的頻域分析[J].航空動力學(xué)報，2007，22（8）：1371-1377. ZHAO Jun，LIU Yangwei，LIU Baojie.Frequency-domain analysis of NASA67 unsteady flow field[J].Journal of Aerospace Power，2007，22（8）：1371-1377.（in Chinese）

[14]李翔.1.5級渦輪非定常流動研究 [D].哈爾濱：哈爾濱工程大學(xué)，2010. LI Xiang.Unsteady flow investigation of a 1.5-stage axial turbine[D]. Harbin：Harbin Engineering University，2010.（in Chinese）

[15]St?ding J，F(xiàn)riedrichs J，Waitz T，et al.The potential of rotor and stator clocking in a 2.5-stage[C]//ASME TURBO EXPO 2012. Copenhagen，Denmark：GT2012-68353，2012.

（編輯：肖磊）

Unsteady Numerical Simulation of High Bypass Ratio Fan

XIE Ya-dong，ZHU Xian，DING Jian-guo
（AVIC Commercial Aircraft Engine Co.Ltd.，Shanghai 201108，China）

The unsteady and steady numerical calculation of a high bypass ratio fan and OGV were obtained by non linear harmonic method to design them better.The bypass flows pressure ratio,stage adiabatic efficiency of the fan which got from the two simulation methods were different,and the reasons were analyzed.Based on the Fast Fourier Transformation (FFT)of the typical profiles of the stage, the rotor-stator interactions between the fan and OGV were analyzed,and the circumferential distribution of typical unsteady perturbation was studied to give out the influences reason in the frequency domain.The phenomena hardly observed in the time domain can be discovered simultaneously.The results show that the stage adiabatic efficiency of the high bypass ratio fan is a sinusolid according to the time,and the amplitude value is about 0.59 percent according to unsteady numerical calculation.The maximum and minimum value of the stage adiabatic efficiency is located on 3/4 and 1/4 OGV that the fan wake passed through.

high bypass ratio fan； non linear harmonic method；unsteady calculation；time domain analysis；frequency domain analysis

V 211.45

A

10.13477/j.cnki.aeroengine.2015.05.004

2014-07-24

解亞東（1987），男，碩士，從事風(fēng)扇/增壓級氣動性能設(shè)計工作；E-mail：xieyadong163@163.com。

解亞東，朱賢，丁建國.大涵道比風(fēng)扇轉(zhuǎn)子的非定常數(shù)值模擬[J].航空發(fā)動機，2015，41（4）：20-23.XIE Yadong，ZHUXian，DING Jianguo.Unsteady numerical simulation ofhigh bypass ratiofan [J].Aeroengine，2015，41（4）：20- 23.