蔣 聰，劉 冕，熊 欣

(中國航空工業(yè)集團公司金城南京機電液壓工程研究中心，江蘇 南京 211106)

進氣畸變是導致航空發(fā)動機氣動穩(wěn)定性惡化的重要因素之一，亟需探索能高效且準確預判壓氣機在畸變流場下的性能和穩(wěn)定性的計算方法[1]?；兞鲌鐾ǔＪ遣粚ΨQ的，因此計算域需要以壓氣機整環(huán)為目標，如采用三維非定常的雷諾平均Navier-Stokes方程(以下簡稱RANS)[2-3]，但當前的硬件水平難以承受其數(shù)量宏大的網格，所需耗費的時間成本對日常工程設計工作而言十分可觀。前人為解決RANS計算量大的情況，提出了徹體力模型[4-6]。

徹體力模型是根據無窮多葉片理論推定的，不再拘泥于葉片的具體結構和形狀，將葉片排對流體產生的壓升及轉向影響簡化為分布式徹體力，并將其作為源項代入到三維歐拉(Euler)方程中針對葉片通道內的畸變流場進行求解，很好地表征了其大尺度特征。

徹體力模型無需直接針對葉片力[7]和黏性細節(jié)進行耗時較大的計算，大大縮短了計算時長，因此將其應用于工程分析具有較大優(yōu)勢。與以前的二維計算模型不同[8]，徹體力模型直接針對三維流動控制方程進行求解，可以更準確地將畸變流場的三維特征反映出來[9]。

本文根據徹體力模型的基本理念，編制出一款可以預測壓氣機在進氣畸變條件下氣動特性及穩(wěn)定性變化的數(shù)值計算程序，使用該程序研究了周向畸變流場中大涵道比風扇轉子特性，并得出畸變流場的三維特征。

1 計算程序建模

1.1 控制方程

本文使用的控制方程是絕對直角坐標系下具有源項的三維非定常Euler方程：

(1)

式中：t為時間；ρ為氣流密度；e為氣流總能量；Vx,Vy,Vz分別為沿x,y,z軸的速度；x,y,z分別為3個維度坐標軸；P為氣流壓力；Fx,Fy,Fz分別為沿x,y,z軸的徹體力；Lu為輪緣功。

1.2 源項構造

徹體力的構造選擇在絕對圓柱坐標系下進行，然后依據圖1所示的坐標關系演變?yōu)橹苯亲鴺讼迪碌男问健?/p>

圖1 徹體力演變

徹體力源項在葉片區(qū)域分解為垂直于相對速度V′的非耗散項φ和平行于相對速度V′的耗散項f，得到式(2)～式(4)：

(2)

(3)

(4)

(5)

式中：α為如圖2所示角度；Vm為子午速度。

圖2 子午面示意圖

在式(5)的基礎上，最終得到式(6)：

(6)

假定子午面網格線為流線，則α可以由網格幾何參數(shù)求得。如此只要求出式(3)中的|f|以及式(6)中的φθ，再結合三維Euler求解器的瞬時速度場即可求得徹體力源項。

基于軸對稱假設及熱力學第一、第二定律，Marble提供如下解法求解|f|及φθ：

(7)

(8)

式中：r為半徑；m為子午坐標；T為氣流靜溫；s為氣流熵。式(7)中?rVθ/?m表示葉片區(qū)域速度環(huán)量沿子午面流線的導數(shù)，離散為Δ(rVθ)/Δm，Δ(rVθ)及Δm本質上應該是葉片區(qū)子午面每一網格單元的參數(shù)。如果假設葉片區(qū)任一葉高進出口的速度環(huán)量變化量Δ(rVθ)沿葉片進出口的子午面流線線性分布，則Δ(rVθ)及Δm就可以表示為葉片區(qū)任一葉高進出口的參數(shù)，而式(7)中的剩余變量ρ,Vm,r仍為網格單元的參數(shù)。對式(8)中?s/?m作相同處理。

根據已知的葉排特性求得葉片進出口的Δ(rVθ)及Δs(氣流熵的變化量)，并結合三維Euler求解器的瞬時流場，就可以利用式(3)及式(8)求得fr，fθ，fz，再由式(7)求得φθ，繼而求出φr，φz，能量方程中的源項Lu=Fθωr，至此得到葉片區(qū)域徹體力源項的分布。

Δ(rVθ)及Δs的計算式如下：

Δ(rVθ)=r2Vθ2-r1Vθ1=r2Vz2tanβ2-r1Vθ1

(9)

(10)

而對于轉子部件總溫比τ*與出氣角β2，它們之間有如下關系：

r2Vz2tanβ2-r1Vθ1=U2Vθ2-U1Vθ1=

(11)

1.3 特性關聯(lián)

2 總壓畸變條件下風扇流場的數(shù)值模擬

2.1 計算對象

本文以某大涵道比渦扇發(fā)動機的進口風扇轉子為算例開展計算與分析，部分參數(shù)見表1。

表1 進口風扇轉子部分參數(shù)

本文計算程序根據流場特性模擬劃分的網格如圖3所示，周向、徑向和軸向的網格單元數(shù)分別為60,20和79。沿壓氣機全環(huán)通道生成網格(不含椎體)，并以其作為計算域，模擬進氣畸變對壓氣機的影響。由于對于徹體力模型而言，網格密度對計算結果精確性的影響不大[10-11]，因此本文所劃分的網格密度是完全可行的。葉片的具體形狀在徹體力模型計算時無需考慮，因此可以靠兩條徑向網格線來確定葉片區(qū)域的前后緣位置，且周向網格單元對于旋轉軸線是對稱的。由此可知，計算網格簡單且均勻，大大加快了收斂和數(shù)值計算的速度。

圖3 網格單元

2.2 計算程序可靠性對比驗證

為驗證本文程序的計算準確性，將其計算所得數(shù)據與NUMECA計算結果進行對比。特性計算針對均勻進氣條件下的風扇開展，在不同背壓條件下得出包含總溫和總壓隨折合流量變化的流場特性，對比情況如圖4所示。

圖4 特性對比

從曲線來看，本程序與NUMECA計算得到的風扇特性結果基本相符，其中最大質量流量點、總溫比-流量特性和總壓比-流量特性與NUMECA計算結果基本吻合，特性線的趨勢也基本一致。

根據以上對比分析結果，可以認為本計算程序采用的計算模型能夠對風扇轉子特性曲線進行精度較高的預估，并且可以可靠地描繪風扇的相關流場。

2.3 穩(wěn)態(tài)進氣條件下畸變流場模擬計算

采用2.1節(jié)給出的計算對象，并按表2數(shù)據設定畸變條件。

表2 畸變條件

進氣總壓畸變條件下畸變特性點在均勻進氣特性圖上的位置如圖5所示。

圖5 進氣畸變特性點

在風扇進口設定范圍為90°的畸變區(qū)，如圖6所示。

圖6 進口總壓畸變區(qū)設定

經過本文計算程序計算后得出的總壓畸變如圖7所示，可以清晰地看出畸變流動的強三維特征，且低壓區(qū)在葉根和葉尖處有明顯的相位移動?？倝旱牟痪栽诮涍^風扇轉子后被削弱了，這是由于低壓區(qū)經過轉子后會產生較高的壓升。在葉片尾緣截面上，轉子葉片旋出低壓畸變區(qū)的區(qū)域附近出現(xiàn)高壓區(qū)，沿著葉片旋轉方向該高壓區(qū)逐漸衰減，最低數(shù)值出現(xiàn)在轉子葉片旋入畸變區(qū)時。

圖7 總壓畸變分布

從圖8可知，風扇進口處的總壓畸變將在轉子葉片尾緣處誘發(fā)產生總溫畸變，且在畸變區(qū)總溫最高，產生這種現(xiàn)象的原因是由于在轉子進口的周向方向上總壓參數(shù)不均勻，繼而導致不同轉子葉片通道內的加功量不同，低總壓區(qū)域的加功量大且表現(xiàn)在葉片出口就是總溫較高。

圖8 葉片尾緣截面總溫畸變分布

3 結束語

本文給出的基于三維徹體力的計算模型具有較高的預估精度，能夠較好地刻畫出畸變流場流動的大尺度和三維特征，為預測壓氣機在進氣畸變條件下的氣動特性提供了一種新的解決方法。但是，目前的研究工作僅通過計算實例來驗證模型的準確性，要真正在實踐中體現(xiàn)其價值，還需要通過試驗結果對比校正，才能對其進行進一步的完善。