CN104298857A - 一种多因素耦合作用下的机构可靠度计算方法 - Google Patents

Abstract

本发明一种多因素耦合作用下的机构可靠度计算方法，首先基于多刚体动力学、间隙碰撞模型和柔性体离散化方法对机构进行建模，获得机构输出的数值计算，从而实现对杆件尺寸误差、装配误差、间隙、摩擦、载荷、速度以及变形等多种影响因素的考虑，然后，在机构建模的模型中，对机构输出有影响的上述多种因素进行参数化；最后，基于本发明提出的最小抽样方法，进行机构可靠度高效高精度计算。本发明提出的机构可靠度计算方法，考虑因素更多，因此更加符合实际工程应用。

Description

一种多因素耦合作用下的机构可靠度计算方法

技术领域

本发明涉及一种多因素耦合作用下的机构可靠度计算方法，尤其综合考虑杆件尺寸误差、装配误差、间隙、摩擦、载荷、速度以及变形等多种因素的作用，更加符合工程实际情况，适用于工程机构设计的可靠性评估，可靠性验证等相关领域。

背景技术

在实际工程机构设计中，机构实际功能的实现与理想设计之间存在误差。这是由于机构受到杆件尺寸制造误差、装配误差、间隙、摩擦、载荷、速度等随机因素的影响，同时，机构在载荷作用下，存在变形。综合上述影响，机构的运动输出具有不确定性，即机构输出在一定误差范围内存在概率分布，如图1所示。因此，有必要进行机构可靠性评估，量化其不确定性。这对于航空航天等重要领域的机构安全性评估和应用具有重要意义。

通常，机构运动可靠度可以用公式表示如下：

$p_f=P[g\left(x\right)<0]\&ap;\frac{N_1}{N_2}---\left(1\right)$

R＝1-p_f

向量x＝[x₁,x₂,...,x_n]中的x₁,x₂,...,x_n为各种影响因素，g(x)＝△-ε(x)为机构功能实现的极限状态函数，通常ε(x)为机构输出，如位移、阻力等，△为机构输出的限定值,通常由机构设计目标给定，N₂为向量x＝[x₁,x₂,...,x_n]的抽样样本总数，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，p_f为失效概率，R为可靠度。例如：当ε(x)为机构位移输出误差时，△为给定的机构运动误差精度，此时P[g(x)<0]表示机构位移输出误差小于给定机构运动误差精度的概率；当ε(x)为机构阻力输出时，△为机构驱动力，此时P[g(x)<0]表示机构阻力输出小于驱动力的概率。

目前机构可靠性计算的研究中，存在如下不足：

(1)主要从运动学的角度，着重考虑杆长误差、装配误差、间隙等对机构运动学运动精度影响，并提出相关的研究方法。但是对于机构设计来讲，既有运动学相关影响因素(如杆件尺寸制造误差、装配误差等)，也有动力学相关的影响因素(如间隙、摩擦、载荷、速度、变形等)。而目前文献考虑因素有限。

(2)间隙从动力学和运动学两方面共同对机构输出造成的影响。由于机构可靠性计算本身比较复杂，在对机构输出的建模中，大多采用相对简单的机构建模方法，从而难以从动力学和运动学两方面同时考虑间隙对机构输出的影响。大多数机构可靠性研究方法在对间隙影响的处理时或者不考虑，或者将间隙视为与杆长作用相同的独立随机变量，或者采用等效连接模型等简化方法。上述处理方法与机构间隙碰撞影响机构输出的实际情况存在较大差异。

(3)同样由于机构可靠性计算本身的复杂性，采用的机构建模方法简单，难以考虑变形对机构输出的影响。

因此，探索能够综合考虑杆件尺寸误差、装配误差、间隙、摩擦、载荷、速度以及变形等多种因素影响的机构可靠性计算方法，更符合实际工程应用。这对于工程实际中机构可靠性分析与设计也具有十分重要的意义。

发明内容

本发明的目的在于提供一种多因素耦合作用下的机构可靠度计算方法，在机构可靠度计算中，综合考虑杆件尺寸误差、装配误差、摩擦、载荷、速度以及变形等多种因素的影响，尤其从运动学和动力学两方面考虑间隙对机构输出的影响，有效选择计算抽样样本，减少极限状态函数的计算次数，从而提高机构可靠度的计算效率，有利于实际应用。

一种多因素耦合作用下的机构可靠度计算方法，其中的机构运动可靠度计算公式表示如下：

$p_f=P[g\left(x\right)<0]\&ap;\frac{N_1}{N_2}---\left(1\right)$

R＝1-p_f

向量x＝[x₁,x₂,...,x_n]中的x₁,x₂,...,x_n为各种影响因素，g(x)＝△-ε(x)为机构功能实现的极限状态函数，ε(x)为机构输出，△为机构输出的限定值，由机构设计目标给定，N₂为向量x＝[x₁,x₂,...,x_n]的抽样样本总数，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，p_f为失效概率，R为可靠度，其特征在于包括如下计算步骤：

步骤1、机构建模：

(1)基于多刚体动力学对机构进行建模，并在机构模型中，将对机构输出有影响的因素进行参数化建模，该影响因素包括杆件长度、装配位置、摩擦、载荷和速度；

(2)在机构模型中，引入间隙碰撞模型，建立间隙碰撞的运动学模型、建立间隙碰撞的力学描述、建立碰撞力模型描述；

(3)杆件变形建模：在机构模型中，首先预判对受载变形相对较大的杆件，然后基于柔性体离散化方法对这些杆件重新进行建模，从而实现对上述杆件受载变形的描述；

步骤2、通过步骤1建立完整考虑多因素耦合作用下的机构模型后，对上述机构模型中的杆件长度、装配位置、摩擦、载荷、速度的多个因素视为随机变量，这里假设随机变量的总数为n个，并用随机变量x₁,x₂,...,x_n表示，同时组成随机向量x＝[x₁,x₂,...,x_n]，按照预设的策略获得高效的抽样样本，即在随机向量x各自分量x₁,x₂,...,x_n的分布范围内，获得一组抽样值x^*＝[x₁ ^*,x₂ ^*,...,x_n ^*]，然后作为输入代入步骤1的机构模型，再通过数值计算获得机构输出ε(x^*)及其对应的极限状态函数输出g(x^*)，利用式(1)即可计算失效概率p_f和可靠度R，具体为：

①应用蒙特卡洛方法在随机向量x的各自分量x_i(i＝1～n)的分布范围内，随机抽样N(＝n)个初始样本点形成初始样本集X'＝[x'₁,x'₂,...,x'_N]^T，其中x_j'＝[x_j1',x_j2',...,x_jn'](j＝1～N),然后将N个初始样本点逐一作为输入代入步骤1的机构模型中，获得机构输出ε(x_j')(j＝1～N)及其对应的极限状态函数输出g(x_j')(j＝1～N)，并组成如下矩阵G'＝[g'₁,g'₂,...,g'_N]^T,，

这里将g(x_j')简写为g'_j(j＝1～N)，上述X'和G'如式(2)所示：

$X^{\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;}\\ x_2^{\&prime;}\\ ...\\ x_N^{\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;}& x_{\mathrm{1,2}}^{\&prime;}& ...& x_{1,n}^{\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;}& x_{\mathrm{2,2}}^{\&prime;}& ...& x_{2,n}^{\&prime;}\\ ...& ...& ...& ...\\ x_{N,1}^{\&prime;}& x_{N,2}^{\&prime;}& ...& x_{N,n}^{\&prime;}\end{array}\right],G^{\&prime;}=\left[\begin{array}{c}{g}_1^{\&prime;}\\ g_2^{\&prime;}\\ ...\\ g_N^{\&prime;}\end{array}\right]---\left(2\right)$

②基于Kriging模型构建X'与G'的映射关系，可得：

G'＝f_kri(X')   (3)

③重新再次生成随机向量x的N₂个抽样样本，N₂远大于N，N₂为随机向量x＝[x₁,x₂,...,x_n]的抽样样本总数，如式(4)所示：

$X^{\&prime;\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;\&prime;}\\ x_2^{\&prime;\&prime;}\\ ...\\ x_{N_2}^{\&prime;\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;\&prime;}& x_{\mathrm{1,2}}^{\&prime;\&prime;}& ...& x_{1,n}^{\&prime;\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;\&prime;}& x_{\mathrm{2,2}}^{\&prime;\&prime;}& ...& x_{2,n}^{\&prime;\&prime;}\\ ...& ...& ...& ...\\ x_{N_2,1}^{\&prime;\&prime;}& x_{N_2,2}^{\&prime;\&prime;}& ...& x_{N_2,n}^{\&prime;\&prime;}\end{array}\right]---\left(4\right)$

以Kriging模型f_kri作为代理模型代替步骤1建立的机构模型，将样本X″代入式(3)，可得N₂个G″＝f_kri(X″)，并计算G″<0的个数，即获得N₁，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，最后利用式(1)即可计算失效概率p_f和可靠度R；

④在步骤①已经生成初始样本集X'的前提下，按照预设的策略，利用成熟优化算法求解式(5)，获得新的样本点x_new，具体如下：

$\begin{array}{c}\max {\mathrm{\&sigma;}}_g\left(x\right)*p\left(x\right)*r{\left(x\right)}^n\\ s.t.f_{\mathrm{kri}}^i\left(x\right)=0\\ s.t.x_{\mathrm{down}}\&le;x\&le;x_{\mathrm{up}}\\ r\left(x\right)=\frac{1}{2}*\max \left[\underset{{x^{\&prime;}}_i\&Element;X^{\&prime;}}{\min }\right||x-{x^{\&prime;}}_i|\left|\right]\\ p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}p\left(x_i\right)\\ p\left(x_i\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{xi}}}{e}^{-{(x-u_{\mathrm{xi}})}^2/\left(2{\mathrm{\&sigma;}}_{\mathrm{xi}}^2\right)}\end{array}---\left(5\right)$

其中σ_g(x)为预测在任意随机向量x输入时，对应极限状态函数g(x)输出时的标准差，σ_g(x)可以利用上一次构建的Kriging模型进行预测；x'_i为初始样本集X'中的已知样本，x_down和x_up为随机向量x的上下极限，n为随机向量x中影响因素的个数，u_xi、σ_xi和p(x_i)分别为对应随机变量x_i的均值、标准差和正态分布概率密度函数，p(x)为随机变量x₁,x₂,...,x_n的联合概率密度函数；

⑤将新样本点x_new加入到初始样本集X'中，增加初始样本集X'的样本数，返回步骤②利用式(3)重新构建逼近精度更高的Kriging模型；同时重复步骤③，并比较相邻2次计算的p_f，若||p_f ⁱ-p_f ^i-1||/p_f ^i-1<δ，取δ＝0.1，则失效概率计算结果基本收敛，停止计算，得到机构可靠度R，否则，重复步骤④至步骤⑤。

本发明的一种多因素耦合作用下的机构可靠度计算方法，首先基于多刚体动力学、间隙碰撞模型和柔性体离散化方法对机构进行建模，获得机构输出的数值计算，从而实现对杆件尺寸误差、装配误差、间隙、摩擦、载荷、速度以及变形等多种影响因素的考虑，然后，在机构建模的模型中，对机构输出有影响的上述多种因素进行参数化；最后，基于本发明提出的最小抽样方法，进行机构可靠度高效高精度计算。本发明提出的机构可靠度计算方法，考虑因素更多，因此更加符合实际工程应用。

附图说明

图1为多因素耦合作用下导致机构输出在一定范围内概率分布的示意图；

图2为本发明实施例的影响因素为随机变量的四杆机构；

图3为本发明实施例的刚体i与刚体j之间的运动副间隙接触模型；

图4为本发明实施例的间隙碰撞模型中刚体i和刚体j无侵入碰撞情形；

图5为本发明实施例的间隙碰撞模型中刚体i和刚体j发生一定深度的侵入碰撞情形；

图6为本发明实施例的间隙碰撞模型中刚体i和刚体j之间产生碰撞力的情形；

图7为本发明实施例考虑的各个影响因素及其分布范围；

图8为本发明实施例中影响因素为确定性数值的理想四杆机构计算实例；

图9为本发明实施例四杆机构运动可靠度计算结果比较。

以下结合附图和具体实施例对本发明做进一步详述。

具体实施方式

本发明一种多因素耦合作用下的机构可靠度计算方法，其中的机构运动可靠度计算公式表示如下：

$p_f=P[g\left(x\right)<0]\&ap;\frac{N_1}{N_2}---\left(1\right)$

R＝1-p_f

向量x＝[x₁,x₂,...,x_n]中的x₁,x₂,...,x_n为各种影响因素，g(x)＝△-ε(x)为机构功能实现的极限状态函数，ε(x)为机构输出，△为机构输出的限定值，由机构设计目标给定，N₂为向量x＝[x₁,x₂,...,x_n]的抽样样本总数，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，p_f为失效概率，R为可靠度，包括如下计算步骤：

步骤1、机构建模：

(1)基于多刚体动力学对机构进行建模，并在机构模型中，将对机构输出有影响的因素进行参数化建模，该影响因素包括杆件长度、装配位置、摩擦、载荷和速度，以便在后续的机构可靠度计算中考虑上述因素对机构输出的影响；

(2)在机构模型中，引入间隙碰撞模型，建立间隙碰撞的运动学模型、建立间隙碰撞的力学描述、建立碰撞力模型描述，以便精确描述机构运动中，间隙碰撞的过程，进而从运动学和动力学两方面考虑间隙对机构输出造成的影响；

(3)杆件变形建模：在机构模型中，首先预判对受载变形相对较大的杆件，然后基于柔性体离散化方法对这些杆件重新进行建模，从而实现对上述杆件受载变形的描述；

步骤2、通过步骤1建立完整考虑多因素耦合作用下的机构模型后，对上述机构模型中的杆件长度、装配位置、摩擦、载荷、速度的多个因素视为随机变量，这里假设随机变量的总数为n个，并用随机变量x₁,x₂,...,x_n表示，同时组成随机向量x＝[x₁,x₂,...,x_n]，然后将随机向量x的某一样本抽样值，代入步骤1的机构模型中，通过数值计算，获得该样本对应的机构输出ε(x)及其对应的极限状态函数输出g(x)。但是上述随机抽样效率低下，导致循环抽样计算次数多，计算量大。这里提出最小抽样计算方法，即按照预设的策略获得高效的抽样样本，即在随机向量x各自分量x₁,x₂,...,x_n的分布范围内，获得一组抽样值x^*＝[x₁ ^*,x₂ ^*,...,x_n ^*]，然后作为输入代入步骤1的机构模型，再通过数值计算获得机构输出ε(x^*)及其对应的极限状态函数输出g(x^*)，从而提高机构可靠度计算效率，并获得高精度的机构可靠度计算结果：

①应用蒙特卡洛方法在随机向量x的各自分量x_i(i＝1～n)的分布范围内，随机抽样N(＝n)个初始样本点形成初始样本集X'＝[x'₁,x'₂,...,x'_N]^T，其中x_j'＝[x_j1',x_j2',...,x_jn'](j＝1～N),然后将N个初始样本点逐一作为输入代入步骤1的机构模型中，获得机构输出ε(x_j')(j＝1～N)及其对应的极限状态函数输出g(x_j')(j＝1～N)，并组成如下矩阵G'＝[g'₁,g'₂,...,g'_N]^T(这里将g(x_j')简写为g'_j(j＝1～N))，上述X'和G'如式(2)所示：

$X^{\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;}\\ x_2^{\&prime;}\\ ...\\ x_N^{\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;}& x_{\mathrm{1,2}}^{\&prime;}& ...& x_{1,n}^{\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;}& x_{\mathrm{2,2}}^{\&prime;}& ...& x_{2,n}^{\&prime;}\\ ...& ...& ...& ...\\ x_{N,1}^{\&prime;}& x_{N,2}^{\&prime;}& ...& x_{N,n}^{\&prime;}\end{array}\right],G^{\&prime;}=\left[\begin{array}{c}{g}_1^{\&prime;}\\ g_2^{\&prime;}\\ ...\\ g_N^{\&prime;}\end{array}\right]---\left(2\right)$

②基于Kriging 模型，构建X'与G'的映射关系，可得：

G'＝f_kri(X')   (3)

现有Kriging 模型建模方法比较成熟，也有现成的计算代码，具体使用方法不再熬述。

③重新再次生成随机向量x的N₂个抽样样本，N₂远大于N，N₂为随机向量x＝[x₁,x₂,...,x_n]的抽样样本总数，如式(4)所示：

$X^{\&prime;\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;\&prime;}\\ x_2^{\&prime;\&prime;}\\ ...\\ x_{N_2}^{\&prime;\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;\&prime;}& x_{\mathrm{1,2}}^{\&prime;\&prime;}& ...& x_{1,n}^{\&prime;\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;\&prime;}& x_{\mathrm{2,2}}^{\&prime;\&prime;}& ...& x_{2,n}^{\&prime;\&prime;}\\ ...& ...& ...& ...\\ x_{N_2,1}^{\&prime;\&prime;}& x_{N_2,2}^{\&prime;\&prime;}& ...& x_{N_2,n}^{\&prime;\&prime;}\end{array}\right]---\left(4\right)$

以Kriging模型f_kri作为代理模型代替步骤1建立的机构模型，将样本X″代入式(3)，可得N₂个G″＝f_kri(X″)，并计算G″<0的个数，即获得N₁，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，最后利用式(1)即可计算失效概率p_f和可靠度R；该步骤中计算N₂个G″＝f_kri(X″)时，由于以Kriging模型f_kri作为代理模型，不必数值求解步骤1建立的机构模型，因而计算速度可以大大节省。

但是由此计算失效概率p_f和可靠度R的精度取决于Kriging代理模型f_kri与机构模型的逼近程度，由于构建Kriging代理模型f_kri时样本数较少(为N＝n个)，需要进一步按照预设的策略，选择新的抽样样本，进一步对构建的Kriging代理模型f_kri进行修正，从而提高失效概率p_f和可靠度R的精度。

④在步骤①已经生成初始样本集X'的前提下，按照预设的策略，利用成熟优化算法求解式(5)，获得新的样本点x_new，具体如下：

$\begin{array}{c}\max {\mathrm{\&sigma;}}_g\left(x\right)*p\left(x\right)*r{\left(x\right)}^n\\ s.t.f_{\mathrm{kri}}^i\left(x\right)=0\\ s.t.x_{\mathrm{down}}\&le;x\&le;x_{\mathrm{up}}\\ r\left(x\right)=\frac{1}{2}*\max \left[\underset{{x^{\&prime;}}_i\&Element;X^{\&prime;}}{\min }\right||x-{x^{\&prime;}}_i|\left|\right]\\ p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}p\left(x_i\right)\\ p\left(x_i\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{xi}}}{e}^{-{(x-u_{\mathrm{xi}})}^2/\left(2{\mathrm{\&sigma;}}_{\mathrm{xi}}^2\right)}\end{array}---\left(5\right)$

其中σ_g(x)为预测在任意随机向量x输入时，对应极限状态函数g(x)输出时的标准差，σ_g(x)可以利用上一次构建的Kriging模型进行预测；x'_i为初始样本集X'中的已知样本，x_down和x_up为随机向量x的上下极限，n为随机向量x中影响因素的个数，u_xi、σ_xi和p(x_i)分别为对应随机变量x_i的均值、标准差和正态分布概率密度函数，p(x)为随机变量x₁,x₂,...,x_n的联合概率密度函数；

按照式(5)获得的新样本点x_new具有如下特点：

(a)新样本点x_new基本处于极限状态曲面上，即

(b)新样本点x_new位于已知样本点(即初始样本集X'内的样本点)之间较稀疏的区域。

(c)新样本点x_new的预测响应值的误差较大，即标准差σ_g(x_new)较大。

(d)新样本点x_new的出现概率较高，即p(x)大。

将具有上述特点的新样本点x_new加入初始样本集X'，有利于减少后续新增样本点的数目，同时最大化提高Kriging模型的逼近精度。

⑤将新样本点x_new加入到初始样本集X'中，增加初始样本集X'的样本数，返回步骤②利用式(3)重新构建Kriging模型f_kri的逼近精度将进一步提高，从而提高失效概率p_f和可靠度R的精度。因此，新样本点x_new选择得合理，有利于减少后续新增样本点的数目，从而可以减少计算机数值求解步骤1所建立机构模型的次数，从而大大减少计算时间。重复步骤③，并比较相邻2次计算的p_f，若||p_f ⁱ-p_f ^i-1||/p_f ^i-1<δ，取δ＝0.1，则失效概率计算结果基本收敛，停止计算，得到机构可靠度R，否则，重复步骤④至步骤⑤。

以图2所示的四杆机构为例，计算该机构中P点的运动误差可靠度，也就是P点运动误差ε(x)小于给定运动精度△＝1mm的概率，即可靠度值R＝1-P[g(x)＝△-ε(x)<0]，具体步骤如下：

步骤(1)、机构建模：

1)对于任意机构，可以应用通用多刚体动力学模型对机构进行建模描述，该多刚体动力学模型如下：

$\left\{\begin{array}{c}M\stackrel{..}{q}+{{\mathrm{\&Phi;}}_q}^T\mathrm{\&lambda;}=Q\\ {\mathrm{\&Phi;}}_q\stackrel{..}{q}=\mathrm{\&gamma;}-2\mathrm{\&alpha;}{\mathrm{\&Phi;}}_q-{\mathrm{\&beta;}}^2{\mathrm{\&Phi;}}_q\\ \mathrm{\&gamma;}=-{\left({\mathrm{\&Phi;}}_q\stackrel{.}{q}\right)}_q\stackrel{.}{q}-{\mathrm{\&Phi;}}_{\mathrm{tt}}-2{\mathrm{\&Phi;}}_{\mathrm{qt}}\stackrel{.}{q}\end{array}\right.---\left(6\right)$

其中，代表系统广义坐标，M＝diag(M₁,M₂,...,M_n)代表系统质量矩阵，Q＝[Q₁ ^T,Q₂ ^T,...,Q_n ^T]^T代表系统广义力矩阵，Φ＝[Φ₁ ^T,Φ₂ ^T,...,Φ_m ^T]^T代表系统约束方程，λ代表拉格朗日乘数向量，n表示刚体的数目，M代表约束方程的数目，α、β为正常数，分别代表速度和位置约束的回馈控制参数，γ代表加速度约束方程，t代表矩阵的转置；

2)间隙碰撞模型建模：如图2所示，机构运转中，始终存在间隙碰撞运动，为全面描述实际间隙接触碰撞过程，这里采用图4所示的间隙接触模型，图3使用两个平面圆(由销轴圆，轴套套孔圆表示)来代表轴套套孔和销轴，分别附加到刚体i和刚体j：

①建立间隙碰撞的运动学模型：式(7)为刚体i和刚体j无侵入碰撞情形下(如图4所示)的间隙运动描述：

$\begin{array}{c}{\stackrel{\&RightArrow;}{r}}_k^P={\stackrel{\&RightArrow;}{r}}_k+{\stackrel{\&RightArrow;}{s}}_k^P(k=i,j)\\ {\stackrel{\&RightArrow;}{e}}_{\mathrm{ij}}={\stackrel{\&RightArrow;}{r}}_j^P-{\stackrel{\&RightArrow;}{r}}_i^P\\ e_{\mathrm{ij}}=\sqrt{{\stackrel{\&RightArrow;}{e}}_{\mathrm{ij}}^T{\stackrel{\&RightArrow;}{e}}_{\mathrm{ij}}}\\ \stackrel{\&RightArrow;}{n}={\stackrel{\&RightArrow;}{e}}_{\mathrm{ij}}/e_{\mathrm{ij}}\end{array}---\left(7\right)$

式(3)中，当k＝i时,表示刚体i的质心O_i在全局坐标系X-Y下的位置向量,表示与刚体i固联的轴套中心点P_i到刚体i的质心O_i的位置向量,为在局部坐标系ξ_i—η_i下与刚体i固联的轴套和刚体i质心O_i之间的位置矢量；当k＝j时,表示刚体j的质心O_j在全局坐标系X-Y下的位置向量,表示与刚体j固联的销轴中心点P_j到刚体j的质心O_j的位置向量,为在局部坐标系ξ_j—η_j下与刚体j固联的销轴和刚体j质心O_j之间的位置矢量；为轴套和销轴的偏心向量；偏心距e_ij为偏心向量的模；为偏心向量的单位方向矢量；

式(8)为刚体i和刚体j发生一定深度的侵入碰撞情形下(如图5所示)的间隙运动描述：

$\begin{array}{c}c=R_i-R_j\\ \mathrm{\&delta;}=e_{\mathrm{ij}}-c\\ {\stackrel{\&RightArrow;}{r}}_k^Q={\stackrel{\&RightArrow;}{r}}_k+{\stackrel{\&RightArrow;}{s}}_i^Q+R_k\stackrel{\&RightArrow;}{n}(k=i,j)\\ {\stackrel{\stackrel{\&CenterDot;}{\&RightArrow;}}{r}}_k^Q={\stackrel{\stackrel{\&CenterDot;}{\&RightArrow;}}{r}}_k+{\stackrel{\&RightArrow;}{s}}_k^Q+R_k\stackrel{\stackrel{\&CenterDot;}{\&RightArrow;}}{n}(k=i,j)\\ {\mathrm{\&nu;}}_n={({\stackrel{\stackrel{\&CenterDot;}{\&RightArrow;}}{r}}_j^Q-{\stackrel{\stackrel{\&CenterDot;}{\&RightArrow;}}{r}}_i^Q)}^T\stackrel{\&RightArrow;}{n}\end{array}---\left(8\right)$

式(8)中，R_i和R_j为轴套和销轴的半径；c为轴套和销轴的间隙；δ表示刚体i和刚体j在碰撞时的侵入深度；当k＝i时,表示刚体i上的碰撞点Q_i到刚体i的质心O_i的位置向量,为刚体i碰撞点Q_i在全局坐标系X-Y下的位置向量，为碰撞点Q_i在全局坐标系X-Y下的速度；当k＝j时,表示刚体j上的碰撞点Q_j到刚体j的质心O_j的位置向量,为刚体j碰撞点Q_j在全局坐标系X-Y下的位置向量，为碰撞点Q_j在全局坐标系X-Y下的速度；ν_N和ν_T分别为碰撞点Q_i和Q_j之间相对运动速度的法向分量和切向分量；单位向量由单位向量逆时针旋转90°获得；

②建立间隙碰撞的力学描述：式(9)对应图6间隙碰撞的力学描述：

$\begin{array}{c}{f}_T=\mathrm{\&mu;}{f}_N\\ f_i=f_N+f_T\\ m_i=-(y_i^Q-y_i){f_i}^x+(x_i^Q-x_i){f_i}^y\\ f_j=-f_i\\ m_j=-(y_j^Q-y_j)f_j^x+(x_j^Q-x_j)f_j^y\end{array}---\left(9\right)$

式(9)中f_N和f_T分别为间隙碰撞中的法向力和切向力；μ为摩擦系数；f_i和m_i为作用在刚体i上的力和力矩；f_j和m_j为作用在刚体j上的力和力矩。在计算上述力和力矩时，将式(9)表示的力和力矩一起叠加到式(6)的广义力向量Q中即可；

③建立碰撞力模型描述：

$\begin{array}{c}{F}_N=K{\mathrm{\&delta;}}^{1.5}(1+\frac{3(1-e^2)}{4}\frac{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}{{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}^{(-)}})\\ K=\frac{4}{3\mathrm{\&pi;}(h_i+h_j)}{\left(\frac{R_i{R}_j}{R_i-R_j}\right)}^{\frac{1}{2}}\\ h_k=(1-v_k^2)/\left(\mathrm{\&pi;}{E}_k\right),k=i,j\end{array}---\left(10\right)$

式(10)中：δ是侵入深度，是相对侵入速度，是冲击碰撞速度e是材料的弹性恢复系数，ν_k和E_k分别为泊松系数和杨氏模量，K表示刚度，F_N表示碰撞力，E_k和ν_k为材料的杨氏模量和泊松比，π＝3.14；h_k是中间代替量，无任何意义，这样式(10)中K的表示式比较简洁，即没有出现E_k，ν_k和π。

3)基于柔性体离散化方法进行杆件变形建模：对于机构中变形较大的杆件，可以采用柔性体离散化方法近似描述杆件的变形，例如对于图2中的杆件2将其转化成多个集中质量单元串联而成，每个质量单元可看作一个刚体，同片相邻的两个集中质量单元之间用无质量的Timoshenko梁连接起来，作为承载元件，式(11)为具体描述：

$\mathrm{\&omega;}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{N}_i{\mathrm{\&omega;}}_i,\mathrm{\&theta;}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{N}_i{\mathrm{\&theta;}}_i$

$K=\underset{e}{\mathrm{\&Sigma;}}{K}^e,a=\underset{e}{\mathrm{\&Sigma;}}{a}^e,P=\underset{e}{\mathrm{\&Sigma;}}{P}^e$

$K^e=K_b^e+K_s^e$

$K_b^e=\frac{\mathrm{EIl}}{2}{\&Integral;}_{-1}^1{B}_b^T{B}_b\mathrm{d\&xi;},K_s^e=\frac{\mathrm{GAl}}{2k}{\&Integral;}_{-1}^1{B}_s^T{B}_s\mathrm{d\&xi;}$

B_b＝[B_b1  B_b2...  B_bn],B_s＝[B_s1  B_s2...  B_sn]

$B_{\mathrm{bi}}=\left[\begin{array}{cc}0& -\frac{{\mathrm{dN}}_i}{\mathrm{dx}}\end{array}\right],B_{\mathrm{si}}=\left[\begin{array}{cc}\frac{{\mathrm{dN}}_i}{\mathrm{dx}}& -N_i\end{array}\right](i=\mathrm{1,2},...,n)$

$P^e=\frac{l}{2}{\&Integral;}_{-1}^1{N}^T\left[\begin{array}{c}q\\ 0\end{array}\right]\mathrm{d\&xi;}+\underset{j}{\mathrm{\&Sigma;}}{N}^T\left({\mathrm{\&xi;}}_j\right)\left[\begin{array}{c}{p}_j\\ 0\end{array}\right]-\underset{m}{\mathrm{\&Sigma;}}N\left({\mathrm{\&xi;}}_m\right)\left[\begin{array}{c}0\\ M_m\end{array}\right]$

N＝[N₁  N₂...  N_n]

$N_i=\left[\begin{array}{cc}{N}_i& \\ & N_i\end{array}\right](i=\mathrm{1,2},...,n)$

$a^e={\left[\begin{array}{cccc}{a}_1^T& a_2^T& ...& a_n^T\end{array}\right]}^T$                                                          (11)

$a_i=\left[\begin{array}{c}{\mathrm{\&omega;}}_i\\ {\mathrm{\&theta;}}_i\end{array}\right](i=\mathrm{1,2},...,n)$

式(11)中，ω和θ为多个Timoshenko梁总的扰度和截面转角，N_i为Hermite插值函数，ω_i为节点i的扰度，θ_i为节点i的转角，n为总节点数K为总刚度矩阵，K^e为单元刚度矩阵，由单元弯曲刚度阵和单元剪切刚度矩阵组成，a为位移矩阵，a^e为单元位移，P为总力矩阵，P^e为单元力矩阵，B_b为总弯曲形函数矩阵，B_bi为节点i的弯曲形函数，B_s为总剪切形函数矩阵，B_si为节点i的剪切形函数，E为弹性模量，I为截面转动惯量，l为梁的长度，G为剪切弹性模量，A为梁的截面面积，k为校正因子，ξ是单元内的自然坐标，x表示形函数沿梁的长度方向，q为均匀分布力，p_j为节点j的集中载荷，M_m为节点m的力矩,N为由Hermite插值函数N_i组成的矩阵，a_i为节点i的位移。

步骤(2)、通过步骤(1)建立完整考虑多因素耦合作用下的机构模型后，即针对随机向量x的某一样本抽样值x'，将其代入机构模型中，通过数值计算，获得机构输出ε(x')及其对应的极限状态函数输出g(x')；

计算图8所示的理想四杆机构的运动输出：图8所示理想四杆机构中在起始运动角度变化到300°过程中，P点随x轴方向和y轴方向的位移输出和这里所谓的理想四杆机构指的是不考虑变形的影响；杆件尺寸没有误差(即l₁～l₅为确定性值)；装配位置没有误差(即x_A,y_A，α,β为确定性值)；不考虑间隙的影响(即C₁～C₃均为0)；摩擦系数没有变化(即u₁～u₃为确定值)；载荷F和速度V为确定性值。上述机构参数l₁～l₅，x_A,y_A，α,β以及为图8所示的设计值。四杆机构在机构参数为上述确定性输入的情况下，其水平运动输出和垂直运动输出也为确定性输出。这里采用成熟的多刚体动力学方法，即式(6)，对四杆机构进行建模，即可获得理想四杆机构的运动输出和

步骤(3)、为克服蒙特卡洛抽样效率低下导致循环抽样计算次数多，计算量大的缺点，提出最小抽样计算方法，即按照预设策略获得高效的抽样样本，该抽样样本为在向量值x各自分量x₁,x₂,...,x_n的分布范围内，获得一组抽样值x^*＝[x₁ ^*,x₂ ^*,...,x_n ^*]，从而提高机构可靠度计算效率，并获得高精度的机构可靠度计算结果。

图7为图2四杆机构可靠性计算中考虑的影响因素。这里将图7所示的18个影响因素视为工程中最为常用的正态分布随机变量，同时将其名义值当做各个影响因素的均值，并且其对应的标准差按照“6σ”原则确定，如图7所示。现进一步介绍综合考虑上述影响因素的四杆机构可靠度计算方法。包括如下：

①将上述18个影响因素用随机变量x₁,x₂,...,x₁₈表示，并组成随机矩阵向量x＝[x₁,x₂,...,x₁₈]。然后应用蒙特卡洛方法在随机矩阵向量x的各自分量的分布范围内，随机抽样N＝18个初始样本集X'＝[x'₁,x'₂,...,x'₁₈]^T(如式(12)所示)，并将N(＝18)个样本点逐一作为输入代入四杆机构模型中，从而获得四杆机构中P点的水平运动输出和垂直运动输出这里i＝1～N(＝18)。

$X^{\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;}\\ x_2^{\&prime;}\\ ...\\ x_{18}^{\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;}& x_{\mathrm{1,2}}^{\&prime;}& ...& x_{\mathrm{1,18}}^{\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;}& x_{\mathrm{2,2}}^{\&prime;}& ...& x_{\mathrm{2,18}}^{\&prime;}\\ ...& ...& ...& ...\\ x_{\mathrm{18,1}}^{\&prime;}& x_{\mathrm{18,2}}^{\&prime;}& ...& x_{\mathrm{18,18}}^{\&prime;}\end{array}\right]---\left(12\right)$

进一步计算考虑多因素影响下的四杆机构运动输出与理想四杆机构运动输出的误差，即

进一步可以得到对应的极限状态函数G'＝[g'₁,g'₂,...,g'_N]^T，其中g'_j＝g(x'_j)＝△-ε(x'_j)。该步计算样本X'的输出(即g'₁,g'₂,...,g'_N)，共计N个。因此，需要计算机N次数值计算求解步骤(6)～(11)所建立的多因素影响的四杆机构模型。通常，工程实际中，机构模型较为复杂，其数值计算求解耗时较长。

②基于Kriging模型，构建X'与G'的映射关系，如式(3)所示。现有Kriging模型建模方法比较成熟，也有现成的计算代码，具体使用方法不再熬述。

③生成随机矩阵向量x＝[x₁,x₂,...,x₁₈]的N₂(假设N₂＝1000)个抽样样本。以Kriging模型f_kri作为代理模型代替步骤(6)～(11)建立的多因素影响的四杆机构模型，将样本X″带入式(3)，可得N₂个G″＝f_kri(X″)。并计算G″<0的个数(即N₁)，最后利用式(1)即可计算失效概率p_f和可靠度R。该步计算N₂个G″＝f_kri(X″)时，由于以Kriging模型f_kri作为代理模型，不必数值求解步骤(6)～(11)建立的多因素影响的四杆机构模型，因而计算速度可以大大节省。

但是由此计算失效概率p_f和可靠度R的精度取决于Kriging代理模型f_kri与多因素影响的四杆机构模型的逼近程度。由于构建Kriging代理模型f_kri时样本数较少(为N＝n个)，需要按照一定策略，选择新的抽样样本，进一步对构建的Kriging代理模型f_kri进行修正，从而提高失效概率p_f和可靠度R的精度。

④在前面已经生成的样本集X'的前提下，按照一定策略合理选择新的样本，如式(5)所示。利用成熟优化算法求解式(5)，可以获得新的样本点x_new。

⑤将新样本点x_new加入到样本集X'中，增加样本集X'的样本数，由此利用式(3)构建Kriging模型f_kri的逼近精度将进一步提高，从而提高失效概率p_f和可靠度R的精度。因此，新样本选择得合理，有利于减少后续新增样本点的数目，从而可以减少计算机数值求解步骤(6)～(11)所建立的多因素影响的四杆机构模型的次数，从而大大减少计算时间。同时重复②和③，并比较相邻2次计算的p_f。若||p_f ⁱ-p_f ^i-1||/p_f ^i-1<δ(这里取δ＝0.1)，则失效概率计算结果基本收敛，最后停止计算。否则，重复步骤④。

⑥如图9所示，按照本发明计算本实例的失效概率为0.0015，与通常作为衡量标准的蒙特卡洛计算结果相当。但是本发明共进行了158次极限状态函数计算，而蒙特卡洛法抽样次数为2000，极限状态函数计算次数也为2000次。由此可见本发明在保证可靠度计算结果的同时，可以大大降低多因素影响的机构模型数值求解的计算次数，计算效率极为高效。

本发明的重点在于：(1)将间隙碰撞模型和柔性体离散化方法整合到多刚体动力学系统模型中，实现了在机构可靠度计算中，应用间隙碰撞模型考虑间隙对机构运动的影响，应用柔性体离散化方法考虑变形的影响，并实现在多刚体动力学系统模型中综合考虑杆件尺寸误差、装配误差、摩擦、载荷、速度的影响，从而更加符合工程实际情况。(2)由于考虑的因素较多，多因素影响的机构模型复杂，需要花费大量时间数值求解，为减少机构可靠度计算中，数值求解极限状态函数的次数，提出了高效获取新样本的方法。该新样本的获取不是随机选择的，而是基于之前样本的信息(包括样本的位置分布和样本的响应情况)，按照一定策略获取，从而提高可靠度计算效率和结果。该方法对于工程实际应用具有重要意义。

以上所述，仅是本发明较佳实施例而已，并非对本发明的技术范围作任何限制，故凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (1)

1.一种多因素耦合作用下的机构可靠度计算方法，其中的机构运动可靠度计算公式表示如下：

$p_f=P[g\left(x\right)<0]\&ap;\frac{N_1}{N_2}---\left(1\right)$

R＝1-p_f

向量x＝[x₁,x₂,...,x_n]中的x₁,x₂,...,x_n为各种影响因素，g(x)＝△-ε(x)为机构功能实现的极限状态函数，ε(x)为机构输出，△为机构输出的限定值，由机构设计目标给定，N₂为向量x＝[x₁,x₂,...,x_n]的抽样样本总数，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，p_f为失效概率，R为可靠度，其特征在于包括如下计算步骤：

步骤1、机构建模：

(1)基于多刚体动力学对机构进行建模，并在机构模型中，将对机构输出有影响的因素进行参数化建模，该影响因素包括杆件长度、装配位置、摩擦、载荷和速度；

(2)在机构模型中，引入间隙碰撞模型，建立间隙碰撞的运动学模型、建立间隙碰撞的力学描述、建立碰撞力模型描述；

(3)杆件变形建模：在机构模型中，首先预判对受载变形相对较大的杆件，然后基于柔性体离散化方法对这些杆件重新进行建模，从而实现对上述杆件受载变形的描述；

步骤2、通过步骤1建立完整考虑多因素耦合作用下的机构模型后，对上述机构模型中的杆件长度、装配位置、摩擦、载荷、速度的多个因素视为随机变量，这里假设随机变量的总数为n个，并用随机变量x₁,x₂,...,x_n表示，同时组成随机向量x＝[x₁,x₂,...,x_n]，按照预设的策略获得高效的抽样样本，即在随机向量x各自分量x₁,x₂,...,x_n的分布范围内，获得一组抽样值x^*＝[x₁ ^*,x₂ ^*,...,x_n ^*]，然后作为输入代入步骤1的机构模型，再通过数值计算获得机构输出ε(x^*)及其对应的极限状态函数输出g(x^*)，利用式(1)即可计算失效概率p_f和可靠度R，具体为：

①应用蒙特卡洛方法在随机向量x的各自分量x_i(i＝1～n)的分布范围内，随机抽样N(＝n)个初始样本点形成初始样本集X'＝[x'₁,x'₂,...,x'_N]^T，其中x_j'＝[x_j1',x_j2',...,x_jn'](j＝1～N),然后将N个初始样本点逐一作为输入代入步骤1的机构模型中，获得机构输出ε(x_j')(j＝1～N)及其对应的极限状态函数输出g(x_j')(j＝1～N)，并组成如下矩阵G'＝[g'₁,g'₂,...,g'_N]^T,，

这里将g(x_j')简写为g'_j(j＝1～N)，上述X'和G'如式(2)所示：

$X^{\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;}\\ x_2^{\&prime;}\\ ...\\ x_N^{\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;}& x_{\mathrm{1,2}}^{\&prime;}& ...& x_{1,n}^{\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;}& x_{\mathrm{2,2}}^{\&prime;}& ...& x_{2,n}^{\&prime;}\\ ...& ...& ...& ...\\ x_{N,1}^{\&prime;}& x_{N,2}^{\&prime;}& ...& x_{N,n}^{\&prime;}\end{array}\right],G^{\&prime;}=\left[\begin{array}{c}{g}_1^{\&prime;}\\ g_2^{\&prime;}\\ ...\\ g_N^{\&prime;}\end{array}\right]---\left(2\right)$

②基于Kriging模型构建X'与G'的映射关系，可得：

G'＝f_kri(X')   (3)

③重新再次生成随机向量x的N₂个抽样样本，N₂远大于N，N₂为随机向量x＝[x₁,x₂,...,x_n]的抽样样本总数，如式(4)所示：

$X^{\&prime;\&prime;}=\left[\begin{array}{c}{x}_1^{\&prime;\&prime;}\\ x_2^{\&prime;\&prime;}\\ ...\\ x_{N_2}^{\&prime;\&prime;}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}^{\&prime;\&prime;}& x_{\mathrm{1,2}}^{\&prime;\&prime;}& ...& x_{1,n}^{\&prime;\&prime;}\\ x_{\mathrm{2,1}}^{\&prime;\&prime;}& x_{\mathrm{2,2}}^{\&prime;\&prime;}& ...& x_{2,n}^{\&prime;\&prime;}\\ ...& ...& ...& ...\\ x_{N_2,1}^{\&prime;\&prime;}& x_{N_2,2}^{\&prime;\&prime;}& ...& x_{N_2,n}^{\&prime;\&prime;}\end{array}\right]---\left(4\right)$

以Kriging模型f_kri作为代理模型代替步骤1建立的机构模型，将样本X″代入式(3)，可得N₂个G″＝f_kri(X″)，并计算G″<0的个数，即获得N₁，N₁为向量x＝[x₁,x₂,...,x_n]的抽样样本中，g(x)<0的数目，最后利用式(1)即可计算失效概率p_f和可靠度R；

④在步骤①已经生成初始样本集X'的前提下，按照预设的策略，利用成熟优化算法求解式(5)，获得新的样本点x_new，具体如下：

$\begin{array}{c}\max {\mathrm{\&sigma;}}_g\left(x\right)*p\left(x\right)*r{\left(x\right)}^n\\ s.t.f_{\mathrm{kri}}^i\left(x\right)=0\\ s.t.x_{\mathrm{down}}\&le;x\&le;x_{\mathrm{up}}\\ r\left(x\right)=\frac{1}{2}*\max \left[\underset{{x^{\&prime;}}_i\&Element;X^{\&prime;}}{\min }\right||x-{x^{\&prime;}}_i|\left|\right]\\ p\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}p\left(x_i\right)\\ p\left(x_i\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{xi}}}{e}^{-{(x-u_{\mathrm{xi}})}^2/\left(2{\mathrm{\&sigma;}}_{\mathrm{xi}}^2\right)}\end{array}---\left(5\right)$

其中σ_g(x)为预测在任意随机向量x输入时，对应极限状态函数g(x)输出时的标准差，σ_g(x)可以利用上一次构建的Kriging模型进行预测；x'_i为初始样本集X'中的已知样本，x_down和x_up为随机向量x的上下极限，n为随机向量x中影响因素的个数，u_xi、σ_xi和p(x_i)分别为对应随机变量x_i的均值、标准差和正态分布概率密度函数，p(x)为随机变量x₁,x₂,...,x_n的联合概率密度函数；

⑤将新样本点x_new加入到初始样本集X'中，增加初始样本集X'的样本数，返回步骤②利用式(3)重新构建逼近精度更高的Kriging模型；同时重复步骤③，并比较相邻2次计算的p_f，若||p_f ⁱ-p_f ^i-1||/p_f ^i-1<δ，取δ＝0.1，则失效概率计算结果基本收敛，停止计算，得到机构可靠度R，否则，重复步骤④至步骤⑤。