基于多场耦合的涡轮叶片静强度计算分析

徐 健,霍施宇,孙 鹏

(中国飞机强度研究所强度与结构完整性全国重点实验室, 西安 710065)

0 引言

航空发动机涡轮冷却叶片在高温高压的燃气环境中以高转速工作,服役过程中受到离心力、气动力等载荷作用,叶片的关键区域(榫齿与伸根部位)在载荷作用下产生了高应力,导致叶片易在关键区域发生破坏[1]。航空发动机涡轮叶片设计首先要满足其静强度条件,保障叶片的正常工作。陈立杰等[2]考虑离心负荷与气流力,对低压涡轮叶片进行了静强度分析,其中气流力沿叶高呈线性分布简化。王婧超[3]等对全三维涡轮叶片进行参数化建模,采用多岛遗传算法及二次序列规划法优化涡轮叶片性能。李立州等[4]基于耦合量曲面拟合的协作优化方法,对巡航状态涡轮转子叶片进行分析,比较了结构和气动子系统中流固耦合面上的温度、压力和变形。Velikanova等[5]采用ANSYS软件中的有限元法和基于杆理论的经典方法对涡轮叶片静强度进行了分析,其中作用在涡轮叶片上的压力场通过 Fluent程序计算获得。郁大照等[6]采用热-应力耦合对两级涡轮叶片开展静强度分析,考虑了叶片榫头与盘外缘榫槽接触的影响。王荣桥等[7]提出了涡轮叶片的多学科优化设计策略,采用流-固-热松散耦合分析方法,从保证结构静强度设计和提高气动等熵效率2个方面,对涡轮带冠叶片进行优化,优化设计目标(叶片总质量和效率)得到明显改善。于复磊等[8]采用单向流-热-固耦合的计算方法获得了在设计工况下离心、气动、热载荷对叶片应力、变形的影响,并按照斯贝发动机EGD-3应力标准对叶片强度进行了校核。彭茂林等[9]通过热-流-固耦合分析对涡轮叶片开展三维气动设计优化和结构可靠度改进。Plesiutschnig等[10]通过金相分析、有限元分析和断裂力学分析,并结合试验数据,计算得出涡轮叶片根部的腐蚀受到离心载荷的影响,导致根部裂纹扩展。沈伟等[11]采用单向流固耦合计算,对涡轮叶盘和叶片的强度分别进行校核。段红燕等[12]采用热力耦合计算,得到了涡轮叶片在离心、温度及气动3种载荷联合作用下的应力、应变分布规律。褚召丰等[13]采用流-热-固单向耦合方法带冠涡轮叶片进行了结构静力分析和应力模态分析。王文瑞等[14]考虑气动力、热应力以及离心力的情况下,对涡轮叶片动态环境下的应力分布进行分析。张筠松等[15]采用流固耦合计算方法研究了涡轮动叶叶顶凹槽喷气冷却的流动与传热特性。梅志恒等[16]通过气-热-固多场耦合仿真计算研究了涡轮叶片典型工况下的应力、应变及形变。杨杰等[17]研究了热载荷和离心载荷耦合工况下涡轮叶片的静强度。

综上所述,现有研究大多在对叶片进行流-热-固耦合分析过程中,流场计算采用理想气体,叶片表面温度分布采用绝热壁面下的流场计算结果,未考虑对流换热和叶片冷却的作用。然而,对于涡轮叶片,不同的气动载荷和热载荷分布规律对叶片变形的影响不尽相同。基于此,本文在考虑工质变比热和叶片在发动机工作过程中表面实际温度场分布规律的情况下,采用单向流-热-固耦合方法建立了单级低压涡轮叶片数值模型,对涡轮叶片进行了强度分析,为开展叶片试验提供强度方法。

1 基于热-流-固耦合分析的涡轮叶片静强度分析方法

基于ANSYS Workbench平台,首先进行单级涡轮的气动定常数值模拟,然后气动载荷和稳态热载荷分别从气动定常计算模块(CFX)和稳态热计算模块(Static-State Thermal)传递到叶片表面,静强度计算模块(Static Structural)用于求解叶片的变形、等效应力等,具体计算流程如图1所示。

通过流-热-固耦合分析能够了解涡轮叶片结构强度的薄弱危险点,从而得到涡轮叶片的关键考核部位,有利于对叶片几何构型进行早期的设计改进。基于热-流-固耦合方法的涡轮叶片静强度分析流程如图2所示。

图2 基于热-流-固耦合方法的涡轮叶片静强度分析流程Fig.2 Static strength analysis process of turbine blades based on thermal-fluid-solid coupling method

涡轮叶片强度计算分析的基本流程如下: 建立涡轮叶片的几何模型,作为强度计算分析的对象;建立叶片温度场分布模型,研究特定工况下叶片温度场的分布;建立叶片强度分析有限元模型,由叶片服役工作状态,设定其气动载荷、边界条件,并将叶片温度场计算结果作预定义场赋给有限元模型;分析叶片整体应力应变状态,确定叶片上的危险点;对几何结构复杂的危险点,采用子模型技术细化网格,求其精确解;基于强度计算准则(安全系数),分析叶片整体的强度。

2 流场仿真分析

2.1 流场计算模型

本研究主要对单级涡轮的通流部分,即一级叶栅流道进行气动分析。单级涡轮动静叶栅模型包括一排静叶栅和一排动叶栅。其中静叶栅包括65个静叶片,动叶栅包括50个动叶片,为了减少计算量,采用叶片约化法,建立了13个静叶和10个静叶组成的1/5个圆周的扇形叶栅模型,边界处设置周期边界条件。这样既不改变动静叶片的相对位置,又能大大减少计算量,计算区域如图3所示。

图3 气动分析计算域Fig.3 Computational domain of aerodynamic analysis

2.2 网格划分

为了减少计算域的边界效应对计算结果的影响,在静叶入口和动叶出口均需要流出一定长度的进出口段。划分网格时,只需要划分一个叶片通道的网格,其他叶片网格在计算中进行圆周阵列。在近壁面处由于流动过程中边界层的存在,垂直壁面方向上流动参数变化剧烈,因此对于壁面处的网格要进行加密,以减少边界层计算的误差。近壁面处采用O型网格,方便在垂直壁面方向网格加密,得到高质量的边界层网格。远壁面处采用H型网格,以提高流道网格的质量。对于叶片的前缘尾缘,由于几何形状尖锐,容易发生流动分离,所以同样进行网格的加密,网格形式如图4所示。

为了保证计算结果的准确性,对网格无关性进行了验证。综合考虑计算资源与结果准确性,采用单流道模型进行分析。网格划分标准与边界条件和上一节所述一致。

2.3 求解设置

采用CFX软件进行数值模拟分析。计算设置需要给定工质参数和各面的边界条件,湍流模型采用SST模型,转速为3 000 r/min,入口和出口的温度压力设置参照表1。

考虑流动工质变比热的影响,定压比热Cp与温度的关系是用以下多项式来确定:

Cp=829.2+0.506 8·T-1.925 4×10-4·T2+

2.736 4×10-8·T3

(1)

气体的粘度由萨瑟兰方程来描述:

(2)

如图5所示,计算区域模型为全周进气,采用了1/5个圆周的通流区域,包含13个静叶和10个动叶,设置周期性边界条件,使前后周期性流动参数对应。在定常计算中,动静交界面采用混合面法。

图5 计算区域的边界条件Fig.5 Boundary conditions of computation region

2.4 网格无关性分析

为了保证计算结果的准确性,本文对网格无关性进行了验证。综合考虑计算资源与结果准确性,采用扇形流道模型进行分析。网格划分标准与边界条件与3.3节所述一致,表2展示了不同网格数的计算结果。表2中可以看出,当网格由389万增加至596万时,流量变化为0.39%,因此网格数389万能够满足网格无关性要求,后续网格均按照此标准进行划分。

表2 网格无关性验证Table 2 Grid independence verification

2.5 流场计算分析

气动载荷是通过流场的气动定常计算获得的,为叶片进行静强度分析中施加气动载荷边界条件提供初始数据。图6为计算得到的叶片表面压力分布云图。

图6 叶片表面压力分布Fig.6 Distribution of static pressure on the blade surface

3 叶片流固耦合分析

3.1 固体计算模型

叶片主要分为以下结构,叶身、缘板及伸根,榫齿。利用NX-UG软件,建立叶片几何模型。如图7所示,三叶片三维模型由叶片和对应的叶轮模型组成。通过对该模型进行三维接触有限元强度分析,获得叶片-叶轮结构的变形及应力分布状况。

图7 三叶片三维模型Fig.7 3D model of three blades

3.2 网格划分

如图8所示,为3个叶片-叶轮模型网格图,除了在叶根平台下方与叶根齿上方之间的过渡区域采用四面体及四棱锥网格划分,模型其他区域均采用8节点六面体单元进行网格划分,模型节点总数为741 011,单元总数为822 049。

图8 固体域网格划分Fig.8 Meshing in the solid domain

3.3 求解设置

完成模型的网格划分以后,下面对模型边界条件的设置进行介绍:

1) 建立接触对。在工作转速下,叶根和轮缘齿面、相邻叶片围带接触面间相互接触,采用面-面接触单元对叶根接触齿面与轮缘接触齿面。

2)设置位移约束。对叶轮进气侧端面上节点施加切向和轴向位移约束;对叶轮出气侧端面上节点施加切向位移约束,并耦合轴向自由度;对叶轮轴线处上节点施加径向位移约束。

3)周期对称约束。对叶轮切向面上节点施加周期对称边界条件。

4)离心力载荷:对模型施加工作转速3 000 r/min下的离心力载荷,即沿轴向施加大小为314.16 rad/s的角速度载荷。

5)气动力载荷:将流场计算得到的压力场加载到叶片表面。

6)温度载荷:热应力对于叶片在高温条件下运行的静强度影响是不可忽视的。

图9为通过实验数据拟合得到的叶片表面温度分布曲线。为了简化计算,在稳态热计算模块(static-state thermal)中采用多项式表达式的形式将温度场数据加载到叶片表面有限元网格的节点上,叶片表面温度分布如图10所示。

图9 叶片表面温度梯度随叶高分布曲线Fig.9 The distribution curve showing the gradient of blade surface temperature with blade height

图10 叶片表面温度分布Fig.10 Temperature distribution on the blade surface

3.4 强度计算分析

采用三维非线性有限元接触方法分析了叶片强度,获得了叶片-轮盘模型的等效应力分布。如图11和图12所示,最大等效应力位于叶片与叶根平台过渡圆角处,其值为401.78 MPa。涡轮叶片的屈服强度为582.5 MPa,计算得到叶身安全系数约为1.45,说明工作转速下该涡轮叶片静应力满足设计要求。

图11 中间叶片Von Mises等效应力分布Fig.11 Von Mises equivalent stress distribution of intermediate blade

图12 叶身圆角处Von Mises等效应力分布Fig.12 Von Mises equivalent stress distribution of blade rounding angle

4 结论

通过对涡轮叶片在温度载荷、离心载荷和气动载荷共同作用下的三维接触有限元应力-应变分析,得到的结论如下:

1) 建立了涡轮叶片流-热-固单向耦合的静强度计算流程,获得了工作转速3 000 r/min时涡轮叶片表面压力场分布规律。

2) 采用三维接触有限元方法对涡轮叶片三叶片模型进行了静强度计算和分析,获得了工作转速3 000 r/min时三叶片-叶轮结构的变形及等效应力分布状况。

3) 在进行叶片强度校核过程中,将气动载荷和热载荷的影响包含在叶身的安全系数计算中,可以减少叶片强度计算带来的误差。