基于流熱固耦合的航空發動機渦輪葉片仿真分析

徐建新，許立敬

中國民航大學，天津 300300

渦輪葉片作為航空發動機核心機重要的熱端部件之一，工作環境比較惡劣，不僅要承受高溫高壓的燃氣作用，還要承受振動等多種載荷的作用[1]。渦輪葉片的可靠性及壽命問題影響著整個航空發動機甚至整架飛機的安全運營。因此，為后續對航空發動機壽命進行預測，對航空發動機高壓渦輪葉片進行仿真分析是至關重要的。但是通過直接利用實物葉片進行試驗的經濟性不好，周期長，需要花費大量的人力、物力及時間成本，并且試驗條件與實際飛行條件有一定的差別，有時甚至差別很大，這直接影響試驗的準確性。現在計算機技術的發展尤為迅速，使得有限元方法成為應用比較廣泛的一種計算機輔助工程（CAE）方法。本文使用有限元方法利用實際的快速存取記錄器（QAR）飛行數據，借助計算機技術和有限元理論對渦輪葉片進行建模、分析研究，周期性較短，精度較高，具有一定的參考價值和使用意義。

目前，大量研究學者對渦輪葉片強度問題進行了仿真分析研究，梅志恒[2]利用飛行循環的載荷譜，結合有限元分析技術，使用有限元分析軟件對渦輪葉片進行氣、熱、固多場耦合仿真計算，獲得葉片的應力應變圖，分析造成渦輪葉片失效的主要原因，得到渦輪葉片需要考核的部位即應力集中疲勞危險點，求得渦輪葉片的壽命并與試驗數據進行對比，結果證明仿真分析得出的結果是正確的。肖力偉[3]針對某燃氣高壓輪機渦輪動葉進行了流熱固耦合數值模擬研究，并在此基礎上對其靜強度和蠕變壽命進行了計算分析。李廣新[4]通過ANSYS Workbench建立流熱固耦合計算平臺，對某新型燃氣輪機跨聲速壓氣機葉輪進行了強度分析，為此類葉輪強度分析和壓氣機的可靠性設計提供了工程設計計算算法。王小宏[5]為真實反映渦輪葉片的受力情況，利用ANSYS軟件進行耦合分析，得到了渦輪葉片的應力應變分布，結果表明葉身根部的吸力面為葉片的疲勞失效危險點。侯甲棟[6]等針對某型航空發動機風扇葉片，考慮離心力和氣動力共同作用的影響，進行了有限元建模和強度分析。Liu Donghuan 等[7]通過有限元方法對渦輪葉片進行仿真分析，并以Lemaitre-Chaboche損傷模型為基礎，通過修改彈性模量考慮蠕變損傷效應，求得了渦輪葉片的壽命并與傳統方法θ進行對比，結果表明了仿真分析結果的正確性和可行性。

本文基于某型飛機實際QAR飛行數據，利用有限元方法并通過有限元軟件ANSYS Workbench 對該型航空發動機渦輪葉片進行流熱固耦合仿真分析，得到渦輪葉片在工作狀況下應力、應變及變形的情況，并對其進行分析，能夠為后續航空發動機渦輪葉片壽命預測提供參考。

1 流熱固耦合分析理論依據

1.1 流體控制方程

渦輪葉片處于高溫高壓燃氣的工作環境中，與高溫高壓燃氣進行熱量交換，是典型的流熱固耦合問題。假設氣體對航空發動機渦輪葉片是不可壓縮的。在實際運行中，葉片的穩定流動遵循質量守恒、動量守恒、能量守恒三大物理守恒定律[8]，三種基本守恒定律對應相應的控制方程。質量守恒定律[9]即由單位時間內流出控制體的流體凈質量等于同時間間隔控制體內因密度變化而減少的質量，質量守恒對應的連續性方程為

式中：ux，uy，uz為x、y、z三個方向的速度分量；t為時間；ρ為密度。

動量守恒方程滿足牛頓第二定律，即對于給定的流體微元，其動量對時間的變化率等于外界作用在該微元體上的各種力之和。葉片流體域中氣體在每個速度方向上的分量都滿足動量守恒方程。其x、y、z三個方向的動量守恒方程為

式中，τxx、τyx、τzx、τxy、τyy、τzy、τxz、τyz、τzz是因分子黏性作用而產生的作用在微元體表面上的黏性應力τ的分量；fx、fy、fz為三個方向上的單位質量力。

能量守恒定律本質是熱力學第一定律，無論氣體在葉片流體域的耗散如何，能量守恒定律都能得到滿足。微體中能量的增加率等于進入微元體的凈熱流通量加上質量力與表面力對微元體所做的功，其表達式為

式中，E為流體微團的總能；hj為組分j的焓；keff為有效熱傳導系數；Jj為組分j的擴散通量；sh為包括了化學反應熱及其他體積熱源項。

1.2 固體控制方程

由流體誘發固體振動、位移的控制方程[10]為

式中，Ms為質量矩陣；Cs為阻尼矩陣；Ks為剛度矩陣；r為固體的位移；τs為固體受到的應力。

1.3 傳熱控制方程

傳熱基本方程[11]為

式中，k為傳熱系數；A為傳熱面積；Δtm為傳熱的平均溫差。

1.4 耦合控制方程

流熱固交界面處應滿足流體與固體的應力、位移、熱流量、溫度等相等，其表達式為

式中，q為熱流量；T為溫度。

2 渦輪葉片有限元分析

2.1 渦輪葉片三維幾何模型

某型航空發動機渦輪葉片主要由葉身、緣板、榫頭三部分組成，其中葉身為主要受力部分，起能量轉換作用，其表面特殊曲面設計可以改善氣流流向；緣板能夠有效減小葉身根部的應力集中，也可以避免渦輪盤受到高溫燃氣的直接沖刷；榫頭為樅樹形，能將葉片固定到渦輪盤中，同時榫頭側面加大與渦輪盤的接觸面積可以減小接觸應力[12]。由于渦輪葉片幾何形狀比較復雜且曲面較多，故首先通過3D掃描儀得到高壓渦輪葉片.stl格式模型，然后在Magics軟件中重新定位擺正模型，最后在SolidWorks 軟件中擬合描線做出曲面得到有限元模型。渦輪葉片模型如圖1所示。

圖1 渦輪葉片模型（單位：m）Fig.1 Turbine blade model

2.2 渦輪葉片流體分析

對渦輪葉片進行流場分析是完成熱流固耦合分析的第一步，采用Fluid Flow（CFX）模塊進行流體分析，得到流固耦合交界處的溫度分布載荷和氣動載荷分布。首先通過ANSYS Workbench中Geometry模塊將渦輪葉片模型導入，并在ANSYS Design Modeler建立渦輪葉片的流場區域，然后采用四面體網格劃分方法對流體區域進行網格劃分，加密葉片臨界處的網格，優化網格質量，四面體網格的數量是3172641個，節點數是581547個，如圖2所示。最后在CFXpre中進行邊界條件的設定，選擇進口總壓、進口總溫、出口壓強（壓力）作為初始邊界條件，湍流模型選擇k-εSST 模型[13]，此模型被公認為標準的工業模型，計算量和精度都符合葉片所處的物理場環境，高速流體考慮流體動能造成的熱量變化，故選取Total Energy全熱模型[14]。

圖2 渦輪葉片流體域網格劃分（單位：mm）Fig.2 Turbine blade fluid domain meshing

本文以巡航最大功率狀態為例，邊界條件的獲取主要從QAR飛行數據結合熱力計算的方式得到，該型號民用航空發動機設計轉速為15183r/min，最大功率為97.7%設計轉速，選取該型號發動機的某次飛行循環數據，在此工況下，高壓壓氣機出口溫度T3為517℃，高壓壓氣機出口壓強p3為196PSIA，通過這些QAR 數據并結合熱力公式計算渦輪葉片進口溫度T4、進口壓強p4以及出口壓強p4.5，其中高壓渦輪進口溫度T4計算公式[15-16]為

式中，f為油氣比；ηb為燃燒效率；Hμ為燃油低熱值；cp為空氣比定壓熱容；cpg為燃氣比定壓熱容。其中查閱參考文獻[17]可知，f= 0.03；Hμ=42900kJ/kg；ηb=0.98；cp=1.005kJ/(kg·K)；cpg=1.224kJ/(kg·K)。

渦輪進口壓強p4計算公式為

式中，σb為燃燒室總壓恢復系數，查閱參考文獻[18]可知，σb=0.97；其中1PSIA=6.89kPa。

渦輪出口壓強p4.5計算公式為

式中，πTH為渦輪落壓比，查閱參考文獻[19]可知，πTH=2.207。

因此，根據計算結果可得到有限元仿真所需的邊界條件參數，具體數值見表1。

表1 有限元仿真邊界條件參數Table 1 Parameters of finite element simulation boundary conditions

最后在CFX-post 后處理模塊查看渦輪葉片溫度和壓力云圖，如圖3、圖4所示。由圖3、圖4可知，渦輪葉片前緣的溫度最高，從前緣到后緣溫度逐漸減小，這是由燃氣流動的方向決定的，葉片前緣附近的燃氣溫度最高，同時由于流道內燃氣速度不斷增大，其溫度逐漸降低導致葉片表面溫度從前緣到后緣依次降低；壓力分布特點與溫度分布基本一致，渦輪葉片壓力從前緣到后緣逐漸減小，但渦輪葉背出現溫度、壓力驟減的現象，是因為氣體流經渦輪葉片葉背處出現了氣流分離的現象，致使葉背部分區域未與高溫高壓的燃氣接觸。

圖3 渦輪葉片溫度云圖（單位：K）Fig.3 Temperature cloud map of turbine blade

圖4 渦輪葉片壓力云圖（單位：Pa）Fig.4 Pressure cloud map of turbine blade

2.3 渦輪葉片熱分析

在對渦輪葉片進行熱分析之前，需要在Engineering Date 中對渦輪葉片材料屬性進行設定，該型號民用航空發動機渦輪葉片使用定向凝固鎳基鑄造高溫合金材料，材料牌號為DZ125，密度為8.48g/cm3，屬于正交各向異性材料，查閱中國航空材料手冊[17]可獲得彈性模量E、泊松比μ、剪切模量G、膨脹系數α及熱導率λ，見表2～表6。

表2 DZ125彈性模量Table 2 Modulus of elasticity of DZ125

表6 DZ125熱導率Table 6 Thermal conductivity of DZ125

然后進行網格劃分，由于只針對渦輪葉片進行熱分析，與流場無關，故需要將流場進行抑制，并對渦輪葉片前緣、尾緣以及榫頭區域進行網格加密，最后將流體分析中得到的流固耦合交界處溫度載荷分布通過采用面載荷的方式加載到Steady-State Thermal 模塊進行穩態熱分析，得到渦輪葉片的溫度載荷。結果如圖5所示。

圖5 渦輪葉片熱分析溫度云圖（單位：mm）Fig.5 Thermal analysis temperature cloud map of turbine blade

觀察渦輪葉片的溫度云圖，可以看出有以下特點：渦輪葉片葉盆溫度明顯大于葉背且其主流流動穩定，這表明熱傳遞的大小決定于邊界層的流動情況和溫度分布；葉身前緣溫度大于后緣溫度。

2.4 渦輪葉片流熱固耦合分析

完成熱分析以后，對渦輪葉片進行靜力學分析完成流熱固耦合，將流體分析中得到的壓力載荷和穩態熱分析中的溫度載荷施加到Static-Structural 中，并施加約束條件和離心力，繼而完成該型號發動機渦輪葉片流熱固耦合分析，其流程如圖6所示。

圖6 渦輪葉片流熱固耦合分析流程Fig.6 Flow-thermal-solid coupling analysis flow chart of turbine blade

其中，離心力是由于葉片繞軸高速轉動而產生的，所帶來的損傷也占總損傷的很大部分，通過施加繞渦輪盤中心轉動速度實現，最后得到渦輪葉片的變形云圖、等效應變云圖以及等效應力云圖，如圖7～圖9所示。

3 流熱固耦合結果分析

渦輪葉片在離心力、氣動力以及熱應力的共同作用下的變形如圖7 所示，最大位移變形位于渦輪葉片的葉尖部位，最大變形是1.034mm，葉片位移變形從葉根往葉身部位變形逐漸增大，主要是離心力作用的結果，離旋轉軸心越大，離心力越大，故變形越大。

圖7 渦輪葉片變形云圖（單位：mm）Fig.7 Deformation cloud map of turbine blade

表3 DZ125泊松比Table 3 Poisson’s ratio of DZ125

渦輪葉片在離心力、氣動力以及熱應力的共同作用下的等效應變和等效應力分布如圖8 和圖9 所示。從圖9 中可以看出，該型航空發動機渦輪葉片的應力集中區域位于葉身下部向緣板過渡區域，且最大等效應力值出現在渦輪葉片的前緣和尾緣部位，最大等效應力值為4601.4MPa，其他區域等效應力值較小，應力集中區域主要集中在這個區域，主要是氣動力和離心力共同作用的結果，此處的等效應力主要是氣動力引起的彎曲應力和離心力引起的拉應力。此外，本文計算出的應力集中區域與參考文獻所確定的應力集中區域位置比較接近[20]。

圖9 渦輪葉片等效應力云圖（單位：mm）Fig.9 Equivalent stress cloud of turbine blade

表4 DZ125剪切模量Table 4 Shear modulus of DZ125

表5 DZ125膨脹系數Table 5 Expansion coefficient of DZ125

由圖8 可知，渦輪葉片的等效應變變化規律與等效應力分布情況大致一樣，最大應變在渦輪葉片根部靠近前后緣的部位，最大應變值為0.026，應變是造成渦輪葉片低周疲勞的主要原因，解釋了渦輪葉片在服役過程中葉片根部容易發生各種低周疲勞失效行為的原因，同時也證明了仿真結果的正確性。

圖8 渦輪葉片等效應變云圖（單位：mm）Fig.8 Equivalent strain cloud map of turbine blade

4 結論

本文以某型民用航空發動機渦輪葉片為例，基于SolidWorks 軟件建立渦輪葉片三維實體模型，并結合飛機實際QAR飛行數據，對該型航空發動機渦輪葉片進行流熱固耦合有限元仿真分析，得出以下結論：

（1）渦輪葉片在最大巡航工作狀態下應力集中區域以及應力最大值區域在葉根部位，最大等效應力值為4601.4MPa，且與其他部位應力值相差較大，主要原因是渦輪葉片受到離心力和氣動力共同作用。

（2）渦輪葉片葉尖變形量最大，最大等效應變值為0.026，并且由葉根向葉尖變形量依次增大，也是造成渦輪葉片疲勞的主要原因。

（3）通過模擬仿真找出了渦輪葉片危險部位以及危險部位的應力應變值，為以后該型航空發動機渦輪葉片壽命預測提供了參考。