步珊珊,陳波,田興,孫皖,馬在勇,張盧騰

(重慶大學(xué)低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,400044,重慶)

球床堆在核能領(lǐng)域中有較多的應(yīng)用,比如高溫氣冷球床堆[1]、熔鹽球床堆[2]和聚變?cè)鲋城虼舶鼘覽3]等。球床堆的有效導(dǎo)熱系數(shù)表征了球床堆芯的宏觀傳熱能力,是反應(yīng)堆熱工設(shè)計(jì)和安全分析程序中的基本參數(shù),其主要由固體表面間的輻射、球體顆粒的導(dǎo)熱和流體的導(dǎo)熱、相鄰顆粒接觸點(diǎn)之間的接觸導(dǎo)熱[4]等3個(gè)部分組成。在事故工況下,當(dāng)冷卻劑失去強(qiáng)制流動(dòng)后,堆芯冷卻劑的雷諾數(shù)很小,對(duì)流傳熱不考慮。

目前,關(guān)于有效導(dǎo)熱系數(shù)的研究主要集中于理論分析[5]和實(shí)驗(yàn)研究[6],基于孔隙尺度的數(shù)值研究還較少。孔隙尺度的數(shù)值方法可以獲得球床堆局部孔隙速度和溫度分布,通過(guò)逆求解導(dǎo)熱微分方程得到有效導(dǎo)熱系數(shù)。球床堆顆粒之間的接觸導(dǎo)熱是球床堆有效導(dǎo)熱系數(shù)中的重要組成部分,基于孔隙尺度的數(shù)值分析中,顆粒之間的接觸區(qū)域網(wǎng)格控制比較困難,通常需要對(duì)相鄰顆粒接觸區(qū)域做一定的簡(jiǎn)化[7],目前常采用間隙模型[8]、面接觸[9]和短圓柱模型[10]這3種方法簡(jiǎn)化接觸區(qū)域。因此,在數(shù)值分析中,相鄰顆粒接觸區(qū)域的簡(jiǎn)化處理會(huì)影響有效導(dǎo)熱系數(shù)的計(jì)算。

為了數(shù)值分析接觸區(qū)域的簡(jiǎn)化方法對(duì)球床堆有效導(dǎo)熱系數(shù)的影響,本文首先采用間隙模型分析了不考慮接觸導(dǎo)熱下的球床堆有效導(dǎo)熱系數(shù),然后采用面接觸和短圓柱模型分析了考慮接觸導(dǎo)熱下的球床堆有效導(dǎo)熱系數(shù),最后分析了不同溫度下接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)的影響。

1 物理模型

本文首先研究了不考慮接觸導(dǎo)熱下3種有序球床堆的有效導(dǎo)熱系數(shù),然后研究了接觸區(qū)域處理方式對(duì)簡(jiǎn)單立方結(jié)構(gòu)球床堆中有效導(dǎo)熱系數(shù)的影響。3種有序球床堆的結(jié)構(gòu)包括簡(jiǎn)單立方結(jié)構(gòu)(SC)球床堆、體心立方結(jié)構(gòu)(BCC)球床堆和面心立方結(jié)構(gòu)(FCC)球床堆。對(duì)于有序堆積結(jié)構(gòu)的球床堆而言,堆積結(jié)構(gòu)本質(zhì)上是周期性的,因此為了簡(jiǎn)化模型和減少計(jì)算量,可以采用部分幾何堆積單元結(jié)構(gòu)來(lái)研究球床堆的有效導(dǎo)熱系數(shù),如圖1所示。

(a)SC (b)BCC (c)FCC圖1 3種有序球床堆示意圖

本文中球顆粒的直徑dp為2 mm。當(dāng)不考慮接觸導(dǎo)熱時(shí),球體直徑減小為原來(lái)直徑的99%,可以看作是間隙模式;當(dāng)考慮接觸導(dǎo)熱時(shí),由于相鄰球形顆粒的接觸點(diǎn)區(qū)域網(wǎng)格生成很難控制,采用面接觸模式和短圓柱接觸模式來(lái)簡(jiǎn)化接觸點(diǎn)區(qū)域。3種不同的接觸模式如圖2所示。面接觸模式的生成方法分為兩種:方法一是球體位置不變,將球體直徑擴(kuò)大使得球體接觸;方法二是球體直徑不變,改變球體位置使得球體接觸。考慮到接觸半徑rc對(duì)傳熱的影響[11-12],又分別考慮了兩種不同接觸半徑下的面接觸模式和3種不同接觸半徑下的短圓柱接觸模式,因此本文計(jì)算模型一共有8種,如表1所示。

(a)間隙模型 (b)短圓柱接觸模式 (c)面接觸模式圖2 3種接觸模式示意圖

模型結(jié)構(gòu)dp/mm孔隙率/%接觸模式rc·d-1p/%SC1.9849.2間隙BCC1.9834.0間隙FCC1.9828.2間隙SC2.0047.6短圓柱8SC2.0047.6短圓柱2.5SC2.0047.6短圓柱1.75SC2.0246.1面接觸7.1SC2.0046.4面接觸6.3

2 數(shù)值計(jì)算方法

本文以ANSYS Fluent14.5為計(jì)算平臺(tái),基于有限容積法對(duì)能量方程進(jìn)行求解。邊界條件如下:在模型z方向(熱流傳遞方向)的前后兩個(gè)面是定溫邊界條件,前后面的溫差為10 ℃。垂直于xy平面的4個(gè)表面是絕熱的。球體材料為鈦酸鋰,孔隙中充滿氦氣,球體和氦氣的接觸面實(shí)施無(wú)滑移條件,假定熱流和溫度的連續(xù)性。求解器選擇基于壓力求解,氣體與固體接觸面設(shè)置為COUPLED,求解算法為SIMPLE。離散格式采用二階迎風(fēng)格式。固體和流體的導(dǎo)熱系數(shù)隨溫度變化[4],通過(guò)一維傅里葉導(dǎo)熱定律來(lái)計(jì)算有效導(dǎo)熱系數(shù)

(1)

式中:q是平均表面熱流量;ΔT是沿?zé)崃鞣较虻臏囟炔?。本?jì)算中流體無(wú)流動(dòng),只需要求解固體區(qū)域和流體區(qū)域的能量方程。計(jì)算工況是穩(wěn)態(tài),且球床內(nèi)無(wú)內(nèi)熱源,物性參數(shù)是溫度T的函數(shù),因此流體和固體區(qū)域的能量方程為

(2)

式中k是固體或者氣體的導(dǎo)熱系數(shù)。

采用ICEM對(duì)幾何模型進(jìn)行非結(jié)構(gòu)化網(wǎng)格劃分,接觸區(qū)域進(jìn)行適當(dāng)?shù)募用堋Ｃ娼佑|模式下固體顆粒表面、固體區(qū)域和氣體區(qū)域的網(wǎng)格如圖3所示。網(wǎng)格無(wú)關(guān)性考核如表2所示,4套網(wǎng)格單元總數(shù)分別為6 979 463、5 986 267、4 794 363和2 735 315。結(jié)果表明,當(dāng)球床平均溫度為500 ℃、網(wǎng)格數(shù)為6 979 463和5 986 267時(shí)模型計(jì)算所得的有效導(dǎo)熱系數(shù)相差很小,可以認(rèn)為網(wǎng)格總數(shù)為5 986 267時(shí)獲得網(wǎng)格無(wú)關(guān)解,其在球體或流體區(qū)域和球體接觸區(qū)

(a)顆粒表面網(wǎng)格 (b)固體區(qū)域和氣體區(qū)域網(wǎng)格圖3 面接觸模式下網(wǎng)格劃分示意圖

網(wǎng)格數(shù)平均熱流密度/W·m-1有效導(dǎo)熱系數(shù)/W·m-1·℃-1相對(duì)誤差/%6 979 4632 1470.859 05 986 2672 1460.858 60.044 794 3632 1420.856 80.252 735 3152 1400.856 10.34

的最大網(wǎng)格特征長(zhǎng)度分別小于0.1 mm和0.05 mm,其余計(jì)算模型也采用類(lèi)似網(wǎng)格進(jìn)行計(jì)算。

3 計(jì)算結(jié)果與分析

3.1 不考慮接觸導(dǎo)熱時(shí)的有效導(dǎo)熱系數(shù)分析

采用間隙模型分析不考慮接觸導(dǎo)熱時(shí)的計(jì)算結(jié)果,先將計(jì)算結(jié)果與不考慮接觸導(dǎo)熱的ZS關(guān)聯(lián)式[13]結(jié)果比較,如圖4所示。ZS關(guān)聯(lián)式如下

(3)

(4)

(5)

式中:ks是固體的導(dǎo)熱系數(shù);kg是氣體的導(dǎo)熱系數(shù);B是球床的形變參數(shù),與孔隙率有關(guān);ε是球床堆的孔隙率。

從圖4可以看出:球床堆的有效導(dǎo)熱系數(shù)隨著溫度升高而增大,3種堆積結(jié)構(gòu)下的球床堆計(jì)算結(jié)果與ZS關(guān)聯(lián)式結(jié)果的差別很小,而且趨勢(shì)走向基本一致。與ZS關(guān)聯(lián)式相比:SC球床堆的有效導(dǎo)熱系數(shù)最大誤差為1.5%;BCC球床堆的有效導(dǎo)熱系數(shù)最大誤差為3.9%;FCC球床堆的有效導(dǎo)熱系數(shù)最大誤差為2.3%。同時(shí),可以看出相同溫度下,FCC球床堆的有效導(dǎo)熱系數(shù)最大,SC球床堆的有效導(dǎo)熱系數(shù)最小。這是由于固體導(dǎo)熱系數(shù)比氣體導(dǎo)熱系數(shù)大,孔隙率越小,固體填充率越大,球床堆的有效導(dǎo)熱系數(shù)越大。綜上所述,在不考慮接觸導(dǎo)熱時(shí),間隙模型下的球床堆計(jì)算模型的計(jì)算結(jié)果能較好地預(yù)測(cè)球床堆的有效導(dǎo)熱系數(shù)。

圖4 不考慮接觸導(dǎo)熱時(shí)計(jì)算結(jié)果與ZS關(guān)聯(lián)式結(jié)果對(duì)比

3.2 考慮接觸導(dǎo)熱時(shí)的有效導(dǎo)熱系數(shù)分析

本文考慮接觸導(dǎo)熱時(shí),采用面接觸和短圓柱接觸模式來(lái)簡(jiǎn)化接觸點(diǎn)區(qū)域的計(jì)算模型。采用Kaviany關(guān)聯(lián)式[5]計(jì)算球體接觸區(qū)域?qū)嵯禂?shù)模型,與ZS關(guān)聯(lián)式疊加后可以獲得考慮接觸導(dǎo)熱下的球床堆有效導(dǎo)熱系數(shù)ZSK關(guān)聯(lián)式。Kaviany使用赫茲彈性形變?cè)?提出的計(jì)算球體接觸區(qū)域?qū)嵯禂?shù)模型如下

(6)

(7)

圖5 面接觸模式計(jì)算結(jié)果與ZSK關(guān)聯(lián)式結(jié)果對(duì)比

(a)大接觸半徑

(b)小接觸半徑圖6 短圓柱接觸模式計(jì)算結(jié)果與ZSK關(guān)聯(lián)式結(jié)果對(duì)比

3.3 不同溫度下接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)的影響

圖7 接觸導(dǎo)熱對(duì)球床堆有效導(dǎo)熱系數(shù)的貢獻(xiàn)

在不考慮輻射的情況下,不同溫度下接觸導(dǎo)熱對(duì)球床堆有效導(dǎo)熱系數(shù)的貢獻(xiàn)如圖7所示,隨著溫度的升高,接觸導(dǎo)熱占有效導(dǎo)熱系數(shù)的份額不斷降低,但是降低的趨勢(shì)逐漸變緩。接觸導(dǎo)熱對(duì)有效導(dǎo)熱系數(shù)的貢獻(xiàn)最大為13%,當(dāng)溫度達(dá)到800 ℃時(shí),接觸導(dǎo)熱的貢獻(xiàn)只占8%左右。這是由于當(dāng)流體靜止時(shí),球床堆的導(dǎo)熱系數(shù)主要受固體顆粒和氣體顆粒的導(dǎo)熱性能影響,氦氣的導(dǎo)熱系數(shù)隨溫度升高而增大;當(dāng)溫度小于550 ℃左右時(shí),固體顆粒材料Li2TiO3的導(dǎo)熱系數(shù)隨溫度的升高而減小,然后隨溫度升高而增大,如圖8所示。在高溫下,如果考慮到輻射,接觸導(dǎo)熱對(duì)有效導(dǎo)熱系數(shù)的貢獻(xiàn)還會(huì)進(jìn)一步降低,因此隨著溫度升高,接觸模式對(duì)有效導(dǎo)熱系數(shù)的影響會(huì)減弱。低溫下接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)的計(jì)算結(jié)果影響會(huì)比較顯著,高溫下接觸模式對(duì)球床堆中有效導(dǎo)熱系數(shù)計(jì)算結(jié)果的影響會(huì)顯著減小。

圖8 球床堆堆積顆粒和氣體的導(dǎo)熱系數(shù)

4 結(jié) 論

本文主要數(shù)值分析了顆粒接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)的影響,主要工作和結(jié)論如下。

(1)在不考慮接觸導(dǎo)熱時(shí),間隙模型下的球床堆計(jì)算模型的計(jì)算結(jié)果與ZS關(guān)聯(lián)式結(jié)果吻合很好。

(2)當(dāng)考慮接觸導(dǎo)熱時(shí),由于受網(wǎng)格質(zhì)量限制,面接觸計(jì)算模型的接觸半徑較大,因此計(jì)算的球床堆有效導(dǎo)熱系數(shù)與ZSK關(guān)聯(lián)式結(jié)果相比誤差較大;在網(wǎng)格生成質(zhì)量允許的情況下,推薦使用較小接觸半徑的短圓柱接觸計(jì)算模型。

(3)低溫下,接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)的計(jì)算結(jié)果影響不能忽略;高溫下,接觸模式對(duì)球床堆有效導(dǎo)熱系數(shù)計(jì)算結(jié)果的影響會(huì)顯著減小。