姚強 郭雪巖 楊帆

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

關(guān)鍵詞：熱浮力;球床反應(yīng)堆;大渦模擬;面心立方;面接觸

中圖分類號： TK11+2??? 文獻標(biāo)志碼： 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

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

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

本文將采用 LES 方法、結(jié)構(gòu)網(wǎng)格研究熱浮力對球床反應(yīng)堆內(nèi) FCC 面接觸單元的燃料球表面及局部流動與傳熱的影響。

1幾何模型與數(shù)值方法

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

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

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

2網(wǎng)格適用性驗證

為了提高計算結(jié)果的準(zhǔn)確性，對所采用的 FCC 面接觸模型進行結(jié)構(gòu)網(wǎng)格創(chuàng)建，采用 ICEM 軟件將流場區(qū)域劃分成多個六面體幾何結(jié)構(gòu)并建立對應(yīng)的六面體塊，最后通過將幾何體與對應(yīng)的六面體塊進行關(guān)聯(lián)，成功完成了結(jié)構(gòu)網(wǎng)格的劃分，并對近壁面網(wǎng)格進行了加密。 FCC 面接觸單元網(wǎng)格如圖3所示。

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

為驗證網(wǎng)格無關(guān)性，分別對網(wǎng)格數(shù)為720萬、503萬、325萬的模型進行了計算，并整理對比了進、出口速度，溫度差。以720萬網(wǎng)格下的計算結(jié)果為標(biāo)準(zhǔn)，計算出另外兩種網(wǎng)格數(shù)下得到的數(shù)據(jù)誤差。網(wǎng)格無關(guān)性驗證結(jié)果如表3所示，結(jié)果顯示3種網(wǎng)格數(shù)的計算結(jié)果無明顯差異，這說明325萬網(wǎng)格已可滿足計算要求，

網(wǎng)格數(shù)的增加對計算結(jié)果無明顯影響。由于采用了 LES 方法，為了能更好地捕捉內(nèi)部流場信息，本文中用于計算的網(wǎng)格數(shù)取為503萬。

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

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

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

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

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

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

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

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

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

4結(jié)論

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

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

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

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

參考文獻：

[1] 雷鳴澤.高溫氣冷堆產(chǎn)業(yè)推廣及應(yīng)用前景[J].中國核電， 2018， 11（1）：26-29.

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

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

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

[5] 孟現(xiàn)珂， 孫中寧， 徐廣展，等.含內(nèi)熱源堆積球床對流換熱特性的實驗研究[J].哈爾濱工程大學(xué)學(xué)報， 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]蔣旭， 郭雪巖.球床反應(yīng)堆流動與傳熱的 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.