杜 蕓，劉曉晶，程 旭

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

1988年，美國NRC提出了最佳估算模型的安全分析程序，針對最佳估算程序的不確定性分析方法也應(yīng)運而生。不確定性分析成為如今進行核能熱工安全分析必不可少的工作。然而，由于核能領(lǐng)域的不確定性分析方法最初是針對系統(tǒng)程序展開的，針對子通道程序開發(fā)的不確定性分析方法還較少。子通道分析程序是計算反應(yīng)堆堆芯熱工水力過程現(xiàn)象的通用工具，子通道程序的模擬計算同樣存在不確定性，對其進行不確定性分析十分必要。

目前針對系統(tǒng)程序開發(fā)的不確定性分析方法[1]，從跟蹤不確定性的方法上可分為輸入?yún)?shù)不確定性的傳播和輸出結(jié)果的誤差傳播兩類。輸入?yún)?shù)不確定性的傳播可通過如下方式得到：首先認定并給出不確定性參數(shù)的范圍和分布，然后通過改變這些輸入?yún)?shù)進行計算。輸出結(jié)果的誤差傳播可通過計算結(jié)果與實驗數(shù)據(jù)的比較直接得到。本工作采用輸入?yún)?shù)不確定性傳播法[2]對COBRA-Ⅳ程序的計算進行不確定性分析。

1 計算對象

圖1 棒束子通道分布示意圖

本文研究的對象為BFBT[3]眾多實驗中的一個。BFBT是由美國NRC與日本金融、貿(mào)易和工業(yè)部一起核準(zhǔn)，最后被經(jīng)濟合作組織(OECD)認可的一國際性工程。本實驗為模擬沸水堆燃料棒束，建立一在高壓、高溫條件下垂直的8×8棒束，棒束橫截面和各子通道分布示意圖示于圖1。組件盒內(nèi)裝有60根燃料棒，呈8×8方式排列，組件中間有一直徑為34.0 mm的不加熱的水棒，其中無流體流動。棒束的軸向功率均勻分布，徑向功率非均勻分布，徑向功率分布示于圖2。定位格架和棒束的一些相關(guān)參數(shù)列于表1。

圖2 徑向相對功率分布

表1 組件的幾何參數(shù)

本工作模擬計算的是該實驗中質(zhì)量含氣率最高的1組穩(wěn)態(tài)工況，該工況的參數(shù)列于表2。

表2 子通道計算工況

2 空泡份額計算值的影響因素

利用子通道分析程序COBRA-Ⅳ對以上實驗對象進行分析計算。如圖1所示，將冷卻劑流通橫截面劃分為80個子通道。為研究影響空泡份額計算值的因素，采用較簡單的模型對實驗進行模擬計算，作為基準(zhǔn)算例進行比較。本工作從邊界條件和計算模型兩方面分析影響空泡份額計算值的因素。

2.1 邊界條件對空泡份額的影響

雖然實驗的邊界條件由實驗設(shè)備控制且由儀器測量，理論上符合工況的設(shè)定，但測量儀器存在誤差[3]。采用單一變量原則，即每次計算時只有1個邊界條件是變量，其余仍為工況設(shè)計值，以此逐一分析邊界條件的偏差對結(jié)果的影響。

為分析每個邊界條件的偏差對結(jié)果的影響并方便比較其對結(jié)果的影響，假設(shè)每個條件均存在同樣的相對偏差，綜合參考實驗設(shè)備對邊界條件的測量誤差[3]，將該相對偏差定為±1%。相對偏差的定義為：

×100%

(1)

其中：d為相對偏差；fa為實際值；fb為基準(zhǔn)值。

分別計算邊界條件發(fā)生偏差時對子通道空泡份額計算結(jié)果造成的影響，得到相對變化量較大的子通道為31、33、48、50。這些子通道均是水棒周圍的子通道。

評估邊界條件的偏差對計算結(jié)果造成偏差的平均效應(yīng)：

εre-d(i)×100%

(2)

其中：εre-d為相對變化量；i為子通道的通道編號。

分別評估邊界條件的偏差對空泡份額計算結(jié)果造成偏差的平均效應(yīng)，結(jié)果列于表3。

表3 邊界條件對空泡份額的影響

由表3計算結(jié)果可見，對空泡份額影響最大的邊界條件是入口流體的焓，其次，出口壓力對空泡份額也有一定影響，而流量及熱流密度對空泡份額的影響較小。雖然邊界條件對空泡份額的影響不大，但邊界條件的不確定性是必然存在的，所以需考慮其對結(jié)果不確定性的影響。

2.2 計算模型對空泡份額的影響

程序中的模型及模型參數(shù)的選擇也會對計算結(jié)果產(chǎn)生影響。對COBRA-Ⅳ中涉及熱工水力計算的主要物理模型及模型參數(shù)進行研究。由于模型眾多，本文僅以空泡份額模型為例。

COBRA-Ⅳ中可供選擇的空泡份額模型[4]有以下幾種。

均相模型：

α=χvG/[χvG+(1-χ)vL]

(3)

Modified Armand模型：

α=(0.833+0.167χ)χvG/((1-χ)vL+χvG)

(4)

Chexal-Lellouche模型：

α=jG/[C0(jG+jL)+vgj]

(5)

滑移模型(滑速比可設(shè)置)：

α=χvG/[(1-χ)vLS+χvG]

(6)

其中：α為空泡份額；vG為氣體比體積；vL為液體比體積；χ為質(zhì)量含氣率；jG為氣相表觀速度；jL為液相表觀速度；vgj為漂移速度；S為滑速比；C0為表示兩種速度關(guān)系的系數(shù)[4]。

基準(zhǔn)計算采用均相模型，其余3種模型的計算結(jié)果與基準(zhǔn)計算進行對比，其中，滑移模型中的滑速比是可設(shè)置的?；谧钚幽芰骷僭O(shè)，理論滑速比為2.7，因此，本文選擇的滑速比為2.0～3.0。

圖3 不同模型的空泡份額計算結(jié)果比較

不同模型空泡份額的計算結(jié)果比較示于圖3。圖3中，H、MA和CL分別代表均相模型、Modified Armand模型和Chexal-Lellouche模型，S2.0、S2.5和S3.0分別代表滑速比為2.0、2.5和3.0的滑移模型。由圖3可見：用不同模型計算的子通道出口空泡份額與基準(zhǔn)計算結(jié)果的趨勢一致，只是大小有異；均相模型的計算結(jié)果最大，這一點由式(3)不難得出。同時，滑移模型的計算結(jié)果受滑速比的影響非常大，由于滑速比定義了流體中氣相與液相流動速度的關(guān)系，由式(6)可知該值在計算空泡份額時起重要作用。

圖3結(jié)果表明，空泡份額模型的選擇對計算結(jié)果的影響相當(dāng)大。在本文計算范圍內(nèi)，平均相對變化量高達20%，這比邊界條件對計算結(jié)果的影響大得多。可見，選擇適合工況的空泡份額模型十分重要。

對模型及模型參數(shù)依次進行研究，綜合所有模型對子通道空泡份額計算值造成的影響，將不同模型的計算結(jié)果與基準(zhǔn)計算結(jié)果進行比較。取80個子通道中相對變化量的最大值的絕對值，結(jié)果列于表4。

由表4可見，與空泡份額的計算直接相關(guān)的空泡份額模型對計算結(jié)果的影響最大，空泡漂移流模型及修正種類的選擇對計算結(jié)果的影響也較大。交混系數(shù)對計算結(jié)果的影響較小，但該影響無法避免。

表4 各模型對子通道出口空泡份額計算值的影響比較

綜上所述，在確定采用滑移模型作為空泡份額模型后，確定6個參數(shù)為不確定性輸入?yún)?shù)。其中，4個參數(shù)為邊界條件，包括出口壓力、入口焓、入口質(zhì)量流密度和熱流密度，2個參數(shù)為模型參數(shù)，包括滑速比和交混系數(shù)。邊界條件的不確定性范圍由儀器的測量精度決定，但不確定性的分布無從知道，所以均按照均勻分布來處理。本文考慮的輸入?yún)?shù)的不確定性范圍及分布列于表5。

表5 輸入?yún)?shù)的不確定性范圍及分布

3 不確定性分析

3.1 數(shù)學(xué)模型及原理

對于給定的計算容忍限，需要確定取樣數(shù)目[5]。文獻[6]提出對1個量進行雙邊容忍限(式(7))及單邊容忍限(式(8))的計算公式：

β=1-αN-N(1-α)αN-1

(7)

β=1-αN

(8)

其中：β為置信度；α為總體空間在兩個限值間的份額；N為最小采樣次數(shù)。

式(7)、(8)即為Wilks公式。

3.2 不確定性分析程序

確定了要求的容忍限，根據(jù)Wilks公式可得到需要的計算次數(shù)。該計算次數(shù)與不確定性輸入?yún)?shù)的個數(shù)、不確定性輸入?yún)?shù)的范圍及分布均無關(guān)，所以該方法適用性很廣。在滿足計算次數(shù)要求的同時，采樣方式須根據(jù)數(shù)學(xué)原理符合計算范圍以及分布，且隨機產(chǎn)生。這樣計算的結(jié)果經(jīng)適當(dāng)處理才能滿足要求的容忍限。

本工作編寫了子通道程序的不確定性分析程序，該程序的計算流程示于圖4。采用排序的方式對N組計算結(jié)果進行處理，將每個子通道的計算結(jié)果按升序排序，其中，排在中間的值作為該子通道的預(yù)測結(jié)果，排在第1位的值為滿足容忍限的不確定性下限，排在最后1位的值為不確定性的上限。這樣，便完成了對該子通道程序計算的不確定性分析。

圖4 不確定性分析程序的流程

3.3 程序驗證

運用本工作編寫的不確定性分析程序產(chǎn)生93組程序輸入?yún)?shù)。由于SUSA方法也使用了Wilks公式的原理，將本文不確定性分析程序產(chǎn)生的數(shù)據(jù)與SUSA方法的數(shù)據(jù)進行比較，驗證本方法的可行性。

1) 驗證各不確定性輸入?yún)?shù)是否符合分布規(guī)律

以入口焓采樣為例，用已排序的數(shù)據(jù)進行對比，結(jié)果示于圖5。由圖5可見，入口焓的采樣值分布近似呈線性，且采樣點疏密較一致，重合率較高，說明兩者采樣均是均勻分布。因此，本文的采樣能基本符合參數(shù)分布的特點。同樣，對于其余5個不確定性輸入?yún)?shù)，本文采樣與SUSA方法的也吻合較好，且符合每個參數(shù)分布的特點。

圖5 入口焓的采樣排序?qū)Ρ?/p>

2) 各參數(shù)組合比較

為比較6個輸入?yún)?shù)的組合方式，將每一種參數(shù)的93個采樣值劃分為大于基準(zhǔn)值和小于等于基準(zhǔn)值兩種，小于等于基準(zhǔn)值的為0類，大于基準(zhǔn)值的為1類。對程序產(chǎn)生的每一組參數(shù)組合按入口焓、入口流量、出口壓力、熱流密度、滑速比和交混系數(shù)的順序編號，表示每組參數(shù)組合的特征。如第1組參數(shù)的編號為011101，說明第1組參數(shù)的特征為：入口焓小于等于基準(zhǔn)值，入口流量大于基準(zhǔn)值，出口壓力大于基準(zhǔn)值，熱流密度大于基準(zhǔn)值，滑速比小于等于基準(zhǔn)值，交混系數(shù)大于基準(zhǔn)值。將編號作為一二進制數(shù)，便可轉(zhuǎn)換為一十進制數(shù)，共有64種組合方式，分別用0到63表示。這樣，第1組參數(shù)的特征可用29來表示。將SUSA和本文方法對于不確定性參數(shù)的組合用數(shù)字表示，并進行比較，結(jié)果示于圖6。

由圖6可知：SUSA方法覆蓋的組合類型為51種，本文方法覆蓋的組合類型為50種，覆蓋率均約為80%，二者均未覆蓋所有組合；本文的組合類型與SUSA方法的組合類型相比，相似率為80%。

綜上可見，本文采樣方法與組合方式是合格的，程序計算出的結(jié)果經(jīng)順序處理，能滿足容忍限的要求。

圖6 不同種類組合數(shù)目的比較

3.4 不確定性分析結(jié)果

由于該實驗屬高空泡份額的實驗，滑移模型較適合計算這類空泡份額，因此選定滑移模型作為本計算的空泡份額模型。其余模型則根據(jù)棒束實驗，選擇較適合實驗的模型即可。將雙邊容忍限定為(95%，95%)，采用自行編寫的不確定性分析程序，得到最終的計算結(jié)果，即每個子通道的出口空泡份額預(yù)測值y，滿足(95%，95%)的不確定性上限y95/95上和不確定性下限y95/95下。

圖7示出子通道空泡份額的計算預(yù)測值與實驗值的對比。圖8示出空泡份額預(yù)測值與實驗值的相對誤差。由圖7、8可知，熱流密度較高的棒束周圍的子通道的空泡份額較高，說明子通道程序?qū)张莘蓊~的預(yù)測符合實驗的規(guī)律以及趨勢。預(yù)測值普遍比實驗值小，平均相對誤差為10%，最大相對誤差在20%以內(nèi)，均在可接受范圍內(nèi)。與實驗值相對誤差較大的是受熱不均勻、不對稱的子通道。

圖7 子通道空泡份額預(yù)測值與實驗值的對比

對于堆芯棒束的熱工水力問題，最關(guān)心的是通道內(nèi)冷卻劑的換熱能力能否得以保證，冷卻劑能否及時有效地將棒束產(chǎn)生的熱量帶走。空泡的存在會削弱冷卻劑的換熱能力，且空泡份額過大時有可能發(fā)生傳熱惡化現(xiàn)象，導(dǎo)致包殼溫度過高，堆芯安全受到威脅。觀察棒束的實驗結(jié)果以及計算結(jié)果(圖8)，發(fā)現(xiàn)程序?qū)τ诔隹诟呖张莘蓊~的通道預(yù)測較準(zhǔn)確。高功率棒束周圍的子通道普遍表現(xiàn)出高的出口空泡份額，包括邊通道和中間通道兩種。按照子通道周圍棒束功率的分布不同，又可將這兩種通道細分。表6列出選擇的具有代表性的子通道。

圖8 空泡份額預(yù)測值與實驗值的相對誤差

表7列出子通道79的不確定性分析結(jié)果。該結(jié)果表示，在考慮了表5所列的不確定性后，該程序?qū)Π羰鴮嶒炛凶油ǖ?9出口空泡份額的計算，其結(jié)果為75.623%的可能性最大，有95%的可能性落在71.787%～79.966%之間，而這個判斷的可信度為95%。由圖2所示的功率分布不難看出，這5種子通道相比較，子通道17周圍的棒束功率最大，導(dǎo)致其空泡份額較高，符合本文分析的規(guī)律。

表6 8×8棒束實驗中典型高空泡子通道

表7 典型高空泡子通道的出口空泡份額計算不確定性范圍

計算了表6中5個子通道的不確定性，結(jié)果示于圖9。由圖9可見，高空泡子通道的出口空泡份額的計算不確定性差別很小，均在約-5.5%～6%之間。

圖9 典型高空泡子通道出口空泡份額的不確定性

子通道空泡份額的計算不確定性示于圖10。由圖10可見，不同子通道的計算不確定性差別很大。其中，邊角子通道的計算不確定性較小，約為±5.5%；水棒周圍不規(guī)則形狀的子通道的不確定性較大，約為±9%。因該類型子通道空泡份額的基礎(chǔ)小，受熱不均勻，且受熱較小，此時功率變化及交混量等因素的變化對其造成的影響不可忽略，因此不確定性較大。由此可見，該程序?qū)Ω呖张葑油ǖ赖某隹诳张萦嬎愕牟淮_定性較小，但由于該類子通道的空泡份額基礎(chǔ)較大，少許不確定性變動均有可能對堆芯安全造成威脅，因此應(yīng)格外重視此類子通道空泡份額的不確定性的計算。

圖10 子通道空泡份額的計算不確定性

圖11 空泡份額的計算不確定帶與實驗值的比較

空泡份額的計算不確定帶與實驗值的比較示于圖11。該不確定帶滿足(95%，95%)的容忍限。然而，其結(jié)果并不能完全包絡(luò)實驗值，說明此次計算的誤差較大。其中，與實驗值差別較大的子通道有10、48、50、56和60，這幾個子通道的共同特征是其周圍棒束的加熱功率均不一致，即子通道周圍的受熱不對稱、不均勻。尤其是水棒周圍的子通道，其一面受到不斷加熱，而另一面又完全沒有受熱，這種極度的不對稱造成計算與實驗偏差較大。而子通道22、30、39、42的形狀規(guī)則，且是受熱對稱的通道，其計算的偏差較小且穩(wěn)定。這說明在計算受熱不均勻的子通道時，子通道程序的計算能力還有待提高。

除此之外，造成不確定性范圍沒有包絡(luò)實驗值的原因還可能是：1) 基準(zhǔn)工況的模型及模型選擇還存在不當(dāng)之處；2) 模型參數(shù)的不確定性范圍無從得知，本文選定范圍并不合適，從而導(dǎo)致計算出的不確定性范圍與客觀事實不符。以上問題還需進一步探討。

4 小結(jié)

本文選擇BFBT的8×8棒束實驗作為計算對象，運用子通道分析程序?qū)Τ隹谔幍目张莘蓊~進行了簡單的分析計算，針對子通道程序進行了不確定性分析，得到的結(jié)論如下。

1) 子通道程序?qū)τ谶吔亲油ǖ阑蚴軣岵痪鶆蜃油ǖ莱隹诳张莘蓊~的計算相對于實驗值的偏差較大。

2) 當(dāng)邊界條件有微小變化時，出口空泡份額的變化量也較小。相比之下，入口焓的變化對結(jié)果影響最大。其次是出口壓力及熱流密度。對于出口空泡份額，子通道分析程序中的空泡漂移模型及修正種類對結(jié)果的影響較大，前者的最大影響高達25.4%，后者的最大影響為6.5%。在計算空泡份額時，這兩種模型的選擇非常關(guān)鍵。

3) 運用Wilks公式的原理確定能滿足容忍限的最少計算次數(shù)，且運用順序統(tǒng)計法對計算結(jié)果進行處理，得到計算值的不確定帶分布范圍。結(jié)果顯示，邊角子通道的計算不確定性較小，約為±5.5%；水棒周圍不規(guī)則形狀的子通道的不確定性較大，約為±9%。對于高空泡子通道的出口空泡份額，其不確定性在-5.5%～6% 之間，與其他子通道相比較小，但由于該類子通道的空泡份額基礎(chǔ)大，少許不確定性變動均有可能對堆芯安全造成威脅，因此應(yīng)格外重視對該類子通道空泡份額的不確定性計算。所有子通道的不確定帶均未完全包絡(luò)實驗值，說明程序計算還存在較大偏差，需進一步改進。

參考文獻：

[1] Best estimate safety analysis for nuclear power plants: Uncertainty evaluation[R]. US: IAEA, 2008.

[2] BOYACK B E, CATTON I, DUFFEY R B, et al. Quantifying reactor safety margins, Part 1: An overview of the code scaling, applicability, and uncertainty evaluation methodology[J]. Nuclear Engineering and Design, 1990, 119(1): 1-15.

[3] NUPEC BWR full-size fine-mesh bundle test benchmark, Volume Ⅰ: Specifications[R]. US: NRC OECD Nuclear Energy Agency, 2005.

[4] 徐濟鋆. 沸騰傳熱和氣液兩相流[M]. 北京：原子能出版社，2001.

[5] WILKS S S. Determination of sample sizes for setting tolerance limits[J]. The Annals of Mathematical Statistics, 1941, 12(1): 91-96.

[6] GUBA A. Statistical aspects of best estimate method-Ⅰ[J]. Reliability Engineering and System Safety, 2003, 80(3): 217-232.