CN104008240B - 一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法 - Google Patents

Abstract

一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，它有七大步骤：一、针对在轨空间柔性齿轮机构中轴是弹性体，建立柔性轴的动力学方程，并采用有限元方法分析；二、建立齿轮动力学方程；三、建立在轨空间柔性齿轮机构中柔性轴与齿轮耦合的非线性动力学方程；四、进行空间柔性齿轮机构传动激励分析；五、采用Newmark算法计算动载荷系数；六、在轨空间柔性齿轮机构中齿轮弯曲应力可靠性极限状态函数的建立；七、空间柔性齿轮机构动态可靠度与故障率分析。本发明避免了传统的传递矩阵法和振型叠加法因非线性因素导致的可靠度计算精度低问题，对于提高空间柔性机构等航天器可靠性具有工程实用价值。本发明同样适用于其他柔性齿轮机构。

Description

一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法

技术领域

本发明涉及一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，它涉及到弹性体动力学及在轨空间柔性齿轮机构时变可靠性仿真分析，属于航天机构产品可靠性设计分析领域。

背景技术

柔性齿轮机构是典型的在轨空间柔性机构，其可靠性水平直接关系到在轨空间卫星等航天器的安全和任务成功。因此，其在轨可靠性设计分析是保证在轨空间机构可靠运转的前提和关键技术之一。柔性齿轮机构中齿轮啮合过程是一个典型的动态响应过程，故障率呈现随时间变化的时变特性，传统采用PHI2方法只能得到近似的柔性齿轮机构时变可靠度，对于可靠性工程中的关键参数故障率并没有明确的给出或者只给出了故障率的计算结果，没能给出故障率的变化规律，这些基于统计学的方法已不适用于柔性齿轮机构的动态故障率分析计算。另外，在轨空间柔性齿轮机构大多采用挠性转子，不考虑轴横向振动和柔性齿轮轴弯曲影响的齿轮啮合动力学模型已不太适用，采用有限元动力学方法建立的齿轮动力学模型，虽然考虑了齿轮轴和齿轮耦合影响，但求解复杂，工程实用性不强，因此需要建立考虑长转子柔性轴弹性形变的齿轮和轴刚柔耦合的动力学方程。其次，考虑柔性齿轮传动机构的刚度矩阵和内外部激励都具有动态特性和非线性，传统的传递矩阵法和振动叠加法忽略了方程中的非线性因素，导致计算精度不够，采用Newmark算法求解动力响应，能一定程度上提高计算解的精度。最后，通过高度非线性刚柔耦合动力学方程建立的可靠性极限状态函数也是一个高维非线性的方程，kriging模型作为一种新型的响应面模型技术对未知信息的模拟精度和整体性要远远优于响应面法模拟技术，在可靠性设计分析中，被视为最优的线性无偏估计。

综上所述，考虑在轨空间柔性齿轮机构中弹性变形的齿轮啮合的强非线性特点，本发明首先建立齿轮啮合过程中柔性轴和齿轮刚柔耦合的非线性动力学方程。其次，针对在轨空间柔性齿轮机构运动过程中故障率呈现的时变特性，本发明提出了一种时变故障率新算法，可得到柔性齿轮机构故障率时变特性规律。然后应用蒙特卡罗算法验证了该方法的精度和可行性，对于提高在轨空间柔性机构产品可靠性具有一定的工程实用价值。本发明的方法同样适用于其他柔性齿轮机构。

发明内容

本发明的目的是为了提供一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，它分析在轨空间柔性齿轮机构啮合过程中故障率和可靠度随时间变化的规律。考虑到空间柔性齿轮机构齿轮啮合过程中的动态特性以及长转子弹性形变的齿轮啮合问题是一个高度非线性问题，本发明首先采用弹性体振动模型和直齿轮平面振动模型建立了齿轮啮合过程中柔性轴和齿轮刚柔耦合的非线性动力学方程；并采用数值的方法分析齿轮动态刚度激励，误差激励和外部激励，在此基础上，采用Newmark法进行仿真计算，得到齿轮的动载荷系数；其次，针对动载荷系数的时变和强非线性特点，应用kriging方法求解某时刻可靠度指标，避免了传统的传递矩阵法和振型叠加法因非线性因素导致的可靠度计算精度低问题；最后，针对在轨空间柔性齿轮机构运动过程中故障率呈现时变特性，本发明提出一种时变故障率新算法，可得到柔性齿轮机构故障率时变特性规律，对于提高在轨空间柔性机构产品的可靠性具有一定的指导意义和工程实用价值。

本发明一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，该方法具体步骤为：

步骤一：针对在轨空间柔性齿轮机构中轴是弹性体，建立柔性轴的动力学方程，采用有限元方法分析，方法如下：

在齿轮径向力和弯矩的作用下，齿轮轴会产生横向和扭转的复合振动。对于其横向位移和扭转角度的位移形函数可以进行独立的插值，类似于Timoshenko梁单元。

齿轮轴的横向位移函数的边界条件：

因此转角度的位移形函数可以采用线性插值：

θ(x,t)＝a₀+a₁x

代入边界条件得到：

计算梁单元的势能和动能，代入拉格朗日方程，进而得到梁单元的振动方程：

$\left[M_{\mathrm{shaft}}\right]\left\{\stackrel{\·\·}{y}\right\}+\left[K_{\mathrm{shaft}}\right]\left\{\stackrel{\·}{y}\right\}=\left\{F_{\mathrm{shaft}}\right\}$

式中，[M_shaft]——单元质量矩阵，[K_shaft]——单元刚度矩阵，{F_shaft}——单元外力列阵；

步骤二：建立齿轮动力学方程，采用附图1所示的二维平面振动模型系统描述齿轮啮合，y₁，y₂表示与齿轮连接的梁单元的端点位移，得到齿轮接触动力学方程：

$\left[M_{\mathrm{gear}}\right]\left\{\stackrel{\·\·}{y}\right\}+\left[C_{\mathrm{gear}}\right]\left\{\stackrel{\·}{y}\right\}+\left[K_{\mathrm{gear}}\right]\left\{y\right\}=\left\{F_{\mathrm{gear}}\right\}$

式中，[M_gear]——齿轮质量矩阵，[C_gear]——阻尼矩阵，[K_gear]——齿轮刚度矩阵，

{F_gear}——外力列阵；

步骤三：建立在轨空间柔性齿轮机构中柔性轴与齿轮耦合的非线性动力学方程。

将所有单元的振动方程和齿轮的振动方程进行叠加，其中各个单元之间以及齿轮和单元之间的相互作用作为内力而抵消，进而得到系统的动力学方程：

$\left[M\right]\left\{\stackrel{\·\·}{y}\right\}+\left[C\right]\left\{\stackrel{\·}{y}\right\}+\left[K\right]\left\{y\right\}=\left\{F\right\}$

式中，[M]——质量矩阵，[C]——阻尼矩阵，[K]——刚度矩阵，{F}——外力列阵；

步骤四：进行空间柔性齿轮机构传动激励分析。

空间柔性齿轮啮合过程中，存在啮合轮齿对数变化以及齿顶到齿根啮合过程弹性的变形，将会导致啮合综合刚度随时间周期变化，从而引起齿轮轮齿啮合力周期变化。采用有限元仿真的方法计算轮齿的弹性形变，利用ANSYS的APDL语言编程建立齿轮啮合模型如附图2所示，进一步计算得到齿轮啮合的动态刚度；同样，齿轮误差激励也是一种周期性激励，适合采用正弦函数进行描述；另外，空间柔性齿轮机构在制造过程中不可避免得会引入制造误差，从而导致转子偏心的现象，偏心转子会产生动载荷激励。

步骤五：采用Newmark算法计算动载荷系数。

采用Newmark法进行仿真计算，进而得到齿轮的动载荷系数，避免了传统的传递矩阵法和振型叠加法因非线性因素导致的计算精度低问题。

由齿轮动载荷计算式：

F_d＝k_ty_t'

以及动载荷系数计算式：

$K_v=\frac{F_d+F_t}{F_t}$

相结合得到齿轮动载荷系数。

式中，F_d——齿轮动载荷，k_t——齿轮啮合刚度，K_v——动载荷系数，

y‘_t——齿轮啮合线相对位移，F_t——圆周力；

步骤六：在轨空间柔性齿轮机构中齿轮弯曲应力可靠性极限状态函数的建立

可知齿轮弯曲应力计算公式为：

${\mathrm{\σ}}_F-\frac{F_t}{\mathrm{bm}}{Y}_S{Y}_b{k}_V$

式中，Y_S——齿形系数，Y_b——应力矫正系数，b——齿轮宽度，m——齿轮模数；

极限状态函数为：

$G={\mathrm{\σ}}_F-{\mathrm{\σ}}_{\mathrm{Flim}}=\frac{F_t}{\mathrm{bm}}{Y}_S{Y}_b{k}_V-{\mathrm{\σ}}_{\mathrm{Flim}};$

步骤七：空间柔性齿轮机构动态可靠度与故障率分析。

针对动载荷系数的时变和强非线性特点，应用kriging方法求解某一时刻可靠度指标β(t)和法向量α(t)，应用基于时变故障率的新算法分析空间柔性齿轮机构可靠度，最后采用蒙特卡罗方法验证该方法的准确性和有效性。

令X(ω,t)代表机械问题的随机变量，t代表时间历程，ω表示样本空间Ω的一个样本点。因此，可以得到时变极限状态函数G(t,X(ω,t))＝0。那么可靠度函数可以表示为：

R(t)＝prob(G(t,X(ω,t)))≤0

本发明定义时变故障率函数为：

$\mathrm{\λ}\left(t\right)\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}\left(A\right|B)}{\mathrm{\Δt}}$

其中：

A＝{G(t+△t,X(ω,t+△t))≤0}

B＝{G(t,X(ω,t))>0}

根据条件概率公式可得：

$\left(\mathrm{\λ}\right)=\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}\left(A\right|B)}{\mathrm{\Δt}}=\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}(A\∩B)}{\mathrm{\Δtprob}\left(B\right)}$

prob(B)＝Φ(β)

引入二维正态分布函数Φ₂和相关系数ρ_t，可得：

prob(A∩B)＝Φ₂(β(t),-β(t+△t),ρ(t,t+△t))

其中ρ(t,t+△t)可以用极限状态面的法向量α表示：

ρ(t,t+△t)＝-α(t)α(t+△t)

Φ₂的计算可以采用如下的方法：

利用可靠度函数可得到整个时间历程的可靠度与故障率时变特性规律。

本发明方法的优点和积极效果在于：

1)本发明为保证在轨空间柔性齿轮机构的动态高可靠性的研制提供了有效的技术途径。

2)本发明避免了传统的传递矩阵法和振型叠加法因非线性因素导致的可靠度计算精度低问题。

3)本发明给出了一种建立齿轮啮合过程中柔性轴与齿轮耦合的非线性动力学方程的方法。

4)本发明提出了一种时变故障率新算法，分析齿轮啮合过程中故障率与可靠度随时间变化的时变特性规律。

5)本发明得到了故障率和可靠度随时间变化的特性规律。

6)本发明同样适用于极限状态函数高度非线性情况。

7)本发明同样适用于其他柔性齿轮机构。

附图说明

图1为直齿轮平面传动模型简图

图2为齿轮啮合有限元模型图

图3为齿轮啮合动态刚度图

图4为动载荷系数图

图5为啮合齿轮针对弯曲应力时间离散可靠度指标图

图6为时变故障率曲线图

图7为数值积分与蒙特卡洛方法比较图

图8为本发明流程框图

图中符号说明如下：

m₁、m₂——当量质量；I₁、I₂——对应质心转动惯量；T₁、T₁——转矩；θ₁、θ₂——转角；r₁、r₂——当量半径；y₁、y₂——端点位移；K——刚度系数；C_h——阻尼系数；e(τ)——电机转子；σ_F——弯曲应力；F_t——圆周力；b——齿轮宽度；m——齿轮模数；Y_S——齿形系数；Y_b——应力矫正系数；K_v——动载荷系数；σ_Flim——允许弯曲应力极限值。

具体实施方式

见图8，本发明一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，该方法具体步骤为：

步骤一：针对在轨空间柔性齿轮机构中轴是弹性体，建立柔性轴的动力学方程，采用有限元方法分析，方法如下：

在齿轮径向力和弯矩的作用下，齿轮轴会产生横向和扭转的复合振动。对于其横向位移和扭转角度的位移形函数可以进行独立的插值，类似于Timoshenko梁单元。

齿轮轴的横向位移函数的边界条件：

因此转角度的位移形函数可以采用线性插值：

θ(x,t)＝a₀+a₁x

代入边界条件得到：

计算梁单元的势能和动能，代入拉格朗日方程，进而得到梁单元的振动方程：

$\left[M_{\mathrm{shaft}}\right]\left\{\stackrel{\·\·}{y}\right\}+\left[K_{\mathrm{shaft}}\right]\left\{\stackrel{\·}{y}\right\}=\left\{F_{\mathrm{shaft}}\right\}$

式中，[M_shaft]——单元质量矩阵；[K_shaft]——单元刚度矩阵；{F_shaft}——单元外力列阵。

步骤二：建立齿轮动力学方程，采用附图1所示的二维平面振动模型系统描述齿轮啮合，y₁，y₂表示与齿轮连接的梁单元的端点位移，得到齿轮接触动力学方程：

$\left[M_{\mathrm{gear}}\right]\left\{\stackrel{\·\·}{y}\right\}+\left[C_{\mathrm{gear}}\right]\left\{\stackrel{\·}{y}\right\}+\left[K_{\mathrm{gear}}\right]\left\{y\right\}=\left\{F_{\mathrm{gear}}\right\}$

式中，[M_gear]——齿轮质量矩阵；[C_gear]——阻尼矩阵；

[K_gear]——齿轮刚度矩阵；{F_gear}——外力列阵。

步骤三：建立在轨空间柔性齿轮机构中柔性轴与齿轮耦合的非线性动力学方程。

将所有单元的振动方程和齿轮的振动方程进行叠加，其中各个单元之间以及齿轮和单元之间的相互作用作为内力而抵消，进而得到系统的动力学方程：

$\left[M\right]\left\{\stackrel{\·\·}{y}\right\}+\left[C\right]\left\{\stackrel{\·}{y}\right\}+\left[K\right]\left\{y\right\}=\left\{F\right\}$

式中，[M]——质量矩阵；[C]——阻尼矩阵；[K]——刚度矩阵；{F}——外力列阵。

步骤四：进行空间柔性齿轮机构传动激励分析。

空间柔性齿轮啮合过程中，存在啮合轮齿对数变化以及齿顶到齿根啮合过程弹性的变形，将会导致啮合综合刚度随时间周期变化，从而引起齿轮轮齿啮合力周期变化。采用有限元仿真的方法计算轮齿的弹性形变，利用ANSYS的APDL语言编程建立齿轮啮合模型如附图2所示，进一步计算得到齿轮啮合的动态刚度(图3齿轮啮合动态刚度图)；同样，齿轮误差激励也是一种周期性激励，适合采用正弦函数进行描述；另外，空间柔性齿轮机构在制造过程中不可避免得会引入制造误差，从而导致转子偏心的现象，偏心转子会产生动载荷激励。

步骤五：采用Newmark算法计算动载荷系数。

采用Newmark法进行仿真计算，进而得到齿轮的动载荷系数，避免了传统的传递矩阵法和振型叠加法因非线性因素导致的计算精度低问题。图4为动载荷系数图。

由齿轮动载荷计算式：

F_d＝k_ty_t'

以及动载荷系数计算式：

$K_v=\frac{F_d+F_t}{F_t}$

相结合得到齿轮动载荷系数。

式中，F_d——齿轮动载荷；k_t——齿轮啮合刚度；K_v——动载荷系数；

y‘_t——齿轮啮合线相对位移；F_t——圆周力。

步骤六：在轨空间柔性齿轮机构中齿轮弯曲应力可靠性极限状态函数的建立

可知齿轮弯曲应力计算公式为：

${\mathrm{\σ}}_F-\frac{F_t}{\mathrm{bm}}{Y}_S{Y}_b{k}_V$

式中，Y_S——齿形系数；Y_b——应力矫正系数；b——齿轮宽度；m——齿轮模数。

极限状态函数为：

$G={\mathrm{\σ}}_F-{\mathrm{\σ}}_{\mathrm{Flim}}=\frac{F_t}{\mathrm{bm}}{Y}_S{Y}_b{k}_V-{\mathrm{\σ}}_{\mathrm{Flim}};$

步骤七：空间柔性齿轮机构动态可靠度与故障率分析。

针对动载荷系数的时变和强非线性特点，应用kriging方法求解某一时刻可靠度指标β(t)和法向量α(t)，应用基于时变故障率的新算法分析空间柔性齿轮机构可靠度，最后采用蒙特卡罗方法验证该方法的准确性和有效性。图5为啮合齿轮针对弯曲应力时间离散可靠度指标图，图6为时变故障率曲线图，图7为数值积分与蒙特卡罗方法比较图。

令X(ω,t)代表机械问题的随机变量，t代表时间历程，ω表示样本空间？的一个样本点。因此，可以得到时变极限状态函数G(t,X(ω,t))＝0。那么可靠度函数可以表示为：

R(t)＝prob(G(t,X(ω,t)))≤0

本发明定义时变故障率函数为：

$\mathrm{\λ}\left(t\right)\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}\left(A\right|B)}{\mathrm{\Δt}}$

其中：

A＝{G(t+△t,X(ω,t+△t))≤0}

B＝{G(t,X(ω,t))>0}

根据条件概率公式可得：

$\left(\mathrm{\λ}\right)=\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}\left(A\right|B)}{\mathrm{\Δt}}=\underset{\mathrm{\Δt}\→0}{\lim }\frac{\mathrm{prob}(A\∩B)}{\mathrm{\Δtprob}\left(B\right)}$

prob(B)＝Φ(β)

引入二维正态分布函数Φ₂和相关系数ρ_t，可得：

prob(A∩B)＝Φ₂(β(t),-β(t+△t),ρ(t,t+△t))

其中ρ(t,t+△t)可以用极限状态面的法向量α表示：

ρ(t,t+△t)＝-α(t)α(t+△t)

Φ₂的计算可以采用如下的方法：

利用可靠度函数可得到整个时间历程的可靠度与故障率时变特性规律。

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

实施案例

实例描述：

根据在轨空间柔性齿轮机构中齿轮和轴相关参数，运用本发明提出的方法建立齿轮与柔性轴耦合的非线性动力学方程，并对齿轮啮合进行时变可靠性分析，得到故障率与可靠度随时间变化的时变特性，最后运用蒙特卡洛方法验证本发明的精度和可行性。

针对在轨空间柔性齿轮机构弹性体弹性形变的齿轮啮合过程，可以总结出空间柔性齿轮机构时变可靠性分析的关键技术与难点如下：

1)齿轮和柔性轴耦合的非线性动力学方程建立问题

齿轮啮合过程中不可避免的会产生径向力，对于齿轮轴较短的情况，其径向力产生的挠度影响可以忽略。但对于长挠性转子，忽略轴的挠度，就会产生较大的误差。因此，在建立齿轮传动动力学方程时，必须考虑转子的挠度影响。因此，本发明采用弹性体振动模型和直齿轮平面振动模型，建立齿轮和柔性轴耦合的非线性动力学方程。

2)动载荷系数求解的难题

空间柔性齿轮机构中传动系统的刚度矩阵和内外部激励都具有时变和高非线性特性。传统的传递矩阵法和振型叠加法需要忽略式中的非线性因素，会导致计算精度不够。本发明采用Newmark算法计算系统的非线性动力学方程能一定程度上提高计算解的精度。

3)空间柔性齿轮机构时变故障率与可靠性分析问题

传统的采用动载荷系数稳定值进行可靠性分析缺少科学性，而且采用上穿率的方法，只能近似估算失效概率的上下界。本发明针对动载荷系数的时变和强非线性特点，应用kriging方法求解某一时刻可靠度指标β(t)和法向量α(t)，应用基于时变故障率的新算法，分析空间柔性齿轮机构的可靠度，最后采用蒙特卡罗方法验证该方法的准确性与有效性。

第一步，齿轮动载荷系数计算

在轨空间柔性齿轮机构传动系统的刚度矩阵和内外部激励都具有时变特性和非线性。传统的传递矩阵法和振型叠加法需要忽略式中的非线性因素，会导致计算精度不够。本文采用Newmark算法计算系统的动力学方程。运用以上方法编写相应的计算程序可以得到动载荷计算结果。结果如附图4所示。

相关参数

表1齿轮相关参数

第二步，空间柔性齿轮机构时变故障率与可靠性分析

采用一种时变可靠性新算法对齿轮弯曲应力进行时变可靠度与故障率计算，并用蒙特卡罗算法验证了本发明方法的精度和可行性。计算出啮合齿轮针对弯曲应力的时间离散可靠性指标，如附图5所示。

附图6描述故障率随时间变化的规律。

附图7描述了蒙特卡洛方法验证的结果，结果表明本发明方法与蒙特卡洛方法最大误差小于0.02，但本发明计算耗时要明显小于蒙特卡洛算法。

分析结论：

1)本发明得到了故障率和可靠度随时间变化的特性规律。

2)本发明方法计算得到的齿轮啮合针对弯曲应力的可靠度与蒙特卡洛相比误差小于0.02。

3)采用本发明方法得到相同的齿轮啮合中针对弯曲应力的时变可靠度值耗时明显比小于蒙特卡洛方法。

4)本发明方法对于强非线性的弹性体动力学齿轮传动系统时变可靠度计算有明显的效果，具有一定的工程应用价值。

5)本发明为保证空间柔性齿轮机构的动态高可靠性的研制提供了有效的技术途径。

6)本发明同样适用于极限状态函数高度非线性情况。

7)本发明同样适用于其他柔性齿轮机构。

Claims (1)

1.一种在轨空间柔性齿轮机构动态耦合时变故障率分析方法，其特征在于：该方法具体步骤为：

步骤一：针对在轨空间柔性齿轮机构中轴是弹性体，建立柔性轴的动力学方程，采用有限元方法分析，方法如下：

在齿轮径向力和弯矩的作用下，齿轮轴产生横向和扭转的复合振动，对于其横向位移和扭转角度的位移形函数进行独立的插值，为Timoshenko梁单元，

齿轮轴的横向位移函数的边界条件：

因此转角度的位移形函数采用线性插值：

θ(x,t)＝a₀+a₁x

代入边界条件得到：

计算梁单元的势能和动能，代入拉格朗日方程，进而得到梁单元的振动方程：

$\left[M_{\mathrm{shaft}}\right]\left\{\stackrel{\·\·}{y}\right\}+\left[K_{\mathrm{shaft}}\right]\left\{\stackrel{\·}{y}\right\}=\left\{F_{\mathrm{shaft}}\right\}$

式中，[M_shaft]——单元质量矩阵，——单元刚度矩阵，{F_shaft}——单元外力列阵；

步骤二：建立齿轮动力学方程，采用二维平面振动模型系统描述齿轮啮合，y₁，y₂表示与齿轮连接的梁单元的端点位移，得到齿轮接触动力学方程：

$\[M_{gear}\]\left\{\stackrel{\·\·}{y}\right\}+\[C_{gear}\]\left\{\stackrel{\·}{y}\right\}+\[K_{gear}\]\left\{y\right\}=\left\{F_{gear}\right\}$

式中，[M_gear]——齿轮质量矩阵，[C_gear]——阻尼矩阵，[K_gear]——齿轮刚度矩阵，

{F_gear}——外力列阵；

步骤三：建立在轨空间柔性齿轮机构中柔性轴与齿轮耦合的非线性动力学方程；

将所有单元的振动方程和齿轮的振动方程进行叠加，其中各个单元之间以及齿轮和单元之间的相互作用作为内力而抵消，进而得到系统的动力学方程：

$\[M\]\left\{\stackrel{\·\·}{y}\right\}+\[C\]\left\{\stackrel{\·}{y}\right\}+\[K\]\left\{y\right\}=\left\{F\right\}$

式中，[M]——质量矩阵，[C]——阻尼矩阵，[K]——刚度矩阵，{F}——外力列阵；

步骤四：进行空间柔性齿轮机构传动激励分析；

空间柔性齿轮啮合过程中，存在啮合轮齿对数变化以及齿顶到齿根啮合过程弹性的变形，将会导致啮合综合刚度随时间周期变化，从而引起齿轮轮齿啮合力周期变化；采用有限元仿真的方法计算轮齿的弹性形变，利用ANSYS的APDL语言编程建立齿轮啮合模型，进一步计算得到齿轮啮合的动态刚度；同样，齿轮误差激励也是一种周期性激励，适合采用正弦函数进行描述；另外，空间柔性齿轮机构在制造过程中不可避免得会引入制造误差，从而导致转子偏心的现象，偏心转子会产生动载荷激励；

步骤五：采用Newmark算法计算动载荷系数；

采用Newmark法进行仿真计算，进而得到齿轮的动载荷系数，避免了传统的传递矩阵法和振型叠加法因非线性因素导致的计算精度低问题；

由齿轮动载荷计算式：

F_d＝k_ty'_t

以及动载荷系数计算式：

$K_v=\frac{F_d+F_t}{F_t}$

相结合得到齿轮动载荷系数；

式中，F_d——齿轮动载荷，k_t——齿轮啮合刚度，K_v——动载荷系数，

y‘_t——齿轮啮合线相对位移，F_t——圆周力；

步骤六：在轨空间柔性齿轮机构中齿轮弯曲应力可靠性极限状态函数的建立

齿轮弯曲应力计算公式为：

${\mathrm{\σ}}_F=\frac{F_t}{bm}{Y}_S{Y}_b{k}_v$

式中，Y_s——齿形系数，Y_b——应力矫正系数，b——齿轮宽度，m——齿轮模数；

极限状态函数为：

$G={\mathrm{\σ}}_F-{\mathrm{\σ}}_{F\lim }=\frac{F_t}{bm}{Y}_S{Y}_b{k}_v-{\mathrm{\σ}}_{F\lim };$

步骤七：空间柔性齿轮机构动态可靠度与故障率分析；

针对动载荷系数的时变和强非线性特点，应用kriging方法求解某一时刻可靠度指标β(t)和法向量α(t)，应用基于时变故障率的新算法分析空间柔性齿轮机构可靠度，最后采用蒙特卡罗方法验证该方法的准确性和有效性；

令X(ω,t)代表机械问题的随机变量，t代表时间历程，ω表示样本空间Ω的一个样本点，因此，得到时变极限状态函数G(t,X(ω,t))＝0，那么可靠度函数表示为：

R(t)＝prob(G(t,X(ω,t)))≤0

定义时变故障率函数为：

$\mathrm{\λ}\left(t\right)=\underset{\mathrm{\Δ}t\→0}{\lim }\frac{prob\left(A\right|B)}{\mathrm{\Δ}t}$

其中：

A＝{G(t+Δt,X(ω,t+Δt))≤0}

B＝{G(t,X(ω,t))＞0}

根据条件概率公式得：

$\left(\mathrm{\λ}\right)=\underset{\mathrm{\Δ}t\→0}{\lim }\frac{prob\left(A\right|B)}{\mathrm{\Δ}t}=\underset{\mathrm{\Δ}t\→0}{\lim }\frac{prob(A\∩B)}{\mathrm{\Δ}tprob\left(B\right)}$

prob(B)＝Φ(β)

引入二维正态分布函数Φ₂和相关系数ρ_t，得：

prob(A∩B)＝Φ₂(β(t),-β(t+Δt),ρ(t,t+Δt))

其中ρ(t,t+Δt)用极限状态面的法向量α表示：

ρ(t,t+Δt)＝-α(t)α(t+Δt)

Φ₂的计算采用如下的方法：

利用可靠度函数得到整个时间历程的可靠度与故障率时变特性规律。