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

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

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

0 引言

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

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

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

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

1.1 葉片模型選取

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

圖1 Rotor37轉子單片葉片幾何模型

1.2 葉片幾何建模流程

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

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

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

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

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

圖3 NASA rotor37轉子的鳥撞有限元模型

表1 葉片材料參數（TC4）

表2 鳥體材料參數

葉片受到鳥撞的數值模擬殘余變形網格模型如圖4所示。

圖4 鳥撞葉片的數值模擬殘余變形網格模型

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

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

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

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

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

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

表3 含變形葉片的Rotor37轉子整級葉片的計算域網格參數

2.2 數值模擬方法

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

圖6 受損Rotor37轉子整級葉片流場計算網格

3 含鳥撞變形葉片的壓氣機氣動性能數值模擬結果分析

3.1 Rotor37轉子鳥撞前后氣動特性對比

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

圖7 Rotor37轉子受到鳥撞受損前后氣動性能對比

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

表4 鳥撞前后轉子氣動特征參數對比

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

3.2 鳥撞后的Rotor37轉子流場分析

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

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

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

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

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

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

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

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

3.2.2 近失穩工況下的流場分析

Rotor37轉子在受到鳥撞后在近失穩工況下不同軸向截面相對馬赫數分布如圖11所示。

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

圖11 Rotor37轉子在受到鳥撞后近失穩工況下不同軸向截面相對馬赫數分布

圖12 Rotor37轉子在受到鳥撞后在近失穩工況下變形葉片附近區域相對馬赫數為0.2的等值面

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

圖13 Rotor37轉子在受到鳥撞后在近失穩工況下不同S1流面相對馬赫數分布

圖14 Rotor37轉子在受到鳥撞后在近失穩工況下遠離變形葉片區域相對馬赫數為0.2的等值面

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

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

4 結論

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

（2）Rotor37轉子氣動試驗數據驗證了本文提出的氣動性能數值模擬方法合理有效；

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

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