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

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

(1.中國(guó)核動(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)用于華龍一號(hào)和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)有效性的評(píng)價(jià)，該模型基于熔池分層結(jié)構(gòu)及材料成分，采用衰變熱內(nèi)熱源驅(qū)動(dòng)的自然對(duì)流計(jì)算得到壓力容器下封頭的表面熱流密度、壓力容器壁厚、氧化殼厚度等參數(shù)，實(shí)現(xiàn)對(duì)壓力容器下封頭完整性的評(píng)價(jià)，并為嚴(yán)重事故緩解策略的制定提供參考。由于熔池形成過程涉及大量的質(zhì)量交換、能量傳遞，使下封頭熔池模型具有關(guān)系式復(fù)雜、輸入?yún)?shù)多且具有很大不確定性的特點(diǎn)，分析輸入因素對(duì)下封頭有效性的影響程度有助于下封頭熔池模型嚴(yán)重事故緩解策略制定的優(yōu)化。

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

1 理論模型

1.1 方差分解法

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

將其稱為對(duì)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ù)的不確定性對(duì)輸出參數(shù)方差的貢獻(xiàn)程度，全部敏感度體現(xiàn)了單個(gè)輸入?yún)?shù)不確定度以及該參數(shù)與其他參數(shù)的相互作用對(duì)輸出參數(shù)方差的貢獻(xiàn)程度。

近年來Saltelli等[8]對(duì)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)化工作，中國(guó)核動(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]已開展了對(duì)CISER軟件在CIS有效性評(píng)價(jià)方面的驗(yàn)證，本文著重對(duì)敏感性分析程序SALib進(jìn)行驗(yàn)證。選擇線性模型和非線性模型分別對(duì)SALib進(jìn)行線性模型和非線性模型敏感性分析能力的驗(yàn)證。隨后基于Khatib-Rahbar等[15]研究中AP1000的熔池初始參數(shù)對(duì)下封頭熔池模型的關(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é)果對(duì)比列于表2。

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

(22)

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

線性模型與理論解的相對(duì)偏差在1%范圍內(nèi)，而非線性模型與理論解的相對(duì)偏差在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等對(duì)壓力容器外部成功再淹沒條件下的事故進(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ù)

壓力容器下封頭外表面的熱流密度與測(cè)量結(jié)果的比值r=qw/qCHF、氧化層硬殼厚度δcr、壓力容器壁面厚度δw常被作為CIS有效性評(píng)價(jià)的指標(biāo)。本文選擇壓力容器下封頭外表面的熱流密度與試驗(yàn)測(cè)量結(jié)果[1,4,12-13]的比值的平均值(qw/qCHF)avg、最大值(qw/qCHF)max和壓力容器壁面厚度的最小值(δw)min及氧化層硬殼厚度的平均值(δcr)avg作為下封頭模型評(píng)價(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)的變化趨勢(shì)。

從圖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ù)的變化趨勢(shì)Fig.4 (qw/qCHF)avg vs. input parameter

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

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

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

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

圖7示出了(δcr)avg的輸入?yún)?shù)敏感性系數(shù)，圖8示出了(δcr)avg隨二氧化鈾質(zhì)量等輸入?yún)?shù)的變化趨勢(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ù)的變化趨勢(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)、壁面厚度對(duì)(δcr)avg有較小的影響(敏感性系數(shù)約0.04～0.05)，下封頭半徑和剩余衰變熱對(duì)(δ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ù)的變化趨勢(shì)。

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

圖8 (δcr)avg隨輸入?yún)?shù)的變化趨勢(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評(píng)價(jià)軟件CISER對(duì)下封頭熔池模型開展了關(guān)鍵結(jié)果參數(shù)與典型輸入?yún)?shù)之間的敏感性分析，配合下封頭熔池模型關(guān)鍵結(jié)果與輸入?yún)?shù)的變化趨勢(shì)圖，得到如下結(jié)論：

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

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

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

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

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