基于Copula函數的多功能機械系統競爭失效分析方法

劉敬一, 莊新臣, 張玉剛, 龐歡

1.長安大學 汽車學院, 陜西 西安 710064; 2.清華大學 機械工程學院, 北京 100091;3.西北工業大學 航空學院, 陜西 西安 710072

隨著機械設計技術的發展和加工制造水平的提高，機械產品通常具有多個功能，與之相關的失效形式也更為多樣化，尤其對應用于航空航天、兵器裝備和重型運載工具等領域的機械系統而言，惡劣的工作條件和極端的載荷環境更易導致多失效模式的發生。對于多失效模式的機械系統，其各個失效模式之間往往是相互作用的，這種相互作用體現于：①失效相關，多個失效模式之間并非獨立，而是因共用組件、同受載荷等原因彼此作用；②相互競爭，某一種失效模式的出現會引起機構工作狀態的變化甚至導致機構系統停止工作，進而造成其他失效模式不會再有機會在同等條件下發生[1]。準確描述多失效模式之間相關且競爭的相互作用關系，成為對機械產品進行可靠性評估、壽命預測和維護策略制定的關鍵。

關于相關競爭失效的理論和建模方法方面，已有大量的文獻開展了相關研究。羅格斯大學的Pham教授較早地把競爭失效的概念推廣應用于可靠性工程領域，提出了可應用于工程上的競爭失效模型，認為對于受突發失效和退化失效影響的系統，任一種失效模式的出現都會導致系統的破壞，從而使另一種失效模式沒有機會再次發生，也就是說2種失效模式以相互競爭的方式綜合決定了系統的可靠性[2-5]。隨后羅格斯大學Coit教授分別考慮變失效閾值、系統多部件、沖擊分類和變退化速率等[6-9]因素，對相關競爭失效的建模方法進行了擴展，以適用于不同工程背景機械系統的研究。同時，實際工程應用中機械系統的相關競爭失效研究也逐漸受到關注。航空航天方面，李偉、王華偉、郭慶和劉曉娟等[10-14]分析了航空發動機的競爭失效，針對可靠性分析、可靠性評估和壽命預測等內容開展研究；車輛工程方面，閆書法等[15]對某型號重型車輛的傳動機構進行了競爭失效的建模，并利用試驗對所建立的模型進行驗證和應用。Yu等[16]提出可用于分析多態退化可修系統的可靠性模型，考慮系統狀態可退化的情況對其可靠性表現進行評估，并把相應模型應用于汽車車架結構的疲勞失效分析；船舶及水下運載工具領域，Bocchetti等[17]考慮磨損失效和熱應力失效，研究了某型輪船燃氣輪機的競爭失效現象，并基于此對該機械系統進行了可靠性評估，鄢偉安[18]則研究了水下魚雷翼面展開機構的競爭失效；兵器裝備領域，王浩偉、蔡忠義、孫麗和逯程等分別研究了導彈貯存[19]、彈翼機構[20]、起爆電容[21]和大型裝備[22]的競爭失效現象及維護保障方法。

綜合現有研究可知，當前研究更多地關注因載荷過大所導致的結構破壞失效模式，卻很少考慮“結構并未破壞但因退化因素累積導致功能無法完成”所造成的機械系統功能失效。在如今的材料科學技術和加工制造水平下，機械系統發生結構破壞失效的比例正逐步降低，更多的失效形式體現于運動卡滯、運動精度不足或構件磨損等形式的退化型功能失效；同時，當前研究往往假設失效模式相互獨立或使用相關系數對失效模式之間的相關性進行描述，獨立假設并不符合工程實際，而通過相關系數無法獲取相關情況下多失效模式的表征量的聯合概率密度函數，難以進行可靠性的計算求解。因此，本文建立了面向多功能機械系統的建立競爭失效分析模型，針對處于突發型結構破壞失效和多個退化型功能失效綜合作用下的機械系統，研究各個失效模式的表征方法和退化規律；基于Copula函數定量化失效模式之間的相關性，獲取多失效模式功能表征量的聯合概率密度函數及演化規律，形成可靠性評估方法；最后利用本文方法針對某型飛機艙門鎖機械系統進行了競爭失效分析，證明了本文所提方法的正確性和工程實用性。

1 機械系統多功能失效描述

功能失效指的是“機械系統在部件尚未發生明顯破壞時因功能無法完成而導致的失效”，機械產品的特定功能可以通過一個正確描述功能現象的指標進行表達，該指標定義為“功能表征量”。常見的功能表征量包括運動精度類、驅動力/阻力類和運動時間類等[23]。隨著機械系統本身的復雜化和綜合化，通常一個機械產品具有多個功能。對于大型機械系統而言，具有精密度高、作用載荷大、功能復雜等特點，需要進行合理的可靠性建模分析。機械產品的失效模式可以概括為突發型結構破壞失效和功能無法完成導致的退化型功能失效。

理想情況下機械系統的各個部件有足夠的能力承受工作載荷。但是由于隨機性的存在，工作載荷本身并非定值，而是一個隨機變量。且在復雜工作環境下的機械系統，一些偶發因素也會引發額外載荷，例如外界物體對機械系統的碰撞，機械系統本身的振動等，都會造成量值較大的載荷作用于各個部件，這種明顯大于工作載荷的偶發載荷，可看作是一種沖擊，當沖擊的量值超過部件的承受能力時，部件發生斷裂，進而引起機械系統突發失效。另一方面，部件參數的退化會引起功能表征量的退化，進而發生功能失效，例如磨損導致摩擦阻力增大而卡滯，運動副間隙過大引起運動精度不足等。機械系統的功能表征量具有明顯的退化特性，且各個功能指標之間具有相關性，需要綜合退化過程和相關性對機械系統的功能失效進行分析。

圖1 多功能機械系統競爭失效分析流程

任一失效模式的發生都會導致機械系統工作特性發生改變，進而其他失效模式不再有機會于同等條件下發生，即各個失效模式處于相互競爭的關系。對處于競爭失效模式下的機械系統進行可靠性分析時，首先需要對機械系統的運動特性進行分析，確定失效模式的具體表現形式(含突發型失效和退化型功能失效)；基于載荷隨機特性和機械系統部件抵抗外界沖擊的能力，確定所承載荷分布和閾值，進行突發型失效建模；根據機械系統各個功能屬性定義功能表征量，基于運動幾何學、運動學和動力學確定功能表征量及其影響因素之間的數學函數關系。對于難以直接表征的功能函數關系，可利用合適的代理模型構造每個功能表征量的功能函數，得到每個功能的累積分布函數；因功能之間的相關性，引入Copula函數構造多功能失效聯合概率分布函數，進而得到考慮失效相關性的多失效模式機械系統的可靠度或失效概率。基于Copula函數的多功能機械系統競爭失效分析的完整流程如圖1所示。

2 多失效模式機械系統失效建模方法

2.1 突發型失效建模方法

機械系統工作的過程中會承受相應的工作載荷，根據對多功能機械系統失效的描述，一些偶然因素會引起量值較大的沖擊載荷，如重型車輛在惡劣路面行駛過程中的意外沖擊，飛機飛行過程中結構顫振導致的沖擊載荷等，從而造成部件的破壞失效。根據定義，這種沖擊載荷在到來時間和每次載荷的量值大小上都具有不確定性。對于每一次突發沖擊載荷而言，一旦超過機械系統構件所能承受的最大閾值，系統就會因構件破壞而失效，故機械系統不發生突發失效的概率表示為公式(1)。

式中：Wi為第i次突發載荷的大小；D0為系統可以承受最大載荷的閾值。突發載荷的到來時刻符合速率為λ的泊松過程，有如下公式成立

式中，N(t)和N(t+τ)分別表示截止時刻t和時刻t+τ突發載荷的發生次數。當t為0時，根據突發載荷的物理意義有N(t)=0，則

即

考慮到時刻t為止外界沖擊載荷次數的所有可能情形，機構不發生突發失效的概率表示為公式(3)。

2.2 退化型功能失效建模方法

機構系統的特定功能可以使用具體的功能表征量進行定量化表征，根據指標判定機構功能是否正常。系統不發生功能失效的概率表示為公式(4)。

R(t)=P(|FS-F0|<δ)(4)

式中：F0為功能指標的理想值；FS為功能指標的實際取值；δ為指標理想值和實際值之間的允許誤差。上述表達式的意義就是如果實際功能指標和功能指標的理想值之間的誤差過大，則機械系統發生功能失效。

機械功能的實現受到多種因素的影響，包括各個部件特性、載荷、連接關系等。功能是否可以完成受到各個部件的共同作用影響。由于各部件在工作中伴隨著不可避免的退化作用(例如磨損、老化、腐蝕等)，進而導致機構的功能也會發生退化。功能表征量與其影響因素存在函數關系，如公式(5)所示。

FS=f(C,L)(5)

式中：C表示機構部件特性和連接關系；L為載荷。函數表達式f的具體形式，需要根據具體機械系統的實際功能相應確定。機械系統的某些參數會隨著工作循環次數和時間發生無法避免的退化(如磨損、老化)，進而通過功能函數，把這種退化傳遞給功能表征量，因此需要對參數的退化進行建模分析。漂移布朗運動是一種典型的用來對退化過程建模的方法。典型的漂移布朗運動的表達式為

X(t)=σB(t)+μt+X0(6)

式中：X(t)表示t時刻的退化量；σ為波動參數；B(t)為標準布朗運動；μ為漂移參數；X0表示退化量的初始值。通常情況下機械系統的退化是加速的，例如磨損的加劇會導致磨損接觸面變得更加粗糙，增大接觸面之間的摩擦力，從而進一步加速磨損；疲勞裂紋的出現使裂紋附近出現應力集中，進而加速裂紋的擴展。為了表示這種加速現象，在典型的漂移布朗運動中引入指數時間項，即把時間項t變為tq，其中q表示加速系數。此時漂移布朗運動的表達式變為公式(7)。

X(t)=σB(t)+μtq+X0(7)

當參數σ,μ,q和退化量的初始值X0為常量時，退化量X(t)～N(μtq+X0,σ2)。利用布朗運動對退化運動建模有幾點優勢[24]：①漂移參數的大小代表著退化過程中整體速率的高低，而引入指數的時間項則有效描述了工程中機械產品退化過程中的加速現象；②波動參數和標準布朗運動的乘積項代表著退化過程中退化量的隨機性，即相同時間段內的退化量并非常值，而是一個有著波動向量的不確定值；③正態分布具有比較高效的計算效率和精確的計算值。

3 基于Copula函數的多功能機械系統可靠度計算方法

對于多功能機械系統而言，分析其可靠性的核心在于建立多功能表征量的聯合分布函數(CDF)，但是復雜的相關關系，使得常規構造CDF的方法存在難度。而Copula函數則可以很好地描述相關性，并建立多變量的聯合分布函數。

“Copula”一詞起源于拉丁詞匯“link”和“tie”，為一組邊緣分布符合標準均勻分布變量的聯合分布函數，表示為公式(8)。

C(u1,u2,…,un)=P(U1

(8)

式中，Ui為標準正態分布,即Ui～U(0,1),i=1,2,…,n。根據Skler定理[25],邊緣分布為F1(x1),F2(x2),…,Fn(xn)的一維分布變量x1,x2,…,xn的n維聯合分布函數F表示為公式(9)。

F(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn))

(9)

式中，C為對應的Copula函數。如果邊緣分布F1(x1),F2(x2),…,Fn(xn)都是連續函數,則存在唯一Copula函數C滿足公式(10)。

(10)

式中，u1=F1(x1),u2=F2(x2),…,un=Fn(xn)。

阿基米德Copula作為一類重要的函數族,因構造形式簡單、計算結果準確而被經常使用。最典型的3種Copula函數分別為Clayton函數、Frank函數和Gumbel函數。3種函數形式如(以二元為例)表1所示。

表1 典型阿基米德Copula函數

當考慮競爭的因素時,如果系統需要保持可靠,則包括突發失效和功能失效在內的任何失效模式都不允許出現。綜合公式(2)、公式(4)和公式(9),可得機械系統的可靠度表達式。

(11)

式中：FSi和F0i分別表示第i個功能失效模式指標的實際值和理想值；δi為第i個功能失效模式許用誤差,其中i=1,2,…m；N(t)表示速率為λd的泊松過程。基于公式(11),即可利用Copula函數計算得到多失效模式相關機械系統可靠度及其演化規律。

4 案例:某飛機艙門鎖機構

4.1 鎖機構組成及工作原理分析

本節針對某型飛機艙門鎖機構進行競爭失效分析。根據外場記錄數據了解到該型號鎖機構存在3種失效模式,失效模式一描述為因外界突發載荷造成的部件破壞,該失效模式屬突發型失效;失效模式二描述為因摩擦因數惡化和彈簧剛度演化引起開鎖過程阻力過大導致的功能失效,該失效模式屬于退化型功能失效;失效模式三描述為關鎖過程中總阻力過大而無法鎖閉導致的功能失效,該失效模式屬于退化型功能失效。

鎖機構的組成結構如圖2所示,主要包括作動筒、活塞、鎖鉤以及進行運動傳遞的3個連桿。作動

圖2 鎖機構組成

圖3 鎖機構工作狀態

筒為驅動源,通過活塞帶動各個連桿按照既定規律運動;鎖鉤為鎖機構運動傳遞的目標構件;鎖環固定在飛機機身,通過與鎖鉤的相對位置關系實現鎖機構的關鎖和開鎖功能。鎖機構的工作狀態如圖3所示。

開鎖過程詳細描述為,鎖機構處于鎖閉狀態,在液壓力的作用下鎖鉤克服鎖環負載和鎖機構的摩擦阻力,從鎖鉤鎖定位置運動至鎖鉤最大旋轉角度位置,這個階段需要最大的液壓驅動力克服運動過程中的阻力;之后從鎖鉤最大旋轉位置運動到鎖鉤鎖環脫離接觸的位置,此時作用在鎖鉤上的載荷成為運動的動力而不是阻力;最后鎖環鎖鉤脫離接觸運動到打開位置,此時鎖機構空載,液壓力只需克服較小的重力。關鎖過程詳細描述為,作動筒驅動活塞和一系列連桿帶動鎖鉤從初始位置開始運動到接觸鎖環,在鎖鉤鎖環接觸之前鎖鉤空載,整個機構的阻力主要來源于各個構件的重力;之后鎖鉤接觸鎖環繼續運動到鎖鉤旋轉到最大角度位置,此階段鎖鉤開始接受來自鎖環傳遞的負載,且負載的作用極大地增加了鎖機構運動的阻力,需要提供較大的液壓驅動力;最后從鎖鉤的最大旋轉角度回復到鎖死位置。這個過程鎖環傳遞給鎖鉤的拉力變成了驅動力,液壓力所需要克服的運動阻力值也較小。鎖閉狀態時,鎖機構始終承受鎖環位置傳遞的作用力,以及飛機運行過程中振動等原因引起的突發沖擊載荷。

4.2 鎖機構競爭失效建模

4.2.1 鎖機構突發失效建模

鎖機構處于鎖閉狀態時,鎖鉤與鎖環緊密連接,受到機身傳遞來的突發沖擊載荷作用,這類載荷的作用時間和量值都具有隨機性。以Wi表示載荷的量值大小,H表示鎖機構構件所能承受的臨界值,根據公式(2),鎖機構不發生突發失效的概率表示為

P(Wi

(12)

式中，Wi是獨立同分布(independent identically distributed)的變量,表示每次出現突發載荷的量值。使用泊松過程表示突發載荷的隨機出現時刻性,則機構不發生突發失效的概率表示為

(13)

式中，N(t)表示泊松過程,用于表示突發沖擊載荷出現時刻的不確定性。不發生突發失效要求沒有突發沖擊載荷出現(即式中i取值為0)或出現的載荷量值皆小于閾值H。

4.2.2 鎖機構退化型功能失效建模

根據對鎖機構工作狀態的分析,失效模式二和失效模式三發生于開/關鎖過程中的驅動力無法克服運動阻力之時。開閉鎖過程中的運動阻力都受到鎖鉤處摩擦因數、鎖鉤與鎖環接觸的角度、鎖鉤處的最大接觸力和彈簧的彈性系數影響,且摩擦因數和彈簧剛度具有明顯的隨時間演化的特征,故此二類失效模式屬于具有退化特征的功能失效,功能表征量分別為“開鎖過程中的最大阻力”和“關鎖過程中的最大阻力”,分別使用Fopen和Fclose表示。因此鎖機構不發生開鎖功能失效和關鎖功能失效的概率分別表示為公式(14)和公式(15)。

Ropen(t)=Fopen(t)

(14)

Rclose(t)=Fclose(t)

(15)

式中：F(t)表示開鎖或關鎖過程中的最大阻力會因退化的存在隨使用時間發生變化；FM表示作動筒所能提供的最大驅動力。為了提高計算效率,使用二階響應面構造功能表征量及其影響因素之間的數學關系,步驟如下:

Step1 依據運動原理搭建參數化機構動力學仿真模型;

Step2 按照各輸入參數的分布(見表1～3)進行拉丁超立方抽樣,將抽樣值代入仿真模型,得到成組的輸入輸出;

Step3 利用Step2中得到的結果構造響應面并進行精度驗證。

構造得到的鎖機構開關鎖過程功能表征量(即所需驅動力)的響應面函數分別為

Fopen(t)=-4 416.10+1 382.69x1+

23.93x1x2+0.32x1x3-0.02x1x4-0.002 4x2x3-

0.000 39x2x4-1.06×10-6x3x4

Fclose(t)=4 557 030.76-33 176.20x1+

115.15x1x2+0.33x1x3+5.07x1x4-

0.009 1x2x3+0.093x2x4+0.000 23x3x4

式中，xi(i=1,2,3,4)為相互獨立的影響因素,詳細含義及取值見表2。當所需驅動力大于作動筒所能提供的最大驅動力時,就發生相應階段的功能失效。根據公式(11)可得,該鎖機構的可靠度表達式進一步推導為公式(16)。

(16)

表2 參數含義及取值

4.3 計算結果

綜合公式(12)～(16)及鎖機構開關鎖過程功能表征量的響應面函數,將各參數取值(見表2)代入,即可得到計算結果。

圖4為開鎖過程所需驅動力的概率密度分布演化示意圖,圖5為閉鎖過程所需驅動力的概率密度分布演化示意圖。從圖中可以看出,隨著開閉鎖循環次數的增加,開鎖過程和閉鎖過程所需驅動力的均值逐步上升,且分散性逐漸增大。

圖4 開鎖過程所需驅動力分布演化圖 圖5 閉鎖過程所需驅動力分布演化圖圖6 可靠度變化曲線

根據公式(16)可以對鎖機構可靠度進行計算,鎖機構可靠度的演化結果如圖6所示。圖中橫坐標為使用鎖機構開閉鎖循環次數表示的等效時間,縱坐標為可靠度。分別將基于3種類型Copula函數、獨立假設情況與蒙特卡洛方法計算所得結果相比,各個方法的均方根誤差(root mean square error,RMSE)如表3所示。根據計算結果可以得到,鎖機構的可靠度隨著使用而逐漸發生退化,以蒙特卡洛方法計算結果為對照,相比于獨立假設,本文所提基于Copula函數的計算結果誤差更小,且其中Gumbel函數形式誤差最小,更適合于本案例的研究。

表3 各類方法所得結果誤差

5 結 論

本文針對處于競爭失效模式作用下的多功能機械系統，基于Copula函數建立了可靠性分析方法。Copula函數可以有效地建立多功能相關情況下功能失效的聯合概率分布函數，進而結合突發失效，計算得到機械系統的總可靠度，并利用所提方法，對某飛機艙門鎖機構進行了可靠性分析。根據分析和計算的結果，可以得到如下結論：

1) 飛機艙門鎖機構開閉鎖所需驅動力隨時間發生演化，且分散性逐漸變大；

2) 本文所提基于Copula函數的考慮失效模式相關性的可靠性分析結果更接近蒙特卡洛方法所得結果，證明本文所提方法的正確性；

3) 針對本文所提案例，Gumbel形式的Copula函數計算結果的RMSE為6.19×10-4，誤差最小，更適合于對此案例進行分析。