趙恩康，李洪雙

(南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016)

0 引言

結(jié)構(gòu)動(dòng)態(tài)可靠性分析的任務(wù)是量化計(jì)算一個(gè)結(jié)構(gòu)或結(jié)構(gòu)系統(tǒng)在工作了一定的時(shí)間t或者時(shí)段[0,t]后的可靠度。盡管針對(duì)結(jié)構(gòu)動(dòng)態(tài)可靠性分析已經(jīng)有了許多研究工作，但是兼顧高精度和高效率仍然是一個(gè)挑戰(zhàn)。傳統(tǒng)的蒙特卡洛法(monte carlo simulation, MCS)可以用于結(jié)構(gòu)動(dòng)態(tài)可靠性分析，但在實(shí)際應(yīng)用中由于需要大量抽樣計(jì)算，所需的計(jì)算成本太高，阻礙了它在實(shí)際工程中的應(yīng)用?，F(xiàn)有的方差縮減技術(shù)同樣面臨大計(jì)算量困難，例如重要抽樣法[1]。因此結(jié)構(gòu)動(dòng)態(tài)可靠性分析領(lǐng)域目前主要采取近似方法，主要包括首次超越法[2-3]和極值法[4-6]。首次超越法認(rèn)為極限狀態(tài)函數(shù)超過(guò)安全閾值的上界或低于其下界，結(jié)構(gòu)就會(huì)發(fā)生失效。通過(guò)假設(shè)不同“超越事件”之間是相互獨(dú)立的，結(jié)構(gòu)動(dòng)態(tài)可靠度就可以通過(guò)在時(shí)間域上對(duì)穿越率進(jìn)行積分來(lái)獲得。大多數(shù)首次超越法的局限性在于忽略了超越時(shí)間之間的相關(guān)性及線性近似帶來(lái)的誤差。極值法的基本思想是考慮極限狀態(tài)函數(shù)在目標(biāo)時(shí)段內(nèi)的極值，如果極值大于給定的閾值，則結(jié)構(gòu)發(fā)生失效。目前，基于極值思想已經(jīng)發(fā)展出多種結(jié)構(gòu)動(dòng)態(tài)可靠性分析方法，例如，WANG提出了一種嵌套極值響應(yīng)面法[6]，DU和HU提出了一種混合高效全局優(yōu)化法[7]，CHEN和WANG則提出了一種自適應(yīng)極值響應(yīng)面法[8]。但是現(xiàn)有極值法的性能依賴于代理模型對(duì)極限狀態(tài)函數(shù)的逼近精度和范圍。若不能高精度地?cái)M合結(jié)構(gòu)極限狀態(tài)函數(shù)及其極值，相應(yīng)的動(dòng)態(tài)可靠性分析結(jié)果是不可信的。

本文基于極值思想，提出一種結(jié)構(gòu)動(dòng)態(tài)可靠性分析的最大熵方法，可以直接處理一般形式的動(dòng)態(tài)極限狀態(tài)函數(shù)。首先采用拓展最優(yōu)線性估計(jì)(expansion optimal linear estimation, EOLE))[9]方法，將極限狀態(tài)函數(shù)中的輸入隨機(jī)過(guò)程離散化為一系列隨機(jī)變量之和的形式。對(duì)所研究結(jié)構(gòu)進(jìn)行結(jié)構(gòu)分析，計(jì)算得出目標(biāo)時(shí)間區(qū)間內(nèi)極限狀態(tài)函數(shù)的極值；進(jìn)而利用極值的前四階統(tǒng)計(jì)矩和最大熵原理擬合極值分布?；谧畲箪胤植?，求解結(jié)構(gòu)動(dòng)態(tài)可靠度。文中利用工程算例驗(yàn)證了所提方法的有效性。

1 結(jié)構(gòu)動(dòng)態(tài)可靠性分析

1.1 動(dòng)態(tài)可靠性定義

結(jié)構(gòu)動(dòng)態(tài)可靠性分析的極限狀態(tài)函數(shù)的一般形式為

g(X,Y(t),t)

(1)

其中X=[X1,X2,…,Xn] 是n個(gè)輸入隨機(jī)變量的向量；Y(t)=[Y1(t),Y2(t),…,Ym(t)] 是m個(gè)與時(shí)間有關(guān)的輸入隨機(jī)過(guò)程。如果在時(shí)間區(qū)間[t1,t2] 內(nèi)存在某一時(shí)間節(jié)點(diǎn)t，極限狀態(tài)函數(shù)值超出安全閾值b，結(jié)構(gòu)發(fā)生失效，即有

?t∈[t1,t2],g(X,Y(t),t)>b

(2)

結(jié)構(gòu)失效概率可以定義為

Pf(t1,t2)=Pr(g(X,Y(t),t)>b,?t∈[t1,t2])

(3)

相應(yīng)的動(dòng)態(tài)可靠度為

R(t1,t2)=1-Pf(t1,t2)

(4)

由于極限狀態(tài)函數(shù)的復(fù)雜性以及時(shí)間參數(shù)的存在，得到式(3)的解析解是非常困難的。蒙特卡洛法能夠得到精確的結(jié)果，但由于需要大量樣本來(lái)評(píng)估極限狀態(tài)函數(shù)，計(jì)算成本極高。因此本文研究基于極值思想的方法。

1.2 極值法原理

極值法關(guān)注極限狀態(tài)函數(shù)一段時(shí)間內(nèi)的極值函數(shù)或者極值分布。那么基于極值思想，式(3)中的結(jié)構(gòu)失效概率可以表示為

Pf(t1,t2)=Pr(gmax>b)

(5)

顯然，式(5)中不含有時(shí)間參數(shù)，這就將時(shí)變問(wèn)題轉(zhuǎn)化為時(shí)不變問(wèn)題。但是由于極限狀態(tài)函數(shù)的復(fù)雜性，在實(shí)際工程中解析地得到其極值往往非常困難，所以目前已有研究工作主要是使用代理模型方法和近似方法得到結(jié)構(gòu)響應(yīng)的極值或者極值分布。

2 動(dòng)態(tài)可靠性分析的最大熵方法

2.1 所提方法

為了解決現(xiàn)有方法的一些缺點(diǎn)，基于最大熵原理，本文提出一種新的結(jié)構(gòu)動(dòng)態(tài)可靠性分析方法。該方法的實(shí)施步驟如下：

1) 采用EOLE方法對(duì)輸入隨機(jī)過(guò)程進(jìn)行離散化；

2) 采用拉丁質(zhì)心Voronoi網(wǎng)格抽樣生成高質(zhì)量抽樣樣本；

3) 計(jì)算每條樣本軌跡的極值，得到響應(yīng)極值的樣本；

4) 計(jì)算極值的前四階矩，分別是均值、標(biāo)準(zhǔn)差、偏度和峰度；

5) 由極值的前四階矩計(jì)算Lagrange乘子，然后擬合極值的最大熵分布；

6) 計(jì)算動(dòng)態(tài)可靠度。

本文提出的動(dòng)態(tài)可靠性分析的最大熵方法流程如圖1所示。

圖1 最大熵方法流程圖

如圖2所示，陰影部分的面積即為結(jié)構(gòu)失效概率。

圖2 極值分布計(jì)算失效概率示意圖

下面針對(duì)步驟中關(guān)鍵理論和技術(shù)進(jìn)行詳細(xì)討論。

2.2 隨機(jī)過(guò)程的離散

在實(shí)際工程應(yīng)用中，許多結(jié)構(gòu)參數(shù)都具有不確定性，例如材料屬性。大多數(shù)情況下，這種不確定性都可以通過(guò)隨機(jī)變量來(lái)表征。但是仍然有一些參數(shù)隨結(jié)構(gòu)服役時(shí)間的累積呈現(xiàn)出時(shí)變隨機(jī)特性，例如結(jié)構(gòu)上的隨機(jī)載荷，需要用隨機(jī)過(guò)程來(lái)表征。若想對(duì)隨機(jī)過(guò)程Y(t)進(jìn)行抽樣計(jì)算，在此之前需要對(duì)該隨機(jī)過(guò)程Y(t)進(jìn)行離散化處理。本文討論的隨機(jī)過(guò)程限定為高斯隨機(jī)過(guò)程，該過(guò)程的均值函數(shù)、標(biāo)準(zhǔn)差函數(shù)和自相關(guān)系數(shù)函數(shù)分別為μY(t)、σY(t) 和ρY(t1,t2)。非高斯隨機(jī)過(guò)程可以通過(guò)轉(zhuǎn)換方法轉(zhuǎn)化為高斯隨機(jī)過(guò)程。

EOLE方法首先將時(shí)間[0,T]離散成s個(gè)離散時(shí)間點(diǎn)ti,i=1,2,…,s?？紤]到計(jì)算上的簡(jiǎn)便，時(shí)間步Δt一般取定值。在時(shí)間點(diǎn)處多維隨機(jī)變量的相關(guān)矩陣可以表示為

(6)

(7)

式中：Ui為相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量；ρY(t)=ρY(t,ti),i=1,2,…,s。

式(7)表明：在每一個(gè)時(shí)間點(diǎn)ti處，隨機(jī)過(guò)程Y(ti)可以表示為一組標(biāo)準(zhǔn)正態(tài)隨機(jī)變量Ui的線性和?；谝陨辖Y(jié)論，結(jié)構(gòu)動(dòng)態(tài)可靠性的極限狀態(tài)函數(shù)可以轉(zhuǎn)化為g(X,Y(t),t)=g(X,U,t) ，使得極限狀態(tài)函數(shù)中不再隱式含有時(shí)間參數(shù)t。

2.3 拉丁質(zhì)心Voronoi網(wǎng)格抽樣

抽樣樣本點(diǎn)的質(zhì)量會(huì)影響擬合的精度。目前常用的拉丁超立方抽樣(latin hypercube sampling, LHS)是一種全空間非重疊填充的隨機(jī)抽樣方法，能夠兼顧可行域與不可行域，但難以保證抽樣的空間均勻性；BURKARDT[10]基于Voronoi網(wǎng)格提出了質(zhì)心Voronoi網(wǎng)格抽樣方法 (centroidal voronoi tessellation, CVT)。它將整個(gè)抽樣空間劃分為指定數(shù)目的Voronoi區(qū)域，Voronoi區(qū)域的質(zhì)心即為抽樣樣本點(diǎn)，具有均勻性與穩(wěn)定性的特點(diǎn)。

ROMERO[11]等人結(jié)合LHS方法和CVT抽樣方法，提出了一種新的抽樣方法，稱為拉丁質(zhì)心Voronoi網(wǎng)格(‘latinized’ centroidal voronoi tessellation，LCVT)抽樣方法。對(duì)質(zhì)心Voronoi網(wǎng)格進(jìn)行特殊處理，使其具有拉丁超立方體的性質(zhì)，此處理過(guò)程稱為拉丁化。拉丁質(zhì)心Voronoi網(wǎng)格抽樣的具體實(shí)施步驟如下：

(8)

式中Uji表示在第j維空間內(nèi)均勻分布的隨機(jī)數(shù)。

(9)

對(duì)質(zhì)心Voronoi網(wǎng)格進(jìn)行拉丁化，獲得的新網(wǎng)格稱為拉丁質(zhì)心Voronoi網(wǎng)格。LCVT抽樣方法兼具CVT方法的均勻性與LHS方法一定的隨機(jī)性，能更好地反映出隨機(jī)變量在設(shè)計(jì)空間的分布情況。

基于極值思想的結(jié)構(gòu)動(dòng)態(tài)可靠性分析方法的關(guān)鍵步驟是準(zhǔn)確地計(jì)算極限狀態(tài)函數(shù)g(X,Y(t),t)的極值。獲得高質(zhì)量抽樣點(diǎn)后，帶入結(jié)構(gòu)分析模型，計(jì)算每個(gè)抽樣點(diǎn)對(duì)應(yīng)的軌跡，并獲得結(jié)構(gòu)響應(yīng)在感興趣時(shí)段內(nèi)的極值。進(jìn)一步計(jì)算結(jié)構(gòu)響應(yīng)極值的前四階統(tǒng)計(jì)矩，為擬合最大熵分布提供約束條件。

2.4 最大熵分布

1948年，SHANNON[12]將熱力學(xué)熵的概念引入到信息科學(xué)中，來(lái)定量描述一個(gè)隨機(jī)事件的不確定性或信息量?；谛畔㈧氐母拍?，JAYNE[13]在1957年提出了最大熵原理。最大熵原理可以表述為：在所有滿足給定的約束條件的概率密度函數(shù)中, 使信息熵最大的概率密度函數(shù)就是最佳，即偏差最小的概率密度函數(shù)。

若X為連續(xù)型隨機(jī)變量，其概率密度函數(shù)為f(x)，則其信息熵HX可定義為

(10)

根據(jù)最大熵原理，應(yīng)該求使信息熵HX最大的概率密度函數(shù)f(x)。通常將隨機(jī)變量的前n階原點(diǎn)矩μi(i=1,…,n)作為約束條件，相應(yīng)的優(yōu)化問(wèn)題為：

(11)

式(11)是一個(gè)等式約束極值問(wèn)題，可以用拉格朗日乘數(shù)法求解，經(jīng)過(guò)簡(jiǎn)單推導(dǎo)，可以得到最大熵分布

(12)

式中λi(i=1,…,n)為L(zhǎng)agrange乘子。根據(jù)概率密度函數(shù)的數(shù)學(xué)定義，求得標(biāo)準(zhǔn)化因子：

在估計(jì)最大熵分布時(shí)，隨著使用的統(tǒng)計(jì)矩μi階次增大，對(duì)樣本中的信息利用也更加充分，最大熵分布擬合的精度也會(huì)提高。文中采用極值數(shù)據(jù)的前四階矩作為約束條件擬合極值分布，它們分別是均值、標(biāo)準(zhǔn)差、偏度和峰度。

由于本文中使用了結(jié)構(gòu)響應(yīng)極值的前四階矩?cái)M合最大熵分布，無(wú)法獲得最大熵分布的顯示解，此時(shí)需要采用諸如標(biāo)準(zhǔn)牛頓法[14]的數(shù)值算法求解Lagrange乘子。但是標(biāo)準(zhǔn)牛頓法需要求解n個(gè)帶有積分的等式，過(guò)程繁瑣，計(jì)算量大，效率較低，因此文獻(xiàn)中使用一種無(wú)約束最小化方法求解Lagrange乘子，基本過(guò)程如下：

定義一個(gè)以λi(i=1,…,n)為自變量的勢(shì)函數(shù)Q(λ1,…,λn)

(13)

勢(shì)函數(shù)Q(λ1,…,λn)對(duì)自變量λi的梯度為：

(14)

基于式(14)，進(jìn)一步計(jì)算Q(λ1,…,λn) 的Hessian矩陣元素：

(15)

3 算例分析

水動(dòng)力渦輪機(jī)可以將水的動(dòng)能轉(zhuǎn)化為葉片的轉(zhuǎn)動(dòng)機(jī)械能。此算例為時(shí)變河流載荷下的水動(dòng)力渦輪葉片[18]。渦輪葉片的簡(jiǎn)化橫截面如圖3所示。

圖3 渦輪葉片根部橫截面

為了表征河流流速的時(shí)變特性，河流流速采用一個(gè)高斯隨機(jī)過(guò)程v(t)表示，其均值μv(t)、標(biāo)準(zhǔn)差σv(t)和自相關(guān)函數(shù)ρv(t1,t2)分別為：

(16)

(17)

ρv(t1,t2)=cos(2π(t2-t1))

(18)

式中系數(shù)a、b和c的值如下：

在河流流速v(t)作用下，葉片根部的彎矩可以表征為

(19)

式中：ρ=1×103kg/m3為水的密度；Cm=0.3422為力矩系數(shù)。

該葉片的極限狀態(tài)函數(shù)定義為

(20)

(21)

表1 渦輪葉片中隨機(jī)變量

在[0,12]個(gè)月區(qū)間上計(jì)算渦輪葉片的失效概率，通過(guò)LCVT方法生成500條樣本曲線，并與樣本數(shù)目為106的MCS計(jì)算結(jié)果進(jìn)行對(duì)比，失效概率結(jié)果如表2和圖4所示，[0,12]個(gè)月區(qū)間極值的概率密度函數(shù)曲線如圖5所示。

表2 渦輪葉片失效概率

圖4 渦輪葉片失效概率

圖5 [0,12]個(gè)月極值概率密度函數(shù)曲線

由表2和圖4可知：隨著時(shí)間的推移，渦輪葉片的失效概率逐漸增大，從第7個(gè)月開始保持不變，這是由于極限狀態(tài)函數(shù)的極值落在[6,7]個(gè)月之間。在整個(gè)時(shí)間區(qū)間內(nèi)，最大熵方法的計(jì)算結(jié)果都和MCS計(jì)算結(jié)果十分接近，并且相比MCS需要106次抽樣，最大熵方法只需要500條樣本曲線，說(shuō)明最大熵方法兼具較好的計(jì)算精度和計(jì)算效率。

4 結(jié)語(yǔ)

本文基于極值思想，提出了一種結(jié)構(gòu)動(dòng)態(tài)可靠性分析的最大熵方法，能夠處理形如g(X,Y(t),t)的一般形式的動(dòng)態(tài)可靠性極限狀態(tài)函數(shù)。該方法首先將隨機(jī)過(guò)程離散為一組隨機(jī)變量，再對(duì)極限狀態(tài)函數(shù)中的隨機(jī)變量進(jìn)行抽樣并計(jì)算樣本點(diǎn)處的極值，最后通過(guò)樣本極值數(shù)據(jù)利用最大熵原理擬合極值分布并計(jì)算動(dòng)態(tài)可靠度。從文中算例的計(jì)算結(jié)果可以看出：該方法在精度上接近蒙特卡洛方法且兼具較高計(jì)算效率，有一定的工程應(yīng)用價(jià)值。