COBRA-Ⅳ對(duì)8×8棒束計(jì)算的不確定性分析

杜 蕓，劉曉晶，程 旭

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

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

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

1 計(jì)算對(duì)象

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

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

圖2 徑向相對(duì)功率分布

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

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

表2 子通道計(jì)算工況

2 空泡份額計(jì)算值的影響因素

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

2.1 邊界條件對(duì)空泡份額的影響

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

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

×100%

(1)

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

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

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

εre-d(i)×100%

(2)

其中：εre-d為相對(duì)變化量；i為子通道的通道編號(hào)。

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

表3 邊界條件對(duì)空泡份額的影響

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

2.2 計(jì)算模型對(duì)空泡份額的影響

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

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)計(jì)算采用均相模型，其余3種模型的計(jì)算結(jié)果與基準(zhǔn)計(jì)算進(jìn)行對(duì)比，其中，滑移模型中的滑速比是可設(shè)置的?；谧钚?dòng)能流假設(shè)，理論滑速比為2.7，因此，本文選擇的滑速比為2.0～3.0。

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

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

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

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

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

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

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

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

3 不確定性分析

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

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

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

(7)

β=1-αN

(8)

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

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

3.2 不確定性分析程序

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

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

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

3.3 程序驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

圖7 子通道空泡份額預(yù)測(cè)值與實(shí)驗(yàn)值的對(duì)比

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

圖8 空泡份額預(yù)測(cè)值與實(shí)驗(yàn)值的相對(duì)誤差

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

表6 8×8棒束實(shí)驗(yàn)中典型高空泡子通道

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

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

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

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

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

圖11 空泡份額的計(jì)算不確定帶與實(shí)驗(yàn)值的比較

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

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

4 小結(jié)

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

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

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

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

參考文獻(xiàn)：

[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] 徐濟(jì)鋆. 沸騰傳熱和氣液兩相流[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.