CN106528982A - 一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法 - Google Patents

Abstract

本发明公开了一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法，结合实际汽轮机叶片失谐状况，首先进行三维接触有限元分析，获得整圈叶片各拉筋和围带接触面各节点处的支反力；其次采用弹簧阻尼单元建立接触面关系，在各弹簧阻尼单元分别考虑正压力、面积、摩擦系数，并基于Mindlin微动滑移摩擦模型和谐波平衡法，建立干摩擦阻尼失谐振动分析模型；然后考虑预应力并采用谐响应分析通过迭代求解获得各叶片振动响应的收敛值；通过循环计算，获得各干摩擦失谐叶片的振动幅值响应曲线，进一步提取可得到固有频率及振型。本发明对叶片干摩擦阻尼失谐分析及改善失谐叶盘振动具有重要意义。

Description

一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法

【技术领域】

本发明属于汽轮机叶片振动分析领域，具体涉及一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法。

【背景技术】

在透平机械叶盘***中，常常由于制造误差、材料磨损等因素而产生失谐现象，***振动能量集中在了少数叶片区域，导致叶片振动响应远大于其他部位，发生振动局部化现象。振动局部化现象使局部叶片应力增大，增加了叶片高周疲劳破坏的危险性。

在汽轮机叶片的设计中常采用干摩擦阻尼结构来提高叶片刚度和阻尼，以降低叶片的振动水平。汽轮机末级长叶片的设计采用了凸台拉筋和整体围带相结合的阻尼结构，由于接触干摩擦阻尼力为非线性，相比于单一形式的摩擦阻尼结构，无论是在实际接触状态还是摩擦减振过程方面都更加复杂。在实际叶盘***中，干摩擦阻尼结构由于安装、磨损等因素，每个叶片拉筋和围带接触面的正压力载荷以及接触状态的不同，会导致整圈叶片摩擦阻尼特性产生失谐，其实际摩擦减振作用与协调叶盘分析结果产生了差异，对叶盘***的振动特性研究和干摩擦阻尼结构的设计带来新的困难，这样就需要一套成熟稳定的方法来准确获得这种复杂整圈干摩擦阻尼失谐叶片的振动响应特性。

【发明内容】

本发明的目的在于克服上述现有技术的缺点，提供一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法。

为达到上述目的，本发明采用以下技术方案予以实现：

一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法，包括以下步骤：

1)在考虑拉筋和围带间隙接触、叶轮轮缘接触以及叶片几何非线性大变形基础上，构造与实际装配情况相一致的整圈叶片-叶轮三维有限元数值模型，并采用能够准确描述复杂几何体特征的三维实体单元对叶片进行网格剖分，得到有限元离散模型，根据所选择的单元形成对应的单元刚度矩阵，集成总刚度矩阵；

在叶根齿与轮缘齿接触面、相邻叶片拉筋接触面和相邻叶片围带接触面间建立面-面接触关系，对叶轮进汽侧轴向端面设置切向位移约束，对叶轮出汽侧轴向端面设置切向和轴向位移约束，对整个模型施加工作转速工况下所对应的角速度载荷，在考虑位移约束、力载荷以及接触约束的情况下加载边界条件，形成有限元方程组；

2)采用有限元数值分析方法在工作转速工况下对叶片进行三维接触有限元分析，算法选用增强拉格朗日法，通过构造修正的势能泛函使有约束问题转化为无约束问题，并选择稀疏矩阵直接消元法对有限元方程进行求解，得到整圈叶片三维接触有限元分析结果，分别获得该转速工况下整圈叶片在拉筋和围带各接触面上各节点的法向支反力，叠加后获得对应工况下接触面正压力载荷；

3)根据实际测得或掌握的数据对数值计算结果各接触面上节点的法向支反力进行对比修正，得到真实情况下整圈干摩擦阻尼失谐叶片在拉筋和围带各接触面上各节点处正压力的分布，然后对叶根部分的所有节点施加完全位移约束，再分别在拉筋和围带的干摩擦接触面各节点间采用弹簧阻尼单元建立接触关系，每个弹簧阻尼单元包括一个刚度矩阵单元和一个阻尼矩阵单元，摩擦接触面通过多个弹簧阻尼单元并联起来，每个弹簧阻尼单元具有各自的刚度值和阻尼值，获得干摩擦阻尼叶片振动特性有限元分析模型；

4)选取该转速下整圈失谐叶片在拉筋和围带接触面各节点处正压力作为叶片振动响应分析的输入参数，同时考虑接触面各节点区域的面积、摩擦系数，通过Mindlin摩擦阻尼分析模型分别计算得到各弹簧阻尼单元上的初始等效刚度系数K_eq和等效阻尼系数C_eq；

根据实际情况测得的气流力数据，对叶片施加对应幅值和频率的简谐气流激振力载荷，采用基于Mindlin微动滑移摩擦模型和谐波平衡法的干摩擦阻尼动态分析方法来描述接触面间的摩擦阻尼特性，在考虑预应力条件下采用基于模态叠加法的谐响应分析方法进行迭代计算获得叶片振动响应的收敛值，更新频率后获得叶片在一系列激振频率下的振动幅值响应曲线，提取共振时各弹簧阻尼单元上的刚度值和阻尼值，得到叶片在对应工况下的固有频率及振型。

本发明进一步的改进在于：

所述步骤4)中获得叶片振动响应收敛值的迭代分析过程为：

4-1)对整圈叶片施加对应幅值和频率ω的简谐气流激振力载荷，假设整圈某两支相邻叶片在工作转速工况下拉筋和围带接触面间的相对运动位移幅值分别为A₁和A₂；

4-2)根据Mindlin微动滑移模型理论，对于两种材料相同的球-平面接触对，接触对在切向力作用下的切向接触刚度为：

式中a为接触半径；ν物体的泊松比；E为物体的弹性模量；

4-3)采用谐波平衡法分析***的响应，对摩擦力f进行傅立叶展开，在权衡计算量和计算精度后取一阶谐波：

f(A,θ)＝f_C(A)cosθ+f_S(A)sinθ (14)

式中f_c(A)为摩擦力的一阶余弦分量，f_s(A)为摩擦力的一阶正弦分量，θ为相位；

4-4)基于Mindlin摩擦模型，推导得到无量纲摩擦力随无量纲相对运动位移变化的迟滞曲线方程表达式；

4-5)根据谐波平衡法的分析，将无量纲位移代入获得无量纲摩擦力接触面间的摩擦力由弹性力和阻尼力叠加来计算，最终得到摩擦接触面间的无量纲等效刚度系数和等效阻尼系数：

4-6)根据干摩擦失谐叶片在拉筋和围带接触面各节点处正压力、摩擦系数、面积大小等因素计算对应弹簧阻尼单元上的等效刚度系数和等效阻尼系数，然后将其分别加载到对应的刚度矩阵单元和阻尼矩阵单元中，在对应幅值和频率的简谐气流激振力载荷作用下考虑预应力，并基于模态叠加法对整圈失谐叶片进行谐响应分析，获得整圈叶片接触面间所有的相对运动位移幅值；

谐响应分析可获得接触面间弹簧阻尼单元对应两节点的位移幅值B₁、B₂以及两幅值间的相位差通过计算得到单个周期内摩擦面间的最大相对滑移量A，这样针对两相邻叶片可以分别得到拉筋和围带接触面的相对运动位移幅值A₁'、A'₂，进一步可以分析得到整圈叶片所有接触面间的相对运动位移幅值；

4-7)对比计算获得的整圈所有接触面间的相对运动位移幅值和前一次计算的整圈所有接触面间的相对运动位移幅值，如果所有差别均在误差允许范围内，则计算收敛，否则更新所有相对运动位移幅值并返回到步骤4-2)继续迭代计算，直到计算收敛；

4-8)在迭代收敛后，获得激振频率为ω时的整圈各叶片振动响应的收敛值，保存当前计算结果，接着返回到步骤4-1)，在正压力等参数不变的情况下改变激振频率计算下一个循环，直到计算完相关的一系列激振频率，获得叶片在一系列激振频率下振动响应，再进一步提取共振时弹簧各阻尼单元的刚度值和阻尼值，计算得到该工况下所对应的叶片固有频率及振型。

与现有技术相比，本发明具有以下有益效果：

本发明结合了实际汽轮机叶片的整圈干摩擦阻尼失谐状况，考虑了干摩擦接触状态的复杂性，采用了Mindlin微动滑移模型描述叶片接触面间的摩擦阻尼特性，在考虑预应力条件下采用了基于模态叠加法的谐响应分析进行求解，基于接触面间运动位移幅值和相位差求得接触面滑移量，进而通过迭代求解得到给定激振频率下叶片振动响应的收敛值，计算中还分别考虑了各节点弹簧阻尼单元的刚度值和阻尼值，分析结果精度高，可靠性强，降低了复杂非线性问题的分析时间，为具有拉筋和围带的整圈汽轮机干摩擦阻尼失谐叶片提供了一种准确有效的分析方法。

【附图说明】

图1是本发明的振动分析流程图；

图2是具有凸台拉筋和整体围带的整圈汽轮机叶片的模型示意图；

图3是各工况下拉筋和围带接触面正压力载荷随转速变化曲线图；

图4是采用弹簧阻尼单元建立的接触面间的局部有限元网格模型示意图；

图5是Mindlin摩擦模型中摩擦力与相对运动位移之间的关系图；

图6是某工况下整圈干摩擦阻尼失谐叶片的振动幅频响应曲线图。

【具体实施方式】

下面结合附图对本发明做进一步详细描述：

参见图1，本发明有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法，包括以下步骤：

1)在考虑拉筋和围带间隙接触、叶轮轮缘接触以及叶片几何非线性大变形基础上，构造与实际装配情况相一致的整圈叶片-叶轮三维有限元数值模型(见图2)，并采用能够准确描述复杂几何体特征的三维实体单元对叶片进行网格剖分，得到有限元离散模型，根据所选择的单元形成对应的单元刚度矩阵，集成总刚度矩阵；

在叶根齿与轮缘齿接触面、相邻叶片拉筋接触面和相邻叶片围带接触面间建立面-面接触关系，对叶轮进汽侧轴向端面设置切向位移约束，对叶轮出汽侧轴向端面设置切向和轴向位移约束，对整个模型施加工作转速工况下所对应的角速度载荷，在考虑位移约束、力载荷以及接触约束的情况下加载边界条件，形成有限元方程组；

2)采用有限元数值分析方法在工作转速工况下对叶片进行三维接触有限元分析，算法选用增强拉格朗日法，通过构造修正的势能泛函使有约束问题转化为无约束问题，接触问题的控制方程为：

[K^*]{u^*}＝{F^*} (1)

其中[K^*]、{u^*}和{F^*}的表达式分别为：

[K_P]＝([N][C])^T[α][N][C] (3)

考虑到拉氏乘子λ的物理意义，用接触力代替λ，使其在迭代计算中作为已知量出现，并选择稀疏矩阵直接消元法对控制方程(1)进行求解，得到整圈叶片三维接触有限元分析结果后，分别获得该转速工况下整圈叶片在拉筋和围带各接触面上各节点的法向支反力，叠加后获得对应工况下接触面正压力载荷(见图3)；

3)根据实际测得或掌握的数据对数值计算的正压力(即各接触面上节点的法向支反力)进行对比修正，得到真实情况下整圈干摩擦阻尼失谐叶片在拉筋和围带各接触面上各节点处正压力的分布，然后对叶根部分的所有节点施加完全位移约束，再分别在拉筋和围带的干摩擦接触面各节点间采用弹簧阻尼单元建立接触关系，每个弹簧阻尼单元包括一个刚度矩阵单元和一个阻尼矩阵单元，摩擦接触面通过多个弹簧阻尼单元并联起来，每个弹簧阻尼单元具有各自的刚度值和阻尼值，这样获得干摩擦阻尼叶片振动特性有限元分析模型(见图4)；

对两点间的刚度矩阵和阻尼矩阵，考虑两自由度***，得到运动微分方程：

只考虑m₁和m₂之间的刚度和阻尼时，为

忽略质量m₁和m₂，得到局部刚度矩阵和阻尼矩阵

将上述两节点之间的刚度矩阵和阻尼矩阵推广至2个节点的12个自由度，得到刚度矩阵单元：

阻尼矩阵单元：

4)选取该转速下整圈失谐叶片在拉筋和围带接触面各节点处正压力作为叶片振动响应分析的输入参数，同时考虑接触面各节点区域的面积、摩擦系数，通过Mindlin摩擦阻尼分析模型分别计算得到各弹簧阻尼单元上的初始等效刚度系数K_eq和等效阻尼系数C_eq；

根据实际情况测得的气流力数据，对叶片施加对应幅值和频率的简谐气流激振力载荷，采用基于Mindlin微动滑移摩擦模型和谐波平衡法的干摩擦阻尼动态分析方法来描述接触面间的摩擦阻尼特性，在考虑预应力条件下采用基于模态叠加法的谐响应分析方法进行迭代计算获得叶片振动响应的收敛值，更新频率后获得叶片在一系列激振频率下的振动幅值响应曲线和振动应力曲线，提取共振时各弹簧阻尼单元上的刚度值和阻尼值，得到叶片在对应工况下的固有频率及振型。

其中，获得叶片振动响应收敛值的迭代分析过程为：

4-1)对整圈叶片施加对应幅值和频率ω的简谐气流激振力载荷，假设整圈某两支相邻叶片在工作转速工况下拉筋和围带接触面间的相对运动位移幅值分别为A₁和A₂；

4-2)根据Hertz弹性接触理论，对于两种材料相同的球-平面接触对，得到接触半径：

式中a为接触半径，N为接触区域的接触压力，R为球面的半径，ν物体的泊松比，E为物体的弹性模量，k为系数，k＝(1-v²)/(πE)；

根据Mindlin微动滑移模型理论，球-平面接触对在切向力作用下的切向接触刚度：

令临界位移幅值为A₀，A₀＝μN/K_d，对叶片摩擦接触面间的相对运动位移幅值A₁和A₂进行无量纲化，即

4-3)采用谐波平衡法分析***的响应，假设：

u＝Acosθ，

式中u为相对运动位移，为无量纲相对运动位移，θ为相对运动位移的相位，对摩擦力f进行傅立叶展开，在权衡计算量和计算精度后取一阶谐波得：

f(A,θ)＝f_C(A)cosθ+f_S(A)sinθ (14)

式中f_c(A)为摩擦力的一阶余弦分量，f_s(A)为摩擦力的一阶正弦分量，θ为相位，对f_c(A)和f_s(A)的表达式进行无量纲化后得：

4-4)基于Mindlin摩擦模型，得到无量纲摩擦力随无量纲相对运动位移变化的迟滞曲线方程表达式(见图5)，当相对运动位移时，曲线方程表达式为：

当相对运动位移时，曲线方程表达式为：

其中无量纲相对运动位移无量纲摩擦力无量纲微动滑移阶段摩擦力幅值无量纲相对运动位移幅值无量纲相对运动位移

4-5)根据谐波平衡法的分析，将无量纲位移代入获得无量纲摩擦力接触面间的摩擦力由弹性力和阻尼力叠加来计算，最终得到摩擦接触面间的无量纲等效刚度系数和等效阻尼系数：

4-6)根据干摩擦失谐叶片在拉筋和围带接触面各节点处正压力、摩擦系数、面积大小等因素计算对应弹簧阻尼单元上的等效刚度系数和等效阻尼系数，然后将其分别加载到对应的刚度矩阵单元和阻尼矩阵单元中，在对应幅值和频率的简谐气流激振力载荷作用下考虑预应力，并基于模态叠加法对整圈失谐叶片进行谐响应分析，获得整圈叶片接触面间所有的相对运动位移幅值；

谐响应分析可获得接触面间弹簧阻尼单元对应两节点的位移幅值B₁、B₂以及两幅值间的相位差通过计算得到单个周期内摩擦面间的最大相对滑移量A，这样针对两相邻叶片可以分别得到拉筋和围带接触面的相对运动位移幅值A₁'、A'₂，进一步可以分析得到整圈失谐叶片所有接触面间的相对运动位移幅值；

4-7)对比计算获得的所有整圈接触面间的相对运动位移幅值和前一次计算的所有整圈接触面间的相对运动位移幅值，如果所有差别均在误差允许范围内，则计算收敛，否则更新所有相对运动位移幅值并返回到步骤4-2)继续迭代计算，直到计算收敛；

4-8)在迭代收敛后，获得激振频率为ω时的整圈叶片振动响应的收敛值，保存当前计算结果，接着返回到步骤4-1)，在正压力等参数不变的情况下改变激振频率计算下一个循环，直到计算完相关的一系列激振频率，获得叶片在一系列激振频率下振动响应(见图6)，再经进一步提取共振时各弹簧阻尼单元的刚度值和阻尼值，计算得到该工况下所对应的叶片固有频率及振型。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (2)

1.一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法，其特征在于，包括以下步骤：

1)在考虑拉筋和围带间隙接触、叶轮轮缘接触以及叶片几何非线性大变形基础上，构造与实际装配情况相一致的整圈叶片-叶轮三维有限元数值模型，并采用能够准确描述复杂几何体特征的三维实体单元对叶片进行网格剖分，得到有限元离散模型，根据所选择的单元形成对应的单元刚度矩阵，集成总刚度矩阵；

在叶根齿与轮缘齿接触面、相邻叶片拉筋接触面和相邻叶片围带接触面间建立面-面接触关系，对叶轮进汽侧轴向端面设置切向位移约束，对叶轮出汽侧轴向端面设置切向和轴向位移约束，对整个模型施加工作转速工况下所对应的角速度载荷，在考虑位移约束、力载荷以及接触约束的情况下加载边界条件，形成有限元方程组；

2)采用有限元数值分析方法在工作转速工况下对叶片进行三维接触有限元分析，算法选用增强拉格朗日法，通过构造修正的势能泛函使有约束问题转化为无约束问题，并选择稀疏矩阵直接消元法对有限元方程进行求解，得到整圈叶片三维接触有限元分析结果，分别获得该转速工况下整圈叶片在拉筋和围带各接触面上各节点的法向支反力，叠加后获得对应工况下接触面正压力载荷；

3)根据实际测得或掌握的数据对数值计算结果各接触面上节点的法向支反力进行对比修正，得到真实情况下整圈干摩擦阻尼失谐叶片在拉筋和围带各接触面上各节点处正压力的分布，然后对叶根部分的所有节点施加完全位移约束，再分别在拉筋和围带的干摩擦接触面各节点间采用弹簧阻尼单元建立接触关系，每个弹簧阻尼单元包括一个刚度矩阵单元和一个阻尼矩阵单元，摩擦接触面通过多个弹簧阻尼单元并联起来，每个弹簧阻尼单元具有各自的刚度值和阻尼值，获得干摩擦阻尼叶片振动特性有限元分析模型；

4)选取该转速下整圈失谐叶片在拉筋和围带接触面各节点处正压力作为叶片振动响应分析的输入参数，同时考虑接触面各节点区域的面积、摩擦系数，通过Mindlin摩擦阻尼分析模型分别计算得到各弹簧阻尼单元上的初始等效刚度系数K_eq和等效阻尼系数C_eq；

根据实际情况测得的气流力数据，对叶片施加对应幅值和频率的简谐气流激振力载荷，采用基于Mindlin微动滑移摩擦模型和谐波平衡法的干摩擦阻尼动态分析方法来描述接触面间的摩擦阻尼特性，在考虑预应力条件下采用基于模态叠加法的谐响应分析方法进行迭代计算获得叶片振动响应的收敛值，更新频率后获得叶片在一系列激振频率下的振动幅值响应曲线，提取共振时各弹簧阻尼单元上的刚度值和阻尼值，得到叶片在对应工况下的固有频率及振型。

2.根据权利要求1所述的一种有拉筋和围带的干摩擦阻尼失谐叶片的振动分析方法，其特征在于，所述步骤4)中获得叶片振动响应收敛值的迭代分析过程为：

4-1)对整圈叶片施加对应幅值和频率ω的简谐气流激振力载荷，假设整圈某两支相邻叶片在工作转速工况下拉筋和围带接触面间的相对运动位移幅值分别为A₁和A₂；

4-2)根据Mindlin微动滑移模型理论，对于两种材料相同的球-平面接触对，接触对在切向力作用下的切向接触刚度为：

$K_d\≈K_{\mathrm{\∞}}=\frac{F}{x_{\mathrm{\∞}}}=\frac{4Ea}{2(1+v)(2-v)}---\left(12\right)$

式中a为接触半径；ν物体的泊松比；E为物体的弹性模量；

4-3)采用谐波平衡法分析***的响应，对摩擦力f进行傅立叶展开，在权衡计算量和计算精度后取一阶谐波：

f(A,θ)＝f_C(A)cosθ+f_S(A)sinθ (14)

式中f_c(A)为摩擦力的一阶余弦分量，f_s(A)为摩擦力的一阶正弦分量，θ为相位；

4-4)基于Mindlin摩擦模型，推导得到无量纲摩擦力随无量纲相对运动位移变化的迟滞曲线方程表达式；

4-5)根据谐波平衡法的分析，将无量纲位移代入获得无量纲摩擦力接触面间的摩擦力由弹性力和阻尼力叠加来计算，最终得到摩擦接触面间的无量纲等效刚度系数和等效阻尼系数：

${\tilde{K}}_{eq}=\frac{K_{eq}}{K_d}=\frac{\frac{f_C}{\mathrm{\μ}N}\·\mathrm{\μ}N}{K_d\·\frac{A}{A_0}\·A_0}=\frac{{\tilde{f}}_C\left(\tilde{A}\right)}{\tilde{A}}=\frac{1}{\mathrm{\π}\tilde{A}}{\∫}_0^{2\mathrm{\π}}f(\tilde{A},\mathrm{\θ})\cos \mathrm{\θ}d\mathrm{\θ}---\left(20\right)$

${\tilde{C}}_{eq}=\frac{C_{eq}}{K_d/\left(\mathrm{\π}\mathrm{\ω}\right)}=-\frac{f_S}{\mathrm{\ω}A\·K_d/\left(\mathrm{\π}\mathrm{\ω}\right)}=-\mathrm{\π}\frac{{\tilde{f}}_S\left(\tilde{A}\right)}{\tilde{A}}=-\frac{1}{\tilde{A}}{\∫}_0^{2\mathrm{\π}}\tilde{f}(\tilde{A},\mathrm{\θ})sin\mathrm{\θ}d\mathrm{\θ}---\left(21\right)$

4-6)根据干摩擦失谐叶片在拉筋和围带接触面各节点处正压力、摩擦系数、面积大小等因素计算对应弹簧阻尼单元上的等效刚度系数和等效阻尼系数，然后将其分别加载到对应的刚度矩阵单元和阻尼矩阵单元中，在对应幅值和频率的简谐气流激振力载荷作用下考虑预应力，并基于模态叠加法对整圈失谐叶片进行谐响应分析，获得整圈叶片接触面间所有的相对运动位移幅值；

谐响应分析可获得接触面间弹簧阻尼单元对应两节点的位移幅值B₁、B₂以及两幅值间的相位差通过计算得到单个周期内摩擦面间的最大相对滑移量A，这样针对两相邻叶片可以分别得到拉筋和围带接触面的相对运动位移幅值A′₁、A'₂，进一步可以分析得到整圈叶片所有接触面间的相对运动位移幅值；

4-7)对比计算获得的整圈所有接触面间的相对运动位移幅值和前一次计算的整圈所有接触面间的相对运动位移幅值，如果所有差别均在误差允许范围内，则计算收敛，否则更新所有相对运动位移幅值并返回到步骤4-2)继续迭代计算，直到计算收敛；

4-8)在迭代收敛后，获得激振频率为ω时的整圈各叶片振动响应的收敛值，保存当前计算结果，接着返回到步骤4-1)，在正压力等参数不变的情况下改变激振频率计算下一个循环，直到计算完相关的一系列激振频率，获得叶片在一系列激振频率下振动响应，再进一步提取共振时弹簧各阻尼单元的刚度值和阻尼值，计算得到该工况下所对应的叶片固有频率及振型。