熱浮力對球床反應堆內燃料球表面對流傳熱的影響

姚強 郭雪巖 楊帆

摘要：采用大渦模擬（LES ）方法研究了低雷諾數（Re =1015）下熱浮力對球床面心立方（FCC ）單元內的局部流動和傳熱的影響。為了精確求解燃料球接觸面附近的流場，球間接觸點表達為面接觸，采用結構化網格處理。研究結果表明：熱浮力的存在會抑制球床各層平均速度的波動;靠近燃料球表面處的相對時均速度和相對時均溫度受熱浮力影響變化較大;在中心流場區域，熱浮力的存在會減小流場中速度分布的不對稱性，使速度最大降低約10%，時均溫度至少升高約20%;球表面上尤其在頂部及接觸面附近受熱浮力影響明顯，努塞爾數 Nu最大降低約6%。

關鍵詞：熱浮力;球床反應堆;大渦模擬;面心立方;面接觸

中圖分類號： TK11+2??? 文獻標志碼： A

Effect of thermal buoyancy on the convective heat transfer on the surface of fuel pebbles in a pebble bed reactor

YAO Qiang，GUO Xueyan，YANG Fan

（School of Energy and Power Engineering， University of Shanghai for Science and Technology， Shanghai 200093， China）

Abstract：In this paper， large eddy simulation （LES） method was adopted to study the influence of thermal buoyancy at low Reynolds number （Re =1015） on local flow and heat transfer in a face- centered cubic （FCC） of a pebble bed. The plane-plane contact among the pebbles and structured grid were adopted to accurately solve the flow field near the interface of fuel pebbles. It is found that the existence of thermal buoyancy could restrain the fluctuation of averaged velocity at each layer. The influence of thermal buoyancy on relative time-averaged velocity and relative time-

averaged temperature near the surface of fuel pebbles was great. In the center of flow field， the existence of thermal buoyancy could reduce the asymmetry of velocity distribution. The velocity decreased by 10% at most and the time-averaged temperature increased by 20% at least. On the spherical ?surface，? especially? at? the? top? and? near? the? interface，? obvious? influence? of thermalbuoyancy was observed. And the Nusselt number Nu number decreased by 6% at most.

Keywords：thermal buoyancy; pebble bed reactor; large eddy simulation; face-centered cubic; plane-plane contact

球床式高溫氣冷堆是第四代核反應堆類型之一。在球床反應堆中，氦氣作為冷卻劑，核燃料嵌在石墨慢化劑中并可以維持很高的溫度。這種結構可以保證燃料不會熔化，使反應堆具有良好的固有安全性[1–4]。該堆型具有結構簡單、燃料元件適合批量生產以及燃料球裝卸方便等特點。

目前，很多學者對球床中的流動和換熱特點進行了研究。孟現柯等[5]用由碳鋼球堆積成球床，以蒸餾水為工質，采用電磁感應加熱方式研究了球床通道內部的換熱特性; Jia 等[6]用黑色玻璃球堆積成球床，以水為工質研究了不同床型對內部工質流動情況的影響?？紤]到實驗成本和復雜性，多數實驗結果很難為球床反應堆提供可靠的參考依據。Lee 等[7]于2007年采用k 一"湍流模型比較了不同球間間隙對流場分布的影響;2014年，Ferng等[8]用雷諾應力湍流模型、非結構網格對體心立方（ BCC ）人工間隙和點接觸兩種接觸方式的共軛模型的模擬結果進行了分析比較;2012—2015年， Shams 等[9–11]采用 q?DNS 方法、多面體網格對面心立方（FCC ）人工間隙模型進行了分析，發現在兩球之間的狹窄區域存在射流現象及速度場、溫度場分布的不對稱性;2017年，蔣旭等[12]采用大渦模擬（ LES ）方法、結構網格研究了 BCC堆積模型的流場和溫度場，通過與k一"湍流模型計算結果對比發現，采用 LES 方法能夠更好地捕捉流場中的渦結構，瞬時溫度變化更具參考價值。采用數值方法可以更容易獲得更多的流場、溫度場的信息，但多數學者采用數值方法時均只考慮球床在強制對流情況下的流動傳熱特性，而忽略了熱浮力的影響。 Huang 等[13]對矩形腔內湍流混合對流特性進行了實驗和數值研究，分析了浮力的影響，通過與實驗結果進行對比，確定了 LES 方法在混合對流計算中的精度;Guardo等[14–15]采用數值方法，對隨機球床內超臨界二氧化碳在高壓下的混合對流傳質進行了分析，研究表明，對于層流狀態下的超臨界流體，可以直觀地看到流量和流動方向對傳質量的影響，在助流狀態下，會得到更大的傳質量。研究表明，混合對流中浮力的影響不可忽略，但對球床內流動和傳熱的研究還不夠充分。

本文將采用 LES 方法、結構網格研究熱浮力對球床反應堆內 FCC 面接觸單元的燃料球表面及局部流動與傳熱的影響。

1幾何模型與數值方法

球的規則排列方式主要有3種[16]：簡單立方（ SC ），體心立方（ BCC ），面心立方（FCC ）。球床內部燃料球堆積結構可以看作是3種方式的隨機組合。兩球的接觸點一般有3種處理方式：人工間隙、點接觸、面接觸。3種規則排列方式和接觸點處理方式如圖1所示。由于3種排列方式中 FCC 結構空隙率最小，而實際球床中燃料球的無規則堆積必然使球床的空隙率最小化，因此本文選擇燃料球的 FCC 排列方式，又考慮到兩球在接觸位置會存在一定程度的變形，對接觸點的處理方式選擇了面接觸，由此得到了一個 FCC 單元。考慮到計算資源有限及為避免受進、出口效應影響，本文中采用4個 FCC 單元堆積模型，幾何模型如圖2所示。主要參數參考 Shams 等[11]中的值。燃料球球面材料為石墨，直徑為60 mm，接觸變形部分是燃料球直徑的1%，內部冷卻工質為氦氣。邊界條件及參數、物性參數分別如表1、2所示。

考慮熱浮力時，密度采用Boussinesq近似，計算 Gr/Re2= g*αv?tl/v2=0.024，其中： g*為本文中采用的重力加速度， g*= gαv（ Tmax一Tmin），Tmax和Tmin分別為;αv為體積膨脹系數;?t為進、出口溫差;l為特征長度;v為流體進口速度; Gr 為格拉曉夫數;Re 為雷諾數。復合對流中，采用Gr/Re2作為判斷自然對流影響程度的判據。一般認為，當Gr/Re2≥0.01時，自然對流的影響不能忽略。本文中重力方向與流動方向相同，熱浮力方向與流動方向相反。

為了能夠捕捉更多的瞬時流場信息，采用 LES 方法進行研究。亞格子模型采用動態Smagorinsky?Lily 模型，其優點在于具有一定的自適應性，得出的結果更可靠。

2網格適用性驗證

為了提高計算結果的準確性，對所采用的 FCC 面接觸模型進行結構網格創建，采用 ICEM 軟件將流場區域劃分成多個六面體幾何結構并建立對應的六面體塊，最后通過將幾何體與對應的六面體塊進行關聯，成功完成了結構網格的劃分，并對近壁面網格進行了加密。 FCC 面接觸單元網格如圖3所示。

為了驗證近壁面網格加密方法的適用性，對 Shams 等[12]采用的 FCC 人工間隙模型，采用相同的結構網格劃分方法進行全結構化網格的劃分。在相同的邊界條件和進、出口參數條件下，采用 LES 方法進行計算，獲取了在垂直來流方向上近球面到中心流場區域的時均速度分布。圖4為近壁面加密網格適用性驗證結果。4（a）為提取數據位置。圖4（ b）為 LES 和 q?DNS 方法結果對比，其中，橫坐標表示的位置與縱坐標表示的時均速度都進行了無量綱化，并與 Shams等[14]的 q?DNS 數據進行對比;為時均速度， U為時均速度，Umax為最大時均速度;（Z 一Zmid）為在中心流場區域所提取線上的數據點 Z 相對中間位置Zmid的坐標;D 為燃料球直徑。結果顯示，與近球面處計算結果高度一致，由此說明近壁面的網格加密方法適用于該球床模型的數值計算。

為驗證網格無關性，分別對網格數為720萬、503萬、325萬的模型進行了計算，并整理對比了進、出口速度，溫度差。以720萬網格下的計算結果為標準，計算出另外兩種網格數下得到的數據誤差。網格無關性驗證結果如表3所示，結果顯示3種網格數的計算結果無明顯差異，這說明325萬網格已可滿足計算要求，

網格數的增加對計算結果無明顯影響。由于采用了 LES 方法，為了能更好地捕捉內部流場信息，本文中用于計算的網格數取為503萬。

3計算結果與分析

在所研究的流場區域沿流動方向取相同的截面，獲取截面上的平均速度，各層截面位置及平均速度、溫度變化如圖5所示。結果顯示，在進、出口處速度變化較大，存在明顯的進、出口效應，擾動較為劇烈。進口處受影響的程度相對較大，并持續到第3層截面之后。從第4層截面開始，流速呈現對稱性波動，表明流動已基本穩定?？梢园l現，自第4層截面之后，熱浮力影響效果明顯，受熱浮力影響的流體流速相對較低，說明熱浮力的存在對流速起到抑制作用。第2層截面之后，每層溫度均有升高，相鄰兩層溫差ΔT 為7 K。為了進一步研究流場內部的流動和傳熱情況，通過以上的分析，選取第5層和第7層截面之間的 FCC 單元作為流場研究區域，選取3號中心球分析球面換熱情況。

3.1熱浮力對速度和溫度場的影響

在所研究的 FCC 單元內選取截面1、2、3，流場區域各截面位置及截面上流線圖如圖6所示。獲取截面1、3上的時均速度和時均溫度分布，結果如圖7、8所示。可以發現熱浮力對局部流場速度與溫度分布影響較為明顯。在截面1的速度場中，受熱浮力影響，低速區域增大，尤其在燃料球底部和頂部對應的尾流區和滯止區位置，低速區域明顯增大。溫度場中，低速區域的增大使高溫區域隨之增大，且在中心區域，尾流區溫度較高。截面3與接觸面相切，對比速度分布，在尾流區和滯止區同樣可以發現與截面1中相同的現象。在接觸面處熱浮力存在時，流體速度會在發生分離的位置略微增加，之后的速度都會相對減小。截面上的平均溫度相對較高，接觸面處較為明顯，高溫區域延伸較長。這再次說明熱浮力的存在對流速起到抑制作用。

為了解速度和溫度在球面和中心流場區域之間的變化情況，在截面1上中心區域沿流動方向和垂直流動方向各取一條截線，如圖9所示。提取線上的時均速度和時均溫度分布，并與不考慮熱浮力的分布進行了對比，結果如圖10所示，其中：Vb、Tb 分別為考慮熱浮力時的時均速度和時均溫度;Vn、Tn 分別為不考慮熱浮力時的時均速度和時均溫度。分析截線1上提取到的時均速度和溫度分布（相對于中心點）：沿 X 正方向，在近球面位置熱浮力對速度和溫度的影響并不明顯;遠離球面后，熱浮力的影響開始增強。在熱浮力影響下，相對中心點的對稱位置處較低的速度被提高，較高的速度被降低，最終減少了速度分布的不對稱性，也同樣減少了溫度分布的不對稱性，截線上的溫度有所提高，考慮到速度差的減小，降低了湍流強度，削弱了換熱。相對于不考慮熱浮力的情況，靠近球面位置，速度和溫度的相對變化較大（速度最大降低70%，相對層間溫差，溫度最大增加43%）。遠離球面后速度和溫度相對變化趨于平緩，在中心流場區域速度相對變化范圍在10%以內，溫度平均增加約20%。

分析截線2上提取到的時均速度分布，沿 Y 負方向（流動方向），不考慮熱浮力時，由上一個中心球底部到下一個中心球頂部，速度先快速增大，后緩慢增大，在中間偏下的位置達到峰值，后緩慢減少，至靠近球面位置急劇降低最后接近滯止，速度受滯止區和尾流區影響明顯。在考慮熱浮力后，由速度分布可以明顯發現，相對中心點呈良好的對稱性，截線上速度相對減小，峰值位置更靠近中心點。靠近球面處速度劇變的情況有所緩和，從一定程度上減少了由于速度變化大引起的湍流強度。從時均溫度分布可以發現，考慮熱浮力后，高低溫差減小，中心區域截線上平均溫度有所增加?？拷蛎嫖恢茫瑫r均速度和溫度的相對變化較大，速度幾乎滯止（相對層間溫差，溫度最大增加50%）?？拷行牧鲌鰠^域，速度和溫度的相對變化有所減小，溫度最少升高約20%。

3.2熱浮力對燃料球表面對流換熱的影響

圖11為3號中心球表面溫度分布。由圖中可見，熱浮力的存在使球面部分位置溫度相對升高，尤其在燃料球的頂部、側面的4個接觸面中間部分以及底部較為明顯。本文中定義努塞爾數 Nu = hD/λ，h = q/（T 一 Tin），其中：h為對流換熱系數;λ為氦氣熱導率;q為熱流密度;T為球面上不同位置溫度; Tin為氦氣進口溫度。球面Nu 分布如圖12所示。對比發現，熱浮力存在時，燃料球頂部及接觸面附近的 Nu 明顯降低，高 Nu 區域明顯減小，且差異較明顯的位置分布在燃料球的上半球面。由時均速度分布情況，頂部附近流場速度相對減小，減弱了對流換熱，使局部 Nu 降低。由圖6截面1、2、3上流線圖可以發現，流場中的渦多集中在接觸面后沿及燃料球底部的尾流區域，擾動較為劇烈。這說明熱浮力的存在削弱了換熱強度，而渦的存在增強了局部擾動，削弱了熱浮力的影響。

為進一步了解頂部和接觸面附近的換熱情況，在球面上3個特殊位置取3條截線，截線3、4、5，取線位置如圖13所示。球面3條截線上 Nu 分布如圖14所示。通過對比每條截線上的 Nu分布可以發現，球面頂部換熱最好，接觸面附近由于接觸面的阻礙作用 Nu較小，換熱較差。在考慮熱浮力的情況下，球面相同位置的 Nu 會相對減少，最大差值約為球面平均 Nu 的6%，即熱浮力的存在對對流換熱起到抑制作用。沿流動方向，位于球面不同位置，熱浮力影響的程度也有所不同。由圖中可以看出，燃料球的上半部分受熱浮力的影響較大，而下半部分由于流體擾動較強受熱浮力影響較小，這與分析球面 Nu 分布的結果一致。

4結論

采用 LES 方法，在考慮熱浮力的情況下，對球床反應堆 FCC 面接觸結構內的流動和傳熱進行了研究。通過對比分析整個球床各層的平均速度變化，局部流場中速度、溫度場的分布及中心燃料球表面的溫度、Nu 分布，得到熱浮力對球床反應堆局部流動和傳熱影響的幾點規律：

（1）低雷諾數下熱浮力的影響不能忽略，熱浮力的存在對球床中各層平均速度變化起到阻礙作用，會抑制速度產生較大的波動。

（2）在局部流場區域，時均速度、時均溫度受熱浮力影響明顯?？拷蛎嫖恢?，相對時均速度和相對時均溫度變化較大。在中心區域，熱浮力的存在會減小流場中速度分布的不對稱性，速度最大降低約10%，時均溫度最少升高約20%。

（3）相對來流方向，燃料球上半部分，尤其在燃料球頂部及接觸面附近，受熱浮力影響 Nu 減小，最大降低6%，削弱了換熱。流體流過接觸面接近球底部，渦結構增加，擾動劇烈，削弱了熱浮力影響。

參考文獻：

[1] 雷鳴澤.高溫氣冷堆產業推廣及應用前景[J].中國核電， 2018， 11（1）：26-29.

[2]吳宗鑫.?我國高溫氣冷堆的發展 [J]. 核動力工程， 2000，?21（1）：39 - 43，80

[3] 符曉銘， 王捷.高溫氣冷堆在我國的發展綜述[J].現代電力， 2006， 23（5）：70-75.

[4] 吳宗鑫， 肖宏才.模塊式高溫氣冷堆的安全特性[J].高技術通訊， 1994（11）：34-38.

[5] 孟現珂， 孫中寧， 徐廣展，等.含內熱源堆積球床對流換熱特性的實驗研究[J].哈爾濱工程大學學報， 2012， 33（9）：1122-1126.

[6] JIA X L， GUI N， YANG X T， et al. Experimental studyof flow field characteristics on bed configurations in thepebble bed reactor[J]. Annalsof Nuclear Energy， 2017， 102：1-10.

[7] LEE? JJ，? PARK? G? C，? KIM? K? Y，? et? al. Numericaltreatment of pebble contact in the flow and heat transferanalysis? of? a? pebble ?bed? reactor? core[J]. Nuclear Engineering and Design， 2007， 237（22）：2183-2196.

[8] FERNG Y M， LIN K Y. CFD investigation of thermal-hydraulic characteristics in a PBR core using differentcontact? treatments? between? pebbles[J]. Annals? of Nuclear Energy， 2014， 72：156-165.

[9] SHAMS? A，? ROELOFS? F，? KOMEN? E? M? J，? et? al.Optimization of a pebble bed configuration for quasi-direct numerical simulation[J]. Nuclear Engineeringand Design， 2012， 242：331-340.

[10] SHAMS A， ROELOFS F， KOMEN E M J. Quasi-directnumerical simulation of a pebble bed configuration. PartⅠ： flow （velocity） field? analysis[J]. Nuclear EngineeringandDesign， 2013， 263：473-489.

[11] SHAMS? A，? ROELOFS? F，? KOMEN? E? M? J，? et? al.Numerical? simulation? of? nuclear? pebble? bedconfigurations[J]. Nuclear Engineering and Design， 2015， 290：51-64.

[12]蔣旭， 郭雪巖.球床反應堆流動與傳熱的 CFD 分析：燃料球尺度[J].能源工程， 2017（6）：8-13，19.

[13] HUANG? YY，? YANG? G，? WU? J? Y. Large? eddysimulation and experimental study of turbulent mixedconvection inside a cavity with large Rayleigh number：effect of buoyancy[J]. BuildingandEnvironment， 2019， 151：268-279.

[14] GUARDO A， COUSSIRAT M， RECASENS F， et al.CFD study on particle-to-fluid heat transfer in fixed bedreactors： convective? heat? transfer? at? low? and? highpressure[J]. Chemical? Engineering? Science， 2006， 61（13）：4341-4353.

[15] GUARDO A， COUSSIRAT M， RECASENS F， et al.CFD studies on particle-to-fluid mass and heat transferin packed beds： free convection effects in supercriticalfluids[J]. Chemical Engineering Science， 2007， 62（18-20）：5503-5511.

[16] OOMS A. Pebble flow in a high temperature reactor[R].Physics? of Nuclear? Reactor. PNR-131-2008-003. TU Delft， The Netherlands， 2008.