姚強(qiáng) 郭雪巖 楊帆

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

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

中圖分類號(hào)： TK11+2??? 文獻(xiàn)標(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)可以保證燃料不會(huì)熔化，使反應(yīng)堆具有良好的固有安全性[1–4]。該堆型具有結(jié)構(gòu)簡(jiǎn)單、燃料元件適合批量生產(chǎn)以及燃料球裝卸方便等特點(diǎn)。

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1熱浮力對(duì)速度和溫度場(chǎng)的影響

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

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

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

3.2熱浮力對(duì)燃料球表面對(duì)流換熱的影響

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

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

4結(jié)論

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

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

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

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

參考文獻(xiàn)：

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

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

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

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

[5] 孟現(xiàn)珂， 孫中寧， 徐廣展，等.含內(nèi)熱源堆積球床對(duì)流換熱特性的實(shí)驗(yàn)研究[J].哈爾濱工程大學(xué)學(xué)報(bào)， 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)堆流動(dò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.