CN101599104A - 一种航空涡轮发动机叶片颤振边界的模拟方法 - Google Patents

Abstract

本发明公开了一种航空涡轮发动机叶片颤振边界的模拟方法，该方法基于能量法原理，采用弱耦合的方法，设计了一套颤振预测方法，通过气动阻尼值得到发动机叶片的颤振边界。该方法首先建立叶片的实体模型，然后对叶片进行模态分析，接着建立叶片的流场模型，获得线性插值叶片的振动位移，采用多层动网格化处理生成流场分析需要的网格文件，调用动网格模块获得流场中各参数，之后获得每个工况下的气动阻尼，最后获得颤振边界。本发明方法基于单向弱耦合，大大简化了计算，节约计算成本，可以实现任意频率和不同流场状态下的气动阻尼计算，具有良好的工程应用性。

Description

一种航空涡轮发动机叶片颤振边界的模拟方法

技术领域

本发明属于航空发动机模拟领域，具体涉及一种基于能量法的航空涡轮发动机叶片颤振边界的模拟方法。

背景技术

20世纪90年代以后，航空发动机的高转速、高推重比、高适用性和耐久性、低成本的发展方向已成为主流趋势。这种趋势使得涡扇发动机压气机的级负荷越来越重，风扇和压气机叶片的工作环境十分严峻。航空发动机中，振动故障占发动机中总故障的60％以上，而叶片故障又占振动故障的70％以上。其中由颤振引起的叶片失效占相当大的比重。如果能够有效的对叶片颤振这种流固耦合现象进行模拟，并提出叶片的结构优化方案，将大大减小实验的周期和材料成本。

对风扇/压气机叶片颤振边界进行模拟，采用计算力学方法来进行模拟，包括能量法和特征值法。能量法的优点是：物理意义明确，指出了颤振的本质是系统内的能量交换。消除颤振就是设法使系统获得的总能量为负值，排除叶片从周围气流中汲取能量的可能性。同时，由于结构动力学分析与非定常气动力的分析各自独立，分析工作大为简化。

采用能量法，需要叶对栅通道中的流场进行分析以确定气动力的分布。以往，对于流固耦合的分析主要依靠非定常气动力计算与模态叠加的方法来完成，这意味着结构计算采用了线性方法处理。然而，实际的气动弹性问题是非线性问题，叠加原理失败，必须探讨全场求解耦合方程的途径。

在固体力学中非线性的有限元分析方法已相当成熟，在流体力学中非线性N-S方程求解的计算流体力学(Computational Fluid Dynamics，简称CFD)技术的进展也十分迅速。两者采用不同的描述运动的方法，固体力学中习惯采用Lagrange坐标系，着眼于质点，而流体力学中更多地使用Euler坐标系，着眼于空间点。这种运动描述方法上的差异，对小运动问题，可不加区分。非定常气动力计算和非线性有限元计算两个系统的耦合技术是实现非线性气动弹性计算的发展方向。

对于叶片颤振这样的流固耦合问题，最大的困难在于固体域和流体域两相界面的协调，如何考虑固体叶片的结构变形与绕叶片的流体运动之间的相互作用，是实现耦合技术的关键。目前还没有一项较为成熟的技术来模拟发动机风扇或压气机叶片的颤振。

发明内容

本发明的目的是模拟发动机叶片的颤振边界。

本发明选择经济性能比较好的弱耦合方法作为模型，它不仅符合叶片小变形假设，并且简单易行，能够降低成本。假设流体对叶片质量和刚度的影响不计，叶片的模态运动体现在使流体作受迫运动上，将叶片振动位移单向传递给流场。

实施步骤如下：

步骤一：建立叶片的实体模型。

建立需要分析的发动机风扇叶片的三维实体模型。

步骤二：对叶片进行模态分析。

将已建立的发动机风扇叶片实体模型导入有限元分析工具ANSYS中，动模态分析获得叶片的固有模态，提取流场分析所需要的频率和模态，导出叶片表面的有限元网格节点信息，即计算结构动力学CSD(Computational Structure Dynamics，简称CSD)节点。

步骤三：建立叶片的流场模型。

将已建立的发动机风扇叶片实体模型导入流体动力学仿真工具CFX中，建立发动机叶片单通道流场模型，划分流体网格。将叶片表面的流体网格节点和单元信息导出，即获得计算流体力学CFD节点。

步骤四：获得线性插值叶片的振动位移。

通过三维线性插值方法，将叶片表面CSD节点上的第P阶振动模态插值到叶片表面的CFD节点上，作为叶片表面所有CFD节点在流场中的振动位移。其中P为正整数，且P∈[1，10]。

步骤五：多层动网格化处理生成流场分析需要的网格文件。

采用多层动网格方法，将叶片表面各个CFD节点的振动位移按各层初始的距离比例，分配到叶片周围M层O型网格区域对应的节点上，由叶片表面向外法向递减，第M层位移为零，最后生成各时刻流场所有节点的网格文件。

步骤六：调用动网格模块获得流场中各参数。

流场分析通过流体动力学仿真工具CFX实现。通道入口给定静温和总压，出口给定静压，循环对称面给定交界面边界条件，叶片表面给定运动网格边界条件。以不加运动网格的定常场作为初始条件。叶片按照第P阶固有频率振动，其运动周期为频率的倒数，每个时间步间隔相同的时间，在每一个时间步上求解N-S方程。输出结果需要设置成每一个时间步输出气动力和对应的节点位移等信息。其中时间步是指叶片一个运动周期的N分之一，N是用户设定的叶片一个周期运动的时间步数，理论上N的值取得越大越好，但实际应用中为保证模拟效率，取N为正整数，且N∈[30，80]。

步骤七：获得每个工况下的气动阻尼。

根据每个时间步叶片表面的气动功和节点位移，按照能量法的概念推导出的等效模态阻尼公式来获得叶片一个振动周期内的气动功和等效模态气动阻尼系数。气动功是按照各节点在各时间步上做的功积分获得，积分区间为叶片的一个振动周期，积分域为叶片表面。

步骤八：获得颤振边界。

获得不同工作状况下的等效模态气动阻尼，获得工作域中等效模态气动阻尼的分布，确定阻尼值为零的等值线为颤振边界。

这套模拟方法适用于各种涡轮风扇发动机的风扇/压气机叶片的颤振边界分析。处理时间短，并具有一定的精度，可以在某种程度上减少实验成本和周期。这套方法中所用到的三维线性插值方法简单并且速度快，而其中的控制因子能够保证插值的精度。多层动网格方法突破了原来CFX中只能运动一层流体网格的限制，使得叶片在更大的范围内可以运动，这样就保证了模拟的准确性，否则很有可能在运动的过程中叶片位移超过了叶片周围的流体网格的尺寸而产生负体积，导致运算无法进行。

下面分别介绍该模拟方法中用到的三维线性插值方法和多层动网格方法。

三维线性插值方法的具体实施步骤如下：

步骤一：寻找距离待插值的流场节点F最近的三个固体域的节点A₁、B₁、C₁。

对叶片表面上的CFD网格节点F，在叶片CSD网格节点中找到与其距离最近的三个点A₁、B₁、C₁。需要注意的是点A₁、B₁、C₁不能在一条直线上或接近一条直线。这是因为当A₁、B₁、C₁发生位移以后，三角形面积会发生变化，如果ΔA₁B₁C₁面积几乎为零或者有一个大钝角，会导致ΔA₂B₂C₂面积变化很剧烈，从而导致CFD网格节点发生剧烈的移动，产生较大误差。引入面积控制因子s，限定ΔA₁B₁C₁的面积不小于s，保证A₁、B₁、C₁三点不接近于直线。

寻找节点A₁、B₁、C₁时，遍历所有的节点以后，将距离节点F最近的w个节点存储在一个向量存储器中，该向量存储器里储存的这w个节点是按照与节点F的距离远近排列的，储存的第一个节点离节点F最近，储存的第二个节点离节点F次近，……，储存的第w个节点离节点F最远；将储存的第一个节点作为节点A₁，将储存的第二个节点作为节点B₁，再通过引入面积控制因子s从剩余的数组元素中找到节点C₁，确保节点C₁既不与节点A₁、节点B₁在同一直线上，又确保节点C₁为向量存储器中除节点A₁、节点B₁外距离节点F最近的节点；

步骤二：获取三个固体域节点A₁、B₁、C₁组成的平面及该平面的法向向量

由A₁、B₁、C₁三点坐标(x_A1，y_A1，z_A1)，(x_B1，y_B1，z_B1)，(x_C1，y_C1，z_C1)组成一个平面方程，根据已知点的坐标确定该平面，平面ΔA₁B₁C₁的法向向量为 $\stackrel{\→}{n}=(a_1,b_1,c_1)$

步骤三：根据三个固体域节点发生的位移确定运动后产生的另一平面ΔA₂B₂C₂。

在点A₁、B₁、C₁上分别作ΔA₁B₁C₁的法线，使得A₁A₂、B₁B₂、C₁C₂分别为点A₁、B₁、C₁的位移值，由此得到A₂、B₂、C₂的坐标，然后再根据坐标确定平面ΔA₂B₂C₂；

步骤四：获取插值得到的节点F的位移。

过F点做ΔA₁B₁C₁的法线，在ΔA₁B₁C₁和ΔA₂B₂C₂上的交点分别为F₁、F₂，则F₂F₂为F点的位移值，其正负由

决定，

与同向时，F₁F₂为正；

与

反向时，F₁F₂为负。

为ΔA₁B₁C₁的外法向向量。

多层动网格方法：多层动网格方法将叶片的振动位移按照一定的比例加到叶片周围的M层O型流体网格节点上，一定程度上加大了叶片可移动的范围，使得流体网格在M层的范围内实现运动。

步骤一：找到叶片表面的M层流体域O型网格节点。

对于O型网格，每个流体单元i中的各个节点都按照一定顺序排列，第1，2，6，5节点分别与4，3，7，8节点对应，而每一层面单元的节点排列顺序为顺时针，即4-3-7-8。由叶片表面的第一层网格节点找到第一层流体单元，通过上述的对应关系找到与结构体表面节点对应的第二层节点。依此类推，找出M层O型网格节点。

步骤二：获得各层网格节点与叶片表面对应点的距离及各层对应点位移比例。

叶片表面节点，即第1层网格节点，其最大运动位移即为插值获得的叶片表面各CFD节点振幅，第M层节点最大运动位移为零；各层按照初始的距离比例分配位移。

步骤三：各层网格节点按比例位移移动，获得运动位移的峰值。

按照步骤二中获得的比例系数分配各层的位移将叶片表面各点的最大位移分配到各层节点上。

步骤四：所有网格节点按正弦运动规律输出网格文件。

令M层网格节点按照正弦运动，一个周期离散为N个时间步，生成网格文件。

本发明的优点在于：

(1)基于能量法的弱耦合手段可以使固体域与流体域分别建模，通过线性插值完成位移载荷的传递，从而使处理过程大大简化，节约成本；

(2)线性插值方法简单，可获得任意固有频率及不同流场状态下的气动阻尼，具有很好的工程应用性；

(3)多层动网格方法的使用使叶片最大振幅不受最小网格尺寸的制约，在更大的范围内可以运动，使动网格的分布趋于合理；

附图说明

图1是本发明航空涡轮发动机叶片颤振边界的模拟方法流程图；

图2是现有NASA67转子流场测量位置；

图3是现有NASA67叶片有限元模型；

图4是现有NASA67叶片流场模型；

图5是现有NASA67叶片单通道流场模型；

图6是本发明三维线性插值方法流程图；

图7是现有三维线性插值几何示意图；

图8是本发明多层动网格方法流程图；

图9是本发明流体域O型网格节点排布示意图；

图10是现有叶尖流场网格图；

图11是本发明等效模态气动阻尼在气动特性图上的三维曲面图；

图12是本发明等效模态气动阻尼在气动特性图上的等值线图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。本发明是一种基于能量法的航空涡轮发动机叶片颤振边界的模拟方法，颤振边界的模拟方法流程如图1所示。

步骤一：建立叶片的实体模型。

建立需要分析的发动机风扇叶片的三维实体模型。

步骤二：对叶片进行模态分析。

将已建立的发动机风扇叶片的三维实体模型导入有限元分析工具ANSYS中，并对其进行有限元单元划分，定义材料和单元属性，加上载荷和约束条件。利用ANSYS中的动模态分析获得叶片的固有模态，提取流场分析所需要的频率和模态，导出叶片表面的有限元网格节点信息，即CSD节点。

步骤三：建立叶片的流场模型。

将已建立的发动机风扇叶片的三维实体模型导入流体动力学仿真工具CFX中，建立发动机叶片单通道流场模型，划分流体网格。将叶片表面的流体网格节点和单元信息导出，获得CFD节点。

步骤四：获得线性插值叶片的振动位移。

采用三维线性插值方法，将叶片表面CSD节点上的第P阶振动模态插值到叶片表面的CFD节点上，作为叶片表面所有CFD节点在流场中的振动位移，其中P为正整数，且P∈[1，10]。小范围选点插值保证了处理速度，引入面积控制因子保证插值的精度。

步骤五：多层动网格化处理生成流场分析需要的网格文件。

采用多层动网格方法，将叶片表面各个CFD节点的振动位移按各层初始的距离比例，分配到叶片周围M层O型网格对应的节点上，由叶片表面向外法向递减，第M层位移为零，最后生成各时刻流场所有节点的网格文件。所述多层动网格方法适用于1～M层，其中M为流场模型中的O型网格的层数。多层动网格方法的应用使模拟结果更贴近于实际，解决了固体域模型在流场中运动受到流体网格尺寸制约的问题。

步骤六：调用动网格模块获得流场中各参数。

流场分析通过CFX实现。在步骤三建立的流场模型中加载边界条件和初始条件。通道入口给定静温和总压，出口给定静压，循环对称面给定交界面边界条件，叶片表面给定运动网格边界条件。以不加运动网格的定常场作为初始条件。叶片按照第P阶固有频率振动，其运动周期为频率的倒数，每个时间步间隔相同的时间，在每一个时间步上求解N-S方程。输出结果需要设置成每一个时间步输出气动力和对应的节点位移信息。其中时间步是指叶片一个运动周期的N分之一，N是用户设定的叶片一个周期运动的时间步数，理论上N的值取得越大越好，但实际应用中，为了保证模拟效率，N的取值范围可以从30到80。

步骤七：获得每个工况下的气动阻尼。

根据每个时间步叶片表面的气动功和节点位移，按照能量法的概念推导出的等效模态阻尼公式来获得叶片一个振动周期内的气动功和等效模态气动阻尼。气动功是按照各节点在各时间步上的做功积分获得，积分区间为一个振动周期，积分域为叶片表面。

其中能量法的概念是指根据叶片在一个振动周期内，从外界所得到的能量的正负来判别颤振发作与否，当从外界所得到的能量为正时，则颤振发作；当从外界所得到的能量为负时，则颤振不发作。按照能量法的概念推导出的等效模态阻尼公式是 ${\mathrm{\ξ}}_i=-\frac{W_{\mathrm{cfd}}}{2\mathrm{\π}{q}_{i0}^2{\mathrm{\ω}}_i^2}$ 其中W_cfd为第i阶模态下状态空间坐标系中非定常气动力做的功，q_i0为第i阶模态的振幅，ω_i为叶片振动的第i阶固有频率。

步骤八：获得颤振边界。

获得不同工作状况下的等效模态气动阻尼，获得工作域中等效模态气动阻尼的分布，确定阻尼值为零的等值线为颤振边界。

下面再通过一个实例来说明此套方法。首先选用NASA67叶片，采用UG建立需要实体模型，其中NASA67是美国宇航局公布的一种发动机风扇叶型，UG是一种三维实体仿真工具。基本的几何参数以及气动参数参见表1。图2则给出了NASA67转子流场的激光测速位置以及气动参数测量站的位置。

表1  NASA67转子基本设计参数

表1给出了实验中NASA67叶片的各项参数，包括几何参数和进出口边界条件，以下的模拟都是按照表中给出的参数进行。

将已建立的发动机风扇叶片实体模型导入ANSYS中，建立了NASA67叶片的有限元模型如图3，参数见表2。X方向为转子周向，Y方向为叶片展向，Z方向为转子轴向。单元类型为solid45，X、Y、Z三个方向共2×20×22个单元，1449个节点，其中solid45是指有限元实体单元。叶片工作转速为16043r/min，约束方式为叶根固支全约束，叶尖自由。可以对比一下叶片表面CFX网格节点，如图4，共10206个节点，可以看到CFX节点要比CSD节点密很多，需要插值获得叶片各阶振动模态。利用ANSYS中的动模态分析获取叶片的固有模态，提取流场分析所需要的频率和模态，导出叶片表面的结构动力学有限元网格节点信息(CSD节点)。

表2  NASA67叶片有限元模型参数

表2给出了有限元模型的参数。将已建立的发动机风扇叶片实体模型导入CFX中，叶片流场通道的CFX模型如图5所示，无叶尖间隙。将叶片表面的流体网格节点和单元信息导出。

通过三维线性插值方法，将叶片表面CSD节点上的一阶振动模态插值到叶片表面的CFD节点上，作为叶片表面所有CFD节点在流场中的振动位移。小范围选点插值保证了处理速度，引入面积控制因子保证插值的精度。三维线性插值方法流程如图6所示。

三维线性插值方法的具体实施步骤如下：

步骤一：寻找距离待插值的流场节点F最近的三个固体域的节点A₁、B₁、C₁。

如图7所示，对叶片表面上的CFD网格节点F，在叶片CSD网格节点中找到与其距离最近的三个节点A₁、B₁、C₁。需要注意的是节点A₁、B₁、C₁不能在一条直线上或接近一条直线。这是因为当A₁、B₁、C₁发生位移以后，三角形面积会发生变化，如果ΔA₁B₁C₁面积几乎为零或者有一个大钝角，会导致ΔA₂B₂C₂面积变化很剧烈，从而导致CFD网格节点发生剧烈的移动，产生较大误差。引入面积控制因子s，限定ΔA₁B₁C₁的面积不小于s，保证A₁、B₁、C₁三个节点不接近于直线。

寻找节点A₁、B₁、C₁时，遍历所有的节点以后，将距离节点F最近的w个节点存储在一个向量存储器中，该向量存储器里储存的这w个节点是按照与节点F的距离远近排列的，储存的第一个节点离节点F最近，储存的第二个节点离节点F次近，……，储存的第w个节点离节点F最远；将储存的第一个节点作为节点A₁，将储存的第二个节点作为节点B₁，再通过引入面积控制因子s从剩余的数组元素中找到节点C₁，确保节点C₁既不与节点A₁、节点B₁在同一直线上，又确保节点C₁为向量存储器中除节点A₁、节点B₁外距离节点F最近的节点。

步骤二：获取三个固体域节点A₁、B₁、C₁组成的平面及该平面的法向向量。

由A₁、B₁、C₁三点坐标(x_A1，y_A1，z_A1)，(x_B1，y_B1，z_B1)，(x_C1，y_C1，z_C1)组成的平面方程为

$\left|\begin{array}{ccc}x-x_{A1}& y-y_{A1}& z-z_{A1}\\ x_{B1}-x_{A1}& y_{B1}-y_{A1}& z_{B1}-z_{A1}\\ x_{Z1}-x_{A1}& y_{Z1}-y_{A1}& z_{Z1}-z_{A1}\end{array}\right|=0---\left(1\right)$

可以化简为

a₁x+b₁y+c₁z+d₁＝0    (2)

式中，

$\left\{\begin{array}{c}{a}_1=(y_{B1}-y_{A1})(z_{C1}-z_{A1})-(y_{C1}-y_{A1})(z_{B1}-z_{A1})\\ b_1=-(x_{B1}-x_{A1})(z_{C1}-z_{A1})+(x_{C1}-x_{A1})(z_{B1}-z_{A1})\\ c_1=(x_{B1}-x_{A1})(y_{C1}-y_{A1})-(x_{C1}-x_{A1})(y_{B1}-y_{A1})\\ d_1=-(a_1{x}_{A1}+b_1{y}_{A1}+c_1{z}_{A1})\end{array}\right.---\left(3\right)$

平面ΔA₁B₁C₁的法向向量为 $\stackrel{\→}{n}=(a_1,b_1,c_1)$

步骤三：根据三个固体域节点发生的位移确定运动后产生的另一平面ΔA₂B₂C₂。

在节点A₁、B₁、C₁上分别作ΔA₁B₁C₁的法线，使得A₁A₂、B₁B₂、C₁C₂分别为节点A₁、B₁、C₁的位移值。

已知A₁、B₁、C₁三节点函数值的变化，也即A₁A₂、B₁B₂和C₁C₂的长度，并且A₁A₂、B₁B₂和C₁C₂的向量与 $\stackrel{\→}{n}=(a_1,b_1,c_1)$ 相同，由此可以根据余弦公式获得A₂、B₂、C₂的坐标。

$\left\{\begin{array}{c}\cos x=\frac{a_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\\ \cos y=\frac{b_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\\ \cos z=\frac{c_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\end{array}\right.---\left(4\right)$

A₂、B₂和C₂的坐标分别为：

$\left\{\begin{array}{c}{x}_{A2}=x_{A1}+A_1{A}_2\cos x\\ y_{A2}=y_{A1}+A_1{A}_2\cos y\\ z_{A2}=z_{A1}+A_1{A}_2\cos z\end{array}\right.,\left\{\begin{array}{c}{x}_{B2}=x_{B1}+B_1{B}_2\cos x\\ y_{B2}=y_{B1}+B_1{B}_2\cos y\\ z_{B2}=z_{B1}+B_1{B}_2\cos z\end{array}\right.,\left\{\begin{array}{c}{x}_{C2}=x_{C1}+C_1{C}_2\cos x\\ y_{C2}=y_{C1}+C_1{C}_2\cos y\\ z_{C2}=z_{C1}+C_1{C}_2\cos z\end{array}\right.---\left(5\right)$

平面ΔA₂B₂C₂的方程为

a₂x+b₂y+c₂z+d₂＝0    (6)

其中 $\left\{\begin{array}{c}{a}_2=(y_{B2}-y_{A2})(z_{C2}-z_{A2})-(y_{C2}-y_{A2})(z_{B2}-z_{A2})\\ b_2=-(x_{B2}-x_{A2})(z_{C2}-z_{A2})+(z_{C2}-x_{A2})(z_{B2}-z_{A2})\\ c_2=(x_{B2}-x_{A2})(y_{C2}-y_{A2})-(x_{C2}-x_{A2})(y_{B2}-y_{A2})\\ d_2=-(a_2{x}_{A2}+b_2{y}_{A2}+c_2{z}_{A2})\end{array}\right.$

步骤四：获取插值得到的节点F的位移。

过节点F做ΔA₁B₁C₁的法线，在ΔA₁B₁C₁和ΔA₂B₂C₂上的交点分别为F₁、F₂，则F₁F₂为节点F的位移值，

为ΔA₁B₁C₁的外法向向量，F₁F₂的正负由

决定：

与

同向时，F₁F₂为正；

与

反向时，F₁F₂为负；

直线F₁F₂的参数方程为

$\left\{\begin{array}{c}x=x_F-a_1{t}_1\\ y=y_F-b_1{t}_1\\ z=z_F-c_1{t}_1\end{array}\right.---\left(7\right)$

获得F₁F₂与ΔA₁B₁C₁的交点F₁(x_F1，y_F1，z_F1)的坐标为：

$\left\{\begin{array}{c}{t}_1=\frac{a_1{x}_F+b_1{y}_F+c_1{z}_F+d_1}{a_1^2+b_1^2+c_1^2}\\ x_{F1}=x_F-a_1{t}_1\\ y_{F1}=y_F-b_1{t}_1\\ z_{F1}=z_F-c_1{t}_1\end{array}\right.---\left(8\right)$

F₂(x_F2，y_F2，z_F2)的坐标为：

$\left\{\begin{array}{c}{t}_2=\frac{a_2{x}_F+b_2{y}_F+c_2{z}_F+d_2}{a_1{a}_2+b_1{b}_2+c_1{c}_2}\\ x_{F2}=x_F-a_1{t}_2\\ y_{F2}=y_F-b_1{t}_2\\ z_{F2}=z_F-c_1{t}_2\end{array}\right.---\left(9\right)$

节点F的位移值F₁F₂为：

$\left|F_1{F}_2\right|=\sqrt{{(x_{F1}-x_{F2})}^2+{(y_{F1}-y_{F2})}^2+{(z_{F1}-z_{F2})}^2}---\left(10\right)$

多层动网格方法流程如图8所示。通过多层动网格方法，将叶片表面各个CFD节点的振动位移按一定的比例分配到叶片周围M层O型网格对应的节点上，由叶片表面向外法向递减，第M层位移为零。以叶片的初始位置作为第一个网格文件，叶片按正弦规律振动，每个周期离散为60个时间步。多层动网格生成的网格文件，给出了叶片在一个周期平均的60个时刻对应的叶片表面节点的空间位置。

多层动网格方法的实施步骤如下：

步骤一：找到叶片表面的M层流体域O型网格节点；

如图9所示，对于O型网格，每个流体单元i中的各个节点都按照如图9所示的顺序排列，1，2，6，5节点分别与4，3，7，8节点对应，而每一层面单元的节点排列顺序为顺时针，即4-3-7-8。由叶片表面的第一层网格节点找到第一层流体单元，通过上述的对应关系找到与结构体表面节点对应的第二层节点。依此类推，找出M层O型网格节点。

步骤二：获得各层网格节点与叶片表面对应点的距离及各层对应点位移比例；

叶片表面节点，即第1层网格节点，其最大运动位移即为插值获得的叶片表面各CFD节点振幅，第M层节点最大运动位移为零；各层按照初始的距离比例分配位移；初始未运动时，第J层网格与第M层网格对应点之间的距离δ_j为：

${\mathrm{\δ}}_j=\sqrt{{(x_m-x_j)}^2+{(y_m-y_j)}^2+{(z_m-z_j)}^2},(j=1,...,m-1)---\left(11\right)$

运动开始以后，各层的位移峰值按照初始的距离比例系数分配，则有

$s_j=\frac{{\mathrm{\δ}}_j}{{\mathrm{\δ}}_1},(j=1,...,m-1)---\left(12\right)$

其中s_j为位移比例系数。

步骤三：各层网格节点按比例位移移动，获得运动位移的峰值；

按照步骤二中获得的位移比例系数分配各层的位移：

d_j＝s_jD，(j＝1，...，m-1)    (13)

式中d_j为第J层网格节点的最大位移，D为插值后叶片表面各CFD节点的最大位移。

步骤四：所有网格节点按正弦运动规律输出网格文件；

令M层网格节点按照正弦运动，一个周期离散为N个时间步，生成并输出网格文件。

图10给出了叶尖的流场网格，从图中可以看出，流场网格由O型网格和H型网格构成，两者有一个比较明显的交界面。叶片的运动只能在O型区域内运动，但在O型区内运动的层数可以由用户自己设置，只要保证在运动过程中各层网格节点不产生过大的变形而出现负体积就可以了。

流场中的气动力的获取通过CFX实现。在已建立的流场模型中加载边界条件和初始条件。叶片表面给定运动网格边界条件，运动网格边界条件是通过CFX中的Junction Box模块实现的。Junction box是CFX中的一个接口模块，主要用于实现指定网格区域的运动规律。需要读入的网格文件为叶片一个周期运动所包含的各时刻的网格节点坐标。50％～60％转速下的初始条件和边界条件如表3所示。叶片按照第一阶固有频率振动，其运动周期为频率的倒数，每个时间步间隔相同的时间，在每一个时间步上求解N-S方程。输出结果需要设置成每一个时间步输出气动力和对应的节点位移等信息。

根据每个时间步叶片表面的气动力和节点位移，按照能量法的概念推导出的等效模态气动阻尼公式为 ${\mathrm{\ξ}}_i=-\frac{W_{\mathrm{cfd}}}{2\mathrm{\π}{q}_{i0}^2{\mathrm{\ω}}_i^2},$ 其中W_cfd为第i阶模态下状态空间坐标系中非定常气动力做的功，q_i0为第i阶模态的振幅，ω_i为叶片振动的第i阶固有频率。由此获取叶片一个振动周期内的气动功和叶片等效模态气动阻尼系数。具体的处理中，气动功即为振荡叶栅上非定常气动力在一个振动周期内对叶片所做的功。叶片单位面积上一个周期内的气动功为： $w(\stackrel{\→}{x},t)={\∫}_0^{t_0+T_{\mathrm{blade}}}p(\stackrel{\→}{x},t)\stackrel{\→}{V}(\stackrel{\→}{x},t)\·\stackrel{\→}{n}(\stackrel{\→}{x},t)\mathrm{dt},\∀x\∈\∂B,$ 式中，

为叶片表面点处的静压，为叶片表面点处的速度矢量，为叶片表面在点

处的单位外法向向量，叶片表面，T_blade为叶片的振动周期。整个叶片的气动功为叶片单位面积上一个周期内气动功在叶片表面上的积分： $W_d^{\mathrm{cfd}}=\underset{\∂}{\∫}\underset{B}{\∫}w\left(\stackrel{\→}{x}\right)\mathrm{dS}.$ dS表示在叶片表面积分。

表3  NASA67叶片各级转速下不同反压下边界条件和初始条件

表3给出了各个工况的初始条件和边界条件，各转速的出口压力包含了整个叶片的工作区域，需要注意的是由于转速不同，叶片的第一阶动频不同，因此，所采用的时间步长需要根据各个转速下的第一阶动频f来确定。叶片表面给定动网格的运动模式，由多层动网格方法生成的各层网格节点在各时刻的坐标来控制网格节点的运动。

获得不同工作状况下的等效模态气动阻尼，在图11上标注出叶片的等转速线以及各等转速线上对应的等效模态阻尼比。流量-压比图上的等值线云图就是所要求的风扇叶片特性图中的等阻尼分布。等效模态阻尼比与流量之间的关系因各级转速的不同而有不同的规律。当转速比较高时，如图12所示100％n、90％n和80％n的等转速线上，等效模态阻尼随流量的增加先增加后减小，变化中间存在一个峰值，整个等转速线上的等效模态阻尼比比较大，其中n表示发动机的工作转速，是一个确定的数值。本例中n为16043r/min。而70％n这一条等转速线可以视为一个过渡，等效模态阻尼比的变化规律从这个转速开始逐渐变得平缓，随着流量增加等效模态气动阻尼比逐渐增加，整个变化过程中出现了拐点，最大值出现在靠近喘振边界的一段。50％n和60％n这两条等转速线比较一致，等效模态阻尼比的变化规律随流量增加变化并不明显，一直呈上升趋势，最大值出现在流量较大的一段。图12给出了比较清晰的流量-压比图上的等值线云图。从图12中可以看出叶片在各工作转速中，通过上述方法得到的阻尼为零的等值线超出了叶片特性图中的失速边界，风扇中没有出现颤振边界。这与实验值的结果是一致的。在本例中，没有考虑叶尖间隙和流场动网格的不确定因素。

对于其他的压气机叶片，可以通过以上的方法，得到压气机特性图上的等阻尼线，以阻尼值为零的线作为颤振边界。

Claims (3)

1、一种航空涡轮发动机叶片颤振边界的模拟方法，其特征在于，该方法包括如下步骤：

步骤一：建立叶片的实体模型；

建立需要分析的发动机风扇叶片的三维实体模型；

步骤二：对叶片进行模态分析；

首先，将已建立的发动机风扇叶片的三维实体模型导入有限元分析工具ANSYS中，并对其进行有限元单元划分，定义材料和单元属性，加上载荷和约束条件；

然后，利用有限元分析工具ANSYS中的动模态分析获得叶片的固有模态，提取流场分析所需要的频率和模态，导出叶片表面的计算结构动力学CSD有限元网格节点信息，得到CSD节点；

步骤三：建立叶片的流场模型；

首先，将已建立的发动机风扇叶片的三维实体模型导入流体动力学仿真工具CFX中，建立发动机叶片单通道流场模型，划分流体网格；

然后，将叶片表面的流体网格节点和单元信息导出，得到计算流体动力学CFD节点；

步骤四：获得线性插值叶片的振动位移；

采用三维线性插值方法，将叶片表面CSD节点上的第阶振动模态插值到叶片表面的CFD节点上，作为叶片表面所有CFD节点在流场中的振动位移，其中P为正整数，且P∈[1，10]；小范围选点插值能够保证处理速度，引入面积控制因子保证插值的精度；

步骤五：多层动网格化处理生成流场分析需要的网格文件；

采用多层动网格方法，将叶片表面各个CFD节点的振动位移按各层初始的距离比例，分配到叶片周围M层O型网格区域对应的节点上，由叶片表面向外法向递减，第M层位移为零，最后生成各时刻流场所有节点的网格文件；

所述多层动网格方法适用于1～M层，其中M为流场模型中的O型网格的层数；

步骤六：调用动网格模块获得流场中各参数；

流场分析通过流体动力学仿真工具CFX实现；

首先，在步骤三建立的流场模型中加载边界条件和初始条件；通道入口给定静温和总压，出口给定静压，循环对称面给定交界面边界条件，叶片表面给定运动网格边界条件；以不加运动网格的定常场作为初始条件；

然后，叶片按照第P阶固有模态振动，其运动周期为频率的倒数，每个时间步间隔相同的时间，在每一个时间步上求解N-S方程；

最后，输出结果设置成每一个时间步输出气动力和对应的节点位移信息；

其中时间步是指叶片一个运动周期的N分之一，N是用户设定的叶片一个周期运动的时间步数，实际应用中为了保证模拟效率，取N为正整数，且N∈[30，80]；

步骤七：获得每个工况下的气动阻尼；

根据每个时间步叶片表面的气动功和节点位移，按照能量法的概念推导出的等效模态阻尼公式来获得叶片一个振动周期内的气动功和等效模态气动阻尼系数；

其中气动功是按照各节点在各时间步上做的功积分获得，积分区间为一个振动周期，积分域为叶片表面；

其中能量法的概念是指根据叶片在一个振动周期内，从外界所得到的能量的正负来判别颤振发作与否，当从外界所得到的能量为正时，则颤振发作；当从外界所得到的能量为负时，则颤振不发作；

按照能量法的概念推导出的等效模态阻尼公式是 ${\mathrm{\ξ}}_i=-\frac{W_{\mathrm{cfd}}}{2\mathrm{\π}{q}_{i0}^2{\mathrm{\ω}}_i^2},$ 其中W_cfd为第i阶模态下状态空间坐标系中非定常气动力做的功，q_i0为第i阶模态的振幅，ω_i为叶片振动的第i阶固有频率；

步骤八：获得颤振边界；

获得不同工作状况下的等效模态气动阻尼，获得工作域中等效模态气动阻尼的分布，确定阻尼值为零的等值线为颤振边界。

2、根据权利要求1所述的一种航空涡轮发动机叶片颤振边界的模拟方法，其特征在于：步骤四中所述的三维线性插值方法的具体实施步骤如下：

步骤一：寻找距离待插值的流场节点F最近的三个固体域的节点A₁、B₁、C₁；

对叶片表面上的CFD网格节点F，在叶片CSD网格节点中找到与其距离最近的三个节点A₁、B₁、C₁，其中节点A₁、B₁、C₁不能在一条直线上或接近一条直线；引入面积控制因子s，限定ΔA₁B₁C₁的面积不小于s，保证A₁、B₁、C₁三个节点不接近于直线；

寻找节点A₁、B₁、C₁时，遍历所有的节点以后，将距离节点F最近的w个节点存储在一个向量存储器中，该向量存储器里储存的这w个节点是按照与节点F的距离远近排列的，储存的第一个节点离节点F最近，储存的第二个节点离节点F次近，……，储存的第w个节点离节点F最远；将储存的第一个节点作为节点A₁，将储存的第二个节点作为节点B₁，再通过引入面积控制因子s从剩余的数组元素中找到节点C₁，确保节点C₁既不与节点A₁、节点B₁在同一直线上，又确保节点C₁为向量存储器中除节点A₁、节点B₁外距离节点F最近的节点；

步骤二：获取三个固体域节点A₁、B₁、C₁组成的平面及该平面的法向向量；

由A₁、B₁、C₁三个节点坐标(x_A1，y_A1，z_A1)，(x_B1，y_B1，z_B1)，(x_C1，y_C1，z_C1)组成的平面方程为：

$\left|\begin{array}{ccc}x-x_{A1}& y-y_{A1}& z-z_{A1}\\ x_{B1}-x_{A1}& y_{B1}-y_{A1}& z_{B1}-z_{A1}\\ x_{Z1}-x_{A1}& y_{Z1}-y_{A1}& z_{Z1}-z_{A1}\end{array}\right|=0$

化简为：

a₁x+b₁y+c₁z+d₁＝0

式中，

$\left\{\begin{array}{c}{a}_1=(y_{B1}-y_{A1})(z_{C1}-z_{A1})-(y_{C1}-y_{A1})(z_{B1}-z_{A1})\\ b_1=-(x_{B1}-x_{A1})(z_{C1}-z_{A1})+(x_{C1}-x_{A1})(z_{B1}-z_{A1})\\ c_1=(x_{B1}-x_{A1})(y_{C1}-y_{A1})-(x_{C1}-x_{A1})(y_{B1}-y_{A1})\\ d_1=-(a_1{x}_{A1}+b_1{y}_{A1}+c_1{z}_{A1})\end{array}\right.$

平面ΔA₁B₁C₁的法向向量为 $\stackrel{\→}{n}=(a_1,b_1,c_1);$

步骤三：根据三个固体域节点发生的位移确定运动后产生的另一平面ΔA₂B₂C₂；

在节点A₁、B₁、C₁上分别作ΔA₁B₁C₁的法线，使得A₁A₂、B₁B₂、C₁C₂分别为节点A₁、B₁、C₁的位移值；

已知A₁、B₁、C₁三节点函数值的变化，也即A₁A₂、B₁B₂和C₁C₂的长度，并且A₁A₂、B₁B₂和C₁C₂的向量与 $\stackrel{\→}{n}=(a_1,b_1,c_1)$ 相同，由此根据余弦公式获得A₂、B₂、C₂的坐标；

$\left\{\begin{array}{c}\cos x=\frac{a_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\\ \cos y=\frac{b_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\\ \cos z=\frac{c_1}{\sqrt{a_1^2+b_1^2+c_1^2}}\end{array}\right.$

A₂、B₂和C₂的坐标分别为：

$\left\{\begin{array}{c}{x}_{A2}=x_{A1}+A_1{A}_2\cos x\\ y_{A2}=y_{A1}+A_1{A}_2\cos y\\ z_{A2}=z_{A1}+A_1{A}_2\cos z\end{array}\right.,$ $\left\{\begin{array}{c}{x}_{B2}=x_{B1}+B_1{B}_2\cos x\\ y_{B2}=y_{B1}+B_1{B}_2\cos y\\ z_{B2}=z_{B1}+B_1{B}_2\cos z\end{array}\right.,$ $\left\{\begin{array}{c}{x}_{C2}=x_{C1}+C_1{C}_2\cos x\\ y_{C2}=y_{C1}+C_1{C}_2\cos y\\ z_{C2}=z_{C1}+C_1{C}_2\cos z\end{array}\right.$

平面ΔA₂B₂C₂的方程为：

a₂x+b₂y+c₂z+d₂＝0

其中 $\left\{\begin{array}{c}{a}_2=(y_{B2}-y_{A2})(z_{C2}-z_{A2})-(y_{C2}-y_{A2})(z_{B2}-z_{A2})\\ b_2=-(x_{B2}-x_{A2})(z_{C2}-z_{A2})+(x_{C2}-x_{A2})(z_{B2}-z_{A2})\\ c_2=(x_{B2}-x_{A2})(y_{C2}-y_{A2})-(x_{C2}-x_{A2})(y_{B2}-y_{A2})\\ d_2=-(a_2{x}_{A2}+b_2{y}_{A2}+c_2{z}_{A2})\end{array}\right.$

步骤四：获取插值得到的节点F的位移；

过节点F做ΔA₁B₁C₁的法线，在ΔA₁B₁C₁和ΔA₂B₂C₂上的交点分别为F₁、F₂，则F₁F₂为节点F的位移值，

为ΔA₁B₁C₁的外法向向量，F₁F₂的正负由

决定：

与同向时，F₁F₂为正；

与

反向时，F₁F₂为负；

直线F₁F₂的参数方程为：

$\left\{\begin{array}{c}x=x_F-a_1{t}_1\\ y=y_F-b_1{t}_1\\ z=z_F-c_1{t}_1\end{array}\right.$

获得F₁F₂与ΔA₁B₁C₁的交点F₁(x_F1，y_F1，z_F1)的坐标为：

$\left\{\begin{array}{c}{t}_1=\frac{a_1{x}_F+b_1{y}_F+c_1{z}_F+d_1}{a_1^2+b_1^2+c_1^2}\\ x_{F1}=x_F-a_1{t}_1\\ y_{F1}=y_F-b_1{t}_1\\ z_{F1}=z_F-c_1{t}_1\end{array}\right.$

F₂(x_F2，y_F2，z_F2)的坐标为：

$\left\{\begin{array}{c}{t}_2=\frac{a_2{x}_F+b_2{y}_F+c_2{z}_F+d_2}{a_1{a}_2+b_1{b}_2+c_1{c}_2}\\ x_{F2}=x_F-a_1{t}_2\\ y_{F2}=y_F-b_1{t}_2\\ z_{F2}=z_F-c_1{t}_2\end{array}\right.$

节点F的位移值F₁F₂为：

$\left|F_1{F}_2\right|=\sqrt{{(x_{F1}-x_{F2})}^2+{(y_{F1}-y_{F2})}^2+{(z_{F1}-z_{F2})}^2}.$

3、根据权利要求1所述的一种航空涡轮发动机叶片颤振边界的模拟方法，其特征在于：步骤五中所述的多层动网格方法的实施步骤如下：

步骤一：找到叶片表面的M层流体域O型网格节点；

对于O型网格，每个流体单元i中的1，2，6，5节点分别与4，3，7，8节点对应，而每一层面单元的节点排列顺序为顺时针，即4-3-7-8；由叶片表面的第一层网格节点找到第一层流体单元，通过上述的对应关系找到与结构体表面节点对应的第二层节点；依此类推，找出M层O型网格节点；

步骤二：获得各层网格节点与叶片表面对应点的距离及各层对应点位移比例；

叶片表面节点，即第1层网格节点，其最大运动位移即为插值获得的叶片表面各CFD节点振幅，第M层节点最大运动位移为零；各层按照初始的距离比例分配位移；初始未运动时，第J层网格与第M层网格对应点之间的距离δ_j为：

${\mathrm{\δ}}_j=\sqrt{{(x_m-x_j)}^2+{(y_m-y_j)}^2+{(z_m-z_j)}^2},(j=1,...,m-1)$

运动开始以后，各层的位移峰值按照初始的距离比例分配，则有：

$s_j=\frac{{\mathrm{\δ}}_j}{{\mathrm{\δ}}_1},(j=1,...,m-1)$

其中s_j为位移比例系数；

步骤三：各层网格节点按比例位移移动，获得运动位移的峰值；

按照步骤二中获得的位移比例系数分配各层的位移：

d_j＝s_jD，(j＝1，...，m-1)

式中d_j为第J层网格节点的最大位移，D为插值后叶片表面各CFD节点的最大位移；

步骤四：所有网格节点按正弦运动规律输出网格文件；

令M层网格节点按照正弦运动，一个周期离散为N个时间步，生成并输出网格文件。

Images