結構可靠度分析的有限單元數(shù)值逼近法

劉旭

河南理工大學土木工程學院 焦作454000

引言

結構可靠性分析［1，2］對工程結構的設計、穩(wěn)定分析以及優(yōu)化設計具有著重要意義。結構可靠性的度量稱為可靠度。當今主流的可靠度理論認為結構體包含的各種參數(shù)如幾何尺寸、抗力、荷載效應等，多為服從某種概率分布的隨機變量，是其諸多不確定性因素的來源［3］。結構可靠度分析是將上述不確定性因素引入概率統(tǒng)計分析的理論［4］，描述結構體功能狀態(tài)的函數(shù)表示為：

式中：X ＝（X1，X2，…，Xn）是n個影響結構體功能的基本隨機變量。當Z＞0、Z＜0、Z＝0 時分別表示結構體處于可靠狀態(tài)、失效狀態(tài)和極限狀態(tài)。設X的聯(lián)合概率密度函數(shù)為fX（x），結構失效概率Pf的表達式為：

這是一個基于失效域的多維積分問題，通常由于極限狀態(tài)面（Z＝0）在X 空間中難以確定，因此直接采用式（2）計算Pf是非常困難的。目前行而有效的方法大體分為兩大類：近似法和模擬法。

近似法是使用平面或曲面替代原來的極限狀態(tài)曲面（失效面），此類方法主要有：國際結構安全性聯(lián)合委員會（JCSS）推薦的JC 法［5］，并以次為基礎發(fā)展了二次二階矩法、高次高階矩法；以及隨著計算機技術的發(fā)展，近些年主流理論中的響應面法［6，7］、支持向量機法［8］、神經(jīng)網(wǎng)絡擬合法［9］等。也可利用函數(shù)近似理論對結構的功能函數(shù)和分布函數(shù)做近似處理的矩法，例如：孟廣偉等人［10－12］提出利用降維算法，建立n個一維函數(shù)替代（1）式中的n維功能函數(shù)，并結合Edgeworth 級數(shù)擬合近似功能函數(shù)的分布函數(shù)從而直接計算失效概率。近似法的不足在于功能函數(shù)或失效面的數(shù)學特性對近似效果的影響較大，如JC法的計算精度很大程度上取決于失效面在設計驗算點附近的非線性程度，且該類方法難以得到準確的誤差范圍。模擬法以蒙特卡洛法（MCS）［13］為主，原理為對已知概率分布的隨機變量產(chǎn)生足夠多的隨機數(shù)，對功能函數(shù)進行重復模擬。雖然此類方法計算精度較高且無求解條件限制，但當處理小概率問題時，需要大量的模擬次數(shù)保證結果的精度［14］，顯得有些捉襟見肘。由于是基于隨機模擬的計算結果，所以該類方法同樣難以對結果做出準確的誤差分析（非概率誤差）。

因此，本文提出一種基于有限單元法［15］處理X分布空間，計算結構失效概率的新方法。該算法分為四個部分：（1）利用功能函數(shù)中基本變量的概率集中性，在預設誤差限的要求下，通過舍棄X的低概率區(qū)間保留高概率區(qū)間，建立有限單元法的計算域。（2）基于有限單元法對計算域進行均勻劃分，得到若干個互不重疊、相互獨立的單元，并規(guī)定以單元的幾何中心點所處位置表征其是否位于失效域。（3）引入功能函數(shù)對每個單元的中心點進行篩選，保留位于失效域的單元并累加其概率得到近似的失效概率。（4）通過MATLAB編程不斷提高劃分單元的精細度，逐次迭代、縮小誤差，實現(xiàn)數(shù)值逼近的計算過程，得到滿足預設誤差限的失效概率。

1 有限單元法處理過程

1.1 建立X空間的計算域

記Xi（Xi∈X，i＝1，2，…，n）的概率密度函數(shù)和分布函數(shù)分別為：fi（xi），F(xiàn)i（xi）。由于Xi的分布具有概率集中性，為提高計算效率，可利用此性質定義Xi的高概率區(qū)間，即包含大部分概率值的區(qū)間。記Xi的高概率區(qū)間為ωi，ωi

內(nèi)包含的概率為Pωi，如：

將式（5）轉變?yōu)閔與xi的關系，并對h取最小值得到：

因此在計算中根據(jù)要求給定Pωi后，即可利用式（6 ～8）求得對應式（3）中ωi的區(qū)間分布范圍。例如：已知X1服從正態(tài)分布（μX1＝20，σX1＝4）；X2服從對數(shù)正態(tài)分布（μX2＝22，σX2＝2）；X3服從極值Ⅰ型分布（μX3＝14，σX3＝3.5）。在給定Pω1＝Pω2＝Pω3＝0.99 后，經(jīng)上述方法分別得到X1、X2、X3的高概率區(qū)間為：ω1＝（9.6967，30.3033）、ω2＝（17.1853，27.4766）、ω3＝（7.2575，25.3613），三者的演示如圖1 所示。

圖1 X1、 X2、 X3 的高概率區(qū)間Fig.1 High probability interval of X1，X2，X3

綜合考慮功能函數(shù)中n個基本變量的高概率區(qū)間，組成了X的高概率區(qū)間，在空間中表現(xiàn)為一種多維超方體，記作Ω，也即是有限單元法的計算域，形式為：

例如：將圖1 中的三個高概率區(qū)間ω1、ω2、ω3共同組成Ω時，演示如圖2 所示。

圖2 ω1、 ω2、 ω3 組成的ΩFig.2 Ω consisting of ω1 and ω2 and ω3

Ω內(nèi)包含的概率記作PΩ，PΩ是一個基于Ω內(nèi)對聯(lián)合概率密度函數(shù)的多重積分，由于隨機變量的獨立性，結合式（3）、式（4）、式（9）PΩ表示為：

在計算域Ω 內(nèi)計算結構失效概率時，為達到目標精度，一般要求PΩ包含大部分概率值，即PΩ的取值近似于1。因此，為了計算方便和平衡每個ωi的壓縮效率，對式（11）作如下規(guī)定：

據(jù)此結合式（11）、式（12）得到Pωi與PΩ的關系為：

因此當有n個隨機變量時，可根據(jù)預設誤差限的要求選取合適的PΩ，利用式（13）計算出Pωi，將其帶入式（6 ～8）得到每個ωi的區(qū)間分布，最后根據(jù)式（9）建立計算域Ω，實現(xiàn)對X 分布空間的高效壓縮。

1.2 計算域Ω內(nèi)的有限單元法

壓縮變量分布空間的過程實現(xiàn)了在區(qū)間范圍較小的計算域Ω內(nèi)，富集了滿足計算要求的概率值。然后介紹均勻劃分Ω的具體方法，采用有限單元法需要提出三點要求：（1）每個單元的相關計算需簡便；（2）每個單元內(nèi)部需要有一點表征該單元是否位于失效域。（3）所有單元應毫無間隔且無重疊的排滿計算域。

據(jù)此三點本文采用垂直坐標軸等距劃分的方法，即將Ω中每個隨機變量的高概率區(qū)間均分N份，得到Nn個尺寸相同的單元，外形上此單元可以視為等比縮小Nn倍的Ω。根據(jù)Ω空間內(nèi)方向的不同分別以j、k、…、l代表單元在X1、X2、…、Xn坐標軸正方向上的位置排序（其中j、k、…、l的范圍為1 到N之間的整數(shù)），將此方法得到的單元記為Aj，k，…，l。例如：對圖2 中由ω1、ω2、ω3組成的Ω設置均分份數(shù)N＝10 時，A1，1，5表示X1坐標軸正方向上的第1 個、X2坐標軸正方向上的第1 個、X3坐標軸正方向上的第5 個單元，演示如圖3 所示。

圖3 均分方法的示意Fig.3 Schematic diagram of uniform segmentation

記Aj，k，…，l在Xi方向上的邊長為hi，由上述均分方法可知其表達式為：記Aj，k，…，l的幾何中心點坐標為Oj，k，…，l，結合式（16）得到其表達式為：

規(guī)定以Oj，k，…，l的坐標表征該單元是否位于失效域，因此幾何尺寸hi和中心點坐標Oj，k，…，l為單元的關鍵參數(shù)。

1.3 由單元計算近似失效概率

記單元Aj，k，…，l內(nèi)部包含的概率為Pj，k，…，l，類似地根據(jù)多重積分的性質和隨機變量的獨立性，使用單元的hi和Oj，k，…，l計算Pj，k，…，l的表達式為：

計算域Ω 被失效面分為可靠域與失效域，將Aj，k，…，l的中心點坐標Oj，k，…，l帶入式（1）功能函數(shù)中，若G（Oj，k，…，l）＜0 則表示該單元位于失效域。因此，累加位于失效域單元的概率可以計算近似失效概率。為了實現(xiàn)有效的篩選，引入指示函數(shù)I（x）。如下：

結合式（16）、式（17）在Ω 內(nèi)計算的近似失效概率P′f可以表示為：

綜上所述，使用式（18）計算的P′f由計算域的概率PΩ和均分份數(shù)N的取值確定。

2 算法的誤差分析

2.1 誤差來源

一般來說，可靠度分析的誤差源有兩個［16］：一個是變量統(tǒng)計誤差，即設計變量的模擬誤差，相當于變量數(shù)學模型的可靠性問題；另一個是數(shù)值誤差，是由于建立算法模型所產(chǎn)生的誤差。本文主要討論算法的數(shù)值誤差，記式（2）中失效概率的解析值

本文算法中使用式（18）計算P′f時，數(shù)值誤差來源有兩個。其一為，建立Ω 在高效壓縮X分布空間的同時，忽略了Ω 以外的低概率區(qū)間，產(chǎn)生了數(shù)值誤差。將Ω內(nèi)計算的失效概率解析值記為其與的誤差記作范圍誤差，范圍誤差限記為Pε，Pε是一種由PΩ 決定的誤差，二者的關系為：

其二為，在Ω 內(nèi)使用有限單元法產(chǎn)生的誤差。位于失效面較近的若干單元會出現(xiàn)失效面穿過這些單元，使其一部分處于可靠域，另一部分處于失效域。由于使用單元的中心點來表征其是否失效是一種非此即彼的前提假設，因此產(chǎn)生了算法誤差。圖4 演示了二維變量的情況，按照中心點表征原則，該單元位于可靠域，從而忽略了位于失效域的部分，產(chǎn)生了相應誤差。

圖4 失效面穿過區(qū)域單元的情況Fig.4 The situation of failure face through the regional unit

2.2 誤差分析

研究發(fā)現(xiàn)，第二種誤差與均分份數(shù)N有關。一般N取值越大，與失效面相交的單元包含概率值越小，該誤差也就越小。隨著N增大，在Ω內(nèi)計算的近似失效概率逐漸收斂因此，為得到一定精度下的失效概率，需要不斷地增加均分份數(shù)N，代入式（18）求算新的近似失效概率，數(shù)值逼近為使每次增加均分份數(shù)N計算的近似失效概率穩(wěn)健收斂，本文采用倍增均分份數(shù)的方法進行迭代。即當初次均分份數(shù)為N0時，第m（m＝1，2，…）次迭代中的均分份數(shù)為：

記第m－1 次迭代計算的近似失效概率為第m次將具體化為：

式中：α1，α2…等為0 到9 中的一個數(shù)字，且α1≠0。當?shù)螖?shù)m足夠大時，的數(shù)值充分接近，若其存在下式關系：

與一般近似法、模擬法等其他方法相比，本文算法的一個優(yōu)點在于根據(jù)式（24）給出了準確的誤差范圍。

3 數(shù)值逼近的迭代計算步驟

在上述計算近似失效概率的框架中，PΩ和N為主要參數(shù)：PΩ指導構建計算域Ω的空間范圍，N為計算域Ω內(nèi)基于有限單元法的計算過程提供劃分精度，二者的取值分別對應算法誤差分析中的Pε與Pε。因此該算法從預設誤差限出發(fā)，得到各種參數(shù)，迭代計算滿足精度要求的失效概率。

綜上所述，迭代步驟如下：

（1）預設計算結果的范圍誤差限Pε，并給式（23）中的c賦值，設定迭代誤差限Pε。

（2）根據(jù)設定的Pε利用式（19）計算PΩ。

（3）根據(jù)PΩ，利用式（13）、式（6 ～8）計算Ω的空間分布。

（4）設定初始均分份數(shù)N0。

（5）根據(jù)設定的Pε將式（22）確定為迭代終止條件。

（6）利用式（20）確定與迭代次數(shù)m對應的均分份數(shù)N。

（7）利用式（14 ～18）計算第m次迭代求算

（8）判斷第m次計算的與第m－1 次計算的是否滿足步驟（5）中的迭代終止條件，若不滿足則令m＝m＋1 重復步驟（6 ～7），計算新的近似失效概率。

（9）若第m次計算的滿足步驟（5）中的迭代終止條件，則停止迭代。

4 數(shù)值算例

算例1 選自文獻［18］，算例2 選自文獻［19］，算例3 選自文獻［20］。本文算法使用MT

ALB編程實現(xiàn)計算過程，使用高精度的MCS 作為對照，檢驗新算法可行性的同時，對其計算精度、效率以及穩(wěn)健性進行驗證。

4.1 算例1

結構中一邊長為b的正方形截面軸壓短柱受到的軸壓為P＝1000kN，設短柱的材料強度為fc，其中b和fc服從正態(tài)分布且相互獨立，它們的均值和方差分別為：μb＝300mm，σb＝6mm；μfc＝22N／mm2，σfc＝5N／mm2。計算柱的失效概率。

此算例中結構的功能函數(shù)為：Z＝b2fc－1000。依照計算步驟，本文算法首先預設計算結果的范圍誤差限和迭代誤差限分別為：Pε＝0.5 ×10－6，Pε＝0.5 ×10－6（c＝－6）；設定初始均分份數(shù)N0＝5；開始迭代。最終迭代次數(shù)m為5 次，均分份數(shù)N為80，取用區(qū)域單元的樣本數(shù)為N2＝6400 個，計算結果為本文算法和MCS 計算結果與使用樣本個數(shù)的分布見圖5。工程中MCS被普遍認為是精確解，以其收斂時（模擬1 ×108次）的計算結果作為準確值，即因為Pε＋Pmε，佐證了該算法提出的誤差范圍在邏輯上是自洽的。該算法與MCS的結果對比見表1 所示。由此可見，在處理失效概率較小的算例時，相同精度下本文算法需要的樣本個數(shù)遠小于MCS，計算效率顯著提高。

圖5 算例1 計算結果的分布Fig.5 The distribution of calculation results in example 1

表1 算例1 中本文方法與MCS計算結果對比Tab.1 Comparative Method and MCS calculation results in example 1

4.2 算例2

設結構的極限狀態(tài)方程為：Z＝X2－8100（X1＋其中X1服從正態(tài)分布，X2和X3服從對數(shù)正態(tài)分布，X4服從極值Ⅰ型分布，均值μX＝（60，2000，24，50）T，標準差σX＝（6.0，74.0，1.2，10.0）T。計算結構的失效概率。

該算例計算方法與算例1 相同，初始參數(shù)選取和迭代結果見表2，與MCS 的比較見表3。該算例表明本文算法在保持計算精度和效率的同時，對非正態(tài)隨機變量和高次非線性功能函數(shù)有很好的兼容性。

表2 算例2 中本文算法初始參數(shù)選取和迭代結果Tab.2 The initial parameter selection and iteration of this product algorithm in example 2

表3 算例2 中本文方法與MCS計算結果對比Tab.3 Comparative Method and MCS calculation results in example 3

4.3 算例3

設結構的功能函數(shù)為：

其中，X1和X2服從標準正態(tài)分布且相互獨立，P為該功能函數(shù)的參數(shù)，隨著P取值的增大，極限狀態(tài)方程的非線性程度也隨之提高。對于算例3，選擇不同的P值進行計算，同樣就計算結果將本文方法與MCS 作對比，數(shù)據(jù)列入表4。為使照組的計算結果達到一定精度，MCS法對該算例采用1 ×107次模擬，平均每次模擬結果的耗時為96s。相同精度下本文方法僅需不到1s的時間。

表4 算例3 的計算結果與對比Tab.4 Calculation results and comparison in example 3

通過分析算例3 發(fā)現(xiàn)，本文算法除上述優(yōu)點外，其對非線性程度各異的功能函數(shù)都有較好的適應性，即穩(wěn)健性較強。

5 結論

結構可靠度分析中，計算結果準確、工作量小是研究的主要方向。為解決MCS 面對小概率事件時求解的困難，本文基于有限單元法提出一種新算法。該算法通過壓縮與均勻劃分隨機變量的多維積分空間實現(xiàn)化整為零、逐個分析；經(jīng)數(shù)值逼近計算結構失效概率。數(shù)值算例檢驗表明：

1.算例1 中在－6.5%的相對誤差下，新算法的計算量和計算耗時可減少至MCS 的0.0001倍以下；算例2 和算例3 中則是以－2.3%與小于6%的相對誤差，分別將工作量減少至0.1 與0.01 以下。因此相同計算精度下，本文算法工作量一般可減小到MCS 方法的0.1 以下，由其當失效概率較小時，該算法的優(yōu)勢更明顯。

2.新算法繼承了MCS方法的諸多優(yōu)點，即：對服從任意分布的基本變量、高次非線性、甚至不連續(xù)的功能函數(shù)都可以很好的求解，并且具有很好的穩(wěn)健性。

3.由于該算法事先預設誤差限，進而計算結構的失效概率，因此就計算結果可以得知其準確的誤差范圍。

本文方法也存在諸多不足之處，在面對“維數(shù)災難”時，區(qū)域單元的取樣將較為冗多，計算效率有待改進。但隨著計算機技術的發(fā)展，該算法在高次非線性的大型復雜結構工程的可靠度分析中應有良好的應用前景。