一种透平机械叶片的振动应力数值分析方法与流程

本发明属于透平机械叶片的设计与计算领域，具体涉及一种透平机械叶片的振动应力数值分析方法。

背景技术：

透平机械中的关键部件叶片在高压和高转速等恶劣的环境中工作，故实践运行经验表明，透平机械的重大事大多与叶片振动产生的高周疲劳损伤有关，而叶片在气流激振力作用下产生的振动应力是导致叶片产生高周疲劳的主要原因。

随着计算机硬件的飞速发展和高效数值计算方法的广泛应用，在叶片设计阶段即采用精确的数值方法对其振动应力进行计算，一方面大大缩短设计周期并节约试验成本，另一方面能有效减小以往采用理论分析和理论-实验相结合的半经验分析等方法对模型过度简化而引起的偏差，从而获得可靠的分析结果，对透平叶片的设计和制造具有重要的意义。然而现有技术中还没有一种相对比较成熟且计算精度高、速度快的振动应力数值分析方法。

技术实现要素：

本发明的目的在于提供一种透平机械叶片振动应力数值分析方法，以解决上述技术问题，该方法主要用来处理透平机械叶片在受到气流激振力交变载荷作用下的振动应力计算。计算过程中采用非定常计算方法获得叶片表面的压强分布，并通过叶片流体区域网格向固体区域网格压强进行插值以获得准确的气流激振力载荷，提高了振动应力计算精度；在位移结果扩展为应力结果的步骤中将叶片划分为不同扇区，在各个扇区采用CPU+GPU异构并行计算的方式，可以大幅度提升计算速度。

为了实现上述目的，本发明采用如下技术方案：

一种透平机械叶片的振动应力数值分析方法，包括以下步骤：

1)、给定透平机械叶片的三维模型和材料参数与运行参数，在网格剖分软件中分别建立整圈叶片固体区域的FEM模型与周期对称的流体区域的CFD模型；并将有限元模型的节点编号以及节点坐标，单元编号以及单元对应的节点编号输出到文件Solid.fem中；将材料参数和运行参数输出到文件Blade.dat中；

2)、对叶片周期对称模型进行气动分析；利用第1)步中获得的周期对称的流体区域的CFD模型，分析采用RNG k-ε湍流模型以及非平衡近壁模型，对模型设置工质类型，进口静压、温度、湍流度和叶片出口静压；首先进行定常计算，静叶栅和动叶栅采用冻结转子法，然后将定常结果作为初始场进行非定常计算，瞬态转子的方式进行耦合，取动叶经过一个静叶栅间距时间作为一个气动周期ΔT，时间步长取一个周期ΔT的计算获得非定常计算收敛后的一个周期各个时间步的叶片表面的压强场分布，并计算稳态压强场和脉动压强场，将获得的稳态压强场和脉动压强场数据存储到文件Pressure.dat中，数据包括点坐标以及对应的压强大小；

3)、压强分布数据转化；读取第1)步获得的solid.fem文件中的单只叶片固体网格数据，计算叶片表面单元中心的坐标；读取第2步获得文件Pressure.dat中的压力场数据，并读取压力场中节点在CFD模型中的网格数据，并对固体表面中心的位置进行插值，获得这些位置在一个气动周期不同时间步的压强；

4)、整圈叶片有限元模态分析；根据第2)步中获得稳态的压强插值施加在固体网格单元表面上，对整圈叶片施加该第1)步获得转速的离心力载荷，在有限元软件中求解其预应力场，并将整圈叶片有限元网格的单元编号与高斯节点的预应力结果输出到文件GSS.dat，节点编号与其对应的预应力结果输出到文件NSS.dat；然后进行考虑预应力场和旋转软化的模态分析，获得整圈叶片阶数m的固有频率f以及振型[φ]，并输出到文件Mode.dat；

5)、整圈叶片激振力施加；根据第3)步中获得的一个气动周期内不同时间步的固体网格的脉动压强，计算透平叶片表面单元每个节点的受力，获得所有时间步叶片表面的节点力载荷向量，并将一个气动周期内所有时间步的载荷向量压缩输出到文件中；

6)、模态叠加法求解振动位移响应；首先设置计算时间步，读取第5)步获得的载荷向量文件，并结合第2)步获得的固有频率以及振型，利用杜哈美积分计算在脉动气流激振力整圈叶片位移响应；

7)、位移响应结果扩展为应力结果；在CPU和GPU上同时将叶片不同扇区的位移响应转化为应力向量，求解各个单元高斯节点的应力分量，并外插获得节点的应力分量，最后将高斯节点和节点的应力分量转化为主应力以及VonMises等效应力；

8)、振动应力结果提取与考核；提取所有扇区叶片在围带圆角，叶根平台圆角，拉筋圆角考核位置在第2)步中预应力分析结果中的最大VonMises等效应力，并提取这些位置在不同时间步的单元或节点的VonMises等效应力，将二者按照线性形式叠加，绘制振动应力时域曲线。

进一步的，步骤1)中材料参数包括叶片和轮盘的密度、弹性模量、泊松比；运行参数包括叶片的转速、工质类型、进口静压、温度、湍流度和叶片出口静压。

进一步的，第2)步中利用一个气动周期各个时间步的压强场结果计算稳态压强场和脉动压强场，计算方式如下：

2.1)提取一个节点在不同时间步的压强值P_(t)；

2.2)取一个周期内P(t)的时均值作为稳态压强

2.3)计算该点各个时间步的脉动压强

2.4)遍历所有节点，获得稳态压强场和脉动压强场P′_(t)，并储存在文件Pressure.dat。

进一步的，第3)中插值的流程如下：

3.1)计算固体网格中一个叶片单元表面中心点的坐标，在CPU和GPU上多线程并行判断该点是否落在流体区域网格的四边形内，记下流体区域网格四边形的四个节点；

3.2)对这四个节点组成的四边形进行剖分：假设四个节点坐标分别为(x1,y1)，(x2,y2)，(x3,y3)，(x4,y4)，在四条边分别均匀地插入9个节点，对应节点连接，将该四边形剖分为10*10的小四边形区域；对经度i和纬度j的小四边形的左下角节点坐标为

依次类推，计算所有剖分节点坐标；

3.3)根据剖分获得的所有小四边形的节点坐标，继续在CPU和GPU上同时判断该点是否落在这10*10的小四边形内，以落在其中的小四边形为新的四边形然后返回第3.2)步继续剖分，如此重复4次，记下这4次剖分过程中，落在其中的四边形的经度(i₁,i₂,i₃,i₄)和纬度(j₁,j₂,j₃,j₄)；

3.4)计算插值坐标

取ξ₁＝-1,ξ₂＝1,ξ₃＝1,ξ₄＝-1，η₁＝-1,η₂＝-1,η₃＝1,η₄＝1；

按照下式进行插值，获得该中心点的压强：

3.5)遍历所有固体网格中单元表面中心点，重复3.1)到3.4)步，获得所有单元中心点的压强；

3.6)按照3.1)—3.5)计算一个周期所有时间步压强场在所有固体网格单元表面中心点的插值，包括第2)步中提取的稳态压强场和脉动压强场。

进一步的，第5)步具体包括：

5.1)选择固体网格一个需要施加压强的单元表面，计算其单位法向矢量与面积S_e，并统计该单元表面的节点数n；

5.2)遍历所有单元，按照下式计算各单元表面各个节点的力载荷：

5.3)按下式计算单元节点载荷F_e在总载荷中的索引值，其中j为节点编号；

5.4)将单元所有节点的载荷按照3)中的索引加到总载荷向量中，并将载荷向量非0位置的索引值及对应元素写入文件Load_i.dat文件中；

5.5)遍历一个稳定周期的时间步，重复5.1)至5.4)步，将所有载荷向量输出。

进一步的，第6)步具体包括：

6.1)设置模态叠加法计算时间步，总时间设置为透平机械的10-20个旋转周期，时间步增量与第2)步中气动分析的时间步增量一致；

6.2)读取第4)步中获得的m阶固有频率f以及振型[φ]；

6.3)读取第5)步获得的所有载荷文件，将非0元素按照索引值填入载荷向量；

6.4)利用GPU计算一个周期内所有时间步的解耦坐标系下的激振力向量同时利用CPU计算阻尼固有频率ω_d

6.5)在CPU上利用梯形法计算所有周期的时间步上的杜哈美积分，获得解耦位移

其中，杜哈美积分中的sinω_d1(t-τ)项只需使用第一个周期的值，此后的周期均是第一个周期的重复；

6.6)利用解耦位移在GPU上计算所有周期的时间步上的实际位移

进一步的，第7)步具体包括：

7.1)将整圈叶片扇区剖分开，扇区总数为R个；CPU计算其中前个扇区，GPU计算后个扇区；其中，K为GPU和CPU计算能力的比值；

7.2)选择被分配扇区中单元，根据单元的节点坐标计算单元高斯节点处的应变矩阵[B]和本构矩阵[D]，其中N_i为单元第i个节点在高斯节点处的插值函数，表示对局部坐标的偏微分，ne为该单元的总节点数；

[B]＝[[B]₁ [B]₂ [B]₃ … [B]_i [B]_i+1 … [B]_ne]

7.3)设需要扩展的时间步为第j到第j+s步，提取该单元节点在时间步第j步到第j+s步的位移向量u^e_i,i+s；利用下式计算单元高斯节点在第j时间步到第j+s时间步的应力分量；

S^e_j,j+s＝[D][B]u^e_j,j+s

u^e_j,j+s＝[u^e_j u^e_j+1 u^e_j+2 … u^e_j+s-1 u^e_j+s]

S^e_j,j+s＝[S^e_j S^e_j+1 S^e_j+2 … S^e_j+s-1 S^e_j+s]

7.4)求第i步到第i+s步应力分量的主应力，即求解矩阵M^e的特征值S₁₁，S₂₂，S₃₃：

S^e＝[σ_x σ_y σ_z τ_xy τ_xz τ_yz]^T。

计算其VonMises等效应力

7.5)遍历单元内所有高斯节点，并计算各个高斯节点的真实坐标；

7.6)遍历所有单元，重复7.2)至7.5)步，对于第i步第j个扇区的应力结果，将获得的该扇区上单元编号及其对应的高斯节点的真实坐标和应力(包括应力分量，主应力以及等效应力)，输出到文件中Stress_i_j.dat。

与传统的计算方法相比，采用本发明的分析方法具有以下有益效果：

本发明在步骤2)中采用了采用非定常计算方法获得叶片表面的压强分布，可以有效提高叶片表面压强场的计算精度；步骤3)中通过将叶片流体区域网格向固体区域网格压强进行插值获得准确的气流激振力载荷，提高了气流激振力的计算精度；步骤7)在位移结果扩展为应力结果的步骤中将叶片划分为不同扇区，在各个扇区采用CPU+GPU异构并行计算的方式，可以大幅度提升计算速度，减少计算时间。

综上所述，本发明旨在提供一种透平机械叶片的振动应力数值分析方法，一方面可以缩短叶片设计周期并节约试验成本，另一方面可以减小以往采用理论分析和理论-实验相结合的半经验分析等方法对模型过度简化而引起的偏差，从而获得可靠的分析结果，对透平叶片的设计和制造具有重要的意义。

附图说明

图1为本发明一种透平机械叶片的振动应力数值分析方法的流程图。

图2为透平叶片的网格剖分示意图；其中图2(a)为叶片固体区域的FEM模型，图2(b)为叶片流体区域的CFD模型。

图3为一个气动周期的示意图。

图4为稳态压强和脉动压强的计算示意图。

图5为叶片固体区域表面网格中心示意图。

图6为叶片固体区域—流体区域网格插值示意图；其中1，2，3，4分别为流体区域网格节点的编号。

图7为叶片不同扇区编号示意图；其中图7(a)为整体示意图；图7(b)为局部示意图。

图8为实例中叶片在第10个时间步时1％叶高位置的压强分布。

图9为实例中叶片在第10个时间步时50％叶高位置的压强分布。

图10为实例中叶片在第10个时间步时99％叶高位置的压强分布。

图11为实例中叶片在第10个时间步时固体区域表面的激振力分布；其中图11(a)中为轴向激振力分布，图11(b)中为切向激振力分布。

图12为实例中整圈叶片模态振型计算结果；其中图12(a)为1阶节圆振动，图12(b)为1阶2节径振动，图12(c)为1阶4节径振动图12(d)为2阶节圆振动，图12(e)为2阶2节径振动，图12(f)为2阶4节径振动，图12(g)为3阶节圆振动，图12(h)为3阶2节径振动，图12(i)为3阶4节径振动。

图13为实例中整圈叶片的振动位移响应和等效应力计算结果；其中图13(a)为第10个时间步时的位移响应分布云图，图13(b)为第50个时间步时的位移响应分布云图，图13(c)为第10个时间步时的振动应力分布云图，图13(d)为第50个时间步时的振动应力分布云图。

图14为实例中叶片在拉筋圆角附近某点的振动应力变化曲线。

具体实施方式

下面结合附图对本发明做进一步详细说明。

请参阅图1所示，本发明一种透平机械叶片的振动应力数值分析方法，包括以下步骤：

1)、给定透平机械叶片的三维模型和材料参数与运行参数，在网格剖分软件中分别建立整圈叶片固体区域的FEM模型与周期对称的流体区域的CFD模型；并将有限元模型的节点编号以及节点坐标，单元编号以及单元对应的节点编号输出到文件Solid.fem中；将材料参数和运行参数输出到文件Blade.dat中。其中，材料参数包括叶片和轮盘的密度，弹性模量，泊松比；运行参数包括叶片的转速，工质类型，进口静压、温度、湍流度和叶片出口静压。

2)、对叶片周期对称模型进行气动分析；利用第1)步中获得的周期对称的流体区域的CFD模型，分析采用RNG k-ε湍流模型以及非平衡近壁模型，对模型设置工质类型，进口静压、温度、湍流度和叶片出口静压；首先进行定常计算，静叶栅和动叶栅采用冻结转子法，然后将定常结果作为初始场进行非定常计算，瞬态转子的方式进行耦合，取动叶经过一个静叶栅间距时间作为一个气动周期ΔT，时间步长取一个周期ΔT的计算获得非定常计算收敛后的一个周期各个时间步的叶片表面的压强场分布，并计算稳态压强场和脉动压强场，，将获得的稳态压强场和脉动压强场的压强场数据存储到文件Pressure.dat中，数据包括点坐标以及对应的压强大小。

利用一个气动周期各个时间步的压强场结果计算稳态压强场和脉动压强场，计算方式如下：

2.1)提取一个节点在不同时间步的压强值P_(t)；

2.2)取一个周期内P(t)的时均值作为稳态压强

2.3)计算该点各个时间步的脉动压强

2.4)遍历所有节点，获得稳态压强场和脉动压强场P′_(t)，并储存在文件Pressure.dat；

3)、压强分布数据转化；读取第1)步获得的solid.fem文件中的单只叶片固体网格数据，计算叶片表面单元中心的坐标，如图2；读取第2步获得文件Pressure.dat中的压力场数据，并读取压力场中节点在CFD模型中的网格数据，并对固体表面中心的位置进行插值，获得这些位置在一个气动周期不同时间步的压强。插值算法的流程如下：

3.1)计算固体网格中一个叶片单元表面中心点的坐标，在CPU和GPU上多线程并行判断该点是否落在流体区域网格的四边形内，记下流体区域网格四边形的四个节点；

3.2)对这四个节点组成的四边形进行剖分：假设四个节点坐标分别为(x1,y1)，(x2,y2)，(x3,y3)，(x4,y4)，在四条边分别均匀地插入9个节点，对应节点连接，将该四边形剖分为10*10的小四边形区域。其中如图标注的称为小四边形的经度，称为小四边形的纬度，对经度i和纬度j的小四边形的左下角节点坐标为

依次类推，可以计算所有剖分节点坐标。

3.3)根据剖分获得的所有小四边形的节点坐标，继续在CPU和GPU上同时判断该点是否落在这10*10的小四边形内，以落在其中的小四边形为新的四边形然后返回第3.2)步继续剖分，如此重复4次，记下这4次剖分过程中，落在其中的四边形的经度(i₁,i₂,i₃,i₄)和纬度(j₁,j₂,j₃,j₄)；

3.4)计算插值坐标

取ξ₁＝-1,ξ₂＝1,ξ₃＝1,ξ₄＝-1，η₁＝-1,η₂＝-1,η₃＝1,η₄＝1。

按照下式进行插值，获得该中心点的压强：

3.5)遍历所有固体网格中单元表面中心点，重复3.1)到3.4)步，获得所有单元中心点的压强；

3.6)按照3.1)—3.5)计算一个周期所有时间步压强场在所有固体网格单元表面中心点的插值，包括第2)步中提取的稳态压强场和脉动压强场；

4)、整圈叶片有限元模态分析；根据第2)步中获得稳态的压强插值施加在固体网格单元表面上，对整圈叶片施加该第1)步获得转速的离心力载荷，在有限元软件中求解其预应力场，并将整圈叶片有限元网格的单元编号与高斯节点的预应力结果输出到文件GSS.dat，节点编号与其对应的预应力结果输出到文件NSS.dat；然后进行考虑预应力场和旋转软化的模态分析，获得整圈叶片阶数m(通常取200-300)的固有频率f以及振型[φ]，并输出到文件Mode.dat。

5)、整圈叶片激振力施加；根据第3)步中获得的一个气动周期内不同时间步的固体网格的脉动压强，计算透平叶片表面单元每个节点的受力，获得所有时间步叶片表面的节点力载荷向量，并将一个气动周期内所有时间步的载荷向量压缩输出到文件中。流程如下：

5.1)选择固体网格一个需要施加压强的单元表面，计算其单位法向矢量与面积S_e，并统计该单元表面的节点数n；

5.2)遍历所有单元，按照下式计算各单元表面各个节点的力载荷：

5.3)按下式计算单元节点载荷F_e在总载荷中的索引值，其中j为节点编号；

5.4)将单元所有节点的载荷按照3)中的索引加到总载荷向量中，并将载荷向量非0位置的索引值及对应元素写入文件Load_i.dat文件中；

5.5)遍历一个稳定周期的时间步，重复1)至4)步，将所有载荷向量输出；

6)、模态叠加法求解振动位移响应；首先设置计算时间步，读取第5)步获得的载荷向量文件，并结合第2)步获得的固有频率以及振型，利用杜哈美积分计算在脉动气流激振力整圈叶片位移响应。具体流程如下：

6.1)设置模态叠加法计算时间步，总时间设置为透平机械的10-20个旋转周期，时间步增量与第2)步中气动分析的时间步增量一致；

6.2)读取第4)步中获得的m阶固有频率f以及振型[φ]；

6.3)读取第5步获得的所有载荷文件，将非0元素按照索引值填入载荷向量；

6.4)利用GPU计算一个周期内所有时间步的解耦坐标系下的激振力向量同时利用CPU计算阻尼固有频率ω_d

6.5)在CPU上利用梯形法计算所有周期的时间步上的杜哈美积分，获得解耦位移

其中，杜哈美积分中的sinω_d1(t-τ)项只需使用第一个周期的值即可，此后的周期均是第一个周期的重复。

6.6)利用解耦位移在GPU上计算所有周期的时间步上的实际位移

7)、位移响应结果扩展为应力结果；在CPU和GPU上同时将叶片不同扇区的位移响应转化为应力向量，求解各个单元高斯节点的应力分量，并外插获得节点的应力分量，最后将高斯节点和节点的应力分量转化为主应力以及VonMises等效应力；具体流程如下

7.1)按照图7示的方式将整圈叶片扇区剖分开，设扇区总数为R个，CPU和GPU分别分配其中一定比例(1：K)的扇区，K值根据CPU和GPU计算能力进行测试确定，测试标准是分配的单元矩阵比例需要保证CPU和GPU计算时间基本相同；则CPU计算其中前个扇区，GPU计算后个扇区；

7.2)选择被分配扇区中单元，根据单元的节点坐标计算单元高斯节点处的应变矩阵[B]和本构矩阵[D]，其中N_i为单元第i个节点在高斯节点处的插值函数，表示对局部坐标的偏微分，ne为该单元的总节点数；

[B]＝[[B]₁ [B]₂ [B]₃ … [B]_i [B]_i+1 … [B]_ne]

7.3)设需要扩展的时间步为第j到第j+s步，提取该单元节点在时间步第j步到第j+s步的位移向量u^e_i,i+s；利用下式计算单元高斯节点在第j时间步到第j+s时间步的应力分量；

S^e_j,j+s＝[D][B]u^e_j,j+s

u^e_j,j+s＝[u^e_j u^e_j+1 u^e_j+2 … u^e_j+s-1 u^e_j+s]

S^e_j,j+s＝[S^e_j S^e_j+1 S^e_j+2 … S^e_j+s-1 S^e_j+s]

7.4)求第i步到第i+s步应力分量的主应力，即求解矩阵M^e的特征值S₁₁，S₂₂，S₃₃：

S^e＝[σ_x σ_y σ_z τ_xy τ_xz τ_yz]^T

计算其VonMises等效应力

7.5)遍历单元内所有高斯节点，并计算各个高斯节点的真实坐标；

7.6)遍历所有单元，重复7.2)至7.5)步，对于第i步第j个扇区的应力结果，将获得的该扇区上单元编号及其对应的高斯节点的真实坐标和应力(包括应力分量，主应力以及等效应力)，输出到文件中Stress_i_j.dat。

8)、振动应力结果提取与考核；提取所有扇区叶片在围带圆角，叶根平台圆角，拉筋圆角考核位置在第2)步中预应力分析结果中的最大VonMises等效应力，并提取这些位置在不同时间步的单元或节点的VonMises等效应力，将二者按照线性形式叠加，绘制振动应力时域曲线。

下面结合具体的实例对本发明涉及的透平机械叶片的振动应力数值分析方法进行说明。

本实例以某带有连接件(围带和拉筋)的整圈叶片为例，三维模型采用四齿枞树形叶根，通过拉筋以及围带发生接触作用将整圈叶片连接起来。叶片在稳定工况下，转速为1500rpm，按照气动数据计算其压强场分布如图8，图9，图10，其中图8为实例中叶片在第10个时间步时1％叶高位置的压强分布，图9为实例中叶片在第10个时间步时50％叶高位置的压强分布，图10为实例中叶片在第10个时间步时99％叶高位置的压强分布。按照计算其固体网格表面节点力的分布如图11，其中图11(a)为实例中叶片在第10个时间步时固体区域表面的轴向激振力分布，图11(b)为实例中叶片在第10个时间步时固体区域表面的切向激振力分布；将稳态气流力施加在固体网格模型上，并计算整圈叶片在离心力和气流力下的固有振型如图12，将瞬态气流力施加在固体网格模型上，并计算其在激振力下的位移响应和等效应力分布云图如图13，其中图13(a)为第10个时间步时的位移响应分布云图，图13(b)为第50个时间步时的位移响应分布云图，图13(c)为第10个时间步时的振动应力分布云图，图13(d)为第50个时间步时的振动应力分布云图；，最后提取其在拉筋圆角附近的振动应力变化曲线如图14。

从图中可以看出叶片的振动响应最大值在叶顶位置分布，而振动应力集中位置主要在叶根平台圆角，拉筋圆角，以及围带圆角附近，计算数值与实验测量结果接近；最大应力在拉筋圆角附近，从其振动应力变化曲线可以看出叶片的振动应力主要受到前四阶气流激振力的影响，与实验测量结果相同；说明了本发明提出的方法可以比较精确地计算透平机械叶片在气流激振力下的振动应力分布。