基于四階矩法的CRDM承壓殼體共因失效可靠性分析

陳 鵬，朱翊洲，謝永誠

(上海核工程研究設(shè)計院有限公司，上海 200233)

0 引言

目前，在設(shè)計核電主設(shè)備時，通常采用分析法設(shè)計進(jìn)行分析，并基于ASME BPVC 或RCC-M等規(guī)范中的相關(guān)設(shè)計準(zhǔn)則進(jìn)行結(jié)構(gòu)完整性評定，盡管評定準(zhǔn)則中會考慮一定的安全系數(shù)，但確定性分析法忽視了結(jié)構(gòu)和運(yùn)行環(huán)境中存在的客觀不確定因素，如加工導(dǎo)致的材料力學(xué)性能和結(jié)構(gòu)尺寸的隨機(jī)性，不同工況下溫度、壓力以及機(jī)械載荷波動的不確定性等。為了研究這些不確定因素對設(shè)備安全性的影響，并為設(shè)備維修和優(yōu)化設(shè)計提供理論指導(dǎo)，有必要在設(shè)計過程中充分考慮潛在失效模式和隨機(jī)因素、開展概率可靠性研究。

近幾十年來，國內(nèi)外學(xué)者對概率可靠性理論和應(yīng)用開展了大量的研究。其中，基于隨機(jī)采樣的Monte-Carlo方法[2-3]被認(rèn)為是最精確的方法，但當(dāng)可靠性模型規(guī)模較大時，該方法的求解效率難以滿足對大樣本的隨機(jī)采樣[1-2]；SONG等[4]基于Markov鏈技術(shù)和重要采樣法，并通過引入失效域子集和條件失效概率提高了Monte-Carlo法在高可靠性問題中的求解效率，但對于工程中的復(fù)雜結(jié)構(gòu)和非正態(tài)分布可靠性問題仍然很難適用。

為了滿足工程需求，學(xué)者們提出了以一次二階矩法[5]為代表的一階可靠性方法(FORM)，并基于對驗算點求解、非正態(tài)分布參數(shù)標(biāo)準(zhǔn)化及功能函數(shù)線性化等關(guān)鍵問題的研究進(jìn)行了修正，但由于FORM中對功能函數(shù)作了正態(tài)分布假設(shè)，導(dǎo)致該方法在求解基本隨機(jī)變量和功能函數(shù)偏度較大的可靠性問題時精度較低。在此基礎(chǔ)上，學(xué)者們將FORM拓展到了高階矩方法。楊周等[6-7]考慮非正態(tài)分布的隨機(jī)變量，用隨機(jī)攝動理論計算功能函數(shù)的前四階中心矩，基于Edgeworth級數(shù)法給出功能函數(shù)的假設(shè)分布，并計算可靠度和靈敏度，該方法在功能函數(shù)和基本變量隨機(jī)部分所做的攝動假設(shè)限制了其在高非線性功能函數(shù)和離散性較大的可靠性問題中的應(yīng)用；李云貴等[8-9]基于最大熵原理，以功能函數(shù)前四階原點矩作為約束條件，通過Lagrange乘子法計算功能函數(shù)的概率分布和可靠度，但對于復(fù)雜結(jié)構(gòu)問題，優(yōu)化求解過程難以收斂；在四階矩理論中，Pearson分布族也被廣泛應(yīng)用于可靠性問題的求解[10]，根據(jù)功能函數(shù)的前四階中心矩，可以直接計算對應(yīng)的分布函數(shù)，然而在實際應(yīng)用中，有12種類型的概率密度函數(shù)，且參數(shù)計算需要通過數(shù)值積分，難以在工程中應(yīng)用；此外，ZHAO等[11]基于三次多項式近似，提出以四階矩為擬合目標(biāo)的標(biāo)準(zhǔn)正態(tài)化技術(shù)，在標(biāo)準(zhǔn)正態(tài)空間中計算結(jié)構(gòu)可靠度，并證明在較大的偏度范圍內(nèi)均有很高的計算精度。

針對核電設(shè)備中失效形式復(fù)雜、非正態(tài)分布變量較多的結(jié)構(gòu)可靠性問題，本文基于四階矩標(biāo)準(zhǔn)正態(tài)化技術(shù)，提出一種適用于復(fù)雜核電設(shè)備多失效模式的結(jié)構(gòu)可靠性分析方法?；谟邢拊托蛄许憫?yīng)面法建立考慮共因失效的多失效模式可靠性功能函數(shù)，用點估計法計算功能函數(shù)四階矩，在此基礎(chǔ)上通過四階矩標(biāo)準(zhǔn)正態(tài)化變換得到標(biāo)準(zhǔn)正態(tài)空間中的可靠度指標(biāo)。通過對控制棒驅(qū)動機(jī)構(gòu)(CRDM)承壓殼體結(jié)構(gòu)多失效模式的可靠度分析，驗證本文方法的適用性，并與Monte-Carlo方法結(jié)果的對比，證明該方法的計算精度。

1 基于響應(yīng)面法的多失效模式功能函數(shù)

對于工程中的結(jié)構(gòu)可靠性問題，描述結(jié)構(gòu)失效的功能函數(shù)通?？梢员硎緸椋?/p>

g(x)=S-σ(x)

(1)

式中，S為基于相關(guān)強(qiáng)度理論和設(shè)計規(guī)范的設(shè)計應(yīng)力強(qiáng)度值；x為影響應(yīng)力分量σ的所有隨機(jī)變量，當(dāng)S也是隨機(jī)變量時，功能函數(shù)對應(yīng)的基本隨機(jī)變量x應(yīng)包含S。

結(jié)構(gòu)可靠度的定義可以用功能函數(shù)描述為：

R=P[g(x)>0]

(2)

對于工程中的大多數(shù)結(jié)構(gòu)的可靠性問題，式(1)所定義的目標(biāo)應(yīng)力分量通常只能通過有限元計算，無法給出顯式的功能函數(shù)，因此難以通過FORM或高階矩方法計算結(jié)構(gòu)可靠度。本文基于響應(yīng)面法(RSM)，用不含交叉項的二次多項式h(x)對隱式功能函數(shù)g(x)進(jìn)行擬合：

(3)

式中，xi為基本隨機(jī)變量(i=1,…,n)；a,bi,ci(i=1,…,n)為2n+1個待定系數(shù)，定義為：υ=[ab1bn…c1…cn]T。

(4)

考慮到基本隨機(jī)變量在每個試驗點的分布概率不同，因此文中通過加權(quán)最小二乘法求解待定系數(shù)向量：

υ=(HTWH)-1HTWg

(5)

其中，W為試驗點權(quán)重矩陣,為對角陣。

(6)

在核電設(shè)備的設(shè)計分析中，通常需要對結(jié)構(gòu)不同位置和不同工況下的應(yīng)力分量進(jìn)行評定。在結(jié)構(gòu)可靠性分析中，同樣需要考慮多失效模式的結(jié)構(gòu)可靠性問題。由于共因失效在功能函數(shù)描述上的復(fù)雜性，以往在多失效模式的結(jié)構(gòu)可靠性分析中，學(xué)者們通常保守地忽略了共因失效的影響，假定各失效模式相互獨立，而直接將各點或各失效模式進(jìn)行串聯(lián)系統(tǒng)處理，即結(jié)構(gòu)可靠度可以表示為：

RInd=R1R2…Rk

(7)

式中，Ri(i=1,…,k)為第i種失效模式單獨計算的可靠度。

張義民等[13]在多失效模式機(jī)械零部件可靠性設(shè)計的結(jié)構(gòu)分析中，也提出了另一種近似可靠度，即考慮最危險環(huán)節(jié)的多失效模式可靠度：

RMin=Min[R1,R2,…,Rk]

(8)

這兩種方法得到的可靠度實際上是兩種假設(shè)的極端，即各失效模式完全相互獨立和完全相關(guān)。

本文在考慮多失效模式結(jié)構(gòu)可靠性問題時，假定各失效模式對應(yīng)功能函數(shù)為gi(x)，i=1,2,…,k，則結(jié)構(gòu)功能函數(shù)可以表示為：

g(x)=min[g1(x),g2(x),…,gk(x)]

(9)

即結(jié)構(gòu)中功能函數(shù)值最小的失效點都可靠時，整個結(jié)構(gòu)才算可靠。顯然，在構(gòu)建功能函數(shù)時，考慮多失效模式比單獨計算各失效模式可靠度之后再考慮相關(guān)性更加合理。由于式(9)所描述的功能函數(shù)是不連續(xù)的，因此在傳統(tǒng)基于功能函數(shù)Taylor展開的可靠性方法中難以適用。本文基于響應(yīng)面法，在有限元基礎(chǔ)上擬合得到各失效模式功能函數(shù)gi(x)的近似功能函數(shù)hi(x)，并以如下功能函數(shù)來描述多失效模式的結(jié)構(gòu)可靠性：

g(x)=min[h1(x),h2(x),…,hk(x)]

(10)

2 可靠性功能函數(shù)的四階矩描述

在結(jié)構(gòu)可靠性問題中，首先認(rèn)為表征結(jié)構(gòu)失效的功能函數(shù)也是一個隨機(jī)變量，結(jié)構(gòu)失效概率或可靠度由功能函數(shù)的分布信息確定，實際工程中，由于功能函數(shù)往往比較復(fù)雜，且多數(shù)情況下基本隨機(jī)變量不能給出具體分布形式，因此得到功能函數(shù)的概率密度函數(shù)(PDF)或累積分布函數(shù)(CDF)通常是不現(xiàn)實的。在此基礎(chǔ)上，矩法被廣泛應(yīng)用于功能函數(shù)的隨機(jī)性描述和可靠性分析中。

當(dāng)給定基本隨機(jī)變量的有限階中心矩時，功能函數(shù)的中心矩通常有兩種求解方式。其一是基于隨機(jī)攝動理論和功能函數(shù)一階Taylor展開的攝動法；其二是標(biāo)準(zhǔn)正態(tài)空間的點估計法。

2.1 隨機(jī)攝動法

隨機(jī)攝動方法假定所有隨機(jī)變量x和功能函數(shù)g(x)都可以表示為確定部分和隨機(jī)部分，且隨機(jī)部分遠(yuǎn)小于確定部分：

x=xd+εxp

(11)

g(x)=gd(x)+εgp(x)

(12)

其中，下標(biāo)d表示變量或函數(shù)的確定部分，下標(biāo)p表示隨機(jī)部分，ε表示一個小量。

當(dāng)隨機(jī)變量的隨機(jī)部分比確定部分小得多時：

(13)

則該功能函數(shù)的前四階矩可以寫為：

(14)

(15)

(16)

(17)

可以看出，攝動法需要知道功能函數(shù)在隨機(jī)變量均值處的偏導(dǎo)數(shù)，此外，功能函數(shù)的各階中心矩完全與基本隨機(jī)變量對應(yīng)階中心矩相關(guān)，對于非線性功能函數(shù)顯然會帶來一定的誤差。

2.2 點估計法

為了提高適用性和求解精度，本文采用另一種可靠性功能函數(shù)四階矩求解方法，即點估計法，采用一系列標(biāo)準(zhǔn)正態(tài)空間的估計點和權(quán)函數(shù)近似計算功能函數(shù)的各階中心矩。

當(dāng)極限狀態(tài)函數(shù)g=g(x1)只含有一個服從任意分布的基本隨機(jī)變量時，g的均值和k階中心矩為：

(18)

(20)

式中，T-1(u)為隨機(jī)變量x1到標(biāo)準(zhǔn)正態(tài)分布變量u的Rosenblatt變換的逆變換過程[14]；uI，PI為估計點和對應(yīng)的權(quán)值，可以通過對權(quán)函數(shù)為exp(-x2)的Hermite迭代得到[15]。

對于7點估計法，估計點和權(quán)值分別為：

u0=0，u1=-u2=1.154 405 4，u3=-u4=2.366 759 4，u5=-u6=3.750 439 7；P0=16/35,P1=P2=0.240 123 3，P3=P4=3.075 71×10-2，P5=P6=5.482 69×10-4。

當(dāng)考慮n個基本隨機(jī)變量時，功能函數(shù)g(x)=g(x1,x2,…,xn)可以近似表示為：

(21)

則g(x)的前四階中心矩可以近似表示為：

(22)

(23)

(24)

(25)

各單變量函數(shù)的前四階中心矩可根據(jù)式(18)～(20)得到?？梢钥闯觯命c估計法計算功能函數(shù)的前四階中心矩不需要功能函數(shù)對基本變量的偏導(dǎo)數(shù)，因此更適用于復(fù)雜功能函數(shù)形式。此外，ZHAO等[14]指出，點估計法在估計功能函數(shù)高階矩時需要更多的估計點才能保證精度，因此對于非正態(tài)性很強(qiáng)的基本隨機(jī)變量，可能需要增加估計點個數(shù)。

本文基于攝動法四階矩計算的可靠度結(jié)果,以后綴“_SD”表示，基于點估計法四階矩計算的可靠度結(jié)果，以后綴“_PE”表示。

3 基于功能函數(shù)四階矩的可靠度

當(dāng)無法給出結(jié)構(gòu)可靠度功能函數(shù)的概率密度函數(shù)，而只能得到前四階中心矩時，要計算結(jié)構(gòu)可靠度和可靠度指標(biāo)，在統(tǒng)計學(xué)中需要通過有限階中心矩對功能函數(shù)及其標(biāo)準(zhǔn)化正態(tài)變換形式做適當(dāng)?shù)暮喕?/p>

首先，定義功能函數(shù)的標(biāo)準(zhǔn)化形式為：

(26)

一種較為常用的簡化形式是采用Edgeworth級數(shù)法，將功能函數(shù)展開為標(biāo)準(zhǔn)正態(tài)分布函數(shù)的多項式形式：

(27)

式中，Hi(·)為i階Hermite多項式；y為標(biāo)準(zhǔn)化功能函數(shù)zu；α3g，α4g為功能函數(shù)的偏度和峰度，α3g=M3g/(M2g)3/2,α4g=M4g/(M2g)2；Φ(·)和φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)和概率密度函數(shù)。

用Edgeworth級數(shù)法可以將結(jié)構(gòu)可靠度近似表示為：

REW=P(g(x)>0)=F(βsm)

(28)

可以看出，Edgeworth級數(shù)法在擬合功能函數(shù)標(biāo)準(zhǔn)化變量的累積分布時作了尾部截斷，試圖用有限階中心矩(如四階矩)直接擬合變量的分布函數(shù)，這也導(dǎo)致了式(27)在很多情況下并不符合一個隨機(jī)變量的累積分布特征。在計算結(jié)構(gòu)可靠度時，可能會出現(xiàn)R>1的情況，對于這種情況，可以采用以下公式進(jìn)行可靠度修正：

(29)

在關(guān)于四階矩結(jié)構(gòu)可靠性的研究中，ZHAO等[16]提出了另一種功能函數(shù)的標(biāo)準(zhǔn)正態(tài)化方法，即高階矩標(biāo)準(zhǔn)化方法(HOMST)，對于標(biāo)準(zhǔn)化后功能函數(shù)分布變量zu，假定存在如下標(biāo)準(zhǔn)正態(tài)變換：

(30)

式中，aj為待定系數(shù)，通過計算Sx(zu)的有限階中心矩，使其分別等于標(biāo)準(zhǔn)正態(tài)分布u的對應(yīng)階中心矩，可以得到各項系數(shù)。

當(dāng)考慮u的三階矩時，將功能函數(shù)的標(biāo)準(zhǔn)化變量zu通過如下關(guān)系式表示為標(biāo)準(zhǔn)正態(tài)變量：

(31)

在標(biāo)準(zhǔn)正態(tài)空間中，四階矩可靠度指標(biāo)可以表示為：

(32)

基于四階矩可靠度指標(biāo)的結(jié)構(gòu)可靠度為：

RFM=Φ(βfm)

(33)

在本文的研究中，首先基于有限元和響應(yīng)面法，根據(jù)結(jié)構(gòu)的潛在失效模式分別建立各失效點的功能函數(shù)，并得到結(jié)構(gòu)的多失效模式功能函數(shù)(見式(10))，采用點估計法計算功能函數(shù)的前四階中心矩，最后用式(32)(33)得到可靠度指標(biāo)和可靠度。

4 CRDM承壓殼體結(jié)構(gòu)可靠性分析

控制棒驅(qū)動機(jī)構(gòu)(CRDM)是核反應(yīng)堆中的關(guān)鍵設(shè)備，某磁力提升型CRDM結(jié)構(gòu)如圖1所示，主要由棒行程殼體、棒位探測器、鉤爪部件、磁軛線圈、驅(qū)動桿部件和鉤爪殼體結(jié)構(gòu)組成。其中，棒行程殼體和鉤爪殼體所組成的承壓殼體結(jié)構(gòu)不僅為CRDM提供支撐，還是整個反應(yīng)堆一回路的壓力邊界，其可靠性對反應(yīng)堆的安全運(yùn)行至關(guān)重要。

圖1 CRDM承壓殼體結(jié)構(gòu)示意

本文研究承壓殼體在設(shè)計工況下的結(jié)構(gòu)可靠性。在Ansys中建立軸對稱模型，其中局部結(jié)構(gòu)如圖2所示。分別考慮承壓殼體與頂蓋貫穿件異種金屬焊縫處的強(qiáng)度失效(評定截面ASN1)、鉤爪殼體中部強(qiáng)度失效(評定截面ASN2)、棒行程殼體和鉤爪殼體Canopy焊縫處的強(qiáng)度失效(評定截面ASN3)。失效準(zhǔn)則基于ASME BPVC第三卷NB分卷，本文考慮設(shè)計工況下的一次薄膜應(yīng)力的規(guī)范限制。

圖2 CRDM承壓殼體有限元模型和評定截面

為了研究3個失效模式在設(shè)計工況下的一次薄膜應(yīng)力強(qiáng)度的概率可靠性，考慮設(shè)計內(nèi)壓P、步躍沖擊載荷F、材料設(shè)計應(yīng)力強(qiáng)度Sm三個變量為基本隨機(jī)變量，則3個失效模式一次薄膜應(yīng)力可靠度的功能函數(shù)可以表示為：

gi(x)=Sm-σi(P,F) (i=1,2,3)

(34)

其中，x=[SmPF]T為基本隨機(jī)變量，并且假定分別服從雙參數(shù)Weibull分布、對數(shù)正態(tài)分布、正態(tài)分布，分布參數(shù)如表1所示。

表1 各基本隨機(jī)變量分布參數(shù)

需要說明的是，設(shè)計內(nèi)壓P和步躍沖擊載荷F的均值取自CRDM設(shè)計規(guī)范書。設(shè)計應(yīng)力強(qiáng)度Sm的均值根據(jù)承壓殼體材料牌號，取ASME規(guī)范中SA-182 F304LN在設(shè)計溫度下的設(shè)計應(yīng)力強(qiáng)度值，各參數(shù)變異系數(shù)均取0.1?？紤]到商密要求，表1中的各分布參數(shù)與真實值之間均有一定的偏差。

為了得到顯式的功能函數(shù)，根據(jù)一次應(yīng)力有限元模型和響應(yīng)面法可以分別擬合得到3個失效模式的近似功能函數(shù)hi(x),i=1,2,3。在實際運(yùn)行過程中，三種失效模式都與基本隨機(jī)變量x相關(guān)，即存在共因失效關(guān)系，則承壓殼體共因失效模式的可靠度功能函數(shù)可以表示為：

g(x)=min[h1(x),h2(x),h3(x)]

(35)

表2列出本文點估計法所計算的功能函數(shù)前四階中心矩，對應(yīng)的功能函數(shù)偏度為-0.384，峰度為3.254，可以看出，當(dāng)忽略三階及以上中心矩時，用傳統(tǒng)二階矩方法得到的可靠度將會偏大。

表2 功能函數(shù)前四階中心矩

為了研究幾種常用可靠度計算方法的精度，并驗證串聯(lián)系統(tǒng)假設(shè)在共因失效問題中的保守性，表3列出了5種可靠度計算結(jié)果。

表3 CRDM承壓殼體多失效模式可靠度

表3中，MC法為基于功能函數(shù)式(35)的106次Monte-Carlo仿真結(jié)果，假定為精確值并作為本文近似方法的驗證對象；EW法為采用Edgeworth級數(shù)法，根據(jù)式(27)(28)得到的四階矩可靠度；SM法為根據(jù)可靠度指標(biāo)βsm得到的二階矩可靠度；SS法為基于傳統(tǒng)串聯(lián)系統(tǒng)假設(shè)，分別以四階矩可靠度指標(biāo)法單獨計算三種失效模式可靠度，并根據(jù)公式(7)得到的可靠度；FM法為本文采用的四階矩可靠度指標(biāo)法，根據(jù)公式(32)(33)得到的可靠度。

從表3可以看出，相對于Edgeworth級數(shù)法和二階矩法，本文所采用的四階矩可靠度指標(biāo)法具有更高的計算精度。在考慮3個基本隨機(jī)變量的情況下，本文方法只需要調(diào)用21次功能函數(shù)就能得到與Monte-Carlo法調(diào)用106次精度相接近的結(jié)果。此外，串聯(lián)系統(tǒng)假設(shè)得到的多失效模式可靠度明顯誤差更大，且比本文考慮共因失效的可靠度保守。

圖3示出了本文考慮共因失效的多失效模式可靠性分析方法(FM)和串聯(lián)系統(tǒng)假設(shè)(SS)的可靠度的對比，可以看出，串聯(lián)系統(tǒng)假設(shè)只有在結(jié)構(gòu)整體可靠度較高時才能適用，結(jié)構(gòu)可靠度越小，用串聯(lián)系統(tǒng)假設(shè)得到的可靠度誤差越大。而本文方法在不同分布參數(shù)下均與Monte-Carlo法的仿真結(jié)果相吻合。

為了研究非正態(tài)分布參數(shù)對二階矩和四階矩方法精度的影響，圖4示出了只考慮失效模式2(ASN2強(qiáng)度失效)的可靠度與應(yīng)力強(qiáng)度Sm形狀參數(shù)間的關(guān)系曲線，圖中虛線為Sm的偏度，實線為不同方法得到的可靠度(其中_SM表示二階矩可靠度、_FM表示四階矩可靠度、_SD表示攝動法求解功能函數(shù)中心矩、_PE表示點估計法求解功能函數(shù)中心矩)。

圖3 不同Sm均值下FM和SS方法可靠度結(jié)果對比

圖4 失效模式2可靠度與Sm形狀參數(shù)間的關(guān)系曲線

由圖4可以看出，對于服從Weibull分布的應(yīng)力強(qiáng)度Sm，當(dāng)形狀參數(shù)逐漸增大時，其偏度也逐漸增大(負(fù)值表示分布向左偏)；當(dāng)Sm的偏度較小時，分別用攝動法和點估計法得到的二階矩可靠度及四階矩可靠度都與Monte-Carlo結(jié)果有相近的精度；當(dāng)偏度達(dá)到一定值時，由于忽略了功能函數(shù)的偏度，二階矩法逐漸偏離Monte-Carlo結(jié)果，此外，基于攝動法的四階矩可靠度與Monte-Carlo法的可靠度誤差也明顯變大，甚至大于二階矩可靠度的誤差。

5 結(jié)語

本文考慮CRDM承壓殼體共因失效條件下的多失效模式可靠性問題，基于有限元和響應(yīng)面法建立了非連續(xù)的多失效模式功能函數(shù)，在考慮任意分布參數(shù)的情況下，用點估計法計算功能函數(shù)前四階中心矩，并依此計算四階矩可靠度指標(biāo)和可靠度。

分析結(jié)果表明，該方法能夠有效求解復(fù)雜結(jié)構(gòu)多失效模式的可靠性問題，相比于Edgeworth級數(shù)四階矩法和傳統(tǒng)二階矩法都有更高的精度，在相同精度的條件下也比Monte-Carlo方法有更小的求解規(guī)模；在求解多失效模式結(jié)構(gòu)可靠性問題時，本文方法也比傳統(tǒng)基于獨立失效假設(shè)的分析方法有更高的精度，適用于核電主設(shè)備的結(jié)構(gòu)可靠性分析。