基于蒙特卡洛重要抽樣法的結(jié)構(gòu)時變可靠性分析

嚴(yán)宇飛，李方義

（1.長沙理工大學(xué) 汽車與機(jī)械工程學(xué)院，長沙 410114；2.長沙理工大學(xué) 橋梁工程安全控制省部共建教育部重點(diǎn)實(shí)驗(yàn)室，長沙 410114；3.廣州大學(xué) 機(jī)械與電氣工程學(xué)院，廣州 510006）

在實(shí)際工程問題中，結(jié)構(gòu)在工況下會受到材料性能退化、隨機(jī)載荷過程等不確定因素的影響，且這些因素往往具有時變效應(yīng)，導(dǎo)致其結(jié)構(gòu)可靠性動態(tài)變化。因此，針對結(jié)構(gòu)動態(tài)可靠性的研究是很有必要的。結(jié)構(gòu)時變可靠性分析就是一種動態(tài)結(jié)構(gòu)可靠性分析，它符合實(shí)際工程所需，能夠客觀、準(zhǔn)確地對某工件或結(jié)構(gòu)的可靠性進(jìn)行預(yù)測或評估。當(dāng)前結(jié)構(gòu)時變可靠分析已有不少研究成果，如：基于跨越率方法［1］演變的“PHI2方法［2－5］”、基于“時間離散”的方法［6］總結(jié)的“準(zhǔn)靜態(tài)法［7－9］”等。然而與MCS相比，當(dāng)用上述方法求解非線性較強(qiáng)的結(jié)構(gòu)可靠性問題時，都有著求解復(fù)雜或求解精度不足的缺陷。

MCS求解結(jié)構(gòu)可靠性問題時，其求解精度僅與模擬次數(shù)有關(guān)，模擬次數(shù)越多，其精度越高，且與問題的維數(shù)無關(guān)。因此，MCS理論上可以求解所有的可靠性問題。MCS在結(jié)構(gòu)時變可靠性問題的應(yīng)用上有：王帥等［10］應(yīng)用MCS求解時域－誤差模型的運(yùn)動時變可靠性；Li［11］基于MCS模擬得到電車軌道連接螺栓工作載荷的概率分布，以分析連接螺栓的結(jié)構(gòu)時變可靠性；Li等［12］利用MCS驗(yàn)證其求解混合變量的時變模型結(jié)果的準(zhǔn)確性；Wang等［13］利用MCS驗(yàn)證其提出的方法求解結(jié)構(gòu)時變可靠性問題的準(zhǔn)確性等。但是由于MCS隨著模擬次數(shù)的增加，會導(dǎo)致計(jì)算量越來越巨大。因此，由于MCS的計(jì)算效率一般來說是比較低的，從而MCS通常被用來驗(yàn)證結(jié)果。

由于MCS具有高的求解精度，因此要將MCS應(yīng)用于結(jié)構(gòu)時變可靠性問題上就一定需要考慮改善并提高M(jìn)CS的計(jì)算效率的問題。當(dāng)前研究者針對這一問題引入了重要抽樣法［14－15］來提高M(jìn)CS的計(jì)算效率?！爸匾闃臃ā敝傅氖窃谀M次數(shù)一定的情況下，通過改變抽樣中心的位置或用新的概率分布對隨機(jī)變量進(jìn)行抽樣，來估計(jì)失效概率的方法。又由于在結(jié)構(gòu)時變可靠性問題上，重要抽樣法的應(yīng)用存在著反應(yīng)量的概率密度函數(shù)未知的問題和抽樣函數(shù)需求的時間多維聯(lián)合概率密度函數(shù)難以獲得的問題［16］，因此當(dāng)前重要抽樣法主要是應(yīng)用在靜態(tài)的結(jié)構(gòu)可靠性問題上，而應(yīng)用在結(jié)構(gòu)時變可靠性問題上的研究成果寥寥無幾，Prablwarter等［17－18］基于俄羅斯輪盤賭、分裂方法對蒙特卡洛重要抽樣法進(jìn)行研究，但該方法對計(jì)算機(jī)運(yùn)算能力的要求比較高；宋曉通和謝紹宇等［19－20］在電力系統(tǒng)可靠性評估上提出了自適應(yīng)重要抽樣算法，但該方法應(yīng)用于其他領(lǐng)域上時，還需要進(jìn)行一定的改善；羅立勝［21］以串聯(lián)理論處理時間變量，提出了一種銹蝕鋼構(gòu)件時變可靠度的重要抽樣法，但對該方法的普遍適用性未進(jìn)行探究。因此，針對結(jié)構(gòu)時變可靠性問題，研究一種精度好、理論簡單、計(jì)算簡便的蒙特卡洛重要抽樣方法是很有必要的。

本文將“蒙特卡洛重要抽樣法”應(yīng)用于結(jié)構(gòu)時變可靠性問題之中，并結(jié)合“準(zhǔn)靜態(tài)法”和“FOSM”，提出一種基于蒙特卡洛重要抽樣法的結(jié)構(gòu)時變可靠性分析方法，使得本文方法不僅能夠達(dá)到接近MCS的求解精度的程度，而且與MCS相比，求解效率也有很大提高。本文給出本文方法的計(jì)算原理和具體的計(jì)算步驟，并以MCS、FOSM以及本文方法對3個算例分別進(jìn)行計(jì)算，通過3種方法的求解結(jié)果對比，證明本文方法的優(yōu)越性，為蒙特卡洛重要抽樣法應(yīng)用于結(jié)構(gòu)動態(tài)可靠性分析上提供新的方法和思路。

1 結(jié)構(gòu)時變可靠性模型

通常來說，大多結(jié)構(gòu)時變可靠性模型的研究都是基于廣義的抗力模型，其極限狀態(tài)方程Z（t）一般由結(jié)構(gòu)抗力隨機(jī)過程R（t）和載荷效應(yīng)隨機(jī)過程S（t）組成的線性或非線性函數(shù)方程，其表達(dá)式為：

式中：t為時間變量；載荷效應(yīng)隨機(jī)過程S（t）一般含有永久載荷A（其通常為n維的隨機(jī)向量）和可變載荷隨機(jī)過程Q（t），即式（1）改寫為：

在實(shí)際工程問題中，結(jié)構(gòu)抗力隨機(jī)過程R（t）的常見的退化形式有指數(shù)退化形式、對數(shù)退化形式等，具體如下：

指數(shù)退化形式表達(dá)式：

對數(shù)退化形式表達(dá)式：

其中：R0為初始抗力；ε為衰減系數(shù)。

由可靠度定義知結(jié)構(gòu)在設(shè)計(jì)使用期T內(nèi)結(jié)構(gòu)可靠度Ps（t）和失效概率Pf（t）為：

2 重要抽樣法理論

MCS生成的隨機(jī)變量樣本點(diǎn)多集中在聯(lián)合概率密度函數(shù)的最大值點(diǎn)附近，但實(shí)際上結(jié)構(gòu)失效是個小概率事件，從而在失效域內(nèi)獲取到樣本機(jī)會小，導(dǎo)致MCS的求解精度較低，因此MCS要大量的樣本來提高算法精度。蒙特卡洛重要抽樣法是通過改變隨機(jī)抽樣的中心，使得樣本點(diǎn)更有可能落在失效域中，使得模擬次數(shù)適量減少也可保證算法求解精度和效率。

設(shè)X＝（X1，X2，…，Xn）是結(jié)構(gòu)的隨機(jī)變量，其概率分布為fX（x），Z＝g（x）為功能函數(shù)，結(jié)構(gòu)失效概率Pf為：

式中I［g（x）］為g（x）的指示函數(shù)，表示為：

隨機(jī)變量X的第i個樣本向量表示為：xi＝（xi1，xi2，…，xin），則結(jié)構(gòu)失效概率Pf的期望估計(jì)值為：

根據(jù)文獻(xiàn)［22］可知，通過改變抽樣方式，抽取的樣本點(diǎn)更可能落入失效域內(nèi)或所謂的重要區(qū)域內(nèi)，使得求解精度提升。因此，式（7）可改寫為：

式中px（x）為：“重要抽樣函數(shù)”，即取代fx（x）對X抽樣的函數(shù)。此時，Pf的期望估計(jì)值為：

3 基于蒙特卡洛重要抽樣法的結(jié)構(gòu)時變可靠性分析方法

蒙特卡洛重要抽樣法中樣本點(diǎn)的選擇通常是在隨機(jī)變量的均值最大可能點(diǎn)或由FOSM得到的MPP點(diǎn)附近選取［23］。

將其應(yīng)用于結(jié)構(gòu)時變可靠性問題時，首先必須先對時間變量進(jìn)行處理，使得動態(tài)可靠性問題變?yōu)殪o態(tài)可靠性問題。因此，利用文獻(xiàn)［6－9］中的方法將設(shè)計(jì)基準(zhǔn)期T劃分為m個均勻相等的時段，即每一個時間段為τ＝T／m，如圖1所示，此時可將抗力隨機(jī)過程R（t）和動載荷隨機(jī)過程Q（t）離散化處理，得到m個隨機(jī)變量Ri和Qi，以第i個時段抗力的中值為Ri，將式（5）轉(zhuǎn)化為：

圖1 抗力R（t）與可變載荷Q（t）隨時間的變化曲線

Qi的值是通過所在時間段內(nèi)的統(tǒng)計(jì)得到，通常情況下，可將Qi視為獨(dú)立同分布的。因此，由串聯(lián)體系可靠度理論可將式（6）改寫為：

從工程意義方面考慮，不妨將Q′取為m個Qi的最大值QT。同時，Qi是獨(dú)立同分布的，根據(jù)統(tǒng)計(jì)學(xué)的原理，設(shè)計(jì)基準(zhǔn)期T內(nèi)最大載荷效應(yīng)QT的概率分布函數(shù)FQT（x）為：

在實(shí)際問題中，通常情況下將載荷效應(yīng)Qi近似為以參數(shù)為ατ和μτ的極值I型分布，因此QT也將服從參數(shù)為αT和μT的極值I型分布，即：

引入新的隨機(jī)變量Q′，其概率密度函數(shù)為fQ′（q′），概率分布函數(shù)為FQ′（q′），式（13）可轉(zhuǎn)化為：

式中：fR1，R2，…，Rm（r1，r2，…，rm）是R1，R2，…，Rm的聯(lián)合概率密度函數(shù)；fA（a）為A的概率密度函數(shù)；FQr為Qi的概率分布函數(shù)；F－1Q′是FQ′的反函數(shù)。

此時，功能函數(shù)轉(zhuǎn)化為：

引入的隨機(jī)變量Q′沒有規(guī)定它的分布的類型，事實(shí)上，在以式（17）為功能函數(shù)利用一次2階矩法（FOSM）的迭代公式計(jì)算可靠指標(biāo)時，公式中只含有FQ′（Q′）項(xiàng)，這時可將其作為一個變量使用，也就是說整個分析過程與Q′的具體分布無關(guān)。但為與現(xiàn)行的結(jié)構(gòu)可靠度設(shè)計(jì)統(tǒng)一標(biāo)準(zhǔn)相協(xié)調(diào)，使Q′具有工程意義，將Q′取為m個Qi的最大值QT，因此可將式（17）轉(zhuǎn)化表示為：

由于最大可變載荷效應(yīng)QT服從極值Ⅰ型分布，因此由式（18）可推導(dǎo)得最終的極限狀態(tài)方程：

FOSM對式（19）進(jìn)行計(jì)算求解，得到各隨機(jī)變量對應(yīng)的MPP點(diǎn)。以對應(yīng)的MPP點(diǎn)進(jìn)行蒙特卡洛重要抽樣，得到各隨機(jī)變量的樣本點(diǎn)，并將其代入式（2）中進(jìn)行計(jì)算，統(tǒng)計(jì)功能函數(shù)失效的次數(shù)nf，然后以Pf＝nf／n得到結(jié)構(gòu)的失效概率，最后求得結(jié)構(gòu)的可靠度指標(biāo)。

本文方法具體步驟及流程（見圖2）如下：

圖2 算法流程框圖

1）以“準(zhǔn)靜態(tài)法”對結(jié)構(gòu)時變可靠性模型進(jìn)行處理，離散時間T為m段，使得隨機(jī)過程等效為隨機(jī)變量，并求出每段中的廣義抗力Ri和載荷QT，使得動態(tài)可靠性問題變?yōu)殪o態(tài)可靠性問題；

2）在ti時間段內(nèi)，利用FOSM的方法對處理后得到的結(jié)構(gòu)時不變可靠性模型進(jìn)行計(jì)算，得到隨機(jī)變量相對應(yīng)的MPP點(diǎn)xi；

3）在xi附近對隨機(jī)變量進(jìn)行重要抽樣法采樣，樣本容量為n；

4）將得到的樣本數(shù)代入式（2）之中計(jì)算，統(tǒng)計(jì)功能函數(shù)失效的次數(shù)nf，以Pf＝nf／n得到ti時間段內(nèi)結(jié)構(gòu)的失效概率，又通過失效概率可求出ti時間段內(nèi)的結(jié)構(gòu)的可靠度指標(biāo)。

4 算例分析

為驗(yàn)證本文方法的精度和效率，將其應(yīng)用于下列3個算例之中，通常情況下的，MCS的樣本容量為：n＝106～108。

4.1 某機(jī)械零件

已知某機(jī)械零件［24］的抗力隨時間的變化規(guī)律為R（t）＝R0exp（－at）MPa，a為衰減系數(shù)，其值為2×10－5，初始抗力R0服從（μR0σR0）MPa的對數(shù)正態(tài)分布，永久載荷效應(yīng)A服從（μA，σA）MPa的正態(tài)分布，如表1所示，以τ＝1 000 h作為一個時間段，動態(tài)載荷Q（t）在第i個時間段內(nèi)的統(tǒng)計(jì)值Qi服從參數(shù)為ατ和μτ的極值I型分布，其中參數(shù)ατ＝0.044 93 MPa－1，μτ＝49.77 MPa。計(jì)算該零件在不同設(shè)計(jì)基準(zhǔn)期時的可靠度指標(biāo)。

表1 某機(jī)械零件的隨機(jī)變量分布情況

當(dāng)設(shè)計(jì)基準(zhǔn)期為5 000 h時，對比MCS、FOSM以及本文方法3種方法求解結(jié)果，其中本文方法模擬次數(shù)分別從104、5×104、105、106進(jìn)行試驗(yàn)計(jì)算，同時為減小偶然誤差，均反復(fù)進(jìn)行了10次運(yùn)算。其中，當(dāng)模擬次數(shù)為104時，最大誤差達(dá)到了0.63%；5×104時，最大誤差達(dá)到了0.51%；105時，最大誤差達(dá)到了0.50%；106時，最大誤差達(dá)到了0.38%，且隨著模擬次數(shù)的增多，結(jié)果值趨近且穩(wěn)定。因此，本文選取模擬次數(shù)為5×104進(jìn)行分析，其結(jié)果值相對誤差在0.29%～0.63%，表2中本文方法可靠度指標(biāo)和相對誤差取10次運(yùn)算后的平均數(shù)，具體如下：

表2 設(shè)計(jì)基準(zhǔn)期5 000 h時某機(jī)械零件的求解結(jié)果

誤差1是指FOSM與MCS之間的誤差；誤差2是指本文方法與MCS之間的誤差。

結(jié)果分析：考慮到參數(shù)的時變特性，求解某機(jī)械零件在10個設(shè)計(jì)基準(zhǔn)期內(nèi)的可靠度指標(biāo)。從表2結(jié)果所示可知：本文方法僅用5×104的樣本數(shù)就能與MCS的結(jié)果接近，模擬次數(shù)僅為MCS的1／20，因此本文方法的求解效率與MCS相比，求解效率有明顯提高。從表3結(jié)果所示可知：隨著設(shè)計(jì)基準(zhǔn)期的延長，結(jié)構(gòu)或產(chǎn)品的可靠度指標(biāo)逐步下降，這完全符合實(shí)際情況。同時，從表4結(jié)果所示可知：在設(shè)計(jì)基準(zhǔn)期內(nèi)，F(xiàn)OSM的最大誤差為0.92%，而蒙特卡洛重要抽樣法的最大誤差為0.91%，均在誤差允許范圍內(nèi)。同時，從誤差1與誤差2的對比中可以發(fā)現(xiàn)：本文方法的精度更高、更準(zhǔn)確，說明本文方法求解結(jié)構(gòu)時變可靠性問題是有效且精度較高的。

表3 某機(jī)械零件各設(shè)計(jì)基準(zhǔn)期的可靠度指標(biāo)

表4 某機(jī)械零件各設(shè)計(jì)基準(zhǔn)期的結(jié)果誤差

4.2 某結(jié)構(gòu)的鋼筋混凝土短柱

某地方某結(jié)構(gòu)［6］的鋼筋混凝土短柱的截面尺寸寬和高分別為b和h，其參數(shù)值分別為300、375 mm?；炷恋膹?qiáng)度等級為C30，柱內(nèi)鋼筋面積As＝1 811.28 mm2，使用的鋼筋是Ⅱ級鋼筋。鋼筋混凝土短柱受抗力隨機(jī)過程為R（t），永久荷載效應(yīng)G以及可變荷載效應(yīng)Q影響，其中抗力隨機(jī)過程為R（t）與混凝土抗壓強(qiáng)度fc和鋼筋屈服強(qiáng)度fy有關(guān)，最大可變荷載效應(yīng)QT服從極值Ⅰ型分布，具體參數(shù)和分布見表5，假設(shè)設(shè)計(jì)基準(zhǔn)期為T＝50年。其他未提及系數(shù)默認(rèn)為確定性參數(shù)，計(jì)算該柱的可靠指標(biāo)和失效概率。功能函數(shù)為：

表5 某結(jié)構(gòu)的鋼筋混凝土短柱的隨機(jī)變量分布情況

鋼筋混凝土柱的抗力R（t）與混凝土抗壓強(qiáng)度fc和鋼筋屈服強(qiáng)度fy的關(guān)系如下：

查得該地方的抗力變化系數(shù)φc（t）和φy（t）為：

當(dāng)設(shè)計(jì)基準(zhǔn)期為25年時，同4.1處理。對比MCS、FOSM以及本文方法3種求解結(jié)果，本文方法結(jié)果值相對誤差在0.19%～1.30%，具體見表6。

表6 設(shè)計(jì)基準(zhǔn)期25年時某結(jié)構(gòu)的鋼筋混凝土短柱的求解結(jié)果

結(jié)果分析：本算例一共預(yù)設(shè)了10個設(shè)計(jì)基準(zhǔn)期，從表7所示可知，隨著設(shè)計(jì)基準(zhǔn)期的延長，某結(jié)構(gòu)的鋼筋混凝土短柱的可靠度指標(biāo)逐步下降。同時，在設(shè)計(jì)基準(zhǔn)期內(nèi)，F(xiàn)OSM的最大誤差為2.03%，而蒙特卡洛重要抽樣法的最大誤差為1.07%，這都是在允許誤差范圍內(nèi)。從表8可知：FOSM結(jié)果誤差偏離的更多，對比算例4.1，算例4.2中的FOSM求解精度降低許多，說明FOSM雖然可以滿足求解的需要，但由于其僅適應(yīng)獨(dú)立正態(tài)分布的變量，使得運(yùn)算中需要有正態(tài)空間的轉(zhuǎn)化以及偏導(dǎo)，從而產(chǎn)生誤差，又由于功能函數(shù)的近似處理和原時變可靠性問題的近似處理，使得FOSM求解精度會有所降低。而MCS由于大的模擬次數(shù)，能夠保證該方法求解精度，因此MCS的精度會高于FOSM。又從誤差1與誤差2的對比中發(fā)現(xiàn)：FOSM精度存在著不穩(wěn)定且精度一般的問題，而本文方法的精度依然較高，且計(jì)算精度比較穩(wěn)定，表6可進(jìn)一步證明本文方法相對MCS，效率明顯有所提高。因此，說明本文方法對于結(jié)構(gòu)時變可靠性分析來說是一個有效的且能夠保證精度的方法。

表7 某結(jié)構(gòu)的鋼筋混凝土短柱各設(shè)計(jì)基準(zhǔn)期的可靠度指標(biāo)

表8 某結(jié)構(gòu)的鋼筋混凝土短柱各設(shè)計(jì)基準(zhǔn)期的結(jié)果誤差

4.3 屋頂?shù)蔫旒芙Y(jié)構(gòu)

圖3所示是根據(jù)文獻(xiàn)［25］改編的一個屋頂?shù)蔫旒芙Y(jié)構(gòu)，底部和承拉的桿件材料為鋼，頂部和承壓的桿件用水泥加固。屋頂承受均勻的分布載荷q（t），q（t）可以轉(zhuǎn)化成等效的節(jié)點(diǎn)力P＝q（t）l／4，節(jié)點(diǎn)C處的垂向位移可以推導(dǎo)為：

圖3 屋頂桁架模型

其中：AC和AS分別表示水泥桿和鋼制桿件的截面積，其數(shù)值分別為0.034 m2和0.000 94 m2；EC和ES分別表示它們的彈性模量。

受載過程中，節(jié)點(diǎn)C位置的垂向位移應(yīng)小于D（t），且隨時間衰減，其衰減規(guī)律為D（t）＝D0［1＋ln（1－0.000 2t）］，其中D0為初始位移。表9列出了此問題中隨機(jī)變量的具體情況。功能函數(shù)定義為：

表9 屋頂?shù)蔫旒芙Y(jié)構(gòu)的隨機(jī)變量分布情況

當(dāng)設(shè)計(jì)基準(zhǔn)期為5年時，同上處理，對比MCS、FOSM以及本文方法3種方法求解結(jié)果，其中本文方法結(jié)果值相對誤差在0.01%～0.37%之間，具體見表10：

表10 設(shè)計(jì)基準(zhǔn)期5年時屋頂?shù)蔫旒芙Y(jié)構(gòu)的求解結(jié)果

結(jié)果分析：本算例一共預(yù)設(shè)了10個設(shè)計(jì)基準(zhǔn)期，由表11可知，隨著設(shè)計(jì)基準(zhǔn)期的延長，屋頂桁架結(jié)構(gòu)的可靠度指標(biāo)逐步下降。同時，在設(shè)計(jì)基準(zhǔn)期內(nèi)，F(xiàn)OSM的最大誤差為2.31%，而蒙特卡洛重要抽樣法的最大誤差為1.11%。很明顯，本文方法的精度遠(yuǎn)高于FOSM，且本算例中FOSM出現(xiàn)了和算例4.2一樣的情況。從表12分析可知：FOSM明顯存在著結(jié)果不穩(wěn)定且精度一般的問題。而本文方法的精度仍然較高，且計(jì)算結(jié)果比較穩(wěn)定。同時，根據(jù)表10中的結(jié)果可證明本文方法相比MCS，效率明顯有所提高。因此，可以說明本文方法對于結(jié)構(gòu)時變可靠性分析是一個有效的方法。

表11 屋頂?shù)蔫旒芙Y(jié)構(gòu)各設(shè)計(jì)基準(zhǔn)期的可靠度指標(biāo)

表12 屋頂?shù)蔫旒芙Y(jié)構(gòu)各設(shè)計(jì)基準(zhǔn)期的結(jié)果誤差

5 結(jié)論

1）本文方法與FOSM對比，雖然效率上有所不及，但在計(jì)算精度上有著明顯的優(yōu)勢且計(jì)算結(jié)果穩(wěn)定；

2）本文方法與MSC對比，兩者求解精度極其接近，且本文方法的模擬次數(shù)僅為MSC的1／20。而對于蒙特卡洛方法來說，模擬次數(shù)的降低意味著計(jì)算效率的提高，因此，本文方法對求解結(jié)構(gòu)時變可靠性問題有較好的精度和效率；

3）本文方法為結(jié)構(gòu)時變可靠性分析近似方法的驗(yàn)證提供了新的途徑，也為工程實(shí)際問題中的結(jié)構(gòu)時變可靠性分析方法提供了新的方法和思路。