高溫球床壁面區域有效導熱系數模型優化

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

(1.重慶大學 低品位能源利用技術及系統教育部重點實驗室，重慶 400044；2.中廣核研究院有限公司，廣東 深圳 518000)

顆粒堆積球床在核能領域中有較多的應用，比如第4代核能系統中的高溫氣冷球床堆[1]、熔鹽球床堆[2]以及聚變增殖球床包層[3]等。球床內部傳熱機制比較復雜，包括固體導熱-表面輻射-固體導熱、固體導熱-氣體導熱-固體導熱和固體導熱-接觸導熱-固體導熱等傳熱路徑[4]，因此通常用有效導熱系數來表征球床的綜合傳熱能力。

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

在前期工作[10，12-14]的基礎上，本文針對無序石墨球床堆的有效導熱系數開展數值研究，分析無序堆積球床主體區域、近壁面區域及壁面區域有效導熱系數的分布特性，針對ZBS模型在球床壁面區域的預測性能進行優化，采用前期球床實驗數據及南非HTTU實驗數據進行對比，分析優化后的ZBS模型在壁面區域的預測能力。

1 數值方法

1.1 計算模型

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

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

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

1.2 數值方法

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

(1)

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

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

(2)

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

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

表1 計算工況Table 1 Calculated case

1.3 網格獨立性驗證

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

表2 網格獨立性驗證Table 2 Grid independence test

2 結果分析與討論

2.1 局部有效導熱系數分析

為計算局部有效導熱系數，沿熱流傳遞方向用13個平面將球床堆均分成14個相等的區域，每個區域的有效導熱系數keff計算公式如下：

(3)

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

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

局部有效導熱系數和局部固體填充率(沿著x方向等分14個區域，每個區域的填充率)分布如圖4所示，在球床中心區域(距離壁面4dp～11dp的區域)，填充率波動比較小，同時這個區域對應的局部有效導熱系數是線性變化的，可認為這個區域是球床主體區域。在距壁面1dp之內，球床局部有效導熱系數和局部填充率的變化比較劇烈，梯度最大，因此可認為是球床壁面區域；在距高溫壁面1dp～4dp范圍或低溫壁面1dp～3dp范圍內，定義為近壁面區域。在中心區域，局部固體填充率基本穩定在0.612，但局部有效導熱系數均勻減小，尤其是在高溫工況，這是因為局部有效導熱系數分布特性不僅受到局部孔隙結構的影響，還會受到溫度影響[18]。當溫度較低時，導熱是主要傳熱機制，隨溫度升高，石墨的導熱系數減小，導熱機制的貢獻減小，而輻射傳熱機制的貢獻增大。當溫度超過690 K左右時，輻射逐漸成為主導傳熱機制[10]，因此對于工況1和2，導熱是主導傳熱機制，主體區域的導熱系數從高溫壁面到低溫壁面呈現上升趨勢；對于工況3，靠近高溫壁面的區域，輻射傳熱是主導傳熱機制；靠近低溫壁面的區域，導熱是主導傳熱機制；因此相對于工況1和2，工況3有效導熱系數的增大幅度呈減小趨勢。對于工況4～9，輻射傳熱成為主導傳熱機制，因此在主體區域，工況4～9有效導熱系數從高溫壁面到低溫壁面呈現下降趨勢。另外，與低溫壁面區域的有效導熱系數相比，工況3～9高溫壁面區域的有效導熱系數的增大更明顯，這也是因為高溫壁面區域的主導傳熱機制是輻射傳熱。

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

圖5 局部有效導熱系數隨溫度的變化Fig.5 Local effective thermal conductivity vs. temperature

2.2 壁面區域有效導熱系數模型優化

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

圖6 有效導熱系數的計算結果和前期改進的ZBS模型結果對比Fig.6 Comparison of effective thermal conductivity calculated result and previous improved ZBS model result

圖7 優化后的ZBS模型和數值計算結果對比Fig.7 Comparison of optimized ZBS model result and numerical result

根據數值計算結果(圖5)，球床壁面區域有效導熱系數相比近壁面區域降低22%，因此，引入修正系數Cw對改進的ZBS模型在球床壁面區域的預測值進行優化，即：

keff,wall=Cwkeff,ZBS

(4)

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

2.3 優化后的壁面區域模型驗證

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

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

3 結論

本文針對無序石墨球床堆有效導熱系數開展了數值研究，分析了無序堆積球床主體區域、近壁面區域及壁面區域有效導熱系數的分布特性，所得主要結論如下。

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

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

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