火箭发动机涡轮叶片动态响应计算及方法

任 众，许开富，朱东华

(西安航天动力研究所，陕西 西安 710100)

0 引言

涡轮叶片作为火箭发动机的关键部件，要求其在高温、高负荷、高转速的工作环境下安全可靠地工作。局部进气的冲击式涡轮为液体火箭发动机中常见的一种涡轮形式，工作过程中叶片受到周期性高速气流激振力的冲击作用，其振动及疲劳问题十分突出和复杂[1-2]。

现有的液体火箭发动机涡轮由于工作时间短，设计过程中主要考虑气动性能及静强度[3-5]，很少考虑由气体激振引起的叶片振动疲劳问题。但随着可重复使用发动机的研制，要求各零部件的服役时间成倍增加，因此，必须对涡轮叶片的高周疲劳问题开展深入的研究。而该问题的核心就是分析叶片在气流激振力作用下的动态响应，获得其振动应力[6-7]。

目前，叶片的动态响应计算分析主要有两种方法。一种是理论模型的方法，将叶片简化成悬臂梁，根据不同的气流激振力的形式获得激振力的简化模型，再建立梁的弯曲振动微分方程，求解方程[8-11]；另一种是三维仿真的方法，首先通过三维非定常流场仿真获得激振力的时域分布，再进行叶片在激振力作用下的振动响应仿真计算[12-13]。其中，第一种方法已在风机叶片及航空发动机压气机叶片中得到了成熟的应用[14-16]，其激振力形式为连续周期性变化的，但局部进气的冲击式涡轮的气流激振力形式有着明显的不同，其激振力形式为周期性脉冲变化，同时气流的冲击马赫数也要高得多，由该激振力产生的动态响应规律会有所不同。目前，对于该种形式的激振力作用下的涡轮叶片动态响应的研究文献较少；而第二种方法理论上可以得到准确的叶片振动响应情况，但涉及到多学科的仿真分析，计算工作量大，占用资源多，计算时间长，不利于工程设计阶段对疲劳问题的快速评估。

本文提出了一套针对局部进气冲击式涡轮叶片动态响应的理论计算模型，并以某液体火箭发动机涡轮为研究对象，采用该方法进行叶片在流体激振力下的动态响应计算，并与三维仿真结果对比，验证结果的准确性，为发动机涡轮高周疲劳设计提供依据。

1 计算方法

1.1 理论模型

1.1.1 受迫振动微分方程

本文的理论模型以悬臂梁受迫振动理论为基础[17-19]。由于涡轮叶片主要是受气动弯矩的作用产生弯曲变形，同时液体火箭发动机涡轮叶片大多数为较短的直叶片，剪切变形以及截面绕中性轴转动惯量对叶片振动的影响较大。因此，将叶片简化成等截面的铁木辛柯悬臂梁，并考虑离心力及阻尼的影响，推导得到梁的受迫振动微分方程为

(1)

初始边界条件

位移边界条件

y(0,t)=0,y′(0,t)=0,y″(l,t)=0,y‴(l,t)=0

其中

f(x,t)=f1(x)·f2(t)

式中：y为梁的位移；x为梁的轴向坐标；t为时间；E为弹性模量；G为剪切模量；ρ为密度；I为截面主惯性矩；J为单位长度梁对截面中性轴的转动惯量；k为截面形状因子，对于矩形截面，k=5/6；A为截面面积；l为梁的长度；c为黏性外阻尼系数；Fa为叶片离心力大小，由于实际工作过程中叶片的离心力沿叶高是变化的，为了使方程有理论解，假设叶片各截面的离心力均为平均离心力；ω为转速；At为叶片截面面积；Rt为叶顶半径；Rh为叶根半径；f(x,t)为叶片单位高度所受的激励力大小；f1(x)为气流激振力沿叶高的分布函数;f2(t)为单位气流激振力随时间的分布函数。

函数f1(x)通过图1中的叶片通道的控制体动量守恒可得[20]

ρg2(x)Ca2(x)Cu2(x)]

(2)

(3)

p=ρgRTg

(4)

式中：下标1为叶片入口截面参数；下标2为叶片出口截面参数；Zb为叶片数；ρg(x)为气流密度的径向分布；Ca为气流轴向速度的径向分布；Cu为气流轴向速度的径向分布；p为气体压力；Tg为气体温度；R为气体常数。

图1 叶片通道控制体示意图Fig.1 Schematic diagram of blade channel control body

函数f2(t)的分布如图2所示，其中局部进气1个周期内流体激振力的作用时间

Δt=εT

(5)

式中：ε为涡轮局部进气度；T为1个旋转周期的时间。

图2 激振力的时间分布函数Fig.2 Time distribution function of excitation force

1.1.2 方程求解方法

采用模态叠加法求解式(1)，求解的形式可表示为[17]

(6)

式中：Yi(x)为梁的第i个主振型函数；qi(t)为第i个广义坐标权函数。在实际求解时，取前3到5阶振型叠加就可得到足够精确的结果。其中，各阶主振型函数即为梁的自由振动的解，令式(1)右侧的激振力为零，再结合悬臂梁的位移边界条件(3)，推导得到梁弯曲振动的特征

(7)

式中λ1、λ2分别为特征方程A1X4+B1X2+C1=0的两个非负解。其中

C1=-ω_i2ρA

求解上述方程的第i个解ω_i，即为梁的第i阶固有频率，对应的振型函数

Yi(x)=[cosh(λ1x)-cos(λ2x)]+

(8)

为了求解各广义坐标下的权函数，将式(6)代入式(1)，并利用振型函数的正交性，推导得

(9)

其中

式(9)左侧为标准的阻尼系统的受迫振动方程，而右侧可看做振幅为D2的周期激振力。将该激励力看做一系列的微脉冲组合，利用杜哈梅积分，得到梁对该激振力的响应

qi(t)=eat[X1cos(bt)+X2sin(bt)]+

(10)

其中

系数X1、X2根据初始条件(2)确定。获得了梁的位移响应分布y(x,t)，可利用材料力学中的关系获得某点A的应力

(11)

式中yA为A点距截面弯曲主轴的距离。

1.2 三维仿真的方法

三维仿真的计算流程见图3，由图可知，整个计算流程包括CFD非定常计算、载荷传递、动态响应计算分析3部分。首先对涡轮整周模型进行三维瞬态流场仿真，获得某个叶片表面的压力分布随时间的变化，即叶片的气体激振力；然后通过插值运算，将流固交界面上的流场压力映射到叶片表面的有限元网格节点上；再进行叶片的动态响应计算计算，获得叶片在激振力作用下的响应情况；在此基础上，结合材料的疲劳性能数据，可以做叶片的高周疲劳分析。整个计算过程在ANSYS Workbench平台下完成。

图3 三维仿真计算方法Fig.3 Calculation method of 3D simulation

2 算例分析

2.1 研究对象

本文以某发动机涡轮为研究对象，计算叶片在气体激振力下的动态响应。该涡轮为单级局部进气冲击式涡轮，工作时高温高压燃气通过单个喷嘴加速，喷出的高速气流冲击涡轮叶片做功。其中，叶片的三维模型见图4(a)，由图可知，叶片为直叶片，主要几何参数如表1所示。

表1 涡轮叶片主要几何参数Tab.1 Main geometric parameters of turbine blades

涡轮叶片材料为高温合金，密度为8 240 kg/m3，工作温度下的弹性模量为146 GPa，泊松比为0.33，屈服极限为860 MPa。

2.2 理论模型求解

根据截面惯性矩相等，将涡轮叶片简化为矩形截面的悬臂梁，如图4所示。简化梁的长L、宽W、高H分别为29.5、9.45、5.44 mm。此外，由于液体火箭发动机涡轮叶片短，刚性大，一般将叶片和盘考虑为一个叶盘系统，因此本文单独将叶片简化为梁结构，未考虑轮盘的阻抗作用，计算结果与实际存在一定的偏差。

图4 叶片模型简化Fig.4 Blade model simplification

首先对梁进行模态计算，计算得到的前三阶固有频率如表2所示，相应的振型如图5所示。由表2可知，理论模型计算的固有频率比三维仿真值偏小，一方面是由于理论模型中假设叶片各截面的离心力均为平均离心力，从而导致旋转刚化作用对固有频率影响的差异；另一方面是由于理论模型将叶片截面简化为矩形截面，截面形状的不同也会导致固有频率的差异。

表2 叶片前三阶频率对比Tab.2 Comparison of the first three frequencies of blades

图5 梁的前三阶振型Fig.5 The first three modes of beam

然后采用模态叠加法求解梁在气流激振力作用下的位移响应及应力响应，其中，由理论方法得到的气体激振力载荷沿叶高的分布规律如图6所示，该激振力在1个旋转周期内的作用时间为0.3 ms；为了对两种方法进行精确比较，理论模型与三维仿真的阻尼系数均取0.2。对前三阶的模态振型叠加，提取了梁的根部A点的应力及顶部B点的位移结果，A点和B点的位置示如图7所示。梁在3个周期内的顶部位移响应如图8中的虚线所示，根部的应力响应如图9中的虚线所示。

图6 梁的激振力分布图Fig.6 Distribution of excitation force on beam

由图8可知，在一次气流激振力的冲击作用下，叶片的响应频率约为4 066 Hz，同时由于阻尼的作用，响应幅值逐渐衰减，经历7个振荡后振幅基本衰减到零，每个旋转周期内该响应重复1次。以该发动机为例，在单次工作时间范围内大约经历40万次上述循环过程；由图9可知，叶片根部的平均应力约为295 MPa，在气流冲击后的第一个响应动应力幅值最大，约为155 MPa，然后由于阻尼的作用动应力幅值逐渐衰减到零。

图7 位置示意图Fig.7 Location diagram

图8 顶部位移响应结果Fig.8 Response results of top displacement

图9 根部应力响应结果Fig.9 Response results of root stress

2.3 三维仿真求解

涡轮的流场三维仿真模型及边界条件见图10。由图可知，计算域由喷嘴通道与动叶通道两部分组成。整个计算域均采用结构网格，单元总数475万。喷嘴入口为总压边界，动叶出口为静压边界，级间采用瞬态转静交接面，其余均为无滑移绝热壁面边界。利用商用软件CFX开展瞬态仿真计算，每个叶片通道设置5个时间步，一共计算10个旋转周期，采用SST湍流模型及高阶差分格式对方程进行离散求解。

图10 流场计算模型及边界Fig.10 CFD model and boundary

计算得到的叶片通道内的马赫数分布如图11所示，由图11可知，一周大概8个叶片受到高速气流的冲击作用。单个涡轮叶片受到的气动弯矩分布如图12所示，由图12可知，单个叶片在一个旋转周期内受到气流的作用时间比例约为0.155，与式(3)计算的作用时间基本吻合，载荷的形状接近图2的方波。

图11 叶片通道马赫数分布Fig.11 Mach number distribution of blade passage

通过上文的非定常流场仿真，获得了叶片的激振力，在此基础上开展动态响应仿真。有限元计算模型及边界如图13所示。由图13可知，模型由叶片及底座两部分组成，均采用六面体结构网格进行划分，网格节点1.25万。首先对叶片进行模态分析，通过对比各阶频率的大小来验证理论梁模型的准确性。然后利用商用软件ANSYS DYNA进行叶片的动态响应仿真，获得叶片的位移及动应力响应分布。

图12 叶片气动弯矩分布Fig.12 Distribution of aerodynamic bending moment of blade

图13 有限元模型及边界Fig.13 FEM model and boundary

模态计算得到的叶片前三阶固有频率如表2所示，由表2可知，三维仿真计算得到的固有频率比理论模型偏大，且随着阶次的增大，偏差逐渐增大，但最大误差不超过10%，说明采用理论梁模型对叶片进行简化的准确性。

计算得到的叶片根部应力响应如图9中的实线所示，由图9可知，三维仿真得到的平均应力约为305 MPa，响应频率为4 000 Hz，响应动应力的最大幅值为175 MPa。与理论模型方法相比，动应力平均值、频率及最大幅值误差分别为3.39%、1.62%、11.4%，从而验证了理论模型进行叶片动态响应计算的准确性。

3 结论

本文建立了一套针对局部进气冲击式涡轮叶片的动态响应的理论计算模型，并以某液体火箭发动机涡轮为研究对象，分别采用理论模型的方法与三维仿真的方法计算了叶片在气流及激振力作用下的动态响应。主要结论如下：

1)根据液体火箭发动机涡轮叶片的特点，将叶片简化成等截面的铁木辛柯悬臂梁，同时考虑离心力及阻尼的影响，推导得到梁的受迫振动微分方程，并采用模态叠加法对方程进行求解。

2)采用流固耦合的方法对叶片动态响应进行三维仿真，首先通过非定常全流场仿真获得气流激振力，在此基础上开展叶片的三维有限元模态分析及动态响应分析。

3)与三维仿真计算方法对比，采用理论模型方法计算得到的叶片固有频率误差不超过10%，同时根部平均动应力低3.39%，响应频率高1.62%，响应最大幅值低11.4%。验证了该理论模型的准确性，可以采用该方法开展发动机涡轮叶片高周疲劳设计。