基于流热固耦合的航空发动机涡轮叶片仿真分析

徐建新，许立敬

中国民航大学，天津 300300

涡轮叶片作为航空发动机核心机重要的热端部件之一，工作环境比较恶劣，不仅要承受高温高压的燃气作用，还要承受振动等多种载荷的作用[1]。涡轮叶片的可靠性及寿命问题影响着整个航空发动机甚至整架飞机的安全运营。因此，为后续对航空发动机寿命进行预测，对航空发动机高压涡轮叶片进行仿真分析是至关重要的。但是通过直接利用实物叶片进行试验的经济性不好，周期长，需要花费大量的人力、物力及时间成本，并且试验条件与实际飞行条件有一定的差别，有时甚至差别很大，这直接影响试验的准确性。现在计算机技术的发展尤为迅速，使得有限元方法成为应用比较广泛的一种计算机辅助工程（CAE）方法。本文使用有限元方法利用实际的快速存取记录器（QAR）飞行数据，借助计算机技术和有限元理论对涡轮叶片进行建模、分析研究，周期性较短，精度较高，具有一定的参考价值和使用意义。

目前，大量研究学者对涡轮叶片强度问题进行了仿真分析研究，梅志恒[2]利用飞行循环的载荷谱，结合有限元分析技术，使用有限元分析软件对涡轮叶片进行气、热、固多场耦合仿真计算，获得叶片的应力应变图，分析造成涡轮叶片失效的主要原因，得到涡轮叶片需要考核的部位即应力集中疲劳危险点，求得涡轮叶片的寿命并与试验数据进行对比，结果证明仿真分析得出的结果是正确的。肖力伟[3]针对某燃气高压轮机涡轮动叶进行了流热固耦合数值模拟研究，并在此基础上对其静强度和蠕变寿命进行了计算分析。李广新[4]通过ANSYS Workbench建立流热固耦合计算平台，对某新型燃气轮机跨声速压气机叶轮进行了强度分析，为此类叶轮强度分析和压气机的可靠性设计提供了工程设计计算算法。王小宏[5]为真实反映涡轮叶片的受力情况，利用ANSYS软件进行耦合分析，得到了涡轮叶片的应力应变分布，结果表明叶身根部的吸力面为叶片的疲劳失效危险点。侯甲栋[6]等针对某型航空发动机风扇叶片，考虑离心力和气动力共同作用的影响，进行了有限元建模和强度分析。Liu Donghuan 等[7]通过有限元方法对涡轮叶片进行仿真分析，并以Lemaitre-Chaboche损伤模型为基础，通过修改弹性模量考虑蠕变损伤效应，求得了涡轮叶片的寿命并与传统方法θ进行对比，结果表明了仿真分析结果的正确性和可行性。

本文基于某型飞机实际QAR飞行数据，利用有限元方法并通过有限元软件ANSYS Workbench 对该型航空发动机涡轮叶片进行流热固耦合仿真分析，得到涡轮叶片在工作状况下应力、应变及变形的情况，并对其进行分析，能够为后续航空发动机涡轮叶片寿命预测提供参考。

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.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所示。

涡轮叶片在离心力、气动力以及热应力的共同作用下的变形如图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

本文以某型民用航空发动机涡轮叶片为例，基于SolidWorks 软件建立涡轮叶片三维实体模型，并结合飞机实际QAR飞行数据，对该型航空发动机涡轮叶片进行流热固耦合有限元仿真分析，得出以下结论：

（1）涡轮叶片在最大巡航工作状态下应力集中区域以及应力最大值区域在叶根部位，最大等效应力值为4601.4MPa，且与其他部位应力值相差较大，主要原因是涡轮叶片受到离心力和气动力共同作用。

（2）涡轮叶片叶尖变形量最大，最大等效应变值为0.026，并且由叶根向叶尖变形量依次增大，也是造成涡轮叶片疲劳的主要原因。

（3）通过模拟仿真找出了涡轮叶片危险部位以及危险部位的应力应变值，为以后该型航空发动机涡轮叶片寿命预测提供了参考。