ABSTRACT

This paper presents a new approach to fitting fatigue datasets to three-parameter Weibull distribution by adjusting location parameter t_o and maximum accumulated fatigue life rate 99% in order to match the intrinsic higher moment statistics such as skewness and kurtosis from sample fatigue data directly. Thereafter, the parameters κ and λ of the Weibull distribution from various manufacture technology can be well estimated. The expected standard deviation σ, mean μ and accumulated failure life L10 can be calculated and compared as well. This new approach has been justified due to available experimental data of McCool.

Keywords:Weibull Distribution, Skewness, Kurtosis, Fitting Indicator η

用于威布尔参数估计的高阶统计量

王桂金

原钢铁研究总院合金钢室，北京

收稿日期：2018年2月6日；录用日期：2018年2月20日；发布日期：2018年2月28日

摘 要

在材料及机械可靠性工程中，威布尔分布是常用的数学工具。然而，不同学者采用的算法通常给出不同的结果。本作者发展一种新算法，通过调整位置参数和名义最大寿命数据(99%)使得三参数威布尔形状参数和疲劳数据内禀高阶统计量斜度/峭度的形状参数κ一致而得到较好的拟合。本文以McCool的高温轴承钢为例，说明在轴承钢寿命分布的形状参数基本相同情况下材料的额定寿命L10如何受到冶炼工艺的影响。

关键词 :威布尔分布，斜度，峭度，拟合指数

Copyright © 2018 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

在材料及机械可靠性工程中，威布尔分布是常用的数学工具。然而，不同学者采用的算法通常给出不同的结果 [1] 。为了提高随机寿命威布尔参数拟合的准确度，本作者发展一种新算法 [2] ，通过调整位置参数t_o和名义最大寿命t_f (99%)使得拟合的三参数威布尔形状参数 $\kappa$ 和疲劳数据内禀的高阶统计量斜度/峭度的形状参数相近而得到更好的拟合。本文采用McCool高温轴承钢数据 [3] 为例说明这种新方法如何用来研究不同冶炼工艺对轴承额定寿命L10的影响。

计算方法

威布尔分布 [1]

威布尔累计失效率有如下三种形式：

三参数威布尔分布；

$f\left(t\right)=1-\exp \left\{-{\left(\frac{t-t_o}{\lambda }\right)}^{\kappa }\right\}\left(\lambda ,t-t_o,\kappa \right)>0$ (1-1)

式中f(t)是实测试样在时刻t的累计失效率。t_o：位置参数，κ：形状参数，λ：尺寸参数，是累计失效率达到63.2%需要的时间，它与形状参数 $\kappa$ 无关。如果t_o = 0，它变成两参数威布尔分布：

$f\left(t\right)=1-\exp \left\{-{\left(\frac{t}{\lambda }\right)}^{\kappa }\right\}\left(\lambda ,t,\kappa \right)>0$ (1-2)

进一步令t以 $\lambda$ 为单位计数，则得到归一化单参数威布尔分布：

$f\left(t\right)=1-\exp \left\{-t^{\kappa }\right\}\left(t,\kappa \right)>0$ (1-3)

威布尔失效密度分布函数就是上述诸式对时间t的导函数。威布尔分布参数的各种估算方法已经被详细研究过了。本文采用的是极大似然法 [1] [3] 。

2. 两参数威布尔分布的拟合

2.1. 极大似然法

设已知一组失效寿命序列(时间或转数) $t_1,t_2,\cdots ,t_f$ ，分别对应于累计失效率 $f\left(t_1\right),f\left(t_2\right),\cdots ,f\left(t_f\right)$ 。采用极大似然法公式(2)进行迭代运算可获得两参数威布尔分布的最佳参数k，然后将k代入式(3)算出尺寸参数 l。为保证结果可靠，本文计算进行到式(2)左边绝对值小于10⁻⁵后停止。

$\frac{1}{\kappa }+\left[\left(\frac{1}{N}\right){\displaystyle {\sum }_{i=1}^N\ln \left(t_i\right)}\right]-\left[{\displaystyle {\sum }_{i=1}^N{t}_i^k\ln \left(t_i\right)}\right]/\left({\displaystyle {\sum }_{i=1}^N{t}_i^k}\right)=0$ (2)

$\lambda ={\left[\left(\frac{1}{N}\right){\displaystyle {\sum }_{i=1}^N{t}_i^k}\right]}^{1/\kappa }$ (3)

由于t趋近无限大时威布尔分布累计失效才等于100%，为了方便，本文设定满足f(t_f) = 99%的t_f为名义全失效寿命。由于极大似然法只保证似然函数极大，为了提高威布尔曲线和试验寿命分布形状(不对称性，峰型)的拟合程度，有必要对拟合曲线和实际寿命曲线的斜度和峭度进行计算和比较。

2.2. 威布尔分布的斜度和峭度

已知形状参数为 $\kappa$ 的威布尔分布具有斜度skewness：

${\gamma }_1=\left[{\Gamma }_3-3{\Gamma }_2{\Gamma }_1+2{\Gamma }_2^3\right]/{\left[{\Gamma }_2-{\Gamma }_1^2\right]}^{3/2}$ (4)

而过盈峭度excess kurtosis (以下简称峭度)可按式(5)计算：

${\gamma }_2=\left\{\left({\Gamma }_4-4{\Gamma }_3{\Gamma }_1+6{\Gamma }_2{\Gamma }_1^2-3{\Gamma }_1^4\right)/{\left({\Gamma }_2-{\Gamma }_1^2\right)}^2\right\}-3$ (5)

式中， ${\Gamma }_i=\Gamma \left(1+i/\kappa \right)$ ，是gamma函数(i = 1, 2, 3, 4)。试样寿命数据的内禀斜度和峭度γ₁和γ₂可由通用统计软件(例如MS OFFICE EXCEL)计算，两者都没有量纲。

2.3. 拟合威布尔形状参数κ的具体步骤

为了评价威布尔分布参数拟合的优劣，本文采用如下拟合指数η：

$\eta =\left({\eta }_1+{\eta }_2\right)/2$ (6)

式中 ${\eta }_1=\kappa \left({\gamma }_1\right)/\kappa$ ， ${\eta }_2=\kappa \left({\gamma }_2\right)/\kappa$ 。并且 $\kappa \left({\gamma }_1\right)$ ， $\kappa \left({\gamma }_2\right)$ 分别是试验数组内禀斜度和峭度的形状参数。 $\kappa$ 是拟合威布尔分布的形状参数。当η₁和η₂皆为1或接近1时威布尔形状参数 $\kappa$ 拟合结果最佳。由于两参数威布尔分布往往达不到最佳拟合，引入可调整的非零位置参数t_o (三参数威布尔分布)以及名义全失效寿命t_f可以提高拟合度 [2] 。

具体步骤如下：

1) 先由小到大排列试样疲劳寿命值t_i，用极大似然法算出两参数威布尔分布的形状参数 $\kappa$ 以及试样寿命的内禀斜度κ，峭度γ₁及它们对应的形状参数 $\kappa \left({\gamma }_1\right)$ ， $\kappa \left({\gamma }_2\right)$ 和拟合指数η。如果η偏离期望值1较远，则

2) 令位置参数t_o为非零值,调整t_o并重复步骤1)，逐步使η靠近1，

3) 如果η仍然偏离1较远，可同时调整t_o和t_f重新计算η直到它接近最佳值1，这就得到拟合形状参数κ。

3.不同冶炼工艺下高温轴承钢的疲劳寿命

3.1. 试验设计

1982年McCool [4] 发表了对三种高温轴承试验钢：M50，T15和高铬不锈钢分别采用真空自耗，真空冶炼 + 真空电弧炉除气或粉末冶金法。总共五种组合的寿命实验代码见表1。

磨损试验设定：

试样尺寸：圆柱直径76.2毫米，长9.53毫米；

转速：10,000转/分，两个面对面191毫米圆盘；

最大接触应力：4826 Mpa即700,000 psi。

3.2. 试验结果及讨论

五种实验方案各有9个试样的疲劳寿命数据见表2。

最下面两行是McCool采用极大似然法(MLE)拟合两参数威布尔分布得到的斜率κ和10%额定累计失效率的寿命。

1) 原始数据的两参数威布尔分布拟合

本作者采用极大似然法(MLE)对表2数据计算两参数威布尔分布参数κ，λ及拟合指数η，结果见表3。其中κ和表2中McCool数值有些差异。按从小到大排列λ值，CVEM M50最小，PP CRB-7最大。参考McCool提供的多组疲劳数据失效机制相近判据，这五组实验的 ${\kappa }_{\max }/{\kappa }_{\min }$ 应小于2.61 ( $n=r=10,k=5$ )。本文数据 ${\kappa }_{\max }/{\kappa }_{\min }=3.4644/1.89622\approx 1.827$ 可认为合格。然而拟合指数全都不接近1。

符号说明(下同)：

t₀，t_f：疲劳失效开始时间和名义全失效时间(10⁶转)；

κ，λ：极大似然法拟合的威布尔形状参数(无量纲)和尺寸参数(10⁶转)；

表1. 五组冶炼工艺的代码

表2. 五组工艺的疲劳寿命数据，单位10⁶转

$\kappa \left({\gamma }_1\right)$ ， $\kappa \left({\gamma }_2\right)$ ：试样数据内禀斜度和内禀峭度对应的威布尔形状参数；

μ和σ：试样寿命数据的均值和标准差(10⁶转)；

δ²：拟合曲线相对于试样实验值的均方偏差(10⁶转)²：

${\delta }^2={\displaystyle {\sum }_{i=1}^9{\left\{f\left(t_i\right)-f\left(t_{i,\text{weib}}\right)\right\}}^2}/9$ (7)

式中 $f\left(t_i\right)$ 是试验寿命t_i的累计失效率， $f\left(t_{i,\text{weib}}\right)$ 是威布尔分布给出的累计失效率。于是，在五种冶炼工艺下的威布尔疲劳寿命累积失效率表达式分别是：

CVEMM50： $f\left(t\right)=1-\exp \left\{{\left(t/7.03549\right)}^{2.30045}\right\}t>0,t=t_1,t_2,\cdots ,t_f$

PP T15： $f\left(t\right)=1-\exp \left\{{\left(t/9.62393\right)}^{2.92675}\right\}t>0,t=t_1,t_2,\cdots ,t_f$

VIMVAR M50： $f\left(t\right)=1-\exp \left\{{\left(t/11.353\right)}^{1.89622}\right\}t>0,t=t_1,t_2,\cdots ,t_f$

PP M50： $f\left(t\right)=1-\exp \left\{{\left(t/11.862\right)}^{2.3971}\right\}t>0,t=t_1,t_2,\cdots ,t_f$

PP CRB-7： $f\left(t\right)=1-\exp \left\{{\left(t/16.407\right)}^{3.4644}\right\}t>0,t=t_1,t_2,\cdots ,t_f$

3.3. 数据经过t_o和t_f修正后按斜度和峭度拟合

五种冶炼工艺的随机疲劳寿命经t_o和t_f修正的最佳拟合指数η见表4。尺寸参数大小顺序略有改变，而且相应数值均比表3的小。因为位置参数t_o均是正值所致。如果把它加进尺寸参数λ，仍然和两参数拟合的尺寸参数λ有差别。这五组三参数威布尔分布 ${\kappa }_{\max }/{\kappa }_{\min }=1.08144/1.0204\approx 1.0598$ 比2.61小得多，而且拟合指数η都很接近1。这表明各种材料虽然成分不同，但材质都较均匀，夹杂物的数量，大小和分布都得到有效控制。此外，三种粉末冶金试样的尺寸参数λ都比较高，值得进一步研究。对冶炼及其他工艺对高温轴承寿命的影响有兴趣的读者请参考综合评述文献 [5] 。

表3. 五组疲劳数据拟合两参数威布尔分布的结果

NA*: No solution of k in normal range (0.5-5.0).

表4. 五组原始疲劳寿命经斜度和峭度修正后拟合三参数威布尔分布的结果

三参数威布尔疲劳寿命累积失效率表达式分别是

CEVM M50： $f\left(t\right)=1-\exp \left\{{\left(\left(t-2.9359\right)/3.7341\right)}^{1.03187}\right\}t-t_0>0,t=t_1,t_2,\cdots ,t_f$

VIMVARM50： $f\left(t\right)=1-\exp \left\{{\left(\left(t-5.46685\right)/4.0021\right)}^{1.0204}\right\}t-t_0>0,t=t_1,t_2,\cdots ,t_f$

PP T15： $f\left(t\right)=1-\exp \left\{{\left(\left(t-3.2237\right)/6.81288\right)}^{1.08144}\right\}t-t_0>0,t=t_1,t_2,\cdots ,t_f$

PP M50： $f\left(t\right)=1-\exp \left\{{\left(\left(t-2.66663\right)/10.279\right)}^{1.04561}\right\}t-t_0>0,t=t_1,t_2,\cdots ,t_f$

PP CRB-7： $f\left(t\right)=1-\exp \left\{{\left(\left(t-6.2439\right)/10.786\right)}^{1.07474}\right\}t-t_0>0,t=t_1,t_2,\cdots ,t_f$

把经斜度和峭度修正的寿命先按 $\left(t-t_0\right)/\lambda$ 归一化，再把45个数据融合成一组，用极大似然法重新拟合两参数威布尔分布得到如下结果，见表5。

融合后的威布尔分布形状参数κ = 1.03414很接近五组κ平均值1.0508，尺寸参数λ = 0.99977很接近1，而且拟合指数 $\eta \left({\gamma }_1,{\gamma }_2\right)=0.97245$ 也表明数据融合是成功的。

下面讨论均值μ，标准差σ和额定寿命L10。

3.4. 标准差σ，均值μ和额定寿命L10

已知形状参数κ和尺寸参数λ的威布尔分布的标准差σ是

$\sigma =\lambda {\left[\Gamma \left(1+2/\kappa \right)-{\Gamma }^2\left(1+1/\kappa \right)\right]}^{1/2}$ (8)

均值μ是：

$\mu =\lambda \left(\Gamma \left(1+1/\kappa \right)\right)$ (9)

随机寿命数据的标准差σ和均值μ可由通用统计软件算出，它们和威布尔分布的σ和μ是否一致也是检验拟合优劣的重要指标。

如表6~表7可见，不管是两参数或三参数拟合，同代码的威布尔参数和疲劳数组得到的均值μ都几乎相同。但是疲劳数据的标准差σ比拟合威布尔分布的σ高，因为σ的随机误差不能互相抵消。此外，有些γ₁，γ₂三参数拟合计算的额定寿命L10比两参数的小，因为它的起始点是从t_o算起。

表5. 五组疲劳寿命经斜度和峭度修正后归一化融合的两参数威布尔分布

表6. 两参数威布尔拟合五组随机疲劳寿命计算的标准差μ，均值σ和额定寿命L10，10⁶转

表7. 按γ₁γ₂拟合五组随机疲劳寿命的威布尔分布的标准差σ，均值μ和额定寿命L10，10⁶转

由于本文采用了新的拟合方法得到了较为可靠的L10排序，为选泽冶炼工艺提供参考。但是由于各厂家的冶炼，热处理，机加工和寿命试验的质量控制不完全一样，它们的实际寿命L10会有差别。因此对重要的轴承，还需要按实际使用环境作疲劳试验核实真正的L10。

4. 结论

本文用新算法对三种轴承钢/五种冶炼工艺的疲劳寿命拟合威布尔分布有如下结论：

1) 和现有威布尔参数拟合方法比较，通过调整位置参数t_o以及名义全失效时间t_f，五种冶炼工艺下拟合三参数威布尔斜率κ相对集中，拟合指数η更接近。各个冶炼工艺的疲劳寿命分布的形状参数κ，尺寸参数λ，标准差σ，均值μ和额定寿命L10都比较合理。这是它优于现有拟合方法之处。但它也对寿命数据质量提出了更高要求。如果经过多次调整t_o和t_f还得不到良好的拟合指数η，可能是样本数量不足，存在系统偏差或者疲劳数据服从其他分布所致。

2) 研究冶炼工艺对额定寿命L10的影响时，尺寸参数λ和形状参数κ都很重要。只有当不同数据组的κ值变化不大时，尺寸参数λ起更大作用。

3) 当需要了解在不同环境条件(诸如应力，温度，湿度，酸碱度等等)的额定寿命L10，最好按本文方法做系统的试验研究。

4) 本文只研究随机疲劳数据的拟合处理，未涉及失效微观机理，例如裂纹萌生及扩展的问题。本文结论只供参考。

文章引用

王桂金. 用于威布尔参数估计的高阶统计量
Higher Moment Statistics for Fitting Random Fatigue Life Data to Weibull Distribution Parameters[J]. 统计学与应用, 2018, 07(01): 65-71. http://dx.doi.org/10.12677/SA.2018.71009

参考文献 (References)