含鳥撞變形葉片的壓氣機氣動性能數(shù)值模擬

劉雙麗，陳 偉，3，王志強，3，羅 剛

（1.南京航空航天大學能源與動力學院，2.機械結(jié)構(gòu)力學及控制國家重點實驗室：南京 210016；3.遼寧省航空發(fā)動機沖擊動力學重點實驗室，沈陽 110015）

0 引言

航空發(fā)動機吸入飛鳥后，不但可能發(fā)生嚴重的機械損傷，同時極易導致整機性能惡化，推力降低甚至停車，危及飛行安全。為避免上述危險，各國適航管理機構(gòu)頒布了適航規(guī)章，如美國FAR33部，歐洲的CSE800，中國的CCAR33部等，以法規(guī)形式規(guī)定了航空發(fā)動機吸入飛鳥后應具備的最低安全工作能力，如吸入中鳥后需要維持75%推力，吸入大型群鳥后需要維持50%推力等標志性要求。長期以來，圍繞航空發(fā)動機吸鳥后風扇/壓氣機葉片的損傷、部件或整機的動態(tài)響應、適航符合性設計與驗證等開展了大量研究，逐漸意識到發(fā)動機吸鳥后的結(jié)構(gòu)損傷引發(fā)的氣動性能變化，可能是導致發(fā)動機整機性能惡化和失去推力或工作能力的重要原因，因此針對變形受損葉片或受損部件的氣動性能持續(xù)進行了相關研究。Imregun等采用穩(wěn)態(tài)CFD分析方法進行了計算流體力學/有限元法（CFD/FEM）耦合計算，對人為假設的鳥撞擊損傷葉片的氣彈穩(wěn)定性問題進行了研究；Kim等研究了70%轉(zhuǎn)速下受到不同形式損傷的風扇轉(zhuǎn)子葉片的氣動穩(wěn)定性；Bohari等通過在Rotor67轉(zhuǎn)子葉片的前緣部分假設模擬鳥撞損傷變形，采用計算流體動力學方法獲取了風扇在不同轉(zhuǎn)速下的特性變化；Li等研究了Rotor67轉(zhuǎn)子葉尖鳥撞損傷對旋轉(zhuǎn)失速與喘振裕度的影響；Muir等對某渦扇發(fā)動機風扇在起飛階段鳥撞的氣動問題進行了研究；楊杰等基于SPH法鳥撞擊發(fā)動機風扇數(shù)值模擬技術，應用NUMACA流體計算軟件開展了受損葉片的流場數(shù)值計算；羅剛等建立了含有3聯(lián)受損葉片的整級風扇CFD模型，對整級風扇氣動性能變化進行了預估；陸嘉華等在平面葉柵試驗驗證的基礎上，采用數(shù)值模擬方法研究了中鳥對典型風扇轉(zhuǎn)子葉片造成損傷后導致的風扇氣動性能變化。

綜上所述，國外以高精度數(shù)值模擬方法為主，開展了較為深入的研究；國內(nèi)開展了部分建模、數(shù)值模擬與平面葉柵試驗驗證工作，其他研究絕大部分仍停留在抗鳥撞結(jié)構(gòu)分析與設計領域，對于發(fā)動機吸鳥后導致壓氣機轉(zhuǎn)子部件的氣動性能變化和整機性能降低過程缺乏了解，對發(fā)動機性能惡化的作用機理不夠明確，在發(fā)動機吸鳥后氣動性能變化的高精度虛擬分析手段亟需掌握。

本文針對含有鳥撞變形葉片風扇/壓氣機氣動性能分析的需求，建立了NASA Rotor37轉(zhuǎn)子在鳥撞前后的氣動性能CFD分析模型，開展了數(shù)值模擬，采用試驗數(shù)據(jù)進行了方法驗證，并開展了對比分析與方法流程研究工作。

1 鳥撞變形葉片的幾何建模

1.1 葉片模型選取

為研究壓氣機轉(zhuǎn)子葉片受到鳥撞發(fā)生變形后其氣動性能的變化，本文選擇了NASA公開的Rotor37轉(zhuǎn)子葉片作為研究對象，Rotor37轉(zhuǎn)子是NASA Glenn Research Center于20世紀70年代設計的低展弦比的核心壓氣機跨聲速第1級轉(zhuǎn)子，其整級具有36片葉片，葉尖相對馬赫數(shù)為1.48，轉(zhuǎn)速為17188 r/min，該轉(zhuǎn)子的性能參數(shù)在跨聲速壓氣機中具有典型性，并有詳細的試驗數(shù)據(jù)，可作為數(shù)值模擬方法的有效驗證依據(jù)。Rotor37轉(zhuǎn)子單片葉片幾何模型如圖1所示。

圖1 Rotor37轉(zhuǎn)子單片葉片幾何模型

1.2 葉片幾何建模流程

采用基于跡線提取的半自動建模方法建立適用于氣動性能分析的變形葉片幾何輸入模型。首先建立鳥撞葉片的有限元模型，在建模分網(wǎng)時排布葉片跡線的節(jié)點編號；隨后在沖擊動力學分析軟件LS-DYNA中開展鳥撞葉片瞬態(tài)動力學分析，輸出含有上述跡線節(jié)點坐標的文件，通過對節(jié)點編號進行編程排布，生成符合氣動分析軟件NUMECA適用的GEOMTURBO文件，導入NUMECA中生成幾何模型。鳥撞變形葉片的幾何建模流程如圖2所示。

圖2 鳥撞變形葉片的幾何建模流程

1.3 鳥撞變形葉片的幾何模型

建立了小鳥撞擊含2聯(lián)葉片的Rotor37轉(zhuǎn)子葉片有限元模型，將Rotor37轉(zhuǎn)子撞擊前的2聯(lián)葉片幾何模型進行網(wǎng)格劃分，選用Solid164實體單元，采用SPH方法進行鳥體建模，將鳥體設置為長徑比為2的圓柱體模擬鳥，由于葉片尺寸較小，因此取鳥體質(zhì)量為9 g。鳥體撞擊位置設置在50%葉高，鳥體與葉片間采用點面侵徹接觸方式，對葉片進行根部約束及旋轉(zhuǎn)狀態(tài)設置，鳥體以100 m/s的入射速度被100%轉(zhuǎn)速旋轉(zhuǎn)的葉片切割。在建模時，對葉片網(wǎng)格加密區(qū)域的跡線進行了節(jié)點編號，所建立的有限元模型如圖3所示。

葉片材料采用PLASTIC_KINEMATIC模型，具體參數(shù)見表1。鳥體采用彈塑性水動力本構(gòu)模型與多項式方程描述其大變形特征。鳥體材料參數(shù)見表2。

圖3 NASA rotor37轉(zhuǎn)子的鳥撞有限元模型

表1 葉片材料參數(shù)（TC4）

表2 鳥體材料參數(shù)

葉片受到鳥撞的數(shù)值模擬殘余變形網(wǎng)格模型如圖4所示。

圖4 鳥撞葉片的數(shù)值模擬殘余變形網(wǎng)格模型

從圖中可見，模擬圓柱鳥體沿發(fā)動機進氣方向撞擊葉片，在撞擊位置使葉片前緣發(fā)生了波浪形鼓包變形，其中受到鳥體撞擊的葉片中間部位的殘余變形幅值最大，對該模型的變形跡線節(jié)點進行編程排布，生成GEOMTURBO文件并導入NUMECA中建立含變形葉片的整級葉片幾何模型，如圖5所示。

圖5 含變形葉片的整級葉片幾何模型

從圖中可見，模型包含了1片鳥撞后變形的葉片和35片完好葉片，構(gòu)成了整級葉片氣動分析用幾何模型。

2 鳥撞變形葉片的氣動性能數(shù)值模擬方法

2.1 含變形葉片的整級葉片流場計算模型的CFD網(wǎng)格劃分

對整級葉片幾何模型（圖5）進行全3維的黏性流場計算用CFD網(wǎng)格劃分，針對受損轉(zhuǎn)子的網(wǎng)格劃分采用NUMECA軟件中的IGG-Autogrid5模塊。葉片區(qū)域網(wǎng)格的拓撲結(jié)構(gòu)采用O4H型，葉片上、下游網(wǎng)格為H型，葉尖徑向間隙內(nèi)網(wǎng)格為蝶形。第1層網(wǎng)格厚度為0.003 mm，徑向設置了73個網(wǎng)格節(jié)點，葉尖徑向間隙內(nèi)設置了13個網(wǎng)格節(jié)點。該受損轉(zhuǎn)子整級36片葉片通道的總網(wǎng)格量約為4166萬。網(wǎng)格劃分結(jié)果見表3并如圖6所示。

表3 含變形葉片的Rotor37轉(zhuǎn)子整級葉片的計算域網(wǎng)格參數(shù)

2.2 數(shù)值模擬方法

數(shù)值模擬和分析目的在于計算含變形葉片的壓氣機穩(wěn)態(tài)氣動性能，為保證計算的效率和穩(wěn)定性，采用定常計算的方法，開展含變形葉片的全通道全3維流場數(shù)值模擬。受到鳥撞受損后的Rotor37轉(zhuǎn)子整級全3維流場數(shù)值模擬采用NUMECA軟件中的Fine/Turbo模塊完成。采用時間推進法求解3維雷諾平均的NS方程，時間項采用4階Runge-Kutta法迭代求解，空間離散采用2階精度的中心差分格式進行，湍流模型選用S-A模型。Rotor37整級全3維流場定常數(shù)值模擬結(jié)果的可靠性采用文獻[22]中的試驗數(shù)據(jù)進行驗證。

圖6 受損Rotor37轉(zhuǎn)子整級葉片流場計算網(wǎng)格

3 含鳥撞變形葉片的壓氣機氣動性能數(shù)值模擬結(jié)果分析

3.1 Rotor37轉(zhuǎn)子鳥撞前后氣動特性對比

Rotor37轉(zhuǎn)子受到鳥撞受損前后的氣動性能對比如圖7所示。圖中黑色方點為文獻[22]給出的試驗測量得到的完好Rotor37轉(zhuǎn)子特性，藍線為計算得到的完好轉(zhuǎn)子的氣動特性，橙線為計算得到的受到鳥撞后轉(zhuǎn)子的氣動特性。

圖7 Rotor37轉(zhuǎn)子受到鳥撞受損前后氣動性能對比

從圖中可見，鳥撞導致葉片變形后，該轉(zhuǎn)子的氣動性能明顯降低，在所有能穩(wěn)定工作的流量狀態(tài)下，受到鳥撞后轉(zhuǎn)子的壓比和效率都要明顯小于鳥撞前的，相關參數(shù)變化結(jié)果見表4。

表4 鳥撞前后轉(zhuǎn)子氣動特征參數(shù)對比

從表中可見，數(shù)值模擬獲得的轉(zhuǎn)子受到鳥撞前的特性與試驗結(jié)果比較接近，說明本文所采用的數(shù)值模擬方法是可靠的。而從鳥撞前后數(shù)值模擬獲得的轉(zhuǎn)子氣動特性的對比可見，鳥撞后的轉(zhuǎn)子氣動性能較撞擊前的雖然明顯降低，但是鳥撞前后的特性變化規(guī)律一致，說明該數(shù)值模擬結(jié)果能夠較合理地反映出含有變形葉片的轉(zhuǎn)子氣動性能影響的一般規(guī)律，說明本文針對受損葉片所采用的氣動性能計算方法較為有效。

3.2 鳥撞后的Rotor37轉(zhuǎn)子流場分析

3.2.1 堵塞和最大效率工況下的流場分析

在堵塞工況和最大效率工況下，Rotor37轉(zhuǎn)子受到鳥撞后不同軸向截面相對馬赫數(shù)分布分別如圖8、9所示。圖中的4個軸向截面分別位于轉(zhuǎn)子上游、葉排中間、轉(zhuǎn)子下游以及計算域出口。

圖8 Rotor37轉(zhuǎn)子受到鳥撞后不同軸向截面相對馬赫數(shù)分布（堵塞工況）

圖9 Rotor37轉(zhuǎn)子受到鳥撞后不同軸向截面相對馬赫數(shù)分布（最大效率工況）

從圖8中可見，整體而言，在該狀態(tài)下轉(zhuǎn)子通道內(nèi)的流動狀態(tài)較好。在遠離變形葉片的區(qū)域，流場結(jié)構(gòu)較為一致，流動的周向不均勻性并不顯著。但是在變形葉片區(qū)域，流場結(jié)構(gòu)還是呈現(xiàn)出了明顯不同。主要表現(xiàn)為在變形葉片的吸力面與相鄰葉片壓力面所形成的通道內(nèi)存在大范圍的低速區(qū)，其產(chǎn)生原因是葉片變形后，相應徑向位置處的葉型幾何進口角大幅增大，有的甚至超過90°（圖6（d）的54%葉高位置），導致氣流攻角也大幅增大，從而引發(fā)了受損葉片吸力面附面層的流動分離，堵塞了葉片通道，降低了整個轉(zhuǎn)子的流動能力，進而使轉(zhuǎn)子在堵塞工況下的流量減小，這與圖7中的結(jié)果一致。進一步觀察受損葉片通道中低速區(qū)的發(fā)展可見，在該狀態(tài)下隨著氣流向下游的發(fā)展，在主流的摻混作用下，該低速區(qū)的影響逐漸減弱，在轉(zhuǎn)子下游截面，其速度已得到明顯提升。再觀察計算域出口截面發(fā)現(xiàn)，出口截面的相對馬赫數(shù)分布已比較均勻，上游低速區(qū)的影響已基本消失。

對比圖8、9可見，在最大效率工況以及堵塞工況下，流場結(jié)構(gòu)并沒有明顯變化，主要特點還是表現(xiàn)為在變形葉片附近的流動受到了葉片變形帶來的影響；而在遠離變形葉片的區(qū)域，各通道內(nèi)流動的一致性較好。由于受損葉片的變形使來流攻角增大，導致氣流在受損葉片吸力面附近形成了大范圍的分離區(qū)。

在堵塞工況下不同S流面的相對馬赫數(shù)分布如圖10所示。結(jié)合圖8、10可見，由于該受損葉片的變形位置主要在40%葉高以上，在葉根附近沒有明顯變形。從圖10（a）中可見，在10%葉高位置，受損葉片通道的流動狀態(tài)相對較好，未出現(xiàn)明顯的附面層分離現(xiàn)象；而從圖10（b）、（c）中可見，在54%葉高以及95%葉高處由于葉片變形嚴重，都出現(xiàn)了由附面層分離導致的低速流動區(qū)域。從圖10中還可見，3個不同葉高截面上等熵馬赫數(shù)分布由于受到受損葉片通道堵塞的影響，受損葉片前緣上游弓形激波的形態(tài)和位置與其它葉片也有明顯區(qū)別，通道的堵塞作用類似于提高了激波后的反壓，使得激波強度增大，激波后的馬赫數(shù)出現(xiàn)了明顯降低；該葉片弓形激波的位置相比于其它葉片前緣的激波更靠近上游；與受損葉片左側(cè)相鄰的約5個葉片通道（與轉(zhuǎn)子旋轉(zhuǎn)方向相反）的流動都受到了該葉片上游弓形激波的左支——外伸激波的影響，其上游的來流馬赫數(shù)也相應偏低；由于受損葉片外伸激波減弱，在受損葉片左側(cè)第6個葉片通道上游，激波后的馬赫數(shù)迅速提高，形成了局部高馬赫數(shù)區(qū)域；雖然受損葉片吸力面附面層的分離發(fā)生在通道上部，但是其形成的堵塞作用卻是3維的，因為在10%葉高截面上也可以看到，其前緣的外伸激波強度較大，激波后的馬赫數(shù)出現(xiàn)了明顯降低，說明通道內(nèi)附面層分離引起的堵塞作用是空間3維的。

圖10 在堵塞工況下不同S1流面的相對馬赫數(shù)分布

3.2.2 近失穩(wěn)工況下的流場分析

Rotor37轉(zhuǎn)子在受到鳥撞后在近失穩(wěn)工況下不同軸向截面相對馬赫數(shù)分布如圖11所示。

Rotor37轉(zhuǎn)子在受到鳥撞后，在近失穩(wěn)工況下不同S流面相對馬赫數(shù)分布、變形葉片區(qū)域和遠離變形葉片區(qū)域相對馬赫數(shù)為0.2的等值面分別如圖12~14所示。

圖11 Rotor37轉(zhuǎn)子在受到鳥撞后近失穩(wěn)工況下不同軸向截面相對馬赫數(shù)分布

圖12 Rotor37轉(zhuǎn)子在受到鳥撞后在近失穩(wěn)工況下變形葉片附近區(qū)域相對馬赫數(shù)為0.2的等值面

與前2種工況相比，在近失穩(wěn)工況下，Rotor37轉(zhuǎn)子內(nèi)的流動發(fā)生了較明顯變化。主要表現(xiàn)為：在該工況下，變形葉片通道內(nèi)的低速區(qū)范圍明顯增大，并且在其相鄰通道內(nèi)的近葉尖區(qū)域也出現(xiàn)了貫穿整個通道的低速區(qū)。從圖12中的變形葉片區(qū)域的低馬赫數(shù)等值面的位置和形狀可見，2個通道的低速區(qū)在葉片下游連成了較大區(qū)域的低速團。從圖13（c）中可見，出現(xiàn)這一現(xiàn)象的原因應該是，隨著該轉(zhuǎn)子的節(jié)流，轉(zhuǎn)子下游的背壓逐漸增大，轉(zhuǎn)子上游的激波強度也逐漸增大，波后馬赫數(shù)逐漸降低，攻角進一步增大，導致分離區(qū)也進一步增大。

圖13 Rotor37轉(zhuǎn)子在受到鳥撞后在近失穩(wěn)工況下不同S1流面相對馬赫數(shù)分布

圖14 Rotor37轉(zhuǎn)子在受到鳥撞后在近失穩(wěn)工況下遠離變形葉片區(qū)域相對馬赫數(shù)為0.2的等值面

另外，從圖11、13（c）以及圖14中可見，在近失穩(wěn)工況下，與變形葉片間隔4個葉片通道（與轉(zhuǎn)子旋轉(zhuǎn)方向相反）的葉尖區(qū)域也出現(xiàn)了大范圍的低速區(qū)。該低速區(qū)占據(jù)的周向范圍較大，約6個葉片通道。這可能是因為受到變形葉片上游外伸激波影響的區(qū)域來流馬赫數(shù)較小，而該外伸激波影響開始消失的區(qū)域的流速又突然增大，導致這些位置處的葉片前緣激波強度增大，在葉尖泄漏流的共同影響下，這些葉片通道內(nèi)的流動狀態(tài)惡化，引起了流動的分離，形成了大范圍的低速區(qū)。同樣，這些區(qū)域的低速區(qū)也會引起流道的堵塞，其作用類似于受損葉片通道內(nèi)的低速區(qū)，從而進一步導致在間隔了若干個葉片通道后，又形成了1個葉尖局部低速區(qū)。這些葉尖局部低速區(qū)與失速團的形式非常接近，其存在導致該轉(zhuǎn)子很快進入失穩(wěn)工況。此外，這些低速區(qū)的存在也導致總壓損失和熵的增大，使得該轉(zhuǎn)子的壓比和效率明顯降低，這與圖7給出的該受損轉(zhuǎn)子的性能變化情況吻合。

以上針對Rotor37轉(zhuǎn)子受到鳥撞前后的氣動性能以及3種典型工況下的流場結(jié)構(gòu)分析表明，本文針對含鳥撞變形葉片的Rotor37轉(zhuǎn)子所采用的全通道全3維黏性流場數(shù)值模擬方法可行，其結(jié)果合理，利用該方法可以較準確地計算出含變形葉片轉(zhuǎn)子內(nèi)部的流場細節(jié)，并獲得其氣動性能的變化特性。

4 結(jié)論

（1）本文發(fā)展的鳥撞變形葉片的幾何建模流程可有效建立含有變形葉片的風扇/壓氣機氣動分析幾何模型；

（2）Rotor37轉(zhuǎn)子氣動試驗數(shù)據(jù)驗證了本文提出的氣動性能數(shù)值模擬方法合理有效；

（3）Rotor37轉(zhuǎn)子鳥撞前后的氣動性能特征參數(shù)變化以及3種典型工況下的流場結(jié)構(gòu)分析表明，在本文構(gòu)建的鳥撞葉片變形條件下，該轉(zhuǎn)子的壓比、效率等氣動性能特征參數(shù)明顯降低，穩(wěn)定工作邊界明顯縮減；

（4）數(shù)值模擬結(jié)果表明，攻角增大導致的局部氣流分離及并發(fā)的低速流動行為的耦合是轉(zhuǎn)子氣動性能惡化與轉(zhuǎn)子進入失穩(wěn)工況的主要原因。