熊青文，黃 濤，茍軍利，杜 鵬，鄧 堅，袁 鵬，周佳樾，胡文楨

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

最佳估算加不確定性(BEPU)分析方法是國際原子能機構(gòu)(IAEA)推薦用于核電廠安全分析和執(zhí)照申請的先進方法[1]，該類方法的一般步驟為：首先基于確定的目標輸出(FOM)建立現(xiàn)象識別排序表(PIRT)識別重要現(xiàn)象、模型或參數(shù)，并以帶有真實不確定性分布的輸入表示這些不確定性源；而后使用最佳估算程序建立電廠的模型，并執(zhí)行多次程序計算將輸入?yún)?shù)的不確定性傳播至FOM；最后基于程序計算結(jié)果量化FOM的不確定性，并執(zhí)行相應(yīng)的不確定性分析和敏感性分析。由于核電廠的結(jié)構(gòu)十分復(fù)雜，因此在模擬電廠時會涉及數(shù)量龐大的輸入?yún)?shù)，如初始/邊界條件、材料物性、狀態(tài)參數(shù)、本構(gòu)模型等[2]。在BEPU分析中考慮所有參數(shù)的不確定性是不現(xiàn)實的，因此需要基于一定的方法篩選和識別在分析工況中比較重要的參數(shù)以簡化分析。其中，PIRT是最常使用的一種方法。通常而言，PIRT的建立需要使用一定的方法，如AHP[3]或MCDA[4]等，并在專家經(jīng)驗的指導(dǎo)下使用定性的方法表征現(xiàn)象或參數(shù)的重要度，如使用H、M和L分別表征高、中、低重要度[5]。PIRT能初步識別事故分析過程中的重要現(xiàn)象、模型或參數(shù)，但是在實際應(yīng)用中會存在一些問題。首先，PIRT建立的成本很大。由于不同反應(yīng)堆在不同工況下的事故序列和熱工水力現(xiàn)象存在區(qū)別，因此理論上針對不同的反應(yīng)堆堆型和不同的事故工況均要建立對應(yīng)的PIRT，導(dǎo)致很多特定的反應(yīng)堆不存在完善的PIRT，難以開展BEPU分析而使用保守參數(shù)進行保守分析。此外，PIRT的建立高度依賴于專家經(jīng)驗，在專家經(jīng)驗不足的情況下往往采取的措施是參考類似反應(yīng)堆的已有PIRT，這種方式存在可能將不重要的參數(shù)考慮在內(nèi)，而將重要參數(shù)遺漏的問題。

核電廠安全分析中失水事故是最為關(guān)注的工況之一，相比于大破口失水事故而言，小破口失水事故(SBLOCA)的BEPU分析研究較少，相關(guān)的PIRT也存在不多。為了能在缺乏PIRT的情況下執(zhí)行準確的參數(shù)重要度排序計算，本文以典型三回路壓水堆(PWR)不同破口尺寸下的SBLOCA為對象開展相關(guān)工作，選定的FOM為堆芯活性區(qū)最低坍塌液位(MCL)。首先，分別在冷管段橫截面積的0.2%、0.4%、0.6%、0.8%和1.0%的5個破口尺寸下開展SBLOCA的瞬態(tài)基準計算，確定工況中的階段劃分，并確定堆芯MCL會出現(xiàn)于噴放和環(huán)路水封存在及清除階段兩個階段。隨后，使用優(yōu)化矩獨立全局敏感性分析方法對5個破口工況中兩個階段分別執(zhí)行定量敏感性分析，并將參數(shù)的重要度排序轉(zhuǎn)換為Savage分數(shù)，按照Savage分數(shù)將所有輸入?yún)?shù)進行重要度分組，從而得到PWR SBLOCA的參數(shù)重要度排序表。

1 優(yōu)化矩獨立方法介紹

本文使用的敏感性分析方法為矩獨立全局敏感性分析方法[6]，其旨在評估輸入?yún)?shù)對目標輸出概率密度函數(shù)(PDF)的影響，通過固定輸入?yún)?shù)的取值并計算目標輸出的無條件和條件PDF間的偏移量來定量評估參數(shù)的敏感性度量。假設(shè)某個模型Y=g(X)存在k個輸入?yún)?shù)，第i個輸入?yún)?shù)的矩獨立敏感性度量δi可表示為：

(1)

(2)

式中：fXi(xi)為第i個輸入?yún)?shù)的概率分布，xi為第i個輸入?yún)?shù)Xi在計算中的取值；s(Xi)為固定第i個輸入?yún)?shù)情況下輸出PDF的偏移量；fY(y)、fY|Xi(y)分別為Y的無條件和條件PDF，y為輸出參數(shù)Y在計算中的取值。

由于計算參數(shù)的矩獨立敏感性定量涉及到兩個嵌套積分計算，因此總體的計算量十分巨大。考慮到核電廠SBLOCA工況的模擬特別耗時，一次計算需耗費數(shù)小時，因此難以采用常規(guī)抽樣方法直接計算各輸入?yún)?shù)的矩獨立敏感性度量。為了能使計算成本特別小且又可相對精確地計算輸入?yún)?shù)的敏感性度量，本文使用一種優(yōu)化的矩獨立全局敏感性分析方法。首先，為了降低積分計算的計算量，使用五點高斯求積方案替代積分計算[7]。基于高斯求積方案，可將矩獨立敏感性度量表示為：

(3)

式中：ωij為第i個輸入?yún)?shù)按照其分布類型確定的5個高斯權(quán)重值中的第j個值；Xij為第i個輸入?yún)?shù)的第j個高斯點取值。因此求解度量δi的關(guān)鍵在于求解s(Xij)，即輸出的無條件PDF與將第i個輸入?yún)?shù)固定于其第j個高斯點時輸出條件PDF間的偏移量。

Liu等[8]提出了一種基于輸出Y的條件和無條件累計分布函數(shù)(CDF)計算s(Xij)的方法，計算得到輸出的CDF即可快速計算s(Xij)。根據(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)}

(4)

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

由于采用五點高斯求積方案簡化計算，因此每個參數(shù)僅需計算5個高斯點處對應(yīng)的函數(shù)輸出值，同時執(zhí)行1次所有參數(shù)的名義值計算，因此該方法需要調(diào)用的模型計算次數(shù)最多為5k+1次。傳統(tǒng)計算矩獨立敏感性度量所需的計算次數(shù)常為輸入?yún)?shù)數(shù)量的數(shù)萬至數(shù)百萬倍，相比于傳統(tǒng)基于隨機抽樣的矩獨立敏感性度量計算方法而言，本文提出的方法能夠以極小的計算量達到類似的計算精度。使用失水事故試驗(LOFT)中的LP-02-6工況對提出的方法進行驗證[11]，驗證計算中使用了15個重要輸入?yún)?shù)，分別使用提出方法和傳統(tǒng)抽樣方法計算輸入對包殼峰值溫度的影響，其中抽樣方法抽樣計算次數(shù)為1.001×106，而優(yōu)化方法總計算次數(shù)僅為81。兩種方法的計算結(jié)果如圖1所示，結(jié)果表明，在計算量縮小近1萬倍的情況下，對于大部分輸入?yún)?shù)該提出方法與傳統(tǒng)方法計算結(jié)果的相對偏差在0.8%～5.2%內(nèi)。

圖1 LOFT LP-02-6工況敏感性分析結(jié)果

2 SBLOCA瞬態(tài)分析

本文考慮的工況為典型PWR 5個不同破口尺寸下的SBLOCA，首先執(zhí)行所有參數(shù)取名義值的基準工況計算，并基于計算結(jié)果分析SBLOCA瞬態(tài)。計算結(jié)果表明，SBLOCA的事故進程可根據(jù)現(xiàn)象分為4個階段，這4個階段及其主要現(xiàn)象介紹如下。

1)噴放階段

小破口出現(xiàn)后，冷卻劑快速從破口處流失，導(dǎo)致一回路壓力快速下降。但由于反應(yīng)堆此時還未觸發(fā)停堆狀態(tài)，主泵還在持續(xù)運轉(zhuǎn)，一回路的冷卻劑處于受迫循環(huán)狀態(tài)。當反應(yīng)堆停堆后，保守假設(shè)主泵立即停止運行，主泵在惰轉(zhuǎn)一段時間后停止，此時一回路內(nèi)冷卻劑進入自然循環(huán)階段。

2)自然循環(huán)階段

主泵停止轉(zhuǎn)動后，一回路冷卻劑處于自然循環(huán)階段。此時堆內(nèi)冷卻劑依靠密度差進行流動，自然循環(huán)階段持續(xù)時間往往較短。

3)環(huán)路水封存在及清除階段

當一回路壓力降至飽和壓力時，一回路內(nèi)冷卻劑開始沸騰，產(chǎn)生大量蒸汽。由于破口尺寸較小，蒸汽不能及時從破口排出，導(dǎo)致積聚在上腔室和蒸汽發(fā)生器(SG)U型管內(nèi)，蒸汽的積聚又阻礙了冷卻劑的流動，進一步導(dǎo)致?lián)Q熱惡化，這個現(xiàn)象被稱為環(huán)路水封。環(huán)路水封的出現(xiàn)會導(dǎo)致堆芯液位快速下降。隨著一回路冷卻劑從破口的不斷流出，環(huán)路水封被逐漸清除，一回路壓力恢復(fù)下降，且堆芯液位開始被安注系統(tǒng)恢復(fù)。

4)長期冷卻階段

環(huán)路水封清除后，堆芯液位開始恢復(fù)，直至完全淹沒，過冷冷卻劑能持續(xù)穩(wěn)定地帶走堆內(nèi)的衰變熱。

對于破口尺寸較小的SBLOCA工況(如0.2%破口工況)而言，其堆芯活性區(qū)MCL出現(xiàn)在噴放階段，主要原因是閃蒸現(xiàn)象導(dǎo)致的堆芯局部沸騰及破口的冷卻劑流出。而對于破口尺寸較大的SBLOCA工況而言，其堆芯活性區(qū)MCL出現(xiàn)在環(huán)路水封存在及清除階段，主要原因是堆內(nèi)傳熱惡化導(dǎo)致的堆芯強烈沸騰及破口的冷卻劑流出。由于兩種情況下堆芯活性區(qū)MCL出現(xiàn)的原理不同，后續(xù)敏感性分析應(yīng)分為兩個階段獨立開展，以評估各輸入?yún)?shù)在噴放階段和環(huán)路水封存在及清除階段對堆芯活性區(qū)MCL的影響。

3 敏感性分析

由于缺乏PWR SBLOCA相關(guān)PIRT，因此擬參考類似研究中建立的PIRT或使用的重要參數(shù)，并結(jié)合目標電廠補充相關(guān)參數(shù)。諸多針對SBLOCA的PIRT中，西屋公司針對AP600開發(fā)的PIRT包含了典型能動式壓水堆和非能動式壓水堆中的主要現(xiàn)象，具有較高的可信度，因此本文主要參考的PIRT為AP600核電廠的SBLOCA PIRT[12]，但不包含非能動設(shè)備相關(guān)的現(xiàn)象或參數(shù)。此外，還參考了文獻[13-15]中對SBLOCA工況開展BEPU分析時使用的重要參數(shù)。基于參考類似PIRT共識別了34個輸入?yún)?shù)，此外，通過分析基準工況計算中出現(xiàn)的重要現(xiàn)象，補充了可能存在重要影響的模型或參數(shù)共20個輸入，因此在敏感性分析計算中共考慮了54個輸入?yún)?shù)。根據(jù)提出的優(yōu)化矩獨立全局敏感性分析方法，54個輸入?yún)?shù)執(zhí)行全局敏感性分析需進行217次程序計算。因此，使用所有54個不確定性輸入?yún)?shù)執(zhí)行全局敏感性分析，于5個破口尺寸工況各執(zhí)行217次程序計算，并進一步計算54個輸入?yún)?shù)在噴放階段和環(huán)路水封存在及清除階段對堆芯活性區(qū)MCL的矩獨立敏感性度量。計算結(jié)果如圖2所示，由于0.2%破口尺寸下沒有環(huán)路水封存在及清除階段，因此沒有相應(yīng)的敏感性分析結(jié)果。

圖2 0.2%、0.4%、0.6%、0.8%、1.0%破口尺寸工況敏感性分析結(jié)果

根據(jù)Pareto法則[16]，復(fù)雜模型可被小部分重要的參數(shù)和大部分影響較小的參數(shù)表征。圖2示出了不同破口尺寸和不同事故階段中各輸入?yún)?shù)對堆芯活性區(qū)MCL的影響，表明大部分影響集中于少數(shù)重要輸入?yún)?shù)，該結(jié)果符合Pareto法則?？紤]到缺乏SBLOCA的PIRT，因此可基于敏感性分析結(jié)果建立參數(shù)重要度排序分組表，可為PWR SBLOCA工況BEPU分析提供參考和指導(dǎo)。需要說明的是，由于大部分影響集中于少數(shù)重要輸入?yún)?shù)，若按照敏感性度量進行參數(shù)排序和重要度分組，會出現(xiàn)重要度排序為高的參數(shù)數(shù)量較少的情況。而且不同破口尺寸工況計算得到的各參數(shù)的敏感性度量無法直接比較，因此可將參數(shù)的重要度排序轉(zhuǎn)換為Savage分數(shù)，然后將所有參數(shù)進行重要度定性均勻分組，分別為高(H)、中(M)、低(L)和無影響(N)。Savage分數(shù)的定義為：

(5)

式中：N為所有輸入?yún)?shù)的數(shù)目；Ri為第i個輸入?yún)?shù)的秩排序。

根據(jù)計算結(jié)果，噴放階段有16個參數(shù)沒有影響，其重要度排序為54，分組為N。剩余38個參數(shù)進行均勻分組，即排序1～13的參數(shù)分組為H，排序14～26的參數(shù)分組為M，排序27～38的參數(shù)分組為L。同理，在環(huán)路水封存在及清除階段有7個參數(shù)對MCL沒有影響，其分組為N，并對剩余47個參數(shù)進行均勻分組。根據(jù)最終對結(jié)果的處理，可得到PWR在SBLOCA工況下的參數(shù)重要度排序表(表1)。由表1可知，噴放階段對MCL影響較大的參數(shù)主要為初始功率、初始一回路壓力、臨界流模型以及一回路中的流動阻力。而在環(huán)路水封存在及清除階段，除破口冷卻劑損失外，堆芯液位下降還與堆內(nèi)沸騰有關(guān)，因此界面阻力模型、各種換熱模型以及輔助給水系統(tǒng)延遲也對MCL有較大的影響。對比發(fā)現(xiàn)，提出方法計算得到的敏感性分析結(jié)果與參考PIRT較為相似，此外亦與現(xiàn)象機理具有高度一致性，因此可說明該敏感性分析結(jié)果的可信度。需要補充說明的是，該計算得到的PIRT能適用于典型能動式PWR SBLOCA工況的現(xiàn)象識別，亦可用作非能動式PWR SBLOCA工況的現(xiàn)象識別參考。

表1 PWR SBLOCA工況參數(shù)重要度排序分組表

4 結(jié)論

本文研究了在缺乏相應(yīng)PIRT的情況下確定各輸入?yún)?shù)的重要度，選取的研究工況是典型三回路PWR在冷管段發(fā)生破口尺寸分別為0.2%、0.4%、0.6%、0.8%和1.0%下的SBLOCA。首先基于基準計算結(jié)果，結(jié)合已有的SBLOCA PIRT，篩選了可能對FOM，即堆芯活性區(qū)MCL，具有影響的54個不確定性輸入?yún)?shù)。然后，使用一種優(yōu)化矩獨立全局敏感性分析方法，計算得到了不同破口尺寸下各輸入?yún)?shù)對FOM的敏感性度量和重要度排序。將參數(shù)的重要度排序轉(zhuǎn)換為Savage分數(shù)，按照Savage分數(shù)定性地將所有輸入?yún)?shù)進行重要度分組，從而得到了PWR SBLOCA的參數(shù)重要度排序表。本文提供了一種在缺乏PIRT時識別重要參數(shù)的方法。計算結(jié)果表明，該文使用的優(yōu)化矩獨立方法能以較小的計算成本完成全局敏感性分析，并能準確地識別對FOM沒有影響的輸入?yún)?shù)?；诰鬲毩⒚舾行苑治龅挠嬎憬Y(jié)果，建立了PWR在SBLOCA下的參數(shù)重要度排序表，該表能夠在SBLOCA相關(guān)分析缺乏PIRT的情況下提供一定的指導(dǎo)和參考價值。敏感性分析結(jié)果表明，噴放階段對MCL影響較大的參數(shù)主要為初始功率、初始一回路壓力、臨界流模型以及一回路中的流動阻力。而在環(huán)路水封存在及清除階段，除破口冷卻劑損失外，堆芯液位下降還與堆內(nèi)沸騰有關(guān)，因此界面阻力模型、各種換熱模型以及輔助給水系統(tǒng)延遲也對MCL有較大的影響。