基于流熱固耦合的航空發(fā)動(dòng)機(jī)渦輪葉片仿真分析

徐建新，許立敬

中國(guó)民航大學(xué)，天津 300300

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

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

本文基于某型飛機(jī)實(shí)際QAR飛行數(shù)據(jù)，利用有限元方法并通過(guò)有限元軟件ANSYS Workbench 對(duì)該型航空發(fā)動(dòng)機(jī)渦輪葉片進(jìn)行流熱固耦合仿真分析，得到渦輪葉片在工作狀況下應(yīng)力、應(yīng)變及變形的情況，并對(duì)其進(jìn)行分析，能夠?yàn)楹罄m(xù)航空發(fā)動(dòng)機(jī)渦輪葉片壽命預(yù)測(cè)提供參考。

1 流熱固耦合分析理論依據(jù)

1.1 流體控制方程

渦輪葉片處于高溫高壓燃?xì)獾墓ぷ鳝h(huán)境中，與高溫高壓燃?xì)膺M(jìn)行熱量交換，是典型的流熱固耦合問(wèn)題。假設(shè)氣體對(duì)航空發(fā)動(dòng)機(jī)渦輪葉片是不可壓縮的。在實(shí)際運(yùn)行中，葉片的穩(wěn)定流動(dòng)遵循質(zhì)量守恒、動(dòng)量守恒、能量守恒三大物理守恒定律[8]，三種基本守恒定律對(duì)應(yīng)相應(yīng)的控制方程。質(zhì)量守恒定律[9]即由單位時(shí)間內(nèi)流出控制體的流體凈質(zhì)量等于同時(shí)間間隔控制體內(nèi)因密度變化而減少的質(zhì)量，質(zhì)量守恒對(duì)應(yīng)的連續(xù)性方程為

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

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

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

能量守恒定律本質(zhì)是熱力學(xué)第一定律，無(wú)論氣體在葉片流體域的耗散如何，能量守恒定律都能得到滿(mǎn)足。微體中能量的增加率等于進(jìn)入微元體的凈熱流通量加上質(zhì)量力與表面力對(duì)微元體所做的功，其表達(dá)式為

式中，E為流體微團(tuán)的總能；hj為組分j的焓；keff為有效熱傳導(dǎo)系數(shù)；Jj為組分j的擴(kuò)散通量；sh為包括了化學(xué)反應(yīng)熱及其他體積熱源項(xiàng)。

1.2 固體控制方程

由流體誘發(fā)固體振動(dòng)、位移的控制方程[10]為

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

1.3 傳熱控制方程

傳熱基本方程[11]為

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

1.4 耦合控制方程

流熱固交界面處應(yīng)滿(mǎn)足流體與固體的應(yīng)力、位移、熱流量、溫度等相等，其表達(dá)式為

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

2 渦輪葉片有限元分析

2.1 渦輪葉片三維幾何模型

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

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

2.2 渦輪葉片流體分析

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

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

本文以巡航最大功率狀態(tài)為例，邊界條件的獲取主要從QAR飛行數(shù)據(jù)結(jié)合熱力計(jì)算的方式得到，該型號(hào)民用航空發(fā)動(dòng)機(jī)設(shè)計(jì)轉(zhuǎn)速為15183r/min，最大功率為97.7%設(shè)計(jì)轉(zhuǎn)速，選取該型號(hào)發(fā)動(dòng)機(jī)的某次飛行循環(huán)數(shù)據(jù)，在此工況下，高壓壓氣機(jī)出口溫度T3為517℃，高壓壓氣機(jī)出口壓強(qiáng)p3為196PSIA，通過(guò)這些QAR 數(shù)據(jù)并結(jié)合熱力公式計(jì)算渦輪葉片進(jìn)口溫度T4、進(jìn)口壓強(qiáng)p4以及出口壓強(qiáng)p4.5，其中高壓渦輪進(jìn)口溫度T4計(jì)算公式[15-16]為

式中，f為油氣比；ηb為燃燒效率；Hμ為燃油低熱值；cp為空氣比定壓熱容；cpg為燃?xì)獗榷▔簾崛荨Ｆ渲胁殚唴⒖嘉墨I(xiàn)[17]可知，f= 0.03；Hμ=42900kJ/kg；ηb=0.98；cp=1.005kJ/(kg·K)；cpg=1.224kJ/(kg·K)。

渦輪進(jìn)口壓強(qiáng)p4計(jì)算公式為

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

渦輪出口壓強(qiáng)p4.5計(jì)算公式為

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

因此，根據(jù)計(jì)算結(jié)果可得到有限元仿真所需的邊界條件參數(shù)，具體數(shù)值見(jiàn)表1。

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

最后在CFX-post 后處理模塊查看渦輪葉片溫度和壓力云圖，如圖3、圖4所示。由圖3、圖4可知，渦輪葉片前緣的溫度最高，從前緣到后緣溫度逐漸減小，這是由燃?xì)饬鲃?dòng)的方向決定的，葉片前緣附近的燃?xì)鉁囟茸罡撸瑫r(shí)由于流道內(nèi)燃?xì)馑俣炔粩嘣龃螅錅囟戎饾u降低導(dǎo)致葉片表面溫度從前緣到后緣依次降低；壓力分布特點(diǎn)與溫度分布基本一致，渦輪葉片壓力從前緣到后緣逐漸減小，但渦輪葉背出現(xiàn)溫度、壓力驟減的現(xiàn)象，是因?yàn)闅怏w流經(jīng)渦輪葉片葉背處出現(xiàn)了氣流分離的現(xiàn)象，致使葉背部分區(qū)域未與高溫高壓的燃?xì)饨佑|。

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

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

2.3 渦輪葉片熱分析

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

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

表6 DZ125熱導(dǎo)率Table 6 Thermal conductivity of DZ125

然后進(jìn)行網(wǎng)格劃分，由于只針對(duì)渦輪葉片進(jìn)行熱分析，與流場(chǎng)無(wú)關(guān)，故需要將流場(chǎng)進(jìn)行抑制，并對(duì)渦輪葉片前緣、尾緣以及榫頭區(qū)域進(jìn)行網(wǎng)格加密，最后將流體分析中得到的流固耦合交界處溫度載荷分布通過(guò)采用面載荷的方式加載到Steady-State Thermal 模塊進(jìn)行穩(wěn)態(tài)熱分析，得到渦輪葉片的溫度載荷。結(jié)果如圖5所示。

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

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

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

完成熱分析以后，對(duì)渦輪葉片進(jìn)行靜力學(xué)分析完成流熱固耦合，將流體分析中得到的壓力載荷和穩(wěn)態(tài)熱分析中的溫度載荷施加到Static-Structural 中，并施加約束條件和離心力，繼而完成該型號(hào)發(fā)動(dòng)機(jī)渦輪葉片流熱固耦合分析，其流程如圖6所示。

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

其中，離心力是由于葉片繞軸高速轉(zhuǎn)動(dòng)而產(chǎn)生的，所帶來(lái)的損傷也占總損傷的很大部分，通過(guò)施加繞渦輪盤(pán)中心轉(zhuǎn)動(dòng)速度實(shí)現(xiàn)，最后得到渦輪葉片的變形云圖、等效應(yīng)變?cè)茍D以及等效應(yīng)力云圖，如圖7～圖9所示。

3 流熱固耦合結(jié)果分析

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

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

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

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

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

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

表5 DZ125膨脹系數(shù)Table 5 Expansion coefficient of DZ125

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

圖8 渦輪葉片等效應(yīng)變?cè)茍D（單位：mm）Fig.8 Equivalent strain cloud map of turbine blade

4 結(jié)論

本文以某型民用航空發(fā)動(dòng)機(jī)渦輪葉片為例，基于SolidWorks 軟件建立渦輪葉片三維實(shí)體模型，并結(jié)合飛機(jī)實(shí)際QAR飛行數(shù)據(jù)，對(duì)該型航空發(fā)動(dòng)機(jī)渦輪葉片進(jìn)行流熱固耦合有限元仿真分析，得出以下結(jié)論：

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

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

（3）通過(guò)模擬仿真找出了渦輪葉片危險(xiǎn)部位以及危險(xiǎn)部位的應(yīng)力應(yīng)變值，為以后該型航空發(fā)動(dòng)機(jī)渦輪葉片壽命預(yù)測(cè)提供了參考。