基于蒙特卡洛重要抽樣法的結構時變可靠性分析

嚴宇飛，李方義

（1.長沙理工大學 汽車與機械工程學院，長沙 410114；2.長沙理工大學 橋梁工程安全控制省部共建教育部重點實驗室，長沙 410114；3.廣州大學 機械與電氣工程學院，廣州 510006）

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

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

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

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

1 結構時變可靠性模型

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

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

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

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

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

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

由可靠度定義知結構在設計使用期T內結構可靠度Ps（t）和失效概率Pf（t）為：

2 重要抽樣法理論

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

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

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

隨機變量X的第i個樣本向量表示為：xi＝（xi1，xi2，…，xin），則結構失效概率Pf的期望估計值為：

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

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

3 基于蒙特卡洛重要抽樣法的結構時變可靠性分析方法

蒙特卡洛重要抽樣法中樣本點的選擇通常是在隨機變量的均值最大可能點或由FOSM得到的MPP點附近選?。?3］。

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

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

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

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

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

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

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

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

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

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

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

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

圖2 算法流程框圖

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

2）在ti時間段內，利用FOSM的方法對處理后得到的結構時不變可靠性模型進行計算，得到隨機變量相對應的MPP點xi；

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

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

4 算例分析

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

4.1 某機械零件

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

表1 某機械零件的隨機變量分布情況

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

表2 設計基準期5 000 h時某機械零件的求解結果

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

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

表3 某機械零件各設計基準期的可靠度指標

表4 某機械零件各設計基準期的結果誤差

4.2 某結構的鋼筋混凝土短柱

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

表5 某結構的鋼筋混凝土短柱的隨機變量分布情況

鋼筋混凝土柱的抗力R（t）與混凝土抗壓強度fc和鋼筋屈服強度fy的關系如下：

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

當設計基準期為25年時，同4.1處理。對比MCS、FOSM以及本文方法3種求解結果，本文方法結果值相對誤差在0.19%～1.30%，具體見表6。

表6 設計基準期25年時某結構的鋼筋混凝土短柱的求解結果

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

表7 某結構的鋼筋混凝土短柱各設計基準期的可靠度指標

表8 某結構的鋼筋混凝土短柱各設計基準期的結果誤差

4.3 屋頂?shù)蔫旒芙Y構

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

圖3 屋頂桁架模型

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

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

表9 屋頂?shù)蔫旒芙Y構的隨機變量分布情況

當設計基準期為5年時，同上處理，對比MCS、FOSM以及本文方法3種方法求解結果，其中本文方法結果值相對誤差在0.01%～0.37%之間，具體見表10：

表10 設計基準期5年時屋頂?shù)蔫旒芙Y構的求解結果

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

表11 屋頂?shù)蔫旒芙Y構各設計基準期的可靠度指標

表12 屋頂?shù)蔫旒芙Y構各設計基準期的結果誤差

5 結論

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

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

3）本文方法為結構時變可靠性分析近似方法的驗證提供了新的途徑，也為工程實際問題中的結構時變可靠性分析方法提供了新的方法和思路。