高溫球床壁面區(qū)域有效導(dǎo)熱系數(shù)模型優(yōu)化

步珊珊，陳 波，楊光超，陳德奇，馬在勇，張盧騰

(1.重慶大學(xué) 低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，重慶 400044；2.中廣核研究院有限公司，廣東 深圳 518000)

顆粒堆積球床在核能領(lǐng)域中有較多的應(yīng)用，比如第4代核能系統(tǒng)中的高溫氣冷球床堆[1]、熔鹽球床堆[2]以及聚變?cè)鲋城虼舶鼘覽3]等。球床內(nèi)部傳熱機(jī)制比較復(fù)雜，包括固體導(dǎo)熱-表面輻射-固體導(dǎo)熱、固體導(dǎo)熱-氣體導(dǎo)熱-固體導(dǎo)熱和固體導(dǎo)熱-接觸導(dǎo)熱-固體導(dǎo)熱等傳熱路徑[4]，因此通常用有效導(dǎo)熱系數(shù)來(lái)表征球床的綜合傳熱能力。

目前球床有效導(dǎo)熱系數(shù)模型的研究比較廣泛[5-7]，其中ZBS模型[8]考慮了接觸導(dǎo)熱以及輻射，有較為廣泛的應(yīng)用，如Enoeda等[9]測(cè)量了鈦酸鋰球床、鋯酸鋰球床以及硅酸鋰球床的有效導(dǎo)熱系數(shù)，發(fā)現(xiàn)在400～800 ℃溫度范圍內(nèi)，ZBS模型預(yù)測(cè)結(jié)果和實(shí)驗(yàn)測(cè)量數(shù)據(jù)吻合較好。Bu等[10]采用數(shù)值計(jì)算和實(shí)驗(yàn)測(cè)量研究了簡(jiǎn)單立方結(jié)構(gòu)球床有效導(dǎo)熱系數(shù)，發(fā)現(xiàn)ZBS模型和計(jì)算結(jié)果以及實(shí)驗(yàn)測(cè)量結(jié)果吻合都很好。De Beer等[11]發(fā)現(xiàn)ZBS模型可很好地預(yù)測(cè)球床中心區(qū)域的有效導(dǎo)熱系數(shù)，但壁面附近區(qū)域的預(yù)測(cè)值和實(shí)驗(yàn)測(cè)量結(jié)果偏差較大。作者前期通過(guò)實(shí)驗(yàn)結(jié)果和4種不同的有效導(dǎo)熱系數(shù)模型的對(duì)比分析也發(fā)現(xiàn)，ZBS模型可更好地預(yù)測(cè)球床主體區(qū)域的有效導(dǎo)熱系數(shù)，但在壁面附近區(qū)域的預(yù)測(cè)性能需要進(jìn)一步提高[12]。因此，高溫球床壁面附近的有效導(dǎo)熱系數(shù)特性需進(jìn)一步研究，同時(shí)ZBS模型在球床壁面區(qū)域的預(yù)測(cè)性能也需優(yōu)化。

在前期工作[10，12-14]的基礎(chǔ)上，本文針對(duì)無(wú)序石墨球床堆的有效導(dǎo)熱系數(shù)開展數(shù)值研究，分析無(wú)序堆積球床主體區(qū)域、近壁面區(qū)域及壁面區(qū)域有效導(dǎo)熱系數(shù)的分布特性，針對(duì)ZBS模型在球床壁面區(qū)域的預(yù)測(cè)性能進(jìn)行優(yōu)化，采用前期球床實(shí)驗(yàn)數(shù)據(jù)及南非HTTU實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，分析優(yōu)化后的ZBS模型在壁面區(qū)域的預(yù)測(cè)能力。

1 數(shù)值方法

1.1 計(jì)算模型

基于PFC3D軟件，采用膨脹法構(gòu)建無(wú)序球床堆結(jié)構(gòu)：在目標(biāo)孔隙率下，直徑為1/2dp(dp為目標(biāo)球徑)的球形顆粒在指定長(zhǎng)方體空間內(nèi)隨機(jī)分布，然后再將球徑擴(kuò)大至目標(biāo)球徑，通過(guò)球體顆粒間的相互擠壓和碰撞運(yùn)動(dòng)達(dá)到平衡。為消除垂直熱流方向上的壁面效應(yīng)以及減少計(jì)算量，首先構(gòu)建了15dp×15dp×14dp的無(wú)序堆積球床(圖1a)，然后在中部區(qū)域選取6dp×6dp×14dp的球床模型(圖1b)作為本文的計(jì)算模型。構(gòu)建的15dp×15dp×14dp的無(wú)序堆積球床徑向孔隙率分布如圖2所示。模型在x方向上的孔隙率分布與De Klerk[15]提出的關(guān)聯(lián)式吻合較好，平均相對(duì)誤差小于5%，驗(yàn)證了構(gòu)建的無(wú)序堆積球床結(jié)構(gòu)的可靠性。

a——無(wú)序堆積球床，15dp×15dp×14dp；b——球床計(jì)算模型，6dp×6dp×14dp圖1 球床堆計(jì)算模型Fig.1 Computational model of pebble bed reactor

圖2 球床堆孔隙率分布及驗(yàn)證Fig.2 Porosity distribution and verification of pebble bed reactor

1.2 數(shù)值方法

由于球床孔隙中的氦氣是靜止的，計(jì)算中只考慮導(dǎo)熱和輻射，因此只需要求解固體區(qū)域和流體區(qū)域的能量方程：

(1)

其中：ks為固體石墨材料的導(dǎo)熱系數(shù)，取值參考文獻(xiàn)[16]，W/(m·K)；kf為氦氣的導(dǎo)熱系數(shù)，參考NIST標(biāo)準(zhǔn)，W/(m·K)；T為溫度，K；x、y、z為坐標(biāo)，m。

本文計(jì)算基于ANSYS FLUENT 20.0軟件平臺(tái)，計(jì)算的邊界條件如圖3所示：左端壁面和右端壁面均為定溫邊界條件；另外4個(gè)壁面絕熱，因此熱流沿著x方向由高溫壁面(左端壁面)向低溫壁面(右端壁面)傳遞。顆粒表面之間的輻射換熱采用S2S(surface-to-surface)模型計(jì)算，顆粒表面從網(wǎng)格單元i離開的熱流qout,i計(jì)算如下：

(2)

圖3 邊界條件及網(wǎng)格劃分示意圖Fig.3 Boundary condition and grid system

其中：ε為發(fā)射率，本計(jì)算中壁面和球體顆粒表面均假設(shè)為漫灰體，ε均為0.8[17]；σ為Stefan-Boltzmann常數(shù)；Ti為溫度，K；ρi為反射率，0.2；Fij為網(wǎng)格單元i和j之間的角系數(shù)，角系數(shù)的計(jì)算是基于固體顆粒表面的網(wǎng)格單元；qout,j為從網(wǎng)格單元j離開的熱流。離散格式采用二階迎風(fēng)格式，當(dāng)溫度的殘差小于10-10判定為計(jì)算收斂。當(dāng)獲得球床溫度分布后，通過(guò)傅里葉導(dǎo)熱定律逆求解可獲得有效導(dǎo)熱系數(shù)。本文計(jì)算工況列于表1。

表1 計(jì)算工況Table 1 Calculated case

1.3 網(wǎng)格獨(dú)立性驗(yàn)證

采用FLUENT Meshing進(jìn)行網(wǎng)格劃分，采用多面體結(jié)合四面體網(wǎng)格的方式劃分計(jì)算域，并在球體表面以及球間短圓柱接觸區(qū)域進(jìn)行適當(dāng)加密，網(wǎng)格模型如圖3所示。通過(guò)改變加密區(qū)域的最小尺寸，獲得了3套質(zhì)量較好的網(wǎng)格系統(tǒng)，其對(duì)應(yīng)的最小網(wǎng)格尺寸分別為2.8%dp、2%dp及1%dp，對(duì)應(yīng)網(wǎng)格的數(shù)量分別為816萬(wàn)、1 317萬(wàn)及2 654萬(wàn)。將上述3套網(wǎng)格在工況4下進(jìn)行敏感性分析，計(jì)算結(jié)果列于表2，第2套網(wǎng)格與第3套的相對(duì)偏差僅為0.5%，可認(rèn)為第2套網(wǎng)格達(dá)到了獨(dú)立性要求，因此選擇該套網(wǎng)格用于數(shù)值計(jì)算分析。

表2 網(wǎng)格獨(dú)立性驗(yàn)證Table 2 Grid independence test

2 結(jié)果分析與討論

2.1 局部有效導(dǎo)熱系數(shù)分析

為計(jì)算局部有效導(dǎo)熱系數(shù)，沿?zé)崃鱾鬟f方向用13個(gè)平面將球床堆均分成14個(gè)相等的區(qū)域，每個(gè)區(qū)域的有效導(dǎo)熱系數(shù)keff計(jì)算公式如下：

(3)

其中：q為熱流密度，本計(jì)算中熱量沿x方向傳遞，每個(gè)平面的熱流密度相同，可以取高溫壁面或低溫壁面的熱流密度，W/m2；ΔT為每個(gè)區(qū)域熱流上下游平面平均溫度的差值，K；Δx為每個(gè)區(qū)域的長(zhǎng)度，1個(gè)球徑。

圖4 局部有效導(dǎo)熱系數(shù)和局部固體填充率分布Fig.4 Distribution of local effective thermal conductivity and local solid fraction

局部有效導(dǎo)熱系數(shù)和局部固體填充率(沿著x方向等分14個(gè)區(qū)域，每個(gè)區(qū)域的填充率)分布如圖4所示，在球床中心區(qū)域(距離壁面4dp～11dp的區(qū)域)，填充率波動(dòng)比較小，同時(shí)這個(gè)區(qū)域?qū)?yīng)的局部有效導(dǎo)熱系數(shù)是線性變化的，可認(rèn)為這個(gè)區(qū)域是球床主體區(qū)域。在距壁面1dp之內(nèi)，球床局部有效導(dǎo)熱系數(shù)和局部填充率的變化比較劇烈，梯度最大，因此可認(rèn)為是球床壁面區(qū)域；在距高溫壁面1dp～4dp范圍或低溫壁面1dp～3dp范圍內(nèi)，定義為近壁面區(qū)域。在中心區(qū)域，局部固體填充率基本穩(wěn)定在0.612，但局部有效導(dǎo)熱系數(shù)均勻減小，尤其是在高溫工況，這是因?yàn)榫植坑行?dǎo)熱系數(shù)分布特性不僅受到局部孔隙結(jié)構(gòu)的影響，還會(huì)受到溫度影響[18]。當(dāng)溫度較低時(shí)，導(dǎo)熱是主要傳熱機(jī)制，隨溫度升高，石墨的導(dǎo)熱系數(shù)減小，導(dǎo)熱機(jī)制的貢獻(xiàn)減小，而輻射傳熱機(jī)制的貢獻(xiàn)增大。當(dāng)溫度超過(guò)690 K左右時(shí)，輻射逐漸成為主導(dǎo)傳熱機(jī)制[10]，因此對(duì)于工況1和2，導(dǎo)熱是主導(dǎo)傳熱機(jī)制，主體區(qū)域的導(dǎo)熱系數(shù)從高溫壁面到低溫壁面呈現(xiàn)上升趨勢(shì)；對(duì)于工況3，靠近高溫壁面的區(qū)域，輻射傳熱是主導(dǎo)傳熱機(jī)制；靠近低溫壁面的區(qū)域，導(dǎo)熱是主導(dǎo)傳熱機(jī)制；因此相對(duì)于工況1和2，工況3有效導(dǎo)熱系數(shù)的增大幅度呈減小趨勢(shì)。對(duì)于工況4～9，輻射傳熱成為主導(dǎo)傳熱機(jī)制，因此在主體區(qū)域，工況4～9有效導(dǎo)熱系數(shù)從高溫壁面到低溫壁面呈現(xiàn)下降趨勢(shì)。另外，與低溫壁面區(qū)域的有效導(dǎo)熱系數(shù)相比，工況3～9高溫壁面區(qū)域的有效導(dǎo)熱系數(shù)的增大更明顯，這也是因?yàn)楦邷乇诿鎱^(qū)域的主導(dǎo)傳熱機(jī)制是輻射傳熱。

球床局部有效導(dǎo)熱系數(shù)隨溫度的變化如圖5所示，近壁面區(qū)域有效導(dǎo)熱系數(shù)存在一定的分散性特點(diǎn)，但其隨溫度的變化趨勢(shì)與主體區(qū)域一致；相同溫度下，壁面區(qū)域有效導(dǎo)熱系數(shù)相對(duì)于主體區(qū)域及近壁面區(qū)域明顯降低，最大降低為23.8%，最小降低為16.8%，平均降低約為22%。在相同溫度下，壁面區(qū)域有效導(dǎo)熱系數(shù)的減小主要是由于固體填充率較小導(dǎo)致的。高溫壁面區(qū)域和低溫壁面區(qū)域的局部固體填充率分別為0.510和0.534，與主體區(qū)域平均填充率0.612相比，分別降低了14%和20%。壁面區(qū)域局部固體填充率小，石墨導(dǎo)熱系數(shù)遠(yuǎn)大于氦氣導(dǎo)熱系數(shù)，因此局部有效導(dǎo)熱系數(shù)也會(huì)更??；另一方面壁面區(qū)域局部固體填充率小，固體顆粒表面之間的輻射傳熱也會(huì)減弱，會(huì)進(jìn)一步減小局部有效導(dǎo)熱系數(shù)。因此，壁面區(qū)域的填充率降低會(huì)同時(shí)削弱導(dǎo)熱和輻射傳熱機(jī)制，使壁面區(qū)域的有效導(dǎo)熱系數(shù)平均下降22%。

圖5 局部有效導(dǎo)熱系數(shù)隨溫度的變化Fig.5 Local effective thermal conductivity vs. temperature

2.2 壁面區(qū)域有效導(dǎo)熱系數(shù)模型優(yōu)化

有效導(dǎo)熱系數(shù)的計(jì)算結(jié)果和前期改進(jìn)的ZBS模型[14]對(duì)比如圖6所示，前期改進(jìn)的ZBS模型對(duì)球床主體區(qū)域及近壁面區(qū)域的有效導(dǎo)熱系數(shù)的預(yù)測(cè)性較好，其預(yù)測(cè)相對(duì)偏差大部分在15%以內(nèi)，最大相對(duì)偏差僅在20%左右，且隨著球床溫度升高，ZBS模型的平均預(yù)測(cè)相對(duì)偏差降低；然而對(duì)于壁面區(qū)域，ZBS模型的預(yù)測(cè)值明顯偏高。這說(shuō)明ZBS模型未能預(yù)測(cè)壁面區(qū)域有效導(dǎo)熱系數(shù)降低的特性。

圖6 有效導(dǎo)熱系數(shù)的計(jì)算結(jié)果和前期改進(jìn)的ZBS模型結(jié)果對(duì)比Fig.6 Comparison of effective thermal conductivity calculated result and previous improved ZBS model result

圖7 優(yōu)化后的ZBS模型和數(shù)值計(jì)算結(jié)果對(duì)比Fig.7 Comparison of optimized ZBS model result and numerical result

根據(jù)數(shù)值計(jì)算結(jié)果(圖5)，球床壁面區(qū)域有效導(dǎo)熱系數(shù)相比近壁面區(qū)域降低22%，因此，引入修正系數(shù)Cw對(duì)改進(jìn)的ZBS模型在球床壁面區(qū)域的預(yù)測(cè)值進(jìn)行優(yōu)化，即：

keff,wall=Cwkeff,ZBS

(4)

其中，對(duì)于球床主體區(qū)域及近壁面區(qū)域，修正系數(shù)Cw=1；對(duì)于壁面區(qū)域，修正系數(shù)Cw=0.78。優(yōu)化后ZBS模型和數(shù)值結(jié)果擬合曲線的對(duì)比如圖7所示，當(dāng)溫度小于1 600 K時(shí)，ZBS模型和計(jì)算結(jié)果在主體區(qū)域的最大相對(duì)偏差為8.7%，在壁面區(qū)域的最大相對(duì)偏差約為8.0%。需要注意的是，當(dāng)溫度大于1 600 K時(shí)，ZBS模型與數(shù)值計(jì)算結(jié)果的趨勢(shì)有所不同。因此，本文提出的修正系數(shù)適用范圍在溫度小于1 600 K，與數(shù)值計(jì)算結(jié)果的相對(duì)偏差在8%以內(nèi)。

2.3 優(yōu)化后的壁面區(qū)域模型驗(yàn)證

為驗(yàn)證優(yōu)化后的ZBS模型在壁面區(qū)域有效導(dǎo)熱系數(shù)的預(yù)測(cè)能力，采用前期無(wú)序球床實(shí)驗(yàn)結(jié)果[12]和南非HTTU實(shí)驗(yàn)數(shù)據(jù)[11]分別進(jìn)行對(duì)比，結(jié)果如圖8所示。由圖8a可見，與前期實(shí)驗(yàn)結(jié)果相比，優(yōu)化后的ZBS模型在主體區(qū)域、近壁面區(qū)域及壁面區(qū)域的有效導(dǎo)熱系數(shù)預(yù)測(cè)能力均較好，優(yōu)化后的ZBS模型可預(yù)測(cè)95%以上的實(shí)驗(yàn)數(shù)據(jù)，其相對(duì)誤差小于15%；和前期實(shí)驗(yàn)結(jié)果的最大相對(duì)誤差為20.1%。由圖8b可見，優(yōu)化后的ZBS模型在主體區(qū)域的預(yù)測(cè)值與實(shí)驗(yàn)結(jié)果吻合較好，最大相對(duì)誤差為13.1%。對(duì)于壁面及近壁面區(qū)域，ZBS模型的預(yù)測(cè)值雖仍有一定程度的偏低，但與球床有效導(dǎo)熱系數(shù)的分布規(guī)律基本吻合。上述驗(yàn)證結(jié)果表明優(yōu)化后的ZBS模型對(duì)球床壁面區(qū)域有效導(dǎo)熱系數(shù)的預(yù)測(cè)能力較好。

a——與前期實(shí)驗(yàn)結(jié)果[12]對(duì)比；b——與南非HTTU實(shí)驗(yàn)數(shù)據(jù)[11]對(duì)比(82 kW)圖8 優(yōu)化后的ZBS模型和實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.8 Comparison of optimised ZBS model and experimental data

3 結(jié)論

本文針對(duì)無(wú)序石墨球床堆有效導(dǎo)熱系數(shù)開展了數(shù)值研究，分析了無(wú)序堆積球床主體區(qū)域、近壁面區(qū)域及壁面區(qū)域有效導(dǎo)熱系數(shù)的分布特性，所得主要結(jié)論如下。

1) 在距壁面1dp之內(nèi)，球床局部有效導(dǎo)熱系數(shù)的梯度最大，被定義為球床壁面區(qū)域；在距離高溫壁面1dp～4dp范圍或低溫壁面1dp～3dp范圍內(nèi)，定義為近壁面區(qū)域。而球床主體區(qū)域則在距離壁面大于4dp～11dp的區(qū)域。

2) 近壁面區(qū)域有效導(dǎo)熱系數(shù)存在一定的分散性特點(diǎn)，但其隨溫度的變化趨勢(shì)與主體區(qū)域一致；由于受到局部固體填充率分布的影響，壁面區(qū)域有效導(dǎo)熱系數(shù)相對(duì)于主體區(qū)域及近壁面區(qū)域明顯降低約為22%。

3) 針對(duì)前期改進(jìn)后的ZBS模型，通過(guò)引入修正系數(shù)Cw來(lái)優(yōu)化其在球床壁面區(qū)域的預(yù)測(cè)能力。對(duì)于球床主體區(qū)域及近壁面區(qū)域，修正系數(shù)Cw=1，對(duì)于壁面區(qū)域，修正系數(shù)Cw=0.78。通過(guò)與前期無(wú)序球床實(shí)驗(yàn)數(shù)據(jù)及南非HTTU實(shí)驗(yàn)數(shù)據(jù)的對(duì)比，驗(yàn)證了優(yōu)化后的ZBS模型能較好地預(yù)測(cè)球床壁面區(qū)域有效導(dǎo)熱系數(shù)。