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

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

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

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

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

1 不確定性量化分析方法

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

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

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

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

(1)

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

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

(2)

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

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

y=h(x,u)

(3)

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

yw=Sxx+Suu+e

(4)

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

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

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

L∝

(5)

(6)

(7)

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

(8)

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

(9)

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

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

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

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

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

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

(10)

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

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

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

圖3為實驗組件截面圖和子通道的具體劃分情況.實驗數(shù)據(jù)為位于距離加熱底端 2 216、2 669、3 177 mm共3個軸向位置的4個中心子通道的空泡份額平均值.在距離加熱底端2 216 mm的底端測量點的空泡份額較低，空泡聚集在受熱表面，因此實驗確定的空泡份額將低于實際的空泡份額，具有較大誤差.對此，僅選取 2 669、3 177 mm 兩個軸向位置的測量值.根據(jù)計算，所選用穩(wěn)態(tài)實驗工況的數(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組工況進行輸入模型參數(shù)的不確定性量化，利用另外5組對計算得到的不確定性分析結(jié)果進行評估，以驗證計算結(jié)果的有效性.

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

①②③

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

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

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

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

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

3 計算結(jié)果分析

將壓力、入口溫度、質(zhì)量流量、熱流密度當(dāng)作已知分布參數(shù)，將滑速比和湍流交混系數(shù)作為待評估的模型參數(shù)，利用25組工況共50個實驗樣本點計算得到滑速比的均值和標(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)，而在進行模型參數(shù)的不確定性分析時所選用的工況壓力范圍為7.36～14.71 MPa，因而滑速比也存在較大的不確定性.為了驗證正態(tài)分布假設(shè)是否成立，使用卡方檢驗來驗證殘差是否滿足標(biāo)準(zhǔn)正態(tài)分布，將殘差劃分為10個區(qū)間，有：

(11)

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

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

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

4 結(jié)論

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

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

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

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