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

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

顆粒堆積球床在核能領(lǐng)域中有較多的應(yīng)用，比如第4代核能系統(tǒng)中的高溫氣冷球床堆[1]、熔鹽球床堆[2]以及聚變增殖球床包層[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ù)來表征球床的綜合傳熱能力。

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

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

1 數(shù)值方法

1.1 計算模型

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

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

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

1.2 數(shù)值方法

由于球床孔隙中的氦氣是靜止的，計算中只考慮導(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。

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

(2)

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

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

表1 計算工況Table 1 Calculated case

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

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

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

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

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

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

(3)

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

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

局部有效導(dǎo)熱系數(shù)和局部固體填充率(沿著x方向等分14個區(qū)域，每個區(qū)域的填充率)分布如圖4所示，在球床中心區(qū)域(距離壁面4dp～11dp的區(qū)域)，填充率波動比較小，同時這個區(qū)域?qū)?yīng)的局部有效導(dǎo)熱系數(shù)是線性變化的，可認(rèn)為這個區(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)的影響，還會受到溫度影響[18]。當(dāng)溫度較低時，導(dǎo)熱是主要傳熱機(jī)制，隨溫度升高，石墨的導(dǎo)熱系數(shù)減小，導(dǎo)熱機(jī)制的貢獻(xiàn)減小，而輻射傳熱機(jī)制的貢獻(xiàn)增大。當(dāng)溫度超過690 K左右時，輻射逐漸成為主導(dǎo)傳熱機(jī)制[10]，因此對于工況1和2，導(dǎo)熱是主導(dǎo)傳熱機(jī)制，主體區(qū)域的導(dǎo)熱系數(shù)從高溫壁面到低溫壁面呈現(xiàn)上升趨勢；對于工況3，靠近高溫壁面的區(qū)域，輻射傳熱是主導(dǎo)傳熱機(jī)制；靠近低溫壁面的區(qū)域，導(dǎo)熱是主導(dǎo)傳熱機(jī)制；因此相對于工況1和2，工況3有效導(dǎo)熱系數(shù)的增大幅度呈減小趨勢。對于工況4～9，輻射傳熱成為主導(dǎo)傳熱機(jī)制，因此在主體區(qū)域，工況4～9有效導(dǎo)熱系數(shù)從高溫壁面到低溫壁面呈現(xiàn)下降趨勢。另外，與低溫壁面區(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)，但其隨溫度的變化趨勢與主體區(qū)域一致；相同溫度下，壁面區(qū)域有效導(dǎo)熱系數(shù)相對于主體區(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ù)也會更小；另一方面壁面區(qū)域局部固體填充率小，固體顆粒表面之間的輻射傳熱也會減弱，會進(jìn)一步減小局部有效導(dǎo)熱系數(shù)。因此，壁面區(qū)域的填充率降低會同時削弱導(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ù)的計算結(jié)果和前期改進(jìn)的ZBS模型[14]對比如圖6所示，前期改進(jìn)的ZBS模型對球床主體區(qū)域及近壁面區(qū)域的有效導(dǎo)熱系數(shù)的預(yù)測性較好，其預(yù)測相對偏差大部分在15%以內(nèi)，最大相對偏差僅在20%左右，且隨著球床溫度升高，ZBS模型的平均預(yù)測相對偏差降低；然而對于壁面區(qū)域，ZBS模型的預(yù)測值明顯偏高。這說明ZBS模型未能預(yù)測壁面區(qū)域有效導(dǎo)熱系數(shù)降低的特性。

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

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

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

keff,wall=Cwkeff,ZBS

(4)

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

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

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

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

3 結(jié)論

本文針對無序石墨球床堆有效導(dǎo)熱系數(shù)開展了數(shù)值研究，分析了無序堆積球床主體區(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)，但其隨溫度的變化趨勢與主體區(qū)域一致；由于受到局部固體填充率分布的影響，壁面區(qū)域有效導(dǎo)熱系數(shù)相對于主體區(qū)域及近壁面區(qū)域明顯降低約為22%。

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