嚴重事故下核電廠安全殼結構概率性能評價

金 松，李鑫波，貢金鑫

(1.大連理工大學建設工程學部，遼寧，大連 116024；2.海岸和近海工程國家重點實驗室，大連理工大學，遼寧，大連 116024；3.河北工業大學土木與交通學院，天津 300401)

安全殼結構作為核電廠最重要的結構，是事故下的防泄漏屏障，在保證核安全方面起到關鍵作用[1?3]。美國三里島核事故以后，安全殼結構的重要性得到重新認識，研究安全殼結構在超設計基準事故下的行為和相關性能成為熱點。目前，關于安全殼結構在事故壓力下的力學性能開展了一系列的試驗和數值模擬研究[4]。安全殼結構試驗研究主要包括桑迪亞實驗室的1∶4模型試驗[5]，印度BARCOM1∶4模型試驗[6]、法國MAEVA模型試驗[7]、英國Sizewell-B1∶10安全殼縮尺模型試驗[8]，以及國內張心斌等[9]進行的1∶10預應力混凝土安全殼模型試驗。通過模型試驗可以深入了解安全殼結構內壓承載力、破壞過程及破壞形態。但由于安全殼模型試驗成本高，制作難度大，試驗周期長，很難大量重復進行。陽濤等[10]開展了采用碳纖維布加固的安全殼模型試驗研究，試驗發現采用CFRP加固能夠顯著提高安全殼結構的抗壓承載力，并且有效控制安全殼的變形和裂縫發展。復雜非線性有限元分析技術的飛速發展使得模擬安全殼結構成為可能。Huang等[11]在考慮時變因素影響的前提下，分析了安全殼結構在超設計基準工況下的力學性能。Shokoohfar和Rahai[12]采用有限元軟件ABAQUS對預應力混凝土安全殼在壓力和溫度共同作用下的力學性能進行分析。臺灣學者Hu等[13]采用有限元軟件ABAQUS對安全殼結構的極限承載力進行分析，詳細考慮了相關材料的非線性以及預應力損失的影響。Sara等[14]采用數值模擬方法詳細評估了預應力損失、混凝土老化、鋼筋銹蝕、鋼襯里銹蝕四種不同老化情況下安全殼極限承載力，研究發現鋼筋銹蝕對安全殼極限承載力影響最大。Wang[15]對安全殼結構穹頂和筒體交接部位在溫度和壓力共同作用下的力學性能進行理論分析，并將理論分析結果與有限元計算結果進行了比較。Luu等[16]開發了循環軟化膜單元，并且利用循環軟化膜單元對安全殼結構的滯回性能進行分析。

安全殼結構有多處幾何不連續區域(例如設備閘門、人員閘門、截椎體、環梁)，這些幾何不連續區域的存在造成安全殼結構受力機理復雜。其次，安全殼結構在運營期間承受各種不同類型的荷載作用，因此分析安全殼結構力學性能，會涉及各種相關不確定性。此外，概率安全性能評價(PSA)方法是評估核電廠結構安全性能評估的主流分析方法[17]。因此采用概率方法評估安全殼結構性能更加合適。在分析與設計安全殼結構的過程中需要考慮各種內部災害和外部災害工況[18],諸如主蒸汽管道破裂等事故工況會造成安全殼結構內部產生超壓[19]。然而，目前關于安全殼結構超壓下失效概率研究主要集中在易損性研究方面。例如Hoseyni等[20]采用中值法和傳統Monte Carlo方法評價了安全殼結構的超壓易損性。Kim等[21]采用簡化的安全殼結構全截面開裂功能函數，對安全殼結構不同位置的可靠度進行分析。Prinja等[22]采用傳統一次二階矩中心點法和驗算點法對安全殼結構在內壓荷載作用下的可靠度進行分析，該分析同樣也是基于簡化的顯示功能函數。

上述關于安全殼結構在嚴重事故工況下可靠性分析或易損性分析，主要基于簡化的顯示功能函數，或者在易損性分析中考慮的不確定性不夠全面。本文基于安全殼結構三維實體精細化有限元模型，對安全殼結構在嚴重事故內壓的易損性進行分析，開發了相應自動化前后處理程序。然后采用統計推斷法和bootstrap法分析了安全殼結構易損性參數的置信區間。在此基礎上采用bootstrap法計算了安全殼在嚴重事故內壓下的可靠指標和總失效概率的統計特性，最后采用中值法和置信法評估了安全殼結構的安全裕度。本文的相關研究成果可以為安全殼結構嚴重事故工況下概率安全評價(PSA)提供技術參考和研究指導。

1 安全殼結構有限元模型

1.1 安全殼幾何模型

本文研究的是采用改進型百萬千瓦級壓水堆核電技術的某核電廠預應力混凝土安全殼結構。該安全殼結構由底板、筒體、穹頂、環梁、鋼襯里和預應力系統組成。除底板外，安全殼結構的總高度為61.7 m。安全殼結構筒體內半徑為18.5 m，外半徑為19.4 m，壁厚為0.9 m。穹頂的曲率半徑為30.2 m，厚度為0.8 m。在安全殼結構外壁上設置有四個扶壁柱。為保證結構完整性，在安全殼結構內設置了由水平預應力鋼束、純豎向鋼束和Gamma鋼束組成的預應力系統。此外，為保證安全殼結構密封功能，在安全殼結構的內壁上設置了6 mm厚的鋼襯里。安全殼結構筒體上也設置了一些孔洞，用來滿足工藝需要(主要包括設備閘門，人員閘門和應急閘門)。安全殼結構的幾何簡圖如圖1所示。

圖1 安全殼結構簡圖Fig.1 Sketch of containment structure

1.2 材料特性及本構關系

有限元軟件ABAQUS提供了多種材料本構關系模型。其中混凝土材料本構關系包含脆性開裂模型、彌散開裂模型和塑性損傷模型。文獻[23]指出塑性損傷模型最適合模擬混凝土材料的非線性行為。因此，本研究的混凝土本構關系選用經典的塑性損傷模型。混凝土單軸受壓本構關系采用Saenz[24]提出的模型，其表達形式如下：

混凝土受拉本構關系采用文獻[26]建議的本構關系。混凝土塑性損傷模型需要定義損傷因子來描述非線性變形產生的損傷。混凝土受拉和受壓損傷因子的計算方法和相關參數選取參考文獻[27]。普通鋼筋、鋼襯里、預應力筋本構關系采用文獻[28]中建議的理想彈塑性模型。此外，普通鋼筋、鋼襯里的彈性模量和泊松比分別取200 000 MPa和0.3。預應力筋的彈性模量取為195 000 MPa，泊松比同樣取0.3。混凝土受壓本構關系如圖2(a)所示，普通鋼筋、鋼襯里、預應力筋本構關系如圖2(b)所示。

圖2 安全殼結構材料本構關系Fig.2 Constitutiverelationship of materials used in containment structure

1.3 接觸約束

在有限元分析中鋼襯里與混凝土采用共節點方式建模(蒙皮法)，不考慮兩者的相對滑移[29]。普通鋼筋、預應力筋與混凝土采用分離式建模，普通鋼筋和預應力筋都嵌入到混凝土內，不考慮它們與混凝土的相對滑移[30]。此外，Chakraborty等[31]研究發現底板對安全殼結構整體力學行為無明顯影響。因此在本文的有限元模擬中忽略底板，采用約束安全殼底部所有自由度的方法模擬底板對安全殼的約束作用。

1.4 分析方法選取與網格劃分

安全殼結構在發生事故的情況下，其內部會產生高溫和高壓，本文分析中考慮的溫度和壓力時程曲線參考文獻[32](溫度壓力時程曲線如圖3所示)。從圖3可以看出，安全殼的溫度約100 ℃～200℃，最高溫度不超過200℃。文獻[33? 34]研究成果表明：對應在溫度不超過200℃范圍內，溫度作用對安全殼極限內壓承載力影響很小。因此在本文有限元分析中忽略溫度作用。安全殼結構非線性有限元分析可以分為兩個分析步：第一個分析步，施加重力荷載和預應力；第二個分析步，沿安全殼結構內表面施加線性增長的內壓荷載。

圖3 溫度和壓力時程曲線Fig.3 Pressure and temperaturetime history curve

本文采用適中的網格尺寸(0.8 m)，對一些洞口等不連續區域的網格做了加密處理。混凝土采用以8節點六面體線性減縮積分實體單元(C3D8R單元)為主，洞口等不連續區域模擬采用6節點六面體線性減縮積分實體單元(C3D6R單元)模擬。普通鋼筋采用三維4節點面單元(SFM3D4單元)，洞口等不連續區域采用三維3節點面單元(SFM3D3單元)模擬，這種單元只能傳遞平面內薄膜力，無彎曲和剪切剛度。預應力筋采用2節點三維桁架單元(T3D2單元)模擬。劃分完網格的安全殼結構部件如圖4所示。

圖4 安全殼結構網格劃分Fig.4 Finite element mesh of containment structure

1.5 安全殼結構失效準則

安全殼結構在內壓荷載作用下通常有兩種典型的失效模式：即功能性失效和結構性失效。功能性失效表示安全殼結構的功能喪失，通常表現為泄漏率超過允許的最大值，結構性失效表示主要承受內壓材料(通常指的是普通鋼筋和預應力筋)發生失效導致安全殼結構快速泄壓[35]。準確描述安全殼的失效準則十分困難。工程上，通常采用基于材料應變的簡化失效準則。本文采用美國核管會導則RG1.216建議的安全殼結構失效準則[36]：

1)遠離不連續區域的預應力筋的應變達到0.8%(包含加壓前的應變，通常加壓前的應變約為0.4%)。

2)自由場位置處的鋼襯里和普通鋼筋的應變達到0.4%。

1.6 安全殼概率有限元分析

安全殼在設計，施工和交付使用過程存在各種不確定性。這些不確定性主要包含幾何不確定性、模型不確定性、材料不確定性以及荷載不確定性。文獻[37]認為安全殼結構施工質量控制嚴格，結構幾何偏差較小，因此可以忽略安全殼結構的幾何不確定性的影響。本文重點關注材料不確定性、模型不確定性、以及荷載作用不確定性。首先分析材料不確定性對安全殼結構內壓承載力的影響。安全殼結構各種組成材料的概率統計特性參考文獻[26]。采用概率有限元方法對安全殼這種復雜結構進行分析，計算工作量巨大。此外，ABAQUS本身是一個確定性有限元計算平臺，無法自動完成抽樣及前后處理[38]。因此，十分有必要開發自動化抽樣和前后處理程序。本文采用Python平臺和Matlab軟件相結合的方式開發了自動化抽樣和前后處理程序(自動化抽樣和前后處理計算程序開發原理如圖5(a)所示)。安全殼結構易損性自動化分析流程如圖5(b)所示。

圖5 自動化計算程序開發原理及易損性分析流程Fig.5 Principle for developing automatic calculation program and fragility analysis process

2 有限元計算結果

2.1 安全殼內壓承載力計算結果

根據上述破壞準則可以得到對應安全殼樣本的內壓承載力。圖6給出了100個安全殼結構樣本對應的內壓承載力分布的箱線圖。

圖6 安全殼內壓承載力分布Fig.6 Pressurecapacity distribution of containment structure

2.2 易損性曲線計算

安全殼結構的易損性曲線通常采用對數正態分布描述[39]，具體表達形式如下：

利用矩估計法可以得到對應安全殼結構易損性參數如下：

式中，pi為安全殼樣本i的內壓承載力。

根據式(3)和式(4)可以得到安全殼結構對應易損性曲線如圖7所示。

圖7 安全殼易損性曲線Fig.7 Fragility curve of containment structure

3 安全殼可靠性及概率性能評價

3.1 易損性參數統計不確定性及其置信區間

結構易損性分析通常采用有限的樣本進行有限元分析并擬合對應易損性曲線。采用有限數量的樣本數量分析得到的結果就會涉及統計不確定性[40]。目前有多種方法可以計算由于統計不確定性造成結構易損性參數的變化，主要包含統計推斷法和bootstrap法。本文采用統計推斷法和bootstrap兩種方法計算易損性參數的置信區間。

a)統計推斷法

b)bootstrap法

bootstrap是一種統計推斷方法，按照抽樣方式的不同可以分為參數化bootstrap和非參數化bootstrap。本文采用Matlab軟件分別編制了參數化bootstrap和非參數化bootstrap程序。利用程序可以計算得到100次、500次、1000次重抽樣對應安全殼結構的內壓易損性曲線族如圖8和圖9所示。

圖8 不同非參數bootstrap重抽樣次數對應安全殼易損性曲線Fig.8 Fragility curvesof nuclear containment structure with different number of nonp-arametric bootstrap re-samplings

圖9 不同參數bootstrap重抽樣次數對應安全殼易損性曲線Fig.9 Fragility curves of containment structure with different number of parametric bootstrap re-samplings

表1列出了采用兩種bootstrap法(參數化bootstrap法和非參數化bootstrap法)對應100次、500次和1000次重抽樣計算的安全殼結構內壓易損性參數pm和 βS的統計特性。從表1可以看出，易損性參數pm和 βS對重抽樣次數和重抽樣方法不敏感。因此，本文采用100次參數化bootstrap重抽樣結果作為易損性參數置信區間構造的依據。圖10給出了100次參數bootstrap重抽樣對應90%置信水平易損性曲線的置信區間。

圖10 100次參數bootstrap重抽樣及其對應90%置信水平易損性曲線Fig.10 Fragility curves of 100 parametric bootstrap resamplings with 90%confidence level

100次參數化bootstrap重抽樣得到易損性參數pm對應90%置信水平的置信區間為[1.596 MPa,1.603 MPa]，對應易損性參數βS的置信區間為[0.0191,0.0195]。對比統計推斷法和bootstrap法預測的置信區間可以看出,兩種方法對易損性參數pm置信區間的估計幾乎相同。對于易損性參數βS，兩種方法估計的置信區間差別較大，統計推斷法往往會高估易損性參數 βS的置信區間。需要說明的是，可以認為100次參數化bootstrap方法估計的易損性參數置信區間為準確結果。

表1不同bootstrap次數計算易損性參數統計特性Table 1 Statistical characteristics of fragility parameters with different number of parametric bootstrap re-samplings

3.2 安全殼可靠指標與總失效概率分析

為了更加合理的評估安全殼結構的內壓承載力，需要考慮模型不確定性。不同堆型的安全殼結構的模型不確定性取值大小差異[41?42]。文獻[20]指出采用先張法和后張法施工的預應力混凝土安全殼模型不確定性 βM可取0.12。此外安全殼結構的事故內壓服從均值μpa=0.663 MPa，對數標準差βpa=0.3的對數正態分布。安全殼結構的可靠指標β和總失效概率CCFP可以采用下式進行計算[26]：

根據上述100次參數化bootstrap重抽樣計算結果，可以采用數理統計法得到對應可靠指標的均值μβ和變異系數 δβ。具體計算公式如下：

其中：

總失效概率CCFP 的統計特性可以采用同樣方法確定。安全殼結構的可靠指標β和總失效概率CCFP的統計特性匯總于表2。從表2可以看出，可靠指標 β和總失效概率CCFP的變異性都很小，且總失效概率CCFP的變異性明顯高于可靠指標β。

表2 安全殼可靠指標和總失效概率統計特性Table 2 Statistical characteristicsof reliability index and cumulative failure probability for containment structure

3.3 安全殼結構概率安全裕度及概率安全性能評價

3.3.1安全殼結構概率安全裕度評價

安全殼結構安全裕度評價是安全殼結構概率安全評價(PSA)工作中的一項重要任務。同時，安全殼結構的安全裕度是核電廠安全審查中關注的重要指標。本文采用中值法和置信法這兩種方法評價安全殼結構的安全裕度。中值法的具體計算思路如下，根據文獻[27]中給出不確定性SRSS組合方法，可以得到安全殼結構總的不確定性參數如下：

安全殼對應95%保證率對應的平均安全裕度如下[26]：

此外，也可采用置信法來評價安全殼結構的安全裕度。具體計算思路如下：置信水平為Q的安全殼結構易損性曲線可以采用如下的形式表達：

對比式(16)和式(18)可知，只須將式(16)中的pm替換成pm·exp(?βM·Φ?1(Q))即可計算得到具有置信水平的安全裕度。定義：

于是可得置信水平為Q的安全裕度表達式如下：

采用中值法和置信法計算得到安全殼結構的安全裕度列在表3。從表3可以看出，采用中值法計算的安全裕度與95%置信水平的置信安全裕度相差不大。置信安全裕度隨著置信水平的提高不斷減小。總體來看，本文研究的安全殼結構可以滿足安全裕度不小于2.5的要求。

表3 安全殼結構安全裕度Table 3 Safety margin of containment structure

3.3.2安全殼概率安全性能評價

安全殼結構在嚴重事故下的性能主要包含確定性能目標和概率性能目標兩個方面。確定性性能目標指的是在堆芯毀損24 h之內，安全殼結構應保證其結構完整性。概率性能指標指的是安全殼結構在嚴重事故下的總失效概率不超過10%[43]。本文采用概率性能目標作為安全殼結構性能評價的依據，從表2的計算結果可以看出，安全殼結構總失效概率基本上穩定在0.0033，遠小于10%，說明本文分析的安全殼結構滿足嚴重事故下的性能要求。

4 結論

本文基于安全殼三維精細化實體有限元模型，對安全殼結構嚴重事故工況下的易損性和概率安全性能進行詳細分析。通過上述研究主要得到以下結論：

(1)統計推斷法和bootstrap法對易損性參數pm置信區間估計幾乎相同，就易損性參數βS而言，兩種方法估計的置信區間差別較大，統計推斷法往往會高估易損性參數 βS的置信區間。

(2)考慮各種不確定性以后，安全殼結構在內壓作用下的總失效概率基本穩定在0.0033。說明本文分析的安全殼結構能滿足嚴重事故下的性能要求。此外，安全殼結構內壓易損性參數的統計不確定性對可靠指標 β和總失效概率CCFP的影響不大，并且易損性參數的統計不確定性對總失效概率的影響明顯大于對可靠指標的影響。

(3)中值法計算的安全裕度與置信法計算的置信水平為95%的置信安全裕度相差不大。同時，置信安全裕度隨著置信水平的增大不斷減小。總體來看，安全殼結構可以滿足安全裕度不小于2.5的要求。