劉旭

河南理工大學(xué)土木工程學(xué)院 焦作454000

引言

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

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

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

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

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

1 有限單元法處理過(guò)程

1.1 建立X空間的計(jì)算域

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

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

將式（5）轉(zhuǎn)變?yōu)閔與xi的關(guān)系，并對(duì)h取最小值得到：

因此在計(jì)算中根據(jù)要求給定Pωi后，即可利用式（6 ～8）求得對(duì)應(yīng)式（3）中ωi的區(qū)間分布范圍。例如：已知X1服從正態(tài)分布（μX1＝20，σX1＝4）；X2服從對(duì)數(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個(gè)基本變量的高概率區(qū)間，組成了X的高概率區(qū)間，在空間中表現(xiàn)為一種多維超方體，記作Ω，也即是有限單元法的計(jì)算域，形式為：

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

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

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

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

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

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

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

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

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

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

記Aj，k，…，l在Xi方向上的邊長(zhǎng)為hi，由上述均分方法可知其表達(dá)式為：記Aj，k，…，l的幾何中心點(diǎn)坐標(biāo)為Oj，k，…，l，結(jié)合式（16）得到其表達(dá)式為：

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

1.3 由單元計(jì)算近似失效概率

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

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

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

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

2 算法的誤差分析

2.1 誤差來(lái)源

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

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

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

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

2.2 誤差分析

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

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

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

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

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

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

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

（1）預(yù)設(shè)計(jì)算結(jié)果的范圍誤差限Pε，并給式（23）中的c賦值，設(shè)定迭代誤差限Pε。

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

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

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

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

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

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

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

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

4 數(shù)值算例

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

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

4.1 算例1

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

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

圖5 算例1 計(jì)算結(jié)果的分布Fig.5 The distribution of calculation results in example 1

表1 算例1 中本文方法與MCS計(jì)算結(jié)果對(duì)比Tab.1 Comparative Method and MCS calculation results in example 1

4.2 算例2

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

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

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

表3 算例2 中本文方法與MCS計(jì)算結(jié)果對(duì)比Tab.3 Comparative Method and MCS calculation results in example 3

4.3 算例3

設(shè)結(jié)構(gòu)的功能函數(shù)為：

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

表4 算例3 的計(jì)算結(jié)果與對(duì)比Tab.4 Calculation results and comparison in example 3

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

5 結(jié)論

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

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

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

3.由于該算法事先預(yù)設(shè)誤差限，進(jìn)而計(jì)算結(jié)構(gòu)的失效概率，因此就計(jì)算結(jié)果可以得知其準(zhǔn)確的誤差范圍。

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