王博, 蔣平, 郭波
(國(guó)防科技大學(xué) 系統(tǒng)工程學(xué)院, 湖南 長(zhǎng)沙 410073)
當(dāng)前,國(guó)家對(duì)國(guó)防軍工裝備提出更高層次的要求,對(duì)重點(diǎn)關(guān)注的航天產(chǎn)品也加大了投入與重視。航天閥門作為液體火箭發(fā)動(dòng)機(jī)的關(guān)鍵部件,要具備很高的可靠性水平,在研制階段就需準(zhǔn)確評(píng)估其可靠性。航天閥門的工作原理與結(jié)構(gòu)復(fù)雜,評(píng)估難度大,同時(shí),由于小樣本、高可靠性的特點(diǎn),在研制階段試驗(yàn)中僅能得到極少失效甚至無失效的試驗(yàn)數(shù)據(jù)。
航天閥門的失效模式和影響分析已經(jīng)趨向成熟,針對(duì)閥門啟動(dòng)(關(guān)閉)失靈、卡滯、閥門污染(泄漏)等故障模式的狀態(tài)識(shí)別技術(shù)、失效機(jī)理分析以及設(shè)計(jì)改進(jìn),已有大量研究成果。但是,在發(fā)動(dòng)機(jī)的整機(jī)試車時(shí),仍然會(huì)出現(xiàn)卡滯、泄漏等問題,這說明現(xiàn)有研究仍未充分考察實(shí)際工作環(huán)境,也難以準(zhǔn)確評(píng)估閥門的可靠性。液體火箭發(fā)動(dòng)機(jī)工作環(huán)境惡劣,閥門可靠性要求閥門在工況溫度和振動(dòng)環(huán)境下,額定的時(shí)間內(nèi)完成規(guī)定的開閉次數(shù)。因此航天閥門在定型之前要進(jìn)行溫度、振動(dòng)和動(dòng)作等研制試驗(yàn),保證閥門在真實(shí)工作環(huán)境下的可靠性。在現(xiàn)有閥門可靠性評(píng)估成果中,這些研制試驗(yàn)數(shù)據(jù)未得到充分利用。
本文提出融合多類型研制試驗(yàn)數(shù)據(jù)的方法來評(píng)估航天閥門的可靠性,對(duì)航天閥門的溫度、振動(dòng)和動(dòng)作試驗(yàn)數(shù)據(jù)進(jìn)行預(yù)處理,融合成統(tǒng)一的先驗(yàn)分布,參與閥門多源信息融合的可靠性評(píng)估過程。該方法利用Bayes理論,對(duì)現(xiàn)場(chǎng)數(shù)據(jù)的評(píng)估結(jié)果進(jìn)行修正,減少對(duì)試驗(yàn)樣本量的需求。
機(jī)械中的疲勞強(qiáng)度、疲勞壽命、磨損壽命、腐蝕壽命大多服從Weibull分布。根據(jù)工程經(jīng)驗(yàn),閥門的可靠工作次數(shù)為壽命型參數(shù),其可靠性分布規(guī)律符合Weibull分布。因此,本文假設(shè)閥門服從兩參數(shù)的Weibull分布——(,),其中為形狀參數(shù),為尺度參數(shù)(也稱特征壽命)。航天閥門屬于高可靠性產(chǎn)品,受到研制時(shí)間和研制經(jīng)費(fèi)的限制,研制試驗(yàn)常常采用截尾試驗(yàn),試驗(yàn)終止時(shí)僅出現(xiàn)極少失效甚至無失效。
閥門的溫度試驗(yàn)采用不等定時(shí)截尾的方式,在試驗(yàn)進(jìn)行到截尾時(shí)間t時(shí),進(jìn)行第次截尾,去掉t個(gè)樣本。已知投入試驗(yàn)的8個(gè)樣本在截尾時(shí)均未失效,截尾時(shí)間與對(duì)應(yīng)的無失效樣本量t如表1所示。
表1 溫度試驗(yàn)的不等定時(shí)截尾數(shù)據(jù)Tab.1 Random censored data from temperature test
一組不等定時(shí)截尾數(shù)據(jù)可以看作多組定時(shí)截尾數(shù)據(jù)。采用傳統(tǒng)統(tǒng)計(jì)方法中的極大似然估計(jì)(MLE)法進(jìn)行可靠性評(píng)估時(shí),若樣本數(shù)目足夠多,收斂性質(zhì)會(huì)表現(xiàn)良好,且此方法簡(jiǎn)單、結(jié)果無偏或漸近無偏。但是,MLE法并不適用于無失效的情況。王玲玲等提出修正MLE(MMLE)法,引入修正常數(shù)(0<<1)作為單組無失效數(shù)據(jù)的偽失效數(shù),再采用MLE方法來開展無失效數(shù)據(jù)的可靠性評(píng)估方法。因此,溫度試驗(yàn)數(shù)據(jù)的似然函數(shù)為
(1)
式中:和為溫度試驗(yàn)數(shù)據(jù)的Weibull參數(shù);為定時(shí)截尾組數(shù)。
形狀參數(shù)和特征壽命的MMLE法估計(jì)結(jié)果如(2)式所示:
(2)
(3)
顯著性水平可以依據(jù)可靠性要求取05、04或03,則特征壽命的評(píng)估結(jié)果為
(4)
根據(jù)相似產(chǎn)品或歷史數(shù)據(jù),結(jié)合專家經(jīng)驗(yàn)給出形狀參數(shù)的取值范圍,通過數(shù)值迭代法確定(2)式中和
振動(dòng)試驗(yàn)以服從高斯分布的隨機(jī)振動(dòng)形式進(jìn)行,樣本在不同的功率譜密度下,分別進(jìn)行定時(shí)截尾試驗(yàn)。將6個(gè)受試樣本分3組進(jìn)行恒應(yīng)力加速壽命試驗(yàn),試驗(yàn)過程中無失效產(chǎn)生,試驗(yàn)結(jié)果如表2所示。已知工況條件下,功率譜密度為2.5/Hz.
在隨機(jī)振動(dòng)下,閥門所受應(yīng)力與壽命之間的關(guān)系滿足逆冪律模型:
=-,
(5)
式中:為振動(dòng)應(yīng)力下的閥門壽命;為系統(tǒng)黏性阻尼常數(shù);為振動(dòng)應(yīng)力水平;為材料結(jié)構(gòu)常數(shù)。
表2 振動(dòng)試驗(yàn)加速應(yīng)力壽命試驗(yàn)數(shù)據(jù)Tab.2 Accelerated stress life test data from vibration test
(5)式兩邊取對(duì)數(shù),得
ln=0+1ln,
(6)
式中:0和1分別為加速方程系數(shù)、常數(shù)。
定義環(huán)境C對(duì)環(huán)境B(工況環(huán)境)的環(huán)境因子為
(7)
式中:為產(chǎn)品在環(huán)境B即工況環(huán)境下的壽命;為產(chǎn)品在環(huán)境C下的壽命。
由(6)式聯(lián)合(7)式推導(dǎo)得到加速方程為逆冪律模型下的環(huán)境因子:
=exp (1(ln-ln)),
(8)
式中:、分別為產(chǎn)品在環(huán)境B和環(huán)境C下的振動(dòng)應(yīng)力水平。
根據(jù)(8)式將表2中的恒應(yīng)力加速試驗(yàn)數(shù)據(jù)換算為工況應(yīng)力水平下的數(shù)據(jù),進(jìn)一步地,參考溫度試驗(yàn)數(shù)據(jù)的MMLE法評(píng)估完成振動(dòng)試驗(yàn)數(shù)據(jù)的可靠性評(píng)估。
動(dòng)作試驗(yàn)同樣采用不等定時(shí)截尾的方式,在進(jìn)行a次開閉后,第次去掉a個(gè)樣本。已知投入試驗(yàn)的10個(gè)樣本在截尾時(shí)均無失效,截尾次數(shù)a與無失效樣本量a如表3所示。
表3 動(dòng)作試驗(yàn)的不等定時(shí)截尾數(shù)據(jù)Tab.3 Random censored data from action test
閥門正常開閉次數(shù)是衡量閥門可靠性的重要指標(biāo)。本文將開閉次數(shù)處理為連續(xù)型數(shù)據(jù),以作為現(xiàn)場(chǎng)數(shù)據(jù)參與多源信息融合的可靠性評(píng)估過程。
參考文獻(xiàn)[13]的Bayes可靠性評(píng)估方法,首先給出單獨(dú)對(duì)動(dòng)作試驗(yàn)數(shù)據(jù)進(jìn)行可靠性評(píng)估的結(jié)果。不等定時(shí)截尾數(shù)據(jù)同樣可以看作多組定時(shí)截尾數(shù)據(jù)。假設(shè):形狀參數(shù)服從均勻分布(,),、分別為形狀參數(shù)的上、下限;任務(wù)可靠度服從均勻分布(,1),為任務(wù)可靠度的下限值。任務(wù)工作次數(shù)為,則動(dòng)作試驗(yàn)數(shù)據(jù)的似然函數(shù)為
(9)
式中:和為動(dòng)作試驗(yàn)數(shù)據(jù)的Weibull參數(shù);為定時(shí)截尾試驗(yàn)的組數(shù);(a)為動(dòng)作試驗(yàn)數(shù)據(jù)下的可靠度函數(shù)。
根據(jù)表3中的數(shù)據(jù),定時(shí)截尾試驗(yàn)組數(shù)=4,在任務(wù)工作次數(shù)時(shí)的可靠度函數(shù)為
(10)
變換可得
(11)
(11)式代入似然函數(shù),由Bayes定理得到形狀參數(shù)和任務(wù)可靠度的聯(lián)合驗(yàn)后分布,并提取核函數(shù):
(12)
馬爾可夫鏈蒙特卡洛(MCMC)抽樣方法的基礎(chǔ)理論是馬爾可夫過程。在MCMC方法中,為在一個(gè)指定的分布上抽樣,根據(jù)馬爾可夫過程,首先從任一狀態(tài)出發(fā),模擬馬爾可夫過程,不斷進(jìn)行狀態(tài)轉(zhuǎn)移,最終收斂到平穩(wěn)分布。Metropolis-Hastings(M-H)抽樣算法和Gibbs抽樣算法是MCMC方法中使用較為廣泛的兩種形式。M-H抽樣算法不要求已知分布是對(duì)稱的,且算法所構(gòu)造的馬爾可夫過程滿足細(xì)致平穩(wěn)條件。因此考慮使用M-H抽樣算法計(jì)算Bayes驗(yàn)后分布。
航天閥門具備所有航天產(chǎn)品的共性,可投入試驗(yàn)的樣本量很小、試驗(yàn)經(jīng)常出現(xiàn)無失效。充分利用研制試驗(yàn)數(shù)據(jù),可在一定程度上增加可靠性評(píng)估的準(zhǔn)確性。在融合壽命相關(guān)數(shù)據(jù)的場(chǎng)合,Bayes方法應(yīng)用最多, 因?yàn)锽ayes方法結(jié)合現(xiàn)場(chǎng)數(shù)據(jù)與各類先驗(yàn)信息,可以到達(dá)修正參數(shù)估計(jì)結(jié)果的效果,同時(shí)需要的樣本量相對(duì)MLE方法來說也更少。將不同類型試驗(yàn)的數(shù)據(jù)轉(zhuǎn)化為先驗(yàn)分布,與現(xiàn)場(chǎng)數(shù)據(jù)進(jìn)行一致性檢驗(yàn),融合成統(tǒng)一的先驗(yàn)分布,可以有效提高閥門Bayes可靠性評(píng)估的精度。
試驗(yàn)數(shù)據(jù)的預(yù)處理工作主要是對(duì)數(shù)據(jù)的可靠性評(píng)估。將溫度和動(dòng)作試驗(yàn)數(shù)據(jù)的預(yù)處理結(jié)果,即可靠度,作為先驗(yàn)信息,可進(jìn)一步轉(zhuǎn)化為先驗(yàn)分布。
已知失效率的共軛先驗(yàn)分布為Gamma分布,根據(jù)變換下的不變性原則可得,閥門任務(wù)工作時(shí)間(min)時(shí)的可靠度先驗(yàn)分布為負(fù)對(duì)數(shù)Gamma分布。取任務(wù)可靠度的先驗(yàn)分布為
(13)
先驗(yàn)分布(13)式中存在兩個(gè)超參數(shù)和,且>0,>0本文采用最大熵法利用可靠度點(diǎn)估計(jì)值確定多源信息的先驗(yàn)分布參數(shù)。
若已知在任務(wù)時(shí)間時(shí)的可靠度點(diǎn)估計(jì)值為,使用最大熵法(14)式可以確定負(fù)對(duì)數(shù)Gamma分布(13)式中的超參值:
(14)
進(jìn)一步簡(jiǎn)化(14)式,可得
(15)
這樣,問題就被轉(zhuǎn)化為一維規(guī)劃問題,有利于未知參數(shù)、的求解。
使用Bayes方法進(jìn)行信息融合,首先要進(jìn)行一致性檢驗(yàn),確保所有先驗(yàn)信息與現(xiàn)場(chǎng)數(shù)據(jù)是屬于同一整體的。 秩和檢驗(yàn)、Mood檢驗(yàn)和假設(shè)檢驗(yàn)都可以直接利用樣本數(shù)據(jù)進(jìn)行檢驗(yàn),但是要求有一定的失效數(shù),這些方法對(duì)無失效數(shù)據(jù)是行不通的。因此,考慮Bayes置信區(qū)間法進(jìn)行檢驗(yàn),由單個(gè)信息源結(jié)合現(xiàn)場(chǎng)數(shù)據(jù)得到參數(shù)的置信區(qū)間,如果現(xiàn)場(chǎng)數(shù)據(jù)在無信息先驗(yàn)下的參數(shù)點(diǎn)估計(jì)值落在置信區(qū)間內(nèi),則認(rèn)為該信息源與現(xiàn)場(chǎng)數(shù)據(jù)是相容的,即通過一致性檢驗(yàn)。
(16)
為融合多源信息,需要先將可靠度的先驗(yàn)分布(13)式轉(zhuǎn)化為失效率的先驗(yàn)分布,可靠度函數(shù)為
=exp (-)
(17)
對(duì)可靠度函數(shù)關(guān)于失效率求導(dǎo),可得
(18)
(18)式中,導(dǎo)函數(shù)恒小于0,因此可靠度函數(shù)是關(guān)于的減函數(shù),有
(19)
等式兩邊關(guān)于求導(dǎo),可得
(20)
假設(shè)形狀參數(shù)服從均勻分布(,),其中,參數(shù)上、下限和根據(jù)相似產(chǎn)品、歷史數(shù)據(jù)由專家經(jīng)驗(yàn)給出,則形狀參數(shù)與失效率的聯(lián)合先驗(yàn)分布為
(21)
21節(jié)解出的超參數(shù)和代入(21)式,分別得到溫度和振動(dòng)試驗(yàn)信息所形成的先驗(yàn)分布(,)和(,),則融合后的先驗(yàn)分布為
(,)=(,)+(,),
(22)
式中:、分別為溫度試驗(yàn)、振動(dòng)試驗(yàn)的先驗(yàn)信息在融合中所占的權(quán)重。根據(jù)第Ⅱ類極大似然估計(jì)(ML-Ⅱ)法可以確定權(quán)重值。
假設(shè)某現(xiàn)場(chǎng)數(shù)據(jù)樣本服從分布(|),為未知參數(shù)。已知多源試驗(yàn)信息形成的先驗(yàn)分布為ρ(),其中=1,2,…,由先驗(yàn)分布ρ()得到的邊緣分布為
(23)
將(|)看作是的函數(shù)(極大似然函數(shù)),它與的先驗(yàn)分布無關(guān),只反映了樣本的信息。ML-Ⅱ法實(shí)質(zhì)上就是根據(jù)在不同的先驗(yàn)分布下,現(xiàn)場(chǎng)樣本出現(xiàn)的似然性大小來確定不同的融合權(quán)重。
假設(shè)現(xiàn)場(chǎng)數(shù)據(jù),,…,,將其看作是由邊緣分布(,ρ)產(chǎn)生的,得到類試驗(yàn)先驗(yàn)分布下的現(xiàn)場(chǎng)數(shù)據(jù)似然函數(shù)為
(24)
則權(quán)重的計(jì)算方法為
(25)
動(dòng)作試驗(yàn)數(shù)據(jù)作為現(xiàn)場(chǎng)數(shù)據(jù),則閥門的極大似然函數(shù)為
(|,)=∏(a;,)=exp (-),
(26)
依據(jù)貝葉斯原理,融合后的驗(yàn)后分布為
(27)
驗(yàn)后分布計(jì)算使用MCMC抽樣方法中的M-H算法,采用MATLAB軟件實(shí)現(xiàn)的抽樣過程如圖1所示。
已得形狀參數(shù)與失效率的點(diǎn)估計(jì)值,由(28)式計(jì)算得到Weibull分布中特征壽命的估計(jì)值,至此完成航天閥門多源信息融合的可靠性評(píng)估。
(28)
溫度試驗(yàn)數(shù)據(jù)依據(jù)11節(jié)方法,給定形狀參數(shù)的區(qū)間范圍為[1,6],數(shù)值迭代法得到溫度試驗(yàn)下的閥門Weibull分布參數(shù):=327,=70724 6
根據(jù)表2中的加速應(yīng)力(振動(dòng))試驗(yàn)數(shù)據(jù)擬合加速模型,獲得加速模型系數(shù):
ln=4470 3-0484 8ln,
(29)
使用環(huán)境因子將加速應(yīng)力數(shù)據(jù)折算到工況水平,結(jié)果如表4所示。
表4 振動(dòng)試驗(yàn)數(shù)據(jù)折算結(jié)果Tab.4 Conversion results of vibration test data
折算后的數(shù)據(jù)仍具備無失效特性,MMLE方法進(jìn)行參數(shù)估計(jì)得到振動(dòng)試驗(yàn)下的閥門Weibull分布參數(shù):形狀參數(shù)=3,特征壽命=154.662 1.
單獨(dú)對(duì)動(dòng)作試驗(yàn)無失效數(shù)據(jù)進(jìn)行可靠性評(píng)估時(shí),采用Bayes方法,給定形狀參數(shù)取值范圍為 [2,3],任務(wù)工作次數(shù)為20次,任務(wù)可靠度下限為0.8,得到動(dòng)作試驗(yàn)下的閥門Weibull分布參數(shù):=2.284 7,=831.818 9.
假設(shè)閥門的任務(wù)工作時(shí)間為10 min,動(dòng)作試驗(yàn)任務(wù)開閉次數(shù)為20次。給定的形狀參數(shù)取值范圍為[2,3]。根據(jù)以上可靠性評(píng)估結(jié)果可以得到各類試驗(yàn)的先驗(yàn)分布參數(shù),如表5所示。
3類試驗(yàn)形成的先驗(yàn)信息已知,對(duì)其進(jìn)行一致 性檢驗(yàn),顯著性水平取0.1,計(jì)算得到溫度、振動(dòng)試驗(yàn)數(shù)據(jù)下的可靠度區(qū)間估計(jì)和動(dòng)作試驗(yàn)作為現(xiàn)場(chǎng)數(shù)據(jù)的無信息先驗(yàn)點(diǎn)估計(jì)值,如表6所示。
表5 3類試驗(yàn)的先驗(yàn)分布參數(shù)Tab.5 Prior distribution parameters for three types of tests
表6 一致性檢驗(yàn)計(jì)算結(jié)果Tab.6 Calculated results of consistency test
現(xiàn)場(chǎng)數(shù)據(jù)的點(diǎn)估計(jì)值落入?yún)^(qū)間估計(jì)范圍內(nèi),滿足Bayes置信區(qū)間法對(duì)多源信息通過一致性檢驗(yàn)的要求。因此,3類試驗(yàn)信息可以進(jìn)行多源信息融合。結(jié)合表1與表4中的數(shù)據(jù),依據(jù)(24)式和(25)式計(jì)算得到融合權(quán)重為==0.5.
經(jīng)推導(dǎo),驗(yàn)后分布的核函數(shù)為
(|)∝()(|)=
((,)+(,))exp (-)
(30)
通過MCMC抽樣方法得到閥門分布參數(shù)的點(diǎn)估計(jì)值為=1595 355×10、=255,抽樣迭代結(jié)果頻數(shù)分布分別如圖2和圖3所示。將抽樣得到的參數(shù)估計(jì)值代入(28)式,得特征壽命=466007 7
圖2 MCMC抽樣失效率λ頻數(shù)分布圖Fig.2 Frequency distribution of failure rate λ by MCMC sampling method
圖3 MCMC抽樣形狀參數(shù)m頻數(shù)分布圖Fig.3 Frequency distribution of shape parameter mby MCMC sampling method
為了更加直觀地展示可靠性評(píng)估的結(jié)果,將3類試驗(yàn)數(shù)據(jù)單組的可靠性評(píng)估結(jié)果與融合后的可靠性評(píng)估結(jié)果歸納于表7,并給出任務(wù)可靠度,在圖4中畫出對(duì)應(yīng)的可靠度函數(shù)曲線。
表7 不同類型試驗(yàn)數(shù)據(jù)和多源信息融合的 可靠性評(píng)估結(jié)果Tab.7 Reliability evaluation results of different types of test data and multi-source information fusion
圖4 可靠度函數(shù)曲線Fig.4 Reliability curves
已知某型號(hào)閥門進(jìn)行3個(gè)批次的生產(chǎn),共組裝37個(gè)閥門。這些閥門經(jīng)過了典型試驗(yàn),包括常溫性能檢查試驗(yàn)、高低溫性能檢查試驗(yàn)、振動(dòng)后性能檢查試驗(yàn)和動(dòng)作壽命試驗(yàn),又進(jìn)行了延壽試驗(yàn)。在系列試驗(yàn)中共出現(xiàn)3個(gè)失效,失效時(shí)的累積動(dòng)作次數(shù)分別為2 050、5 450和4 400. 為了簡(jiǎn)化計(jì)算過程,假設(shè)該批次閥門的壽命服從指數(shù)分布。通過MLE法計(jì)算任務(wù)動(dòng)作次數(shù)為20次時(shí)的可靠度為0.995.
本文方法計(jì)算出的任務(wù)可靠度0.999 7與大樣本極大似然的計(jì)算結(jié)果0.995相比較,誤差在0.5%以內(nèi)。此外,在第2 050次動(dòng)作試驗(yàn)失效的閥門實(shí)際上經(jīng)歷了1 000次的高溫動(dòng)作壽命試驗(yàn),因此,閥門的真實(shí)任務(wù)可靠度大于0.995. 閥門設(shè)計(jì)要求為0.999 83,專家一致認(rèn)為該型號(hào)閥門達(dá)到任務(wù)可靠度指標(biāo)。通過以上對(duì)比可以有效說明本文評(píng)估方法的準(zhǔn)確性和精確度。
為了進(jìn)一步驗(yàn)證本文方法的穩(wěn)定性,計(jì)算機(jī)仿真運(yùn)行10 000次,計(jì)算任務(wù)可靠度的平均值和平均相對(duì)誤差,結(jié)果如表8所示。
表8 方法穩(wěn)定性驗(yàn)證結(jié)果Tab.8 Stability verification results of the proposed method
由表8可知,平均相對(duì)誤差相較可靠度平均值而言非常小,說明本文方法具有較好的穩(wěn)定性。
在航天閥門研制過程中,定量的可靠性評(píng)估還缺乏理論方法的支持。航天閥門投入可靠性試驗(yàn)的樣本量小,而且屬于高可靠性產(chǎn)品,在試驗(yàn)過程中往往僅能得到有限的無失效試驗(yàn)數(shù)據(jù)。
針對(duì)航天閥門不同試驗(yàn)數(shù)據(jù)的特點(diǎn),需要采用合適的方法對(duì)數(shù)據(jù)分別進(jìn)行預(yù)處理。將溫度試驗(yàn)和振動(dòng)試驗(yàn)數(shù)據(jù)信息融合為統(tǒng)一的先驗(yàn)分布,動(dòng)作數(shù)據(jù)作為現(xiàn)場(chǎng)數(shù)據(jù),基于Bayes原理完成對(duì)航天閥門的可靠性評(píng)估?;诙嘣葱畔⑷诤系暮教扉y門可靠性評(píng)估方法,既解決了小子樣、無失效數(shù)據(jù)的可靠性評(píng)估問題,又充分利用已有數(shù)據(jù),提高了閥門可靠性評(píng)估的精度和可信度。