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

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

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

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

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

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

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

1.1 質(zhì)量守恒方程

(1)

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

1.2 軸向動量守恒方程

(2)

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

1.3 橫向動量守恒方程

(3)

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

1.4 能量守恒方程

(4)

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

1.5 計算流程

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

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

1.6 物性模型

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

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

(5)

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

1.7 摩擦阻力模型

程序采用達西公式計算壓降：

(6)

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

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

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

(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 湍流交混模型

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

(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)系式計算的交混系數(shù)與實驗結(jié)果對比

1.9 對流換熱模型

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

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

(19)

2 程序驗證

2.1 子通道程序流動驗證

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

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

圖5 7棒束子通道編號及劃分

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

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

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

為比較結(jié)果，將間隙橫流速度u在整個流向上進行積分，表示為坐標(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)格進行無關(guān)性驗證。圖7為網(wǎng)格計算結(jié)果與LES結(jié)果(ANL-LES)在A-A、B-B、C-C間隙處的橫向速度對比。網(wǎng)格3與網(wǎng)格4的結(jié)果相近且與LES結(jié)果吻合，證明基于大渦模擬流場數(shù)據(jù)的CFD模擬結(jié)果較好，可基于網(wǎng)格3劃分子通道進一步驗證子通道程序。

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

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

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

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

2.2 子通道程序傳熱驗證

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

圖10 19棒束的子通道編號和劃分

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

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

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

3 結(jié)論

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

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