基于貝葉斯估計的核電廠安全殼內壓概率安全評估

田澳楠,鄭 志,潘曉蘭,蘇春陽,王 勇

(太原理工大學 土木工程學院,山西 太原 030024)

近些年來,頻發的地震災害給核電廠的安全性能帶來了嚴重威脅,這些超出設計基準的地震事故可能會對核電廠的主要系統和部件帶來不可預測的損壞。為此,許多研究人員對核電廠的安全性能進行了概率評估,以抵御核電廠在嚴重事故下的潛在影響。Prinja等[1]使用一階可靠度分析法預測了安全殼在超壓荷載下的失效概率,并分析了帶鋼襯里的安全殼和不帶鋼襯里的安全殼的結構可靠性。Kim等[2]對考慮長期鋼筋束退化的安全殼結構進行了內壓荷載下的可靠性分析,結果表明鋼筋束數量與安全殼的可靠性呈正相關關系。Hoseyni等[3]采用中值壓力法對安全殼在超壓作用下的失效概率進行了綜合評估,并給出了安全殼在特定壓力下的失效概率和極限承壓能力。Granger等[4]采用置信度的方法對安全殼進行了在內壓荷載下的概率安全評估,并對安全殼的抗泄漏能力進行了分析。金松等[5]考慮材料不確定性、統計不確定性以及建模不確定性等因素并對安全殼開展了概率安全評估,最后使用泰勒展開方法得到了安全殼的總失效概率。Zhao等[6]在考慮結構不確定性的情況下提出了一種基于PE-IDA的地震易損性分析方法,并通過單自由度體系對該方法進行了驗證。結果表明,結構參數的不確定性對易損性的影響與結構損傷程度呈正比關系。Kwag等[7]提出了一種基于完全采樣的地震概率評估量化方法,該方法能夠準確地表示安全殼結構地震易損性信息的相關性。王曉磊等[8-9]對我國核電廠安全殼進行了地震易損性研究,并基于地震風險解析函數和風險卷積函數評估了我國某核電廠安全殼地震風險,提出了適用于我國核電廠地震易損性及風險評估的方法。上述學者在對安全殼進行事故荷載或地震作用下的易損性分析時均假定易損性函數服從對數正態分布模型,然而該模型未能綜合考慮認知、經驗、統計等多種不確定性因素的影響,這可能會極大影響安全殼概率安全評估的準確性。當前,Gardoni等[10-11]提出了完整的易損性評估理論,其通過采用貝葉斯估計建立結構的概率能力模型和需求模型,并綜合考慮模型所有相關的不確定性,最大限度地提高了易損性評估和概率安全評估的準確性,然而該理論還未應用于核電廠安全殼的易損性和概率安全評估。應用貝葉斯估計理論進行安全殼的易損性和概率安全評估面臨兩大難點:1) 建立安全殼內壓概率需求模型時首先需要確定能夠反映結構整體反應變化的參數;2) 準確的概率需求模型還需要兼顧模型的簡潔性,這需要對模型進行修正。

為此,本文提出一種新的混凝土損傷參數來評估安全殼混凝土的整體損傷性能,同時利用該參數建立確定性的內壓需求模型。然后采用貝葉斯估計方法建立安全殼的概率需求模型,并對模型參數進行優化。在此基礎上對不同損傷狀態下的安全殼開展易損性評估及概率安全評估。

1 安全殼結構有限元模型

1.1 安全殼幾何模型介紹

預應力混凝土安全殼由穹頂、筒體、鋼襯里、扶壁柱、普通鋼筋、預應力筋等組成。安全殼結構總高度為69.0 m,其中筒體高度為48.0 m,穹頂高度為21.0 m。穹頂可以減小安全殼在極高內壓下的不利影響,穹頂的跨度和曲率半徑分別為40.0 m和20.0 m;筒體的內外徑分別為40.0 m和42.2 m。雙層普通鋼筋分別為豎直和水平分布。為了保證安全殼結構的完整性,在安全殼筒壁上布置了預應力筋系統。預應力筋系統由330根預應力筋組成,其中190根水平預應力筋錨固在扶壁柱上,140根倒“U”型預應力筋錨固在底板上。為了保證安全殼的密封性,在安全殼的內表面設置了6 mm厚的鋼襯里。此外,為了滿足安全殼的功能要求,在筒壁上留有一個直徑為7.0 m的設備洞口。安全殼的幾何簡圖如圖1所示。

圖1 安全殼幾何簡圖

1.2 安全殼材料屬性及其本構關系

混凝土的強度等級為C50。一般來說,混凝土被認為是各向異性的材料,其本構關系是影響裂縫演化的關鍵因素。為了盡可能準確地模擬安全殼的非線性行為,采用了經典的混凝土塑性損傷(CDP)模型[12]。該模型采用各向同性損傷和各向同性拉壓彈性的概念來模擬混凝土的非線性行為,可以很好地表示混凝土在受力情況下的損傷演變。CDP模型通過式(1)定義拉伸破壞和壓縮破壞;非彈性應變以式(2)表示;等效塑性應變可用式(3)來實現。損傷參數的確定采用了Sidoroff破壞模型[13]。

(1)

(2)

(3)

普通鋼筋、鋼襯里及預應力筋采用理想彈塑性模型。泊松比均取0.3,鋼襯里和預應力筋的彈性模量為2×105MPa,鋼筋的彈性模量為1.95×105MPa。鋼襯里和鋼筋的屈服強度分別為320 MPa、400 MPa,預應力筋的極限抗拉強度為1 860 MPa。

1.3 安全殼有限元分析方法

本研究通過3個加載步驟來獲得安全殼在內壓作用下的易損性評估結果。第1步施加重力荷載;第2步通過降溫法[14]給預應力筋施加預應力,如式(4)所示;第3步施加線性增加的內壓。

(4)

式中:ΔT為預應力筋的降溫溫度;Δσ為控制應力;λ為線膨脹系數;E為預應力筋的彈性模量。

1.4 有限元網格劃分及邊界條件設定

安全殼有限元模型計算時混凝土容易產生局部應力集中現象,網格單元過密則會產生頻繁的不收斂問題。因此在劃分網格時,在保持精度的前提下選取適當的網格尺寸。表1列出不同網格尺寸下完成安全殼預應力筋張拉時的相關統計信息。

表1 不同網格尺寸的統計信息

由表1可看出,采用0.8 m的網格尺寸可以在保證精度的前提下節約模型運行時間。為精確評估安全殼洞口的破壞表現,在洞口區域對網格進行加密處理。在劃分網格單元時,混凝土采用實體單元(C3D8R)、鋼襯里采用殼單元(S4R)模擬,兩者通過共節點連接,且不考慮它們之間的相對滑移[15]。為減小鋼筋建模的復雜性,普通鋼筋采用表面單元(SFM3D4)模擬,預應力筋通過桁架單元(T3D2)模擬。安全殼的結構網格劃分如圖2所示。鋼筋、預應力筋與混凝土采用分離式建模,其中鋼筋和預應力筋完全嵌入混凝土中,不考慮它們與混凝土之間的相對滑移[16]。

a——混凝土;b——普通鋼筋;c——鋼襯里;d——預應力筋

1.5 安全殼損傷參數定義

安全殼結構的完整性和密閉性一直是眾多學者關注的重點。當前,工程界常采用應力、應變或位移來定義安全殼的功能或結構失效,這類定義能夠有效反映結構的局部開裂或損傷情況,卻不能直接反映結構整體的開裂或損傷情況。在持續增長的內壓作用下,安全殼各部件的強度、剛度會發生退化和損傷,結構損傷不斷累積,最終產生破壞。綜合來說,安全殼混凝土的損傷行為是一個不斷累積的過程,在承壓過程中會產生貫穿裂縫并不斷擴展至整個筒壁,已經產生貫穿裂縫的區域會發生強度和剛度的退化,這會對安全殼整體損傷和失效產生影響?；诖?本研究提出了“損傷面積比”的參數來量化安全殼混凝土的損傷行為,以完整地反映安全殼在損傷演化過程中的失效情況。

(5)

式中:DRc為混凝土損傷面積比;AreaD為混凝土產生貫穿裂縫的區域面積;AreaT為混凝土的總面積。DRc的取值范圍為[0,1],當DRc=0時代表安全殼混凝土尚未產生貫穿裂縫,當DRc=1時代表安全殼混凝土均已產生貫穿裂縫。

當混凝土拉應變超過峰值拉應變時,混凝土損傷超過損傷限值,即代表安全殼混凝土層發生受拉破壞,關于計算混凝土產生貫穿裂縫的區域面積在ABAQUS程序的后處理模塊完成,詳細計算步驟如下:1) 依據《混凝土設計規范》(GB50010—2010),C50混凝土的標準抗拉強度為2.64 MPa,彈性模量為34 500 MPa,由此可計算得到C50混凝土的峰值拉應變和相應的損傷因子;2) 在ABAQUS中的Contour Plot模塊中定義該損傷因子為損傷限值,當混凝土的損傷超過損傷限值時即代表安全殼混凝土層發生受拉破壞;3) 在ABAQUS中選定混凝土損傷超過損傷限值的區域,定義為AreaD。

圖3示出安全殼混凝土產生貫穿裂縫的判別方式,左側區域混凝土損傷穿過整個混凝土層,其為混凝土貫穿裂紋區域;右側區域混凝土損傷未穿過整個混凝土層,因此不能判別為混凝土貫穿裂紋區域。

圖3 安全殼混凝土裂縫貫穿

1.6 安全殼概率有限元分析

開展安全殼概率安全評估時,有必要綜合考量多種誤差因素對結果造成的影響。這些誤差因素主要包括安全殼幾何不確定性、模型不確定性、材料不確定性以及荷載不確定性等。文獻[17]認為安全殼是施工質量受到嚴格控制的結構,可以忽略由幾何不確定性帶來的影響。因此,本文在進行安全殼概率安全評估中重點考察了材料不確定性、統計不確定性、荷載不確定性以及模型不確定性的影響。安全殼各種材料的統計特性[15,18]列于表2。表2中:fc為混凝土單軸抗壓強度;fy400為鋼筋屈服強度;Ey400為鋼筋彈性模量;fl為鋼襯里屈服強度;El為鋼襯里彈性模量;fp為預應力筋屈服強度;Ep為預應力筋彈性模量。利用表2給出的7個材料參數的分布形式和范圍,本文采用拉丁抽樣方法在Matlab中生成了100組樣本,依次導入ABAQUS中運行得到計算結果。

表2 隨機變量的統計特性

2 易損性分析理論

傳統的易損性分析通常假定易損性函數服從對數正態分布的形式,如式(6)所示。這種假定使得進行易損性評估時僅需要估計關鍵的易損性參數即可,大大簡化了易損性評估的計算量。然而,這種簡化使得易損性評估缺乏嚴格的理論分析,并且不能全面考慮概率分析中存在的認知、經驗、統計等相關不確定性因素,這可能會對易損性評估結果產生較大影響。

(6)

式中:F為失效概率;Φ為標準正態累積分布函數;pm為安全殼結構承壓能力中值;βs為對數標準差;p為安全殼事故壓力。

貝葉斯估計是一種廣泛應用于結構易損性評估和可靠性分析的統計學方法。它通過建立結構的概率能力和概率需求模型綜合考慮概率分析中的相關不確定性因素,從而得到更可靠的易損性分析結果。

不論安全殼材料、承受的荷載如何變化,安全殼混凝土的損傷面積比情況反映了混凝土整體的損傷情況以及可能的泄漏水平,因此本研究中假定安全殼概率能力模型是確定的損傷面積比。而安全殼概率需求模型可通過下式[11]確定:

(7)

(8)

由于安全殼在內壓作用下的力學行為較為復雜,很難建立損傷與內壓的確切物理力學關系。因此,確定性模型通過擬合100組樣本的損傷中值與內壓的數學關系得到,如圖4所示?？梢园l現,安全殼混凝土的損傷中值與內壓近似服從均值為1.351 8、方差為0.112 8的累積正態分布,如式(9)所示。

(9)

圖4 混凝土損傷面積比與內壓的數學關系

式中,P為內壓。

在構建安全殼概率需求模型時,混凝土抗壓強度(h1(x)=fc)、鋼筋強度(h2(x)=fy400)、鋼筋彈性模量(h3(x)=Ey400)、預應力筋強度(h4(x)=fp)、預應力筋彈性模量(h5(x)=Ep)、鋼襯里強度(h6(x)=fl)、鋼襯里彈性模量(h7(x)=El)及內壓(h8(x)=P)被選擇作為確定性需求模型的修正項。然而,修正因子過多可能會導致模型過擬合,模型變得過于復雜也會導致其無法很好地推廣到新樣本,因此有必要對修正項進行參數優化和校正。

模型優化和校正過程分為以下3步:1) 采用貝葉斯估計方法計算所有模型參數和標準差的后驗分布;2) 對比解釋函數的變異系數,去除變異系數最大的一項,該項被認為是解釋函數中包含信息量最少的一項;3) 去除某項后如果后驗均值未發生明顯變化,則認為去除該項是合理的,重復此流程,直到后驗均值明顯增大,認為去除該項影響了模型的精度,則該項不能去除,因此,返回上一步,模型優化停止。

按照上述的參數優化方法對參數進行取舍,如表3所列。第1次去除了參數θ3(1.867),第2次去除了參數θ7(21.419),隨后依次去除了θ5(2.044)、θ2(5.864)、θ6(1.848),當去除參數θ1(0.646)后,模型標準差顯著增大,這表明去掉該項會使模型精度降低,因此,參數優化過程停止,取上一次優化的結果作為最優模型。

表3 模型參數的統計信息

經貝葉斯估計方法優化后,參數θ1、θ4、θ8最終保留為修正項,因此,安全殼混凝土的損傷-內壓概率需求模型可簡化為式(10)。計算模型每個參數的標準差和變異系數,以及后驗分布的均值,如表4所列。

(10)

表4 模型參數的后驗統計

極限狀態方程定義了安全殼的安全狀態與失效狀態的邊界,這在進行易損性評估中起著關鍵作用。通常,極限狀態方程的形式表達為能力與需求之差,形式如下:

Zk=Ck(x)-Dk(x)

(11)

式中:Zk為安全殼混凝土的第k個失效狀態的安全余量;Ck為第k個狀態下的能力函數;Dk為第k個狀態下的需求函數。

需求模型是包含隨機變量的函數,因此具有概率分布,而能力模型被認為是安全殼特定狀態下的恒定損傷值,它反映了安全殼的整體損傷與泄漏水平,因此極限狀態方程也具有概率分布。當Zk<0時表明該狀態下安全殼的損傷能力值小于損傷需求值,即安全殼達到該狀態下的失效水平。安全殼的易損性評估可采用極限狀態方程的積分形式來表達:

(12)

式中:Fk為第k個狀態下的安全殼的失效概率;P為第k個狀態下安全殼損傷能力值小于損傷需求值的概率;fzk(Z)為概率密度函數;Z為功能函數。

概率需求模型由多個具有不確定性因素的隨機變量組成,這導致易損性評估存在不確定性。因此,有必要給出具有置信邊界的易損性曲線,該曲線可以表征結構在特定損傷狀態下的失效概率范圍,獲得更加精確的易損性評估結果。易損性曲線的置信邊界可通過1階分析法獲得[10],已知失效概率和可靠性指標存在如下關系:

β(c,d,Θ)=Φ-1[1-Fk(c,d,Θ)]

(13)

式中:β為可靠性指標;c為能力模型中的變量;d為需求模型中的變量;Θ為模型參數。

可靠性指標β的方差可近似由下式獲得:

(14)

置信邊界可被定為{Φ[-β-σβ]、Φ[-β+σβ]},該置信邊界對應于15%和85%的置信水平。

3 有限元計算結果

3.1 安全殼樣本差異性分析

圖5示出考慮材料不確定性情況下混凝土損傷演化的差異。根據混凝土的抽樣分布情況將結果樣本分為Mean-1Std、Mean和Mean+1Std 3個區域,Mean表示混凝土抗拉強度均值,Std表示混凝土抗拉強度標準差。從每個區域內隨機挑選出樣本,比較其損傷演化和樣本間的差異性。在0.8 MPa內壓下,混凝土裂縫集中在設備洞口上下兩側,此時,由材料不確定性引起的變異性并不明顯。當內壓增長至1.0 MPa時,選自Mean-1Std區域的安全殼樣本的混凝土裂縫已經擴展至穹頂頂部,而選自Mean及Mean+1Std區域的安全殼樣本尚未發生明顯變化。隨著內壓繼續增長至1.2 MPa,選自Mean-1Std區域的安全殼混凝土裂縫廣泛開展至穹頂、洞口周圍及底部,樣本混凝土損傷已經達到18.61%,而選自Mean及Mean+1Std區域的安全殼樣本損傷僅達到4.63%和1.99%。內壓增長至1.4 MPa時,選自Mean-1Std區域的安全殼混凝土裂縫擴展至整個筒體,損傷程度遠超過選自Mean及Mean+1Std區域的安全殼樣本。當內壓增長至1.6 MPa時,安全殼混凝土絕大部分開裂,樣本差異性則明顯減小。結果表明,由材料不確定性引起的安全殼樣本損傷的差異性不可忽略,尤其當內壓處于[1.2 MPa,1.6 MPa]時,這種差異性更為明顯。采用混凝土損傷面積比的形式可以反映出安全殼混凝土在內壓作用下的損傷演化行為,并為研究其在不同損傷水平下的失效概率奠定了基礎。

圖5 安全殼混凝土損傷演化

3.2 易損性曲線分析結果

圖6示出安全殼在不同損傷狀態下的概率密度曲線。隨著混凝土損傷面積比增大,概率密度曲線的峰值逐漸向右移動,這表明樣本損傷中值逐漸增加。曲線的峰值逐漸減小且曲線的輪廓逐漸變寬,這表明樣本分布的標準差逐漸增大,樣本的可變性逐漸增加,由非彈性行為引起的不確定性增加。

圖6 不同損傷狀態下的概率密度曲線

表5列出相同損傷下采用傳統方式和貝葉斯估計獲得的易損性參數,并繪制了相應的易損性曲線,如圖7所示。可以看出,在相同損傷下,采用貝葉斯估計和傳統方式獲得的易損性參數中值相差不大,采用貝葉斯估計獲得的易損性參數變異系數大于傳統方式,使得貝葉斯易損性曲線相對平緩。這是由于貝葉斯估計在構建易損性曲線時經過了嚴格的理論分析,并且綜合考慮了各項材料參數對需求模型的影響,最終獲得的易損性曲線具有較高的準確度。

表5 傳統方式和貝葉斯估計易損性參數對比

DRc:a——0.1;b——0.2;c——0.3;d——0.4;e——0.5;f——0.6;g——0.7;h——0.8;i——0.9

圖8示出不同損傷狀態下的貝葉斯易損性曲線及相應的置信邊界。隨著混凝土損傷面積比增大,易損性曲線的能力中值逐漸增大,這表明在內壓一定時,安全殼越難以發生損傷程度嚴重的失效。如在內壓增長至1.2 MPa時,安全殼混凝土有49.7%的可能發生損傷程度為0.1的失效,而僅有7.4%的可能發生損傷程度為0.5的失效。通過量化不同損傷狀態下安全殼的失效概率有利于掌握安全殼在內壓增長過程中的損傷及失效情況,這有利于監管者針對不同的損傷情況做出相應的防護措施。

DRc:a——0.1;b——0.2;c——0.3;d——0.4;e——0.5;f——0.6;g——0.7;h——0.8;i——0.9

4 安全殼概率安全評價

概率安全評估被認為是核工程界評估易損性的重要標準。通常情況下,安全殼結構的事故壓力遵循對數正態分布,中值為0.663 MPa,對數標準差為0.3[19]。安全殼結構的總失效概率(CCFP)可采用下式進行計算:

(15)

式中:F(fail|p)為安全殼結構在特定損傷狀態下的易損性曲線;fp(p)為事故壓力的概率密度函數。根據式(15)可計算得到安全殼在任一損傷狀態下的總失效概率,表6列出采用傳統易損性評估方法和貝葉斯估計方法得到的安全殼在不同損傷狀態下的總失效概率的相關信息。表6中:μCCFP為總失效概率的均值;σCCFP為總失效概率的標準差;δCCFP為總失效概率的變異系數。

表6 安全殼結構的總失效概率的統計信息

由表6可見,總體而言安全殼的總失效概率隨著混凝土損傷面積比的增大而減小,且安全殼在任一損傷狀態下的總失效概率均小于0.1[20],這表明本文分析的安全殼結構在損傷演化的全過程中均滿足嚴重事故下的性能要求。不同之處在于,采用貝葉斯估計方法獲得的總失效概率的均值大于傳統易損性評估方法,而變異系數小于傳統易損性評估方法。這表明傳統易損性評估可能低估了安全殼在損傷狀態下的總失效概率,采用貝葉斯估計方法可以取得更加保守的概率安全評估結果。

5 結論

本文采用貝葉斯估計方法建立了安全殼的損傷-內壓概率需求模型,從損傷演化的角度對安全殼在內壓作用下的易損性和概率安全性能進行了全面的評估,得到了以下結論。

1) 由材料不確定性引起的安全殼樣本間的損傷變異性不可忽略,不同的樣本在承受相同內壓時往往表現出不同的損傷情況,這為進行安全殼的易損性評估奠定了基礎。

2) 在綜合考慮經驗、認知、材料、統計等相關不確定因素的基礎上,本文采用貝葉斯估計建立了安全殼內壓-損傷概率需求模型,進而實現了安全殼在不同損傷狀態下的易損性和概率安全評估結果。這有助于從概率角度掌握安全殼增壓全過程中的整體損傷及失效情況,有利于監管者針對不同損傷情況做出相應的防護措施。

3) 安全殼的總失效概率隨著混凝土損傷面積比的增大而逐漸減小。采用貝葉斯估計方法獲得的總失效概率的均值大于傳統易損性評估方法,而變異系數小于傳統易損性評估方法。本文建議采用貝葉斯估計方法,其可以取得更加保守的概率安全評估結果,且評估結果更加穩定。