摘要： 設(shè)計(jì)階段的大型礦井提升機(jī)主軸承隨機(jī)振動(dòng)系統(tǒng)的設(shè)計(jì)可靠性分析因試驗(yàn)樣本不足，無(wú)法獲取振動(dòng)加速度響應(yīng)的完備概率信息。提出了概率信息不完備下大型礦井提升機(jī)主軸承設(shè)計(jì)可靠性分析的技術(shù)路線：滾動(dòng)軸承多尺度耦合系統(tǒng)動(dòng)力學(xué)模型?振動(dòng)加速度概率密度演化?概率密度演化路徑隨機(jī)過(guò)程建模?振動(dòng)功率譜密度的概率分布?基于條件概率的設(shè)計(jì)可靠度計(jì)算。運(yùn)用已采集的工況數(shù)據(jù)驅(qū)動(dòng)建立滾動(dòng)軸承多尺度耦合系統(tǒng)動(dòng)力學(xué)模型，開(kāi)展了振動(dòng)加速度概率密度演化研究；基于概率密度演化和Karhunen?Loève展開(kāi)，提出了滾動(dòng)軸承振動(dòng)加速度的非平穩(wěn)隨機(jī)過(guò)程建模方法，獲得了滾動(dòng)軸承振動(dòng)加速度隨機(jī)序列的孿生數(shù)據(jù)；研究了大型礦井提升機(jī)主軸承振動(dòng)加速度功率譜密度的概率分布，并計(jì)算了大型礦井提升機(jī)主軸承在服役時(shí)間上的設(shè)計(jì)可靠度指標(biāo)。

關(guān)鍵詞： 隨機(jī)振動(dòng)； 設(shè)計(jì)可靠性分析； 孿生數(shù)據(jù)； 大型礦井提升機(jī)主軸承； 概率信息不完備

中圖分類號(hào)： O324； TD534 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)： 1004-4523（2024）06-0915-13

DOI： 10.16385/j.cnki.issn.1004-4523.2024.06.002

1 概 述

設(shè)計(jì)階段的大型礦井提升機(jī)主軸承將在低轉(zhuǎn)速、重載荷和強(qiáng)沖擊的復(fù)雜環(huán)境中工作，其運(yùn)行穩(wěn)定性受到宏尺度系統(tǒng)動(dòng)態(tài)特性和微尺度材料磨損的共同作用。由軸承內(nèi)圈沖擊激勵(lì)產(chǎn)生的系統(tǒng)振動(dòng)會(huì)引起外圈滾道與滾動(dòng)體的局部接觸應(yīng)力變大。外圈滾道在局部高接觸應(yīng)力作用下會(huì)因疲勞磨損而導(dǎo)致其表面材料剝落，形成點(diǎn)蝕。同時(shí)，外圈滾道的點(diǎn)蝕故障激勵(lì)會(huì)引起滾動(dòng)軸承較大的振動(dòng)響應(yīng)［1］。實(shí)際中，滾動(dòng)軸承的宏尺度系統(tǒng)動(dòng)態(tài)特性和微尺度材料磨損還受到制造水平的影響，進(jìn)而導(dǎo)致滾動(dòng)軸承的振動(dòng)響應(yīng)遵循某種隨機(jī)過(guò)程［2］。此外，設(shè)計(jì)階段的大型礦井提升機(jī)主軸承因其價(jià)格昂貴和工況復(fù)雜，只能開(kāi)展小樣本試驗(yàn)，甚至無(wú)法開(kāi)展試驗(yàn)。雖然隨機(jī)過(guò)程的協(xié)方差函數(shù)可以從小樣本隨機(jī)序列中統(tǒng)計(jì)，但是概率分布無(wú)法從小樣本隨機(jī)序列中恢復(fù)，特別是四階矩［3］都隨時(shí)間變化的概率分布。因此，信息不完備下大型礦井提升機(jī)主軸承的設(shè)計(jì)可靠性分析是一個(gè)工程難題。

數(shù)字孿生技術(shù)［4］提供了很好的解決思路，在狀態(tài)量無(wú)法監(jiān)測(cè)的前提下，可以通過(guò)監(jiān)測(cè)環(huán)境變量來(lái)驅(qū)動(dòng)物理模型，進(jìn)而獲取滾動(dòng)軸承的振動(dòng)響應(yīng)并開(kāi)展可靠性研究。孿生數(shù)據(jù)驅(qū)動(dòng)的設(shè)計(jì)可靠性分析通過(guò)綜合考慮設(shè)計(jì)產(chǎn)品的制造水平和復(fù)雜工況，把未來(lái)服役的設(shè)計(jì)產(chǎn)品的運(yùn)行狀態(tài)全部提前映射出來(lái)，并基于孿生的狀態(tài)數(shù)據(jù)給出設(shè)計(jì)可靠度。

圖1詳細(xì)描述了孿生數(shù)據(jù)驅(qū)動(dòng)的大型礦井提升機(jī)主軸承隨機(jī)振動(dòng)系統(tǒng)的設(shè)計(jì)可靠性分析的技術(shù)路線。主要包括：第一，分析的對(duì)象是設(shè)計(jì)階段的主軸承；第二，復(fù)雜工況隨機(jī)分配給每個(gè)確定的主軸承；第三，主軸承的運(yùn)行狀態(tài)被提前孿生?；谠O(shè)計(jì)院確定的設(shè)計(jì)參數(shù)，建立共性的滾動(dòng)軸承多尺度系統(tǒng)動(dòng)力學(xué)模型；根據(jù)制造商的不完備制造水平信息，模型參數(shù)隨機(jī)化并確定個(gè)性的動(dòng)力學(xué)模型，且每個(gè)軸承產(chǎn)品的個(gè)性動(dòng)力學(xué)模型具備概率屬性；考慮到經(jīng)銷商隨機(jī)發(fā)貨給不同的煤礦企業(yè)，即每個(gè)軸承產(chǎn)品的服役工況具有不確定性，采集的工況數(shù)據(jù)被隨機(jī)分配給個(gè)性的動(dòng)力學(xué)模型；通過(guò)工況數(shù)據(jù)驅(qū)動(dòng)個(gè)性的動(dòng)力學(xué)模型，孿生出大量的振動(dòng)響應(yīng)，進(jìn)而建立振動(dòng)響應(yīng)的非平穩(wěn)非高斯隨機(jī)模型并計(jì)算這批軸承的可靠度。

以上孿生數(shù)據(jù)驅(qū)動(dòng)的大型礦井提升機(jī)主軸承隨機(jī)振動(dòng)系統(tǒng)的設(shè)計(jì)可靠性分析技術(shù)路線涉及概率分配的技術(shù)問(wèn)題，即個(gè)性的動(dòng)力學(xué)模型以多少概率出現(xiàn)和工況數(shù)據(jù)以多少概率分配給個(gè)性的動(dòng)力學(xué)模型，在這兩個(gè)概率的耦合下產(chǎn)生的振動(dòng)響應(yīng)也具備某種概率屬性。

本文以礦井提升機(jī)主軸承為研究對(duì)象，建立滾動(dòng)軸承宏微尺度耦合系統(tǒng)的動(dòng)力學(xué)模型；針對(duì)信息不完備的概率空間，提出一種基于線性矩的GF偏差代表點(diǎn)選取策略，并開(kāi)展?jié)L動(dòng)軸承振動(dòng)響應(yīng)的概率密度演化研究；針對(duì)概率密度演化缺失隨機(jī)序列信息，提出一種基于概率密度演化路徑的非平穩(wěn)隨機(jī)過(guò)程建模方法，進(jìn)而獲取隨機(jī)振動(dòng)功率譜密度的概率密度函數(shù)；最后計(jì)算礦井提升機(jī)主軸承的可靠度。

2 主軸承多尺度系統(tǒng)動(dòng)力學(xué)模型

2.1 宏尺度系統(tǒng)動(dòng)力學(xué)方程

如圖2所示，礦井提升機(jī)主軸承的內(nèi)圈同時(shí)受到X和Y軸方向上的徑向力和。滾動(dòng)軸承的非線性動(dòng)力學(xué)模型［1］包括內(nèi)、外圈水平和豎直方向的4個(gè)自由度以及單元諧振器豎直方向的1個(gè)自由度。其中，，和分別為內(nèi)圈、外圈和單元諧振器的質(zhì)量；，和分別為內(nèi)圈、外圈和單元諧振器的阻尼；，和分別為內(nèi)圈、外圈和單元諧振器的剛度；和分別為接觸力f在X和Y方向的分量。

根據(jù)牛頓第二定律，滾動(dòng)軸承的動(dòng)力學(xué)方程為：

（1）

式中 和分別為內(nèi)圈和外圈在X方向的位移；和分別為內(nèi)圈和外圈在Y方向的位移；為軸承的振動(dòng)響應(yīng)。

式（1）還可以寫(xiě)成矩陣的形式，其表達(dá)式為：

（2）

式中 M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；X為位移向量；F為接觸力向量；Q（t）為徑向力向量。

設(shè)滾動(dòng)軸承節(jié)圓直徑為，滾動(dòng)體直徑為，內(nèi)圈角速度為w，則保持架角速度為：

（3）

表示第i個(gè)滾動(dòng)體在時(shí)間t的角位置，其公式為：

；（4）

式中 為滾動(dòng)體的數(shù)目；為滾動(dòng)體的初始角位置。

滾動(dòng)軸承受內(nèi)圈徑向載荷的作用，第i個(gè)滾動(dòng)體在角位置的總接觸變形量為：

（5）

式中 為滾動(dòng)軸承的徑向間隙；為外圈滾道上點(diǎn)蝕的深度；為判斷函數(shù)，如下式所示：

（6）

式中 為外圈滾道上點(diǎn)蝕所在的位置角，為點(diǎn)蝕的角寬度。

根據(jù)非線性Hertz接觸理論，第i個(gè)滾動(dòng)體的接觸力為：

（7a）

（7b）

式中 為接觸剛度；為判斷函數(shù)，其表達(dá)式為：

。

2.2 微尺度材料滑動(dòng)磨損方程

假設(shè)外圈滾道的磨損為屈服線性磨損過(guò)程，即體積磨損率與法向載荷成正比。由Holm?Archard方程［5］可知：

（8）

式中 V為體積磨損量；s為滑動(dòng)距離；K為無(wú)量綱磨損系數(shù)；H為材料硬度；為法向載荷。

將式（8）的等號(hào)兩邊同時(shí)除以表觀接觸面積，可得：

（9）

式中 h為磨損深度；p為法向接觸壓力。

如圖3所示，滾動(dòng)體與滾道間的接觸面積近似為長(zhǎng)方形，其長(zhǎng)為 mm，寬為 mm，施加在第i個(gè)滾動(dòng)體上的法向載荷為，則在第i個(gè)滾動(dòng)體處施加在外圈滾道的平均法向接觸壓力為：

（10）

假設(shè)H為常數(shù)，則K和H可組合為無(wú)量綱磨損系數(shù)k，式（9）還可以表達(dá)為：

（11）

當(dāng)滾動(dòng)軸承旋轉(zhuǎn)一周時(shí)，每個(gè)滾動(dòng)體在外圈滾道上點(diǎn)蝕處的滑動(dòng)距離為，則外圈滾道上點(diǎn)蝕處的磨損深度為：

（12）

式中 h即為式（5）的點(diǎn)蝕深度。

2.3 多尺度耦合系統(tǒng)動(dòng)力學(xué)仿真

為了獲得礦井提升機(jī)主軸承宏微尺度耦合系統(tǒng)的振動(dòng)加速度響應(yīng)，利用Runge?Kutta法求解式（1）。礦井提升機(jī)主軸承多尺度耦合系統(tǒng)動(dòng)力學(xué)模型的主要參數(shù)如表1所示，其中，，，和分別為和的均值。

圖4（a）～（c）分別為礦井提升機(jī)主軸承多尺度耦合系統(tǒng)的振動(dòng)加速度時(shí)間歷程。隨著外圈滾道磨損深度的增加，滾動(dòng)軸承加速度的振幅逐漸增大。在服役時(shí)間 s內(nèi)，滾動(dòng)軸承加速度的最大振幅為0.998 m/s2；在服役時(shí)間 s內(nèi)，滾動(dòng)軸承加速度的最大振幅為2.769 m/s2；在服役時(shí)間 s內(nèi)，滾動(dòng)軸承加速度的最大振幅為4.441 m/s2。

3 主軸承隨機(jī)振動(dòng)的不確定性量化

3.1 信息不完備下概率空間選點(diǎn)的線性矩法

概率空間由d維隨機(jī)向量X構(gòu)成，且各分量的概率分布未知。實(shí)際中，只有少量樣本信息能被知曉?；谟邢迋€(gè)（ns>200）隨機(jī)樣本，的前4個(gè)線性矩［6?7］分別定義為：

（13a）

（13b）

（13c）

（13d）

式中 為概率權(quán)重矩，其表達(dá)式為：

（14）

式中 為的樣本容量；為的第j個(gè)次序統(tǒng)計(jì)量。

為了標(biāo)準(zhǔn)化的高階線性矩，它的偏度和峰度用線性矩比來(lái)描述，如下式所示：

（15a）

（15b）

式中 和分別為的線性偏度和線性峰度。

根據(jù)三次正態(tài)變換多項(xiàng)式［6］，可近似為：

（16）

式中 為三次多項(xiàng)式函數(shù)；為標(biāo)準(zhǔn)正態(tài)隨機(jī)變量的樣本；，，和為多項(xiàng)式系數(shù)，其解析表達(dá)式見(jiàn)附錄A。

按照變量的等概率變換原則，的累積分布函數(shù)表達(dá)為：

（17）

式中 和分別為和的累積分布函數(shù)，和的對(duì)應(yīng)關(guān)系如式（16）所示。

相較傳統(tǒng)的Pearson相關(guān)系數(shù)，Kendall（或Spearman）秩相關(guān)系數(shù)能在非線性變換的前提下保持不變。由于和的對(duì)應(yīng)關(guān)系是非線性變換，因此，秩相關(guān)系數(shù)被用來(lái)度量和未知變量的關(guān)聯(lián)程度。

假設(shè)和的Kendall秩相關(guān)系數(shù)為，那么，中間標(biāo)準(zhǔn)正態(tài)隨機(jī)變量和的Kendall秩相關(guān)系數(shù)也為。與對(duì)應(yīng)的Pearson相關(guān)系數(shù)為：

（18）

對(duì)應(yīng)的Pearson相關(guān)矩陣描述為：

（19）

基于正態(tài)隨機(jī)變量的線性變換性質(zhì)，Z可以用U來(lái)表示，其表達(dá)式為：

（20）

式中 A為R的Cholesky分解；為A的轉(zhuǎn)置，其表達(dá)式為：

。

根據(jù)式（16）和（20），和最終可以分別表示為：

（21a）

（21b）

3.2 概率空間剖分與賦得概率計(jì)算

代表點(diǎn)選取的基本原則是在均勻抽樣的基礎(chǔ)上確保經(jīng)驗(yàn)累積分布函數(shù)的偏差很小。針對(duì)信息不完備的概率空間，本節(jié)提出了一種改進(jìn)的GF偏差代表點(diǎn)［8］選取策略。

一個(gè)完整的概率空間被剖分成一系列子域（），即滿足：（1）；（2），。在子域內(nèi)的代表點(diǎn)用表示，其中，是一個(gè)d維隨機(jī)向量。的Voronoi區(qū)域用表示，其表達(dá)式為：

（22）

式中 為概率空間的代表點(diǎn)集；為N個(gè)d維隨機(jī)向量組成的空間。

顯然，的體積需要用Monte Carlo方法獲取，其代表點(diǎn)的賦得概率與的體積有關(guān)。

圖5為利用Voronoi區(qū)域計(jì)算代表點(diǎn)的賦得概率的示意圖。以的概率空間為例，基于外接圓法計(jì)算Voronoi區(qū)域的面積，并采用微元法近似其代表點(diǎn)的賦得概率。每個(gè)Voronoi區(qū)域的外接圓的最大半徑可以表達(dá)為：

（23）

如圖5所示，以代表點(diǎn)為圓心，以為半徑作圓，圓的面積為。采用Monte Carlo方法在圓內(nèi)均勻撒入個(gè)測(cè)試點(diǎn) ，計(jì)算每個(gè)測(cè)試點(diǎn)與代表點(diǎn)的距離。若測(cè)試點(diǎn)滿足公式（23），則：

（24）

式中 為與相鄰的代表點(diǎn)，則屬于代表點(diǎn)的Voronoi區(qū)域。

記測(cè)試點(diǎn) 中屬于的Voronoi區(qū)域的數(shù)目為，則Voronoi區(qū)域的面積可以近似為：

（25）

位于Voronoi區(qū)域的每個(gè)測(cè)試點(diǎn)代表的面積可以近似為：

（26）

采用微元法，位于Voronoi區(qū)域的代表點(diǎn)的賦得概率計(jì)算為：

（27）

式中 為測(cè)試點(diǎn)的聯(lián)合概率密度函數(shù)。

根據(jù)式（16）和（17），可以表達(dá)為：

（28）

式中 為標(biāo)準(zhǔn)正態(tài)分布密度函數(shù)；為測(cè)試點(diǎn)的等概率變換點(diǎn)，它們對(duì)應(yīng)關(guān)系如式（16）所示；為的第j個(gè)分量。若考慮的兩個(gè)分量間的相關(guān)性，可以參照式（21a）和（21b）。

對(duì)代表點(diǎn)（）的概率進(jìn)行歸一化處理，如下式所示：

（29）

式中 為代表點(diǎn)的最終賦得概率。

實(shí)際上，概率空間只能被剖分成有限個(gè)且盡可能少的子域。因此，GF偏差通常被用來(lái)作為概率空間最佳剖分的指標(biāo)，其表達(dá)式為：

（30）

式中 為第j個(gè)變量的邊際分布函數(shù)，其表達(dá)式如式（17）所示；為概率空間剖分影響下第j個(gè)變量的經(jīng)驗(yàn)邊際分布函數(shù)，可以表達(dá)為：

（31）

式中 為指示函數(shù)；為代表點(diǎn)的第j個(gè)分量。

代表點(diǎn)集的GF偏差越小，概率空間的剖分越合理。因此，概率空間的剖分可以轉(zhuǎn)化為最優(yōu)化問(wèn)題。這種最優(yōu)化問(wèn)題可以通過(guò)遺傳算法或其他優(yōu)化方法解決，但要占據(jù)大量的計(jì)算資源。本節(jié)采用點(diǎn)集重排法實(shí)現(xiàn)快速的概率空間剖分。該方法包括兩個(gè)步驟：（1）生成初始點(diǎn)集；（2）重新排列以減少GF偏差。

由Sobol序列［9］生成在單位超立方體的d維均勻點(diǎn)集，如下式所示：

（32）

根據(jù)式（19a）和（19b），在物理空間被第一次均勻化的初始點(diǎn)集可以定義為：

（33）

式中 為在標(biāo)準(zhǔn)正態(tài)空間的點(diǎn)集，其表達(dá)式為：

（34）

式中 為的逆函數(shù)。

為了確保初始點(diǎn) 所處的Voronoi區(qū)域的彼此接近，在每個(gè)維度上被第二次均勻化，如下式所示：

（35）

式中 為在標(biāo)準(zhǔn)正態(tài)空間被第二次均勻化的點(diǎn)集，定義為：

（36）

式中 為在物理空間內(nèi)第j個(gè)維度上被第一次均勻化的點(diǎn)集，當(dāng)時(shí)，。

根據(jù)式（29）和（31），的最終賦得概率和經(jīng)驗(yàn)邊際分布函數(shù)被分別確定。為了減少的GF偏差，在經(jīng)驗(yàn)邊際分布函數(shù)上被第三次均勻化，如下式所示：

（37）

式中 為在標(biāo)準(zhǔn)正態(tài)空間被第三次均勻化的點(diǎn)集，定義為：

（38）

式中 為在物理空間內(nèi)第j個(gè)維度上被第二次均勻化的點(diǎn)集；為代表點(diǎn)的最終賦得概率。當(dāng)時(shí)，，。

在信息不完備的概率空間中，由以上的點(diǎn)集重排法所確定的代表點(diǎn)集 通常只是次佳的。但是，的精度在實(shí)際應(yīng)用中足夠滿足大部分工程的要求。表2為影響礦井提升機(jī)主軸承多尺度耦合系統(tǒng)振動(dòng)響應(yīng)的隨機(jī)變量。

如圖6所示，代表點(diǎn)數(shù)量q是影響選點(diǎn)策略精度的主要因素，而測(cè)試點(diǎn)數(shù)量對(duì)選點(diǎn)策略精度的影響較小。當(dāng)時(shí)，選點(diǎn)策略能達(dá)到較好的精度效果，其GF偏差的范圍為。當(dāng)時(shí)，采用獲得的選點(diǎn)策略的GF偏差為0.0271；而采用獲得的選點(diǎn)策略的GF偏差為0.0303。因此，的增加并不能顯著地提高選點(diǎn)策略的穩(wěn)定性。在工程應(yīng)用中，可以選取較大的q和較小的，在多次嘗試的情況下選擇滿足工程要求的代表點(diǎn)集。

3.3 隨機(jī)振動(dòng)的概率密度演化

礦井提升機(jī)主軸承多尺度耦合系統(tǒng)動(dòng)力學(xué)方程的參數(shù)用一系列基本隨機(jī)變量表示。根據(jù)式（2），主軸承的隨機(jī)動(dòng)力學(xué)方程可以表示為：

（39）

基于隨機(jī)動(dòng)力系統(tǒng)的概率守恒原理，式（39）的廣義概率密度演化方程［10］為：

（40）

式中 為隨機(jī)變量A和的聯(lián)合概率密度函數(shù)；為所關(guān)注的物理量，即主軸承的振動(dòng)加速度響應(yīng)；為所關(guān)注的物理量的狀態(tài)變化速率；a和x分別為A和對(duì)應(yīng)的樣本。

式（33）的初始條件為：

（41）

式中 為的聯(lián)合概率密度函數(shù)；為感興趣的物理量的初始值；為Dirac?Delta函數(shù)。

采用有限差分法求解式（40）和（41），的概率密度函數(shù)可以計(jì)算為：

（42）

圖7為礦井提升機(jī)主軸承宏微尺度耦合系統(tǒng)在 s時(shí)間歷程的振動(dòng)加速度概率密度演化。圖7（a）～（c）反映了滾動(dòng)軸承振動(dòng)加速度的概率分布隨時(shí)間的變化趨勢(shì)。滾動(dòng)軸承振動(dòng)加速度的變化范圍為-6～6 m/s2，其中，振動(dòng)加速度以較大的概率出現(xiàn)在-2～2 m/s2。礦井提升機(jī)主軸承在全服役時(shí)間歷程的振動(dòng)加速度概率密度演化見(jiàn)附錄B。

4 主軸承動(dòng)態(tài)可靠性分析

4.1 概率密度演化路徑非平穩(wěn)隨機(jī)過(guò)程建模

概率密度演化缺失隨機(jī)序列信息，無(wú)法從頻譜分析的角度對(duì)主軸承開(kāi)展可靠性分析。本節(jié)提出了一種基于概率密度演化路徑的非平穩(wěn)隨機(jī)過(guò)程建模方法，其基本思想是使用潛在的非平穩(wěn)高斯隨機(jī)過(guò)程去近似目標(biāo)的非平穩(wěn)非高斯隨機(jī)過(guò)程。

為定義在一個(gè)概率空間上的隨機(jī)過(guò)程，表示t時(shí)刻的A在樣本點(diǎn)處的值，它的前4個(gè)中心矩被表示為，，和。采用梯形法對(duì)3.3節(jié)的概率密度演化進(jìn)行數(shù)值積分，的前4個(gè)中心矩［11］能夠被確定。

根據(jù)三次正態(tài)變換多項(xiàng)式［12］，可近似為：

（43）

式中 ，，和為三次多項(xiàng)式的系數(shù)；為非平穩(wěn)高斯隨機(jī)過(guò)程。

已知的協(xié)方差函數(shù)為，它的自相關(guān)函數(shù)定義為：

（44）

根據(jù)式（43）和（44），的自相關(guān)函數(shù)與的自相關(guān)函數(shù)的關(guān)系可以表達(dá)為：

（45）

式中 為的自相關(guān)函數(shù)，其解析表達(dá)式見(jiàn)附錄C。

由于和，的協(xié)方差函數(shù)與它的自相關(guān)函數(shù)相等，即

（46）

基于Mercer定理［13］，的譜分解可以表達(dá)為：

（47）

式中 和分別為第二類Fredholm積分方程［14］的特征值和特征函數(shù)，如下式所示：

（48）

式中 為的時(shí)間域，有一個(gè)完整的正交集，并滿足下式：

（49）

式中 為Kronecker delta函數(shù)。

基于Karhunen?Loève展開(kāi)［15］，表達(dá)為：

（50）

式中 m為展開(kāi)項(xiàng)的數(shù)目；為一系列獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量。

將式（50）帶入式（43），即得非平穩(wěn)非高斯隨機(jī)過(guò)程。

圖8為譜分解的特征函數(shù)和特征值?；谑剑?8）求解產(chǎn)生個(gè)特征函數(shù)和特征值。

圖9（a）的為給定的協(xié)方差函數(shù)。由3.2節(jié)概率空間剖分產(chǎn)生個(gè)滾動(dòng)軸承振動(dòng)加速度隨機(jī)序列樣本，并統(tǒng)計(jì)得到的協(xié)方差函數(shù)為?；谑剑?4）～（46）求解得到的協(xié)方差函數(shù)。圖9（b）為的Karhunen?Loève近似，其近似表達(dá)式見(jiàn)式（47）。

圖10為由概率密度演化路徑非平穩(wěn)隨機(jī)過(guò)程產(chǎn)生的106個(gè)振動(dòng)加速度隨機(jī)序列。其振動(dòng)加速度幅值和協(xié)方差函數(shù)分別與圖4（c）的振動(dòng)加速度幅值和協(xié)方差函數(shù)存在著某種“相似性”，即設(shè)計(jì)階段主軸承加速度狀態(tài)量的“孿生”。

4.2 隨機(jī)振動(dòng)頻譜概率密度函數(shù)分布

如圖10所示，基于4.1節(jié)建立的隨機(jī)過(guò)程，采用Monte Carlo方法產(chǎn)生106個(gè)加速度隨機(jī)序列。對(duì)每個(gè)確定的隨機(jī)序列分別進(jìn)行傅里葉變換，如下式所示：

（51）

式中 f為頻率。

對(duì)于離散的頻率 ，其對(duì)應(yīng)的傅里葉變換定義為：

（52）

式中 為采樣間隔；為時(shí)間段的采樣點(diǎn)數(shù)。

根據(jù)維納?欣欽定理，的功率譜密度定義為：

（53）

對(duì)于落在間隔內(nèi)的每個(gè)確定的隨機(jī)序列的功率譜密度的片段，統(tǒng)計(jì)其功率譜密度的平均值，如下式所示：

（54）

式中 為落在間隔內(nèi)的數(shù)量；為間隔的中心頻率，其表達(dá)式為：

。

以1/2倍頻程增加用于計(jì)算下一個(gè)間隔內(nèi)的功率譜密度的平均值，。直到間隔包括最后的頻率。

將轉(zhuǎn)換為以dB為單位的數(shù)值，其表達(dá)式為：

（55）

對(duì)于每個(gè)確定的，的概率密度為：

（56）

式中 為106個(gè)隨機(jī)序列的落在某個(gè)長(zhǎng)度為1 dB的間隔內(nèi)的數(shù)量。

圖11為礦井提升機(jī)主軸承振動(dòng)加速度功率譜密度的概率密度函數(shù)。如圖11（a）所示，在加速度頻率34～40 Hz內(nèi)，滾動(dòng)軸承振動(dòng)加速度每間隔兩個(gè)頻率出現(xiàn)一個(gè)較大的峰值，且較大的峰值以特定的概率出現(xiàn)在之間。圖11（b）反映了在頻率 Hz上的滾動(dòng)軸承振動(dòng)加速度較大峰值的概率分布。

4.3 基于條件概率的可靠度計(jì)算

礦井提升機(jī)主軸承外圈的故障頻率定義為：

（57）

式中 和為隨機(jī)變量，其概率分布見(jiàn)表2。

當(dāng)出現(xiàn)在故障頻率上的加速度功率譜密度的峰值超過(guò)失效閾值時(shí)，礦井提升機(jī)主軸承視為失效。因此，礦井提升機(jī)主軸承的失效概率表達(dá)為：

（58）

式中 為失效閾值；為在“”發(fā)生條件下“大于”發(fā)生的概率，如式（56）和圖11（b）所示；為“”發(fā)生的概率。

根據(jù)式（57）和二階Taylor級(jí)數(shù)展開(kāi)［16］，“”的前4個(gè)中心矩可被確定，它們分別為，，和?；谌握龖B(tài)變換多項(xiàng)式，的概率密度函數(shù)可表達(dá)為：

（59）

式中 ，和為三次正態(tài)變換的多項(xiàng)式系數(shù)；和的對(duì)應(yīng)關(guān)系為，其中。

圖12為礦井提升機(jī)主軸承外圈故障頻率的概率密度函數(shù)，的定義域?yàn)?。相較圖11（a），滾動(dòng)軸承振動(dòng)加速度在頻率 內(nèi)共出現(xiàn)5個(gè)峰值，其中較大峰值的概率分布如圖11（b）所示。

將的定義域均等分成5段，每段包括對(duì)應(yīng)的加速度頻率。定義式（58）的功率譜密度的失效閾值為。設(shè)計(jì)階段礦井提升機(jī)主軸承在服役時(shí)間 s內(nèi)的最終失效概率計(jì)算為：

（60）

相應(yīng)的可靠度為：

（61）

5 結(jié) 論

針對(duì)設(shè)計(jì)階段的大型提升機(jī)主軸承隨機(jī)振動(dòng)的設(shè)計(jì)可靠性分析缺少試驗(yàn)樣本的工程難題，本文提出了孿生數(shù)據(jù)驅(qū)動(dòng)的大型礦井提升機(jī)主軸承隨機(jī)振動(dòng)系統(tǒng)的設(shè)計(jì)可靠性分析的技術(shù)路線。該技術(shù)路線總結(jié)如下：

（1）提出了大型礦井提升機(jī)主軸承宏微尺度耦合系統(tǒng)的動(dòng)力學(xué)模型，運(yùn)用已采集的載荷歷程驅(qū)動(dòng)滾動(dòng)軸承的動(dòng)力學(xué)模型，并獲取振動(dòng)加速度響應(yīng)。

（2）針對(duì)概率分配的技術(shù)問(wèn)題：個(gè)性的動(dòng)力學(xué)模型以多少概率出現(xiàn)和工況數(shù)據(jù)以多少概率分配給個(gè)性的動(dòng)力學(xué)模型，提出了信息不完備的概率空間剖分和賦得概率計(jì)算方法。該方法以GF偏差作為概率空間最佳剖分的指標(biāo)，保證了工程的精度要求。

（3）現(xiàn)有的概率密度演化缺失隨機(jī)序列信息，無(wú)法從頻譜分析的角度對(duì)主軸承開(kāi)展可靠性分析。提出了概率密度演化路徑非平穩(wěn)隨機(jī)過(guò)程建模方法。相較傳統(tǒng)的隨機(jī)過(guò)程建模方法，該方法中的分布類型無(wú)需假設(shè)，隨機(jī)序列的每個(gè)樣本由概率密度演化給定的時(shí)間上的分布類型和從物理模型中統(tǒng)計(jì)的協(xié)方差函數(shù)共同決定。

（4）提供了滾動(dòng)軸承隨機(jī)振動(dòng)的加速度功率譜密度的概率分布。該功率譜密度的概率分布反映了發(fā)生在頻率上的加速度功率譜密度的峰值概率。若發(fā)生在故障頻率上的加速度功率譜密度的峰值很小，則滾動(dòng)軸承視為安全。

（5）提供了大型礦井提升機(jī)主軸承在服役時(shí)間歷程上的可靠度指標(biāo)，為考慮批量生產(chǎn)的滾動(dòng)軸承設(shè)計(jì)提供了依據(jù)。

參考文獻(xiàn)：

［1］ 胡愛(ài)軍， 許莎， 向玲， 等. 滾動(dòng)軸承外圈多點(diǎn)故障特征分析［J］. 機(jī)械工程學(xué)報(bào)， 2020， 56（21）： 110-120.

HU Aijun， XU Sha， XIANG Ling， et al. Characteristic analysis of multi-point faults on the outer race of rolling element bearing［J］. Journal of Mechanical Engineering， 2020， 56（21）： 110-120.

［2］ 張文虎， 鄧四二， 陳國(guó)定， 等. 圓柱滾子軸承半物理仿真試驗(yàn)技術(shù)研究［J］. 機(jī)械工程學(xué)報(bào)， 2018， 54（19）： 78-87.

ZHANG Wenhu， DENG Sier， CHEN Guoding， et al. Hardware-in-the-loop simulation technology of cylindrical roller bearing［J］. Journal of Mechanical Engineering， 2018， 54（19）： 78-87.

［3］ 劉宇， 李天翔， 劉闊， 等. 基于四階矩法車削顫振可靠性研究［J］. 機(jī)械工程學(xué)報(bào)， 2016， 52（20）： 193-200.

LIU Yu， LI Tianxiang， LIU Kuo， et al. Chatter reliability of turning processing system based on fourth moment method［J］. Journal of Mechanical Engineering， 2016， 52（20）： 193-200.

［4］ 付洋， 曹宏瑞， 郜偉強(qiáng)， 等. 數(shù)字孿生驅(qū)動(dòng)的航空發(fā)動(dòng)機(jī)渦輪盤(pán)剩余壽命預(yù)測(cè)［J］. 機(jī)械工程學(xué)報(bào)， 2021， 57（22）： 106-113.

FU Yang， CAO Hongrui， GAO Weiqiang， et al. Digital twin driven remaining useful life prediction for aero-engine turbine discs［J］. Journal of Mechanical Engineering， 2021， 57（22）： 106-113.

［5］ 尹家寶， 盧純， 全鑫， 等. 列車制動(dòng)塊磨損行為動(dòng)態(tài)演變數(shù)值分析［J］. 機(jī)械工程學(xué)報(bào)， 2021， 57（18）： 204-213.

YIN Jiabao， LU Chun， QUAN Xin， et al. Analysis of wear behavior dynamic evolution on railway brake pad［J］. Journal of Mechanical Engineering， 2021， 57（18）： 204-213.

［6］ ULRYCH T J， VELIS D R， WOODBURY A D， et al. L-moments and C-moments［J］. Stochastic Environmental Research and Risk Assessment， 2000， 14： 50-68.

［7］ CAO Shuang， LU Hao， PENG Yuxing， et al. A novel fourth-order L-moment reliability method for L-correlated variables［J］. Applied Mathematical Modelling， 2021， 95： 806-823.

［8］ CHEN Jianbing， YANG Junyi， LI Jie. A GF-discrepancy for point selection in stochastic seismic response analysis of structures with uncertain parameters［J］. Structural Safety， 2016， 59： 20-31.

［9］ RADOVI? I， SOBOL I M， TICHY R F. Quasi-Monte Carlo methods for numerical integration： comparison of different low discrepancy sequences［J］. Monte Carlo Methods and Applications， 1996， 2（1）： 1-14.

［10］ LIU Gang， GAO Kai， YANG Qingshan， et al. Improvement to the discretized initial condition of the generalized density evolution equation［J］. Reliability Engineering & System Safety， 2021， 216： 107999.

［11］ LU Hao， CAO Shuang， ZHU Zhencai， et al. An improved high order moment-based saddlepoint approximation method for reliability analysis［J］. Applied Mathematical Modelling， 2020， 82（6）： 836-847.

［12］ LU Zhaohui， CAI Chaohuang， ZHAO Yangang， et al. Normalization of correlated random variables in structural reliability analysis using fourth-moment transformation［J］. Structural Safety， 2020， 82： 101888.

［13］ CHIOU J M， CHEN Y T， YANG Y F. Multivariate functional principal component analysis：a normalization approach［J］. Statistica Sinica， 2014， 24： 1571-1596.

［14］ WANG Xingtao， LI Yuanmin. Numerical solution of fredholm integral equation of the second kind by general legendre wavelets［J］. International Journal of Innovative Computing， Information and Control， 2012， 8： 799-805.

［15］ PHOON K K， HUANG S P， QUEK S T. Simulation of second-order processes using Karhunen?Loève expansion［J］. Computers & Structures， 2002， 80（12）： 1049-1060.

［16］ 鄭宏偉， 孟廣偉， 李鋒， 等. 基于高階矩最大熵方法的結(jié)構(gòu)混合可靠性分析［J］. 機(jī)械工程學(xué)報(bào)， 2021， 57（14）： 282-290.

ZHENG Hongwei， MENG Guangwei， LI Feng， et al. Hybrid reliability analysis for structures based on high-order moments and maximum entropy method［J］. Journal of Mechanical Engineering， 2021， 57（14）： 282-290.

Twin data-driven design reliability analysis of random vibration system for main bearing of a large mine hoist

CAO Shuang1，2， LU Hao1，2， ZHU Zhen-cai1，2， ZHANG Yi-min3

（1. School of Mechanical & Electrical Engineering， China University of Mining and Technology， Xuzhou 221008， China；2. Jiangsu Key Laboratory of Mine Mechanical and Electrical Equipment， China University of Mining and Technology，Xuzhou 221008， China； 3. School of Mechanical and Power Engineering， Shenyang University of Chemical Technology， Shenyang 110142， China）

Abstract： In the design phase of a mine hoist’s main bearing， the reliability analysis of its random vibration cannot obtain complete probability information of the vibration acceleration response due to insufficient experimental samples. This paper proposes a new technical route for the reliability analysis of the main bearing of a mine hoist under incomplete probability information. This route includes the dynamic model of a multi-scale coupling system for rolling element bearings， the probability density evolution of vibration acceleration， stochastic process modeling of the probability density evolution path， the probability distribution of vibration power spectral density， and the calculation of design reliability based on conditional probability. By using the collected condition data to drive the established dynamic model of multi-scale coupling system for rolling element bearings， the probability density evolution of vibration acceleration is carried out. Based on the probability density evolution and Karhunen-Loève expansion， a modeling approach for the non-stationary random process of vibration acceleration of rolling element bearings is proposed. This approach obtains the twin data of the random sequence of vibration acceleration for rolling element bearings. The probability distribution of the vibration power spectral density for the main bearing of a mine hoist is studied， and the reliability index of the main bearing of a mine hoist within the service time is calculated.

Key words： random vibration； design reliability analysis； twin data； main bearing of a large mine hoist； incomplete probability information

作者簡(jiǎn)介： 曹 爽（1993―），男，博士研究生。 E-mail： caoshuangyc@163.com。

通訊作者： 盧 昊（1985―），男，博士，副教授。 E-mail： haolucumt@163.com。

附錄A：基于正態(tài)變換的三次多項(xiàng)式系數(shù)的解析表達(dá)式

設(shè)式（16）中的前4個(gè)線性矩與式（16）中的前4個(gè)線性矩相等，可以得到關(guān)于多項(xiàng)式系數(shù)的方程組，如下式所示：

（A1）

（A2）

（A3）

（A4）

其中，和的表達(dá)式分別為：

，。

求解式（A1）～（A4），多項(xiàng)式系數(shù)，，和可被確定為：

，，。

附錄B：全服役時(shí)間歷程的加速度概率密度演化

附錄C：的解析表達(dá)式

（1） 當(dāng)時(shí)：

（C1）

（2） 當(dāng)且時(shí)：

（C2）

（3） 當(dāng)，且時(shí)：

（C3）

（4） 當(dāng)，且時(shí)：

（C4）

（5） 當(dāng)且時(shí)：

（C5）

（6） 當(dāng)且時(shí)：

（C6）

其中，g，，，B，C和D的表達(dá)式分別為：

（C7-1）

（C7-2）

（C7-3）

（C7-4）

（C7-5）

（C7-6）

其中，，和的表達(dá)式分別為：

，

，。