陳仁宗，王 輝，孫曉暉

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

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

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

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

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

1 計(jì)算模型

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

圖1 BEPU計(jì)算流程Fig.1 Flowchart of BEPU calculation

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

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

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

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

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

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

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

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

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

表2 事故序列Table 2 Accident sequence

1.2 不確定性分析方法

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

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

β=1-γN

(1)

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

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

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

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

結(jié)合用于設(shè)計(jì)基準(zhǔn)事故的LBLOCA PIRT表[17-18]、與鋯包殼破裂相關(guān)的嚴(yán)重事故現(xiàn)象和專家判斷進(jìn)行輸入?yún)?shù)的確定。綜合考慮核電廠設(shè)計(jì)文件、工程經(jīng)驗(yàn)判斷和相關(guān)文獻(xiàn)資料中的參考值進(jì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]實(shí)現(xiàn)，選用Mersenne Twister隨機(jī)數(shù)產(chǎn)生器，種子設(shè)定值為98 765。對(duì)不確定性輸入?yún)?shù)分別進(jìn)行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é)果達(dá)到了預(yù)期的效果。

3 結(jié)果分析

3.1 不確定性量化

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

3.2 敏感性分析

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

圖5 氣隙釋放時(shí)間和堆芯出口溫度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所示，計(jì)算結(jié)果表明，氣隙釋放對(duì)應(yīng)的堆芯出口溫度對(duì)衰變熱系數(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)重事故為研究對(duì)象，MELCOR為研究工具，開展了輸入?yún)?shù)對(duì)堆芯出口溫度的不確定性和敏感性研究，得到的結(jié)論如下：1) 氣隙釋放時(shí)間和堆芯出口溫度單側(cè)統(tǒng)計(jì)容忍限值(95/95)分別為1 905.71 s和430.85 ℃；2) 氣隙釋放對(duì)應(yīng)的堆芯出口溫度對(duì)衰變熱系數(shù)和包殼厚度較為敏感。