CN113705045B - 一种基于代理模型的转静子系统碰摩可靠性分析方法 - Google Patents

Abstract

本发明公开了一种基于代理模型的转静子系统碰摩可靠性分析方法，包括步骤1：建立转静子碰摩的有限元仿真模型；步骤2：找出设计阶段、加工及装配、工作运转环节中存在的不确定参数，通过灵敏度分析确定影响转静子碰摩的主要不确定性来源及各参数的概率分布类型，输入参数；步骤3：利用自适应实验设计算法，在基于极限状态曲面附近添加精心选择的样本，构建代理模型来替代有限元仿真模型；步骤4：利用训练后的代理模型，在给定的初始转速下预测MCS样本点的功能响应，根据功能响应结果计算失效概率；步骤5：在所考虑的转速范围内，迭代增加转速，重复步骤4，依次计算每个转速下的失效概率，最终得到失效概率曲线。

Description

一种基于代理模型的转静子系统碰摩可靠性分析方法

技术领域

本发明涉及可靠性技术领域，特别涉及一种基于代理模型的转静子系统碰摩可靠性分析方法。

背景技术

旋转机械如航空发动机、船用燃气轮机和压缩机等由于高性能要求，使得转子和定子之间的设计间隙越来越小，这很可能导致设备运行过程中发生碰摩故障。碰摩作为次生故障不可避免，旋转机械的设计很难平衡无碰摩和间隙减小这两方面的矛盾。转静子之间的碰摩会产生高冲击力和摩擦力，影响系统安全运行，很可能导致叶片断裂或者机匣包容性问题，引起严重振动甚至灾难性事故。

转静件碰摩引起的振动响应为复杂的非线性振动，会表现出极其丰富的动力学行为，如偏摩、全周碰摩和反向涡动，系统出现准周期甚至混沌与分岔现象，因此转静子碰摩故障的预防一直备受重视。实际上，影响转静子碰摩的参数具有明显的随机性，设计、加工、装配及旋转机械运转过程中的参数具有不确定性，某一转速甚至某一时刻转静子各参数都是动态变化的。

因此，碰摩的发生具有随机性，为一概率事件，需要进行可靠性评估。现有的传统蒙特卡罗及抽样方法是一种用于简单结构响应概率分析的有效技术，当应用于复杂非线性机械系统动态、多目标和多学科的可靠性问题时，其计算效率将不可接受。转静子碰摩可靠性分析即属于这一类复杂问题，蒙特卡罗及抽样方法在转静子碰摩概率分析的计算量将会非常大甚至无法计算。

为了提高计算效率，在可靠性分析领域，先进的代理模型被开发用来解决复杂结构的可靠性问题，其中代理模型结合蒙特卡罗仿真的方法是一种有效的可靠性分析技术，在土木工程、电力系统等领域被广泛采用。目前有关转静子碰摩失效概率的分析尚无快速有效的评估方法，如何针对转静子碰摩进行不确定分析及量化，以度量碰摩发生的概率这一问题亟待需要得到解决，进而有利于稳健的产品设计及旋转机械的可靠性提升。

发明内容

为了克服现有技术中的不足，本发明提供一种基于代理模型的转静子系统碰摩可靠性分析方法，可实现不同转频不确定概率参数区间内转静子碰摩概率的快速有效计算。

为了达到上述发明目的，解决其技术问题所采用的技术方案如下：

一种基于代理模型的转静子系统碰摩可靠性分析方法，包括以下步骤：

步骤1：建立转静子碰摩的有限元仿真模型；

步骤2：找出设计阶段、加工及装配、工作运转环节中存在的不确定参数，通过灵敏度分析确定影响转静子碰摩的主要不确定性来源及各参数的概率分布类型，输入参数；

步骤3：利用自适应实验设计算法，在基于极限状态曲面附近添加精心选择的样本，构建代理模型来替代有限元仿真模型；

步骤4：利用训练后的代理模型，在给定的初始转速下预测MCS样本点的功能响应，根据功能响应结果计算失效概率；

步骤5：在所考虑的转速范围内，迭代增加转速，重复步骤4，依次计算每个转速下的失效概率，最终得到失效概率曲线。

进一步的，步骤1中所述的转静子碰摩模型通过有限元方法建立的动力学运动微分方程如下：

其中，M、C、G及K分别为系统质量、阻尼、陀螺及刚度矩阵，H为Heaviside函数，F_u为转子不平衡激励，q、及/>分别为位移、速度和加速度向量，ω为转子转速，x为转子接触点在水平方向振动位移，g为转静子初始间隙，t为时间；

根据赫兹接触与库伦摩擦定律，转静子单圆盘位置单处单点碰摩力表达式为：

式中，F_n、F_t分别为法向碰撞力与切向摩擦力，k_c为接触刚度，μ为摩擦系数；

在有限元模型对应的圆盘节点位置，转子圆盘的不平衡激励计算表达式为：

式中，m_ei、e_i(i＝1,2)分别为两圆盘偏心质量及偏心距。

进一步的，步骤2中，考虑转子系统自身材料性能中的接触刚度k_c、装配误差影响的转静子初始间隙g、不平衡量u这三个参数的不确定性，假设接触刚度k、初始间隙g为高斯分布，初始间隙变化范围随机性更强，假设为beta分布，不确定参数的概率描述为：

其中，和/>为接触刚度、不平衡量与间隙的平均值，/>δ_X和δ_J分别为三参数的不确定度，不确定度与平均值的乘积则为随机参数的标准差，ξ₁、ξ₂和ξ₃则为三个随机变量，其中变量ξ₁和ξ₂服从高斯分布，变量ξ₃服从beta分布。然后指定参数均值、不确定度，将随机样本点代入对应参数的概率发分布中生成对应的概率分布。

进一步的，在步骤3中，利用自适应算法采样构建训练Kriging代理模型的步骤包括：

步骤3.1：通过LHS生成一个小的初始样本X＝{x₁,...,x_N0}，N₀为向量的维数,样本能够有效覆盖足够广泛的设计变量,通过有限元模型计算得到真实响应，评估极限状态函数响应Y＝{y₁,...,y_N0}；

步骤3.2：基于步骤3.1的仿真实验设计{X,Y}训练一个初始克里金模型

克里金模型获得由4个子步确定：

选择克里金回归模型的基函数；

选择一个相关函数；

如果超参数θ未知，通过最大似然估计优化确定；

使用θ的最佳值，计算其余未知克里金参数(如β)，如有必要，计算σ²；克里金模型由以下等式描述：

G(x)＝f^T(x)β+z(x，ω) (7)

式中，第一项f^T(x)β为回归模型，即克里金模型的均值或者趋势，β为回归系数向量，f(x)为随机变量的多项式，提供模拟的全局近似，可以为0阶、1阶或是2阶多项式，一般将其设为常数，常数值的大小并不影响模型的近似精度，回归模型包含P个任意函数{f_j；j＝1,...,P}和相应的系数{β_j；j＝1,...,P}；第二项z(x，ω)为平稳高斯过程，反映局部偏差的近似，决定模型的准确性，其方差为σ²和均值为0，基本概率空间用ω表示，并根据相关函数R及其超参数θ定义，z(x，ω)的协方差矩阵表示为：

cov[z(x_i)，z(x_j)＝σ²R(x_i，x_j；θ) (8)

式中，R(x_i，x_j)是带有参数θ的相关函数，描述了输出空间中两个样本点x_i和x_j之间的相关性，相关函数选用高斯核函数，表达式为：

式中，与/>分别为样本x_i和x_j的第k个分量，θ_k为相关性参数，由相关函数的极大似然估计确定，即：

式中，|R|表示为相关函数R的行列式；

通过θ便可确定克里金模型表达式中β和σ²的估计值：

式中，F是克里金超模型均值的观察(设计)矩阵，即：

F_ij＝f_j(x⁽ⁱ⁾)，i＝1，...，N₀；j＝0，...，P (13)

由此便建立了初始克里金模型；

步骤3.3：生成N_MC个待预测候选参数样本，通过训练的初始克里金模型预测相应的碰摩诱导的代理模型响应/>

新的待测点x的预测响应值为：

式中，为是未知点x和所有已知实验点之间的相关向量；

得到的克里金模型均值和方差的预测值为：

其中，u(x)＝F^TR^-1r(x)-f(x)；

步骤3.4：基于适当的学习函数，选择下一个最佳样本s*添加到实验设计X中；

基于错误分类概念的U-函数定义如下：

当代理模型的符号和基础极限状态函数的符号不匹配时，就会发生错误分类，相应的错误分类概率是:

P_m(x)＝Φ(-U(x)) (18)

其中，Φ是标准高斯变量的累计分布函数，从集合S中选择下一候选样本作为最有可能被代理极限状态函数错误分类为安全/失败的样本：

步骤3.5：检查是否满足给定的收敛标准，如果是，请跳到步骤3.7，否则继续步骤；

收敛标准终止了克里金元模型实验设计中样本的添加，从而终止了故障概率估计准确性的提高，与U函数相关的标准收敛标准定义为当min_iU(s⁽ⁱ⁾)>2时，迭代停止；

步骤3.6：在代理模型的实验设计中加入s*和相应的极限状态函数响应y*＝g(s*)，返回步骤2；

步骤3.7：用最终极限状态函数替代该代理模型便可以通过蒙特卡罗模拟估算失效概率。

进一步的，步骤4：利用训练后的克里金Kriging代理模型，在给定的初始转速下预测MCS中所有样本点的功能响应，根据功能响应结果计算失效概率；

在给定的参数不确定空间内，求解转子稳态振动响应；当盘位置的响应幅值大于初始间隙g时，发生碰摩，转子系统被认为是失效的；相反，如果盘位置的响应幅值小于初始间隙g，则不会发生碰摩，转子系统被认为是安全的；因此，转子碰摩失效概率可以定义为：

P_f＝∫_s(x，g)≤0f(Z)dZ (20)

式中，Z代表系统随机振动响应，f(Z)表示随机响应的联合概率函数，S(x,g)为间隙判断函数，定义为极限状态函数：

S(x，g)＝g-x (21)

定义通过MCS对自适应克里金模型进行仿真模拟，根据中心极限则失效概率可近似为：

式中，N_MC为所述蒙特卡罗模拟的样本点数目。

本发明由于采用以上技术方案，使之与现有技术相比，具有以下的优点和积极效果：

本发明提出了一种快速计算转静子碰摩可靠度的计算方法，考虑到转子系统自身材料性能、加工装配公差、载荷等产生的不确定性因素，首先确定随机参数的概率分布类型，采用初始采样点策略和主动学习函数，构建主动学习代理模型替代原始极限状态函数，结合蒙特卡罗模拟来评估转静子碰摩失效概率，同时给出随机参数对碰摩诱导的非线性振动响应的不确定量化，该方法为转静子碰摩失效概率计算提供了新的思路，且代理模型的采用大大减少了计算样本点数，有效的提高了计算效率。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍。显而易见，下面描述中的附图仅仅是本发明的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。附图中：

图1为本发明提出的基于克里金模型的转静子碰摩可靠度分析流程图；

图2是双盘单转子有限元碰摩仿真示意图；

图3是转静子单点碰摩示意图；

图4是Case1转静子碰摩失效概率曲线图；

图5是Case2转静子碰摩失效概率曲线图；

图6是Case3转静子碰摩失效概率曲线图；

图7是三种Beta分布频率分布直方图；

图8是间隙分布为Beta(5,3)转静子碰摩失效概率曲线图；

图9是间隙分布为Beta(5,5)转静子碰摩失效概率曲线图；

图10是间隙分布为Beta(3,5)转静子碰摩失效概率曲线图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1-10所示，本实施例公开了一种基于代理模型(以克里金Kriging代理模型为例)的转静子系统碰摩可靠性分析方法，包括以下步骤：

步骤1：建立转静子碰摩的有限元仿真模型；

步骤2：找出设计阶段、加工及装配、工作运转等环节中存在的不确定参数，通过灵敏度分析确定影响转静子碰摩的主要不确定性来源及各参数的概率分布类型，输入参数；

步骤3：利用自适应实验设计算法，在基于极限状态曲面(S(x)＝0)附近添加精心选择的样本，构建代理模型来替代有限元仿真模型；

步骤4：利用训练后的代理模型，在给定的初始转速下预测MCS样本点的功能响应，根据功能响应结果计算失效概率；

步骤5：在所考虑的转速范围内，迭代增加转速，重复步骤4，依次计算每个转速下的失效概率，最终得到失效概率曲线。

在上述技术方案中，步骤1的具体步骤如下：

图2为双盘单转子有限元碰摩仿真示意图，只考虑左盘单点碰摩，步骤1中所述的转静子碰摩模型通过有限元方法建立的动力学运动微分方程如下：

其中，M、C、G及K分别为系统质量、阻尼、陀螺及刚度矩阵，H为Heaviside函数，F_u为转子不平衡激励，q、及/>分别为位移、速度和加速度向量，ω为转子转速，x为转子接触点在水平方向振动位移，g为转静子初始间隙，t为时间；

图3为转静子单点碰摩示意图，根据赫兹接触与库伦摩擦定律，转静子单圆盘位置单处单点碰摩力表达式为：

式中，F_n、F_t分别为法向碰撞力与切向摩擦力，k_c为接触刚度，μ为摩擦系数；

在有限元模型对应的圆盘节点位置，转子圆盘的不平衡激励计算表达式为：

式中，，m_ei、e_i(i＝1,2)分别为两圆盘偏心质量及偏心距。

进一步的，步骤2中，考虑转子系统自身材料性能中的接触刚度k_c、装配误差影响的转静子初始间隙g、不平衡量u这三个参数的不确定性，三个参数的概率分布已知且可以为任意的分布类型，如均匀分布、高斯分布及beta分布等。假设接触刚度k、初始间隙g为高斯分布，初始间隙变化范围随机性更强，假设为beta分布，不确定参数的概率描述为：

其中，和/>为接触刚度、不平衡量与间隙的平均值，/>δ_X和δ_J分别为三参数的不确定度，不确定度与平均值的乘积则为随机参数的标准差，ξ₁、ξ₂和ξ₃则为三个随机变量，其中变量ξ₁和ξ₂服从高斯分布，变量ξ₃服从beta分布。然后指定参数均值、不确定度，将随机样本点代入对应参数的概率发分布中生成对应的概率分布。

进一步的，在步骤3中，利用自适应算法采样构建训练Kriging代理模型的步骤包括：

步骤3.1：通过LHS生成一个小的初始样本X＝{x₁,...,x_N0}，N₀为向量的维数，样本能够有效覆盖足够广泛的设计变量,通过有限元模型计算得到真实响应，评估极限状态函数响应Y＝{y₁,...,y_N0}；

步骤3.2：基于步骤3.1的仿真实验设计{X,Y}训练一个初始克里金模型

克里金模型获得由4个子步确定：

选择克里金回归模型的基函数；

选择一个相关函数；

如果超参数θ未知，通过最大似然估计优化确定；

使用θ的最佳值，计算其余未知克里金参数(如β)，如有必要，计算σ²；克里金模型由以下等式描述：

G(x)＝f^T(x)β+z(x，ω) (7)

式中，第一项f^T(x)β为回归模型，即克里金模型的均值或者趋势，β为回归系数向量，f(x)为随机变量的多项式，提供模拟的全局近似，可以为0阶、1阶或是2阶多项式，一般将其设为常数，常数值的大小并不影响模型的近似精度，回归模型包含P个任意函数{f_j；j＝1,...,P}和相应的系数{β_j；j＝1,...,P}；第二项z(x，ω)为平稳高斯过程，反映局部偏差的近似，决定模型的准确性，其方差为σ²和均值为0，基本概率空间用ω表示，并根据相关函数R及其超参数θ定义，z(x，ω)的协方差矩阵表示为：

cov[z(x_i)，z(x_j)＝σ²R(x_i，x_j；θ) (8)

式中，5[_L[_M是带有参数θ的相关函数，描述了输出空间中两个样本点x_i和x_j之间的相关性，相关函数选用高斯核函数，表达式为：

式中，与/>分别为样本x_i和x_j的第k个分量，θ_k为相关性参数，由相关函数的极大似然估计确定，即：

式中，|R|表示为相关函数R的行列式；

通过θ便可确定克里金模型表达式中β和σ²的估计值：

式中，F是克里金超模型均值的观察(设计)矩阵，即：

F_ij＝f_j(x⁽ⁱ⁾)，i＝1，...，N₀；j＝0，...，P (13)

由此便建立了初始克里金模型；

步骤3.3：生成N_MC个待预测候选参数样本，通过训练的初始克里金模型预测相应的碰摩诱导的代理模型响应/>

新的待测点x的预测响应值为：

式中，为是未知点x和所有已知实验点之间的相关向量；

得到的克里金模型均值和方差的预测值为：

其中，u(x)＝F^TR^-1r(x)-f(x)；

步骤3.4：基于适当的学习函数，选择下一个最佳样本s*添加到实验设计X中；

基于错误分类概念的U-函数定义如下：

当代理模型的符号和基础极限状态函数的符号不匹配时，就会发生错误分类，相应的错误分类概率是:

P_m(x)＝Φ(-U(x)) (18)

其中，Φ是标准高斯变量的累计分布函数，从集合S中选择下一候选样本作为最有可能被代理极限状态函数错误分类为安全/失败的样本：

步骤3.5：检查是否满足给定的收敛标准，如果是，请跳到步骤3.7，否则继续步骤；

收敛标准终止了克里金元模型实验设计中样本的添加，从而终止了故障概率估计准确性的提高，Echard等人提出的与U函数相关的标准收敛标准定义为当min_iU(s⁽ⁱ⁾)>2时，迭代停止；

步骤3.6：在代理模型的实验设计中加入s*和相应的极限状态函数响应y*＝g(s*)，返回步骤2；

步骤3.7：用最终极限状态函数替代该代理模型便可以通过蒙特卡罗模拟估算失效概率。

进一步的，步骤4：利用训练后的克里金Kriging代理模型，在给定的初始转速下预测MCS中所有样本点的功能响应，根据功能响应结果计算失效概率；

在给定的参数不确定空间内，求解转子稳态振动响应；当盘位置的响应幅值大于初始间隙g时，发生碰摩，转子系统被认为是失效的；相反，如果盘位置的响应幅值小于初始间隙g，则不会发生碰摩，转子系统被认为是安全的；因此，转子碰摩失效概率可以定义为：

P_f＝∫_s(x，g)≤0f(Z)dZ (20)

式中，Z代表系统随机振动响应，f(Z)表示随机响应的联合概率函数，S(x,g)为间隙判断函数，定义为极限状态函数：

S(x，g)＝g-x (21)

定义通过MCS对自适应克里金模型进行仿真模拟，根据中心极限则失效概率可近似为：

式中，N_MC为所述蒙特卡罗模拟的样本点数目。

实施例：

为了进一步说明本申请的可靠性分析方法，利用案例进行具体介绍。附图说明：图1为本发明提出的基于克里金模型的转静子碰摩可靠度分析流程。

图2所示，为双盘单转子有限元碰摩仿真示意图，只考虑左盘单点碰摩，步骤1中所述的转静子碰摩模型通过有限元方法建立的动力学运动微分方程如下：

其中，M、C、G及K分别为系统质量、阻尼、陀螺及刚度矩阵，H为Heaviside函数，F_u为转子不平衡激励，T、及/>分别为位移、速度和加速度向量，ω为转子转速，[为转子接触点在水平方向振动位移，J为转静子初始间隙，W为时间；

图3为转静子单点碰摩示意图，根据赫兹接触与库伦摩擦定律，转静子单圆盘位置单处单点碰摩力表达式：

式中，F_n、F_t分别为法向碰撞力与切向摩擦力，k_c为接触刚度，μ为摩擦系数；

在有限元模型对应的圆盘节点位置，转子圆盘的不平衡激励计算表达式为：

式中，m_ei、e_i(i＝1,2)分别为两圆盘偏心质量及偏心距。

接触刚度k_c、不平衡量u及转静子初始间隙g单一参数不确定统计信息见表1。

案例	不确定参数	均值	不确定度	分布类型
					Case1	u	180g·mm	3％	高斯
Case2	kc	2×105N/m	1.5％	高斯
					Case3	g	3×10-4N/m	10％	高斯

表1转子单一不确定参数统计信息

在三组单参数的不确定度下，失效概率随频率变化曲线如图4-图6所示。AKM的结果与MCS结果在全局频率范围内十分接近。在表1三种案例下，当转子通过临界转速时，其不平衡振动的响应幅度均会大于转子/定子的初始间隙。因此，碰摩失效概率为1。当不平衡激励频率远离共振区时，不平衡响应幅值始终小于初始间隙，碰摩失效概率为0。除第2种案例外，当不平衡激励频率接近临界转速时，才会发生一定概率的碰摩。图5所示的第2种案例，接触刚度的变化不会影响碰摩的出现条件，这是由公式(23)定义的极限状态函数决定的。多参数不确定统计信息见表2。

表2转子多源不确定参数统计信息

Beta分布定义在区间[0，1]上，区间由两个形参数参数化，表示为α以及β。在本实施例中，β分布的形状参数分别为(5，3)、(5，5)、(3，5)。图7中绘制了相应的概率密度函数。值得注意的是，不同的形状参数可以使beta分布模拟其他先验分布，如均匀分布。

图8-图10显示了根据AKM与MCS两种方法的统计结果计算出的碰摩失效概率曲线。从这三种情况下的失效概率曲线可以得出相同的规律：第一，失效概率在第一临界转速附近达到峰值。第二，当转子频率远离第一临界转速(低于或高于第一临界转速)时，转子振动幅值不会超过初始随机间隙，失效概率减小到0，不同的是，随着间隙概率分布趋于一个较小的值，共振转速范围内的碰摩概率增大，最大失效概率由0.16增大到0.71，碰摩故障发生的频率范围也增大。这表明间隙对碰摩失效概率至关重要。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (4)

1.一种基于代理模型的转静子系统碰摩可靠性分析方法，其特征在于，包括以下步骤：

步骤1：建立转静子碰摩的有限元仿真模型；

步骤2：找出设计阶段、加工及装配、工作运转环节中存在的不确定参数，通过灵敏度分析确定影响转静子碰摩的主要不确定性来源及各参数的概率分布类型，输入参数；

步骤3：利用自适应实验设计算法，在基于极限状态曲面附近添加精心选择的样本，构建代理模型来替代有限元仿真模型；

在步骤3中，利用自适应算法采样构建训练Kriging代理模型的步骤包括：

步骤3.1：通过LHS生成一个小的初始样本X＝{x₁,...,x_N0}，N₀为向量的维数，样本能够有效覆盖足够广泛的设计变量,通过有限元模型计算得到真实响应，评估极限状态函数响应Y＝{y₁,...,y_N0}；

步骤3.2：基于步骤3.1的仿真实验设计{X,Y}训练一个初始克里金模型

克里金模型获得由4个子步确定：

选择克里金回归模型的基函数；

选择一个相关函数；

如果超参数θ未知，通过最大似然估计优化确定；

使用θ的最佳值，计算其余未知克里金参数β，计算σ²；

克里金模型G(x)由以下等式描述：

G(x)＝f^T(x)β+z(x,ω) (7)

式中，第一项f^T(x)β为回归模型，即克里金模型的均值或者趋势，β为回归系数向量，f(x)为随机变量的多项式，提供模拟的全局近似，可以为0阶、1阶或是2阶多项式，将其设为常数，常数值的大小并不影响模型的近似精度，回归模型包含P个任意函数{f_j；j＝1,...,P}和相应的系数{β_j；j＝1,...,P}；第二项z(x,ω)为平稳高斯过程，反映局部偏差的近似，决定模型的准确性，其方差为σ²和均值为0，基本概率空间用ω表示，并根据相关函数R及其超参数θ定义，z(x,ω)的协方差矩阵表示为：

cov[z(x_i),z(x_j)]＝σ²R(x_i,x_j；θ) (8)

式中，R(x_i,x_j)是带有参数θ的相关函数，描述了输出空间中两个样本点x_i和x_j之间的相关性，相关函数选用高斯核函数，表达式为：

式中，与/>分别为样本x_i和x_j的第k个分量，θ_k为相关性参数，由相关函数的极大似然估计确定，即：

式中，|R|表示为相关函数R的行列式，m为矢量维数；

通过θ便可确定克里金模型表达式中β和σ²的估计值：

式中，F是克里金超模型均值的观察矩阵，即：

F_ij＝f_j(x⁽ⁱ⁾),i＝1,…,N₀；j＝0,…,P (13)

由此便建立了初始克里金模型；

步骤3.3：生成N_MC个待预测候选参数样本，通过训练的初始克里金模型预测相应的碰摩诱导的代理模型响应/>

新的待测点x的预测响应值为：

式中，为是未知点x和所有已知实验点之间的相关向量；

得到的克里金模型均值和方差的预测值为：

其中，u(x)＝F^TR^-1r(x)-f(x)；

步骤3.4：基于适当的学习函数，选择下一个最佳样本s*添加到实验设计X中；

基于错误分类概念的U-函数定义如下：

当代理模型的符号和基础极限状态函数的符号不匹配时，就会发生错误分类，相应的错误分类概率是:

P_m(x)＝Φ(-U(x)) (18)

其中，Φ是标准高斯变量的累计分布函数，从集合S中选择下一候选样本作为最有可能被代理极限状态函数错误分类为安全/失败的样本：

步骤3.5：检查是否满足给定的收敛标准，如果是，请跳到步骤3.7，否则继续步骤；

收敛标准终止了克里金元模型实验设计中样本的添加，从而终止了故障概率估计准确性的提高，与U函数相关的标准收敛标准定义为当min_iU(s⁽ⁱ⁾)>2时，迭代停止；

步骤3.6：在代理模型的实验设计中加入s*和相应的极限状态函数响应y*＝g(s*)，返回步骤2；

步骤3.7：用最终极限状态函数替代该代理模型便可以结合蒙特卡罗模拟估算失效概率；

步骤4：利用训练后的代理模型，在给定的初始转速下预测MCS样本点的功能响应，根据功能响应结果计算失效概率；

步骤5：在所考虑的转速范围内，迭代增加转速，重复步骤4，依次计算每个转速下的失效概率，最终得到失效概率曲线。

2.根据权利要求1所述的一种基于代理模型的转静子系统碰摩可靠性分析方法，其特征在于，步骤1中所述的转静子碰摩模型通过有限元方法建立的动力学运动微分方程如下：

其中，M、C、G及K分别为系统质量、阻尼、陀螺及刚度矩阵，H为Heaviside函数，F_u为转子不平衡激励，q、及/>分别为位移、速度和加速度向量，ω为转子转速，x为转子接触点在水平方向振动位移，g为转静子初始间隙，t为时间；

根据赫兹接触与库伦摩擦定律，转静子单圆盘位置单处单点碰摩力表达式为：

式中，F_n、F_t分别为法向碰撞力与切向摩擦力，k_c为接触刚度，μ为摩擦系数；

在有限元模型对应的圆盘节点位置，转子圆盘的不平衡激励计算表达式为：

式中，m_ei、e_i(i＝1,2)分别为两圆盘偏心质量及偏心距。

3.根据权利要求1所述的一种基于代理模型的转静子系统碰摩可靠性分析方法，其特征在于，步骤2中，考虑转子系统自身材料性能中的接触刚度k_c、装配误差影响的转静子初始间隙g、不平衡量u这三个参数的不确定性，假设接触刚度k、初始间隙g为高斯分布，初始间隙变化范围随机性更强，假设为beta分布，不确定参数的概率描述为：

式中，和/>为接触刚度、不平衡量与间隙的平均值，/>δ_u和δ_g分别为三参数的不确定度，不确定度与平均值的乘积则为随机参数的标准差，ξ₁、ξ₂和ξ₃则为三个随机变量，其中变量ξ₁和ξ₂服从高斯分布，变量ξ₃服从beta分布,然后指定参数均值、不确定度，将随机样本点代入对应参数的概率发分布中生成对应的概率分布。

4.根据权利要求1所述的一种基于代理模型的转静子系统碰摩可靠性分析方法，其特征在于，步骤4：利用训练后的克里金Kriging代理模型，在给定的初始转速下预测MCS中所有样本点的功能响应，根据功能响应结果计算失效概率；

在给定的参数不确定空间内，求解转子稳态振动响应；当盘位置的响应幅值大于初始间隙g时，发生碰摩，转子系统被认为是失效的；相反，如果盘位置的响应幅值小于初始间隙g，则不会发生碰摩，转子系统被认为是安全的；因此，转子碰摩失效概率可以定义为：

P_f＝∫_S(x,g)≤0f(Z)dZ (20)

式中，Z代表系统随机振动响应，f(Z)表示随机响应的联合概率函数，S(x,g)为间隙判断函数，定义为极限状态函数：

S(x,g)＝g-x (21)

定义通过MCS对自适应克里金模型进行仿真模拟，根据中心极限则失效概率可近似为：

式中，N_MC为所述蒙特卡罗模拟的样本点数目。