CN110610049A - 一种叶片和机匣系统在碰摩故障下的力学特性分析方法 - Google Patents

Abstract

本发明属于转定子系统机械模拟技术领域，尤其涉及一种叶片和机匣系统在碰摩故障下的力学特性分析方法。该方法包括如下步骤：A1、对叶片进行离散化，获得叶片的梁单元模型和叶片的弹簧单元模型；A2、对机匣进行离散化，获得机匣的梁单元模型；A3、基于所述叶片的梁单元模型、叶片的弹簧单元模型和机匣的梁单元模型，构建叶片和机匣系统的有限元模型；A4、基于接触算法，并结合所述构建叶片和机匣系统的有限元模型，得到碰摩故障下的叶片和机匣系统的有限元模型；A5、对所述碰摩故障下的叶片和机匣系统的有限元模型进行降维处理，获得叶片和机匣系统的力学特性。该分析方法考虑了叶片根部的弹性支撑问题，得到的力学特性更为准确。

Description

一种叶片和机匣系统在碰摩故障下的力学特性分析方法

技术领域

本发明属于转定子系统机械模拟技术领域，尤其涉及一种叶片和机匣系 统在碰摩故障下的力学特性分析方法。

背景技术

在航空发动机中，以转子叶片和定子机匣为代表的一类系统会存在一些 结构故障，如错位、质量偏心和机匣变形等，会使转子叶片尖端和机匣产生 碰摩。

目前，很多学者对此进行了研究，提出了一些有关碰摩力学特性的分析 方法。这些方法使设计人员可以模拟出叶片和机匣出现碰摩故障时的力学特 性，从而对设计高性能、稳定运行的航空发动机具有一定的指导意义。

但是，现有的方法中叶片大多采用直板叶片，用欧拉梁单元建模，没有 考虑叶片根部的弹性支撑问题；在计算碰摩力时，也是采用基于欧拉梁的接 触算法。

发明内容

(一)要解决的技术问题

针对现有存在的技术问题，本发明提供一种叶片和机匣系统在碰摩 故障下的力学特性分析方法，该分析方法考虑了叶片根部的弹性支撑问 题，得到的力学特性更为准确。

(二)技术方案

本发明提供一种叶片和机匣系统在碰摩故障下的力学特性分析方 法，包括如下步骤：

A1、对叶片进行离散化，获得叶片的梁单元模型和叶片的弹簧单元 模型；

A2、对机匣进行离散化，获得机匣的梁单元模型；

A3、基于所述叶片的梁单元模型、叶片的弹簧单元模型和机匣的梁 单元模型，构建叶片和机匣系统的有限元模型；

A4、基于接触算法，并结合所述构建叶片和机匣系统的有限元模型， 得到碰摩故障下的叶片和机匣系统的有限元模型；

A5、对所述碰摩故障下的叶片和机匣系统的有限元模型进行降维处 理，获得叶片和机匣系统的力学特性。

进一步地，所述叶片的弹簧单元模型包括三个线弹簧和 两个扭簧

进一步地，所述叶片的弹簧单元模型满足如下公式：

式中：M_b、C_b、K_b、分别为叶片的质量矩阵、阻尼矩阵、刚度矩 阵和外力矩阵；u_b、分别为叶片的位移、速度和加速度。

进一步地，所述机匣的梁单元模型满足如下公式：

式中：M_c、C_c、K_c、分别为机匣的质量矩阵、阻尼矩阵、刚度矩 阵和外力矩阵；u_c、分别为机匣的位移、速度和加速度。

进一步地，所述步骤A3中，叶片和机匣系统的有限元模型满足如下 公式：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩 阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和 加速度。

进一步地，所述步骤A4中，接触算法为改进的拉格朗日乘子法。

进一步地，所述改进的拉格朗日乘子法中，包括如下步骤：

S1、将参数θ_C作为输入数据，输入预先训练的第一模型，获取接触 区域；

所述第一模型为基于历史参数θ与对应的接触区域训练而获得的接 触区域判定模型；

S2、在所述接触区域内，将预先获得的参数μ、 g_N作为输入数据，输入预先训练的第二模 型，获取沿s向和t向的碰摩力类型和滑动摩擦接触面与s/t向的夹角；

其中，μ为摩擦系数，分别为叶尖点A在t向的位移、速度 和加速度；分别为叶尖点A在s向的位移、速度和加速度； 分别为A’在t向的位移、速度和加速度；g_N为法向接触间隙；

所述第二模型为基于历史参数μ、u^t、u^s、u_， ^t、g_N与 对应的沿s向和t向的碰摩力类型训练而获得的碰摩力类型判定模型；

S3、将获得的沿s向和t向的碰摩力类型、滑动摩擦接触面与s/t向 的夹角和预先获得的作为输入数据，输入预先训练 的第三模型，获取拉格朗日乘子；

其中，分别为A’相对于A在t向的速度、加速 度和A’相对于A在s向的速度、加速度；

所述第三模型为基于历史参数沿s向和t向的碰摩力类型、滑动摩擦 接触面与s/t向的夹角、训练而获得的切向摩擦力判 定模型。

进一步地，所述步骤A5中，基于对叶片和机匣的收敛性分析，对碰 摩故障下的叶片和机匣系统的有限元模型进行降维处理。

进一步地，所述碰摩故障下的叶片和机匣系统的有限元模型满足如 下公式：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩 阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和 加速度；B_c,n,B_c,t和B_c,s分别为法向、周向和轴向接触约束矩阵；λ_N为法向拉 格朗日乘子。

(三)有益效果

本发明提供的方法，采用自编弹簧单元和铁木辛柯梁单元，建立了 柔性支承叶片和机匣的有限元模型。在不考虑叶片扭转和机匣面外运动 的情况下，采用自编的点面接触算法，结合向前增量拉格朗日乘子法和 向前增量罚函数法，得到了系统的碰摩动力学特性。同时，为了提高计 算效率，在求解时中采用了模型降维技术。该分析方法得到的力学特性 更为准确。

附图说明

图1为本发明提供的叶片和机匣系统在碰摩故障下的力学特性分析 方法的流程图；

图2本发明中叶片和机匣系统的整体结构示意图；

图3为本发明中叶片和机匣系统中叶片的结构示意图；

图4为本发明中叶片和机匣系统中机匣的结构示意图；

图5为本发明中拉格朗日乘子法的接触区域判断示意图；

图6为本发明中拉格朗日乘子法的接触区判断流程图；

图7为本发明中拉格朗日乘子法的A和A’的空间分布示意图；

图8为本发明中叶片和机匣接触示意图；

图9为本发明中拉格朗日乘子法中切向摩擦力的预判流程图；

图10为本发明中拉格朗日乘子法中切向摩擦力的计算流程图；

图11(a)为本发明中五个弹性支撑下的叶片有限元模型；

图11(b)-11(d)分别为本发明中叶片模态截断数(MTNB)对X^r-、 Y^r-和Z^r-振动响应的影响示意图；

图12(a)为本发明中给出了转速n∈[0,15000]rpm时，6倍转频激 励下临界共振点的示意图；

图12(b)为本发明中n＝15000rpm下叶片的前五阶振型；

图13(a)为本发明中机匣的有限元模型；

图13(b)为本发明中机匣的前11阶振型；

图13(c)、13(d)分别为机匣模态截断数(MTNC)对机匣上节 点1处X向和Y向振动响应的影响示意图；

图14为本发明中偏心量e对系统位移影响示意图；

图15为本发明中偏心量e对系统碰摩力响应示意图；

图16(a)为本发明中叶片转速为1900转每分时叶片和机匣的法向 接触间隙随时间的变化；

图16(b)为本发明中叶片转速为1900转每分时摩擦区域的对比；

图17(a)为本发明中叶尖点的Z向位移频谱图；

图17(b)为本发明中机匣节点1处Y向位移频谱图；

图18(a)为本发明中A₆共振带处的碰摩区域图；

图18(b)为本发明中A₆共振带处法向、切向和轴向的碰摩时域波 形；

图19为本发明中摩擦系数μ对位移幅频响应的影响示意图；

图20为本发明中摩擦系数μ对碰摩力幅频响应的影响示意图；

图21为本发明中叶片和机匣的频谱图；

图22为本发明中基于罚函数法和拉格朗日乘子法得到的机匣各阶模 态应变能的分布示意图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实 施方式，对本发明作详细描述。

以下为本发明中所涉及的字母及其对应含义：

A 叶尖点

A₀ 叶片截面面积

A′ 投影到XOY平面BC和OA线的交点

A″ 消除浸入后的A的最终位置

b₀,b_c 叶根宽度；机匣宽度

B,C 与A接触的机匣单元BC的节点

B_c,n,B_c,t,B_c,s 分别在法向、切向和轴向接触约束矩阵

C_b，C_c 叶片和机匣的阻尼矩阵

D 叶片-机匣系统的阻尼矩阵包括瑞利阻尼矩阵和科氏力矩阵

降维后的系统阻尼矩阵

E 弹性模量

e 偏心率

f_r 转频

降维后的接触力向量

F_ext 系统外力向量

叶片和机匣受的外力

降维后的系统外力向量

F^e 梁单元在o^ex^ey^ez^e的外力向量

g_N 法向接触间隙

A相对于A′的切向位移和轴向位移

h₀,h_c 叶片根部厚度；机匣厚度

i,j 梁单元的端节点

I₀ 绕z^e轴的截面惯性矩

I 描述θ_c的元素按升序排列的索引向量

k 第k个单元

节点i处的切向和法向支承刚度

在O^rX^rY^rZ^r中，三个线性支撑弹簧和两个角支撑弹簧

K 系统整体刚度矩阵

K_b，K_c 叶片和机匣的刚度矩阵

降维后的系统刚度矩阵；采用中心差分法时的等效刚度矩阵

K_c,eq 在R_c,min处的等效刚度

K^e, 结构刚度,离心刚化矩阵和旋转软化矩阵

L 叶片长度

m 叶片旋转前的障碍物数

M 系统整体质量矩阵

M_b，M_c 叶片和机匣的质量矩阵

降维后的系统质量矩阵

M^e 梁单元的质量矩阵

n 转速

n 法向接触单位矢量

n_x,n_y n的x向和y向的分量

N_v,N_u 弯曲方向和轴向的形函数矩阵

p_b,p_b,0；p_c,p_c,0 分别为叶片和机匣的气动力及其振幅

q, 分别为广义位移向量、速度向量和加速度向量

q_b,q_c 分别为叶片和机匣的广义位移向量

r_A, A的坐标、速度和加速度向量

A在t和s方向上的速度和加速度

A′的坐标、速度和加速度向量

A′相对于A在t和s方向上的速度和加速度

R_c,min 最小机匣半径

R_d 圆盘半径

s 在θ_c,m指示θ_A等于θ_i的索引号

s 轴向接触单位向量

t 时间

t 切向接触单位向量

t_x,t_y t在X向和Y向的分量

u_A,v_A,w_A； A在O^rX^rY^rZ^r中的位移，速度和加速度

A在t和s方向的位移，速度和加速度

根据λ_N修正

A′沿t向的位移，速度和加速度

u_B,v_B； B在OXY中的位移，速度和加速度

u_C,v_C； C在OXY中的位移，速度和加速度

u, 系统的位移，速度和加速度

u_b,u_c 叶片和机匣的位移矢量

u_ABC 节点A,B和C的位移，速度和加速度

x_k o^e在o^bx^by^bz^b中的x的坐标

X₀ 节点A,B和C的初始位置向量

X_A,Y_A,Z_A； A的坐标，速度和加速度分量

X_A′,Y_A′； A′的坐标，速度和加速度分量

X_B,0,X_C,0；Y_B,0,Y_C,0, 变形前B和C的X坐标和Y坐标

X_i,X_j；Y_i,Y_j 在OXY中第i个和第j个节点的X和Y坐标

希腊字母

Φ,Φ_b,Φ_c 系统，叶片和机匣的降维特征向量

Γ_i sin(θ_i)和cos(θ_i)的乘积

滑动摩擦接触面与t/s的夹角

Λ 直线Oi和X轴所夹的锐角

Θ 机匣椭圆度

α, 旋转角；旋转角速度；角加速度

β₀ 叶片安装角

β₁,β₂ 机匣单元在位置ξ处的差值位移矩阵

δ 变分号

γ B_c,t的权向量

γ_k 第k个梁单元的预扭角

γ(L) 叶片的预扭角

η B_c,t的权向量

B和C相对于Z轴的角位移，速度和加速度

κ_x,κ_y 剪切因子

λ_N, 法向，切向和轴向的拉格朗日乘子

沿s向和t向的碰摩力类型

μ 摩擦系数

υ 泊松比

θ_c 机匣所有节点的相位角向量

θ_c,m 按升序排列θ_c

θ_A,θ_B,θ_C 节点A,B和C在OXY中的相位角

θ_i 节点i在OXY中的相位角

ρ 密度

τ_b,τ_h 宽度和厚度锥度比

本发明提供一种叶片和机匣系统在碰摩故障下的力学特性分析方 法，如图1所示，包括如下步骤：

A1、对叶片进行离散化，获得叶片的梁单元模型和叶片的弹簧单元 模型；

A2、对机匣进行离散化，获得机匣的梁单元模型；

A3、基于所述叶片的梁单元模型、叶片的弹簧单元模型和机匣的梁 单元模型，构建叶片和机匣系统的有限元模型；

A4、基于接触算法，并结合所述构建叶片和机匣系统的有限元模型， 得到碰摩故障下的叶片和机匣系统的有限元模型；

A5、对所述碰摩故障下的叶片和机匣系统的有限元模型进行降维处理， 获得叶片和机匣系统的力学特性。

具体地，受安装条件、热梯度和气动载荷的影响，机匣构型往往不是圆 形的，而是呈椭圆形的。基于此，椭圆形机匣将是本发明的研究对象。单叶 片和机匣系统的示意图如图2所示，在图2中，假设圆盘内径固定，叶根与 圆盘的连接处引入五个支承弹簧来模拟叶片的弹性支撑。如图3所示，为叶 片的结构示意图，包括三个线弹簧和两个扭簧如图4 所示，为机匣的结构示意图，一系列径向弹簧和切向弹簧均匀分布在机 匣圆周上。此外，本发明采用现有技术中的Timoshenko梁单元结合弹簧单 元构建了叶片和机匣系统的有限元模型。需要说明的是，Timoshenko梁单元 矩阵是在单元坐标系o^ex^ey^ez^e中得到的(如图3所示)。因此，在构建叶片 和机匣系统的有限元模型时，需进行必要的坐标变换。本发明中，叶片的单 元矩阵是转化到旋转坐标系O^rX^rY^rZ^r中，而机匣的单元矩阵则是转化到全局 坐标系OXYZ中。此外，本发明研究中忽略了叶片扭转和机匣面外振动。

进一步地，叶片的弹簧单元模型满足如下公式：

式中：M_b、C_b、K_b、分别为叶片的质量矩阵、阻尼矩阵、刚度矩 阵和外力矩阵；u_b、分别为叶片的位移、速度和加速度。

机匣的梁单元模型满足如下公式：

式中：M_c、C_c、K_c、分别为机匣的质量矩阵、阻尼矩阵、刚度矩 阵和外力矩阵；u_c、分别为机匣的位移、速度和加速度。

综合上述，在外力F_ext的作用下，叶片和机匣系统的有限元模型满足 如下公式：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩 阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和 加速度。

进一步地，基于接触算法并结合叶片和机匣系统的有限元模型，得 到碰摩故障下的叶片和机匣系统的有限元模型。其中，接触算法包括拉 格朗日乘子法和罚函数法。

拉格朗日乘子法

接触区判断

在用拉格朗日乘子法推导含接触冲击的单叶片和机匣系统的有限元 模型之前，需首先确定叶尖点A与哪个机匣段发生接触。例如，在图5 中，θ_B≤θ_A<θ_C表示叶尖点A与机匣段BC发生接触，而A′点可以视作直 线OA在XOY面上的投影与BC的交点。接触区判断的流程图如图6所 示。图中，(θ_A，r_A)、(θ_B，r_B)和(θ_C，r_C)分别是点A，点B和点C 在OXY面内的瞬时位置信息。

接触约束矩阵

图7中n，t和s分别表示法向、周向和轴向接触单位向量；g_N为法 向接触间隙；和分别表示A″(消除浸入后的A的最终位置)相对于 A′的周向和轴向位移。

如图8所示，为叶片和机匣接触示意图。如图9、10所示，分别为 拉格朗日乘子法中切向摩擦力的预判流程图和切向摩擦力的计算流程 图。

采用拉格朗日乘子法将接触边界引入系统中，并利用哈密顿原理推 导出相应的接触约束矩阵B_c,n,B_c,t和B_c,s：

式中：u_ABC和X₀分别为点A，点B和点C的位移向量和初始位置向 量；B_c,n,B_c,t和B_c,s分别为法向、周向和轴向接触约束矩阵。

(1)法向接触约束矩阵

公式(4)中，g_N在推导B_c,n时具有非常重要的作用，因此应先计算g_N。 假设A点的径向变形为u_A，横向变形为v_A和w_A，则点A在OXY面内的 坐标表示为：

式中：e、R_d、L和α分别为偏心率、轮盘半径、叶片长度和转角。

令ξ＝BA′/BC，则A'在OXY面内的坐标表示为：

式中：N_v＝[N_v1(ξ),N_v2(ξ),N_v3(ξ),N_v4(ξ)]^T和N_u＝[N_u1(ξ),N_u2(ξ)]^T分别为弯曲和轴向的插值形函数。(n_x，n_y)和(t_x，t_y)分别是单位向量n和t在OXY 面内的分量；[X_B,0,Y_B,0]^T和[X_C,0,Y_C,0]^T分别为t＝0时点B和点C在OXY面内的 初始位置向量；为点B和点C在OXY面内的位移向量。

基于公式(5)、(6)，g_N可以表示为：

g_N＝[n_x,n_y,-n_x,-n_y][X_A,Y_A,X_A′,Y_A′]^T (7)

基于哈密顿原理对公式(7)进行变分可得：

(2)周向(t)和轴向(s)接触约束矩阵

根据公式(4)结合图7，周向和轴向接触约束矩阵与切向摩擦力密 切相关，而摩擦力的方向你取决于A′相对于A的运动。因此，需首先判 断其相对运动。

公式(5)、(6)分别对时间t求一阶和二阶导数，A和A′的速度和 加速度可以表示为：

根据公式(9)-(12)，A′相对于A的沿t和s向的速度和加速度如 下：

需要注意的是，在法向接触力存在的情况下，叶尖与机匣之间可能 存在无摩擦和摩擦两种状态。对于摩擦状态，本发明均假设为滑动摩擦。

B_c,t和B_c,s可以由下面公式计算：

式中：表示总界面摩擦力与t向间的夹角；总界面摩擦力与s向 间的夹角。

罚函数法

罚函数不需要控制叶尖对机匣的侵入量控制，其法向拉格朗日乘子 的表达式如下：

式中：

Γ₁＝Γ₀/K_c,eq,

A₀＝b₀h₀,Η＝(R_d+L)/L

由于碰摩主要发生在叶尖点和机匣最小半径处(如图2所示)，K_c,eq可以通过力除以位于最小机匣半径处节点的静变形来计算。和的计 算可以参照拉格朗日乘子法中的相关内容进行。

基于上述接触算法，并结合叶片和机匣系统的有限元模型，得到碰 摩故障下的叶片和机匣系统的有限元模型，如下式所示：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩 阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和 加速度；B_c,n,B_c,t和B_c,s分别为法向、周向和轴向接触约束矩阵；λ_N为法向拉 格朗日乘子。

式中，节点位移u与广义位移q之间的转换关系为：

将公式(18)代入公式(17)并将公式(17)左乘Φ^T可以得到：

式中：

数值仿真

叶片和机匣的仿真参数如表1所示。

表1叶片和机匣的几何和材料参数

叶片的收敛性分析

本发明中五个弹性支撑下的叶片有限元模型如图11(a)所示。图11 (b)-11(d)分别为叶片模态截断数(MTNB)对X^r-、Y^r-和Z^r-振动响 应的影响。气动力的表达式其中p_b,0＝150N，转速n＝15000rpm。从图11(b)-11(d)中可以看出，当MTNB≥10时， X^r-、Y^r-和Z^r-振动响应数值解已基本收敛，因此，在接下来的分析中均假 设MTNB＝10。图12(a)中给出了转速n∈[0,15000]rpm时，6倍转频激 励下临界共振点。在图12(b)中，绘制了n＝15000rpm下叶片的前五阶 振型。

机匣的收敛性分析

机匣的有限元模型如图13(a)所示。在机匣内表面上作用如图13 (b)中的2(ND2)下的广义压力p_c，其表达式为其 中，p_c,0＝150Pa。机匣模态截断数(MTNC)对机匣上节点1处X向和Y 向振动响应的影响如图13(c)和13(d)所示。由图可知，当MTNC＝11 时，数值解已基本收敛，因而在接下来的分析中均假设MTNC＝11。

基于两种接触算法的碰摩动力学特性

基于拉格朗日乘子法和罚函数法讨论了系统偏心和摩擦系数对系统 碰摩动力学特性的影响,相关仿真参数如表2所示。

表2仿真参数

工况1：偏心影响

偏心量e对系统位移和碰摩力响应的影响如图14和图15所示。图中 反映的一些主要结论如下：

(1)基于拉格朗日乘子法和罚函数法得到的X^r方向、Y^r方向和Z^r方向的位移幅频响应，其在振幅上存在明显的差异，尤其是在共振点处 (如A₁、A₂、A₃、A₄、A₅和A₆)的差异性更为明显，但在变化趋势上基 本一致相同(如图14所示)。此外，法向碰摩力，周向碰摩力和轴向碰 摩力的幅频响应与位移幅频响应的变化规律相似(如图15所示)。造成 这些幅值差异的原因主要在于拉格朗日乘子法需要消除叶片对机匣的侵 入量，从而导致比罚函数法更大的摩擦力和振动位移(如图16(a)所示)。 此外，拉格朗日乘子法和罚函数法在机匣上的摩擦区域也存在明显的差 异性(如图16(b)所示)。

(2)当碰摩特别严重时(如e＝1.75mm)，基于拉格朗日乘子法得 到的位移/碰摩力幅频响应会产生共振带。而基于罚函数法得到的结果不 存在类似现象，这主要是由于罚函数法得到的碰摩力远小于拉格朗日乘 子法得到的摩擦力(如图15所示)。为了进一步研究共振点A₁、A₂、A₃、A₄和共振带内的点A₅和点A₆在振动响应上的差异性(如图14所示)， 由拉格朗日乘子法得到的频谱图及其对应的频率分量分别如图17(a)、 17(b)和表3所示。从图17可以看出，组合频率的出现(见图17和表3 中的A₆)是共振带出现的必要条件。图18(a)为机匣上对应的碰摩区， 同时这也说明对于A₆工况下的碰摩不再是单点碰而是局部碰。需要指出 的是，图17中的边频fr/4是由碰摩力导致的(如图18(b)所示)。

表3不同谐振点处的转频及其倍频

注:表中所列值的单位为Hz；表中的黑体表示谐振频率。

工况2：摩擦系数的影响

摩擦系数μ对位移幅频响应和碰摩力幅频响应的影响见图19和图20 所示。图中反映的主要规律如下：

(1)如图19所示，Y^r向和Z^r向位移幅频响应表明μ的增加会导致 系统呈现出一定的软式非线性。法向、周向和轴向碰摩力的幅频响应在 变化规律上是相似的，但幅值大小上是不同的，这主要是由于拉格朗日 乘子法和罚函数法对叶片侵入量的控制不同所导致的(图16)。此外，图 20的碰摩力幅频响应表明存在两个明显的碰摩区间A和B，但其产生的 物理机制却是不一样的。碰摩区A的是由于转速靠近叶片的一阶临界转 速诱发的共振(见图12(a)中的P1)，而区域B的存在则源自于高转 速下叶片因离心力作用导致的伸长效应诱发的共振。

(2)为了进一步识别位于共振点B₁(n＝12700rpm)和B₂(n＝12950rpm) 的主导模态(图19)，叶片和机匣的频谱图如图21所示。表4列出了相 应的固有频率。在图21中，叶片频谱图中存在多处幅值放大现象(见f_r， 3f_r，5f_r(6f_r)和11f_r)，这主要是由于表4所示的叶片固有频率接近叶片固 有频率从而引发共振。在图21中，基于罚函数法得到的幅值放大现象位 于频率9f_r(1942.5Hz)处，而基于拉格朗日乘子法得到的幅值放大现象位 于频率10f_r(2116.7Hz)处。从表4中可以看出，机匣各阶模态的固有频率 均较为接近，因此，共振频率9f_r(罚函数法)和10f_r(拉格朗日乘子法) 会被误认为分别接近f_n1和f_n11。为了进一步阐明这个问题，从叶片和机匣 的振动响应中分别提取叶片和机匣各阶模态的应变能占比，结果如图22 所示。在图22的叶片各阶的模态应变能占比中，罚函数法以第3阶模态 和第10阶模态为主。但需要注意的是，叶片的第10阶模态是轴向拉伸 模态，其应变能的很大一部分来源于离心力引起的静应变能。因此，第3 阶模态才是真正的主导模态，这也与拉格朗日乘子法得到的模态识别结 果一致。在图22的机匣各阶的模态应变能占比中，罚函数法的第1阶模 态应变能(见表4中的f_n1)最大，而第6阶模态次之，这与拉格朗日乘子法 得到的结果是相反的。因此，对于拉格朗日乘子法来说，共振频率10f_r (2116.7Hz)更接近f_n6，而不是f_n11，相应的振型图如图13(b)所示。总 的来说，基于罚函数法和拉格朗日乘子法得到的机匣各阶模态应变能的 分布是一致的(图22)。

表4叶片和机匣的固有频率

本发明采用自编弹簧单元和铁木辛柯梁单元，建立了柔性支承叶片 和机匣的有限元模型。在不考虑叶片扭转和机匣面外运动的情况下，采 用自编的点面接触算法，结合向前增量拉格朗日乘子法和向前增量罚函 数法，得到了系统的碰摩动力学特性。同时，为了提高计算效率，在求 解时中采用了模型降维技术。基于本发明所建立的动力学模型，详细讨 论了偏心量和摩擦系数对系统碰摩特性的影响。相关结论如下：

(1)与罚函数法相比，拉格朗日乘子法可以较好地消除叶尖与机匣 之间的浸入量，但会导致较大的碰摩力和振动位移；

(2)碰摩严重程度的增加使得局部共振点演变成共振带甚至导致边 频带的出现；

(3)摩擦系数的增大使局部共振点左移，呈现出软式非线性现象。 此外，应变能结合频谱图能对系统中的主导模态进行较好的识别。

以上结合具体实施例描述了本发明的技术原理，这些描述只是为了 解释本发明的原理，不能以任何方式解释为对本发明保护范围的限制。 基于此处解释，本领域的技术人员不需要付出创造性的劳动即可联想到 本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之内。

Claims (9)

1.一种叶片和机匣系统在碰摩故障下的力学特性分析方法，其特征在于，包括如下步骤：

A1、对叶片进行离散化，获得叶片的梁单元模型和叶片的弹簧单元模型；

A2、对机匣进行离散化，获得机匣的梁单元模型；

A3、基于所述叶片的梁单元模型、叶片的弹簧单元模型和机匣的梁单元模型，构建叶片和机匣系统的有限元模型；

A4、基于接触算法，并结合所述构建叶片和机匣系统的有限元模型，得到碰摩故障下的叶片和机匣系统的有限元模型；

A5、对所述碰摩故障下的叶片和机匣系统的有限元模型进行降维处理，获得叶片和机匣系统的力学特性。

2.根据权利要求1所述的力学特性分析方法，其特征在于，所述叶片的弹簧单元模型包括三个线弹簧和两个扭簧

3.根据权利要求2所述的力学特性分析方法，其特征在于，所述叶片的弹簧单元模型满足如下公式：

式中：M_b、C_b、K_b、分别为叶片的质量矩阵、阻尼矩阵、刚度矩阵和外力矩阵；u_b、分别为叶片的位移、速度和加速度。

4.根据权利要求3所述的力学特性分析方法，其特征在于，所述机匣的梁单元模型满足如下公式：

式中：M_c、C_c、K_c、分别为机匣的质量矩阵、阻尼矩阵、刚度矩阵和外力矩阵；u_c、分别为机匣的位移、速度和加速度。

5.根据权利要求4所述的力学特性分析方法，其特征在于，所述步骤A3中，叶片和机匣系统的有限元模型满足如下公式：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和加速度。

6.根据权利要求5所述的力学特性分析方法，其特征在于，所述步骤A4中，接触算法为改进的拉格朗日乘子法。

7.根据权利要求6所述的力学特性分析方法，其特征在于，所述改进的拉格朗日乘子法中，包括如下步骤：

S1、将参数θ_C作为输入数据，输入预先训练的第一模型，获取接触区域；

所述第一模型为基于历史参数θ与对应的接触区域训练而获得的接触区域判定模型；

S2、在所述接触区域内，将预先获得的参数μ、g_N作为输入数据，输入预先训练的第二模型，获取沿s向和t向的碰摩力类型和滑动摩擦接触面与s/t向的夹角；

其中，μ为摩擦系数，分别为叶尖点A在t向的位移、速度和加速度；分别为叶尖点A在s向的位移、速度和加速度；分别为A’在t向的位移、速度和加速度；g_N为法向接触间隙；

所述第二模型为基于历史参数μ、u^t、u^s、u_' ^t、g_N与对应的沿s向和t向的碰摩力类型训练而获得的碰摩力类型判定模型；

S3、将获得的沿s向和t向的碰摩力类型、滑动摩擦接触面与s/t向的夹角和预先获得的作为输入数据，输入预先训练的第三模型，获取拉格朗日乘子；

其中，分别为A’相对于A在t向的速度、加速度和A’相对于A在s向的速度、加速度；

所述第三模型为基于历史参数沿s向和t向的碰摩力类型、滑动摩擦接触面与s/t向的夹角、训练而获得的切向摩擦力判定模型。

8.根据权利要求7所述的力学特性分析方法，其特征在于，所述步骤A5中，基于对叶片和机匣的收敛性分析，对碰摩故障下的叶片和机匣系统的有限元模型进行降维处理。

9.根据权利要求8所述的力学特性分析方法，其特征在于，所述碰摩故障下的叶片和机匣系统的有限元模型满足如下公式：

式中：M、D、K、分别为叶片和机匣系统的质量矩阵、阻尼矩阵、刚度矩阵和外力矩阵；u、分别为叶片和机匣系统的位移、速度和加速度；B_c,n,B_c,t和B_c,s分别为法向、周向和轴向接触约束矩阵；λ_N为法向拉格朗日乘子。