王洋洋，孫曉暉，石興偉

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

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

圖1 基于輸入傳遞的不確定性評(píng)估流程

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

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

1 模型建立

1.1 程序介紹

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

1.2 模型建立

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

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

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

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

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

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

1.4 瞬態(tài)計(jì)算

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

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

圖3 堆芯功率

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

圖5 破口流量

2 理論方法

2.1 簡(jiǎn)單隨機(jī)抽樣

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 拉丁超立方抽樣

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

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

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

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

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

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

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

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

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

2）通過(guò)公式（9）將0 到1 之間的隨機(jī)浮點(diǎn)數(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分位點(diǎn)wi；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 不確定性量化

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

圖6 基于SRS 的124 組PCT 散點(diǎn)

圖7 基于LHS 的124 組PCT 散點(diǎn)圖

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

3.3 敏感性分析

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

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

基于簡(jiǎn)單隨機(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ù)

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

4 結(jié)論

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

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

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

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