一种转子-叶片耦合系统固有频率的确定方法

摘要

本发明一种转子-叶片耦合系统固有频率的确定方法，属于机械动力学技术领域，本发明节省了实验测试需要的传感器、放大器以及显示或记录仪表所用的成本费用；能够获得较高阶次的固有频率以及旋转状态下的固有频率；本发明无需重复建模，仅需修改系统的结构尺寸后即可得到不同转子-叶片耦合系统的固有频率；本发明考虑了耦合系统复杂阶梯转轴的弯曲和扭转的影响、轴和盘的陀螺效应，以及离心刚化、旋转软化和叶片的科氏力影响，能够得到更加准确的结果；此外，本发明还能进行转子-叶片耦合系统不平衡响应、叶尖碰摩等故障分析，从而实现系统结构的优化。

说明书

技术领域

本发明属于机械动力学技术领域，具体涉及一种转子-叶片耦合系统固有频率的确定方 法。

背景技术

目前，现有的转子-叶片耦合系统的固有频率确定方法主要有以下几种途径：

1.基于实验测试

固有频率测试就是通过传感器、放大仪器以及显示或记录仪表，获得固有频率等特性参 数，为机械或工程结构的动力设计服务。

力锤激励法具有测试方便、所要求的仪器设备简单等特点，已成为测试固有频率的首选 方法；利用力锤激励法，可由响应信号和锤头所激励的宽频力信号的比值，即频响函数来辨 识结构的各阶的固有频率和模态振型；力锤激励法实验方便，花费少，但其能量有限且不易 控制，对于大型结构，往往因能量不足造成远端响应太小，使测量信噪比小，造成较大误差； 不同传感器位置及激励点的位置对测试精度也有较大影响；且实验获得的固有频率有限，不 易获得系统的高阶固有频率以及旋转状态下的固有频率。

2.基于大型商业有限元分析软件

将CAD三维模型导入有限元分析软件中或者在有限元分析软件中建立的三维模型；并选择 合适的单元及各零部件的材料参数，对三维模型进行网格划分，建立有限元模型；设置约束 及添加载荷后选择适宜的方法进行模态分析，得到系统结构的固有频率；但在利用成熟的有 限元分析软件对较复杂结构进行动力特性分析时，建模(软件的前处理部分)过程非常复杂且 繁重，而且对于具有相类似结构，需要进行近乎重复的、繁重的建模过程，不同方式的建模得 到的固有频率也会有较大差距。

3.基于理论计算

通过建立转子-叶片耦合系统的动力学模型，从理论上推导系统的固有频率；然而现有的 转子-叶片的耦合动力学模型主要考虑了转轴的扭转和叶片的弯曲变形的影响，且现有的动力 学模型仍旧适用于光轴这种较为理想的系统结构，对于复杂阶梯轴的横向、扭转运动以及叶 片的耦合影响的技术处于空白状态。

发明内容

针对现有技术的缺点，本发明提出一种转子-叶片耦合系统固有频率的确定方法，以达到 在考虑转轴的弯曲和扭转，转轴和盘的陀螺效应，叶片的离心刚化、旋转软化和科氏力影响 的情况下，准确确定旋转态下转子-叶片耦合系统固有频率的目的。

一种转子-叶片耦合系统固有频率的确定方法，包括以下步骤：

步骤1、构建转子-叶片耦合系统所需三维坐标系，包括：整体坐标系、第i个叶片局部 坐标系、圆盘坐标系和第i个旋转坐标系；

具体如下：

整体坐标系OXYZ：以静止态下转子-叶片耦合系统的圆盘中心点为原点构建整体坐标系， 该坐标Z轴平行于转子-叶片耦合系统的转轴方向；

第i个叶片局部坐标系ox^by^bz^b：以转子-叶片耦合系统第i个叶片根部中心点为原点建立 叶片局部坐标系，该坐标的x^b轴沿着叶片长度方向，y^b轴沿着叶片厚度方向，z^b轴沿着叶片 宽度方向；

圆盘坐标系ox^dy^dz^d：以旋转态下转子-叶片耦合系统的圆盘中心点为原点建立圆盘坐标 系，该坐标z^d轴垂直于圆盘，该坐标的x^d、y^d轴分别平行于整体坐标系的X、Y轴；

第i个旋转坐标系ox^ry^rz^r：以旋转态下转子-叶片耦合系统的圆盘中心点为原点建立旋转 坐标系，该坐标z^r轴垂直于圆盘，该坐标x^r轴与圆盘坐标系x^d轴的夹角为

其中，θ(t)为圆盘运动的角位移，表示第i个叶片在叶片组中 的位置；N_b为叶片数；

步骤2、在转子-叶片耦合系统的旋转态下，确定系统的总动能，具体如下：

步骤2-1、确定转子-叶片耦合系统中第i个叶片的动能，具体步骤如下：

步骤2-1-1、根据圆盘处转轴的扭转角，确定第i个叶片局部坐标系与第i个旋转坐标系 的转换矩阵A₁；

步骤2-1-2、根据圆盘运动的角位移和每个叶片在所有叶片中的位置，确定第i个旋转坐 标系与圆盘坐标系的转换矩阵A₂；

步骤2-1-3、在转子-叶片耦合系统旋转态下，根据旋转时圆盘所产生的位移、圆盘半径 和叶片的剪切角，确定第i个叶片上任意一点Q在整体坐标系中的位移向量；

计算公式如下：

其中，x表示Q点在第i个叶片局部坐标系ox^by^bz^b中沿叶片长度方向的位置；y表示Q 点在第i个叶片局部坐标系ox^by^bz^b中沿叶片厚度方向的位置；表示第i个叶片局部坐标系 ox^by^bz^b中叶片剪切角；x_d表示圆盘旋转时在整体坐标系OXYZ中的X方向位移；y_d表示圆盘 旋转时在整体坐标系OXYZ中的Y方向位移；z_d表示圆盘旋转时在整体坐标系OXYZ中的Z 方向位移；R_d表示圆盘的半径；u为Q点在第i个叶片局部坐标系ox^by^bz^b中叶片长度方向位 移；v为Q点在第i个叶片局部坐标系ox^by^bz^b中叶片厚度方向位移；w为Q点在第i个叶片 局部坐标系ox^by^bz^b中叶片宽度方向位移；

步骤2-1-4、根据叶片的横截面积、叶片密度、叶片长度和Q点的速度，获得叶片动能；

第i个叶片动能T_blade计算公式如下：

$T_{\mathrm{blade}}=\frac{1}{2}{\int }_0^L{\rho }_b{A}_b{\stackrel{.}{r}}_Q^2\mathrm{dx}=\frac{1}{2}(T_1+T_2)---\left(2\right)$

其中，A_b表示叶片的横截面积；ρ_b表示叶片密度；L表示叶片长度；为第i个叶片上 任意一点Q的速度；T₁表示第i个叶片振动产生的动能项；T₂表示第i个叶片与转子耦合振 动产生的动能项；T₁的表达式如下：

公式(3)中，表示圆盘运动的角速度；I_b表示叶片任意位置处的截面惯性矩；表示v 对时间t的一阶导数；表示u对时间t的一阶导数；表示对时间t的一阶导数；

T₂的表达式如下：

公式(4)中，表示Ψ对时间t的一阶导数；

步骤2-2、确定转子-叶片耦合系统中转轴的动能，具体步骤如下：

步骤2-2-1、根据圆盘质心与圆盘形心位移之间的关系，获得圆盘质心在整体坐标系中X 方向的位移x_c和Y方向的位移y_c；

步骤2-2-2、根据获得的圆盘质心位移，结合圆盘的极转动惯量、轴段的极转动惯量、圆 盘的直径转动惯量、轴段的直径转动惯量、圆盘在整体坐标系中的摆角和轴段在整体坐标系 中的摆角，确定系统中转轴的动能；

转轴的动能T_rotor计算公式如下：

$\left(\begin{array}{c}{T}_{\mathrm{rotor}}=\frac{1}{2}{J}_{\mathrm{pd}}{(\stackrel{.}{\theta }+\stackrel{.}{\Psi })}^2+\frac{1}{2}{m}_d({\stackrel{.}{x}}_c^2+{\stackrel{.}{y}}_c^2)-J_{\mathrm{pd}}(\stackrel{.}{\theta }+\stackrel{.}{\Psi }){\stackrel{.}{\theta }}_{\mathrm{yd}}{\theta }_{\mathrm{xd}}+\frac{1}{2}{J}_{\mathrm{dd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}^2+{\stackrel{.}{\theta }}_{\mathrm{yd}}^2)\\ +\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{m}_j({\stackrel{.}{x}}_j^2+{\stackrel{.}{y}}_j^2)-\underset{j=1}{\overset{N_d}{\Sigma }}{J}_{\mathrm{pj}}\stackrel{.}{\theta }{\stackrel{.}{\theta }}_{\mathrm{yj}}{\theta }_{\mathrm{xj}}+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{J}_{\mathrm{dj}}({\stackrel{.}{\theta }}_{\mathrm{xj}}^2+{\stackrel{.}{\theta }}_{\mathrm{yj}}^2)\end{array}\right)---\left(5\right)$

式中，表示圆盘质心在整体坐标系中X方向位移x_c对时间t的一阶导数；表示圆盘 质心在整体坐标系中Y方向位移y_c对时间t的一阶导数；J_pd表示圆盘处集中质量点的极转动 惯量；J_dd表示圆盘处集中质量点的直径转动惯量；J_pj表示第j个集中质量点的极转动惯量； J_dj表示第j个集中质量点的直径转动惯量；θ_xd表示圆盘绕整体坐标系X轴的摆角；θ_yd表示圆 盘绕整体坐标系Y轴的摆角；m_d表示圆盘的质量；m_j表示转轴上的第j个集中质量点的质 量；x_j表示转轴上的第j个集中质量点在整体坐标系中X方向位移；y_j表示转轴上的第j个集 中质量点在整体坐标系中Y方向位移；θ_xj表示转轴上的第j个集中质量点绕整体坐标系X轴 的转角；θ_yj表示转轴上的第j个集中质量点绕整体坐标系Y轴的转角；表示θ_xd对时间t的 一阶导数；表示θ_yd对时间t的一阶导数；表示θ_xj对时间t的一阶导数；表示θ_yj对时 间t的一阶导数；N_d表示集中质量点数；j≠p_disc，p_disc表示圆盘处集中质量点的编号；

根据圆盘质心和形心之间的位置关系，获得转轴动能最终计算公式如下：

$\left(\begin{array}{c}{T}_{\mathrm{rotor}}=\frac{1}{2}{J}_{\mathrm{pd}}{(\stackrel{.}{\theta }+\stackrel{.}{\Psi })}^2+\frac{1}{2}{m}_d({({\stackrel{.}{x}}_d-e\sin (\theta +\Psi )(\stackrel{.}{\theta }+\stackrel{.}{\Psi }))}^2+{({\stackrel{.}{y}}_d+e\cos (\theta +\Psi )(\stackrel{.}{\theta }+\stackrel{.}{\Psi }))}^2)\\ -(\stackrel{.}{\theta }+\stackrel{.}{\Psi })J_{\mathrm{pd}}{\stackrel{.}{\theta }}_{\mathrm{yd}}{\theta }_{\mathrm{xd}}+\frac{1}{2}{J}_{\mathrm{dd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}^2+{\stackrel{.}{\theta }}_{\mathrm{yd}}^2)+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{m}_j({\stackrel{.}{x}}_j^2+{\stackrel{.}{y}}_j^2)-\underset{j=1}{\overset{N_d}{\Sigma }}\stackrel{.}{\theta }{J}_{\mathrm{pj}}{\stackrel{.}{\theta }}_{\mathrm{yj}}{\theta }_{\mathrm{xj}}+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{J}_{\mathrm{dj}}{({\stackrel{.}{\theta }}_{\mathrm{xj}}+{\stackrel{.}{\theta }}_{\mathrm{yj}})}^2\end{array}\right)---\left(6\right)$

公式(6)中，e表示圆盘质心与形心不对中时的偏心距；

步骤2-3、将转子-叶片耦合系统的叶片动能与转轴动能进行求和，获得系统总的动能；

步骤3、在转子-叶片耦合系统旋转态下，确定系统的总势能，具体如下：

步骤3-1、根据叶片的弹性模量、叶片的剪切模量、剪切修正系数、离心力和法向力，确 定转子-叶片耦合系统中第i个叶片的势能；

第i个叶片的势能V_blade计算公式如下：

其中，E_b表示叶片的弹性模量；G_b表示叶片的剪切模量；κ表示叶片的剪切修正系数； f_c表示叶片的离心力；F_n表示叶尖受到的法向力；

步骤3-2、根据圆盘处转轴的扭转刚度和转轴的各集中质量点位移向量，确定转子-叶片 耦合系统中转轴的势能；

转轴的势能计算公式如下：

$V_{\mathrm{rotor}}=\frac{1}{2}{k}_{\Psi }{\Psi }^2+\frac{1}{2}{q}_r^T{K}_r{q}_r---\left(8\right)$

其中，k_Ψ表示圆盘处转轴的扭转刚度；q_r表示转轴的各集中质量点位移向量， q_r＝[…x_j θ_yj……y_j θ_xj…]^T；K_r表示转轴刚度矩阵；

步骤3-3、将转子-叶片耦合系统的叶片势能与转轴势能进行求和，获得系统总的势能；

步骤4、根据获得的转子-叶片耦合系统的总的动能和总的势能，结合哈密顿原理通过变 分运算确定转子-叶片耦合系统的运动情况，具体如下：

步骤4-1、根据获得的系统总的动能和势能，结合哈密顿原理进行变分运算，确定第i个 叶片的运动微分方程；

步骤4-2、根据哈密顿原理进行变分运算，确定圆盘位置处的扭转运动微分方程；

步骤4-3、根据哈密顿原理进行变分运算，确定转轴的横向运动微分方程；

步骤4-4、根据正则坐标对第i个叶片的运动微分方程中的叶片的长度方向位移、厚度方 向位移和剪切角进行离散化处理，获得第i个叶片的质量矩阵、阻尼矩阵和刚度矩阵；根据 正则坐标对圆盘位置处的扭转运动微分方程中的叶片的长度方向位移、厚度方向位移和剪切 角进行离散化处理，获得系统质量矩阵耦合项、阻尼矩阵耦合项和刚度矩阵耦合项；

步骤4-5、对叶片的质量矩阵、叶片科氏力矩阵、叶片刚度矩阵、质量矩阵耦合项、阻尼 矩阵耦合项、刚度矩阵耦合项、圆盘位置处转轴的扭转质量、圆盘位置处转轴的扭转阻尼、 圆盘位置处转轴的扭转刚度、转轴的质量矩阵、转轴的陀螺矩阵和转轴的刚度矩阵进行组集， 获得转子-叶片耦合系统的质量矩阵、阻尼矩阵和刚度矩阵，结合转子-叶片耦合系统的广义 坐标向量和外激振力向量构建转子-叶片耦合系统的运动微分方程；

步骤5、设定外激振力向量为零，根据获得的转子-叶片耦合系统的质量矩阵、阻尼矩阵 和刚度矩阵，确定转子-叶片耦合系统的固有频率；

步骤6、根据获得的转子-叶片耦合系统的固有频率，确定系统的工作转速，避免共振， 使系统运行稳定。

步骤2-1-4所述的叶片的横截面积和叶片任意位置处的截面惯性矩，计算公式如下：

$\left(\begin{array}{c}{A}_b=A_0(1-{\tau }_b\frac{x}{L})(1-{\tau }_h\frac{x}{L}),{\tau }_b=1-\frac{b_1}{b_0}\\ I_b=I_0(1-{\tau }_b\frac{x}{L}){(1-{\tau }_h\frac{x}{L})}^3,{\tau }_h=1-\frac{h_1}{h_0}\end{array}\right)---\left(9\right)$

其中，b₁表示叶尖处的叶片宽度；h₁表示叶尖处的叶片厚度；b₀表示叶根处的叶片宽度； h₀表示叶根处的叶片厚度；A₀表示叶根处的横截面积，A₀＝b₀h₀；I₀表示叶根处的惯性矩， $I_0=\frac{1}{12}{b}_0{h}_0^3.$

步骤4-1所述的确定每个叶片的运动微分方程，具体如下：

将系统总动能和总势能，代入哈密顿原理表达式(10)中：

$\delta {\int }_{t_1}^{t_2}(T_{\mathrm{total}}-V_{\mathrm{total}})\mathrm{dt}=0---\left(10\right)$

其中，T_total表示转子-叶片耦合系统总动能；V_total表示转子-叶片耦合系统总的势能；t₁表 示初始时间；t₂表示终止时间；

以δu为独立变量进行变分运算，获得第i个叶片的第一个运动微分方程为：

其中，ü表示u对时间t的二阶导数；表示θ对时间t的二阶导数；A′_b表示A_b对积分 域x的一阶导数；u′表示u对积分域x的一阶导数；u″表示u对积分域x的二阶导数；表 示Ψ对时间t的二阶导数；表示x_d对时间t的二阶导数；表示y_d对时间t的二阶导数；

以δv为独立变量进行变分运算，获得第i个叶片的第二个运动微分方程为：

其中，表示v对时间t的二阶导数；f′_c(x)表示叶片离心力对积分域x的一阶导数；v′表 示v对积分域x的一阶导数；v″表示v对积分域x的二阶导数；表示对积分域x的一阶 导数；F′_n表示叶片的法向力对积分域x的一阶导数；

以为独立变量进行变分运算，获得第i个叶片的第三个运动微分方程为：

其中，表示对时间t的二阶导数；表示对积分域x的二阶导数；I′_b表示I_b对积 分域x的一阶导数。

步骤4-2所述的确定圆盘位置处的扭转运动微分方程，具体如下：

以δΨ为独立变量进行变分运算，获得圆盘位置处的扭转运动微分方程为：

$\left(\begin{array}{c}(J_{\mathrm{pd}}+e^2{m}_d+C_1)\stackrel{..}{\Psi }+C_2\stackrel{.}{\Psi }+(k_{\Psi }-e^2{m}_d{\stackrel{.}{\theta }}^2+C_3)\Psi \\ -{\mathrm{em}}_d\sin \theta {\stackrel{..}{x}}_d+{\mathrm{em}}_d\cos \theta {\stackrel{..}{y}}_d+e^2{m}_d\stackrel{..}{\theta }-J_{\mathrm{pd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}{\stackrel{.}{\theta }}_{\mathrm{yd}}+{\theta }_{\mathrm{xd}}{\stackrel{..}{\theta }}_{\mathrm{yd}})\\ +\stackrel{..}{\theta }{J}_{\mathrm{pd}}+C_4-C_5-C_6=0\end{array}\right)---\left(14\right)$

其中，

步骤4-3所述的确定转轴的横向运动微分方程，具体如下：

以δq_r为独立变量进行变分运算，获得转轴的横向运动微分方程为：

$M_r{\stackrel{..}{q}}_r+(G_r+C_b){\stackrel{.}{q}}_r+K_r{q}_r=F_r---\left(15\right)$

其中，M_r表示转轴的质量矩阵；G_r表示转轴的陀螺矩阵；F_r表示转轴的外激振力向量； C_b为转轴轴承的阻尼矩阵。

步骤4-4所述的叶片的长度方向位移、厚度方向位移和剪切角进行离散化处理，具体如 下：

公式(16)中，U_ξ(t)表示离散处理后叶片的长度方向位移；V_ξ(t)表示离散处理后叶片的 厚度方向位移；ψ_ξ(t)表示离散处理后叶片的剪切角；φ_1ξ(x)表示叶片长度方向振动的第ξ阶振 型函数；φ_2ξ(x)表示叶片厚度方向振动的第ξ阶振型函数；φ_3ξ(x)表示叶片剪切角的第ξ阶振 型函数；

$\left(\begin{array}{c}{\phi }_{1\xi }\left(x\right)=\frac{\sin \left({\alpha }_{\xi }x\right)}{{\alpha }_{\xi }}\\ {\phi }_{2\xi }\left(x\right)=\frac{1-\cos \left({\alpha }_{\xi }x\right)}{{\alpha }_{\xi }}\\ {\phi }_{3\xi }\left(x\right)=\sin \left({\alpha }_{\xi }x\right)\end{array}\right)---\left(17\right)$

公式(17)中，ξ＝1，2，3，…，N_mod，L为叶片的长度，N_mod为截断的模态数。

步骤4-5所述的构建转子-叶片耦合系统的运动微分方程，具体如下：

转子-叶片耦合系统的运动微分方程如下：

$M_{\mathrm{RB}}{\stackrel{..}{q}}_{\mathrm{RB}}+D_{\mathrm{RB}}{\stackrel{.}{q}}_{\mathrm{RB}}+K_{\mathrm{RB}}{q}_{\mathrm{RB}}=f_{\mathrm{RB}}---\left(18\right)$

其中，M_RB表示转子-叶片耦合系统的质量矩阵，

$M_{\mathrm{RB}}=\left(\begin{array}{ccc}{M}_b& M_{c1}& M_{c2}\\ M_{c1}^T& M_{\Psi }& 0\\ M_{c2}^T& 0& M_r\end{array}\right)---\left(19\right)$

公式(19)中，M_b表示叶片的质量矩阵；M_c1为扭转运动的质量耦合项；M_c2为横 向运动的质量耦合项；M_Ψ为圆盘位置处转轴的扭转质量；M_r转轴的质量矩阵；

D_RB表示转子叶片耦合系统的阻尼矩阵；

D_RB的表达式如下：

D_RB＝G_RB+D_b       (20)

公式(20)中，D_b为系统轴承的阻尼矩阵；

公式(20)中，G_RB的表达式如下：

$G_{\mathrm{RB}}=\left(\begin{array}{ccc}{G}_b& G_{c1}& 0\\ -G_{c1}^T& G_{\Psi }& 0\\ 0& 0& G_r\end{array}\right)---\left(21\right)$

公式(21)中，G_b表示叶片的科氏力矩阵；G_c1表示阻尼耦合项；G_Ψ表示圆盘位置处转 轴的扭转阻尼；G_r表示系统转轴的陀螺矩阵；

K_RB表示转子-叶片耦合系统的刚度矩阵；

$K_{\mathrm{RB}}=\left(\begin{array}{ccc}{K}_b& 0& 0\\ 0& K_{\Psi }& 0\\ 0& 0& K_r\end{array}\right)---\left(22\right)$

公式(22)中，K_b表示叶片刚度矩阵；K_Ψ表示圆盘位置处转轴的扭转刚度；K_r表示转 轴的刚度矩阵；

q_RB表示转子-叶片耦合系统的广义坐标向量，q_b表示叶片的广义坐标 向量；表示q_RB对时间t的二阶导数；表示q_RB对时间t的一阶导数；

f_RB表示转子-叶片耦合系统的外激振力向量。

步骤5所述的确定转子-叶片耦合系统的固有频率，具体如下：

步骤5-1、设置外激振力向量f_RB为零，构建特征方程；

特征方程如下：

$\left(\begin{array}{cc}{D_{\mathrm{RB}}}^T& {M_{\mathrm{RB}}}^T\\ {M_{\mathrm{RB}}}^T& 0\end{array}\right)\left(\begin{array}{c}{\stackrel{.}{q}}_{\mathrm{RB}}\\ {\stackrel{..}{q}}_{\mathrm{RB}}\end{array}\right)+\left(\begin{array}{cc}{K_{\mathrm{RB}}}^T& 0\\ 0& -{M_{\mathrm{RB}}}^T\end{array}\right)\left(\begin{array}{c}{q}_{\mathrm{RB}}\\ {\stackrel{.}{q}}_{\mathrm{RB}}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)---\left(23\right)$

步骤5-2、对所构建的特征方程进行简化；

即设置如下关系：

$A_{\mathrm{mod}}=\left(\begin{array}{cc}{K_{\mathrm{RB}}}^T& 0\\ 0& -{M_{\mathrm{RB}}}^T\end{array}\right),B_{\mathrm{mod}}=\left(\begin{array}{cc}{D_{\mathrm{RB}}}^T& {M_{\mathrm{RB}}}^T\\ {M_{\mathrm{RB}}}^T& 0\end{array}\right),Z_{\mathrm{mod}}=\left(\begin{array}{c}{q}_{\mathrm{RB}}\\ {\stackrel{.}{q}}_{\mathrm{RB}}\end{array}\right)---\left(24\right)$

将A_mod、B_mod、Z_mod代入公式(23)，获得如下方程：

$B_{\mathrm{mod}}{\stackrel{.}{Z}}_{\mathrm{mod}}+A_{\mathrm{mod}}{Z}_{\mathrm{mod}}=0---\left(25\right)$

并设置Z_mod＝Ze^-λt，代入到式(25)中化简获得如下方程：

(A_mod-λB_mod)Z＝0     (26)

公式(26)中，λ表示特征方程的特征值；Z表示特征方程的特征向量；

步骤5-3、将化简后方程中的系数行列式设置为零；

即：

|A_mod-λB_mod|＝0      (27)

步骤5-4、根据该系数行列式的特征值，确定转子-叶片耦合系统的固有频率；

即计算获得系数行列式的一组特征值λ，取其虚部的绝对值除以2π，并进行从小到大排 序，获得一组固有频率ω_k，其中，k表示转子-叶片耦合系统模态的第k阶，k＝1，2，…。

本发明优点：

本发明为一种转子-叶片耦合系统固有频率的确定方法，节省了实验测试需要的传感器、 放大器以及显示或记录仪表所用的成本费用；能够获得较高阶次的固有频率以及旋转状态下 的固有频率；本发明无需重复建模，仅需修改系统的结构尺寸后即可得到不同转子-叶片耦合 系统的固有频率；本发明考虑了耦合系统复杂阶梯转轴的弯曲和扭转的影响、轴和盘的陀螺 效应，以及离心刚化、旋转软化和叶片的科氏力影响，能够得到更加准确的结果；此外，本 发明还能进行转子-叶片耦合系统不平衡响应、叶尖碰摩等故障分析，从而实现系统结构的优 化。

附图说明

图1为本发明一种实施例的转子-叶片试验台结构示意图；

图2为本发明一种实施例的转子-叶片耦合系统固有频率的确定方法流程图；

图3为本发明一种实施例的转子-叶片耦合系统示意图；

图4为本发明一种实施例的盘片系统示意图；

图5为本发明一种实施例的圆盘的瞬时位置示意图；

图6为本发明一种实施例的转子-叶片耦合系统矩阵组集示意图，其中，图(a)为质量矩 阵，图(b)为阻尼矩阵，图(c)为刚度矩阵；

图7为本发明一种实施例的转子-盘片耦合系统的两个有限元模型，其中，图(a)为有限 元模型1示意图，图(b)为有限元模型2示意图；

图8为本发明一种实施例的转子叶片耦合系统在4、6、8叶片情况下的Campbell图，其 中，图(a)为4叶片情况示意图，图(b)为6叶片情况示意图，图(c)为8叶片情况示意图。

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。

本发明实施例中系统实验台如图1所示，1表示电机(型号YVP100L2-3)、2表示弹性 联轴器(型号YLD7)、3表示左轴承、4表示右轴承、5表示阶梯轴、6表示轮盘和7表示叶 片。轴的几何参数，如表1所示。盘的半径和厚度分别为140mm和58mm。叶片的长度、 宽度和厚度分别为82mm、44mm和3mm。转子和叶片的材料参数均为：杨氏模量E＝200GPa、 密度ρ＝7800kg/m³和泊松比υ＝0.3。两个滚珠轴承完全相同，均通过线性弹簧和阻尼来模拟， 忽略交叉耦合项，水平和竖直方向刚度假定为k_bx1＝k_by1＝k_bx2＝k_by2＝1.5×10⁷N/m。

本发明实施例基于表1所示数据，转轴划分为15个集中质量单点，圆盘所处的集中质量 点编号为9。为了研究系统的无阻尼固有频率，忽略了转子-叶片耦合系统轴承的阻尼矩阵D_b， 在计算固有频率时，忽略转轴扭转和叶片变形的高阶耦合项。

其轴段具体参数如表1所示：

表1

本发明实施例中转子-叶片耦合系统固有频率的确定方法，方法流程图如图2所示，包括 以下步骤：

步骤1、构建转子-叶片耦合系统所需三维坐标系，包括：整体坐标系、第i个叶片局部 坐标系、圆盘坐标系和第i个旋转坐标系；

本发明实施例中，考虑转子的弯曲和扭转变形、叶片的纵向和横向振动的转子-叶片耦合 系统示意图如图3所示，图中OXYZ为整体坐标系、ox^dy^dz^d为圆盘的坐标系，由于转子在运 动中产生涡动，所以圆盘的坐标系与整体坐标系并不重合；ox^ry^rz^r为旋转坐标系；ox^by^bz^b为 叶片的局部坐标系；坐标系建立具体如下：

整体坐标系OXYZ：以静止态下转子-叶片耦合系统的圆盘中心点为原点构建整体坐标系， 该坐标Z轴平行于转子-叶片耦合系统的转轴方向；

第i个叶片局部坐标系ox^by^bz^b：以转子-叶片耦合系统第i个叶片根部中心点为原点建立 叶片局部坐标系，该坐标的x^b轴沿着叶片长度方向，y^b轴沿着叶片厚度方向，z^b轴沿着叶片 宽度方向；

圆盘坐标系ox^dy^dz^d：以旋转态下转子-叶片耦合系统的圆盘中心点为原点建立圆盘坐标 系，该坐标z^d轴垂直于圆盘，该坐标的x^d、y^d轴分别平行于整体坐标系的X、Y轴；

第i个旋转坐标系ox^ry^rz^r：以旋转态下转子-叶片耦合系统的圆盘中心点为原点建立旋转 坐标系，该坐标z^r轴垂直于圆盘，如图4所示，该坐标x^r轴与圆盘坐标系x^d轴的夹角为； 其中，θ(t)为圆盘运动的角位移，表示第i个叶片在叶片组中的位 置；N_b为叶片数；

图3中，x_b1表示左轴承水平方向、y_b1表示左轴承竖直方向、m_b1表示左轴承质量、c_bx1表 示左轴承水平方向阻尼、c_by1表示左轴承竖直方向阻尼、k_bx1表示左轴承水平方向刚度、k_by1表 示左轴承竖直方向刚度、x_b2表示右轴承水平方向、y_b2表示右轴承竖直方向、m_b2表示右轴承 轴承质量、c_bx2表示右轴承水平方向阻尼、c_by2表示右轴承竖直方向阻尼、k_bx2表示右轴承水 平方向刚度、k_by2表示右轴承竖直方向刚度。

本发明实施例中，对系统做如下假设：

(1)材料各向同性，本构关系满足Hooke定律；

(2)忽略叶-盘-轴之间的接触问题，认为三者之间是固定连接的；

(3)圆盘简化为刚性盘，不考虑圆盘的柔性变形；

(4)轴和刚性盘均采用质点模型，其轴向变形很小，可以忽略不计；

(5)轴承部分采用线性刚度和阻尼模型(弹簧阻尼模型)。

步骤2、在转子-叶片耦合系统的旋转态下，确定系统的总动能，具体如下：

步骤2-1、确定转子-叶片耦合系统中第i个叶片的动能，具体步骤如下：

步骤2-1-1、根据圆盘处转轴的扭转角，确定第i个叶片局部坐标系与第i个旋转坐标系 的转换关系：

根据图4确定第i个叶片局部坐标系与第i个旋转坐标系的转换矩阵A₁如下：

$A_1=\left(\begin{array}{ccc}\cos \Psi & -\sin \Psi & 0\\ \sin \Psi & \cos \Psi & 0\\ 0& 0& 1\end{array}\right)=\left(\begin{array}{ccc}1& -\Psi & 0\\ \Psi & 1& 0\\ 0& 0& 1\end{array}\right)---\left(28\right)$

其中，Ψ为圆盘处转轴的扭转角；

步骤2-1-2、根据圆盘运动的角位移和第i个叶片在所有叶片中的位置，确定第i个旋转 坐标系与圆盘坐标系的转换关系：

根据图4确定第i个旋转坐标系与圆盘坐标系的转换矩阵A₂如下：

步骤2-1-3、在转子-叶片耦合系统旋转态下，根据旋转时圆盘所产生的位移、圆盘半径 和叶片的剪切角，确定第i个叶片上任意一点Q在整体坐标系中的位移向量；

计算公式如下：

其中，x表示Q点在第i个叶片局部坐标系ox^by^bz^b中沿叶片长度方向的位置；y表示Q 点在第i个叶片局部坐标系ox^by^bz^b中沿叶片厚度方向的位置；表示第i个叶片局部坐标系 ox^by^bz^b中叶片剪切角；x_d表示圆盘旋转时在整体坐标系OXYZ中的X方向位移；y_d表示圆盘 旋转时在整体坐标系OXYZ中的Y方向位移；z_d表示圆盘旋转时在整体坐标系OXYZ中的Z 方向位移；R_d表示圆盘的半径；u为Q点在第i个叶片局部坐标系ox^by^bz^b中叶片长度方向位 移；v为Q点在第i个叶片局部坐标系ox^by^bz^b中叶片厚度方向位移；w为Q点在第i个叶片 局部坐标系ox^by^bz^b中叶片宽度方向位移；

步骤2-1-4、根据叶片的横截面积、叶片密度、叶片长度和Q点的速度，获得叶片动能；

第i个叶片动能T_blade计算公式如下：

$T_{\mathrm{blade}}=\frac{1}{2}{\int }_0^L{\rho }_b{A}_b{\stackrel{.}{r}}_Q^2\mathrm{dx}=\frac{1}{2}(T_1+T_2)---\left(2\right)$

其中，A_b表示叶片的横截面积；ρ_b表示叶片密度；L表示叶片长度；为第i个叶片上 任意一点Q的速度；T₁表示第i个叶片振动产生的动能项；T₂表示第i个叶片与转子耦合振 动产生的动能项；T₁的表达式如下：

公式(3)中，表示圆盘运动的角速度；I_b表示叶片任意位置处的截面惯性矩；表示v 对时间t的一阶导数；表示u对时间t的一阶导数；表示对时间t的一阶导数；

T₂的表达式如下：

公式(4)中，表示Ψ对时间t的一阶导数；

所述的叶片任意位置处的横截面积A_b和截面惯性矩I_b，计算公式如下：

$\left(\begin{array}{c}{A}_b=A_0(1-{\tau }_b\frac{x}{L})(1-{\tau }_h\frac{x}{L}),{\tau }_b=1-\frac{b_1}{b_0}\\ I_b=I_0(1-{\tau }_b\frac{x}{L}){(1-{\tau }_h\frac{x}{L})}^3,{\tau }_h=1-\frac{h_1}{h_0}\end{array}\right)---\left(9\right)$

其中，b₁表示叶尖处的叶片宽度；h₁表示叶尖处的叶片厚度；b₀表示叶根处的叶片宽度； h₀表示叶根处的叶片厚度；A₀表示叶根处的横截面积，A₀＝b₀h₀；I₀表示叶根处的惯性矩， $I_0=\frac{1}{12}{b}_0{h}_0^3.$

步骤2-2、确定转子-叶片耦合系统中转轴的动能，具体步骤如下：

步骤2-2-1、根据圆盘质心与圆盘形心位移之间的关系，获得圆盘质心在整体坐标系中X 方向的位移x_c和Y方向的位移y_c；

如图5所示，圆盘质心位移与圆盘形心位移之间的关系满足：

$\left(\begin{array}{c}{x}_c=x_d+e\cos (\theta +\Psi )\\ y_c=y_d+e\sin (\theta +\Psi )\end{array}\right)---\left(30\right)$

其中，e表示圆盘质心与形心不对中时的偏心距；x_d、y_d依次表示转子圆盘形心处的水平 和垂直方向上的位移；

步骤2-2-2、根据获得的圆盘质心位移，结合圆盘的极转动惯量、轴段的极转动惯量、圆 盘的直径转动惯量、轴段的直径转动惯量和圆盘在圆盘坐标系中的摆角，确定系统中转轴的 动能；

转轴的动能T_rotor计算公式如下：

$\left(\begin{array}{c}{T}_{\mathrm{rotor}}=\frac{1}{2}{J}_{\mathrm{pd}}{(\stackrel{.}{\theta }+\stackrel{.}{\Psi })}^2+\frac{1}{2}{m}_d({\stackrel{.}{x}}_c^2+{\stackrel{.}{y}}_c^2)-J_{\mathrm{pd}}(\stackrel{.}{\theta }+\stackrel{.}{\Psi }){\stackrel{.}{\theta }}_{\mathrm{yd}}{\theta }_{\mathrm{xd}}+\frac{1}{2}{J}_{\mathrm{dd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}^2+{\stackrel{.}{\theta }}_{\mathrm{yd}}^2)\\ +\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{m}_j({\stackrel{.}{x}}_j^2+{\stackrel{.}{y}}_j^2)-\underset{j=1}{\overset{N_d}{\Sigma }}{J}_{\mathrm{pj}}\stackrel{.}{\theta }{\stackrel{.}{\theta }}_{\mathrm{yj}}{\theta }_{\mathrm{xj}}+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{J}_{\mathrm{dj}}({\stackrel{.}{\theta }}_{\mathrm{xj}}^2+{\stackrel{.}{\theta }}_{\mathrm{yj}}^2)\end{array}\right)---\left(5\right)$

式中，表示圆盘质心在整体坐标系中X方向位移x_c对时间t的一阶导数；表示圆盘 质心在整体坐标系中Y方向位移y_c对时间t的一阶导数；J_pd表示圆盘处集中质量点的极转动 惯量；J_dd表示圆盘处集中质量点的直径转动惯量；J_pj表示第j个集中质量点的极转动惯量； J_dj表示第j个集中质量点的直径转动惯量；θ_xd表示圆盘绕整体坐标系X轴的摆角；θ_yd表示圆 盘绕整体坐标系Y轴的摆角；m_d表示圆盘的质量；m_j表示转轴上的第j个集中质量点的质量； x_j表示转轴上的第j个集中质量点在整体坐标系中X方向位移；y_j表示转轴上的第j个集中质 量点在整体坐标系中Y方向位移；θ_xj表示转轴上的第j个集中质量点绕整体坐标系X轴的转 角；θ_yj表示转轴上的第j个集中质量点绕整体坐标系Y轴的转角；表示θ_xd对时间t的一阶 导数；表示θ_yd对时间t的一阶导数；表示θ_xj对时间t的一阶导数；表示θ_yj对时间t 的一阶导数；N_d表示集中质量点数；j≠p_disc，p_disc表示圆盘处集中质量点的编号；

获得转轴动能最终计算公式如下：

$\left(\begin{array}{c}{T}_{\mathrm{rotor}}=\frac{1}{2}{J}_{\mathrm{pd}}{(\stackrel{.}{\theta }+\stackrel{.}{\Psi })}^2+\frac{1}{2}{m}_d({({\stackrel{.}{x}}_d-e\sin (\theta +\Psi )(\stackrel{.}{\theta }+\stackrel{.}{\Psi }))}^2+{({\stackrel{.}{y}}_d+e\cos (\theta +\Psi )(\stackrel{.}{\theta }+\stackrel{.}{\Psi }))}^2)\\ -(\stackrel{.}{\theta }+\stackrel{.}{\Psi })J_{\mathrm{pd}}{\stackrel{.}{\theta }}_{\mathrm{yd}}{\theta }_{\mathrm{xd}}+\frac{1}{2}{J}_{\mathrm{dd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}^2+{\stackrel{.}{\theta }}_{\mathrm{yd}}^2)+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{m}_j({\stackrel{.}{x}}_j^2+{\stackrel{.}{y}}_j^2)-\underset{j=1}{\overset{N_d}{\Sigma }}\stackrel{.}{\theta }{J}_{\mathrm{pj}}{\stackrel{.}{\theta }}_{\mathrm{yj}}{\theta }_{\mathrm{xj}}+\underset{j=1}{\overset{N_d}{\Sigma }}\frac{1}{2}{J}_{\mathrm{dj}}{({\stackrel{.}{\theta }}_{\mathrm{xj}}+{\stackrel{.}{\theta }}_{\mathrm{yj}})}^2\end{array}\right)---\left(6\right)$

步骤2-3、将转子-叶片耦合系统的叶片动能与转轴动能进行求和，获得系统总的动能；

$T_{\mathrm{total}}=\underset{i=1}{\overset{N_b}{\Sigma }}{T}_{\mathrm{blade}}+T_{\mathrm{rotor}}---\left(31\right)$

步骤3、在转子-叶片耦合系统旋转态下，确定系统的总势能，具体如下：

步骤3-1、根据叶片的弹性模量、叶片的剪切模量、剪切修正系数、离心力和法向力，确 定转子-叶片耦合系统中第i个叶片的势能；

第i个叶片的势能V_blade计算公式如下：

其中，E_b表示叶片的弹性模量；G_b表示叶片的剪切模量；κ表示叶片的剪切修正系数； f_c表示叶片的离心力；F_n表示叶尖受到的法向力；

所述的离心力，其表达式为：

$f_c\left(x\right)={\int }_x^{L}d{f}_c\left(x\right)=\frac{1}{2}{\rho }_b{A}_b{\stackrel{.}{\theta }}^2(L^2+2R_{d}L-2R_{d}x-x^2)---\left(32\right)$

步骤3-2、根据圆盘处转轴的扭转刚度和转轴的各集中质量点位移向量，确定转子-叶片 耦合系统中转轴的势能；

转轴的势能计算公式如下：

$V_{\mathrm{rotor}}=\frac{1}{2}{k}_{\Psi }{\Psi }^2+\frac{1}{2}{q}_r^T{K}_r{q}_r---\left(8\right)$

其中，k_Ψ表示圆盘处转轴的扭转刚度；q_r表示转轴的各集中质量点位移向量， q_r＝[…x_j θ_yj……y_j θ_xj…]^T；K_r表示转轴刚度矩阵；

所述的转轴刚度矩阵，包括结构刚度矩阵和轴承支承刚度矩阵，其表达式为：

其中，K_x表示转轴结构刚度矩阵X方向刚度矩阵；K_y为转轴结构刚度矩阵Y方向刚度 矩阵；K_bx表示轴承在X方向的刚度矩阵；K_by表示轴承在Y方向的刚度矩阵：

$\left(\begin{array}{c}{K}_{\mathrm{bx}}=\mathrm{diag}\left(\begin{array}{ccccc}0...0& k_{\mathrm{bx}1}& 0...0& k_{\mathrm{bx}2}& 0...0\end{array}\right)\\ K_{\mathrm{by}}=\mathrm{diag}\left(\begin{array}{ccccc}0...0& k_{\mathrm{by}1}& 0...0& k_{\mathrm{by}2}& 0...0\end{array}\right)\end{array}\right)---\left(34\right)$

步骤3-3、将转子-叶片耦合系统的叶片势能与转轴势能进行求和，获得系统总的势能；

$V_{\mathrm{total}}=\underset{i=1}{\overset{N_b}{\Sigma }}{V}_{\mathrm{blade}}+V_{\mathrm{rotor}}---\left(35\right)$

步骤4、根据获得的转子-叶片耦合系统的总的动能和总的势能，结合哈密顿原理通过变 分运算确定转子-叶片耦合系统的运动情况，具体如下：

步骤4-1、根据获得的系统总的动能和势能，结合哈密顿原理进行变分运算，确定第i个 叶片的运动微分方程；

所述的确定第i个叶片的运动微分方程，具体如下：

将系统总动能和总势能，代入哈密顿原理表达式(10)中：

$\delta {\int }_{t_1}^{t_2}(T_{\mathrm{total}}-V_{\mathrm{total}})\mathrm{dt}=0---\left(\text{10}\right)$

其中，T_total表示转子-叶片耦合系统总动能；V_total表示转子-叶片耦合系统总的势能；t₁表 示初始时间；t₂表示终止时间；

以δu为独立变量进行变分运算，获得第i个叶片的第一个运动微分方程为：

其中，ü表示u对时间t的二阶导数；表示θ对时间t的二阶导数；A′_b表示A_b对积分 域x的一阶导数；u′表示u对积分域x的一阶导数；u″表示u对积分域x的二阶导数；表 示Ψ对时间t的二阶导数；表示x_d对时间t的二阶导数；表示y_d对时间t的二阶导数；

以δv为独立变量进行变分运算，获得第i个叶片的第二个运动微分方程为：

其中，表示v对时间t的二阶导数；f′_c(x)表示叶片离心力对积分域x的一阶导数；v′表 示v对积分域x的一阶导数；v″表示v对积分域x的二阶导数；表示对积分域x的一阶 导数；F′_n表示叶片的法向力对积分域x的一阶导数；

以为独立变量进行变分运算，获得第i个叶片的第三个运动微分方程为：

其中，表示对时间t的二阶导数；表示对积分域x的二阶导数；I′_b表示I_b对积 分域x的一阶导数。

步骤4-2、根据哈密顿原理进行变分运算，确定圆盘位置处的扭转运动微分方程；

所述的确定圆盘位置处的扭转运动微分方程，具体如下：

以δΨ为独立变量进行变分运算，获得圆盘位置处的扭转运动微分方程为：

$\left(\begin{array}{c}(J_{\mathrm{pd}}+e^2{m}_d+C_1)\stackrel{..}{\Psi }+C_2\stackrel{.}{\Psi }+(k_{\Psi }-e^2{m}_d{\stackrel{.}{\theta }}^2+C_3)\Psi \\ -{\mathrm{em}}_d\sin \theta {\stackrel{..}{x}}_d+{\mathrm{em}}_d\cos \theta {\stackrel{..}{y}}_d+e^2{m}_d\stackrel{..}{\theta }-J_{\mathrm{pd}}({\stackrel{.}{\theta }}_{\mathrm{xd}}{\stackrel{.}{\theta }}_{\mathrm{yd}}+{\theta }_{\mathrm{xd}}{\stackrel{..}{\theta }}_{\mathrm{yd}})\\ +\stackrel{..}{\theta }{J}_{\mathrm{pd}}+C_4-C_5-C_6=0\end{array}\right)---\left(14\right)$

其中，

步骤4-3、根据哈密顿原理进行变分运算，确定转轴的横向运动微分方程；

所述的确定转轴的横向运动微分方程，具体如下：

以δq_r为独立变量进行变分运算，获得转轴的横向运动微分方程为：

$M_r{\stackrel{..}{q}}_r+(G_r+C_b){\stackrel{.}{q}}_r+K_r{q}_r=F_r---\left(15\right)$

其中，M_r表示转轴的质量矩阵；G_r表示转轴的陀螺矩阵；F_r表示转轴的外激振力向量； C_b为转轴轴承的阻尼矩阵。

所述转轴的质量矩阵，其表达式为：

其中，M_rx表示转轴质量矩阵X方向质量矩阵；M_ry表示转轴质量矩阵Y方向质量矩阵； $M_{\mathrm{rx}}=M_{\mathrm{ry}}=\mathrm{diag}{\left(\begin{array}{ccccccccc}{m}_1& J_{d1}& ...& m_d+N_b{\int }_0^L{\rho }_b{A}_b\mathrm{dx}& J_{\mathrm{dd}}& ...& m_j& J_{\mathrm{dj}}& ...\end{array}\right)}^T.$

所述转轴的陀螺矩阵，其表达式为：

其中，J₁＝diag[0 J_p1 … 0 J_pd … 0 J_pj …]^T。

所述转轴轴承的阻尼矩阵，其表达式为：

其中，C_bx表示转轴轴承的阻尼矩阵X方向阻尼矩阵；C_by表示转轴轴承的阻尼矩阵Y方 向阻尼矩阵；C_bx＝diag[0… 0c_bx10 … 0c_bx20 …0]；C_by＝diag[0… 0c_by10 … 0c_by20 …0]。

步骤4-4、根据正则坐标对第i个叶片的运动微分方程中的叶片的长度方向位移、厚度方 向位移和剪切角进行离散化处理，获得第i个叶片的质量矩阵、阻尼矩阵和刚度矩阵；根据 正则坐标对圆盘位置处的扭转运动微分方程中的叶片的长度方向位移、厚度方向位移和剪切 角进行离散化处理，获得系统质量矩阵耦合项、阻尼矩阵耦合项和刚度矩阵耦合项；

本发明实施例中采用假设模态法对叶片的每个方向变形进行离散化，通过引入正则坐标 U_ξ(t)，V_ξ(t)和ψ_ξ(t)，纵向位移u(x，t)、横向位移v(x，t)和剪切角可以写作：

    $\left(\begin{array}{c}{\phi }_{1\xi }\left(x\right)=\frac{\sin \left({\alpha }_{\xi }x\right)}{{\alpha }_{\xi }}\\ {\phi }_{2\xi }\left(x\right)=\frac{1-\cos \left({\alpha }_{\xi }x\right)}{{\alpha }_{\xi }}\\ {\phi }_{3\xi }\left(x\right)=\sin \left({\alpha }_{\xi }x\right)\end{array}\right)---\left(17\right)$

公式(16)中，U_ξ(t)表示离散处理后叶片的长度方向位移；V_ξ(t)表示离散处理后叶片的 厚度方向位移；ψ_ξ(t)表示离散处理后叶片的剪切角；φ_1ξ(x)表示叶片长度方向振动的第ξ阶振 型函数；φ_2ξ(x)表示叶片厚度方向振动的第ξ阶振型函数；φ_3ξ(x)表示叶片剪切角的第ξ阶振 型函数；公式(17)中，ξ＝1，2，3，…，N_mod，L为叶片的长度，N_mod为截断的模态 数。

步骤4-5、对叶片的质量矩阵、叶片科氏力矩阵、叶片刚度矩阵、质量矩阵耦合项、阻尼 矩阵耦合项、刚度矩阵耦合项、圆盘位置处转轴的扭转质量、圆盘位置处转轴的扭转阻尼、 圆盘位置处转轴的扭转刚度、转轴的质量矩阵、转轴的陀螺矩阵和转轴的刚度矩阵进行组集， 获得转子-叶片耦合系统的质量矩阵、阻尼矩阵和刚度矩阵，结合转子-叶片耦合系统的广义 坐标向量和外激振力向量构建转子-叶片耦合系统的运动微分方程；

根据上面得到的第i个叶片的运动方程以及转子弯扭耦合运动方程，经过离散化，写成 矩阵形式。由于第i个叶片都有相应的质量、阻尼和刚度矩阵，需要将各个叶片的矩阵与转 子系统矩阵组集起来，形成转子-多叶片耦合系统整体矩阵。图6中图(a)、图(b)和图(c) 为转子-多叶片耦合系统矩阵组集示意图。图中N_dof为第i个叶片的自由度数，N_rdof为转轴自 由度数。

所述的构建转子-叶片耦合系统的运动微分方程，具体如下：

转子-叶片耦合系统的运动微分方程如下：

$M_{\mathrm{RB}}{\stackrel{..}{q}}_{\mathrm{RB}}+D_{\mathrm{RB}}{\stackrel{.}{q}}_{\mathrm{RB}}+K_{\mathrm{RB}}{q}_{\mathrm{RB}}=f_{\mathrm{RB}}---\left(18\right)$

公式(18)中，M_RB表示转子-叶片耦合系统的质量矩阵，

$M_{\mathrm{RB}}=\left(\begin{array}{ccc}{M}_b& M_{c1}& M_{c2}\\ M_{c1}^T& M_{\Psi }& 0\\ M_{c2}^T& 0& M_r\end{array}\right)---\left(19\right)$

公式(19)中，M_b表示叶片的质量矩阵；M_c1为扭转运动的质量耦合项；M_c2为横 向运动的质量耦合项；M_Ψ为圆盘位置处转轴的扭转质量；M_r转轴的质量矩阵；

所述叶片的质量矩阵，其表达式为：

$M_b=\mathrm{diag}\left(\begin{array}{ccccc}{M}_b^1& ...& M_b^i& ...& M_b^{N_b}\end{array}\right)---\left(39\right)$

其中，为第i个叶片的质量矩阵，i＝1，2，3，…，N_b，矩阵各元素如下：

$M_b^i(m,n)=(1+{\Psi }^2){\rho }_b{\int }_0^L{A}_b{\phi }_{1m}{\phi }_{1n}\mathrm{dx},M_b^i(m,n)$表示矩阵的第m行，n列元素，后文中 矩阵元素表示方式同理，m，n＝ξ＝1，2，…，N_mod；

$M_b^i(m+N_{\mathrm{mod}},n+N_{\mathrm{mod}})=(1+{\Psi }^2){\rho }_b{\int }_0^L{A}_b{\phi }_{2m}{\phi }_{2n}\mathrm{dx},$

$M_b^i(m+{2N}_{\mathrm{mod}},n+{2N}_{\mathrm{mod}})=(1+{\Psi }^2){\rho }_b{\int }_0^L{I}_b{\phi }_{3m}{\phi }_{3n}\mathrm{dx};$

所述扭转运动的质量耦合项，其表达式为：

$M_{c1}={\left(\begin{array}{ccccc}{M}_{c1}^1& ...& M_{c1}^i& ...& M_{c1}^{N_b}\end{array}\right)}^T---\left(40\right)$

所述横向运动的质量耦合项，其表达式为：

$M_{c2}={\left(\begin{array}{ccccc}{M}_{c2}^1& ...& M_{c2}^i& ...& M_{c2}^{N_b}\end{array}\right)}^T---\left(41\right)$

其中，表示第i个叶片的扭转运动的质量耦合项，表示第i个叶片的横向运 动的质量耦合项，矩阵各元素如下：$M_{c1}^i(m,1)=0,M_{c1}^i(m+N_{\mathrm{mod}},1)={\rho }_b{\int }_0^L(R_d+x)A_b{\phi }_{2m}\mathrm{dx},$

$M_{c2}^i(m+2N_{\mathrm{mod}},2p_{\mathrm{disc}}-1+\frac{N_{\mathrm{rdof}}}{2})=0;$

所述圆盘位置处转轴的扭转质量，其表达式为：

其中，$U={[U_1,...,U_{N_{\mathrm{mod}}}]}^T,$$V={[V_1,...,V_{N_{\mathrm{mod}}}]}^T,\psi ={[{\psi }_1,...,{\psi }_{N_{\mathrm{mod}}}]}^T.$

D_RB表示转子-叶片耦合系统的阻尼矩阵；D_RB的表达式如下：

D_RB＝G_RB+D_b     (20)

公式(20)中，D_b为系统轴承的阻尼矩阵；

G_RB的表达式如下：

$G_{\mathrm{RB}}=\left(\begin{array}{ccc}{G}_b& G_{c1}& 0\\ -G_{c1}^T& G_{\Psi }& 0\\ 0& 0& G_r\end{array}\right)---\left(21\right)$

公式(21)中，G_b表示叶片的科氏力矩阵；G_c1表示阻尼耦合项；G_Ψ表示圆盘位置处转 轴的扭转阻尼；G_r表示转轴的陀螺矩阵；

所述叶片的科氏力矩阵，其表达式为：

$G_b=\mathrm{diag}\left(\begin{array}{ccccc}{G}_b^1& ...& G_b^i& ...& G_b^{N_b}\end{array}\right)---\left(43\right)$

其中，为第i个叶片的科氏力矩阵，矩阵各元素如下：$G_b^i(m,n+N_{\mathrm{mod}})=-2((1+{\Psi }^2)\stackrel{.}{\theta }+\stackrel{.}{\Psi }){\rho }_b{\int }_0^L{A}_b{\phi }_{2n}{\phi }_{1m}\mathrm{dx},G_b^i(m+N_{\mathrm{mod}},n)=2((1+{\Psi }^2)\stackrel{.}{\theta }+\stackrel{.}{\Psi }){\rho }_b{\int }_0^L{A}_b{\phi }_{1n}{\phi }_{2m}\mathrm{dx},$$G_b^i(m+N_{\mathrm{mod}},n+N_{\mathrm{mod}})=2\Psi \stackrel{.}{\Psi }\stackrel{.}{\theta }{\rho }_b{\int }_0^L{A}_b{\phi }_{2n}{\phi }_{2m}\mathrm{dx},G_b^i(m+2N_{\mathrm{mod}},n+2N_{\mathrm{mod}})=2\Psi \stackrel{.}{\Psi }\stackrel{.}{\theta }{\rho }_b{\int }_0^L{I}_b{\phi }_{3n}{\phi }_{3m}\mathrm{dx};$

所述阻尼耦合项矩阵，其表达式为：

$G_{c1}={\left(\begin{array}{ccccc}{G}_{c1}^1& ...& G_{c1}^i& ...& G_{c1}^{N_b}\end{array}\right)}^T---\left(44\right)$

其中，为第i个叶片的阻尼耦合项矩阵，矩阵各元素如下：

$G_{c1}^i(m,1)=-2\stackrel{.}{\theta }{\rho }_b{\int }_0^L{A}_b(R_d+x){\phi }_{1m}\mathrm{dx},G_{c1}^i(m+N_{\mathrm{mod}},1)=G_{c1}^i(m+2N_{\mathrm{mod}},1)=0;$

所述圆盘位置处转轴的扭转阻尼，其表达式为：

其中，表示U对时间t的一阶导，表示V对时间t的一阶导。

K_RB表示转子-叶片耦合系统的刚度矩阵；

$K_{\mathrm{RB}}=\left(\begin{array}{ccc}{K}_b& 0& 0\\ 0& K_{\Psi }& 0\\ 0& 0& K_r\end{array}\right)---\left(22\right)$

公式(22)中，K_b表示叶片刚度矩阵；K_Ψ表示圆盘位置处转轴的扭转刚度；K_r表示转轴 的刚度矩阵；

所述叶片刚度矩阵，其表达式为：

$K_b=\mathrm{diag}\left(\begin{array}{ccccc}{K}_b^1& ...& K_b^i& ...& K_b^{N_b}\end{array}\right),---\left(46\right)$

其中，为第i个叶片的刚度矩阵，矩阵各元素如下：

$K_b^i(m,n)=(\Psi \stackrel{..}{\Psi }-(1+{\Psi }^2){\stackrel{.}{\theta }}^2-2\stackrel{.}{\Psi }\stackrel{.}{\theta }){\rho }_b{\int }_0^L{A}_b{\phi }_{1n}{\phi }_{1m}\mathrm{dx}+E{A}_b{\phi }_{1n}^{\prime }{{\phi }_{1m}|}_{x=L}-{\int }_0^{L}E(A_b^{\prime }{\phi }_{1n}^{\prime }+A_b{\phi }_{1n}^{\prime \prime }){\phi }_{1m}\mathrm{dx},$

$K_b^i(m,n+N_{\mathrm{mod}})=-((1+{\Psi }^2)\stackrel{..}{\theta }+2\Psi \stackrel{.}{\Psi }\stackrel{.}{\theta }+\stackrel{..}{\Psi }){\rho }_b{\int }_0^L{A}_b{\phi }_{2n}{\phi }_{1m}\mathrm{dx},$

$K_b^i(m+N_{\mathrm{mod}},n)=(\stackrel{..}{\Psi }+(1+{\Psi }^2)\stackrel{..}{\theta }+2\Psi \stackrel{.}{\Psi }\stackrel{.}{\theta })\rho {\int }_0^L{A}_b{\phi }_{1n}{\phi }_{2m}\mathrm{dx},$

$\left(\begin{array}{c}{K}_b^i(m+N_{\mathrm{mod}},n+N_{\mathrm{mod}})=\kappa G_b{A}_b{\phi }_{2n}^{\prime }{{\phi }_{2m}|}_{x=L}-{\int }_0^L\kappa G_b(A_b^{\prime }{\phi }_{2n}^{\prime }+A_b{\phi }_{2n}^{\prime \prime }){\phi }_{2m}\mathrm{dx}\\ -{\int }_0^L(f_c^{\prime }\left(x\right){\phi }_{2n}^{\prime }+f_c\left(x\right){\phi }_{2n}^{\prime \prime }){\phi }_{2m}\mathrm{dx}\\ +F_n{\phi }_{2n}^{\prime }{{\phi }_{2m}|}_{x=L}-{\int }_0^L(F_n{\phi }_{2n}^{\prime \prime }+F_n^{\prime }{\phi }_{2n}^{\prime }){\phi }_{2m}\mathrm{dx}\\ +(\Psi \stackrel{..}{\Psi }+{\stackrel{.}{\Psi }}^2-{\Psi }^2{\stackrel{.}{\theta }}^2-{(\stackrel{.}{\Psi }+\stackrel{.}{\theta })}^2){\rho }_b{\int }_0^L{A}_b{\phi }_{2n}{\phi }_{2m}\mathrm{dx}\end{array}\right),$

$K_b^i(m+N_{\mathrm{mod}},n+2N_{\mathrm{mod}})=-\kappa G_b{A}_b{\phi }_{3n}{{\phi }_{2m}|}_{x=L}+{\int }_0^L\kappa G_b(A_b^{\prime }{\phi }_{3n}+A_b{\phi }_{3n}^{\prime }){\phi }_{2m}\mathrm{dx},$

$K_b^i(m+3N_{\mathrm{mod}},n+2N_{\mathrm{mod}})=-{\int }_0^L\kappa G_b{A}_b{\phi }_{2n}^{\prime }{\phi }_{3m}\mathrm{dx},$

$\left(\begin{array}{c}{K}_b^i(m+3N_{\mathrm{mod}},n+3N_{\mathrm{mod}})=E{I}_b{\phi }_{3n}^{\prime }{{\phi }_{3m}|}_{x=L}-{\int }_0^{L}E(I_b^{\prime }{\phi }_{3n}^{\prime }+I{\phi }_{3n}^{\prime \prime }){\phi }_{3m}\mathrm{dx}+{\int }_0^f\kappa G_b{A}_b{\phi }_{3n}{\phi }_{3m}\mathrm{dx}\\ +(\Psi \stackrel{..}{\Psi }-{\Psi }^2{\stackrel{.}{\theta }}^2-2\stackrel{.}{\Psi }\stackrel{.}{\theta }-{\stackrel{.}{\theta }}^2)\rho {\int }_0^L{I}_b{\phi }_{3n}{\phi }_{3m}\mathrm{dx}\end{array};\right)$其中，φ′_1n表示φ_1n对 积分域x的一阶导，φ″_1n表示φ_1n对积分域x的二阶导，φ′_2n表示φ_2n对积分域x的一阶导，φ″_2n表 示φ_2n对积分域x的二阶导，φ′_3n表示φ_3n对积分域x的一阶导，φ″_3n表示φ_3n对积分域x的二阶导；

所述圆盘位置处转轴的扭转刚度，其表达式为：

其中，表示U对时间t的二阶导，表示V对时间t的二阶导。

q_RB表示转子-叶片耦合系统的广义坐标向量，q_b表示叶片的广义坐标 向量；表示q_RB对时间t的二阶导数；表示q_RB对时间t的一阶导数；

所述叶片的广义坐标向量，其表达式为：

$q_b={\left(\begin{array}{ccccc}{q}_b^1& ...& q_b^i& ...& q_b^{N_b}\end{array}\right)}^T---\left(48\right)$

其中，为第i个叶片的广义坐标向量，

f_RB表示转子-叶片耦合系统的外激振力向量。

步骤5、设定外激振力向量为零，根据获得的转子-叶片耦合系统的质量矩阵、阻尼矩阵 和刚度矩阵，确定转子-叶片耦合系统的固有频率；

步骤5-1、设置外激振力向量f_RB为零，构建特征方程；

特征方程如下：

$\left(\begin{array}{cc}{D_{\mathrm{RB}}}^T& {M_{\mathrm{RB}}}^T\\ {M_{\mathrm{RB}}}^T& 0\end{array}\right)\left(\begin{array}{c}{\stackrel{.}{q}}_{\mathrm{RB}}\\ {\stackrel{..}{q}}_{\mathrm{RB}}\end{array}\right)+\left(\begin{array}{cc}{K_{\mathrm{RB}}}^T& 0\\ 0& -{M_{\mathrm{RB}}}^T\end{array}\right)\left(\begin{array}{c}{q}_{\mathrm{RB}}\\ {\stackrel{.}{q}}_{\mathrm{RB}}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)---\left(23\right)$

步骤5-2、对所构建的特征方程进行简化；

即设置如下关系：

$A_{\mathrm{mod}}=\left(\begin{array}{cc}{K_{\mathrm{RB}}}^T& 0\\ 0& -{M_{\mathrm{RB}}}^T\end{array}\right),B_{\mathrm{mod}}=\left(\begin{array}{cc}{D_{\mathrm{RB}}}^T& {M_{\mathrm{RB}}}^T\\ {M_{\mathrm{RB}}}^T& 0\end{array}\right),Z_{\mathrm{mod}}=\left(\begin{array}{c}{q}_{\mathrm{RB}}\\ {\stackrel{.}{q}}_{\mathrm{RB}}\end{array}\right)---\left(24\right)$

将A_mod、B_mod、Z_mod代入公式(23)，获得如下方程：

$B_{\mathrm{mod}}{\stackrel{.}{Z}}_{\mathrm{mod}}+A_{\mathrm{mod}}{Z}_{\mathrm{mod}}=0---\left(25\right)$

并设置Z_mod＝Ze^-λt，代入到式(25)中化简获得如下方程：

(A_mod-λB_mod)Z＝0     (26)

公式(26)中，λ表示特征方程的特征值；Z表示特征方程的特征向量；

步骤5-3、将化简后方程中的系数行列式设置为零；

即：

|A_mod-λB_mod|＝0     (27)

步骤5-4、根据该系数行列式的特征值，确定转子-叶片耦合系统的固有频率；

即计算获得系数行列式的一组特征值λ，取其虚部的绝对值除以2π，并进行从小到大排 序，获得一组固有频率ω_k，其中，k表示转子-叶片耦合系统模态的第k阶，k＝1，2，…。

步骤6、根据获得的转子-叶片耦合系统的固有频率，确定系统的工作转速，避免共振， 使系统运行稳定；此外，还可根据公式(18)进行转子-叶片耦合系统不平衡响应、叶尖碰摩 等故障分析，从而指导系统结构的设计和优化。

模型验证和数值仿真分析

为了更好的验证本发明所提出的方法，如图7所示建立了两个有限元模型来分析本发明 的有效性。图7中图(a)为限元模型1，其中转轴划分为14个Timoshenko梁单元(ANSYS 中BEAM188单元)和15个节点，除了最左端节点约束其轴向自由度外，其余节点均采用6 个自由度；盘采用壳单元来模拟(ANSYS中SHELL181单元)，划分了96个单元和335个节 点。每个节点有6个自由度；四个叶片也采用Timoshenko梁单元(ANSYS中BEAM188单 元)，每个叶片划分4个单元和5个节点。轴和盘采用刚性连接，盘和叶片通过共用节点实现 连接。轴承采用弹簧-阻尼单元(ANSYS中COMBI214单元)；图7中图(b)为限元模型2， 其中盘采用ANSYS中MASS21单元来模拟，轴上的节点9和叶片根部通过命令“CERIG” 实现刚性连接，其它处理过程同有限元模型1。

本发明实施例中，对两个子系统即转子-轴承系统和叶片的固有频率进行分析，如表2所 示；对于转子-轴承系统，叶片的影响被忽略，其它部件的处理方式同转子-叶片耦合系统。 对于叶片，叶片根部固支，仅仅第1阶弯曲固有频率被计算。表3为转子-叶片耦合系统固有 频率，展示了4种方法(本发明、有限元模型1、有限元模型2和实验)的固有频率对比结 果，本发明和实验结果的误差对比，如表3第7列所示。

表2

表3

注：“----“代表无值

由表2和表3对比，获得以下结论：

(1)由于叶片附加质量和惯性的影响，转子-叶片耦合系统的大部分固有频率略小于转子- 轴承系统。

(2)本发明与有限元模型和测定的试验结果均吻合很好，这均表明本发明所提的动力学 模型是合理的。可能存在少许误差的原因：叶片和盘实际连接形式为榫连接触，然而仿真忽 略这种连接影响；测试过程中轻质加速传感器放在叶尖位置会有一定的影响对于测试结果； 转轴和叶片的几何参数和材料参数也影响仿真结果。

(3)由于均为刚性盘的假定，本发明更接近于有限元模型2的结果，而有限元模型1的 结果更接近于实测结果，这是因为有限元模型2中将圆盘考虑为弹性体，而本发明和有限元 模型1将圆盘考虑为刚性体。表2与表3中本发明与有限元模型和实验结果的误差均不超过 5％，说明本发明满足计算精度，而且由于本发明采用较少的自由度，因而具有较高的计算效 率，在对于计算转子-叶片耦合系统复杂的非线性振动响应方面更具前景。

对于转子-叶片系统旋转态下的模态进行实际测试是十分困难的，而在ANSYS软件中， 对于非轴对称的转子-叶片系统而言，在静态坐标系中叶片的离心刚化效应和转轴的陀螺效应 不能同时考虑。分别考虑转轴(轴和盘)的陀螺效应和叶片的旋转软化、离心刚化效应，计 算转子-叶片系统的动频。其中计算转子系统的扭转固有频率和叶片的弯曲固有频率时仅考虑 叶片的离心刚化和旋转软化的影响，计算转子有关的俯仰、摆动和弯曲固有频率时仅考虑了 轴和刚性盘陀螺效应的影响。为了方便对比，这些动频均绘制在一个Campbell图中。在4、6 和8叶片情况下，三种仿真模型计算的Campbell图，如图8中图(a)、图(b)和图(c)所 示。图中A表示转轴扭转模态，B₁表示系统水平俯仰模态，B₂表示系统竖直俯仰模态，C₁表示圆盘反进动摆动模态，C₂表示圆盘正进动摆动模态，D叶片一阶弯曲模态，E₁表示转轴 反进动弯曲模态，E₂表示转轴正进动弯曲模态；“-“代表本发明，“----“代表有限元模型1， “-·-“代表有限元模型2。图8所展示的动力学现象如下：

(1)动频曲线的对比表明，本发明所提的方法对于所研究的转子-叶片系统是非常准确的。

(2)扭转固有频率(A)随着转速的增加而减小；扭转固有频率的缩减程度随着叶片数的 增加而增大；其中缩减程度最大的为本发明模型，其次为有限元模型2，最后为有限元模型1， 这也表明盘的柔性对于系统的扭转固有频率有一定影响。由于陀螺力矩的影响，系统俯仰振 动模态(B₁和B₂)有轻微的改变。

(3)随着转速的增加，三种类型的振动模态改变较大。由于陀螺效应的影响，转速对摆 动模态(C₁和C₂)有最大的影响；由于叶片离心刚化和旋转软化的影响，转速对叶片第1阶 弯曲模态(D)有第二大影响；第三大影响为转子的弯曲模态(E₁和E₂)。

(4)随着叶片数的增加，耦合振动模态数(D附近)也增加。在4、6、8叶片情况下，分 别出现了和叶片数相同的耦合振动模态(参见表3所列4叶片情况)。由于存在较小的误差， 这些动频曲线机会重叠在一起(参见图8中D)。

本发明实施例中旋转柔性叶片采用Timoshenko梁来模拟，其考虑旋转导致的离心刚化、 旋转软化和科氏力影响；刚性盘采用集中质量点来模拟，考虑陀螺效应；转轴通过多个无质 量弹性轴段相连的多个质量点来模拟，考虑转轴的横向运动和陀螺力矩影响。本发明实施例 中转轴的扭转振动仅在盘的位置被考虑。转子-叶片系统的运动方程通过Hamilton原理结合假 设模态法来进行离散化。为了计算系统的无阻尼固有频率，忽略了轴承的阻尼矩阵，和非线 性方程中转轴扭转、叶片变形的高阶耦合项。

通过两个有限元模型计算和实验测定的固有频率对本发明进行验证。结果表明本发明计 算的固有频率和有限元模型结果及实测结果均吻合很好；系统扭转固有频率随着转速的增加 而减小，其减幅程度随着叶片数的增加而增大；在4、6和8叶片情况下会出现和叶片数相同的 多叶片耦合模态。