基于不同抽樣方法的核電廠大破口失水事故BEPU 分析

王洋洋，孫曉暉，石興偉

（1.電力規(guī)劃設(shè)計總院，北京 100120；2.中國核電工程有限公司，北京 100840；3.生態(tài)環(huán)境部核與輻射安全中心，北京 100082）

最佳估算加不確定性（BEPU）分析方法同傳統(tǒng)的保守性分析方法相比，可獲取相對現(xiàn)實的分析結(jié)果和安全裕量，能在保證核電廠安全性的前提下，提高核電廠的經(jīng)濟(jì)性和運(yùn)行靈活性。根據(jù)國際原子能機(jī)構(gòu)（IAEA）的定義[1]，最佳估算加不確定性分析方法應(yīng)滿足3 個條件：對于選定的驗收準(zhǔn)則，事故分析中不有意引入悲觀性；使用最佳估算程序；對計算結(jié)果進(jìn)行不確定性分析。不確定性分析方法依據(jù)不確定性傳輸主體的不同可分為“輸入不確定性傳播”方法和“輸出不確定性傳播”方法[1]。典型的“輸入不確定性傳播”方法有美國的CSAU 和ASTRUM、德國的GRS、法國的IPSN 和西班牙的ENUSA 等。典型的“輸出不確定性傳播”方法為意大利的UMAE/CIAU。目前相比“輸出不確定性傳播”方法，“輸入不確定性傳播”方法被世界更多國家研究和使用?！拜斎氩淮_定性傳播”方法以熱工水力系統(tǒng)程序（最佳估算程序）為不確定性傳輸媒介，將輸入不確定性傳播到輸出不確定性，其具體流程如圖1 所示?；谔囟ǖ某闃臃椒▽Σ淮_定性重要輸入?yún)?shù)進(jìn)行隨機(jī)抽樣是不確定性評估流程的重要環(huán)節(jié)，隨機(jī)抽樣參數(shù)作為熱工水力系統(tǒng)程序的輸入，將直接關(guān)系到目標(biāo)響應(yīng)參數(shù)的不確定性量化結(jié)果。

圖1 基于輸入傳遞的不確定性評估流程

美國核管會（NRC） 在1974 年發(fā)布的10CFR 50.46 輕水堆核電廠應(yīng)急堆芯冷卻系統(tǒng)（ECCS）驗收準(zhǔn)則[2]，成為輕水堆保守性失水事故（LOCA）分析的國際通用規(guī)范。失水事故驗收準(zhǔn)則為：包殼峰值溫度不能超過限值（1 477 K）以防止包殼脆化；包殼總氧化率不超過氧化前包殼總厚度的17%；燃料包殼與水或蒸汽發(fā)生化學(xué)反應(yīng)的產(chǎn)氫量不應(yīng)超過堆芯所有鋯產(chǎn)氫量的1%；堆芯保持可冷卻的幾何形狀；保持對堆芯的長期冷卻。NRC 在1988 年對10CFR 50.46 及其附錄K 進(jìn)行了修訂，LOCA 驗收準(zhǔn)則保持不變，允許采用BEPU分析方法進(jìn)行LOCA 分析，需將不確定性分析結(jié)果同驗收準(zhǔn)則進(jìn)行比較，確保驗收準(zhǔn)則規(guī)定限制有很高的概率不被超過。盡管許多國家采用TRACE、CATHARE和RELAP5 等最佳估算程序?qū)Υ笃瓶谑鹿剩↙BLOCA）進(jìn)行了不確定性分析，但系統(tǒng)地基于不同抽樣方法開展響應(yīng)參數(shù)不確定性量化與敏感性分析工作尚無。

本文采用RELAP5 程序建立百萬千瓦級壓水堆核電廠大破口失水事故模型，并對模型進(jìn)行穩(wěn)態(tài)調(diào)試和瞬態(tài)分析；采用Python3.7 基于簡單隨機(jī)抽樣（SRS）和拉丁超立方抽樣（LHS）對重要不確定性輸入?yún)?shù)進(jìn)行隨機(jī)抽樣產(chǎn)生模型輸入矩陣，通過模型對輸入矩陣實施計算獲取關(guān)鍵安全參數(shù)；采用Wilks 非參數(shù)統(tǒng)計方法對關(guān)鍵安全參數(shù)進(jìn)行不確定性量化以獲取單側(cè)統(tǒng)計容忍上限，采用Spearman 秩相關(guān)系數(shù)法對關(guān)鍵安全參數(shù)進(jìn)行敏感性分析以識別主要影響的輸入?yún)?shù)。并評估2 種抽樣方法在不確定性量化和敏感性分析中的計算結(jié)果的差異。

1 模型建立

1.1 程序介紹

RELAP5/MOD3.4 程序基于非均勻、非平衡、兩流體動力學(xué)模型，采用快速半隱式數(shù)值技術(shù)求解，能夠模擬反應(yīng)堆在幾乎各種運(yùn)行瞬態(tài)或假想事故下的行為，廣泛應(yīng)用于事故分析、規(guī)程制定、審評計算和事故緩解措施評價等各個方面。是目前國際通用的輕水堆熱工水力瞬態(tài)分析程序，可用于反應(yīng)堆破口類事故分析。

1.2 模型建立

某百萬千瓦級壓水堆核電廠采用我國具有完整自主知識產(chǎn)權(quán)的第三代壓水堆核電技術(shù)。核電廠主系統(tǒng)具有3 個環(huán)路，每個環(huán)路有1 臺主泵和1 臺蒸汽發(fā)生器，系統(tǒng)共用1 臺穩(wěn)壓器。反應(yīng)堆熱功率為3 180 MW，冷卻劑平均溫度為310℃，系統(tǒng)正常運(yùn)行壓力為15.5MPa。

采用RELAP5/MOD3.4 程序針對百萬千瓦級壓水堆核電廠的設(shè)計特點，對大破口失水事故建立最佳估算分析模型，系統(tǒng)節(jié)點劃分如圖2 所示。建模對象主要包括反應(yīng)堆壓力容器、穩(wěn)壓器、蒸汽發(fā)生器、主泵、主蒸汽系統(tǒng)、安注系統(tǒng)、主給水系統(tǒng)、輔給水系統(tǒng)、主管道和破口等。其中，壓力容器針對入口、下降環(huán)腔、下封頭、下腔室、堆芯通道、堆芯旁通通道、上腔室和上封頭等進(jìn)行了模擬。安全殼采用時間相關(guān)控制體進(jìn)行模擬，大破口失水事故的發(fā)生通過觸發(fā)閥打開進(jìn)行模擬。

圖2 RELAP5 系統(tǒng)節(jié)點

1.3 穩(wěn)態(tài)調(diào)試

在進(jìn)行大破口失水事故瞬態(tài)計算之前，需要對建立的模型進(jìn)行穩(wěn)態(tài)調(diào)試，以確保各主要參數(shù)同電廠名義值一致。表1 為模型的穩(wěn)態(tài)計算值和電廠名義值的相對偏差，從模型穩(wěn)態(tài)計算結(jié)果可見，主要參數(shù)模型穩(wěn)態(tài)計算值和名義值符合較好。

表1 模型穩(wěn)態(tài)計算值相對偏差%

1.4 瞬態(tài)計算

模型在穩(wěn)態(tài)運(yùn)行1 000 s 后發(fā)生冷管段雙端剪切斷裂大破口失水事故，事故瞬態(tài)計算時間為300 s，大破口失水事故事件序列見表2。圖3—5 分別為大破口失水事故下歸一化堆芯功率、穩(wěn)壓器壓力、破口總流量的模型計算曲線和參考曲線。從瞬態(tài)計算結(jié)果可見，計算曲線與參考曲線大小趨勢基本一致。

表2 大破口失水事故事件序列s

圖3 堆芯功率

圖4 穩(wěn)壓器壓力

圖5 破口流量

2 理論方法

2.1 簡單隨機(jī)抽樣

2.1.1 均勻分布參數(shù)的簡單隨機(jī)抽樣

對于均勻分布參數(shù)，簡單隨機(jī)抽樣實現(xiàn)步驟如下。

1）通過公式（1）采用Python3.7 程序中的random 模塊調(diào)用random（）方法生成0 到1 之間的隨機(jī)浮點數(shù)u。

式中：random 為Python3.7 程序中的模塊，random（）為模塊中的方法。

2）通過公式（2）將隨機(jī)浮點數(shù)u映射為參數(shù)分布范圍內(nèi)的隨機(jī)數(shù)。

式中：a為均勻分布參數(shù)的下限，b為均勻分布參數(shù)的上限。

2.1.2 正態(tài)分布參數(shù)的簡單隨機(jī)抽樣

對于正態(tài)分布參數(shù)，簡單隨機(jī)抽樣實現(xiàn)步驟如下。

1）通過公式（1）采用Python3.7 程序中的random 模塊調(diào)用random（）方法生成0 到1 之間的隨機(jī)浮點數(shù)u；

2）將u視為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)值，采用公式（3）借助標(biāo)準(zhǔn)正態(tài)分布反函數(shù)映射到上1-u分位點z；

式中：z為標(biāo)準(zhǔn)正態(tài)分布的上1-u分位點，?-1為標(biāo)準(zhǔn)正態(tài)分布的反函數(shù)。

3）采用線性變換公式（4）將z映射為正態(tài)分布抽樣散點。

式中：σ 為正態(tài)分布的標(biāo)準(zhǔn)差，μ 為正態(tài)分布的期望值。

2.2 拉丁超立方抽樣

2.2.1 均勻分布參數(shù)的拉丁超立方抽樣

對于均勻分布參數(shù)，拉丁超立方抽樣[3]實現(xiàn)步驟如下。

1）將均勻分布參數(shù)的取值范圍等概率分割為滿足公式（5）的N個子區(qū)間；

式中：X為隨機(jī)變量，Ii為等概率區(qū)間坐標(biāo)值。

2）采用公式（6）將0 到1 之間的隨機(jī)浮點數(shù)u映射為歸一化區(qū)間的各子區(qū)間的隨機(jī)數(shù)；

3）通過公式（7）將隨機(jī)數(shù)zi映射為參數(shù)分布范圍內(nèi)的隨機(jī)數(shù)；

2.2.2 正態(tài)分布參數(shù)的拉丁超立方抽樣[3]

對于正態(tài)分布參數(shù)，拉丁超立方抽樣實現(xiàn)步驟如下。

1）依據(jù)公式（8）將正態(tài)分布參數(shù)區(qū)間進(jìn)行等概率劃分

2）通過公式（9）將0 到1 之間的隨機(jī)浮點數(shù)u映射為歸一化區(qū)間的各子區(qū)間的隨機(jī)數(shù)；

3）將νi視為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)值，采用公式（10）借助標(biāo)準(zhǔn)正態(tài)分布反函數(shù)映射到上1-νi分位點wi；

4）采用線性變換公式（11）將wi映射為正態(tài)分布抽樣散點。

2.3 Wilks 非參數(shù)統(tǒng)計

Wilks 非參數(shù)統(tǒng)計方法[4]廣泛應(yīng)用于BEPU 統(tǒng)計分析中?；舅枷胧峭ㄟ^Wilks 公式根據(jù)選定的置信水平和概率水平獲取最小的計算次數(shù)，通過對響應(yīng)參數(shù)的有序統(tǒng)計理論來獲取統(tǒng)計容忍區(qū)間。

對于概率密度分布函數(shù)g（y），確定一統(tǒng)計容忍區(qū)間（L，U），此區(qū)間對y的覆蓋率至少為γ 的概率為β，通過式（12）表示為

式中：β 為置信水平，γ 為概率水平，L為容忍下限，U為容忍上限，P為概率。

依據(jù)有序統(tǒng)計理論和多項式分布定律可推導(dǎo)出單側(cè)統(tǒng)計容忍區(qū)間Wilks 公式為

式中：N為計算次數(shù)。

Wilks 非參數(shù)統(tǒng)計方法具有黑箱統(tǒng)計模型的特點，由公式（13）可知：計算次數(shù)僅與置信水平和概率水平有關(guān)，與輸入變量的數(shù)量及其概率密度分布無關(guān)。Guba[5]在2003 年將單側(cè)統(tǒng)計容忍區(qū)間Wilks 公式推廣到下式

式中：p為輸出變量數(shù)量。

依據(jù)Wilks 非參數(shù)統(tǒng)計方法，采取95%置信水平和95%概率水平（簡記為95/95），若只有一個輸出變量，確定95/95 單側(cè)統(tǒng)計容忍上限至少需進(jìn)行59 次運(yùn)算，取59 次運(yùn)算最大值；若計算次數(shù)為93 次，取次大值作為95/95 單側(cè)統(tǒng)計容忍上限；若計算次數(shù)為124次，取第三大值作為95/95 單側(cè)統(tǒng)計容忍上限。

3 不確定性量化與敏感性分析

3.1 輸入?yún)?shù)抽樣

主要依據(jù)大破口失水事故PIRT 表[6]，綜合考慮大破口失水事故過程中燃料棒、反應(yīng)堆堆芯、換料水箱和安注箱等重要部件重要現(xiàn)象，識別大破口失水事故下對關(guān)鍵安全參數(shù)PCT 有重要影響的模型參數(shù)、電廠參數(shù)（電廠初始條件和電廠邊界條件），并結(jié)合相關(guān)資料等多種方式綜合確定不確定性輸入?yún)?shù)。選用的不確定性輸入?yún)?shù)見表3。對各種不確定性輸入?yún)?shù)的范圍和概率密度分布，結(jié)合核電站設(shè)計文件、已有的相關(guān)文獻(xiàn)資料及工程經(jīng)驗判斷加以判定。在信息不足的情況下，不確定性輸入?yún)?shù)通常假設(shè)為均勻分布。

表3 不確定性輸入?yún)?shù)

采用Python3.7 程序調(diào)用Mersenne Twister 隨機(jī)數(shù)產(chǎn)生器。對不確定性輸入?yún)?shù)分別進(jìn)行124 次簡單隨機(jī)抽樣和拉丁超立方抽樣。

可見抽樣散點反映了輸入?yún)?shù)的分布特征，對于均勻分布的輸入?yún)?shù)，散點在取值范圍內(nèi)分布均勻。對于正態(tài)分布參數(shù)，散點分布呈現(xiàn)出中間聚集效應(yīng)。拉丁超立方抽樣比簡單隨機(jī)抽樣散點分布更加均勻。以堆芯功率散點分布為例，簡單隨機(jī)抽樣相鄰散點距離的標(biāo)準(zhǔn)差為0.517 7，而拉丁超立方抽樣相鄰散點距離的標(biāo)準(zhǔn)差僅為0.252 6。

3.2 不確定性量化

將簡單隨機(jī)抽樣和拉丁超立方抽樣結(jié)果組合作為RELAP5/MOD3.4 程序輸入并執(zhí)行計算，得到基于2 種抽樣的124 組熱點處燃料包殼溫度瞬態(tài)曲線。通過AptPlot 程序提取熱點處包殼峰值溫度數(shù)據(jù)。相應(yīng)的包殼峰值溫度散點分布如圖6 和圖7 所示。

圖6 基于SRS 的124 組PCT 散點

圖7 基于LHS 的124 組PCT 散點圖

將124 組包殼峰值溫度由大到小進(jìn)行排序，根據(jù)Wilks 非參數(shù)統(tǒng)計理論可知，排序第3 位的值即為包殼峰值溫度, 在95%置信水平上具有95%概率水平的單側(cè)統(tǒng)計容忍上限?；诤唵坞S機(jī)抽樣計算得到包殼峰值溫度單側(cè)統(tǒng)計容忍上限95/95PCT 為1 264.057 K，基于拉丁超立方抽樣計算得到包殼峰值溫度單側(cè)統(tǒng)計容忍上限95/95PCT 為1 267.617 K?；? 種抽樣方法得到的包殼峰值溫度，單側(cè)統(tǒng)計容忍上限均低于驗收準(zhǔn)則規(guī)定的限值1 477 K?；? 種抽樣方法得到的包殼峰值溫度單側(cè)統(tǒng)計容忍上限較為接近，僅相差3.56K。

3.3 敏感性分析

Spearman 秩相關(guān)系數(shù)法作為全局敏感性分析方法[7]，可檢驗多個輸入?yún)?shù)對響應(yīng)參數(shù)的全局影響。采用Spearman 秩相關(guān)系數(shù)作為輸入?yún)?shù)對包殼峰值溫度敏感性的度量。Spearman 秩相關(guān)系數(shù)的絕對值大小表示敏感性強(qiáng)弱，其計算公式為

式中：RXi為Xi在X中的大小排序，RYi為Yi在Y中的大小排序。

基于簡單隨機(jī)抽樣和拉丁超立方抽樣的不確定性輸入?yún)?shù)與包殼峰值溫度間的Spearman 秩相關(guān)系數(shù)如圖8 和圖9 所示。

圖8 基于SRS 輸入?yún)?shù)的Spearman 秩相關(guān)系數(shù)

圖9 基于LHS 輸入?yún)?shù)的Spearman 秩相關(guān)系數(shù)

基于簡單隨機(jī)抽樣的敏感性分析結(jié)果表明包殼峰值溫度對衰變熱、換料水箱水溫、安注箱水溫和氣隙導(dǎo)熱率較為敏感；基于拉丁超立方抽樣的敏感性分析結(jié)果表明，包殼峰值溫度對衰變熱、換料水箱水溫、安注箱水溫和氣隙導(dǎo)熱率較為敏感。因此采用Spearman 秩相關(guān)系數(shù)法基于2 種不同的抽樣方法識別出的包殼峰值溫度四大主要影響參數(shù)相一致。

4 結(jié)論

以百萬千瓦級壓水堆核電廠大破口失水事故為分析對象，基于簡單隨機(jī)抽樣和拉丁超立方抽樣對包殼峰值溫度開展不確定性量化和敏感性分析計算，得到的結(jié)論如下。

1）采用簡單隨機(jī)抽樣和拉丁超立方抽樣對不確定性輸入?yún)?shù)開展隨機(jī)抽樣，拉丁超立方抽樣散點比簡單隨機(jī)抽樣散點分布更加均勻，能夠更好地復(fù)現(xiàn)參數(shù)的分布特征。

2）采用Wilks 非參數(shù)統(tǒng)計理論基于簡單隨機(jī)抽樣和拉丁超立方抽樣計算得到的95/95 PCT 均滿足ECCS 驗收準(zhǔn)則，且95/95 PCT 非常接近，僅相差3.56 K。

3）采用Spearman 秩相關(guān)系數(shù)法基于簡單隨機(jī)抽樣和拉丁超立方抽樣均識別出包殼峰值溫度對衰變熱、換料水箱水溫、安注箱水溫和氣隙導(dǎo)熱率較為敏感。