CN103984813B - 一种离心压缩机裂纹叶轮结构的振动建模与分析方法 - Google Patents

Abstract

本发明公开了一种离心压缩机裂纹叶轮结构的振动建模与分析方法，首先，提出了一种离心叶轮的建模方法，通过叶轮的一个扇区的有限元模型旋转变换得到整个模型，采用混合界面模态综合法对系统的自由度进行缩减，并通过在裂纹界面上定义接触的方式模拟裂纹的呼吸效应。其次，提出了一种求解大型对称稀疏矩阵逆矩阵的方法，显著地减少了计算所需的内存空间。最后，提出了一种对裂纹叶轮结构的谐振频率进行统计分析的方法，考虑了实际叶轮存在的制造误差及状态劣化等失谐因素，获得了在随机失谐模式下结构谐振频率的统计规律性。该方法在保证计算精度的同时，计算效率提升效果显著，为叶轮的优化设计以及裂纹的定量诊断提供了一种高效的分析方法。

Description

一种离心压缩机裂纹叶轮结构的振动建模与分析方法

技术领域

本发明属结构动力学建模领域，具体涉及一种离心压缩机裂纹叶轮结构的振动建模与分析方法。

背景技术

离心压缩机是一种广泛应用于石化、钢铁、电力等工业领域的重要设备。叶轮是离心压缩机转子的核心部件，也是最容易发生故障的部件之一。由于作用于叶轮上的复杂的交变应力的影响，结构容易在薄弱位置处萌发疲劳裂纹。当裂纹扩展到一定深度时，就会发生叶片断裂事故，对处于高速运转状态下的转子，断裂的叶片不仅会破坏机组，还威胁着现场人员的生命安全。

离心压缩机转子是一个复杂的系统，而且运行在复杂气流激励的环境下。为开发一项有效的叶轮裂纹诊断技术，需要深入分析裂纹对结构动态响应的影响规律。叶轮结构的有限元模型通常包含数十万甚至数百万自由度，分析如此大的模型对计算机硬件的要求非常高。另外，裂纹的呼吸效应对结构动态响应有显著的影响，是分析中不能忽略的因素。而裂纹的呼吸效应是一种非线性接触问题，求解过程包含很多个迭代步，直接采用结构的有限元模型来模拟裂纹的呼吸效应是非常困难的，不仅计算时间非常长，而且分析的收敛性也难以保证。因此，分析裂纹对系统动态响应的影响规律需要一种高效的、能够定量预测结构响应的建模方法，同时还应能够很容易融合裂纹的呼吸效应。

其次，由于制造误差以及结构状态劣化等因素的影响，实际的叶轮不可避免地出现失谐现象。由于实际叶轮的失谐模式往往是难以测量的，因此，通常引入一些简化的失谐模式来模拟失谐效应。由于叶轮是一种强耦合结构，导致结构的响应幅值对裂纹不敏感，因而传统的针对幅值的统计分析方法无法提供有价值的诊断指标。而叶轮结构的某些谐振频率表现出对裂纹和失谐比较敏感的特性，由于失谐模式具有很强的随机性，分析参数必须建立在统计分析的基础上。然而，由于统计分析往往包含数千甚至上万种随机模式，如果直接采用结构的有限元模型，分析所需的时间是难以接受的。另外，统计分析还需要很高的计算硬件水平，这是普通工作站无法满足的。因此，需要研究在有限的计算资源以及有限的计算时间下实现对响应参数的统计分析的方法。

发明内容

针对上述缺陷或不足，本发明的目的在于提供一种离心压缩机裂纹叶轮结构的振动建模与分析方法，该方法能够融合裂纹的呼吸效应以及随机失谐效应，实现在有限计算资源条件下对裂纹叶轮结构动态响应以及谐振频率统计规律性的分析。

为达到以上目的，本发明的技术方案为：

包括以步骤：

1)首先，通过叶轮一个扇区的有限元模型旋转变换，得到整个结构的有限元模型；然后，运用混合界面模态综合法对有限元模型的系统自由度进行缩减，并在系统自由度进行缩减过程中，采用能够减少计算所需空间的大型稀疏逆矩阵求解方法进行计算，最后，在有限元模型的裂纹界面上定义接触对，以模拟裂纹的呼吸效应；

2)融合实际叶轮存在的制造误差及状态劣化失谐因素，获得裂纹叶轮结构在随机失谐模式下结构谐振频率的统计规律性，建立谐振频率与裂纹深度和失谐水平间的定量关系。

所述步骤1)中通过叶轮一个扇区的有限元模型旋转变换，得到整个结构的有限元模型方法如下：

首先，获得叶轮一个扇区有限元模型的质量和刚度矩阵分别为M_s和K_s。在笛卡尔坐标系下，第n扇区的质量矩阵M_n和刚度矩阵K_n通过旋转变换的方式得到：

式中，R_n为从基本扇区到第n扇区的旋转变换矩阵，且R_n为：

式中，I_S为维数为n_s的一个单位矩阵，且n_s为一个扇区中所包含的节点数，为Kronecker乘积,R_ne为节点的旋转变换矩阵；第n扇区与基本扇区对应节点间的旋转角度θ＝(n-1)2π/N，其中N为叶轮的扇区总数；将各扇区的质量矩阵和刚度矩阵的计算合并起来，则：

式中，M_R、K_R为旋转变换后的结构质量矩阵和刚度矩阵，R为一个以R_n为对角元素的pseudo-block对角矩阵；

其次，采用混合界面模态综合法对系统的自由度进行缩减，该混合界面模态综合方法对应的自由度变换方程为：

式中，u_i和u_b分别为系统的内部和界面自由度，p_k、p_a和u_b分别为保留的正则模态自由度、附加模态自由度和固定界面边界自由度，Φ_k、Φ_a和Φ_c分别为保留的正则模态集、剩余模态和约束模态；

然后，将裂纹界面上的节点自由度保留在系统的缩减模型中，并在界面节点之间定义接触，以模拟裂纹界面的呼吸效应；接触对的非线性接触力是在节点的局部坐标系下定义的，即：

式中，和分别为目标面和接触面上的接触力，k为法向刚度惩罚系数，符号<·>为Macaulay括号，和分别为在局部坐标系下目标节点和接触节点的法向位移，g_i为接触对的法向初始间隙；则包含裂纹呼吸效应的系统非线性动力学方程为：

式中，M、C和K分别为缩减模型的质量、阻尼和刚度矩阵，q为缩减系统的广义自由度向量，b为外载荷向量，f_nl(q)为作用在裂纹界面上的非线性接触力。

所述采用能够减少计算所需空间的大型对称稀疏矩阵逆矩阵的求解方法，该求解方法包括：

在计算大型对称稀疏矩阵K_I的逆矩阵时，将矩阵K_I等分为2×2的子块，即：

则，其逆矩阵可以表示为：

式中，

K_A＝(K₁₁)^-1

K_B＝(K₁₁)^-1K₁₂

其中，子矩阵K₁₁、K₁₂及K₂₂分为2×2的子块后进行计算。

所述步骤2)的具体步骤包括：

首先，计算裂纹叶轮结构的谐振频率，将非线性系统近似为双线性系统；计算中需要的信息为两个状态下的固有频率，状态一为不考虑裂纹界面的呼吸效应情况，状态二为约束了裂纹界面间的法向相对运动，释放了界面间的切向相对运动，则裂纹叶轮结构的第i阶谐振频率f_bi可估计为：

式中，f_1i为系统在状态一的固有频率，f_2i为系统在状态二的固有频率。

其次，引入简化的失谐模式来模拟失谐效应，随机改变各扇区叶片的弹性模量，则第n扇区叶片的弹性模量E_n为：

E_n＝(1+δ_n)E₀

式中，E₀为叶片材料的名义弹性模量，δ_n为服从均匀分布的随机变量，且均值为0，标准差为σ，将该标准差定义为结构的失谐水平；

最后，采用统计分析的方法分析随机失谐对谐振频率的影响规律，所述采用统计分析的方法分析随机失谐对谐振频率的影响规律具体包括：

a、先对叶轮结构的失谐水平进行估计，并随机生成较多数量的随机失谐模式；

b、采用Monte Carlo模拟方法对谐振频率进行统计分析，计算在每一种失谐模式下的结构谐振频率，根据每一种失谐模式下的结果得到各谐振频率的概率分布曲线；

c、分析在均匀分布的失谐模式下的谐振频率的统计分布，并研究分布的均值、方差与失谐水平间的关系；

d、根据谐振频率的统计规律性，建立谐振频率与裂纹深度和失谐水平间的定量关系。

与现有技术比较，本发明的有益效果为：

1)本发明提出的离心叶轮振动建模方法只需要叶轮一个扇区的有限元模型，并采用混合界面模态综合法对系统的自由度进行了显著地缩减，提高了建模的效率，该方法还可以考虑裂纹界面的呼吸效应。

2)通过采用统计分析的方法，考虑了实际叶轮存在的失谐因素，在估计出叶轮的失谐水平之后，即可进而建立起谐振频率与裂纹深度和失谐水平间的定量关系。

3)该方法适用于分析具有失谐和裂纹损伤的叶轮的振动问题，为叶轮的优化设计以及裂纹的定量诊断提供一个高效的分析方法。

进一步的，本发明在求解大型对称稀疏矩阵的逆矩阵时，通过将矩阵分解为若干个子矩阵的形式，使得只有在每一步中参与计算的子矩阵会被导入到内存中，显著地降低了计算所需的内存空间。

附图说明

图1为本发明离心压缩机裂纹叶轮结构振动建模与分析的计算框架示意图；

图2为本发明一个离心压缩机裂纹叶轮结构的有限元模型图；

图3为本发明分析中采用的裂纹呼吸模型图；

图4为本发明谐振频率统计分析的流程图；

图5为本发明裂纹结构的频率响应曲线图；

图6为本发明裂纹叶轮谐振频率随裂纹深度的变化趋势图；

图7为本发明结构第2阶谐振频率的概率密度函数图；

图8为本发明结构第5阶谐振频率的概率密度函数图。

具体实施方式

下面结合附图对本发明做详细描述。

本发明提出了一种离心压缩机裂纹叶轮结构的振动建模与分析方法，通过对系统自由度的缩减，并融合裂纹的呼吸效应和随机失谐效应，可以高效地实现对裂纹叶轮动态响应的分析以及对谐振频率的统计分析，包括以下步骤：

1)通过叶轮一个扇区的有限元模型旋转变换，得到整个结构的有限元模型，运用混合界面模态综合法对有限元模型的系统自由度进行缩减，并在有限元模型的裂纹界面上定义接触对，以模拟裂纹的呼吸效应；在系统自由度进行缩减过程中，采用能够减少计算所需空间的大型稀疏逆矩阵求解方法进行计算；

2)根据有限元模型的裂纹界面上定义接触对以及实际叶轮存在的制造误差及状态劣化产生的失谐因素，获得在随机失谐模式下结构谐振频率的统计规律性，建立谐振频率与裂纹深度和失谐水平间的定量关系。

具体的，采通过将矩阵分解为若干个子矩阵的形式，对有限元模型的系统自由度进行缩减，实现了对大型对称稀疏矩阵逆矩阵的求解，显著降低了计算所需内存空间的要求。

3)运用混合界面模态综合法对有限元模型的系统自由度进行缩减，采用对裂纹叶轮结构的谐振频率进行统计分析，以获得在外部随机失谐因素影响下结构谐振频率的统计规律。

下面结合附图对本发明的内容作进一步详细说明：

本发明所提出的方法的分析流程如图1所示，主要分析步骤详述如下。

(1)扇区模型的旋转变换

首先，基于如图2所示的叶轮结构的有限元模型，从ANSYS中将叶轮一个扇区模型的质量矩阵和刚度矩阵导入到Matlab中，其余的分析过程在Matlab中完成。在笛卡尔坐标系下，第n扇区的矩阵可以通过旋转变换得到，即：

式中，R_n为从基本扇区到第n扇区的旋转变换矩阵，且定义如下：

式中，I_S为维数为n_s的一个单位矩阵，且n_s为一个扇区中所包含的节点数，为Kronecker乘积,R_ne为节点的旋转变换矩阵。第n扇区与基本扇区对应节点间的旋转角度θ＝(n-1)2π/N，其中N为叶轮的扇区总数。将各扇区的质量矩阵和刚度矩阵的计算合并起来，则：

式中，R为一个以R_n为对角元素的对角矩阵。由于相邻扇区共享同一个边界，为得到最终的叶轮模型，需要将相邻边界的重复自由度消去，本质上也对应着一个变换矩阵T_e。因此，通过旋转变换得到的整个叶轮的质量和刚度矩阵为：

此时，系统的动力学方程为：

式中，M_E、C_E和K_E分别为系统的质量、阻尼和刚度矩阵，u为未经缩减的叶轮结构的位移向量，b为作用于结构上的外力向量，f(u)为子结构联接界面间的载荷向量。

(2)结构的自由度缩减

在式(5)模型的基础上，采用混合界面模态综合法对系统自由度进行缩减，对应的自由度变换方程为：

式中，u_i和u_b分别为系统的内部和界面自由度，p_k、p_a和u_b分别为保留的正则模态自由度、附加模态自由度和固定界面边界自由度，Φ_k、Φ_a和Φ_c分别为保留的正则模态集、剩余模态和约束模态。此时，系统的结构动力学方程为：

式中，和分别为自由度缩减后的系统的质量、阻尼和刚度矩阵，p为缩减系统的广义位移向量，为外力向量，为子结构联接界面间的载荷向量。

(3)大型稀疏矩阵的逆矩阵

为计算式(6)所示变换矩阵中的模态，涉及大型对称稀疏矩阵逆矩阵的计算问题，将矩阵等分为若干个子矩阵的形式，对于稀疏矩阵K_I，可将其分为2×2的子块，即：

则，该分块矩阵的逆矩阵为：

为计算式(9)中的各子块，需要先计算：

K_A＝(K₁₁)^-1

K_B＝(K₁₁)^-1K₁₂ (10)

将式(10)中的各元素带入式(9)中，可以得到矩阵K_I的逆矩阵。在实际的计算中，逆矩阵的各子块只需逐个计算。为进一步减少内存需求，可以将子矩阵K₁₁、K₁₂及K₂₂分为2×2的子块后进行计算，计算所需的内存空间可以进一步地降低。

(4)裂纹的建模

为模拟裂纹的呼吸效应，将裂纹界面节点自由度保留在式(7)所示的系统的自由度缩减模型中，并在界面节点之间定义接触对，如图3所示。然后在节点的局部坐标系下定义接触对的非线性接触力，

式中，和分别为目标面和接触面上的接触力，k为法向刚度惩罚系数，符号<·>为Macaulay括号，和分别为在局部坐标系下目标节点和接触节点的法向位移，g_i为接触对的法向初始间隙。在式(7)中考虑裂纹的呼吸效应，则，包含裂纹呼吸效应的系统非线性动力学方程为：

式中，M、C和K分别为裂纹叶轮结构缩减模型的质量、阻尼和刚度矩阵，q为系统的广义自由度向量，b为外载荷向量，f_nl(q)为作用在裂纹界面上的非线性接触力。

(5)谐振频率的估计

为对谐振频率进行统计分析，必须实现裂纹叶轮结构谐振频率的快速计算。因此，将裂纹非线性系统近似为双线性系统，在计算中涉及系统两种状态下的固有频率。首先，计算系统在状态一下系统的固有频率，状态一为不考虑裂纹呼吸效应的情况。其次，计算系统状态二下的固有频率，状态二为约束了裂纹界面间的法向相对运动，释放了界面间的切向相对运动。因此，结构的第i阶谐振频率f_bi可估计为：

式中，f_1i为系统在状态一的固有频率，f_2i为系统在状态二的固有频率。

(6)谐振频率的统计分析

在式(12)的系统非线性动力学方程的基础上考虑结构失谐因素，采用的是一种比例失谐模型，即随机改变各扇区叶片的弹性模量，第n扇区叶片的弹性模量E_n为：

E_n＝(1+δ_n)E₀ (13)

式中，E₀为叶片材料的名义弹性模量，δ_n为服从均匀分布的随机变量，且均值为0，标准差为σ，将该标准差定义为结构的失谐水平。进而，采用Monte Carlo模拟方法分析随机失谐对谐振频率的影响规律，分析流程如图4所示。

首先，在一定失谐水平下，由均匀分布随机产生N_m个失谐模式。根据子结构的失谐模式对系统刚度矩阵中对应于叶片结构的元素进行修改，再将失谐模式分别带入第5步中的状态一和状态二下的刚度矩阵中，并分别计算两种状态下的结构固有频率，进而得到结构的各谐振频率。反复执行上述过程，直至分析完所有失谐模式。其次，将得到的N_m组结构谐振频率进行统计分析，绘出各谐振频率的概率密度函数，进而确定谐振频率所服从的统计分布。最后，根据谐振频率的分布特征建立起结构谐振频率与裂纹深度和失谐水平间的定量关系，进而为离心叶轮裂纹定量诊断提供理论指导。

实施案例

该实施案例结合图2所示的某离心压缩机的叶轮模型验证了该发明的有效性，其中，1为端盖，2为叶片，3为轮盘。

某离心叶轮结构由端盖、叶片和轮盘焊接而成，根据工程经验，疲劳裂纹常出现在进气端外缘的叶片与端盖的焊接位置，该实施案例分析的裂纹即位于该处。裂纹叶轮结构的分析框架主要包含：前处理模块、分析模块和后处理模块等三个部分。

首先，通过旋转变换的方式得到整个叶轮结构的有限元模型，进而验证经过混合界面模态综合法缩减的系统的计算精度。在进行结果对比时，分析的对象为无裂纹叶轮结构，经过缩减后的系统包含1281个自由度，而原始的有限元模型包含302877个自由度，自由度缩减的比例为236倍。将系统缩减模型计算的前10阶固有频率与ANSYS的结果进行对比，计算结果以及相对误差如表1所示。可以看出前10阶固有频率的相对误差在0.1％以下，说明缩减模型具有较好的计算精度。

表1系统自由度缩减模型在固有频率上的计算误差

其次，在对系统进行自由度缩减时，涉及到大型对称稀疏矩阵逆矩阵的计算问题。计算中最大的一个稀疏矩阵的维度为119889，在Matlab中以双精度存储该矩阵的逆矩阵所需要的内存空间为：

S_m＝8×119889²Byte≈115GB (14)

如此巨大的内存需求使得大部分工作站无法满足要求，因此，必须对计算过程进行修改以大幅度降低内存的需求。通过采用本发明提出的大型对称稀疏矩阵逆矩阵的分块计算方法，计算过程中使用的最大内存约为30.2GB，内存的需求下降3.8倍，使得该求解可以在普通的工作站上完成。该方法为采用模态综合法分析工程复杂结构提供了计算上的可能性。

在系统自由度缩减模型的基础上，可以求解系统的频率响应曲线。图5为裂纹叶轮结构在2090至2130频率范围内的频率响应曲线，并包含了无裂纹以及20mm裂纹两种情况，可以看出，裂纹造成了结构谐振峰的分裂。

进而，研究裂纹叶轮结构谐振频率的统计规律性。根据图4所示的分析流程，需要计算结构的谐振频率。由于裂纹的存在，结构的响应本质上是非线性的。因此，将裂纹叶轮结构近似为双线性系统，以快速地计算结构的谐振频率。以具有20mm裂纹的叶轮为例，说明采用该方法的计算精度。图6为结构的第2和第3阶谐振频率随裂纹深度的变化情况。其中，标签中“BFA”表示的是双线性系统计算得到的结果，“Real”为实际的谐振频率，其后的数字为谐振频率的阶数。结果表明，双线性近似方法具有较好的计算精度。

其次，研究裂纹叶轮结构谐振频率的统计规律性。根据随机变量的均匀分布以及失谐水平，随机生成了25000种失谐模式，然后采用Monte Carlo模拟的方法对谐振频率进行统计分析。对分析得到的结果进行统计分析，图7和图8分别为结构的第2以及第5阶谐振频率的概率密度函数。图中红色的曲线表示的是标准正态分布的概率密度曲线，绿色的柱状图为实际的概率密度值。对分析结果进行正态分布检验，确定在均匀分布的失谐水平下，结构的谐振频率服从正态分布。因此，可以得到叶轮结构谐振频率与裂纹深度、失谐水平间的定量关系，进而得到基于谐振频率指标的裂纹诊断数据库，以实现离心叶轮裂纹的定量诊断问题。

本专利提出的离心压缩机裂纹叶轮结构振动建模与分析方法，能够方便地融合裂纹的呼吸效应以及失谐效应，并能够实现在有限计算资源条件下对裂纹叶轮结构动态响应以及谐振频率统计规律性的分析。该方法的分析结果可从理论上指导离心叶轮的裂纹定量诊断问题。

Claims (3)

1.一种离心压缩机裂纹叶轮结构的振动建模与分析方法，其特征在于，包括以下步骤：

1)首先，通过叶轮一个扇区的有限元模型旋转变换，得到整个结构的有限元模型；然后，运用混合界面模态综合法对有限元模型的系统自由度进行缩减，并在系统自由度进行缩减过程中，采用能够减少计算所需空间的大型稀疏逆矩阵求解方法进行计算，最后，在有限元模型的裂纹界面上定义接触对，以模拟裂纹的呼吸效应；

2)融合实际叶轮存在的制造误差及状态劣化失谐因素，获得裂纹叶轮结构在随机失谐模式下结构谐振频率的统计规律性，建立谐振频率与裂纹深度和失谐水平间的定量关系，具体步骤包括：

首先，计算裂纹叶轮结构的谐振频率，将非线性系统近似为双线性系统；计算中需要的信息为两个状态下的固有频率，状态一为不考虑裂纹界面的呼吸效应情况，状态二为约束了裂纹界面间的法向相对运动，释放了界面间的切向相对运动，则裂纹叶轮结构的第i阶谐振频率f_bi可估计为：

$f_{bi}=\frac{2f_{1i}{f}_{2i}}{f_{1i}+f_{2i}}$

式中，f_1i为系统在状态一的固有频率，f_2i为系统在状态二的固有频率；

其次，引入简化的失谐模式来模拟失谐效应，随机改变各扇区叶片的弹性模量，则第n扇区叶片的弹性模量E_n为：

E_n＝(1+δ_n)E₀

式中，E₀为叶片材料的名义弹性模量，δ_n为服从均匀分布的随机变量，且均值为0，标准差为σ，将该标准差定义为结构的失谐水平；

最后，采用统计分析的方法分析随机失谐对谐振频率的影响规律，所述采用统计分析的方法分析随机失谐对谐振频率的影响规律具体包括：

a、先对叶轮结构的失谐水平进行估计，并随机生成较多数量的随机失谐模式；

b、采用Monte Carlo模拟方法对谐振频率进行统计分析，计算在每一种失谐模式下的结构谐振频率，根据每一种失谐模式下的结果得到各谐振频率的概率分布曲线；

c、分析在均匀分布的失谐模式下的谐振频率的统计分布，并研究分布的均值、方差与失谐水平间的关系；

d、根据谐振频率的统计规律性，建立谐振频率与裂纹深度和失谐水平间的定量关系。

2.根据权利要求1所述的一种离心压缩机裂纹叶轮结构的振动建模与分析方法，其特征在于，所述步骤1)中通过叶轮一个扇区的有限元模型旋转变换，得到整个结构的有限元模型方法如下：

首先，获得叶轮一个扇区有限元模型的质量和刚度矩阵分别为M_s和K_s，在笛卡尔坐标系下，第n扇区的质量矩阵M_n和刚度矩阵K_n通过旋转变换的方式得到：

$M_n=R_n^T{M}_s{R}_n,K_n=R_n^T{K}_s{R}_n$

式中，R_n为从基本扇区到第n扇区的旋转变换矩阵，且R_n为：

$R_n=I_S\&CircleTimes;R_{ne}$

式中，I_S为维数为n_s的一个单位矩阵，且n_s为一个扇区中所包含的节点数，为Kronecker乘积,R_ne为节点的旋转变换矩阵；第n扇区与基本扇区对应节点间的旋转角度θ＝(n-1)2π/N，其中N为叶轮的扇区总数；将各扇区的质量矩阵和刚度矩阵的计算合并起来，则：

$M_R=R^T(I_N\&CircleTimes;M_s)R,K_R=R^T(I_N\&CircleTimes;K_s)R$

式中，M_R、K_R为旋转变换后的结构质量矩阵和刚度矩阵，R为一个以R_n为对角元素的pseudo-block对角矩阵；

其次，采用混合界面模态综合法对系统的自由度进行缩减，该混合界面模态综合方法对应的自由度变换方程为：

$u=\left\{\begin{array}{c}{u}_i\\ u_b\end{array}\right\}=\left[\begin{array}{ccc}{\mathrm{\&Phi;}}_k& {\mathrm{\&Phi;}}_a& {\mathrm{\&Phi;}}_c\\ 0& 0& I\end{array}\right]\left\{\begin{array}{c}{p}_k\\ p_a\\ u_b\end{array}\right\}={\mathrm{\&Phi;}}_{H}p$

式中，u_i和u_b分别为系统的内部和界面自由度，p_k、p_a和u_b分别为保留的正则模态自由度、附加模态自由度和固定界面边界自由度，Φ_k、Φ_a和Φ_c分别为保留的正则模态集、剩余模态和约束模态；

然后，将裂纹界面上的节点自由度保留在系统的缩减模型中，并在界面节点之间定义接触，以模拟裂纹界面的呼吸效应；接触对的非线性接触力是在节点的局部坐标系下定义的，即：

$\left\{\begin{array}{c}{f}_T^i=-k\&CenterDot;<u_T^i+u_C^i-g_i>\\ f_C^i=f_T^i\end{array}\right.$

式中，和分别为目标面和接触面上的接触力，k为法向刚度惩罚系数，符号<·>为Macaulay括号，和分别为在局部坐标系下目标节点和接触节点的法向位移，g_i为接触对的法向初始间隙；则包含裂纹呼吸效应的系统非线性动力学方程为：

$M\stackrel{\&CenterDot;\&CenterDot;}{q}+C\stackrel{\&CenterDot;}{q}+Kq=b+f_{nl}\left(q\right)$

式中，M、C和K分别为缩减模型的质量、阻尼和刚度矩阵，q为缩减系统的广义自由度向量，b为外载荷向量，f_nl(q)为作用在裂纹界面上的非线性接触力。

3.根据权利要求1所述的一种离心压缩机裂纹叶轮结构的振动建模与分析方法，其特征在于，所述采用能够减少计算所需空间的大型对称稀疏矩阵逆矩阵的求解方法，该求解方法包括：

在计算大型对称稀疏矩阵K_I的逆矩阵时，将矩阵K_I等分为2×2的子块，即：

$K_I=\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{12}^T& K_{22}\end{array}\right]$

则，其逆矩阵可以表示为：

${\left(K_I\right)}^{-1}=\left[\begin{array}{cc}{K}_A+K_B{K}_C{K}_B^T& -K_B{K}_C\\ -{\left(K_B{K}_C\right)}^T& K_C\end{array}\right]$

式中，

K_A＝(K₁₁)^-1

K_B＝(K₁₁)^-1K₁₂

$K_C={(K_{22}-K_{12}^T{\left(K_{11}\right)}^{-1}{K}_{12})}^{-1}$

其中，子矩阵K₁₁、K₁₂及K₂₂分为2×2的子块后进行计算。