陳仁宗，王 輝，孫曉暉

(1.清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084；2.中國核電工程有限公司，北京 100840)

1988年，美國核管會(NRC)修訂了美國聯(lián)邦法規(guī)10CFR50.46允許使用最佳估算方法進行大破口事故的認證級分析，但必須考慮不確定性，并加以量化計算，以保證分析結(jié)果在驗收準(zhǔn)則之內(nèi)具有較高的概率[1-4]。最佳估算方法應(yīng)符合以下3個條件：應(yīng)根據(jù)選定的驗收準(zhǔn)則，事故分析不引入有意的悲觀性；使用最佳估算程序；包含不確定性分析[5]。最佳估算方法利用盡可能詳細的模型而非簡單模型來保證結(jié)果更為接近物理現(xiàn)實，用最佳估算結(jié)果的不確定性來量化分析結(jié)果與物理現(xiàn)實之間的差距[6]。

目前，最佳估算加不確定性(BEPU)方法主要應(yīng)用于設(shè)計基準(zhǔn)事故的安全分析[7-9]。與設(shè)計基準(zhǔn)事故相比，嚴(yán)重事故現(xiàn)象更加復(fù)雜，涉及堆芯裸露升溫、燃料包殼氧化升溫失效、堆芯熔毀和裂變產(chǎn)物釋放等諸多過程。復(fù)雜的嚴(yán)重事故現(xiàn)象導(dǎo)致了嚴(yán)重事故更大的不確定性和更為復(fù)雜的不確定性分析計算過程。嚴(yán)重事故現(xiàn)象復(fù)雜，其不確定性較大，目前尚無官方發(fā)布的嚴(yán)重事故現(xiàn)象識別排序(PIRT)表。本文采用的方法是結(jié)合設(shè)計基準(zhǔn)事故PIRT表中參數(shù)和與嚴(yán)重事故現(xiàn)象相關(guān)的參數(shù)進行嚴(yán)重事故輸入?yún)?shù)的不確定性分析。

堆芯出口溫度(CET)是核電廠安全運行的重要監(jiān)測參數(shù)，在嚴(yán)重事故管理導(dǎo)則(SAMG)中作為堆芯狀態(tài)的重要表征參數(shù)，可作為堆芯損傷評價和簡化源項評估的整定值，并可作為嚴(yán)重事故管理導(dǎo)則的入口條件。同時，考慮堆芯裸露后復(fù)雜的嚴(yán)重事故現(xiàn)象，文中并未針對包殼破裂后的嚴(yán)重事故進程進行不確定性分析。包殼破裂通常被作為裂變產(chǎn)物氣隙釋放過程的開始，在堆芯損傷評估和簡化源項評估過程中扮演重要的角色。

本文以VVER1000壓水堆核電廠為研究對象，選取大破口失水事故(LBLOCA)始發(fā)嚴(yán)重事故進行包殼破裂對應(yīng)堆芯出口溫度的BEPU分析。

1 計算模型

本文以VVER1000壓水堆核電廠為研究對象，BEPU計算流程如圖1所示。第1步，建立VVER1000的最佳估算模型，并對模型進行驗證；第2步，選取不確定性分析的輸入?yún)?shù)，并對輸入?yún)?shù)進行取樣，根據(jù)取樣結(jié)果進行計算；第3步，整理計算結(jié)果，進行不確定性和敏感性分析。

圖1 BEPU計算流程Fig.1 Flowchart of BEPU calculation

本文采用的最佳估算程序為MELCOR。MELCOR是NRC委托桑迪亞國家實驗室開發(fā)的用于概率風(fēng)險評估的程序[10]，主要用于嚴(yán)重事故分析，如SAMG中時間窗口的計算、氫氣風(fēng)險和源項評估等，目前也是監(jiān)管機構(gòu)主要的安全評審工具。MELCOR的運行需要兩步，第1步執(zhí)行MELGEN進行輸入檢查，第2步啟動MELCOR進行迭代計算。

1.1 VVER1000模型建立與瞬態(tài)計算

VVER1000壓水堆核電廠是由俄羅斯研發(fā)的第3代核電技術(shù)[11]。核電廠主回路系統(tǒng)包括4個環(huán)路，每個環(huán)路有1臺主泵和1臺直流式蒸汽發(fā)生器，共用1臺穩(wěn)壓器。反應(yīng)堆熱功率為3 000 MW，堆芯出口處壓力為15.7 MPa，入口冷卻劑溫度為291 ℃。采用MELCOR對VVER1000壓水堆核電廠的設(shè)計特點建立最佳估算模型，系統(tǒng)節(jié)點圖如圖2所示。建模對象包括反應(yīng)堆壓力容器、穩(wěn)壓器、直流式蒸汽發(fā)生器、主泵、主給水系統(tǒng)和主管道等。其中壓力容器包括下降段、下封頭、下腔室、堆芯區(qū)域、旁通段和上腔室的模擬。

在進行瞬態(tài)計算前要對模型進行穩(wěn)態(tài)調(diào)試，確保模擬值與設(shè)計值一致。冷卻劑流量、堆芯入口溫度、堆芯出口溫度、堆芯壓力的模擬值與設(shè)計值對比如圖3所示，模擬值與設(shè)計值的相對誤差列于表1，由圖3和表1可知，模擬值與設(shè)計值之間的誤差較小，相對誤差均在1%以內(nèi)，說明所建模型可信，可用于進一步的計算。

圖2 主回路系統(tǒng)節(jié)點圖Fig.2 Node diagram of main circuit system

LBLOCA是輕水堆安全設(shè)計的基準(zhǔn)事故[12]。國際原子能機構(gòu)把LBLOCA始發(fā)嚴(yán)重事故作為計算簡化源項的基準(zhǔn)事故[13]，國內(nèi)及美國NRC將NUREG-1465報告中LBLOCA始發(fā)嚴(yán)重事故的源項作為核應(yīng)急的參考源項[14]。本文選取LBLOCA作為始發(fā)事故，事故假設(shè)條件如下：0 s時刻，反應(yīng)堆冷管段發(fā)生當(dāng)量直徑為30 cm的破口；余熱排出系統(tǒng)不可用；安注箱的非能動部分可用，能動部分失效；蒸汽發(fā)生器二次側(cè)輔助給水喪失。

事故發(fā)生后，大量冷卻劑泄漏至安全殼內(nèi)，導(dǎo)致安全殼內(nèi)壓力升高，同時一回路壓力迅速降低，從而觸發(fā)反應(yīng)堆緊急停堆及主泵停運。壓力降低到一定值時，安注箱非能動部分開始向堆芯注水，隨著應(yīng)急冷卻水的不足，堆芯開始裸露升溫，然后包殼破裂，開始氣隙釋放。大破口失水事故序列列于表2。

圖3 不同參數(shù)的模擬值與設(shè)計值對比Fig.3 Comparison of simulation value and design value of different parameters

表1 模擬值與設(shè)計值的相對誤差Table 1 Relative error between simulation value and design value

表2 事故序列Table 2 Accident sequence

1.2 不確定性分析方法

本文采用拉丁超立方抽樣方式對輸入?yún)?shù)進行抽樣。與簡單隨機抽樣相比，拉丁超立方抽樣具有更高的抽樣維度，抽取的樣品參數(shù)更具代表性。拉丁超立方抽樣的實現(xiàn)步驟如下：1) 將參數(shù)按其范圍進行等概率劃分，劃分為很多子范圍；2) 用隨機數(shù)產(chǎn)生器產(chǎn)生指定數(shù)量的[0，1]范圍內(nèi)的隨機數(shù)；3) 根據(jù)均勻分布或正態(tài)分布將隨機數(shù)映射為實際范圍內(nèi)的參數(shù)。

本文采用Wilks非參數(shù)統(tǒng)計方法[15-16]進行不確定性量化計算，此方法廣泛應(yīng)用于BEPU統(tǒng)計分析中。其基本思想是通過Wilks公式根據(jù)選定的置信水平和概率水平獲取最小的計算次數(shù)，通過對響應(yīng)參數(shù)的有序統(tǒng)計理論來獲取統(tǒng)計容忍區(qū)間。Wilks公式為：

β=1-γN

(1)

式中：β為置信水平；γ為概率水平；N為計算次數(shù)。

基于Wilks公式，對于單側(cè)統(tǒng)計容忍限值(95/95)，樣本數(shù)決定了限值的選取原則：樣本為59組，則許用限值為59組計算結(jié)果的最大值；若樣本為93組，則許用限值為93組計算結(jié)果中的第2大值；若樣本為124組，則許用限值為124組計算結(jié)果中的第3大值。本文選取的樣本數(shù)為59。

2 輸入?yún)?shù)選取與抽樣

2.1 輸入?yún)?shù)選取

結(jié)合用于設(shè)計基準(zhǔn)事故的LBLOCA PIRT表[17-18]、與鋯包殼破裂相關(guān)的嚴(yán)重事故現(xiàn)象和專家判斷進行輸入?yún)?shù)的確定。綜合考慮核電廠設(shè)計文件、工程經(jīng)驗判斷和相關(guān)文獻資料中的參考值進行輸入?yún)?shù)范圍和概率密度分布的判定。不確定性輸入?yún)?shù)及其范圍、分布列于表3。

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

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

輸入?yún)?shù)的抽樣采用DAKOTA程序[19]實現(xiàn)，選用Mersenne Twister隨機數(shù)產(chǎn)生器，種子設(shè)定值為98 765。對不確定性輸入?yún)?shù)分別進行59次拉丁超立方抽樣，歸一化結(jié)果如圖4所示。

圖4 不確定性輸入?yún)?shù)歸一化分布Fig.4 Normalized distribution of uncertainty input parameter

由圖4可看出，參數(shù)的歸一化結(jié)果與參數(shù)分布形式有關(guān)，均勻分布參數(shù)的歸一化結(jié)果呈現(xiàn)均勻分布，正態(tài)分布參數(shù)的歸一化結(jié)果呈現(xiàn)中間聚集效應(yīng)。輸入?yún)?shù)抽樣結(jié)果達到了預(yù)期的效果。

3 結(jié)果分析

3.1 不確定性量化

對59組MELCOR計算結(jié)果進行分析整理，得到氣隙釋放時間和堆芯出口溫度，如圖5所示。根據(jù)Wilks非參數(shù)統(tǒng)計方法，對于單側(cè)統(tǒng)計容忍限值(95/95)，樣本為59組，則許用限值為59組計算結(jié)果的最大值。據(jù)此得到的氣隙釋放時間和堆芯出口溫度單側(cè)統(tǒng)計容忍限值分別為1 905.71 s和430.85 ℃。

3.2 敏感性分析

Spearman秩相關(guān)系數(shù)法是全局敏感性方法，可檢驗多個輸入?yún)?shù)對響應(yīng)參數(shù)的全局影響。采用Spearman秩相關(guān)系數(shù)法作為輸入?yún)?shù)對堆芯出口溫度敏感性進行分析，相關(guān)系數(shù)取值范圍為[-1，1]，絕對值的大小表示敏感性的強弱；數(shù)值的正負號表示正負相關(guān)性，正號表示正相關(guān)性，負號表示負相關(guān)性，零表示不相關(guān)，其表達式如式(2)所示。

圖5 氣隙釋放時間和堆芯出口溫度Fig.5 Gas release time and CET

ρ=

(2)

式中：RXi為Xi在X中的大小排序；RYi為Yi在Y中的大小排序；ρ為秩相關(guān)系數(shù)；n為樣本數(shù)。

輸入?yún)?shù)與堆芯出口溫度間的Spearman秩相關(guān)系數(shù)如圖6所示，計算結(jié)果表明，氣隙釋放對應(yīng)的堆芯出口溫度對衰變熱系數(shù)和包殼厚度較為敏感。

圖6 輸入?yún)?shù)與堆芯出口溫度的 Spearman秩相關(guān)系數(shù)Fig.6 Spearman rank relational coefficient of input parameter and CET

4 結(jié)論

以VVER1000壓水堆核電廠LBLOCA始發(fā)嚴(yán)重事故為研究對象，MELCOR為研究工具，開展了輸入?yún)?shù)對堆芯出口溫度的不確定性和敏感性研究，得到的結(jié)論如下：1) 氣隙釋放時間和堆芯出口溫度單側(cè)統(tǒng)計容忍限值(95/95)分別為1 905.71 s和430.85 ℃；2) 氣隙釋放對應(yīng)的堆芯出口溫度對衰變熱系數(shù)和包殼厚度較為敏感。