鉛鉍冷卻快堆含繞絲燃料組件子通道程序開發(fā)與驗(yàn)證

劉佳泰，彭天驥,3，蘇興康，顧 龍,4,*

(1.中國(guó)科學(xué)院 近代物理研究所，甘肅 蘭州 730000；2.中國(guó)科學(xué)院大學(xué)，北京 100049；3.先進(jìn)能源科學(xué)與技術(shù)廣東省實(shí)驗(yàn)室，廣東 惠州 516003；4.蘭州大學(xué)，甘肅 蘭州 730000)

中國(guó)加速器驅(qū)動(dòng)嬗變研究裝置(CiADS)由直線加速器、散裂靶和次臨界反應(yīng)堆組成[1]，其中次臨界反應(yīng)堆選用液態(tài)鉛鉍(LBE)冷卻快堆[2]。鉛鉍冷卻快堆具有較高的安全性和可持續(xù)性等優(yōu)點(diǎn)，然而受限于LBE非線性的湍流普朗特?cái)?shù)以及LBE具有腐蝕性等困難，通過(guò)實(shí)驗(yàn)獲得堆芯的精細(xì)熱工水力參數(shù)所需要滿足的條件十分嚴(yán)格。CFD的幾何刻畫精細(xì)，但對(duì)整個(gè)堆芯或燃料組件進(jìn)行網(wǎng)格劃分和數(shù)值計(jì)算需要大量的計(jì)算資源，計(jì)算時(shí)間較長(zhǎng)。子通道分析方法通過(guò)將燃料組件劃分為內(nèi)通道、邊通道和角通道，將經(jīng)驗(yàn)關(guān)系式代入單個(gè)子通道的守恒方程中使方程組封閉，迭代求解獲得一定精度的熱工水力參數(shù)，可大幅減少精細(xì)計(jì)算時(shí)間。

對(duì)于液態(tài)金屬，國(guó)內(nèi)外開發(fā)了許多適用于含繞絲燃料組件的子通道程序。早期多是基于鈉冷快堆，如ASFRE-Ⅲ[3]基于分布阻力模型對(duì)鈉冷快堆進(jìn)行計(jì)算；MATRA-LMR[4]是基于壓水堆程序MATRA開發(fā)的鈉冷快堆子通道程序；COBRA-WC[5]和COBRA-LM[6]均是基于美國(guó)西北太平洋實(shí)驗(yàn)室開發(fā)的COBRA程序改進(jìn)成的鈉冷快堆程序。近年來(lái)針對(duì)鉛鉍冷卻快堆的研究日益增多，如SACOS-PB[7]可計(jì)算鉛冷快堆穩(wěn)態(tài)條件下的溫度分布；ANTEO+[8]是一種適用于液態(tài)金屬?gòu)?qiáng)制對(duì)流的通用程序。目前已有的液態(tài)金屬子通道程序由于使用的關(guān)系式不同，具有不同的適用條件和限制，是否適用于CiADS的子通道分析需進(jìn)一步驗(yàn)證。

本文開發(fā)適用于鉛鉍冷卻快堆的子通道程序，對(duì)液態(tài)鉛鉍的摩擦阻力模型、湍流交混模型和對(duì)流換熱模型進(jìn)行適用性分析，并與含繞絲燃料組件的LBE大渦模擬(LES)計(jì)算結(jié)果和傳熱實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證。

1 守恒方程及本構(gòu)關(guān)系

本研究以鉛鉍冷卻快堆含繞絲燃料組件為研究對(duì)象，采用子通道分析程序?qū)Ψ€(wěn)態(tài)強(qiáng)迫對(duì)流工況下冷卻劑及包殼表面的溫度分布進(jìn)行計(jì)算。程序基于Fortran平臺(tái)進(jìn)行開發(fā)，由于LBE在正常工況下不會(huì)產(chǎn)生沸騰，因此采用單相流狀態(tài)下的守恒方程。

1.1 質(zhì)量守恒方程

(1)

式中：A為軸向流動(dòng)面積；ρ為冷卻劑密度；t為時(shí)間；m為軸向質(zhì)量流量；z為軸向控制體高度；eik為符號(hào)函數(shù)，i為子通道的編號(hào)，k為子通道i周邊的間隙；w為單位長(zhǎng)度橫向質(zhì)量流量。方程各項(xiàng)分別為質(zhì)量隨時(shí)間的變化、軸向質(zhì)量流量隨空間的變化和間隙處的橫向質(zhì)量流量。

1.2 軸向動(dòng)量守恒方程

(2)

式中：U為軸向流速；p為壓力；g為重力加速度；θ為冷卻劑通道與垂直方向的夾角；f為摩擦阻力系數(shù)；D為燃料棒直徑；K為形阻系數(shù)；fT為湍流動(dòng)量交換系數(shù)；w′為單位長(zhǎng)度湍流交混流量。方程左邊3項(xiàng)分別為軸向動(dòng)量隨時(shí)間的變化、軸向動(dòng)量通量隨空間的變化和橫流導(dǎo)致的橫向動(dòng)量變化。方程右邊分別為壓力項(xiàng)、重力項(xiàng)、摩擦力項(xiàng)和湍流交混項(xiàng)。方程右邊各項(xiàng)共同形成軸向力。

1.3 橫向動(dòng)量守恒方程

(3)

式中：s為間隙寬度；l為子通道i與子通道k之間的質(zhì)心距離；KG為橫向阻力系數(shù)；Δp為相鄰子通道間的橫向壓差。方程左邊分別為橫向流量隨時(shí)間的變化和橫向流量通量隨空間的變化，方程右邊分別為橫向壓差項(xiàng)和間隙阻力項(xiàng)。由于繞絲結(jié)構(gòu)產(chǎn)生橫流，使燃料組件的同一橫截面內(nèi)具有明顯的壓力差。

1.4 能量守恒方程

(4)

式中：h為流體焓值；PW為燃料棒r面向子通道i的加熱周長(zhǎng)；φir為燃料棒r面向子通道i的份額；q為燃料棒的線功率密度；CQ為冷卻劑產(chǎn)熱份額；φin為燃料棒與通道之間的接觸份額；q′為燃料棒傳入冷卻劑的線功率密度；Δh為相鄰?fù)ǖ篱g焓差；Ck為橫向換熱系數(shù)；ΔT為相鄰?fù)ǖ篱g的溫度差。方程左邊分別為焓隨時(shí)間的變化、軸向焓通量隨空間的變化和通道內(nèi)所有間隙處的橫向焓通量。方程右邊分別為燃料棒傳入流體的熱量、直接在冷卻劑中產(chǎn)生的熱量、湍流交混的熱量和相鄰?fù)ǖ篱g的熱量交換。

1.5 計(jì)算流程

子通道程序在軸向節(jié)點(diǎn)的入口沿流向逐步計(jì)算。首先通過(guò)能量方程更新焓值，之后對(duì)質(zhì)量方程和動(dòng)量方程聯(lián)立的矩陣采用高斯-賽德爾迭代計(jì)算得到橫向流量。最后將橫向流量帶入到質(zhì)量方程計(jì)算得到軸向流量，帶入到橫向動(dòng)量方程計(jì)算得到壓力梯度。程序求解流程如圖1所示，主程序在讀取邊界條件和幾何條件后，依次調(diào)用物性模塊、能量方程迭代模塊和動(dòng)量方程迭代模塊完成循環(huán)，直到軸向節(jié)點(diǎn)全部計(jì)算收斂后輸出計(jì)算結(jié)果。

圖1 子通道程序中冷卻劑求解流程

1.6 物性模型

對(duì)程序模型進(jìn)行不確定性分析是評(píng)估和驗(yàn)證模型準(zhǔn)確性的關(guān)鍵環(huán)節(jié)。冷卻劑的物性作為模型關(guān)鍵輸入?yún)?shù)，在不確定性的傳播過(guò)程中起到關(guān)鍵作用。冷卻劑相關(guān)物性的不確定度列于表1[9]，包括密度、比定壓熱容、動(dòng)力黏度及熱導(dǎo)率。程序中還需已知冷卻劑的焓值，通過(guò)焓的定義計(jì)算得到。

表1 LBE的推薦物性關(guān)系式

(5)

式中：h0為熔點(diǎn)溫度處的焓值；T為冷卻劑溫度；TM為熔點(diǎn)；cp為比定壓熱容。

1.7 摩擦阻力模型

程序采用達(dá)西公式計(jì)算壓降：

(6)

式中：L為通道長(zhǎng)度；De為通道水力直徑；v為冷卻劑流速。

求解壓降的關(guān)鍵是得到合適的摩擦阻力系數(shù)f，而摩擦阻力系數(shù)一般是雷諾數(shù)Re和幾何關(guān)系的函數(shù)。Fan等[10]以CiADS棒束為原型，開展了阻力系數(shù)測(cè)量實(shí)驗(yàn)，并將實(shí)驗(yàn)結(jié)果與摩擦關(guān)系式進(jìn)行對(duì)比。如圖2所示，在3 750≤Re≤16 250范圍內(nèi)，Rehme關(guān)系式[11]能在±10%的范圍對(duì)全部54個(gè)實(shí)驗(yàn)點(diǎn)準(zhǔn)確預(yù)測(cè)，均方根誤差為7.8%。因此本程序選用Rehme關(guān)系式計(jì)算摩擦阻力系數(shù)：

圖2 Rehme關(guān)系式與CiADS實(shí)驗(yàn)的摩擦阻力系數(shù)對(duì)比

(7)

式中：Pwb為棒束濕周；Pwt為總濕周；F為幾何因子，用于考慮繞絲的幾何結(jié)構(gòu)。

(8)

式中：P為相鄰燃料棒質(zhì)心的距離；Dw為繞絲直徑；H為繞絲螺距。關(guān)系式的適用范圍為：103≤Re≤3×105，1.125≤P/D≤1.417，5≤H/D≤50。

1.8 湍流交混模型

繞絲的存在會(huì)增加相鄰子通道間的橫向交混，具有展平橫向溫度分布的作用。假設(shè)湍流交混不引起質(zhì)量交換，只引起動(dòng)量和能量上的交換。考慮湍流交混時(shí)的橫向流量，采用以下公式：

(9)

內(nèi)部通道為：

(10)

(11)

(12)

壁面處通道為：

C1L,Lam=0.413(H/D)-0.5(Ar2/A′2)0.5tanθ

(13)

C1L,Tur=0.73(H/D)-0.5(Ar2/A′2)0.5tanθ

(14)

(15)

其中：

ψb=lg(Reb/RebL)/lg(RebT/RebL)

(16)

(17)

(18)

圖3 湍流交混關(guān)系式計(jì)算的交混系數(shù)與實(shí)驗(yàn)結(jié)果對(duì)比

1.9 對(duì)流換熱模型

對(duì)于包殼外表面和冷卻劑之間的對(duì)流換熱，需要通過(guò)努塞爾數(shù)(Nu)來(lái)確定對(duì)流換熱系數(shù)。對(duì)于LBE，Nu通常被認(rèn)為是棒徑比和貝克萊數(shù)(Pe)的函數(shù)。Pacio等[13]在LBE傳熱實(shí)驗(yàn)中選取了3個(gè)加熱段測(cè)量位置ML1(z/H=1/6)、ML2(z/H=11/6)和ML3(z/H=15/6)，并與液態(tài)金屬的對(duì)流換熱關(guān)系式進(jìn)行對(duì)比，最終推薦Kazimi等[14]關(guān)系式用于計(jì)算液態(tài)鉛鉍的努塞爾數(shù)，如式(19)所示。圖4給出經(jīng)驗(yàn)關(guān)系式與實(shí)驗(yàn)測(cè)得努塞爾數(shù)之間的對(duì)比，可看出關(guān)系式低估了ML1處的Nu，這是因?yàn)镸L1在熱發(fā)展區(qū)，但在位于充分發(fā)展區(qū)的ML2和ML3時(shí)對(duì)Nu的預(yù)測(cè)比較準(zhǔn)確，均方根誤差控制在7.1%。因此程序采用Kazimi等關(guān)系式計(jì)算努塞爾數(shù)，適用范圍為：10≤Pe≤5 000，1.1≤P/D≤1.40。

圖4 對(duì)流換熱關(guān)系式計(jì)算的努塞爾數(shù)與實(shí)驗(yàn)結(jié)果對(duì)比

(19)

2 程序驗(yàn)證

2.1 子通道程序流動(dòng)驗(yàn)證

采用CFD可獲得精細(xì)的流場(chǎng)數(shù)據(jù)，為驗(yàn)證CFD結(jié)果可靠，采用大渦模擬計(jì)算結(jié)果作為基準(zhǔn)驗(yàn)證[15]，該結(jié)果為美國(guó)阿貢國(guó)家實(shí)驗(yàn)室通過(guò)開源程序Nek5000計(jì)算得到，具有高保真度。本研究在此模型基礎(chǔ)上進(jìn)一步劃分子通道供程序驗(yàn)證使用。

棒束結(jié)構(gòu)及子通道編號(hào)如圖5所示，繞絲沿流向逆時(shí)針?lè)较蛐D(zhuǎn)。表2列出棒束的具體幾何參數(shù)。幾何建模時(shí)，由于繞絲與燃料棒之間為線接觸，不利于網(wǎng)格劃分，因此在每個(gè)繞絲與燃料棒接觸處做0.25 mm的倒圓角處理，以提高網(wǎng)格質(zhì)量。

圖5 7棒束子通道編號(hào)及劃分

表2 7棒束幾何參數(shù)

采用商業(yè)計(jì)算流體軟件STARCCM+進(jìn)行多面體網(wǎng)格劃分。圖6示出整體網(wǎng)格及間隙處的局部網(wǎng)格細(xì)節(jié)。在近壁面處設(shè)置4層邊界層網(wǎng)格，選取一個(gè)螺距高度進(jìn)行計(jì)算，進(jìn)出口設(shè)置周期性邊界，進(jìn)口雷諾數(shù)為9 457，生成網(wǎng)格后，使用SSTk-ω湍流模型進(jìn)行計(jì)算。

圖6 7棒束網(wǎng)格劃分上視圖

為比較結(jié)果，將間隙橫流速度u在整個(gè)流向上進(jìn)行積分，表示為坐標(biāo)s的函數(shù)：

(20)

分別選取0.8 mm(網(wǎng)格1)、0.6 mm(網(wǎng)格2)、0.3 mm(網(wǎng)格3)和0.2 mm(網(wǎng)格4) 4套不同尺寸的網(wǎng)格進(jìn)行無(wú)關(guān)性驗(yàn)證。圖7為網(wǎng)格計(jì)算結(jié)果與LES結(jié)果(ANL-LES)在A-A、B-B、C-C間隙處的橫向速度對(duì)比。網(wǎng)格3與網(wǎng)格4的結(jié)果相近且與LES結(jié)果吻合，證明基于大渦模擬流場(chǎng)數(shù)據(jù)的CFD模擬結(jié)果較好，可基于網(wǎng)格3劃分子通道進(jìn)一步驗(yàn)證子通道程序。

a——間隙A-A；b——間隙B-B；c——間隙C-C

為驗(yàn)證程序流動(dòng)計(jì)算的準(zhǔn)確性，將子通道程序的計(jì)算結(jié)果與CFD結(jié)果進(jìn)行對(duì)比。圖8為出口處不同子通道的質(zhì)量流量分布，邊通道的出口流量最大，內(nèi)通道其次，角通道出口流量最小，程序能準(zhǔn)確計(jì)算出口質(zhì)量流量。圖9為子通道1的冷卻劑軸向速度沿流向在1個(gè)螺距長(zhǎng)度的分布，可看出子通道程序計(jì)算的軸向速度與CFD計(jì)算結(jié)果具有相同的趨勢(shì)。

圖8 出口處子通道質(zhì)量流量分布

圖9 子通道1軸向速度分布

2.2 子通道程序傳熱驗(yàn)證

子通道程序要求能準(zhǔn)確計(jì)算冷卻劑溫度和燃料棒包殼外表面溫度。為驗(yàn)證程序傳熱計(jì)算的準(zhǔn)確性，選取Pacio等[13]開展的19棒束含繞絲燃料組件的棒束傳熱實(shí)驗(yàn)進(jìn)行驗(yàn)證。燃料組件的具體幾何參數(shù)列于表3。具體邊界條件包括進(jìn)口溫度Tin=473 K、進(jìn)口質(zhì)量流量M=19.18 kg/s及總加熱功率Q=197.04 kW，加熱功率均勻分布在每根棒的加熱段上。在程序中模擬824 mm的流動(dòng)發(fā)展段和870 mm的加熱段，棒束子通道劃分與編號(hào)方式如圖10所示。

圖10 19棒束的子通道編號(hào)和劃分

包殼和冷卻劑溫度計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比如圖11～13所示，分別對(duì)應(yīng)沿加熱段軸向高度的3個(gè)不同測(cè)量位置ML1、ML2和ML3。其中ML1處于熱發(fā)展段，ML2和ML3處于熱平衡段。除實(shí)驗(yàn)數(shù)據(jù)外，另外選取子通道程序SACOS-PB的計(jì)算結(jié)果進(jìn)行對(duì)比[16]。由圖12、13可看出，內(nèi)通道比邊通道和角通道更熱，ML3比ML2的溫度曲線變化更明顯。同一截面上，邊通道和角通道的溫度均低于流體平均溫度。程序計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)在子通道處的最大相對(duì)誤差為4.08%，包殼外表面處的最大相對(duì)誤差為4.36%，且與其他子通道程序的計(jì)算精度類似，總體計(jì)算相對(duì)誤差低于5%，驗(yàn)證結(jié)果較好。

圖11 ML1處包殼外表面和冷卻劑溫度分布

圖12 ML2處包殼外表面和冷卻劑溫度分布

3 結(jié)論

本文開發(fā)了適用于液態(tài)鉛鉍冷卻含繞絲燃料組件的子通道分析程序，并對(duì)不同流動(dòng)和傳熱關(guān)系式進(jìn)行了評(píng)價(jià)?；?棒束大渦模擬流場(chǎng)數(shù)據(jù)驗(yàn)證了程序流動(dòng)模塊的有效性，基于19棒束傳熱實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了程序傳熱模塊的有效性。結(jié)果表明開發(fā)的子通道程序能較好地計(jì)算液態(tài)鉛鉍冷卻繞絲定位燃料組件的流動(dòng)傳熱特性，可代替或輔助CFD計(jì)算，節(jié)約了計(jì)算時(shí)間。

圖13 ML3處包殼外表面和冷卻劑溫度分布