基于方差分解的堆芯下封頭熔池模型敏感性分析

李治剛，安 萍，明平洲，潘俊杰，蘆 韡,*，余紅星

(1.中國核動(dòng)力研究設(shè)計(jì)院，四川 成都 610041；2.核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041)

在發(fā)生堆芯熔化的嚴(yán)重事故中，堆芯熔融物可能會(huì)遷移至壓力容器下封頭，嚴(yán)重威脅壓力容器的完整性[1-2]。熔融物堆內(nèi)滯留(IVR)是一項(xiàng)重要的嚴(yán)重事故緩解措施，已應(yīng)用于華龍一號和AP1000等[3]三代堆中。熔池熔融物的分布及流動(dòng)換熱特性決定了壓力容器壁面熱流密度的分布，直接影響IVR的有效性。Theofanous等[4]在AP600的IVR分析報(bào)告中提出了氧化物層和金屬層的兩層熔池模型。近年來，Parker等[5]發(fā)現(xiàn)二氧化鈾、二氧化鋯與金屬鋯在高溫熔池內(nèi)會(huì)發(fā)生復(fù)雜的化學(xué)反應(yīng)，形成以鈾為主要成分的重金屬層，提出了三層熔池模型。

采用集總參數(shù)法的下封頭熔池模型在工程上廣泛用于堆腔注水冷卻系統(tǒng)(cavity injection and cooling system, CIS)有效性的評價(jià)，該模型基于熔池分層結(jié)構(gòu)及材料成分，采用衰變熱內(nèi)熱源驅(qū)動(dòng)的自然對流計(jì)算得到壓力容器下封頭的表面熱流密度、壓力容器壁厚、氧化殼厚度等參數(shù)，實(shí)現(xiàn)對壓力容器下封頭完整性的評價(jià)，并為嚴(yán)重事故緩解策略的制定提供參考。由于熔池形成過程涉及大量的質(zhì)量交換、能量傳遞，使下封頭熔池模型具有關(guān)系式復(fù)雜、輸入?yún)?shù)多且具有很大不確定性的特點(diǎn)，分析輸入因素對下封頭有效性的影響程度有助于下封頭熔池模型嚴(yán)重事故緩解策略制定的優(yōu)化。

傳統(tǒng)的局部敏感性分析方法在應(yīng)用于復(fù)雜模型的敏感性分析時(shí)具有效率低、計(jì)算量大的缺點(diǎn)，本文采用中國核動(dòng)力研究設(shè)計(jì)院自主研發(fā)的基于方差分解的全局敏感性分析工具SALib(sensitivity analysis library)和下封頭冷卻劑注入系統(tǒng)(coolant inject system, CIS)有效性評價(jià)程序CISER對下封頭熔池模型進(jìn)行輸入?yún)?shù)敏感性分析，系統(tǒng)介紹方差分解法和下封頭熔池模型的基本理論，并對敏感性分析工具SALib進(jìn)行初步驗(yàn)證，針對下封頭熔池模型開展AP1000重要結(jié)果的敏感性分析。

1 理論模型

1.1 方差分解法

敏感性分析是一種定量描述模型輸入?yún)?shù)對輸出結(jié)果的重要性程度的方法，包括局部敏感性分析方法和全局敏感性分析方法。局部敏感性分析只檢驗(yàn)單個(gè)參數(shù)對模型的影響程度；而全局敏感性分析[6]可同時(shí)考慮多個(gè)參數(shù)對模型輸出的影響，并分析各參數(shù)之間的相互作用對模型輸出的影響，廣大研究者對該方法開展了深入的研究和推廣，提出了如方差分解法、神經(jīng)網(wǎng)絡(luò)法、PaD2方法、FAST方法等。

方差分解法是一種典型的全局敏感性分析方法，由數(shù)學(xué)家Sobol[7]提出，其核心是把模型分解為單個(gè)屬性及屬性之間相互組合的函數(shù)。

假設(shè)模型為y=f(x)，其中x=(x1,x2,…,xn)，xi服從[0，1]均勻分布，且f2(x)可積，模型可分解為：

(1)

其中，f0為常數(shù)，且每個(gè)分解項(xiàng)fi1…is(xi1,…,xis)滿足：

1≤s≤n

(2)

將其稱為對f(X)的分解，并且分解是唯一的。

則模型總的方差也可分解為單個(gè)參數(shù)和每個(gè)參數(shù)相互組合的影響：

(3)

將式(3)歸一化，并設(shè)：

Si1,…,in=Di1,…,in/D

(4)

(5)

(6)

顯然D和Di1…is分別是f(x)和fi1…is(xi1,…,xis)的方差。

根據(jù)式(1)，可獲得模型單個(gè)參數(shù)及參數(shù)之間相互作用的敏感度S。

(7)

式中：Si稱為1次敏感度；Sij稱為2次敏感度，依次類推，S1,2,…,n為n次敏感度，共有2n-1項(xiàng)，敏感度即方差分解法計(jì)算的敏感性系數(shù)。參數(shù)xi總敏感度定義為：

STj=∑Si

(8)

它表示所有包含參數(shù)xi的敏感度。1次敏感度Si體現(xiàn)了單個(gè)輸入?yún)?shù)的不確定性對輸出參數(shù)方差的貢獻(xiàn)程度，全部敏感度體現(xiàn)了單個(gè)輸入?yún)?shù)不確定度以及該參數(shù)與其他參數(shù)的相互作用對輸出參數(shù)方差的貢獻(xiàn)程度。

近年來Saltelli等[8]對Sobol方法中敏感度的計(jì)算進(jìn)行了優(yōu)化，如式(9)和(10)所示，使得該方法計(jì)算效率高、易編程實(shí)現(xiàn)。

Si=Vxi(Ex～i(Y|xi))/V(Y)

(9)

STi=Ex～i(Vxi(Y|x～i))/V(Y)

(10)

具體的步驟如下：

1) 生成N×2D的樣本矩陣，該步驟可采用抽樣方法或Sobol sequence[9]生成；

2) 將矩陣的前D列設(shè)置為矩陣A，后D列設(shè)置為矩陣B；

3) 構(gòu)造N×D的矩陣ABi(i=1,2,…,D)，即用矩陣B中的第i列替換矩陣A的第i列；

計(jì)算1次敏感度Si和總體敏感度STi，Saltelli等[10-11]的計(jì)算公式列于表1。

表1 敏感度計(jì)算公式Table 1 Sensitivity calculation formula

基于方差分解敏感性分析方法的基本理論及Saltelli等的優(yōu)化工作，中國核動(dòng)力研究設(shè)計(jì)院開發(fā)了具有通用性的敏感性分析工具SALib，SALib的計(jì)算流程如圖1所示。

圖1 SALib程序計(jì)算流程圖Fig.1 Calculation flow chart of SALib code

1.2 下封頭熔池模型

3層熔池結(jié)構(gòu)[12-14]的傳熱模型如圖2所示。

圖2 熔池傳熱模型示意圖Fig.2 Molten pool heat transfer model

模型關(guān)系式如下(其氧化物層和輕金屬層的關(guān)系式適用于兩層熔池結(jié)構(gòu))。

中間氧化層：

Qo,vVo=qo,upAo,up+qo,dn(Ao,side+Ao,dn)

(11)

(12)

(13)

(14)

頂部輕金屬層：

Ql,vVl+ql,dnAl,dn=ql,upAl,up+ql,sideAl,side

(15)

(16)

(17)

(18)

底部重金屬層：

Qh,vVh+qh,upAh,up=qh,dnAh,dn

(19)

(20)

式中：下標(biāo)o、h、l分別代表氧化層、重金屬層、輕金屬層；Qo,v、Qh,v、Ql,v為熔池各層的體積功率密度；qo,up、qh,up、ql,up為熔池各層向上的熱流密度；qo,dn、qh,dn、ql,dn為熔池各層向下的熱流密度；Tl,bulk、Tl,m分別為輕金屬層主體溫度和熔化溫度；Ts,i、Ts,o分別為堆內(nèi)結(jié)構(gòu)件內(nèi)、外表面溫度；As、Al,up、Al,dn、Al,side為堆內(nèi)結(jié)構(gòu)件、輕金屬層的向上、向下和向側(cè)面的傳熱面積；Ah,up、Ah,dn分別為重金屬層向上、向下的傳熱面積；Ao,up、Ao,dn、Ao,side分別為氧化層向上、向下和向側(cè)面的傳熱面積；Nu為努塞爾數(shù)；εi、εs分別為輕金屬層和結(jié)構(gòu)件的輻射發(fā)射率；σ為斯特藩-玻爾茲曼常數(shù)；Vo、Vh、Vl為各層的總體積。

熔融池中總的衰變熱為停堆t小時(shí)后的衰變熱Pdecay-FP(t)。Pdecay-FP(t)為輸入值，氧化物層的體積功率密度為：

(21)

式中，fo-mUO2/UO2tot為氧化層中二氧化鈾占二氧化鈾總質(zhì)量的比例。

2 驗(yàn)證與應(yīng)用

向清安等[12]已開展了對CISER軟件在CIS有效性評價(jià)方面的驗(yàn)證，本文著重對敏感性分析程序SALib進(jìn)行驗(yàn)證。選擇線性模型和非線性模型分別對SALib進(jìn)行線性模型和非線性模型敏感性分析能力的驗(yàn)證。隨后基于Khatib-Rahbar等[15]研究中AP1000的熔池初始參數(shù)對下封頭熔池模型的關(guān)鍵結(jié)果參數(shù)開展敏感性分析，文中敏感性分析的樣本空間均由Sobol sequence[9]生成。

2.1 SALib敏感性分析能力驗(yàn)證

1) 線性模型

假設(shè)線性模型f(x1,x2)=ax1+bx2，xi在[0,1]內(nèi)均勻分布，該模型的理論解為D=a2/2+b2/2，D1=a2/2，D2=b2/2。a=1、b=1、樣本容量N=1 000時(shí)線性模型的敏感性分析結(jié)果對比列于表2。

表2 a=1、b=1、N=1 000時(shí)線性模型的敏感性分析結(jié)果對比Table 2 Sensitivity analysis result of linear model at a=1, b=1, and N=1 000

(22)

該模型的敏感性系數(shù)計(jì)算結(jié)果及對比情況列于表3、4。

線性模型與理論解的相對偏差在1%范圍內(nèi)，而非線性模型與理論解的相對偏差在5%范圍內(nèi)，采用Saltelli關(guān)系式計(jì)算的敏感性系數(shù)與Homma[16]結(jié)果相當(dāng)，具有較好的精度。

表3 a=7、b=0.1、N=1 000時(shí)非線性模型SiTable 3 Si result of nonlinear model at a=7, b=0.1, and N=1 000

表4 a=7、b=0.1、N=1 000時(shí)非線性模型STiTable 4 STi result of nonlinear model at a=7, b=0.1, and N=1 000

2.2 下封頭熔池模型的敏感性分析

1) 輸入?yún)?shù)

Khatib-Rahbar等[15]利用嚴(yán)重事故系統(tǒng)分析軟件MELCOR、MAAP等對壓力容器外部成功再淹沒條件下的事故進(jìn)程進(jìn)行了分析計(jì)算，得到了IVR分析的初始參數(shù)。IVR分析的初始參數(shù)包括初始質(zhì)量、材料物性、壓力容器的幾何結(jié)構(gòu)、衰變熱等參數(shù)，本文選擇了鋯金屬質(zhì)量、二氧化鈾質(zhì)量等9個(gè)代表性輸入?yún)?shù)，并假設(shè)其在初始參數(shù)10%范圍內(nèi)均勻分布(表5[1,15])。

2) 關(guān)鍵結(jié)果參數(shù)

壓力容器下封頭外表面的熱流密度與測量結(jié)果的比值r=qw/qCHF、氧化層硬殼厚度δcr、壓力容器壁面厚度δw常被作為CIS有效性評價(jià)的指標(biāo)。本文選擇壓力容器下封頭外表面的熱流密度與試驗(yàn)測量結(jié)果[1,4,12-13]的比值的平均值(qw/qCHF)avg、最大值(qw/qCHF)max和壓力容器壁面厚度的最小值(δw)min及氧化層硬殼厚度的平均值(δcr)avg作為下封頭模型評價(jià)的關(guān)鍵輸出結(jié)果。

表5 IVR分析初始參數(shù)Table 5 Initial parameter of IVR analysis

3) 敏感性系數(shù)

Saltelli關(guān)系式具有較好的計(jì)算效率和精度，本文將采用該關(guān)系式計(jì)算壓力容器下封頭壁面熱流密度比值等關(guān)鍵結(jié)果參數(shù)與輸入?yún)?shù)之間的敏感性系數(shù)，并分析不同樣本空間下敏感性系數(shù)的收斂性。

(1) 下封頭壁面熱流密度比值的平均值

圖3示出了(qw/qCHF)avg的輸入?yún)?shù)敏感性系數(shù)，圖4示出了(qw/qCHF)avg隨下封頭半徑(圖4a)和剩余衰變熱(圖4b)的變化趨勢。

從圖3可知，下封頭半徑和剩余衰變熱的敏感性系數(shù)分別為0.408和0.544，是影響下封頭壁面熱流密度比值平均值的關(guān)鍵輸入?yún)?shù)。從圖4可知，下封頭壁面熱流密度比值平均值隨下封頭半徑的增加而減小，隨剩余衰變熱的增加而增加，呈現(xiàn)典型的線性關(guān)系。

圖3 (qw/qCHF)avg的輸入?yún)?shù)敏感性系數(shù)Fig.3 Input parameter sensitivity coefficient of (qw/qCHF)avg

圖4 (qw/qCHF)avg隨輸入?yún)?shù)的變化趨勢Fig.4 (qw/qCHF)avg vs. input parameter

(2) 下封頭壁面熱流密度比值的最大值

圖5示出了(qw/qCHF)max的輸入?yún)?shù)敏感性系數(shù)，圖6示出了(qw/qCHF)max隨下封頭半徑(圖6a)和剩余衰變熱(圖6b)的變化趨勢。

從圖5、6可知，下封頭半徑和剩余衰變熱對(qw/qCHF)max的敏感性系數(shù)分別為0.559和0.426，(qw/qCHF)max隨下封頭半徑的增加而減小，隨剩余衰變熱的增加而增加。由圖3和圖5可知，剩余衰變熱和下封頭半徑對下封頭壁面熱流密度比值平均值和最大值的影響程度是不同的。

(3) 氧化層硬殼厚度平均值

圖7示出了(δcr)avg的輸入?yún)?shù)敏感性系數(shù)，圖8示出了(δcr)avg隨二氧化鈾質(zhì)量等輸入?yún)?shù)的變化趨勢。

圖5 (qw/qCHF)max的輸入?yún)?shù)敏感性系數(shù)Fig.5 Input parameter sensitivity coefficient of (qw/qCHF)max

圖6 (qw/qCHF)max隨輸入?yún)?shù)的變化趨勢Fig.6 (qw/qCHF)max vs. input parameter

圖7 (δcr)avg的輸入?yún)?shù)敏感性系數(shù)Fig.7 Input parameter sensitivity coefficient of (δcr)avg

從圖7可知，二氧化鈾質(zhì)量、壓力容器壁面熔點(diǎn)、壁面厚度對(δcr)avg有較小的影響(敏感性系數(shù)約0.04～0.05)，下封頭半徑和剩余衰變熱對(δcr)avg的敏感性系數(shù)分別為0.4和0.468。由圖8可知，(δcr)avg隨下封頭半徑的增加而增加，隨剩余衰變熱的增加而減小。

(4) 下封頭壁面厚度最小值

圖9示出了(δw)min的輸入?yún)?shù)敏感性系數(shù)，圖10示出了(δw)min隨壁面熔點(diǎn)等輸入?yún)?shù)的變化趨勢。

從圖9可知，壁面熔點(diǎn)和厚度、下封頭半徑、剩余衰變熱對(δw)min有明顯的影響，敏感性系數(shù)分別為0.275、0.336、0.179和0.162。從圖10可知，(δw)min隨壁面熔點(diǎn)的增加而增加、隨下封頭半徑的增加而增加、隨剩余衰變熱的增加而減小。

圖8 (δcr)avg隨輸入?yún)?shù)的變化趨勢Fig.8 (δcr)avg vs. input parameter

圖9 (δw)min的輸入?yún)?shù)敏感性系數(shù)Fig.9 Input parameter sensitivity coefficient of (δw)min

3 結(jié)論

本文采用基于方差分解的全局敏感性分析方法，采用自研的敏感性分析工具SALib和壓力容器下封頭IVR評價(jià)軟件CISER對下封頭熔池模型開展了關(guān)鍵結(jié)果參數(shù)與典型輸入?yún)?shù)之間的敏感性分析，配合下封頭熔池模型關(guān)鍵結(jié)果與輸入?yún)?shù)的變化趨勢圖，得到如下結(jié)論：

1) 參數(shù)樣本容量≥1 000時(shí)敏感性系數(shù)基本趨于一致；

2) 鋯質(zhì)量、二氧化鈾質(zhì)量、鋯氧化份額、結(jié)構(gòu)件的發(fā)射率的變化對下封頭壁面熱流密度比值的平均值等關(guān)鍵結(jié)果影響較?。?/p>

3) 剩余衰變熱和下封頭半徑對熔池下封頭模型中的結(jié)果參數(shù)影響很大，下封頭半徑的增加，有助于氧化層硬殼厚度的增加，降低壁面熱流密度比值;

4) 壓力容器壁面熔點(diǎn)的變化對下封頭壁面平均厚度有較大的影響，而對壁面熱流密度比值和氧化層硬殼厚度平均值的影響較小。

圖10 (δw)min隨輸入?yún)?shù)的變化趨勢Fig.10 (δw)min vs. input parameter