CN112182802B - 计入随机失谐的叶轮机械气动弹性优化设计实现方法 - Google Patents

Abstract

一种计入随机失谐的叶轮机械气动弹性优化设计实现方法，基于影响系数法的颤振计算、基础失谐模型和蒙特卡洛法的强迫振动分析，在约束气动性能以及静应力的前提下，使工业加工的叶片在不同的安装方式下具有相似的振动特性；显著提升了颤振稳定性并降低了在随机失谐的情况下叶片强迫振动响应，本方法对于选取的轴流风扇叶片，在约束气动性能和静应力的前提下，其最小的气动阻尼增加了约19％和随机失谐下最大的强迫振动幅值降低了约18％。

Description

计入随机失谐的叶轮机械气动弹性优化设计实现方法

技术领域

本发明涉及的是一种航空发动机制造领域的技术，具体是一种计入随机失谐的叶轮机械气动弹性优化设计实现方法。

背景技术

叶片作为叶轮机械中的做功部件，旋转中的叶片受到各种载荷(如旋转产生的离心力，上游静叶/导叶尾迹产生的非稳态气动力等等)，不仅会影响气动性能还会引起叶片振动问题。因此，为了保证高效且安全的运行要求，叶片的设计必须兼顾气动及结构性能。

叶片和流体的相互作用而导致的气动弹性振动会导致叶片疲劳失效，其中最常见的两种气动弹性现象是：1.颤振；2.强迫振动。颤振是一种自激现象，叶片不断地从流场中吸收能量，导致叶片振动幅值不断放大；强迫振动是由于外界激励而导致叶片振动，振动的频率通常与转速频率成整数倍关系。这两种振动形式都会使叶片所受应力增加，在有限个振动周期内造成叶片疲劳失效：轻则叶片产生裂纹，重则叶片断裂，造成严重事故。通常叶片受到的应力大小与其振动的振幅成正比。叶片振动的幅值除了与材料属性，激励的大小，摩擦阻尼等因素相关之外，制造加工误差也能对其产生较大的影响。

发明内容

本发明针对现有技术的空白，提出一种计入随机失谐的叶轮机械气动弹性优化设计实现方法，基于影响系数法的颤振计算、基础失谐模型和蒙特卡洛法的强迫振动分析，在约束气动性能以及静应力的前提下，显著提升了颤振稳定性并降低了在随机失谐的情况下叶片强迫振动响应，为叶片设计提供了一个新途径，适用于轴流及离心式叶轮机械，具有良好的实用价值和广阔的应用前景。

本发明是通过以下技术方案实现的：

本发明涉及一种计入随机失谐的叶轮机械气动弹性优化设计实现方法，包括以下步骤：

步骤一：对原始叶片进行参数化处理：提取各截面翼型的弦长、中弧线、前缘与尾缘的半径和安装角，再分别定义构叶片弯、扭和掠的形变方向及大小，完成叶片的三维造型设计。

步骤二：利用拉丁超立方采样，跟据叶片设计变量的数量生成满足精度要求的实验样本。

步骤三：分别建立单通道流场模型、全通道流场模型、单个叶片模型及整体叶盘模型，相应地得到：单通道流场网格、全通道流场网格、叶片有限元网格、整体叶盘有限元网格。

步骤四：利用步骤三得到的单通道流场网格计算气动性能、压比(η)和效率(π)，并将叶片表面压力以文件形式输出；将来自单通道流场网格计算到的叶片压力插值到叶片有限元网格中作为叶片所受气动力，将由叶片旋转产生的离心力和气动力作为预应力进行静力学分析，得到叶片的静应力σ_s。

步骤五：在静力学分析的基础上进行单个叶片的模态分析，得到叶片质量归一化的模态振型

最大振幅(x_fem)和固有频率(f₀)。

步骤六：将步骤五中得到的模态振型插值到全通道流场网格，采用动网格的方式进行叶片振动模拟；

由于采用的是基于线性假设的影响系数法，所以只需要一个叶片振动即可，叶片振动的幅值(x_cfd)一般小于叶尖弦长的1％；叶片振动频率为其固有频率，计算得到叶片所受到的非定常气动力。

步骤七：根据模态力定义，将步骤六得到的非定常气动力投影到全通道流场模型的叶片有限元网格中的模态坐标得到模态力(Q)，根据模态力构建影响系数矩阵(C)，然后采用特征值法求得谐调状态下的气动阻尼(ζ)并取最小气动阻尼(ζ_min)为优化目标。

步骤八：根据步骤三建立的整体叶盘模型，进行模态分析计算得到叶盘各模态簇的固有频率与节径的关系；将步骤六得到的模态力与叶盘固有频率及其对应的节径导入基础失谐模型中，得到谐调状态下的强迫振动响应。

所述的导入，优选采用扫频的方式，在一定的频率范围内计算叶片的动态响应；频率范围由叶片固有频率确定。

步骤九：利用蒙特卡洛法，生成1000组服从正态分布的具有相同失谐量但不同形式的随机失谐方案并利用基础失谐模型分别按步骤一～步骤八计算每个叶片设计在1000组随机失谐的强迫振动响应；采用统计的方法得到叶片归一化振动幅值的正态分布概率密度函数；以其中位置参数μ、尺度参数σ作为优化目标。

步骤十：根据步骤二中生成的实验样本，经步骤七和九的计算得到的对应的目标函数(ζ_min，μ和σ)，采用响应面法(Response surface methodology，RSM)分别构建叶片设计变量(v₁，v₂，v₃，…，v_n)与每个优化目标的关系：ζ_min＝f(v₁，v₂，v₃，…，v_n)，μ＝g(v₁，v₂，v₃，…，v_n)，σ＝h(v₁，v₂，v₃，…，v_n)。

步骤十一：在约束条件下，采用多目标优化算法，得到帕累托前沿，并选择若干预测的帕累托解，即预测的较优的叶片设计。

所述的预测的帕累托解的数量优选为10个。

所述的约束条件包括：η₀≤η，π₀≤π，σ_s≤σ_s0，其中：η₀和π₀分别表示原始叶片的效率以及压比；σ_s0为原始叶片的静应力；η和π分别表示优化后叶片的效率以及压比；σ_s为优化后叶片的静应力。

步骤十二：重复步骤三到步骤九，计算若干预测的叶片设计的效率、压比、静应力、最小气动阻尼和正态分布的位置参数和尺度参数。

所述的预测的叶片的数量优选为10个。

步骤十三：验证若干预测的叶片设计是否某个设计满足约束条件及优化目标；当满足约束条件及优化目标，则输出优化后叶片的参数，气动性能，静应力，颤振分析结果及强迫振动分析结果；否则将预测的设计带入到响应面模型中，重新生成新的响应面模型

然后再采用多目标优化算法预测的更新叶片设计，重复步骤三到十三。

所述的优化目标包括：ζ_min0≤ζ_min，μ≤μ₀，σ≤σ₀，其中：ζ_min0，μ₀和σ₀分别表示原始叶片的最小的气动阻尼，正态分布的位置参数和尺度参数；ζ_min，μ和σ分别表示优化后叶片的最小的气动阻尼，正态分布的位置参数和尺度参数。

技术效果

本发明整体解决了现有的叶片由于气动弹性振动导致的疲劳失效的问题。

与现有技术相比，本发明在保证气动性能的前提下，降低了叶片强迫振动对机械加工误差的敏感性，使工业加工的叶片在不同的安装方式下具有相似的振动特性；本发明降低了叶片由机械加工误差导致的强迫振动幅值，从而提升叶片的寿命。

附图说明

图1为本发明流程示意图；

图2为流场网格划分-单通道示意图；

图3为有限元分析网格划分示意图；

图4为一阶弯曲模态示意图；

图5为模态振动位移插值到流体网格示意图；

图6为相邻5个叶片的模态力随计算步数的变化示意图；

图7为相邻5个叶片的模态力幅值示意图；

图8为不同叶高处的翼型比较示意图；

图中：原始叶片与优化后叶片分别在20％，50％，80％叶高处的翼型截面示意图；

图9a、b为气动性能参数对比示意图；

图中：原始叶片与优化后叶片气动性能参数的比较。优化后的叶片在压比和效率方面都有较明显的提升示意图；

图10a、b为无量纲静应力分布示意图；

图中：原始叶片与优化后叶片在OP3处的静应力分布。优化后的叶片的静应力有较明显的下降示意图；

图11a、b为颤振特性示意图；

图中：原始叶片与优化后叶片的颤振特性的比较。优化后的叶片的最小气动阻尼明显高于原始叶片的最小气动阻尼，即系统在当前工况相对更稳定示意图；

图12为归一化振动幅值的概率密度函数-2％随机失谐示意图；

图中：在2％的随机失谐量下，叶片强迫振动响应的概率密度函数于归一化振动幅值的关系示意图；优化后的叶片的位置参数和尺度参数小于原始叶片对应的值。从统计学的角度上看，优化后的叶片的强迫振动响应相对有所降低示意图；

图13为实施例系统示意图。

具体实施方式

如图1所示，为本实施例涉及一种计入随机失谐的叶轮机械气动弹性优化设计方法，包括以下步骤：

步骤一：选取某轴流风扇叶片，采用Matlab编写的函数(叶片参数提取单元)对原始叶片进行参数化处理。

步骤二：实验设计：根据叶片设计变量的数量及每个变量的上下限，利用拉丁超立方采样，生成合适的实验数量。

步骤三：采用自定义Matlab函数驱动计算流体力学软件CFX中的Turbogrid模块进行流场单通道网格划分，如图2所示。叶片表面为O型网格，且网格尺度满足Y+小于5。将单通道流场网格导入CFX-Pre中，利用Turbo Rotation以旋转轴复制(N-1)份网格即可得到全通道流场网格，其中：N为叶片的数量。采用Matlab驱动Workbench中Static Structural模块进行有限元网格划分，如图3对单个叶片进行网格划分。通常，有限元计算所需网格数量远小于流场计算的网格。同样地，将整体叶盘模型导入Static Structural模块，可以对整体叶盘进行网格划分。

步骤四：采用自定义Matlab函数驱动计算流体力学软件CFX进行稳态计算，具体包括：

4.1)将步骤三中的单通道流场网格导入CFX-pre，分析类型选择Steady State(稳态求解)，湍流模型为K-Epsilon，入口给定总温总压，出口给定质量流量，求得气动性能、压比(η)、效率(π)，并输出叶片表面压力文件。

4.2)为了保证在一定范围内具有较好的气动性能，选取3个工况点(Operatingpoint，OP)进行计算：OP1靠近堵塞点，OP3靠近最高效率点，OP2介于两者之间。

4.3)采用自定义Matlab函数驱动Static Structural模块自动设置材料属性，定义边界条件，转子转速，并将OP3的叶片压力插值到叶片表面，作为叶片所受气动力，进行静力学分析，得到叶片的静应力σ_s。

步骤五：在静力学分析的基础上，利用自定义Matlab函数驱动Workbench中的Modal模块对单叶片进行模态分析，得到叶片的模态振型，并导出质量归一化的模态振型。图4显示了叶片的一阶弯曲模态。

步骤六：将步骤三得到的全通道流场网格导入CFX-Pre中。因为影响系数法要求一个叶片振动即可，所以将叶片质量归一化的模态振型插值到某一个叶片的流场网格，如图5所示，采用动网格技术使其作简谐振动。振动的幅值(x_cfd)通常取小于叶尖弦长的1％，以此满足线性假设。采用CFX求解器进行非稳态计算，分析类型选择Transient blade row，湍流模型为K-Epsilon。非稳态计算的时间步长(τ)与振动频率(f0)相关，

其中：n为一个振动周期内计算步数。以步骤四中OP3的稳态计算结果作为初始解进行计算，得到叶片所受到的非定常气动力。

步骤七：根据模态力(Q)的定义：

其中：A为叶片表面积，

为单元面法向量，p为叶片表面压力，

表示单元面所受到的气动力。因此为了计算模态力，需要在CFX-Pre创建模态力变量表达式：-Pressure*(modeshape1 X*Normal X+modeshape1 Y*Normal Y+modeshape1 Z*Normal Z)，其中：Pressure对应叶片表面压力p，modalshape1 xY Z分别表示一阶模态振型在x，Y，Z三个方向的分量，Normal xY Z分别表示单元法向量在x，Y，Z三个方向的分量。

如图6所示，展示了原始叶片5个明显的模态力随计算步数的变化，因为0号叶片在振动，所以它的模态力最大，离0号叶片越远，模态力越小。取最后一个周期作快速傅里叶变换(FFT)：q_i＝FFT(Q_i)，其中：Q_i表示第i个叶片的模态力。

各个叶片的模态力幅值如图7所示。

利用特征值法进行颤振特性计算，即求解方程：

其中：叶片固有圆频率(单位rad)ω₀＝f₀*2π，I为单位矩阵，λ包含颤振频率(f_c)和气动阻尼信息(ζ)，

z包含各个节径(ND)下每个叶片的振动幅值即相位信息，

，其中：α为相邻两叶片间的相位角之差，单位rad，C影响系数矩阵，是由模态力构成的N*N的循环对称矩阵：

其中：

上述方程中影响系数矩阵和叶片固有圆频率为已知条件，求解可得到叶片的颤振频率及其对应的气动阻尼。

步骤八：利用Matlab自定义函数将叶盘模型导入Workbench的DesignModeler中，并联接Static Structural和modal模块。为了节约计算，可以采用扇形模型，利用循环边界条件求得叶盘振动特性，导出第一阶模态簇节径与其对应的固有频率。

步骤九：利用蒙特卡洛法生成1000组服从正态分布的相同失谐量(ε，2％)，不同的形式的随机失谐，失谐量

其中：f_bi表示第i号叶片的固有频率，采用基础失谐模型求解振动响应：

其中：A为失谐矩阵，ω为叶片振动频率，F_e为由外部激励引起的模态力向量；最后利用统计学方法，求概率密度函数的位置参数μ和尺度参数σ。

步骤十：计算完所有的实验设计之后，利用响应面法分别构建叶片设计变量(v₁，v₂，v₃，…，v_n)与每个目标函数(ζ_min，μ，σ)的关系：ζ_min＝f(v₁，v₂，v₃，…，v_n)，μ＝g(v₁，v₂，v₃，…，v_n)，σ＝h(v₁，v₂，v₃，…，v_n)。

步骤十一：在约束条件下，采用多目标优化算法，得到帕累托前沿，并选择10个预测的帕累托解，即10个预测的较优的叶片设计。

所述的约束条件包括：η₀₁≤η₁，η₀₂≤η₂，η₀₃≤η₃，π₀₁≤π₁，π₀₂≤π₂，π₀₃≤π₃，σ_s≤σ_s0，其中：η₀₁，η₀₂，η₀₃，π₀₁，π₀₂，π₀₃表示原始叶片分别在OP1，OP2，和OP3处的的效率，压比；σ_s0为原始叶片在OP3处的静应力，其中：η₁，η₂，η₃，π₁，π₂，π₃表示优化后叶片分别在OP1，OP2，和OP3处的效率，压比；σ_s为优化后叶片在OP3处的静应力。

步骤十二：重复步骤三到九，计算10个预测的叶片设计的效率，压比，静应力，最小气动阻尼和正态分布的位置参数和尺度参数。

步骤十三：验证叶片设计是否满足约束条件及优化目标：ζ_min0≤ζ_min，μ≤μ₀，σ≤σ₀，其中：ζ_min0，μ₀和σ₀分别表示原始叶片在OP3处的最小的气动阻尼，正态分布的位置参数和尺度参数；ζ_min，μ和σ分别表示优化后叶片在OP3处的最小的气动阻尼，正态分布的位置参数和尺度参数。如果满足约束条件及优化目标，则输出优化后叶片的参数，气动性能，静应力，颤振结果分析及强迫振动结果分析；如果不满足，则将预测的10个设计带入到响应面模型中，重新生成新的响应面模型

然后再采用多目标优化算法预测10个新的叶片设计，重复步骤三到十三。

如图13所示，为本实施例涉及一种实现上述方法的系统，包括：叶片参数提取单元、实验设计及叶片造型单元、网格划分单元，流场及结构计算单元以及优化单元，其中：叶片参数提取单元与实验设计及叶片造型单元相连并传输叶片的设计变量信息，实验设计及叶片造型单元与网格划分单元相连并传输叶片几何造型信息，网格划分单元与流场及结构计算单元相连并传输网格信息，流场及结构计算单元与优化单元相连并传输流场及结构振动特性的信息。所有单元基于Matlab编写的函数实现。

与现有技术相比，本方法对于选取的轴流风扇叶片，在约束气动性能和静应力的前提下，其最小的气动阻尼增加了约19％和随机失谐下最大的强迫振动幅值降低了约18％。

上述具体实施可由本领域技术人员在不背离本发明原理和宗旨的前提下以不同的方式对其进行局部调整，本发明的保护范围以权利要求书为准且不由上述具体实施所限，在其范围内的各个实现方案均受本发明之约束。

Claims (9)

1.一种计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征在于，包括以下步骤：

步骤一：对原始叶片进行参数化处理：提取各截面翼型的弦长、中弧线、前缘与尾缘的半径和安装角，再分别定义构叶片弯、扭和掠的形变方向及大小，完成叶片的三维造型设计；

步骤二：利用拉丁超立方采样，跟据叶片设计变量的数量生成满足精度要求的实验样本；

步骤三：分别建立单通道流场模型、全通道流场模型、单个叶片模型及整体叶盘模型，相应地得到：单通道流场网格、全通道流场网格、叶片有限元网格、整体叶盘有限元网格；

步骤四：利用步骤三得到的单通道流场网格计算气动性能、压比η和效率π，并将叶片表面压力以文件形式输出；将来自单通道流场网格计算到的叶片压力插值到叶片有限元网格中作为叶片所受气动力，将由叶片旋转产生的离心力和气动力作为预应力进行静力学分析，得到叶片的静应力σ_s；

步骤五：在静力学分析的基础上进行单个叶片的模态分析，得到叶片质量归一化的模态振型

最大振幅x_fem和固有频率f₀；

步骤六：将步骤五中得到的模态振型插值到全通道流场网格，采用动网格的方式进行叶片振动模拟，得到叶片所受到的非定常气动力；

步骤七：根据模态力定义，将步骤六得到的非定常气动力投影到全通道流场模型的叶片有限元网格中的模态坐标得到模态力Q，根据模态力构建影响系数矩阵C，然后采用特征值法求得谐调状态下的气动阻尼ζ并取最小气动阻尼ζ_min为优化目标；

步骤八：根据步骤三建立的整体叶盘模型，进行模态分析计算得到叶盘各模态簇的固有频率与节径的关系；将步骤六得到的模态力与叶盘固有频率及其对应的节径导入基础失谐模型中，得到谐调状态下的强迫振动响应；

步骤九：利用蒙特卡洛法，生成1000组服从正态分布的具有相同失谐量但不同形式的随机失谐方案并利用基础失谐模型分别按步骤一～步骤八计算每个叶片设计在1000组随机失谐的强迫振动响应；采用统计的方法得到叶片归一化振动幅值的正态分布概率密度函数；以其中位置参数μ、尺度参数σ作为优化目标；

步骤十：根据步骤二中生成的实验样本，经步骤七和九的计算得到的对应的目标函数ζ_min，μ和σ，采用响应面法分别构建叶片设计变量v₁，v₂，v₃，…，v_n与每个优化目标的关系：ζ_min＝f(v₁，v₂，v₃，…，v_n)，μ＝g(v₁，v₂，v₃，…，v_n)，σ＝h(v₁，v₂，v₃，…，v_n)；

步骤十一：在约束条件下，采用多目标优化算法，得到帕累托前沿，并选择若干预测的帕累托解，即预测的较优的叶片设计；

步骤十二：重复步骤三到步骤九，计算若干预测的叶片设计的效率、压比、静应力、最小气动阻尼和正态分布的位置参数和尺度参数；

步骤十三：验证若干预测的叶片设计是否某个设计满足约束条件及优化目标；当满足约束条件及优化目标，则输出优化后叶片的参数，气动性能，静应力，颤振分析结果及强迫振动分析结果；否则将预测的设计带入到响应面模型中，重新生成新的响应面模型

然后再采用多目标优化算法预测的更新叶片设计，重复步骤三到十三。

2.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤四包括：

4.1)将步骤三中的单通道流场网格导入CFX-pre，分析类型选择稳态求解，湍流模型为K-Epsilon，入口给定总温总压，出口给定质量流量，求得气动性能、压比η、效率π，并输出叶片表面压力文件；

4.2)为了保证在一定范围内具有较好的气动性能，选取3个工况点进行计算，即：OP1靠近堵塞点，OP3靠近最高效率点，OP2介于两者之间；

4.3)采用自定义Matlab函数驱动Static Structural模块自动设置材料属性，定义边界条件，转子转速，并将OP3的叶片压力插值到叶片表面，作为叶片所受气动力，进行静力学分析，得到叶片的静应力σ_s。

3.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤六中，由于采用的是基于线性假设的影响系数法，所以只需要一个叶片振动即可，叶片振动的幅值x_cfd小于叶尖弦长的1％；叶片振动频率为其固有频率，计算得到叶片所受到的非定常气动力。

4.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤八中所述的导入，采用扫频的方式计算叶片的动态响应；频率范围由叶片固有频率确定。

5.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤十一中所述的预测的帕累托解的数量为10个。

6.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤十一中所述的约束条件包括：η₀≤η，π₀≤π，σ_s≤σ_s0，其中：η₀和π₀分别表示原始叶片的效率以及压比；σ_s0为原始叶片的静应力；η和π分别表示优化后叶片的效率以及压比；σ_s为优化后叶片的静应力。

7.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，步骤十二中所述的预测的叶片的数量为10个。

8.根据权利要求1所述的计入随机失谐的叶轮机械气动弹性优化设计实现方法，其特征是，所述的优化目标包括：ζ_min0≤ζ_min，μ≤μ₀，σ≤σ₀，其中：ζ_min0，μ₀和σ₀分别表示原始叶片的最小的气动阻尼，正态分布的位置参数和尺度参数；ζ_min，μ和σ分别表示优化后叶片的最小的气动阻尼，正态分布的位置参数和尺度参数。

9.一种实现上述任一权利要求所述方法的系统，其特征在于，包括：叶片参数提取单元、实验设计及叶片造型单元、网格划分单元，流场及结构计算单元以及优化单元，其中：叶片参数提取单元与实验设计及叶片造型单元相连并传输叶片的设计变量信息，实验设计及叶片造型单元与网格划分单元相连并传输叶片几何造型信息，网格划分单元与流场及结构计算单元相连并传输网格信息，流场及结构计算单元与优化单元相连并传输流场及结构振动特性的信息，所有单元基于Matlab编写的函数实现。

Images