基于高精度子通道程序的棒束臨界熱流密度機(jī)理模型

郭俊良,蔣 理,高 逵,孔煥俊,桂 淼,單建強(qiáng)

(西安交通大學(xué) 動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710049)

為了提高壓水堆核電站的經(jīng)濟(jì)性,過冷沸騰因其高熱流密度和高換熱系數(shù)成為首選的沸騰模式[1]。由于堆芯內(nèi)工質(zhì)的低含汽率和泡狀流流型,容易發(fā)生偏離核態(tài)沸騰(DNB)。當(dāng)發(fā)生DNB型沸騰危機(jī)時(shí),燃料包殼表面被大量汽膜覆蓋導(dǎo)致傳熱能力的惡化,嚴(yán)重威脅著反應(yīng)堆的安全。因此,準(zhǔn)確預(yù)測棒束臨界熱流密度(CHF)并提高CHF是壓水堆設(shè)計(jì)的最終目標(biāo),以提高反應(yīng)堆的經(jīng)濟(jì)性和安全性[2]。

由于壓水堆堆芯燃料組件復(fù)雜且獨(dú)特的幾何形狀(如緊湊的棒束結(jié)構(gòu)、攪混翼格架和導(dǎo)向管)以及寬范圍的運(yùn)行工況,準(zhǔn)確預(yù)測棒束CHF具有一定的挑戰(zhàn)性。目前,工程上常用的棒束CHF預(yù)測方法是針對(duì)感興趣的工況開展大規(guī)模棒束CHF實(shí)驗(yàn)并利用其數(shù)據(jù)開發(fā)特定棒束結(jié)構(gòu)的CHF關(guān)聯(lián)式。使用子通道分析程序建立特定棒束幾何結(jié)構(gòu)的數(shù)學(xué)物理模型,并采用集總參數(shù)的方法將子通道內(nèi)的“平均”參數(shù)作為“局部”參數(shù)輸入到CHF關(guān)聯(lián)式中擬合出最優(yōu)系數(shù)。工程上常用的關(guān)聯(lián)式如EPRI[3]、WRB-2[4]以及FC-2000[5]均采用上述方法得到。然而,現(xiàn)有的子通道分析程序?qū)ψ油ǖ赖膭澐州^為粗糙,不能準(zhǔn)確刻畫燃料棒表面附近局部流場的細(xì)節(jié),特別是攪混翼格架和導(dǎo)向管存在的情況下。而DNB型CHF又是一個(gè)局部現(xiàn)象,因此,需要在CHF關(guān)聯(lián)式中反映攪混翼格架和導(dǎo)向管(冷壁效應(yīng))的影響[5]。CHF關(guān)聯(lián)式的應(yīng)用受限于棒表面局部參數(shù)計(jì)算的準(zhǔn)確性,這也是為什么特定的CHF關(guān)系式只能用于特定棒束的原因。準(zhǔn)確計(jì)算棒表面的局部參數(shù)是準(zhǔn)確預(yù)測CHF的前提。

除了準(zhǔn)確計(jì)算棒表面局部參數(shù)外,棒束CHF預(yù)測面臨的另一大挑戰(zhàn)是CHF的預(yù)測方法,包括機(jī)理模型、關(guān)聯(lián)式以及查詢表[6]。CHF機(jī)理模型以物理機(jī)制為基礎(chǔ),可以適當(dāng)?shù)耐馔撇⑻峁┱_的參數(shù)趨勢,一直是眾多學(xué)者研究的對(duì)象。此外,CHF機(jī)理模型可以為新型燃料組件或新型堆提供適當(dāng)?shù)念A(yù)測。在過去的幾十年中,針對(duì)不同熱工水力條件已經(jīng)建立了許多CHF機(jī)理模型[7-9]。其中Weisman等[7]的模型適用于最廣泛的工況范圍,并且已被應(yīng)用于棒束CHF的預(yù)測。然而,Weisman等的模型在應(yīng)用于棒束CHF預(yù)測時(shí)受限于子通道局部參數(shù)的不準(zhǔn)確性,預(yù)測精度不高,且對(duì)于導(dǎo)向管柵元CHF預(yù)測時(shí)仍需進(jìn)行冷壁效應(yīng)的額外修正[10]。因此,本研究基于高精度子通道分析程序建立一個(gè)棒束CHF機(jī)理模型,并試圖擺脫攪混翼格架和導(dǎo)向管對(duì)CHF模型的修正,使CHF的預(yù)測重新回到局部現(xiàn)象的本質(zhì)。

1 高精度子通道分析程序

ATHAS(Advanced Thermal-Hydraulics Analysis Subchannel)是西安交通大學(xué)開發(fā)的具有自主知識(shí)產(chǎn)權(quán)的熱工水力子通道分析程序。該子通道分析程序基于漂移流模型,采用全隱式壓力-速度修正算法求解,具有較高的預(yù)測模擬精度。最近Dong等[11]基于分布式阻力模型開發(fā)了一種用于ATHAS的攪混翼橫流混合模型;Jiang等[12]同樣基于分布式阻力模型建立了高精度子通道分析程序。本文簡要介紹高精度子通道分析程序的子通道劃分和相關(guān)本構(gòu),具體細(xì)節(jié)參考文獻(xiàn)[12]。

1.1 精細(xì)化子通道劃分

在子通道分析中,人為地將棒束通道劃分成若干個(gè)子通道并假設(shè)每個(gè)子通道內(nèi)的熱工水力行為像圓管一樣,且子通道大小的劃分是任意的。如果子通道劃分得太細(xì),則必須考慮兩相微觀結(jié)構(gòu)的細(xì)節(jié);如果子通道劃分得過于粗糙,則無法準(zhǔn)確預(yù)測感興趣的棒表面的局部流動(dòng)參數(shù)[13],精細(xì)化子通道劃分就是在這種情況下提出的。在傳統(tǒng)子通道劃分的基礎(chǔ)上,根據(jù)攪混翼格架葉片分布的特點(diǎn),將子通道進(jìn)一步細(xì)分。傳統(tǒng)子通道和精細(xì)化子通道的劃分方法分別如圖1所示。

a——傳統(tǒng)子通道劃分;b——精細(xì)化子通道劃分

精細(xì)化子通道劃分后,每個(gè)子通道內(nèi)至多存在1個(gè)攪混翼葉片,因此能夠更好地適用于格架分布式阻力模型。另外,相比于傳統(tǒng)子通道,更細(xì)致的劃分可以相對(duì)更準(zhǔn)確的計(jì)算出棒表面的局部參數(shù)且能得到傳統(tǒng)子通道內(nèi)部的橫流。雖然更細(xì)致的子通道劃分使子通道和間隙數(shù)量都大大增加,但計(jì)算效率仍遠(yuǎn)高于CFD類軟件。

1.2 攪混翼格架分布式阻力模型

分布式阻力模型最早應(yīng)用于鈉冷快堆繞絲形阻的分析,通過受力分析將繞絲對(duì)流體的作用力加入到動(dòng)量方程中[14]。類似的,也可將燃料棒和攪混翼格架對(duì)流體的作用力用表達(dá)式量化以反映其對(duì)流體的作用,這就是攪混翼格架分布式阻力模型的基本思想。流體在經(jīng)過攪混翼格架所在的軸向控制體內(nèi)分別受到來自燃料棒以及攪混翼葉片的作用力。將燃料棒的阻力沿軸向和徑向分解為軸向力Fu和橫向力Fv,攪混翼葉片的作用力沿葉片切向和法向分解為Fu1和Fv1。流體的受力示意圖如圖2所示。用同樣的方式將流體的總速度vall以兩種方式分解(軸向和橫向分解,切向和法向分解),如圖3所示。

圖2 流體受力示意圖

圖3 速度分量示意圖

利用Rehme軸向阻力關(guān)系式[15]和Gunter-Shaw橫向阻力關(guān)系式[14]計(jì)算速度對(duì)應(yīng)方向的阻力并加入到軸向和橫向動(dòng)量方程中,對(duì)應(yīng)的關(guān)系式如下:

(1)

(2)

(3)

(4)

式中:AR、AM分別為燃料棒和攪混翼葉片的潤濕面積;fR、fM分別為燃料棒的軸向摩擦系數(shù)和攪混翼的切向摩擦系數(shù);CR、CM分別為燃料棒的橫向摩擦系數(shù)和攪混翼的法向阻力系數(shù);ρ為流體密度;u、v分別為軸向分速度和橫向分速度;u1、v1分別為切向分速度和法向分速度。

此外,傳統(tǒng)子通道分析方法大多假設(shè)子通道中軸向的流動(dòng)占主要地位,忽略橫流對(duì)流項(xiàng)。相關(guān)研究表明[16],攪混翼所引起的定向橫流可以達(dá)到與軸流相同的數(shù)量級(jí),忽略橫流對(duì)流項(xiàng)顯然是不合理的。因此,高精度子通道分析程序在動(dòng)量方程中加入了橫流對(duì)流項(xiàng),同時(shí)考慮攪混翼所在的子通道與其下游子通道之間的不完全耦合。橫流對(duì)流項(xiàng)表達(dá)式為:

GSK=〈N〉CSK·ρAKvK

(5)

式中:GSK為橫流對(duì)流項(xiàng);〈N〉表示橫流方向;CSK為不完全耦合系數(shù);AK為間隙處的橫向面積;vK為間隙處的流速。

1.3 湍流攪混模型

相鄰子通道之間除了橫向壓差引起的橫流以及攪混翼引起的定向橫流外,湍流攪混也是子通道間質(zhì)量、動(dòng)量以及能量交換的一種重要作用機(jī)制。子通道之間的湍流攪混本質(zhì)上是以渦的形式進(jìn)行傳遞的。但由于壁面阻礙導(dǎo)致徑向渦的擴(kuò)散受到制約,其擴(kuò)散系數(shù)小于周向渦擴(kuò)散系數(shù),渦擴(kuò)散系數(shù)呈各向異性[17]。對(duì)于傳統(tǒng)子通道劃分所形成的第一類間隙而言(圖1a),湍流攪混方向平行于棒表面(周向攪混)。而對(duì)于精細(xì)化子通道劃分所新形成的第二類間隙(圖1b),湍流攪混方向垂直于棒表面(徑向攪混)。因此,這兩類間隙的湍流攪混速率是不同的。ATHAS中的湍流攪混模型分別計(jì)算兩類不同間隙的湍流攪混速率。

第一類間隙的湍流攪混速率采用Carlucci等[18]的模型來計(jì)算?？偟耐牧鲾嚮焖俾视善鄶嚮旌鸵合鄶嚮靸刹糠纸M成,而每一部分又分為均相部分和汽泡引起的附加部分:

w′l=w′l,hom+w′l,inc

(6)

w′g=w′g,hom+w′g,inc

(7)

式中:w′l、w′g分別為液相、汽相攪混速率;w′l,hom、w′g,hom分別為液相、汽相均勻部分?jǐn)嚮焖俾?w′l,inc、w′g,inc分別為液相、汽相附加部分?jǐn)嚮焖俾省?/p>

對(duì)于第二類間隙,考慮到1個(gè)棒與其周圍4個(gè)子通道中第二類間隙的攪混類似于環(huán)管的徑向攪混[19],因此第二類間隙處的湍流攪混速率采用基于Levy[20]提出的環(huán)管渦擴(kuò)散系數(shù)表達(dá)式以及普朗特混合長度理論推導(dǎo)得到:

(8)

式中:w′為單相攪混速率;l為棒表面到間隙的距離;τ為剪應(yīng)力;Sij和zij分別為間隙寬度和相鄰子通道的質(zhì)心距。

兩相湍流攪混倍增因子采用Beus[21]模型計(jì)算,在得到總的兩相湍流攪混速率后用含汽率進(jìn)行加權(quán)從而分別得到汽相和液相的湍流攪混速率。

高精度子通道分析程序ATHAS已經(jīng)過GE 3×3棒束實(shí)驗(yàn)[22]及PSBT基準(zhǔn)題[23]的驗(yàn)證,計(jì)算值與實(shí)驗(yàn)值符合很好,與傳統(tǒng)子通道計(jì)算結(jié)果相比精度大幅度提高。

2 棒束CHF機(jī)理模型開發(fā)

2.1 圓管汽泡壅塞模型

Weisman等[7]將整個(gè)流動(dòng)區(qū)域劃分成汽泡層區(qū)和主流區(qū),如圖4所示。汽泡層區(qū)與主流區(qū)之間的湍流交換被認(rèn)為是CHF觸發(fā)的核心機(jī)制。當(dāng)汽泡層內(nèi)的空泡份額超過臨界空泡份額(α2=0.82)時(shí),大量汽泡壅塞在一起從而阻礙了渦向壁面的徑向輸運(yùn),導(dǎo)致CHF的發(fā)生。

圖4 汽泡層與主流區(qū)輸運(yùn)示意圖

對(duì)汽泡層區(qū)和主流區(qū)分別建立質(zhì)量守恒方程,聯(lián)立可得:

G′(x2-x1)=qb/hfg

(9)

式中:G′為兩個(gè)區(qū)域交界面處的徑向脈動(dòng)質(zhì)量流速;x1為主流區(qū)的真實(shí)含汽率;x2為汽泡層區(qū)的真實(shí)含汽率,由α2=0.82計(jì)算得到;qb為用來產(chǎn)生蒸汽的熱流密度;hfg為汽化潛熱。

qb采用Lahey等[13]的模型計(jì)算:

(10)

式中:q為總熱流密度;hl和hf分別為主流區(qū)過冷液體焓和飽和液焓;hld為汽泡脫離點(diǎn)的焓,由Levy[24]模型計(jì)算得到。

將式(10)代入式(9)并整理成無量綱的形式,即可得到Weisman模型的CHF(qCHF)表達(dá)式:

(11)

交界面處的G′由湍流強(qiáng)度與軸向質(zhì)量流速計(jì)算得到,并考慮壁面產(chǎn)生的蒸汽會(huì)阻礙主流到汽泡層的徑向脈動(dòng)。通過引入?yún)?shù)Ψ得到有效湍流脈動(dòng)到壁面的比例。在交界面處,G′的計(jì)算表達(dá)式為:

G′=GibΨ

(12)

式中:G為軸向質(zhì)量流速;ib為徑向湍流脈動(dòng)強(qiáng)度;Ψ為一個(gè)關(guān)于壁面蒸汽速度的復(fù)雜函數(shù)。

徑向湍流脈動(dòng)強(qiáng)度由圓管單相徑向脈動(dòng)速度分布擬合得到[7],同時(shí)考慮兩相脈動(dòng)倍增因子。徑向湍流脈動(dòng)強(qiáng)度ib為:

[1+a(ρl-ρg)/ρg]

(13)

式中:k為經(jīng)驗(yàn)系數(shù),文獻(xiàn)[7]中為2.4;Db為汽泡直徑;D為圓管直徑;a為與質(zhì)量流速相關(guān)的經(jīng)驗(yàn)系數(shù);ρl和ρg分別為液相和汽相的密度。

2.2 相關(guān)本構(gòu)的修改

1) 棒束近壁區(qū)的湍流強(qiáng)度

圖5 等效圓管示意圖

由于棒表面到零剪應(yīng)力位置的距離不是定值,而子通道分析是采用集總參數(shù)的方法,須對(duì)零剪應(yīng)力位置距離進(jìn)行平均。通過敏感性分析計(jì)算,用不同的權(quán)(包括使用最大距離和最小距離的算數(shù)平均)對(duì)距離進(jìn)行加權(quán),結(jié)果變化最大為2%。因此,仿照子通道的水力直徑的定義,以(rm(θ)+ri(θ))為權(quán)重進(jìn)行加權(quán),推導(dǎo)等效圓管直徑如下:

(14)

從式(14)可看出,等效圓管的平均直徑與子通道的水力直徑定義類似,只不過分母為濕周與第二類間隙長度的總和。以此平均直徑替換式(13)中的圓管直徑D可將圓管表面的湍流強(qiáng)度分布轉(zhuǎn)換成棒表面的湍流強(qiáng)度分布,相比于之前的方法更具有物理意義。

需要強(qiáng)調(diào)的是,子通道的水力直徑和等效圓管直徑分別獨(dú)立,在計(jì)算子通道的局部參數(shù)時(shí)使用的仍是子通道的水力直徑,在CHF模型中計(jì)算湍流強(qiáng)度分布時(shí)使用的是等效圓管直徑,該直徑只是為了將圓管湍流強(qiáng)度分布轉(zhuǎn)成棒束湍流強(qiáng)度分布。CHF是一個(gè)局部現(xiàn)象,子通道分析程序在CHF預(yù)測中的作用是盡可能來提供相對(duì)更準(zhǔn)確的、集總的局部參數(shù),棒表面的湍流強(qiáng)度分布不應(yīng)受子通道的劃分而隨之改變,因此獨(dú)立開來。

2) CHF位置處的汽泡直徑

在Weisman的CHF機(jī)理模型中,采用Levy[24]模型來計(jì)算汽泡脫離直徑,并將此直徑作為CHF發(fā)生位置處的汽泡直徑。Levy模型是基于單個(gè)附著在壁面上的汽泡的受力平衡模型,其計(jì)算表達(dá)式如下:

(15)

使用該模型計(jì)算汽泡脫離點(diǎn)的直徑是合理的,但在CHF發(fā)生位置處由于大量汽泡壅塞在一起,汽泡被壓扁成長徑比為3/1的橢圓形汽泡(依據(jù)Wesiman的假設(shè))。因此,式(13)中的汽泡直徑不能用汽泡脫離直徑代替。考慮到CHF發(fā)生位置處有大量汽泡滯留在壁面,壁面上的汽泡層厚度可被認(rèn)為是粗糙度ε[27]。

按照Lee等[28]的假設(shè),汽泡層邊緣被認(rèn)為當(dāng)渦的大小為汽泡直徑Db的k倍時(shí),汽泡層便終止,因此汽泡層厚度yc可由式(16)計(jì)算得到:

(16)

式中,F2為兩相混合普朗特長度倍增因子,F2=1。

將式(16)中的汽泡層厚度yc作為壁面粗糙度ε代入到式(15)進(jìn)行迭代求解,即可得到CHF發(fā)生位置處的汽泡直徑。

3 棒束CHF機(jī)理模型評(píng)價(jià)

將棒束CHF機(jī)理模型以子程序的方式嵌入到ATHAS中。CHF子程序中所涉及到的局部參數(shù)通過ATHAS的流場計(jì)算得到,CHF子程序在ATHAS計(jì)算流場收斂后被調(diào)用。輸入?yún)?shù)為子通道程序計(jì)算局部參數(shù)值和流體物性值,輸出是預(yù)測的CHF。

3.1 棒束實(shí)驗(yàn)數(shù)據(jù)

用來驗(yàn)證本機(jī)理模型的實(shí)驗(yàn)數(shù)據(jù)來自西安交通大學(xué)棒束CHF實(shí)驗(yàn)數(shù)據(jù)庫。由于本模型是針對(duì)DNB型CHF開發(fā)的,選取其中空泡份額小于0.8的5×5棒束CHF數(shù)據(jù)進(jìn)行評(píng)價(jià)。驗(yàn)證數(shù)據(jù)共有12組,總共601個(gè)棒束CHF數(shù)據(jù)點(diǎn)。實(shí)驗(yàn)數(shù)據(jù)工況范圍列于表1,實(shí)驗(yàn)工況范圍覆蓋了壓水堆的典型工況范圍。表2列出實(shí)驗(yàn)布置情況及數(shù)據(jù)統(tǒng)計(jì)。表2中:N/A表示沒有導(dǎo)向管;D為最后一個(gè)攪混翼格架距加熱長度末端距離;gsp為攪混翼格架跨距;C/P為徑向功率因子(中間棒和外圈棒的功率比值)。

表1 實(shí)驗(yàn)數(shù)據(jù)工況

表2 實(shí)驗(yàn)布置情況及數(shù)據(jù)統(tǒng)計(jì)

3.2 棒束CHF機(jī)理模型評(píng)價(jià)

圖6示出601個(gè)數(shù)據(jù)的預(yù)測CHF與實(shí)驗(yàn)CHF的比較。從圖6可看出,幾乎所有數(shù)據(jù)都在±15%的相對(duì)誤差范圍帶內(nèi)。預(yù)測值與實(shí)驗(yàn)值之比(P/M)的平均值Rav為0.99,均方差RMS為4.69%。由于高精度子通道程序?qū)鹘y(tǒng)子通道一分為四并且增加了格架分布式阻力模型,因此可更準(zhǔn)確地計(jì)算出傳統(tǒng)子通道內(nèi)局部參數(shù)的不均勻性,所計(jì)算出來的局部參數(shù)更能表征棒表面附近的流場。此外,結(jié)合本研究所開發(fā)的機(jī)理模型,均方差RMS可達(dá)到5%以內(nèi),優(yōu)于目前任何一個(gè)棒束CHF機(jī)理模型。

圖6 預(yù)測CHF與實(shí)驗(yàn)CHF的比較

圖7示出P/M隨壓力p、質(zhì)量流速G、臨界含汽率x的分布以檢驗(yàn)其無偏性。從圖7可見,P/M隨質(zhì)量流速和臨界含汽率沒有明顯的參數(shù)傾向性,但似乎本模型在低壓下(6～8 MPa)預(yù)測偏保守。這是由于本模型是基于Weisman模型進(jìn)一步開發(fā)得到的,而在這個(gè)壓力范圍內(nèi)Weisman模型預(yù)測圓管CHF本身就偏保守[29]。

圖7 P/M隨壓力、質(zhì)量流速和臨界含汽率的分布

表3列出不帶導(dǎo)向管(典型柵元)和帶導(dǎo)向管(導(dǎo)向管柵元)的分組誤差統(tǒng)計(jì),可以看出,本模型可很好地適用于帶導(dǎo)向管和不帶導(dǎo)向管的棒束CHF預(yù)測。冷壁效應(yīng)的本質(zhì)為在一個(gè)子通道內(nèi)焓(或含汽率)存在分布,在用傳統(tǒng)子通道計(jì)算棒束CHF時(shí)需進(jìn)行冷壁效應(yīng)修正。在將子通道細(xì)分后,則無需在CHF模型上進(jìn)行任何修正即可成功預(yù)測帶導(dǎo)向管的棒束CHF。本模型從機(jī)理層面上揭示了冷壁效應(yīng)的本質(zhì)。

表3 冷壁無偏性檢驗(yàn)

4 結(jié)論

本文基于高精度子通道程序ATHAS開發(fā)了棒束CHF機(jī)理模型,并采用5×5棒束CHF實(shí)驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行了驗(yàn)證,主要結(jié)論如下。

1) 傳統(tǒng)子通道劃分方式過于粗糙,得到的集總參數(shù)不能表征棒表面的局部參數(shù)。基于CHF是一個(gè)局部現(xiàn)象的假設(shè),使用高精度子通道程序ATHAS進(jìn)一步細(xì)分子通道以獲得相對(duì)更準(zhǔn)確的棒表面局部參數(shù)。

2) 基于Weisman汽泡壅塞模型建立了棒束CHF機(jī)理模型,提出等效圓管的概念將圓管表面湍流脈動(dòng)強(qiáng)度的分布轉(zhuǎn)換成棒束表面的強(qiáng)度分布?？紤]CHF發(fā)生位置處的汽泡壅塞造成的汽泡形變,將汽泡層厚度作為粗糙度迭代求解汽泡直徑,使其更具有物理意義。

3) 使用5×5棒束CHF實(shí)驗(yàn)數(shù)據(jù)對(duì)本文模型進(jìn)行了評(píng)價(jià)和驗(yàn)證,Rav為0.99,RMS為4.69%,預(yù)測精度優(yōu)于目前現(xiàn)有的任何一個(gè)棒束模型。本模型結(jié)合高精度子通道程序,預(yù)測導(dǎo)向管柵元和典型柵元的棒束CHF時(shí)無需進(jìn)行修正,是一個(gè)具有發(fā)展前景的棒束CHF機(jī)理模型。