子通道程序?qū)SBT空泡分布實(shí)驗(yàn)計(jì)算的不確定性量化分析

張俊濤， 劉曉晶， 張滕飛， 柴 翔

(上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240)

1989年美國(guó)發(fā)布管理導(dǎo)則RG 1.157[1]，允許在應(yīng)急堆芯冷卻系統(tǒng)分析中使用最佳估算程序.出于安全考慮，最佳估算程序必須與不確定性分析相結(jié)合，以判斷計(jì)算與實(shí)際情況之間的差異.截至目前，國(guó)際上已有多種不同的不確定性分析方法，例如程序比例、適用性和不確定性分析方法(CSAU)[2]、德國(guó)核安全中心方法(GRS)[3]、基于精度外推的不確定性方法(UMAE)[4]、先進(jìn)統(tǒng)計(jì)學(xué)處理的不確定性方法(ASTRUM)[5]等.這些不確定性分析方法大體可以分為兩大類，一類是基于輸入?yún)?shù)不確定性的傳遞，另一類是基于輸出結(jié)果不確定性的外推[6].基于輸入?yún)?shù)不確定性傳遞的統(tǒng)計(jì)性方法的研究和應(yīng)用相對(duì)更成熟和廣泛，如在由經(jīng)濟(jì)合作與發(fā)展組織(OECD)發(fā)起的不確定性研究項(xiàng)目“最佳估算方法-不確定性和敏感性分析”階段V中，14個(gè)參與者中有12個(gè)采用輸入不確定性傳遞的方法[7].

基于輸入?yún)?shù)不確定性傳遞的方法需要知道輸入?yún)?shù)的不確定性.然而，大多數(shù)熱工水力程序都沒(méi)有提供關(guān)于輸入?yún)?shù)不確定性的足夠信息，尤其是對(duì)于某些無(wú)法通過(guò)實(shí)驗(yàn)直接測(cè)量的輸入模型參數(shù)，通常通過(guò)專家判斷來(lái)進(jìn)行量化.目前一些取代專家判斷的評(píng)估方法被陸續(xù)提出，如比薩大學(xué)開發(fā)的基于快速傅里葉變換的方法(FFTBM)[8]，法國(guó)開發(fā)的輸入?yún)?shù)經(jīng)驗(yàn)確定方法(DIPE)[8]和模型參數(shù)不確定性計(jì)算方法(CIRCé)[9]，韓國(guó)原子力研究院開發(fā)的通過(guò)數(shù)據(jù)同化進(jìn)行模型標(biāo)定方法(MCDA)[10].李冬[11]對(duì)CIRCé方法可評(píng)估參數(shù)有限和MCDA方法求解非線性問(wèn)題耗時(shí)過(guò)長(zhǎng)這兩個(gè)缺點(diǎn)進(jìn)行了改進(jìn)；Xiong等[12]開發(fā)了一種優(yōu)化的最佳估算加不確定性分析方法，該方法包含了一個(gè)結(jié)構(gòu)化的模型不確定性量化方法.

子通道程序已被廣泛應(yīng)用于堆芯熱工水力分析，相較于計(jì)算流體力學(xué)(CFD)，子通道分析消耗的計(jì)算資源更少，能夠快速提供分析所需的熱工參數(shù).由于子通道程序?qū)τ谖锢砟Ｐ痛嬖谝欢ê?jiǎn)化假設(shè)以及數(shù)值計(jì)算誤差等因素的影響，程序直接計(jì)算出的結(jié)果存在一定的不確定性，而目前針對(duì)子通道程序的不確定性分析還較少.本研究使用基于輸入?yún)?shù)不確定性傳遞的統(tǒng)計(jì)學(xué)方法，利用期望最大化的思想來(lái)反推輸入模型參數(shù)不確定性的分布，對(duì)子通道程序COBRA-IV 的計(jì)算進(jìn)行不確定性量化分析.

1 不確定性量化分析方法

不確定性分析方法包括對(duì)影響計(jì)算程序結(jié)果的所有不確定性因素的確認(rèn)和描述，并用統(tǒng)計(jì)學(xué)方法將這些不確定性因素整合得到輸出結(jié)果的總體不確定性，其一般過(guò)程如圖1所示.首先從眾多輸入?yún)?shù)中選取對(duì)計(jì)算結(jié)果影響較大的輸入?yún)?shù)；其次量化這些輸入模型參數(shù)的范圍或不確定性分布；再次通過(guò)蒙特卡洛抽樣或拉丁超立方抽樣對(duì)這些輸入?yún)?shù)進(jìn)行抽樣，并將這些抽樣的輸入?yún)?shù)輸入子通道程序進(jìn)行計(jì)算；最后對(duì)大量目標(biāo)參數(shù)值進(jìn)行統(tǒng)計(jì)，以得到輸出目標(biāo)參數(shù)的不確定性.其中，輸入?yún)?shù)不確定性的確定和對(duì)輸入?yún)?shù)進(jìn)行抽樣是關(guān)鍵步驟.

1.1 輸入模型參數(shù)不確定性量化

在分析過(guò)程中, 認(rèn)為輸出不確定性y是真實(shí)值yr與程序計(jì)算值yc之間的誤差.實(shí)際可以得到的數(shù)據(jù)是實(shí)驗(yàn)值ye與對(duì)應(yīng)的程序計(jì)算值yc，兩者之差可以表示為

yw=ye-yc=(ye-yr)+(yr-yc)=e+y

(1)

式中：e=(ej),j=1, 2，…,J為實(shí)驗(yàn)誤差，其中J為樣本點(diǎn)數(shù).

一般來(lái)說(shuō)，模型參數(shù)的不確定性分析通常要確定模型參數(shù)基準(zhǔn)值的無(wú)量綱乘數(shù)，記為k=(ki),i=1, 2, …,p，其中p為模型參數(shù)的個(gè)數(shù)；將待評(píng)估的模型參數(shù)不確定性記為x=(xi),i=1, 2, …,p.一般情況下，ki與xi存在以下關(guān)系：

(2)

因此可以采用ki=exi來(lái)表示模型參數(shù)系數(shù)和不確定性之間的轉(zhuǎn)化關(guān)系[8].

輸入模型參數(shù)不確定性x和輸出不確定性y之間的傳遞關(guān)系可簡(jiǎn)化為

y=h(x,u)

(3)

近似認(rèn)為x與y之間為線性關(guān)系，得到

yw=Sxx+Suu+e

(4)

式中：Sx和Su為樣本點(diǎn)對(duì)應(yīng)待求參數(shù)的敏感性矩陣，敏感性矩陣中的偏導(dǎo)數(shù)可以通過(guò)有限差分法求得[11].

當(dāng)存在數(shù)據(jù)缺失問(wèn)題時(shí)，期望最大化(EM)算法是極大似然估計(jì)的一種常用迭代算法.但是實(shí)際計(jì)算過(guò)程中EM算法一般迭代次數(shù)過(guò)多，收斂較慢，因此使用EM算法的改進(jìn)算法，即ECME算法[13]進(jìn)行計(jì)算.

將x的均值和方差記為μx和Cx，即待求參數(shù)，統(tǒng)一記為θ={μx,Cx}.完整數(shù)據(jù)的似然函數(shù)表示為L(zhǎng)(θ;x1,x2, …,xJ,yw)，簡(jiǎn)化表示為[11]

L∝

(5)

(6)

(7)

式中：ywi為yw的第j個(gè)分量.最后最大化似然函數(shù)，可以得到

(8)

由于以上是基于正態(tài)分布假設(shè)的，所以為了證明結(jié)果的可靠性需要進(jìn)行假設(shè)檢驗(yàn).對(duì)殘差進(jìn)行檢驗(yàn)，對(duì)于每個(gè)實(shí)驗(yàn)樣本點(diǎn)，殘差為

(9)

如果殘差服從標(biāo)準(zhǔn)正態(tài)分布，則說(shuō)明正態(tài)分布假設(shè)成立[8].

為測(cè)試該方法的可靠性，建立簡(jiǎn)單的數(shù)學(xué)模型對(duì)該方法進(jìn)行測(cè)試.在測(cè)試中給定70個(gè)樣本點(diǎn)，令已知分布的輸入?yún)?shù)不確定性u(píng)的維度為1，待評(píng)估的模型參數(shù)不確定性x的維度為2，敏感性系數(shù)由均勻分布U(5, 50)得到，已知分布參數(shù)服從u1～N(1.5, 0.42)，認(rèn)為實(shí)驗(yàn)誤差為ej～N(0, 0.052)，用于生成yw的真實(shí)分布滿足x1～N(-0.6, 0.32),x2～N(1.2, 0.22)，然后根據(jù)這些數(shù)據(jù)反推x的均值和方差，其計(jì)算結(jié)果如表1所示.可知，計(jì)算結(jié)果與真實(shí)值比較接近，相對(duì)誤差不超過(guò)11%.

表1 評(píng)估所得參數(shù)與真實(shí)值的比較Tab.1 Comparison of evaluated parameters and real values

1.2 輸入?yún)?shù)抽樣

不確定性的正向傳遞可以通過(guò)蒙特卡洛抽樣進(jìn)行.對(duì)確定的所有重要輸入?yún)?shù)同時(shí)抽樣，根據(jù)次序統(tǒng)計(jì)理論，抽樣數(shù)目只與輸出結(jié)果的容忍區(qū)間以及置信水平有關(guān).滿足特定容忍區(qū)間的最小抽樣數(shù)目由Wilks公式[14]確定，在雙側(cè)容忍區(qū)間的情況下，其形式為

β=1-γZ+Z(1-γ)γZ-1

(10)

式中：β為置信水平；γ為概率；Z為最小樣本數(shù)目.對(duì)于統(tǒng)計(jì)方法，美國(guó)核管會(huì)(NRC)建議采用95%的置信水平，不確定性評(píng)估必須證明在β=95%時(shí)，95%的計(jì)算結(jié)果(容忍區(qū)間)在安全裕度內(nèi)[15].因此，令β=γ=95%，對(duì)于雙側(cè)容忍區(qū)間可得Z=93，即抽樣計(jì)算93次后即可滿足“95/95準(zhǔn)則”, 得到的最大值和最小值可以認(rèn)為是雙95%置信區(qū)間的上下界.

2 PSBT空泡份額棒束基準(zhǔn)

壓水堆子通道和棒束實(shí)驗(yàn)(PSBT)是2010年發(fā)布的一套基準(zhǔn)，主要用于驗(yàn)證現(xiàn)有子通道和計(jì)算流體力學(xué)程序模擬流體在單管和棒束的空泡份額、含氣率以及偏離泡核沸騰(DNB)等.本研究計(jì)算對(duì)象為PSBT基準(zhǔn)一階段穩(wěn)態(tài)棒束空泡份額分布問(wèn)題.表2為本文計(jì)算涉及的3種組件具體參數(shù)[16].其中，B5和B6組件的徑向功率分布場(chǎng)為類型A，如圖2(a)所示；B7組件的徑向功率分布為類型B, 如圖2(b)所示.

圖3為實(shí)驗(yàn)組件截面圖和子通道的具體劃分情況.實(shí)驗(yàn)數(shù)據(jù)為位于距離加熱底端 2 216、2 669、3 177 mm共3個(gè)軸向位置的4個(gè)中心子通道的空泡份額平均值.在距離加熱底端2 216 mm的底端測(cè)量點(diǎn)的空泡份額較低，空泡聚集在受熱表面，因此實(shí)驗(yàn)確定的空泡份額將低于實(shí)際的空泡份額，具有較大誤差.對(duì)此，僅選取 2 669、3 177 mm 兩個(gè)軸向位置的測(cè)量值.根據(jù)計(jì)算，所選用穩(wěn)態(tài)實(shí)驗(yàn)工況的數(shù)據(jù)范圍為壓力7.36～14.71 MPa、質(zhì)量流量1 380.6～2 288.9 kg/(m2·s)、加熱功率 1 920～3 536 kW、入口溫度173.5～287.8 ℃.共選取30組工況，并從中選取25組工況進(jìn)行輸入模型參數(shù)的不確定性量化，利用另外5組對(duì)計(jì)算得到的不確定性分析結(jié)果進(jìn)行評(píng)估，以驗(yàn)證計(jì)算結(jié)果的有效性.

表2 組件參數(shù)Tab.2 Parameters of assemblies

①②③

從邊界條件和計(jì)算模型兩方面分析影響空泡份額分布的因素.對(duì)于邊界條件，PSBT基準(zhǔn)報(bào)告中給出了邊界條件的測(cè)量精度，參照文獻(xiàn)[17]的處理方式：壓力、入口溫度、質(zhì)量流量、熱流密度的不確定性的分布均為均值為0的正態(tài)分布，其標(biāo)準(zhǔn)差依次為0.333%、0.133%、0.5%、0.333%.

在子通道程序的計(jì)算過(guò)程中所選的重要模型和參數(shù)如表3所示.空泡份額模型采用滑移模型，滑速比可以設(shè)置.根據(jù)最小熵增原理[18]，在14 MPa下的滑速比計(jì)算結(jié)果為1.86，但是一般情況下計(jì)算值比大多數(shù)的實(shí)驗(yàn)值偏高[19]，而最典型的單相流模型是均相模型，其把兩相流體看成某種單相混合物流體，認(rèn)為氣相與液相之間不存在相對(duì)速度，即滑速比為1，因此滑速比基準(zhǔn)值最終取1.2.湍流交混模型中需要給定湍流交混系數(shù)，湍流交混系數(shù)主要與幾何尺寸、棒束排布方式、流體類型以及雷諾數(shù)(Re)范圍有關(guān)，根據(jù)與計(jì)算工況相近的Castellana等[20]的實(shí)驗(yàn)所得湍流交混系數(shù)關(guān)系式，取湍流交混系數(shù)的基準(zhǔn)值為0.008.

表3 PSBT基準(zhǔn)題計(jì)算模型及參數(shù)選擇

模型參數(shù)使用基準(zhǔn)值進(jìn)行計(jì)算.如圖4所示，橫軸為空泡份額實(shí)驗(yàn)值ve，縱軸為子通道程序空泡份額計(jì)算值vc，越接近黑色實(shí)線代表計(jì)算值與測(cè)量值越符合，上、下兩條虛線代表實(shí)驗(yàn)測(cè)量不確定性2σ(=0.08).可知，較多的計(jì)算值位于實(shí)驗(yàn)值的兩倍標(biāo)準(zhǔn)差(2σ)范圍之外，且整體預(yù)測(cè)值偏高，證明計(jì)算選用的模型過(guò)于保守，仍需要進(jìn)一步改進(jìn).

在進(jìn)行模型參數(shù)不確定性量化過(guò)程中，計(jì)算敏感性矩陣是重要步驟.在使用有限差分法求偏導(dǎo)數(shù)的過(guò)程中，步長(zhǎng)的選擇并不是越小越好.因?yàn)椴介L(zhǎng)太小可能會(huì)導(dǎo)致程序響應(yīng)變化不明顯，也會(huì)受到數(shù)值不穩(wěn)定性的影響，步長(zhǎng)過(guò)大又會(huì)導(dǎo)致誤差過(guò)大.所以為了減小偏導(dǎo)數(shù)計(jì)算過(guò)程中的誤差，需要進(jìn)行平滑處理[11].例如在求解滑速比和湍流交混系數(shù)對(duì)應(yīng)的敏感性系數(shù)時(shí)，分別取步長(zhǎng)為 ±0.01, ±0.03, ±0.1 來(lái)進(jìn)行計(jì)算，對(duì)于每個(gè)樣本點(diǎn)可以得到6個(gè)偏導(dǎo)數(shù)，然后舍掉最大值和最小值，求平均值作為最終的偏導(dǎo)數(shù).

3 計(jì)算結(jié)果分析

將壓力、入口溫度、質(zhì)量流量、熱流密度當(dāng)作已知分布參數(shù)，將滑速比和湍流交混系數(shù)作為待評(píng)估的模型參數(shù)，利用25組工況共50個(gè)實(shí)驗(yàn)樣本點(diǎn)計(jì)算得到滑速比的均值和標(biāo)準(zhǔn)差分別為 0.440 5、0.259 9；湍流交混系數(shù)的均值和標(biāo)準(zhǔn)差分別為-0.103 4、0.048 8.可知，滑速比不確定性的標(biāo)準(zhǔn)差較大.根據(jù)最小熵增理論，滑速比的大小與壓力密切相關(guān)，而在進(jìn)行模型參數(shù)的不確定性分析時(shí)所選用的工況壓力范圍為7.36～14.71 MPa，因而滑速比也存在較大的不確定性.為了驗(yàn)證正態(tài)分布假設(shè)是否成立，使用卡方檢驗(yàn)來(lái)驗(yàn)證殘差是否滿足標(biāo)準(zhǔn)正態(tài)分布，將殘差劃分為10個(gè)區(qū)間，有：

(11)

則有95%的置信度認(rèn)為參數(shù)服從正態(tài)分布.

在得到了輸入?yún)?shù)不確定性的具體分布之后，對(duì)另外5組工況進(jìn)行不確定性的正向傳播.將得到的模型參數(shù)及邊界條件參數(shù)按照其分布進(jìn)行93次抽樣，并通過(guò)子通道程序進(jìn)行計(jì)算，將得到的最大值和最小值當(dāng)作雙95%置信區(qū)間的上下界.圖5為4個(gè)中心子通道的空泡份額平均值計(jì)算值不確定帶對(duì)實(shí)驗(yàn)值的包絡(luò)情況.可知，當(dāng)空泡份額較高時(shí)，不確定帶對(duì)實(shí)驗(yàn)值的包絡(luò)性非常好；當(dāng)空泡份額較低時(shí)，不確定帶對(duì)實(shí)驗(yàn)數(shù)據(jù)的包絡(luò)情況不夠理想，有部分未包絡(luò)，但不確定帶下界也非常接近實(shí)驗(yàn)值.

利用所求出的統(tǒng)計(jì)均值對(duì)模型系數(shù)進(jìn)行修正，圖6為基準(zhǔn)模型修正前后4個(gè)中心子通道的空泡份額平均值計(jì)算結(jié)果對(duì)比.可知，樣本點(diǎn)模型參數(shù)經(jīng)過(guò)均值修正后，計(jì)算結(jié)果與實(shí)驗(yàn)值更接近，計(jì)算精度有了明顯改善.

4 結(jié)論

本文在假設(shè)輸入?yún)?shù)不確定性與輸出參數(shù)不確定性近似為線性關(guān)系的基礎(chǔ)上，將邊界條件壓力、入口溫度、質(zhì)量流量、熱流密度當(dāng)作已知分布輸入?yún)?shù)，對(duì)子通道程序COBRA-IV輸入模型參數(shù)中滑速比和湍流交混系數(shù)的不確定性進(jìn)行評(píng)估；然后進(jìn)行輸入?yún)?shù)不確定性的正向傳遞，得到了空泡份額的不確定帶，并利用求得的模型參數(shù)不確定性統(tǒng)計(jì)均值對(duì)模型進(jìn)行修正，計(jì)算結(jié)果表明：

(1) 在所選取的空泡分布實(shí)驗(yàn)工況范圍內(nèi)，滑速比存在較大的不確定性，而湍流交混系數(shù)的變化范圍較小.

(2) 計(jì)算結(jié)果的不確定帶對(duì)實(shí)驗(yàn)值的包絡(luò)情況較好，尤其是當(dāng)空泡份額較高時(shí).

(3) 利用統(tǒng)計(jì)均值對(duì)模型進(jìn)行修正后，可以得到比原模型更接近實(shí)驗(yàn)值的計(jì)算結(jié)果，對(duì)于輸入模型參數(shù)的不確定性評(píng)估結(jié)果可以作為計(jì)算模型改進(jìn)的重要參考.