壓力畸變試驗進氣道軸向力數(shù)值計算研究

吳志昌，徐國輝，劉 毅，魏利永，高雷兵

(中國航發(fā)沈陽發(fā)動機設(shè)計研究所，沈陽 110015)

1 引言

當(dāng)今世界對戰(zhàn)斗機隱身性及機動性的要求越來越高，隨之帶來的問題是發(fā)動機進口流場品質(zhì)變差，這不僅會降低發(fā)動機性能，嚴重時甚至?xí)?dǎo)致其進入不穩(wěn)定工作狀態(tài)，為此對發(fā)動機的抗畸變能力提出了更高的要求[1-4]。

為了評估航空渦輪風(fēng)扇發(fā)動機的穩(wěn)定性，國內(nèi)外研究者對壓力畸變及相關(guān)試驗進行了大量研究。如Braithwaite 等[5]基于J85-GE-13 發(fā)動機研究了進口180°范圍內(nèi)的總溫總壓畸變對發(fā)動機總壓比及損失的影響，并由此發(fā)展了一個新的計算模型；Mehalic 等[6]基于渦扇發(fā)動機，研究了進口總溫總壓組合畸變對發(fā)動機氣動穩(wěn)定性的影響；Beale 等[7]在發(fā)動機地面試車臺模擬了多種飛行狀態(tài)下進氣畸變的狀態(tài)，并提出了相關(guān)優(yōu)化方法；朱愛迪[2]用CFD 手段模擬插板及板后流場，并在此基礎(chǔ)上研究了模擬板的設(shè)計技術(shù)；丁寧[3]對大涵道比渦扇發(fā)動機進氣總壓畸變進行了多種狀態(tài)的模擬研究，并對其穩(wěn)定性進行了計算分析；李亮等[8]為弄清國軍標(biāo)GJB/Z 64A-2004(K)采用簡單的單個弦月形插板作為進口總壓畸變評定實驗用標(biāo)準插板的緣由，在雙級低速軸流壓氣機試驗器上開展了4 種不同形式插板的進氣總壓畸變實驗驗證研究，并詳細比較了各插板的綜合畸變指數(shù)特性；馬明明等[9]在吊艙進氣道進口安裝擾流板進行試驗，研究了擾流板進氣畸變的影響因素及總壓畸變特征，得到了進氣道出口若干馬赫數(shù)下進氣總壓畸變的定量數(shù)據(jù)；程邦勤等[10]分析了進氣總壓畸變對某型渦扇發(fā)動機氣動穩(wěn)定性的影響，得到了穩(wěn)態(tài)總壓場不均勻度與總壓紊流度隨插板深度及發(fā)動機轉(zhuǎn)速的變化關(guān)系，確定了該型發(fā)動機的最大綜合畸變指數(shù)和總壓畸變敏感系數(shù)；劉作宏等[11]分別采用平插板和角插板方式，在插板式壓力畸變發(fā)生器上開展了某型發(fā)動機進氣畸變試驗，并采用綜合指數(shù)法對數(shù)據(jù)進行分析；彭生紅等[12]基于中等涵道比渦扇發(fā)動機結(jié)構(gòu)及其氣動穩(wěn)定性特點，通過該類型某發(fā)動機高空壓力畸變試驗開展了相關(guān)研究，獲取了研究對象合理可行的高空壓力畸變試驗方案設(shè)計、畸變裝置選取以及數(shù)據(jù)處理規(guī)范。但在公開文獻中，未見針對風(fēng)險系數(shù)較高的發(fā)動機總壓畸變用設(shè)備的安全設(shè)計。

壓力畸變試驗中，帶有插板的進氣道在給發(fā)動機提供進氣壓力畸變的同時，也使得自身受到較大的氣動載荷[13]。為評估插板及固定進氣道臺架設(shè)備的安全系數(shù)，保證發(fā)動機進氣壓力畸變試驗順利進行，有必要計算其所承受的軸向載荷。本文以帶有插板裝置的喇叭口型進氣道[14]為研究對象，利用Fluent 軟件對插板位于進氣道2 個典型位置流場進行數(shù)值模擬，然后對數(shù)值求解得到的云圖進行分析，最后利用積分的方法分別求得2 個典型位置時進氣道唇口和插板所受沿軸向的合力。

2 數(shù)值模型及網(wǎng)格劃分

圖1 為設(shè)計的進氣道型面，除插板插入深度不同外，插板軸向位置(距離型面曲線切點2D，D為進氣道直徑)以及進氣道尺寸等其他幾何條件相同。

圖1 進氣道型面尺寸Fig.1 Air intake dimensions

圖2 為數(shù)值計算的三維模型，插板A 位置較插板B 位置插入進氣道內(nèi)的深度較淺。

圖2 數(shù)值計算三維模型Fig.2 Numerical 3D model

圖3 為整體數(shù)值計算的網(wǎng)格模型。對模型網(wǎng)格局部進行加密，保證能夠模擬比較復(fù)雜的流場細節(jié)。為保證對壁面附近流動的準確模擬，進行了邊界層網(wǎng)格設(shè)置，靠近壁面第一層網(wǎng)格高度為0.01 mm。經(jīng)計算查看，Y+大小約為100，滿足所用湍流模型30～300 的計算要求。

圖3 數(shù)值計算網(wǎng)格模型Fig.3 Numerical gird model

取進氣道內(nèi)某點壓力作為網(wǎng)格無關(guān)性驗證參考點，如圖4 所示經(jīng)過網(wǎng)格無關(guān)性驗證(后3 次計算誤差約為0.04%)。對于插板位于進氣道A 位置時的計算模型，最終得到總的網(wǎng)格節(jié)點數(shù)為797 446，總的網(wǎng)格單元數(shù)為4 537 270；對于插板位于B 位置時的計算模型，總的網(wǎng)格節(jié)點數(shù)為827 822，總的網(wǎng)格單元數(shù)為4 709 033。

圖4 網(wǎng)格驗證Fig.4 Mesh validation

3 數(shù)值計算方法

3.1 湍流模型

本次數(shù)值求解的湍流模型選擇雙方程模型，即標(biāo)準k-ε模型。此模型是在單方程模型的基礎(chǔ)上，再引入了一個關(guān)于湍動能耗散率ε的方程。標(biāo)準k-ε模型適用于初始迭代、設(shè)計選型和參數(shù)研究，并對近壁面選用標(biāo)準壁面函數(shù)進行處理[15]。

3.2 邊界條件

數(shù)值模擬時，需要對計算域的進出口和壁面等邊界進行合理的設(shè)置。插板位于A 位置和B 位置時試車間進口均設(shè)置為總溫總壓進口，試車間出口為靜壓出口，進氣道出口為靜壓出口(對應(yīng)發(fā)動機某一流量狀態(tài))；其余設(shè)置為無滑移壁面邊界。

3.3 求解方法

壓力-速度耦合采用Simple 算法，空間離散中梯度求解選擇基于最小二乘單元算法，壓力求解選擇Standard(標(biāo)準)形式，動量、能量、湍動能和湍流擴散率的求解選擇二階迎風(fēng)格式。

4 插板位于進氣道A 位置時計算結(jié)果

4.1 數(shù)值計算結(jié)果

圖5 示出了插板位于進氣道A 位置時數(shù)值計算得到的壓力分布云圖?？梢钥吹?，進氣道唇口部分逆航向視圖的壓力要小于順航向視圖的壓力，說明唇口部分進氣道受到的軸向力的合力方向應(yīng)該為發(fā)動機順航向；而插板部分與此正好相反，插板部分受到的合力方向應(yīng)該為發(fā)動機逆航向。

圖5 插板位于進氣道A 位置時的壓力分布云圖Fig.5 Pressure cloud of position A

4.2 插板不同半徑位置前后表面壓力分布

圖6 為插板位于進氣道A 位置時的幾何示意圖。圖中，R為插板表面某一位置距離進氣道中心軸線的距離，其變化范圍為148.50～450.00 mm。通過積分的方法求解作用在插板表面的壓力。

圖6 插板位于進氣道A 位置示意圖Fig.6 Schematic of position A

表1 為插板不同半徑位置前后表面壓力分布。表中，p1、p2分別表示插板不同半徑位置前后表面壓力，Δp表示其差值。p1、p2大小通過前面的數(shù)值計算得到。

表1 插板不同半徑位置前后表面壓力分布(A 位置)Table 1 Pressure distribution at different radius of flashboard (location A)

圖7 給出了插板前后表面壓差Δp與位置R的多項式關(guān)系曲線，該關(guān)系式通過圖中曲線擬合得到。

圖7 插板前后表面壓差與半徑的關(guān)系曲線(插板A 位置)Fig.7 Pressure difference vs.radius of flashboard (location A)

插板所受軸向力的計算公式為：

其中：R1=148.50 mm，R2=450.00 mm，p(R)由圖7中擬合曲線的四次多項式方程給出。

由公式(1)計算得到的插板所受軸向力為7 481 N。

4.3 進氣道唇口前后表面壓力分布

表2 示出了進氣道唇口不同半徑位置處的前后表面壓力分布。其中，R′為進氣道唇口表面某一位置距離進氣道中心軸線的距離，其變化范圍為465.00～805.00 mm；、分別表示進氣道唇口不同位置前后表面壓力，Δp′表示其差值。、大小通過前面的數(shù)值計算得到。

表2 進氣道唇口前后表面壓力分布(A 位置)Table 2 Pressure distribution at different radius of intake (location A)

圖8 為進氣道唇口前后表面壓差Δp′與位置R′的多項式關(guān)系曲線。關(guān)系式通過曲線擬合得到。

圖8 進氣道唇口前后表面壓差與半徑的關(guān)系曲線(A 位置)Fig.8 Pressure difference vs.radius of intake (location A)

進氣道唇口所受軸向力的計算公式為：

其中：=465.00 mm，=805.00 mm，p(R′)由圖8中的曲線擬合得到的多項式方程給出。由公式(2)計算得到的進氣道唇口所受的軸向力為7 374 N。

根據(jù)前面計算結(jié)果，插板位于進氣道A 位置時，進氣道唇口和插板部分產(chǎn)生的軸向力的合力大小為107 N，方向為發(fā)動機逆航向。

5 插板位于進氣道B 位置時計算結(jié)果

5.1 數(shù)值計算結(jié)果

圖9 為插板位于進氣道B 位置時數(shù)值計算得到的壓力云圖?？梢钥吹?，進氣道唇口部分逆航向視圖的壓力要小于順航向視圖的壓力，說明唇口部分進氣道受到的軸向力的合力方向應(yīng)該為發(fā)動機順航向；而插板部分與此正好相反，插板部分受到的合力方向應(yīng)該為發(fā)動機逆航向。

圖9 插板位于B 位置壓力分布云圖Fig.9 Pressure cloud of position B

5.2 插板不同半徑位置前后表面壓力分布

圖10 為插板位于進氣道B 位置時的幾何示意圖。此時，R的變化范圍為_13.50～450.00 mm。通過積分的方法求解作用在插板表面的壓力。

圖10 插板位于進氣道B 位置示意圖Fig.10 Schematic of position B

表3 為插板不同半徑位置處的前后表面壓力。其中，p1和p2大小通過前面的數(shù)值計算得到。

表3 插板不同半徑位置前后表面壓力分布(B 位置)Table 3 Pressure distribution at different radius of flashboard (location B)

圖11 為插板前后表面壓差Δp與位置R的多項式關(guān)系曲線。關(guān)系式通過圖中曲線擬合得到。

圖11 插板前后表面壓差與半徑的關(guān)系曲線(B 位置)Fig.11 Pressure difference vs.radius of flashboard (location B)

插板位于進氣道B 位置時，R1=13.50 mm，R2=450.00 mm，p(R)由圖11 中曲線擬合的四次多項式方程給出，由公式(1)計算得到的插板所受的軸向力為15 819 N。

5.3 進氣道唇口前后表面壓力分布

表4 為進氣道唇口不同位置前后表面壓力。其中，和的大小通過前面的數(shù)值計算得到。此時，R′的變化范圍仍為465.00～805.00 mm。

表4 進氣道唇口前后表面壓力分布(B 位置)Table 4 Pressure distribution at different radius of intake (location B)

圖12 為進氣道唇口前后表面壓差Δp′與位置R′的多項式關(guān)系曲線。關(guān)系式通過擬合得到。

圖12 進氣道唇口前后表面壓差與半徑關(guān)系曲線(B 位置)Fig.12 Pressure difference vs.radius of intake (location B)

插板位于進氣道B 位置時，=465.00 mm，=805.00 mm，p(R′)由圖12 中曲線擬合的三次多項式方程給出，由公式(2)計算得到的進氣道唇口所受的軸向力為2 636 N。

根據(jù)前面計算結(jié)果，插板位于進氣道B 位置時，進氣道唇口和插板部分產(chǎn)生的軸向力的合力大小為13 183 N，方向仍為發(fā)動機逆航向。

6 結(jié)論

以帶有插板裝置的喇叭口型進氣道為研究對象，對2 個典型插板位置進行了數(shù)值模擬與數(shù)據(jù)分析，分別求得進氣道唇口和插板所受軸向力?？梢缘玫饺缦陆Y(jié)論：

(1) 帶插板的進氣道受到的軸向力主要來自于進氣道唇口和插板2 個部分，且2 個分力方向相反；

(2) 插板位于進氣道A 位置時，進氣道唇口和插板部分產(chǎn)生的軸向力合力大小為107 N，方向為發(fā)動機逆航向；

(3) 插板位于進氣道B 位置時，進氣道唇口和插板部分產(chǎn)生的軸向力合力大小為13 183 N，方向也為發(fā)動機逆航向；

(4) 插板位于進氣道B 位置時，插板本身和進氣道都將受到比較大的軸向載荷(分別為15 819 N和13 183 N)，在插板、進氣道以及固定進氣道機械結(jié)構(gòu)的強度設(shè)計時應(yīng)加以考慮；

(5) 對于固定在試車臺動架上的進氣道，其軸向受力將通過進氣道支架直接傳遞到試車臺架其他部位，進而會影響到發(fā)動機及吊裝發(fā)動機臺架設(shè)備的安全，本文計算結(jié)果可為航空發(fā)動機試驗臺架設(shè)備的安全設(shè)計提供參考。