適用于核電廠BEPU的高效全局敏感性分析方法開發(fā)及應(yīng)用

熊青文，茍軍利，杜 鵬，鄧 堅，邱志方，黃 濤，申亞歐

(1.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都 610213；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

近年來，隨著世界核電技術(shù)的快速發(fā)展，核安全問題也愈加受到重視。2003年，國際原子能機構(gòu)(IAEA)發(fā)表的報告總結(jié)了各種事故分析方法的分類，其將安全分析方法分為了4個選項，其中選項3，即使用最佳估算程序計算帶有不確定性的真實輸入數(shù)據(jù)的方法，也被稱為最佳估算加不確定性(BEPU)分析方法，是目前核電安全分析的主流方法[1]。

BEPU方法中一個關(guān)鍵的步驟為參數(shù)的敏感性分析，其旨在評估輸入?yún)?shù)對輸出參數(shù)的影響大小。敏感性分析方法可根據(jù)其特征分為兩類，即局部敏感性分析和全局敏感性分析[2]。局部敏感性分析的實現(xiàn)相對簡單，對于簡單的線性系統(tǒng)或單調(diào)系統(tǒng)，局部方法能用于評估輸入?yún)?shù)的敏感性度量，但對于復(fù)雜的大型核反應(yīng)堆系統(tǒng)，由于其存在高度非線性及參數(shù)間的相互作用，使用局部敏感性分析方法可能得到不充分的結(jié)果，因此發(fā)展全局敏感性分析方法很有必要。

常規(guī)用于工程應(yīng)用的全局方法有基礎(chǔ)效應(yīng)方法、方差分解方法、矩獨立方法等[3]。其中矩獨立敏感性度量旨在評估輸入?yún)?shù)對輸出參數(shù)概率密度函數(shù)(PDF)的影響，通過計算輸出參數(shù)的無條件PDF及固定一個輸入?yún)?shù)的條件PDF之間的偏差值來評估輸入?yún)?shù)對輸出參數(shù)的影響。矩獨立重要性測度由于其廣泛的適用性而在近年來得到了大量工程應(yīng)用[4-5]。

實際情況中應(yīng)用矩獨立方法最大的困難在于其需要大量的計算來估計輸出參數(shù)的PDF，其計算過程涉及兩個循環(huán)抽樣計算，涉及的計算成本十分昂貴[6]。因此該方法往往被用于關(guān)系式模型或簡單系統(tǒng)的敏感性分析，很難應(yīng)用于核反應(yīng)堆失水事故分析等大型復(fù)雜且計算耗時的工況。因此，如何降低計算成本是將矩獨立全局敏感性分析方法應(yīng)用至核電廠BEPU分析的一個重要課題。本研究開發(fā)一個優(yōu)化的低成本全局矩獨立敏感性分析方法，通過多個例題對其可靠性進行驗證，并將其應(yīng)用于LOFT(loss-of-fluid test)大破口失水事故(LBLOCA)的敏感性分析。

1 方法優(yōu)化

矩獨立全局敏感性分析方法旨在評估輸入?yún)?shù)對目標輸出PDF的影響。假設(shè)某個模型Y=g(X)存在k個輸入?yún)?shù)，即X=(X1,X2,…,Xk)T，每個輸入?yún)?shù)均服從一定的概率分布fXi(xi)，其不確定性可通過模型計算傳播至輸出Y。將Y的無條件PDF和累積分布函數(shù)(CDF)表示為fY(y)和FY(y)，將第i個輸入?yún)?shù)Xi取某一固定值時得到的Y的條件PDF和CDF表示為fY|Xi(y)和FY|Xi(y)。根據(jù)定義，可將第i個輸入?yún)?shù)的矩獨立敏感性指標表示為：

(1)

其中，s(Xi)為固定第i個輸入?yún)?shù)情況下輸出PDF的偏移量，可表示為：

(2)

如上所示，由于計算參數(shù)的矩獨立敏感性指標涉及到兩個嵌套積分計算，且估計輸出的PDF是一個不適定的問題[7]，因此總體的計算量十分巨大。傳統(tǒng)方法的計算中總計算次數(shù)為R(kN1N2+N1)，其中N1為計算輸出的PDF的抽樣計算次數(shù)，推薦取值N1>1 000，N2為計算每個輸入?yún)?shù)影響期望值所需的抽樣計算次數(shù)，推薦取值N2>1 000，R為計算抽樣樣本方差所需的獨立重復(fù)計算次數(shù)，由于大部分工程例題不存在解析解，因此對于傳統(tǒng)方法，計算量越大，計算結(jié)果越精確。將求解輸出的PDF轉(zhuǎn)換為求解其CDF可在一定程度上優(yōu)化計算結(jié)果。假設(shè)輸出的fY(y)和fY|Xi(y)存在m個交點，表示為a1，a2，…，am，則s(Xi)可表示為m+1個子面積之和，即：

s(Xi)=s1+s2+…+sj+…+sm+sm+1

j=1,2,…,m+1

(3)

其中，fY(y)和fY|Xi(y)的交點和每個子面積sj可根據(jù)輸出的FY(y)和FY|Xi(y)計算得到。為了能以十分小的計算成本且又相對精確地計算輸入?yún)?shù)的敏感性指標，使用多種方法進行優(yōu)化計算。首先，為降低積分計算的計算量，使用五點高斯求積方案替代積分計算[8]。基于高斯求積方案，可將矩獨立敏感性指標表示為：

(4)

式中：ωi,j為第i個輸入?yún)?shù)按照其分布類型確定的5個高斯權(quán)重值中的第j個值；Xi,j為第i個輸入?yún)?shù)的第j個高斯點取值。ωi,j和Xi,j的取值與參數(shù)的分布類型及不確定性區(qū)間有關(guān)，具體取值可參考文獻[8]。

求解s(Xi,j)的關(guān)鍵在于求解輸出的條件和無條件CDF，即FY(y)和FY|Xi(y)。根據(jù)CDF的定義，可將其表示為：

FY(y)=P{Y≤y}=P{g(X)≤y}=

P{g(X)-y≤0}=P{z(X,y)≤0}=

Pf{z(X,y)}

(5)

其中：z(X,y)=g(X)-y為定義的新函數(shù)；Pf為失效概率。因此，可將求解模型g(X)的CDF轉(zhuǎn)換為求解模型z(X,y)的失效概率，可根據(jù)四階矩估計方法[9]和Pearson系統(tǒng)[10]來求解模型的失效概率。

由于采用五點高斯求積方案進行簡化計算，因此每個參數(shù)僅需計算5個高斯點處對應(yīng)的函數(shù)的輸出值，同時還需計算一次參考點處對應(yīng)的模型輸出值，因此該方法需調(diào)用的模型計算次數(shù)最多為5k+1次。此外，對于不確定性分布服從正態(tài)分布、均勻分布或?qū)?shù)正態(tài)分布的參數(shù)，其每個參數(shù)僅需執(zhí)行4次計算。

2 優(yōu)化方法驗證

為驗證所提出優(yōu)化的矩獨立方法的準確性和適應(yīng)性，本研究中將提出的優(yōu)化矩獨立方法應(yīng)用至1個數(shù)值線性關(guān)系例題及兩個工程非線性非單調(diào)例題，以驗證該方法在簡單線性系統(tǒng)及非線性非單調(diào)系統(tǒng)中的適應(yīng)性。由于工程中大多數(shù)問題不存在矩獨立指標的解析解，因此驗證過程中使用基于雙循環(huán)拉丁超立方抽樣(LHS)的方法來求解矩獨立敏感性指標的參考值，比較提出方法的計算值與參考值即可完成方法的驗證。

2.1 數(shù)值例題

首先使用了一個數(shù)值例題來對雙循環(huán)LHS方法和提出的方法進行驗證[11]。該數(shù)值例題的表達式可表示為：

g(X)=1.5X1+1.6X2+1.7X3+

1.8X4+1.9X5+2.0X6

(6)

式中6個輸入?yún)?shù)均服從標準正態(tài)分布N(0,1)。首先根據(jù)Borgonovo等[12]的工作可計算得到每個參數(shù)敏感性指標δi的解析解，而后使用雙循環(huán)LHS方法和提出的方法分別執(zhí)行1.201×109次和25次程序計算，如圖1所示。此外，本算例中亦參考了其他作者的研究結(jié)果，圖1中雙循環(huán)MC和單循環(huán)MC是Wei等的計算結(jié)果[11]，兩種方法的計算量分別為1.200 2×108次和6×104次。

圖1 數(shù)值例題分析結(jié)果Fig.1 Result of analytical test

由圖1可見，雙循環(huán)方法由于使用了LHS方法，且計算量較多，因此相比于Wei等的計算結(jié)果，其結(jié)果更加接近于解析解。就數(shù)值而言，雙循環(huán)LHS方法計算得到的敏感性指標相比于解析解相對誤差在0.03%～0.24%之間，故雙循環(huán)LHS方法的結(jié)果可作為參考解以驗證優(yōu)化方法的準確性。相對于雙循環(huán)LHS方法，提出方法計算得到的矩獨立度量誤差范圍在4.2%～6.4%之間，在計算量大幅減小的情況下依舊具有較高的計算精度。

雙循環(huán)LHS方法和提出的優(yōu)化方法計算結(jié)果間存在一個偏差，該偏差主要是由于使用高斯求積方案以及Pearson系統(tǒng)引入的。相對于敏感性指標δi的取值，該偏差相對較小，且其不影響參數(shù)的敏感性排序，因此提出方法的準確性可得到證明。

2.2 懸梁臂模型

懸梁臂模型是一個典型的結(jié)構(gòu)力學(xué)問題，其經(jīng)常被用于執(zhí)行參數(shù)的敏感性分析，該模型的介紹可參考文獻[11]。懸梁臂模型的表達式可表示為：

(7)

式中，參數(shù)X服從正態(tài)分布N(500, 402)，Y服從N(1 000, 802)，E服從N(2.9×107,2 320 0002)，w服從N(2.448 7,0.195 8962)，t服從N(3.888 4, 0.311 0722)，L服從N(100,82)。使用雙循環(huán)LHS方法和提出的方法分別執(zhí)行1.201×109次和25次程序計算，結(jié)果如圖2所示。

圖2 懸梁臂模型分析結(jié)果Fig.2 Result of cantilever beam structure

驗證結(jié)果顯示兩種方法結(jié)果十分接近，由于參數(shù)本身數(shù)值較小，因此計算得到的敏感性指標誤差在1.15%～17.27%之間，且參數(shù)敏感性排序一致，進一步驗證了優(yōu)化方法的可靠性。

2.3 機翼重量函數(shù)

機翼重量函數(shù)旨在模擬計算輕型飛機的機翼重量，該問題也常被用于參數(shù)篩選的示例，其相關(guān)描述可參考文獻[13]。機翼重量函數(shù)的函數(shù)表達式為：

(8)

式(8)中各變量的意義及分布列于表1，需要說明的是，由于不清楚參數(shù)的具體分布，因此假設(shè)所有參數(shù)均服從均勻分布。

表1 機翼重量函數(shù)輸入?yún)?shù)及其分布Table 1 Input parameter and distribution of wing weight function

使用雙循環(huán)LHS方法和提出的方法分別執(zhí)行2.001×109次和41次程序計算，結(jié)果如圖3所示。

圖3 機翼重量函數(shù)分析結(jié)果Fig.3 Result of wing weight function

由圖3可見，所提出優(yōu)化方法的結(jié)果與雙循環(huán)LHS方法的結(jié)果依舊十分接近，兩種方法敏感性排序都是一致的。由于第2個參數(shù)對輸出的影響小，因此基于降維方法計算得到的敏感性指標非常小，這是由于算法導(dǎo)致的，但不影響敏感性的排序。對于重要參數(shù)，優(yōu)化方法計算得到的敏感性指標相比于LHS方法誤差在0.8%～1.7%之間，由此可見，所提出的方法能以非常少的計算量，很好地估算出輸入?yún)?shù)對輸出的矩獨立敏感性指標。

3 LOFT應(yīng)用

如上所述，優(yōu)化方法能以非常小的計算成本評估兩個矩獨立敏感性度量，且其在簡單線性系統(tǒng)、非線性以及非單調(diào)系統(tǒng)中的適用性和準確性得到了評估和驗證。為進一步驗證優(yōu)化方法是否適用于復(fù)雜的核反應(yīng)堆系統(tǒng)，本研究中以LOFT裝置LP-02-6 LBLOCA試驗為對象開展了相關(guān)研究。

LOFT試驗項目最早由美國NRC負責(zé)執(zhí)行，后期成為由OECD成員國合作的國際化項目[14]。該試驗項目探究了核反應(yīng)堆在發(fā)生各種事故工況時的表現(xiàn)，其中以LOCA為主要分析對象。LOFT裝置被設(shè)計為模擬典型四回路壓水堆，是一個IET縮比試驗裝置，主要包括1個壓力容器，1個用于模擬破口發(fā)生的破裂回路，1個包含有主泵、穩(wěn)壓器以及蒸汽發(fā)生器(SG)的完整回路，以及1個用于模擬安全殼收容作用的馳壓容器。LOFT裝置上開展了諸多試驗，其中LP-02-6試驗為模擬一回路冷管段雙端斷裂的LBLOCA，并模擬了失去廠外電源且只有1個ECCS投入使用的情況。

由于LBLOCA中最重要的驗收準則參數(shù)為包殼峰值溫度(PCT)，因此以PCT為目標輸出，開展LP-02-6 LBLOCA中輸入?yún)?shù)對PCT的矩獨立敏感性度量計算。基于現(xiàn)象識別排序表(PIRT)篩選了15個重要輸入?yún)?shù)[15]，這些參數(shù)及其不確定性分布列于表2。

表2 LOFT LP-02-6輸入?yún)?shù)Table 2 Input parameter of LOFT LP-02-6 scenario

為驗證提出的優(yōu)化方法在核反應(yīng)堆系統(tǒng)中的適用性，分別使用雙循環(huán)LHS方法和優(yōu)化方法執(zhí)行了矩獨立敏感性分析計算。其中，考慮到LOFT LP-02-6試驗的模擬計算較為耗時，因此雙循環(huán)LHS方法最終確定抽樣次數(shù)為1.001×106，使用高性能服務(wù)器5個并行通道耗時約為1個月執(zhí)行了這些計算。而對于優(yōu)化方法，總計算次數(shù)僅為81，總計算耗時不到20 min。

雙循環(huán)LHS方法和優(yōu)化方法的計算結(jié)果如圖4所示。

由計算結(jié)果可知，優(yōu)化方法與雙循環(huán)LHS方法的計算結(jié)果十分接近，對于重要參數(shù)，優(yōu)化方法計算得到的敏感性指標相比于LHS方法誤差在0.8%～31.4%之間，參數(shù)重要度排序亦高度相似，僅有3個重要度較小的參數(shù)敏感性排序不同，即衰變功率、穩(wěn)壓器水溫和完整回路冷管段溫度。因此，優(yōu)化方法在計算影響較大參數(shù)的敏感性度量時相當(dāng)準確，而對于影響較小的參數(shù)則可能存在一定的偏差，該計算精度可通過提高計算次數(shù)提高。考慮到核電廠事故中的敏感性分析主要目的在于確定工況中的重要參數(shù)，且考慮到諸多工況的單次模擬耗時較長，因此大多數(shù)情況下該優(yōu)化方法能適用于核電廠的敏感性分析。

4 結(jié)論

針對局部敏感性分析方法原理不適用，而全局敏感性分析方法計算成本過大難以應(yīng)用至核反應(yīng)堆系統(tǒng)的問題，本研究以矩獨立全局敏感性分析方法為對象，通過使用高斯求積方案、Pearson系統(tǒng)等方法開展了該方法的優(yōu)化研究。研究中提出的優(yōu)化方法的計算成本僅為輸入?yún)?shù)數(shù)目的4～5倍，因此成本十分低，可被用于計算耗時的核反應(yīng)堆系統(tǒng)的事故安全分析。

為驗證優(yōu)化方法的可靠性和適用性，首先使用一個簡單線性系統(tǒng)以及懸梁臂模型、機翼重量函數(shù)兩個非線性系統(tǒng)對其進行驗證。為獲得準確的參考解以評估優(yōu)化方法的可靠性，本研究中使用了基于大量抽樣計算的雙循環(huán)LHS方法。通過對比雙循環(huán)LHS方法和優(yōu)化方法的計算結(jié)果，驗證了優(yōu)化方法在以上系統(tǒng)中的準確性和適用性。而后，使用LOFT LP-02-6 LBLOCA工況驗證優(yōu)化方法在核反應(yīng)堆系統(tǒng)中的可靠性。通過對比雙循環(huán)LHS方法和提出優(yōu)化方法的計算結(jié)果，可得出優(yōu)化方法能以十分小的計算成本準確識別工況中的重要參數(shù)的結(jié)論。