俞 水，汪忠來

(電子科技大學(xué) 機(jī)械與電氣工程學(xué)院，成都 611731)

1 引 言

在機(jī)械產(chǎn)品的壽命周期內(nèi)，產(chǎn)品的材料性能、工作環(huán)境和載荷效應(yīng)等呈現(xiàn)時變性和非線性，這將給產(chǎn)品的可靠性分析和設(shè)計(jì)帶來新的挑戰(zhàn)。由于傳統(tǒng)的靜態(tài)可靠性分析方法[1-3]忽略了時變不確定性對可靠性的影響，從而無法獲得滿足精度的可靠度。因此，針對產(chǎn)品的時變不確定性進(jìn)行時變可靠性分析具有重要的工程意義。

時變可靠性分析是目前可靠性分析和設(shè)計(jì)的熱點(diǎn)問題，已有大量的研究工作，主要可以分為解析法和數(shù)值仿真法兩類。數(shù)值仿真方法利用隨機(jī)過程時間序列離散，將時變可靠性問題轉(zhuǎn)化成高維的靜態(tài)可靠性問題，主要包括蒙特卡洛仿真方法[4]、重要度抽樣技術(shù)[5]和子集模擬[6]。蒙特卡洛仿真雖然可以獲得精確解，但是針對高可靠度評估問題計(jì)算效率極低，因此常將蒙特卡洛仿真方法作為評價其他可靠度分析方法的一個標(biāo)準(zhǔn)。重要度抽樣技術(shù)在處理高維性能函數(shù)時往往無法給定合適的重要性抽樣密度，而且重要度抽樣密度的選擇依賴于時效域的相關(guān)信息。子集模擬是最近出現(xiàn)的高可靠度評估的有效方法，但其計(jì)算效率取決于性能函數(shù)的線性程度。此外，國內(nèi)外學(xué)者從解析法的角度開展了大量的研究工作，主要包括極值響應(yīng)面方法[7]、Gamma過程[8]、復(fù)合隨機(jī)過程法[9]、攝動法[10]、復(fù)合極限狀態(tài)法[11]和穿越率方法[12]。極值響應(yīng)面方法、Gamma過程和復(fù)合隨機(jī)法在評估時變可靠性時需要假設(shè)系統(tǒng)參數(shù)或者性能服從特定的分布，將產(chǎn)生模型近似誤差；攝動法對系統(tǒng)參數(shù)有較強(qiáng)的依賴性；復(fù)合極限狀態(tài)法在處理非線性程度較高的性能函數(shù)時，將產(chǎn)生無法接受的誤差；相對而言，穿越率方法在計(jì)算效率方面優(yōu)勢較大，引起學(xué)術(shù)界和工業(yè)界的廣泛關(guān)注。Rice[12]首先提出了穿越率方法，隨后許多學(xué)者在此基礎(chǔ)上開展了大量的研究，主要包括可微高斯過程法[13]、矩形波更新過程法[14]、拉普拉斯積分法[15]、PHI2方法[16]和PHI2+方法[17]。由于近似時間積分的復(fù)雜性和泊松分布假設(shè)，可微高斯過程法、矩形波更新過程法和拉普拉斯積分法只能處理特定隨機(jī)過程條件下的穿越率計(jì)算問題，且計(jì)算精度不高。PHI2和PHI2+方法利用并聯(lián)系統(tǒng)可靠性框架來提高計(jì)算精度以及擴(kuò)寬穿越率計(jì)算方法的適用范圍。但現(xiàn)有研究表明，在處理非單調(diào)性系統(tǒng)時，基于并聯(lián)系統(tǒng)可靠性框架的穿越率計(jì)算方法顯現(xiàn)出較低的計(jì)算精度[18]。

本文針對產(chǎn)品全壽命周期的時變可靠性問題，建立了首次穿越點(diǎn)的概率分布模型，提出了一種求解時變可靠性的計(jì)算方法，可獲得產(chǎn)品在壽命周期內(nèi)的可靠性累計(jì)概率密度函數(shù)。

2 時變可靠性簡介

時變可靠性是指產(chǎn)品在規(guī)定的時間內(nèi)和條件下承受時變不確定性作用能夠完成規(guī)定功能的概率。因此在壽命周期[0,T]內(nèi)，產(chǎn)品的可靠度R(0,T)可以表示為

R(0,T)=Pr{g(X,Y(t),t)>0,?t∈[0,T]}

(1)

式中X是n維隨機(jī)變量向量，Y(t)是m維隨機(jī)過程向量，g是產(chǎn)品的性能函數(shù)。

因此，產(chǎn)品在壽命周期內(nèi)的失效概率為

P(0,T)=Pr{g(X,Y(t),t)≤0,?t∈[0,T]}

(2)

3 基于首次穿越時間概率模型的時變可靠性分析

穿越率方法已經(jīng)成為產(chǎn)品時變可靠性評估的主流方法，首先基于一些隨機(jī)過程的假設(shè)來計(jì)算穿越率，然后利用穿越率估計(jì)可靠度或者失效概率是穿越率法的一般流程。然而，由于穿越率計(jì)算時的隨機(jī)過程假設(shè)，導(dǎo)致了該方法目前僅能針對少數(shù)幾類特殊形式的性能函數(shù)給出有效的穿越率計(jì)算方法，而對于一般類型的性能函數(shù)，采用蒙特卡洛計(jì)算穿越率效率極低。為此，本文提出了一種基于首次穿越時間概率模型的新型時變可靠性分析方法(F-PTPD方法)。

3.1 基于首次穿越時間概率密度函數(shù)的轉(zhuǎn)化模型

為了簡化計(jì)算過程，輸入隨機(jī)過程向量Y(t)通過離散隨機(jī)過程方法[19]分解成隨機(jī)向量Y*和時間t的耦合函數(shù)。因此，性能函數(shù)變?yōu)殡S機(jī)變量向量和時間t的函數(shù)，不妨設(shè)為g(Z,t)，其中，Z=(X,Y*)。

當(dāng)賦予性能函數(shù)輸入隨機(jī)變量其自身分布的樣本值Z時，性能函數(shù)視為關(guān)于時間t的一元函數(shù)，圖1展示了樣本值Z和其相對應(yīng)的性能g(Z,t)的關(guān)系。

如圖1所示，g(Z,t)與時間軸的三個交點(diǎn)分別是t1，t2和t3，稱第一個交點(diǎn)t1為首次穿越點(diǎn)。首次穿越點(diǎn)與壽命周期[0,T]的位置關(guān)系決定了當(dāng)輸入為Z時產(chǎn)品的失效情況。如圖1(a)所示，t1>T時，產(chǎn)品可靠；反之，如圖1(b)所示，當(dāng)t1∈[0,T]時，產(chǎn)品失效。根據(jù)產(chǎn)品性能函數(shù)g(Z,t)可知，首次穿越點(diǎn)是輸入隨機(jī)變量的函數(shù)，記為t1=t(Z)。則產(chǎn)品在壽命周期[0,T]的失效概率可以表示為

P(0,T)=Pr{t1∈[0,T]}

(3)

對于首次穿越點(diǎn)函數(shù)t1=t(Z)，若已知其概率密度函數(shù)f(τ)，則式(3)所示的產(chǎn)品失效概率可以轉(zhuǎn)化為

(4)

3.2 首次穿越時間的概率密度函數(shù)

由于不同產(chǎn)品性能函數(shù)不同，而且復(fù)雜產(chǎn)品的性能函數(shù)維度和非線性程度較高，導(dǎo)致無法直接利用g(Z,t)=0來求解首次穿越點(diǎn)。因此，利用將性能函數(shù)在性能函數(shù)均值函數(shù)的首次穿越點(diǎn)處二階泰勒展開為關(guān)于時間的二次函數(shù)，利用二次函數(shù)的性質(zhì)來求解首次穿越點(diǎn)與輸入變量的函數(shù)關(guān)系；再利用稀疏網(wǎng)絡(luò)隨機(jī)配置方法[20]得到首次穿越點(diǎn)的四階原點(diǎn)矩；基于首次穿越點(diǎn)的四階原點(diǎn)矩，采用最大熵估計(jì)方法來估計(jì)首次穿越時間的概率密度函數(shù)。

3.2.1 首次穿越點(diǎn)函數(shù)的求解

根據(jù)隨機(jī)過程均值函數(shù)定義，性能函數(shù)g(Z,t)的均值函數(shù)可以表達(dá)為

(5)

圖1Z和其相對應(yīng)的性能函數(shù)g(Z,t)的關(guān)系

Fig.1 Geometrical relationship betweeng(Z,t)andZ

(6)

對于式(6)的多重積分模型，利用稀疏網(wǎng)絡(luò)隨機(jī)配置方法來執(zhí)行均值函數(shù)的計(jì)算。根據(jù)稀疏網(wǎng)絡(luò)隨機(jī)配置方法的特點(diǎn)，本文假設(shè)所有的輸入隨機(jī)變量均服從相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布，否則利用Rosenblatt變換[20]將其他分布類型轉(zhuǎn)化成為相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)分布。基于Smolyak算法[20]，利用稀疏網(wǎng)絡(luò)隨機(jī)配置來計(jì)算式(6)的多重積分：

(7)

(8)

式中q值取決于性能函數(shù)的非線性程度，在實(shí)際工程中為了保證計(jì)算效率和精度，取2≤q≤4。

當(dāng)輸入變量積分點(diǎn)產(chǎn)生后，式(7)就是關(guān)于時間t的函數(shù)，即性能函數(shù)的均值函數(shù)。因此，均值函數(shù)的首次穿越點(diǎn)可由μ(t)=0獲得，記為t*。將性能函數(shù)g(Z,t)在均值函數(shù)的首次穿越點(diǎn)t*處二階泰勒展開為

(9)

3.2.2 首次穿越點(diǎn)的概率密度函數(shù)

當(dāng)獲得首次穿越點(diǎn)關(guān)于輸入隨機(jī)變量的函數(shù)t1=t(Z)后，首次穿越點(diǎn)的四階原點(diǎn)矩可以表達(dá)為

(10)

利用稀疏網(wǎng)絡(luò)隨機(jī)配置來估計(jì)式(10)的多重積分，計(jì)算過程和均值函數(shù)計(jì)算過程相同。于是，得到首次穿越點(diǎn)的四階原點(diǎn)矩M1,…,M4?；谒碾A原點(diǎn)矩利用最大熵方法來估計(jì)首次穿越點(diǎn)的概率密度函數(shù)，首次穿越點(diǎn)的概率密度函數(shù)f(τ)可由最大熵模型(11)獲得。

(11)

(12)

3.3 基于首次穿越點(diǎn)概率模型的時變可靠性分析

由稀疏網(wǎng)絡(luò)隨機(jī)配置法和最大熵概率密度估計(jì)方法可以得到首次穿越點(diǎn)的概率密度函數(shù)，因此在壽命周期[0,T]的可靠度可以表示為

(13)

F-PTPD方法流程如圖2所示，具體步驟如下。

(1)隨機(jī)過程離散。根據(jù)輸入隨機(jī)過程信息，將輸入隨機(jī)過程向量Y(t)通過離散隨機(jī)過程方法轉(zhuǎn)化成隨機(jī)向量Y*和時間t的耦合函數(shù)。

(2)求解均值函數(shù)首次穿越點(diǎn)。采用稀疏網(wǎng)絡(luò)隨機(jī)配置方法求解性能函數(shù)均值函數(shù)μ(t),并通過求解方程μ(t)=0獲得均值函數(shù)首次穿越點(diǎn)t*。

(3)確定性能函數(shù)首次穿越點(diǎn)函數(shù)。通過對性能函數(shù)在t*二階泰勒展開，并利用一元二次函數(shù)求根方法，求得性能函數(shù)首次穿越點(diǎn)函數(shù)t1=t(Z)。

(4)獲得性能函數(shù)首次穿越點(diǎn)概率密度函數(shù)。利用稀疏網(wǎng)絡(luò)隨機(jī)配置方法結(jié)合最大熵概率密度估計(jì)方法，可獲得性能函數(shù)首次穿越點(diǎn)概率密度函數(shù)f(τ)。

(5)時變可靠性計(jì)算。對概率密度函數(shù)f(τ)積分，可得到壽命周期[0,T]的可靠度。

4 數(shù)值算例

算例1考慮性能函數(shù)(14)

(14)

圖2 F-PTPD方法的流程

Fig.2 Flowchart of F-PTPD method

式中X1和X2均服從正態(tài)分布，且X1～N(5,0.25),X2～N(5,0.25)。圖3為利用本文方法得到的首次穿越點(diǎn)的概率密度函數(shù)。圖4為由本文方法得到的概率密度函數(shù)通過對產(chǎn)品服役時間的積分求得產(chǎn)品時變失效概率的變化趨勢，同時也展示了作為參考的蒙特卡洛方法估算的時變失效概率?？梢钥闯?，兩種方法所得結(jié)果非常接近，說明由本文方法能夠精確地估計(jì)產(chǎn)品壽命周期的時變失效概率的變化趨勢。

此外，在利用隨機(jī)網(wǎng)絡(luò)系數(shù)配置法來估計(jì)性能函數(shù)均值函數(shù)和首次穿越點(diǎn)原點(diǎn)矩時，取q=3,于是可產(chǎn)生13個輸入樣本，同時算例1利用前四階原點(diǎn)矩和最大熵優(yōu)化模型來估計(jì)概率密度函數(shù)，求解最大熵優(yōu)化模型需要迭代約20次，因此利用本文方法來估計(jì)首次穿點(diǎn)概率密度函數(shù)的函數(shù)調(diào)用次數(shù)為13+13×20=273。在蒙特卡洛仿真過程中，輸入隨機(jī)變量的樣本個數(shù)為105，此外，在估計(jì)每段時間區(qū)間可靠度時，取時間間隔為0.01。如在利用蒙特卡洛仿真來估計(jì)時間區(qū)間[0,15]的時變可靠度時，需要函數(shù)調(diào)用次數(shù)為105×15×100=1.5×108。本文方法只需273次函數(shù)調(diào)用就可以得到產(chǎn)品壽命周期關(guān)于時間的概率密度函數(shù)，然后通過一次積分運(yùn)算就可以得到任意區(qū)間的失效概率，然而，在利用蒙特卡洛仿真計(jì)算時間區(qū)間[0,T]的失效概率時，需要107×T次函數(shù)調(diào)用?？梢钥闯?，本文方法能夠在保證精度的前提下極大地提高時變失效概率的評估效率。

圖3 算例1首次穿越點(diǎn)概率密度函數(shù)

Fig.3 Probability density of first-passage time point for the example 1

圖4 算例1的時變失效概率

Fig.4 Time-dependent failure probability for the example 1

算例2圖5為一懸臂梁結(jié)構(gòu)，其中，長度L=5 m，界面是寬為b0，高為h0的矩形，且b0和h0均為隨機(jī)變量。材料密度為ρ=7.85 kg/m3,屈服強(qiáng)度為σ。該簡支梁中點(diǎn)承受集中隨機(jī)載荷F(t)和均勻分布載荷p的作用，其中隨機(jī)載荷F(t)是均值為3500 N、變異系數(shù)為0.2、自相關(guān)函數(shù)為exp[-(3τ)2]的高斯過程；均勻分布載荷p可以表示為ρgb0h0。

該懸臂梁從t=0時刻開始腐蝕，且在整個截面上腐蝕程度各向同性，同時腐蝕部分完全失去機(jī)械強(qiáng)度，則未腐蝕部分的面積表示為

A(t)=b(t)×h(t)

(15)

式中b(t)=b0-2κt,h(t)=h0-2κt,κ=0.03 mm/年。當(dāng)應(yīng)力達(dá)到材料的極限應(yīng)力時，該簡支梁發(fā)生失效，于是，性能函數(shù)可以表達(dá)為

g(X,Y(t),t)=σb(t)h(t)2/4-[F(t)L/4+

ρb(t)h(t)L2/8]

(16)

該懸臂梁的隨機(jī)分布信息列入表1。圖6給出了由本文方法所估計(jì)得到的懸臂梁的概率密度函數(shù)，該懸臂梁在服役期間的失效概率變化趨勢如圖7所示，同時圖7包含了蒙特卡洛仿真結(jié)果。可以看出，兩種方法得到的結(jié)果非常接近，這說明本文方法能夠在提高計(jì)算效率的同時保證估計(jì)結(jié)果的精確性。

圖5 懸臂梁

Fig.5 Cantilever beam

表1 簡支梁的隨機(jī)參數(shù)分布

Tab.1 Distribution information of cantilever beam

參數(shù)分布類型均值變異系數(shù)b0正態(tài)分布0.2m0.05h0正態(tài)分布0.04m0.10σ正態(tài)分布180MPa0.10

圖6 算例2的首次穿越點(diǎn)概率密度函數(shù)

Fig.6 Probability density of first-passage time point for the example 2

圖7 算例2的時變失效概率

Fig.7 Time-dependent failure probability for the example 2

5 結(jié) 論

本文基于對產(chǎn)品性能函數(shù)的分解，利用稀疏網(wǎng)絡(luò)隨機(jī)配置方法和最大熵方法估計(jì)性能函數(shù)首次穿越點(diǎn)的概率密度函數(shù)，提出一種求解時變可靠性的方法，可獲得產(chǎn)品整個壽命周期時變失效概率的變化趨勢。通過對性能函數(shù)二階泰勒展開，將性能函數(shù)分解成關(guān)于時間的二次函數(shù)；利用二次方程求根公式得到性能函數(shù)首次穿越點(diǎn)關(guān)于輸入隨機(jī)變量的函數(shù)表達(dá)式；采用稀疏網(wǎng)絡(luò)隨機(jī)配置方法和最大熵方法估算首次穿越點(diǎn)的概率密度函數(shù)。與傳統(tǒng)的時變可靠性評估方法相比，本文方法可獲得產(chǎn)品整個壽命周期失效概率的變化趨勢，而不是某個時間區(qū)間，極大地提高了評估效率，對復(fù)雜產(chǎn)品的可靠性評估設(shè)計(jì)有一定的工程指導(dǎo)意義。同時兩個算例表明，該方法的精度滿足工程需求。