楊 濤，張 鈺 奇，李 凱 旋，韓 明 海，趙 華 東

(1.鄭州大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，河南 鄭州 450001； 2.河南省陸渾水庫(kù)管理局，河南 洛陽(yáng) 471000)

0 引 言

水工閘門隨著運(yùn)行年限的增加，會(huì)出現(xiàn)腐蝕、磨損等情況，從而嚴(yán)重影響閘門的運(yùn)行安全[1-2]。一旦閘門出現(xiàn)破壞，會(huì)帶來(lái)巨大的經(jīng)濟(jì)損失，甚至?xí)斐扇藛T傷亡[3]。目前，隨著數(shù)字化、信息化技術(shù)的發(fā)展，對(duì)閘門進(jìn)行健康狀態(tài)評(píng)估越來(lái)越重要。剩余壽命預(yù)測(cè)作為健康監(jiān)測(cè)的核心部分，對(duì)于指導(dǎo)閘門進(jìn)行預(yù)測(cè)性維護(hù)以保證閘門安全運(yùn)行具有重要意義。

近年來(lái)，閘門壽命預(yù)測(cè)方法越來(lái)越受到學(xué)者們的重視。趙林章等[4]利用分形維數(shù)法對(duì)閘門進(jìn)行了壽命預(yù)測(cè)；Abdul等[5]計(jì)算了閘門側(cè)輪和軌道的疲勞壽命；徐衛(wèi)敏等[6]利用Gamma隨機(jī)過(guò)程對(duì)鋼板銹蝕量進(jìn)行隨機(jī)模擬，并采用Monte-Carlo法對(duì)閘門上某構(gòu)件進(jìn)行了時(shí)變可靠度計(jì)算；李偉康等[7]考慮到閘門銹蝕對(duì)材料性能會(huì)產(chǎn)生影響，基于D-H曲線推導(dǎo)出了閘門構(gòu)件壽命預(yù)測(cè)的公式。但以上研究還存在著一些不足之處，主要體現(xiàn)在壽命預(yù)測(cè)的精度不高，用某個(gè)構(gòu)件代替整個(gè)閘門，缺乏合理性。

基于隨機(jī)過(guò)程的壽命預(yù)測(cè)方法可以很好地描述裝備的不確定性，現(xiàn)已逐漸被應(yīng)用于長(zhǎng)壽命、高可靠度的裝備壽命預(yù)測(cè)[8-9]。Huang等[10]基于逆高斯過(guò)程對(duì)刀具壽命進(jìn)行預(yù)測(cè)，達(dá)到了較高的精度；Lin等[11]基于Gamma過(guò)程對(duì)鋰電池剩余壽命進(jìn)行了預(yù)測(cè)；Wang等[12]利用維納過(guò)程建立了狀態(tài)空間模型，并對(duì)軸承進(jìn)行了壽命預(yù)測(cè)；Wu等[13]基于維納過(guò)程對(duì)復(fù)雜機(jī)電系統(tǒng)進(jìn)行了壽命預(yù)測(cè)。但是可以發(fā)現(xiàn)，以往研究未充分考慮閘門在運(yùn)行過(guò)程中的復(fù)雜工況，且單個(gè)指標(biāo)難以準(zhǔn)確反映出其健康狀態(tài)。Copula函數(shù)可以合理刻畫(huà)出不同性能指標(biāo)間的相關(guān)特性。Zhou等[14]運(yùn)用Copula函數(shù)表征了腐蝕管道不同缺陷間的相關(guān)關(guān)系；Song等[15]利用Copula函數(shù)分析了軸承多種失效模式的相關(guān)性，并進(jìn)行了壽命預(yù)測(cè)。因此本文采用多個(gè)指標(biāo)進(jìn)行壽命預(yù)測(cè)。

基于上述分析，本文選取閘門在腐蝕情況下變形和應(yīng)力兩個(gè)性能指標(biāo)來(lái)表征其健康狀態(tài)。首先，基于維納過(guò)程對(duì)閘門變形和應(yīng)力進(jìn)行退化性建模；其次利用Copula函數(shù)分析其相關(guān)特性，得到閘門剩余壽命概率密度函數(shù)；通過(guò)分步極大似然估計(jì)法對(duì)維納過(guò)程和Copula函數(shù)中的參數(shù)進(jìn)行估計(jì)；最后對(duì)比分析了基于二元維納過(guò)程和只考慮應(yīng)力影響下的一元維納過(guò)程進(jìn)行壽命預(yù)測(cè)的精度，同時(shí)利用剩余壽命預(yù)測(cè)結(jié)果指導(dǎo)閘門確定了合理的維修時(shí)間。

1 基于二元維納過(guò)程的退化建模

1.1 維納過(guò)程

當(dāng)采用維納過(guò)程進(jìn)行閘門退化建模時(shí)，第k個(gè)性能指標(biāo)的退化過(guò)程可表示為

X(k)(t)=X(k)(0)+λ(k)t+σ(k)B(t)

(1)

式中：X(k)(t)為t時(shí)刻第k個(gè)指標(biāo)退化值；X(k)(0)為退化初始值；λ(k)為漂移系數(shù)；σ(k)為擴(kuò)散系數(shù)；B(t)為標(biāo)準(zhǔn)維納過(guò)程。

1.2 基于Copula函數(shù)的相關(guān)性分析

由于閘門的變形和應(yīng)力兩個(gè)指標(biāo)之間不是相互獨(dú)立，也不是簡(jiǎn)單的線性關(guān)系，所以需要選擇出合理的方法刻畫(huà)不同指標(biāo)間的關(guān)系。Copula函數(shù)能有效地將聯(lián)合分布函數(shù)與邊緣分布函數(shù)聯(lián)合起來(lái)，常用于分析不同指標(biāo)之間的相關(guān)特性。因此，本文利用Copula函數(shù)來(lái)分析閘門上多元性能指標(biāo)之間的相關(guān)性，如式(2)所示：

F(x1，x2，…，xn)=C(F1(x1)，F(xiàn)2(x2)，…，F(xiàn)n(xn))

=C(u1，u2，…，un；θ)

(2)

式中：F(x1，x2，…，xn)為各個(gè)變量聯(lián)合分布函數(shù)；F1(x1)，F(xiàn)2(x2)，…，F(xiàn)n(xn)為各個(gè)變量的邊緣分布函數(shù)；C為Copula函數(shù)；θ為Copula中的參數(shù)。

Copula函數(shù)常見(jiàn)的4種類型如表1所列，其中u1，u2表示邊緣分布函數(shù)。

由Sklar定理得，閘門變形和應(yīng)力的相關(guān)特性可表示為[16]

F(x1，x2)=C[F1(x1)，F(xiàn)2(x2)；θ]

(3)

由于不同的Copula函數(shù)對(duì)于描述變形和應(yīng)力兩個(gè)退化指標(biāo)的相關(guān)特性具有較大影響，所以應(yīng)結(jié)合閘門實(shí)際退化情況選擇合適的Copula函數(shù)。AIC信息準(zhǔn)則鼓勵(lì)數(shù)據(jù)擬合的優(yōu)良性，被廣泛運(yùn)用于評(píng)價(jià)模型擬合數(shù)據(jù)效果的優(yōu)劣。AIC值越小，則說(shuō)明擬合的效果越好。本文采用AIC準(zhǔn)則選擇Copula函數(shù)：

AIC=-2lnA+2m

(4)

式中：A為模型對(duì)應(yīng)的似然函數(shù)；m為模型中參數(shù)的個(gè)數(shù)。

2 剩余壽命預(yù)測(cè)模型

2.1 剩余壽命分布

當(dāng)閘門兩個(gè)性能指標(biāo)均服從維納過(guò)程時(shí)，其剩余壽命服從逆高斯分布[17]，剩余壽命的概率密度函數(shù)為

(5)

式中：ωk為第k個(gè)性能指標(biāo)的失效閾值；xjk為第k個(gè)性能指標(biāo)在tj時(shí)刻的性能退化值。

當(dāng)獲得閘門上變形和應(yīng)力各自的剩余壽命概率密度函數(shù)后，由Copula函數(shù)可以得到應(yīng)力和變形剩余壽命分布函數(shù)的聯(lián)合概率密度函數(shù)：

(6)

式中：F1(t)、F2(t)為變形和應(yīng)力剩余壽命邊緣分布函數(shù)；c(F1(t)，F(xiàn)2(t)；θ)為C(F1(t)，F(xiàn)2(t)；θ)的概率密度函數(shù)。

定義變形和應(yīng)力中任一性能指標(biāo)超過(guò)失效閾值ωk，便認(rèn)為裝備失效。則閘門的剩余壽命T=min(T1，T2)。

因此，閘門剩余壽命的分布函數(shù)為

FRUL=P{min(T1，T2)≤t}

=F1(t)+F2(t)-C(F1(t)，F(xiàn)2(t))

(7)

閘門剩余壽命的概率密度函數(shù)為

(8)

2.2 參數(shù)估計(jì)

由于要估計(jì)的參數(shù)較多，分步極大似然估計(jì)法具有計(jì)算簡(jiǎn)單、流程清晰等優(yōu)點(diǎn)，所以采用分步極大似然估計(jì)法更新模型中的參數(shù)。

(1) 由維納過(guò)程的性質(zhì)可知，增量服從正態(tài)分布：

(9)

式中：ΔXk為第k個(gè)性能指標(biāo)退化量增量。

所以參數(shù)(λk，σk)的極大似然函數(shù)為

(10)

對(duì)λk，σk求偏導(dǎo)，可得λk，σk的極大似然估計(jì)值：

(11)

(12)

(13)

3 工程實(shí)例分析

陸渾水庫(kù)溢洪道露頂式弧形鋼閘門構(gòu)件采用Q235鋼，閘門面板中心弧面半徑為12 m，面板、縱梁厚度為12 mm，支臂厚度為16 mm，主橫梁厚度為18 mm，作用水頭為9.84 m。閘門變形和應(yīng)力作為閘門上重要的性能指標(biāo)，在運(yùn)行過(guò)程中會(huì)發(fā)生退化，可反映出閘門的健康狀態(tài)狀況。變形會(huì)影響閘門安全運(yùn)行[18]，本文定義變形閾值為10 mm。根據(jù)SL 74-2019《水利水電工程鋼閘門設(shè)計(jì)規(guī)范》，容許應(yīng)力與閘門運(yùn)行情況和鋼材的厚度有關(guān)[19]。該弧形閘門最大應(yīng)力處構(gòu)件厚度均不大于16 mm，屬于鋼材尺寸分組中的第一組，同時(shí)取修正系數(shù)K=0.9，所以定義閘門上最大應(yīng)力閾值為144 MPa。

本文采用文獻(xiàn)[20]中所論述的Guedes Soares銹蝕深度非線性變化公式，如式(14)所示，并根據(jù)文獻(xiàn)[21]中的統(tǒng)計(jì)數(shù)據(jù)確定參數(shù)，計(jì)算時(shí)段取13 a，同時(shí)利用有限元仿真得到閘門不同年份上的變形和應(yīng)力變化。

(14)

式中：d∞為腐蝕量的長(zhǎng)期厚度，全面腐蝕取d∞=1.6 mm，局部腐蝕取d∞=2.5 mm；d(t)為t時(shí)刻的腐蝕量厚度；τc為涂層壽命，取τc=5 a；τt為過(guò)渡時(shí)間，取τt=4.5 a。

由于ANSYS中的APDL模塊有著簡(jiǎn)單的格式，可以方便對(duì)有限元模型進(jìn)行計(jì)算，所以利用APDL對(duì)弧形鋼閘門進(jìn)行參數(shù)化建模和求解。由于閘門屬于空間薄壁結(jié)構(gòu)，因此采用shell 181單元進(jìn)行建模。閘門上銹蝕主要分為全面銹蝕和局部銹蝕：對(duì)于全面銹蝕，采用平均蝕余厚度法進(jìn)行模擬；對(duì)于局部銹蝕，通過(guò)改變銹蝕坑所在位置節(jié)點(diǎn)的厚度來(lái)實(shí)現(xiàn)對(duì)銹蝕坑的模擬[22]。根據(jù)閘門實(shí)際運(yùn)行情況，局部銹蝕一般情況下出現(xiàn)在閘門的面板右下角漏水的部位、右下角邊梁處以及下主梁經(jīng)常性積水區(qū)域，所以通過(guò)改變上述區(qū)域的節(jié)點(diǎn)厚度實(shí)現(xiàn)對(duì)局部銹蝕坑的模擬。閘門的局部銹蝕區(qū)域如圖1所示。

圖1 閘門局部銹蝕區(qū)域Fig.1 Local corrosion area of gate

利用APDL對(duì)閘門進(jìn)行有限元分析，可以得到閘門的變形和等效應(yīng)力。閘門初始狀態(tài)的變形和等效應(yīng)力云圖如圖2～3所示。

圖2 閘門變形云圖Fig.2 Gate deformation nephogram

圖3 閘門應(yīng)力云圖Fig.3 Gate stress nephogram

根據(jù)腐蝕非線性規(guī)律改變閘門對(duì)應(yīng)部分的厚度，同時(shí)利用APDL仿真可以得出閘門在退化過(guò)程中不同運(yùn)行年限的變形和等效應(yīng)力，將不同運(yùn)行年限的閘門變形和應(yīng)力值連接成線，可以得到閘門上變形和應(yīng)力退化曲線，如圖4所示。

圖4 閘門性能指標(biāo)退化曲線Fig.4 Degradation curve of gate performance indices

為了驗(yàn)證使用維納過(guò)程模型對(duì)閘門因腐蝕而產(chǎn)生的退化過(guò)程建模是否合理，需要對(duì)變形和應(yīng)力退化數(shù)據(jù)進(jìn)行檢驗(yàn)。由于K-S檢驗(yàn)法能夠利用樣本數(shù)據(jù)推斷樣本的總體是否服從某一理論分布，因此利用K-S檢驗(yàn)法對(duì)閘門退化數(shù)據(jù)進(jìn)行檢驗(yàn)，檢驗(yàn)結(jié)果如表2所列。其中：h=0表示在置信水平為α=0.05時(shí)接受退化量增量為正態(tài)分布的原假設(shè)，h=1時(shí)表示拒絕原假設(shè)；當(dāng)p>α?xí)r，表示不拒絕退化增量為正態(tài)分布的原假設(shè)。通過(guò)分析可以看出閘門上的變形和應(yīng)力退化量均滿足正態(tài)分布假設(shè)，故可以采用維納過(guò)程模型來(lái)描述閘門退化過(guò)程。

表2 K-S檢驗(yàn)結(jié)果

通過(guò)式(5) 得到變形和應(yīng)力的剩余壽命邊緣概率密度函數(shù)后，利用AIC準(zhǔn)則選擇合適的Copula函數(shù)，不同Copula函數(shù)得到的AIC值如表3所列。由表3可以看出：Gumbel函數(shù)的AIC值最小，所以選擇Gumbel函數(shù)來(lái)分析變形和應(yīng)力的相關(guān)特性。由Gumbel函數(shù)得到變形和應(yīng)力的剩余壽命聯(lián)合概率密度函數(shù)圖如圖5所示。由圖5可以看出Gumbel函數(shù)具有不對(duì)稱的尾部結(jié)構(gòu)，上尾高，下尾低，這說(shuō)明在分布的上尾部變形和應(yīng)力之間具有較強(qiáng)的相關(guān)性，在分布的下尾部?jī)烧咧g則為漸進(jìn)獨(dú)立的。

表3 4種Copula函數(shù)的AIC值

圖5 二元Gumbel-Copula概率密度函數(shù)Fig.5 Density function of 2D Gumbel-Copula function

在選擇好合適的Copula函數(shù)之后，由式(9)～(13)可得到不同運(yùn)行年限的參數(shù)值，如表4所列。之后結(jié)合式(8)可以得到不同時(shí)刻的閘門剩余壽命概率密度函數(shù)，同時(shí)計(jì)算出基于應(yīng)力退化的一元維納過(guò)程預(yù)測(cè)的剩余壽命概率密度函數(shù)，將閘門概率密度函數(shù)最大值對(duì)應(yīng)的時(shí)刻作為剩余壽命的預(yù)測(cè)值，如圖6所示。由圖6可以看出，隨著閘門運(yùn)行年限的增加，數(shù)據(jù)的增多，剩余壽命預(yù)測(cè)結(jié)果的分布范圍越來(lái)越小，說(shuō)明其預(yù)測(cè)結(jié)果的精確度越來(lái)越高，閘門發(fā)生失效的概率越來(lái)越高，閘門壽命減少，符合工程實(shí)際。

表4 不同運(yùn)行年限的參數(shù)估計(jì)值

為了能夠清晰地看出基于一元維納過(guò)程和二元維納過(guò)程預(yù)測(cè)閘門剩余壽命的變化趨勢(shì)，將閘門實(shí)際壽命取SL 226-98《水利水電金屬結(jié)構(gòu)報(bào)廢標(biāo)準(zhǔn)》中所規(guī)定的大型閘門折舊年限30 a[23]，將兩種方法預(yù)測(cè)的閘門剩余壽命與真實(shí)剩余壽命相比較，如圖7所示。從圖7可以看出隨著閘門運(yùn)行年限的增加，基于二元維納過(guò)程的閘門剩余壽命預(yù)測(cè)與真實(shí)壽命差距較小。

圖6 閘門剩余壽命概率密度Fig.6 Probability density of remaining life of gate

圖7 不同方法預(yù)測(cè)壽命變化趨勢(shì)Fig.7 Predicting trends of life expectancy by differet methods

為了客觀比較兩種方法預(yù)測(cè)閘門剩余壽命的精度，引入均方根誤差(RMSE)和相對(duì)誤差指標(biāo)進(jìn)行評(píng)價(jià)[24]。定義相對(duì)誤差指標(biāo)E：

(15)

式中：S為閘門剩余壽命預(yù)測(cè)值；T為閘門實(shí)際剩余壽命。

一元維納過(guò)程和二元維納過(guò)程預(yù)測(cè)結(jié)果的RMSE分別為4.72和3.46，這說(shuō)明相較于單一指標(biāo)，兩個(gè)指標(biāo)更能全面準(zhǔn)確反映出閘門的壽命狀態(tài)。不同預(yù)測(cè)方法的相對(duì)誤差E如表5所列，由表5可以看出二元維納過(guò)程的E<20%占比達(dá)到85.8%，而一元維納過(guò)程E<20%僅占比28.6%，同時(shí)，基于二元維納過(guò)程的E≥25%，占比為0，這說(shuō)明基于二元維納過(guò)程進(jìn)行閘門壽命預(yù)測(cè)的方法可信度較高。

表5 不同預(yù)測(cè)方法的相對(duì)誤差指標(biāo)

通過(guò)對(duì)閘門剩余壽命預(yù)測(cè)可以指導(dǎo)閘門確定合理的維修時(shí)間，既避免了經(jīng)常性維修造成的資金浪費(fèi)，又避免了閘門因過(guò)久不維修而發(fā)生破壞，造成巨大的安全事故?；陔S機(jī)過(guò)程預(yù)測(cè)閘門維修時(shí)間，可以根據(jù)閘門的重要性程度和維修的費(fèi)用，確定閘門上的維修閾值ωm[25]。

本文假設(shè)閘門上的變形維修閾值為9.5 mm，應(yīng)力維修閾值為143 MPa，當(dāng)閘門到達(dá)維修閾值后，便進(jìn)行一次全面大修，以延長(zhǎng)閘門的剩余使用壽命。由于數(shù)據(jù)越多，預(yù)測(cè)越精確，因此根據(jù)13 a的數(shù)據(jù)，將維修閾值代入式(5)～(13)可以得到閘門的維修時(shí)間，如圖8所示。由圖8可知，維修閾值下剩余壽命還剩6 a，說(shuō)明至少應(yīng)在閘門運(yùn)行19 a的時(shí)候進(jìn)行一次大修，以確保閘門的運(yùn)行安全。

圖8 閘門維修時(shí)間預(yù)測(cè)Fig.8 Forecast on gate maintenance time

4 結(jié) 論

(1) 利用二元維納過(guò)程描述腐蝕情況下閘門退化過(guò)程，預(yù)測(cè)閘門剩余壽命的均方根誤差RMSE為3.46，低于單一應(yīng)力指標(biāo)的RMSE值，說(shuō)明用二元維納過(guò)程方法來(lái)描述閘門退化過(guò)程具有合理性。

(2) 考慮到單一指標(biāo)難以準(zhǔn)確反映閘門的健康狀態(tài)，利用變形和應(yīng)力兩個(gè)退化性能指標(biāo)來(lái)預(yù)測(cè)閘門壽命，運(yùn)用Gumbel函數(shù)描述其相關(guān)特性，預(yù)測(cè)壽命相對(duì)誤差指標(biāo)E≥25%占比為0，說(shuō)明基于二元維納過(guò)程的閘門剩余壽命預(yù)測(cè)具有較高的精度。

(3) 基于閘門剩余壽命預(yù)測(cè)結(jié)果，確定維修閾值，確定了閘門至少在運(yùn)行19 a進(jìn)行一次大修，以延長(zhǎng)閘門使用壽命。同時(shí)，閘門上的失效形式眾多，如何考慮更多的性能指標(biāo)和維修過(guò)程對(duì)閘門的退化過(guò)程的影響，以達(dá)到更高的預(yù)測(cè)精度有待進(jìn)一步研究。