王寶生，唐秀歡，朱 磊，包利紅

(西北核技術(shù)研究所，陜西 西安 710024)

西安脈沖堆(XAPR)作為池式研究堆，具有停堆后的非能動安全性，即當反應(yīng)堆發(fā)生事故時，依靠堆池水的自然循環(huán)冷卻將堆芯余熱非能動導(dǎo)出。然而這樣的B型非能動系統(tǒng)[1]在緩解嚴重事故時易受到不確定性的影響，從而引起系統(tǒng)的功能失效[2]。因此，有必要對XAPR堆芯自然循環(huán)冷卻能力的功能可靠性進行研究，發(fā)展一種評價可靠性的方法，為XAPR PSA提供重要的信息。

近年來，Marques等[3]和Jafari等[4]分別提出了功能可靠性評估方法RMPS和REPAS，并用于兩相流自然循環(huán)系統(tǒng)。謝國鋒等[5-6]分別采用響應(yīng)面(RS)法和蒙特卡羅抽樣(MCS)法計算了高溫氣冷堆余熱排出系統(tǒng)功能失效概率。夏少雄等[7]和潘曉磊等[8]分別采用神經(jīng)網(wǎng)絡(luò)法和響應(yīng)面擬合法計算了中國鉛基研究堆非能動系統(tǒng)的可靠性。研究表明，功能可靠性評估可采用MCS、減小方差的MCS及RS法。然而上述方法對于小功能失效概率需進行大量的抽樣，計算效率不高[9]。

本文結(jié)合中破口失水事故(MLOCA)對XAPR堆芯自然循環(huán)冷卻能力進行可靠性評價。通過改進響應(yīng)面(IRS)法獲得系統(tǒng)輸出與輸入變量的近似解析表達式，利用重要抽樣子集模擬(ISSS)法獲得系統(tǒng)輸出響應(yīng)的累積分布函數(shù)和置信區(qū)間，進而對系統(tǒng)進行功能可靠性評估與分析。

1 功能可靠性評估數(shù)學模型

假設(shè)影響系統(tǒng)功能的隨機輸入向量為X=(x1,x2,…,xn)，系統(tǒng)的功能失效概率Pf表示為：

(1)

式中：n為不確定參數(shù)個數(shù)；g(X)為功能函數(shù)，g(X)=A-y(X)，A為失效準則，y(X)表示熱工水力學程序，g(X)<0表示系統(tǒng)處于安全狀態(tài)，g(X)>0表示系統(tǒng)處于失效狀態(tài)，g(X)=0表示系統(tǒng)處于臨界狀態(tài)；fX(X)為變量X的聯(lián)合概率密度函數(shù)。

為便于計算，式(1)可改寫為如下形式：

(2)

式中：I[g(X)]為指示性函數(shù)，若g(X)<0，則I[g(X)]=1，若g(X)≥0，則I[g(X)]=0；Ω為整個積分區(qū)域。

對于復(fù)雜的自然循環(huán)系統(tǒng)，難以用快速積分的方法求解式(2)。因為g(X)是一個復(fù)雜的熱工水力學程序，其功能函數(shù)都是模擬程序，需用數(shù)值方法求解。由于功能失效概率較小，采用MCS計算需進行多次抽樣，且每次抽樣都需運行熱工水力學程序，計算效率較低。為此，本文提出提高計算效率的兩種方法：一是建立響應(yīng)面替代回歸模型，以加快運算速度；二是使用高效蒙特卡羅方法即ISSS法以較少的運算次數(shù)獲得更加穩(wěn)健的估計結(jié)果。

2 功能可靠性評估分析方法

2.1 IRS法

XAPR堆芯自然循環(huán)過程由熱工水力學程序模擬，其功能函數(shù)均為隱式，可以構(gòu)建失效響應(yīng)面的方法來近似表達隱式功能函數(shù)[5]。若用響應(yīng)面法直接求解，需運行熱工水力學程序k(2n+1)次，當?shù)螖?shù)k較大時，計算耗時[8]。為減少計算次數(shù)，本文提出IRS法構(gòu)造失效響應(yīng)面功能函數(shù)，并用該函數(shù)代替熱工水力學程序。IRS法的流程如圖1所示，具體步驟如下。

1) 獲得樣本點

本文樣本點分為兩部分：一是通過直接MCS在不確定性輸入?yún)?shù)區(qū)間內(nèi)隨機抽取M組輸入向量，代入熱工水力學程序計算相應(yīng)輸出向量；二是在失效域附近進行重要密度函數(shù)抽取N組條件樣本點，并由熱工水力學程序計算相應(yīng)輸出量。

2) 估計響應(yīng)面功能函數(shù)系數(shù)

采用bootstrap方法從步驟1獲得的樣本點中產(chǎn)生t組bootstrap自舉樣本作為輸入樣本，并由熱工水力學程序計算相應(yīng)輸出量。利用最小二乘法計算出每組響應(yīng)面功能函數(shù)系數(shù)，并計算均方根誤差，以均方根誤差最小為目標函數(shù)，確定響應(yīng)面功能函數(shù)系數(shù)。

圖1 IRS法流程Fig.1 Flowchart of IRS method

3) 功能失效概率計算

從不確定性輸入?yún)?shù)的分布中進行ISSS抽樣，代入步驟2的響應(yīng)面功能函數(shù)中，獲得系統(tǒng)輸出響應(yīng)的累積分布函數(shù)和置信區(qū)間，進行統(tǒng)計計算得到系統(tǒng)的功能失效概率。

2.2 ISSS法

通過引入合理的中間失效事件，將小失效概率表達為一系列較大的條件概率的乘積形式[10-11]，而較大條件失效概率通過重要抽樣來高效估計。給定功能函數(shù)g(X)的失效域F，引入一系列中間失效事件臨界值0

=…

(3)

式(3)中P(F1)可由MCS計算，條件失效概率P(Fi|Fi-1)可由重要抽樣密度函數(shù)通過自動分層計算。ISSS法抽樣流程如圖2所示，具體過程如下。

a——直接蒙特卡羅模擬；b——A1的選取；c——重要抽樣模擬；d——A2的選取圖2 ISSS法抽樣流程Fig.2 Illustration of importance sampling subset simulation procedure

i=1,2,…,m-1

(4)

5) 重復(fù)步驟3和4，直到滿足條件Am≥A，此時自動分層結(jié)束。

6) 計算功能失效概率估算值為：

(5)

3 XAPR堆芯自然循環(huán)功能可靠性分析

3.1 系統(tǒng)安全功能

XAPR采用池式反應(yīng)堆結(jié)構(gòu)，堆芯采用自然循環(huán)冷卻。當反應(yīng)堆停堆后，堆芯余熱可依靠池水的自然循環(huán)冷卻，池水與環(huán)境間的熱量交換即可將堆芯余熱非能動導(dǎo)出。當發(fā)生失水事故時，事故進程可劃分為3個階段：1) 水冷自然循環(huán)冷卻階段；2) 堆芯半裸露階段；3) 空氣冷卻自然循環(huán)階段。這樣即可在失水事故時建立堆芯長期自然循環(huán)冷卻，將堆芯余熱導(dǎo)出，保證燃料元件的完整性。

3.2 系統(tǒng)建模與驗證

采用RELAP5模擬失水事故堆芯自然循環(huán)冷卻過程。XAPR熱工水力建模計算節(jié)點示于圖3。圖3中控制體101～107為堆池，120～124為堆池混凝土壁外空氣邊界，150～163為一回路管道，170～172為簡化二回路邊界，109為堆池頂部空氣邊界。為驗證模型的準確性，采用穩(wěn)態(tài)運行工況作為邊界條件進行計算，得到堆芯自然循環(huán)流量為12.21 kg/s，該流量與實際自然循環(huán)流量12.13 kg/s基本一致。

3.3 系統(tǒng)不確定性參數(shù)

XAPR自然循環(huán)運行時涉及到2種不確定性：第1種是“偶然性”所導(dǎo)致的不確定性，它與模型的幾何性質(zhì)有關(guān)；第2種是“認知局限”所導(dǎo)致的不確定性，它與運行和實驗數(shù)據(jù)的缺乏而導(dǎo)致的對有關(guān)現(xiàn)象認知的局限性有關(guān)[10]。本文主要研究因“認知局限”所導(dǎo)致的不確定性。

圖3 RELAP5計算節(jié)點圖Fig.3 Nodalization scheme of RELAP5 calculation

通過層次分析法選出14個不確定性輸入?yún)?shù)。在目前處理中，當可得到的數(shù)據(jù)有限時，各輸入?yún)?shù)的概率分布及取值區(qū)間主要以工程設(shè)計標準為基礎(chǔ)，并結(jié)合專家評價的主觀方法得出。表1列出輸入?yún)?shù)的概率分布及特征參數(shù)。這些參數(shù)的取值及分布考慮了所有可能的事故工況，是一種較為保守的方法。

3.4 系統(tǒng)失效準則

定義功能準則為：當包殼溫度低于或等于500 ℃時，燃料芯體最高溫度應(yīng)低于1 150 ℃；包殼溫度高于500 ℃時，燃料芯體最高溫度應(yīng)低于970 ℃。從概率安全角度出發(fā)，假設(shè)燃料芯體最高溫度超過安全限值970 ℃則認為功能失效。設(shè)向量X為系統(tǒng)不確定性參數(shù)(如表1所列)，To,max(X)為燃料芯體最高溫度，則功能函數(shù)g(X)可表示為g(X)=970-To,max(X)，當g(X)<0時，即認為失效。

4 功能可靠性評價與結(jié)果分析

4.1 燃料芯體最高溫度置信水平

用于擬合響應(yīng)面功能函數(shù)的樣本點包括：第1部分采用MCS抽取M=120組隨機樣本點；第2部分采用重要密度函數(shù)抽取N=24組樣本點；第3部分在接近失效域附件采用bootstrap方法抽取12組樣本點。共156組樣本點，其RELAP5運行計算耗時為159.8 h，其中tM=120.1 h、tN=27.3 h及ttb=12.4 h。利用156組樣本點求解響應(yīng)面功能函數(shù)的系數(shù)，結(jié)果如下：

1.72Tg+0.18p-31.19ξ1+1.34×102ξ2-

6.3ξ3+6.21ξ4+0.09dx+13.48kin+

12.98kout-4.79td+8.54tm+12.89tp-

(6)

圖4示出燃料芯體最高溫度對比與殘差散點圖。由圖4a可看出，響應(yīng)面函數(shù)的擬合值與REALP5的輸出值吻合很好。由圖4b可看出，各工況的殘差均在±6 ℃以內(nèi)。因此可將響應(yīng)面函數(shù)作為RELAP5替代模型進行功能失效概率計算。

為研究IRS在抽樣效率上的提升，計算了156組抽樣在RS和IRS抽樣中的比例。在利用156組抽樣得到的輸出結(jié)果中，RS大約有2組數(shù)據(jù)導(dǎo)致To,max大于970 ℃，所占比例約為1.28%；IRS得到15組數(shù)據(jù)導(dǎo)致To,max大于970 ℃，所占比例約為9.61%，失效事件出現(xiàn)的概率提高了將近8倍，抽樣效率明顯提升。圖5示出RS和IRS對20組樣本點的To,max的預(yù)測結(jié)果。

表1 輸入?yún)?shù)的概率分布及特征參數(shù)Table 1 Probability distribution and characteristic parameter of input parameter

圖4 燃料芯體最高溫度對比(a)與殘差散點圖(b)Fig.4 Comparison of fuel peak temperature (a) and residual scatter plot (b)

圖5 20組樣本點的To,max預(yù)測結(jié)果Fig.5 Prediction result of To,max of 20 samples

圖6示出α=0.9、0.95、0.99及0.999等4種不同置信水平下燃料芯體最高溫度的累積分布函數(shù)(CDF)。由圖6可看出，不同置信水平下，ISSS均有相當數(shù)量的樣本點落在計算區(qū)域內(nèi)。但MCS抽取400個樣本點時，在高置信水平下，很難得到精確值。這是由于當抽樣數(shù)為400時，只有平均400×(1-0.99)=4個樣本點落在0.99～1區(qū)間內(nèi)，如圖6c所示；而只有400×(1-0.999)=0.4個點落入0.999～1區(qū)域，如圖6d所示。然而，對于ISSS能保證有100、10個樣本點落入相關(guān)區(qū)域，這就給可靠性分析提供了足夠的信息。注意到MCS抽樣10 000個樣本點，分別也有100、10個樣本點落入目標區(qū)域，如圖6c、d所示，但抽樣次數(shù)較多，計算效率較低。

4種不同置信水平下，燃料芯體最高溫度均值與標準差列于表2。由表2可看出，ISSS的結(jié)果與MCS抽樣10 000次的結(jié)果基本一致，即便是在高置信水平下，兩者的標準差也相差不大，但ISSS的抽樣數(shù)明顯低于MCS的抽樣數(shù)，僅約為后者的1/25。

α：a——0.9；b——0.95；c——0.99；d——0.999圖6 不同置信水平下燃料芯體最高溫度的累積分布函數(shù)Fig.6 CDF for fuel peak temperature at different confidence levels

表2 不同置信水平下燃料芯體最高溫度均值與標準差Table 2 Mean value and standard deviation of fuel peak temperature at different confidence levels

4.2 功能失效概率計算及比對

不同層次下燃料芯體最高溫度的自動分層結(jié)果如圖7所示。由圖7可看出，各中間失效事件F1、F2、F3、F4的臨界值A(chǔ)1=843.52 ℃、A2=910.25 ℃、A3=948.79 ℃及A4=999.08 ℃，由于A4=999.08 ℃>970 ℃(圖7d)，自動分層結(jié)束。根據(jù)式(5)計算自然循環(huán)冷卻能力的功能失效概率為Pf=3.796×10-3。為便于各種計算方法的比較，還同時采用了MCS、IRS、重要抽樣(IS)法進行計算，計算結(jié)果列于表3。

由表3可知，ISSS+IRS抽樣400次，即可達到MCS+RELAP5抽樣104時的計算精度。MCS+IRS抽樣400次時，其計算結(jié)果為1.875×10-3，導(dǎo)致這樣的結(jié)果是因為當失效概率為10-3數(shù)量級時，只有約400×10-3=0.4個樣本點落入失效域內(nèi)。IRS實質(zhì)是曲線擬合近似方法，其計算精度依賴于功能函數(shù)的線性程度，XAPR自然循環(huán)的功能函數(shù)非線性程度較高，所以結(jié)果欠佳。MCS+RELAP5抽樣效率低、計算費時，在抽樣次數(shù)為104的情況下計算耗時約為ISSS+IRS的6.25倍，因而不適用。在相同抽樣次數(shù)的情況下，IS+IRS計算耗時與ISSS+IRS接近，但其計算精確度略差。

4.3 參數(shù)靈敏度分析

靈敏度分析可反映輸入?yún)?shù)的不確定性對系統(tǒng)功能可靠性的影響，從而為提高系統(tǒng)功能可靠性提供指導(dǎo)性建議，有效地減小不確定性。由于非能動物理過程為隱式非線性非單調(diào)關(guān)系，傳統(tǒng)表征參數(shù)靈敏度的相關(guān)系數(shù)(CCS)、標準回歸系數(shù)(SRCS)、秩相關(guān)系數(shù)(RCCS)及標準秩回歸系數(shù)(SRRCS)等具有局限性[12]。本文將基本變量xi(i=1,2,…,n)均值μxi的變化引起的功能失效概率變化的比率表征為靈敏度，在數(shù)學上表達為由失效概率Pf對基本變量均值μxi的偏導(dǎo)數(shù)，即?Pf/?μxi。為消除量綱影響，引入了無量綱靈敏度系數(shù)Sμxi=(?Pf/?μxi)×(σxi/Pf)[13]，其中σxi為輸入變量標準差。靈敏度分析結(jié)果如圖8所示。

a——層次1；b——層次2；c——層次3；d——層次4圖7 不同層次下的燃料芯體最高溫度Fig.7 Fuel peak temperature at different conditional levels

表3 不同計算方法計算的功能失效概率Table 3 Functional failure probability for different calculation methods

由圖8可見，空氣入口溫度Tg、衰變熱功率相對偏差Q、流道進口阻力系數(shù)kin、堆池水入口溫度Tw及排水泵開啟延遲時間tp對系統(tǒng)輸出有較大影響，其中Tg、Q及kin最為敏感。因為這些輸入?yún)?shù)直接影響事故發(fā)生后堆芯自然循環(huán)的輸熱能力和建立時間，對堆芯余熱排出功能可靠性影響較大，因此，減小上述輸入?yún)?shù)的不確定性可有效降低堆芯自然循環(huán)的功能失效概率，提高系統(tǒng)功能的可靠性。

圖8 輸入?yún)?shù)的無量綱靈敏度系數(shù)Fig.8 Normalized sensitivity coefficient of input parameter

5 結(jié)論

1) 本文方法將IRS和ISSS相結(jié)合，基于失效響應(yīng)面函數(shù)將功能失效概率表述為一系列較大條件失效概率的乘積，通過子集模擬使得失效域內(nèi)樣本計算效率大為提高。相比于其他方法，該方法具有較高的抽樣效率，同時又能保證較好的計算精度，且對于非線性程度較高的自然循環(huán)系統(tǒng)有很強的適應(yīng)性。

2) XAPR中破口失水事故時，由“認知局限”不確定性導(dǎo)致的堆芯自然循環(huán)失效總有非零的發(fā)生概率，其失效概率為3.796×10-3。為此，在涉及非能動自然循環(huán)的可靠性評估中，功能失效是重要的，其失效概率應(yīng)基于熱工水力計算和可靠性方法進行全面評價。

3) 減小Tg、Q、kin、Tw及tp等5個關(guān)鍵參數(shù)的不確定性，可更有效地降低XAPR堆芯自然循環(huán)的失效概率，提高XAPR系統(tǒng)的功能可靠性。