球床式高溫氣冷堆球流及球床輻射的理論與實驗研究進展

姜勝耀，桂 南，楊星團，屠基元,2

(1.清華大學 核能與新能源技術研究院，北京 100084； 2.墨爾本皇家理工大學，墨爾本 VIC 3083，澳大利亞)

高溫氣冷堆滿足第4代先進反應堆系統的各種要求，是21世紀美國甚至全世界最有前途的堆型[1]。我國在2006年將高溫氣冷堆項目列入《國家中長期科學和技術發展規劃綱要(2006—2020)》，是16個國家科技重大專項之一。高溫氣冷堆的發展路線主要有球床式和柱狀式，柱狀式氣冷堆主要包括美國與俄羅斯合作的GT-MHR(Gas Turbine-Modular Helium Reactor)[2-3]和日本HTTR(High Temperature Engineering Test Reactor)[4-5]等。球床式高溫氣冷堆主要包括美國MPBR(Modular Pebble Bed Reactor)[6]、南非PBMR(Pebble Bed Modular Reactor)[7]以及中國10 MW高溫氣冷堆(HTR-10)[8-10]及模塊化球床高溫氣冷堆商用核電站示范工程(HTR-PM)[8,11]。

球床式高溫氣冷堆最突出的優點之一就是其固有安全性。球床式高溫氣冷堆的燃料元件是將全陶瓷型包覆顆粒彌散在石墨球基體中制成的，在1 600 ℃以下能保持完整性[10]。由清華大學核能與新能源技術研究院研發的HTR-10在2000年實現臨界。目前HTR-PM[12]正在建設調試中，預計于2020年實現并網。

球床堆中隨機密集堆積的燃料球在重力驅動下極緩慢地流過堆芯，形成了一種極緩慢的特殊顆粒流，稱之為球流。球流運動特性是球床式堆芯幾何設計、物理設計及核燃料循環等的基礎，在球床式高溫氣冷堆設計改進、熱工安全分析中具有基礎性的地位。德國、美國、南非、中國等均積極開展了相關研究[13-15]，包括實驗的和理論的，但研究尚不充分。

另一方面，堆芯球床中的熱傳導包括球體材料的導熱、球球間接觸導熱、球表面的輻射傳熱、球與氣體間的對流傳熱等多種傳熱機制，工程上可用一個綜合的等效導熱系數將其概括為簡單的傅里葉導熱[16]。因此，堆芯球床等效導熱系數是綜合反映球床堆傳熱能力的宏觀參數。雖然德國研究中心安全分析和反應堆技術研究所進行了不同工況的SANA-1實驗，給出了球床等效導熱系數最早的實驗數據，但其實驗給出的球床等效導熱系數溫度范圍在1 000 ℃以下，未能覆蓋高溫堆失冷失壓事故工況極限溫度范圍(1 640 ℃)，而高溫堆堆芯局部溫度可達1 600 ℃，高溫下顆粒間的輻射換熱主導著整個球床傳熱，該溫度范圍內的數據恰是最需要的[17]。南非西北大學和M-Tech公司聯合為PBMR項目建立了用于研究堆芯球床內的熱傳導及流體流動現象的球床導熱實驗裝置(HTTU)，但1 200 ℃以上時，實驗曲線向下掉，失真度較高[18]。因此，鑒于堆芯球床等效導熱系數是直接影響高溫氣冷堆燃料最高溫度和堆芯溫度分布的關鍵參數，在余熱導出過程中起主導作用，開展球床等效導熱系數的實驗研究對于反應堆分析程序的完善、研究提高高溫氣冷堆單堆功率的可能性以及工程的安全分析具有重要意義[16]。

本文對近年來的球流研究進行總結，同時給出下一階段球流研究的方向。

1 實驗研究

1.1 球流實驗

1) 唯象實驗研究

為更好地掌握球流運動特征，根據球床式高溫氣冷堆堆芯幾何設計，建造了縮比的球流實驗裝置[19-21](圖1a)。球流唯象實驗分析成功獲得了球停留時間、滯留區及滯留特性[22]、球流流線及球流紡(在獲得足夠多顆粒運動軌跡信息的基礎上得到約95%軌跡線所包絡的區域稱為球流紡)、雙區堆芯及球交混及擴散特性[23]的定性規律，并給后續的定量分析、機理研究和模擬對比驗證提供了基礎數據。基于球流唯象實驗，提出了標志球法[24]來分析球流流場，標志球法是通過直接測量標志球(彩色球)在特定時刻的空間位置獲得其運動軌跡，在獲得大量標志球(圖1b中id0～id8)的運動軌跡的基礎上得到球的坐標和時間信息并計算得到球床內顆粒球流動的定量結果(圖1b)。

2) 定量化實驗研究

首先發展了二維球流運動界面分析法(SAM)[25]，在控制加球和卸球方式的工況下達到穩定分區球流態(即卸球中不同顏色球所占的比例與加球中相應顏色球占比相同；均值流線的曲面的球流通量處處為0，通過均值線面的瞬時球流通量一段時間內累積為0)。達到穩定分區時，可得到不同加球比例下的均值流線，并可進一步計算不同高度處的中心區寬度，給出得到的平均速度的分布等。

圖1 球流實驗裝置(a)和標志球法得到的球流流線(b)Fig.1 Pebble flow setup (a) and measured pebble streamlines (b)

灰點表示單顆粒速度，綠點代表局部平均速度，R=0處顯示該處的標準差圖3 不同高度處的VyFig.3 Vertical velocity Vy on different heights

進而對球流進行了顆粒速度測量學(PTV)分析[26]，采用了動態灰度閾值方法(DTB)和新型匹配幾率法(RM法)，并經顆粒識別方法檢驗[27]。實驗用黑色玻璃球，循環球速率為150 min-1，每秒卸球的個數約占總裝球量的0.003 5%，按平均值計算顆粒球每秒移動的位移小于球直徑的0.01倍。利用高速攝影儀進行連續拍照，對圖片進行預處理、顆粒識別、顆粒匹配及數據分析，如圖2所示。各工況重復實驗3～5次用于檢驗實驗的可重復性，以保證統計可靠性。

通過PTV技術對球速度定量分析，可得不同高度條件下平均垂直速度的水平方向定量分布(圖3)。圖3中，Vy為垂直方向的顆粒速度，R為徑向坐標(以顆粒直徑d為單位)，h為數據提取點所在高度。局部平均速度通過灰色散點表示并呈現高斯分布特征。

圖2 實驗用顆粒球(a)、顆粒球位置標記(b) 及顆粒球軌跡(c)Fig.2 Experimental images of particles (a),particles marked by red crosses (b) and trajectories of particles (c)

基于PTV實驗數據獲得的球床內的極緩慢稠密顆粒流動信息，采用自相關時間τc和間歇性指數C兩種特征參數在不同均值流線上的分布，可研究不同流線上的球流速度相關性、間歇性，進而對球床分區進行分析。如球床上部不同流線間自相關性較高，說明流動的協調性較高，體現整體流特征，顆粒隨機運動及橫向擴散小。隨著高度的下降，球流隨機性增強，橫向擴散增強，τc也逐漸減小，兩者關系呈現冪律分布τc～Hα(H為高度，α為整體流水平指數)，該規律尚需進一步研究。在徑向方向上，越靠近側壁，球流受壁面影響協調性增強，自相關性增強[25]。引入流動均勻性指數(σ)和α，通過單個顆粒球速度矢量與球床流場合速度矢量方向夾角的統計分布衡量流場運動一致性水平。

同時，基于PTV數據對球床內的搭橋現象進行了分析，發現搭橋概率分布函數在半對數坐標系下，呈現二階多項式分布特性，在近壁區，搭橋顆粒呈現較為規則的堆積結構，而在中心區，搭橋顆粒排列取向的交錯性較強。搭橋顆粒球的連續相鄰顆粒球的夾角的概率分布函數表征了搭橋的穩定性。通過互相關函數方法，證明了搭橋崩塌和顆粒空間拓撲結構變化存在一定相關性。通過計算球床內典型窗口區域內的搭橋顆粒平均速度與窗口平均速度的互相關函數，發現搭橋顆粒速度和窗口平均速度呈現較高的相關性。搭橋持續特征時間和速度自相關時間呈現正相關特征。搭橋動力學特征應是球床內速度脈動的主要原因[28]。

1.2 球床等效導熱系數實驗

清華大學HTR-PM三維堆芯球床等效導熱系數測量實驗(TS-HTTF)分兩個階段實施。第1階段是建立一個比例縮小的前期結構驗證和材料篩選試驗裝置，構建一個模擬高溫氣冷堆溫度的環境，測試溫度要求達到1 600 ℃，具有與高溫氣冷堆相同性質的氣體氛圍，進行加熱材料、測溫元件材料和保溫材料的篩選，以及加熱元件結構、爐膛密封結構、水冷結構和絕緣結構等關鍵結構的驗證[29]。第2階段是建立HTR-PM三維堆芯全尺寸球床等效導熱系數實驗系統[28]。該裝置直徑5 m，床高7 m，內部含70 000個隨機堆積的石墨球，分布于內徑1.2 m、外徑4.2 m、高度1 m的環形區域中，并保持20 Pa壓力。布置90個熱電偶，其中周向5個(C1～C5)，高度方向3個，徑向6個(T1～T6)。目前完成了真空條件下靜態填充球床內最高達1 200 ℃的瞬態法等效導熱系數測試，并與南非的HTTU數據進行了比較[29-30]。圖4示出了球床等效導熱系數實驗臺測點布置、瞬態床溫變化及等效導熱系數對比。通過真空工況下純堆積造成的球床等效導熱系數與有冷卻劑存在的情況下球床等效導熱系數的比較，可深入分析球床等效導熱系數的關鍵影響因素。按照實驗計劃，還將進一步開展氦氣和氮氣氛圍下1 200 ℃和1 600 ℃的球床等效導熱系數實驗測試。

2 理論與數值研究

2.1 球流模型及模擬

1) 球流紡及橫向擴散

借助于球流實驗及離散元數值模擬，研究了球床內球流紡結構及球流紡內的橫向擴散分布函數[31]。實驗和計算均發現球流紡分布區域較窄，球床本體部分幾乎與球流紡相互平行，在基底收縮向泄流口(圖5a)。對球流紡的統計分析發現內部球的橫向分布呈現近高斯分布。由于球床內為極緩慢的球流，采用準靜態分析發現球流的概率分布結構為楊輝三角結構(圖5b)，沿高度方向向下，其高概率區間(95%)分布于非常窄的帶形區域內，即球流紡。由于楊輝三角橫向分布為二項分布，其漸進分布為高斯分布，因而解釋了球流紡內的徑向擴散為近高斯形式。

圖4 球床等效導熱系數實驗臺測點布置(a)、瞬態床溫變化(b)及等效導熱系數對比(c)Fig.4 Schematic of vertical cut through test facility and thermocouples distribution among whole pebble bed (a),average temperature of five sets in middle level in whole test period (b), and comparison between three effective thermal conductivities derived by inverse method and steady-state method (c)

圖5 實驗和計算中的球流紡(a)、球流紡內的徑向概率分布 及楊輝三角分布(b)以及高概率區間(95%)的正視和側視圖(c)Fig.5 Experimental snapshots of particles in granular bed and particle trajectories originated from nine fixed points at bed top (a), comparison of normalized PDFs between numerical and experimental results and probability distribution of diffusion process on all levels and positions (b), and ideal diffusion process on the first 20 levels (colored by probability on each level) and the probability view of ideal diffusion process (c)

2) 球流均勻性、滯留區(率)及床幾何構型優化

為減小滯留區并分析球床基底結構對球流均勻度和滯留區的影響，分析了給定基底高度和出口直徑的情況下，過渡曲線為圓心在斜面中垂線上的{Rn，n=1,2,3}圓弧曲線基底的球床及與直線母線圓錐{Rn，n=∞}和最速降線形狀的比較(圖6a)[32]。發現在給定循環球數下，最速降線內的球流分布最均勻，每個球的循環概率最集中，滯留率最小(圖6b)。而從流型上看，球流幾乎不存在滯留區(圖6c)。僅次于最速降線的是R2構型。目前采用的直線母線的圓錐構型較R2構型更差一些，說明在加工條件和成本允許的情況下應盡量采用最速降線構型或R2構型。

進一步分析可發現，當采用更大的底面傾角時，如采用60°傾角替代目前的30°傾角，即使采用直線母線圓錐構型，整個滯留區亦將消失[33]。但60°傾角帶來的問題是球床底部泄流口的受力要重新校核。

此外，對球床內施加內錐結構進行導流也可對床內流型進行調制，內錐應分布于基底內不超過基底高度且接近泄流口的位置[34]。但內錐的存在需要支撐，而支撐可能會對流型施加不利影響，因此還需進一步研究。

球床優化目前已采用的手段包括用近壁防結晶凹槽防止出現球流結晶。對三角、半圓及梯形結構凹槽的分析和比較可發現，有凹槽結構相對于光滑無凹槽結構能避免近壁球堆積結晶的形成，而三角凹槽在1d(d為球徑)深度和1.67d間距的布置下相對于其他形狀的凹槽的防結晶效果最好[35]。

3) 物性參數對球流影響

數值模擬可分析球-球摩擦與球-壁摩擦(圖7)，以及不同部位的球-壁摩擦系數對球流流型和流動性的影響。

研究中采用不同的流場評價指標來定量地描述球流流場，包括整體流指數、滯留率、垂直標準差和平均次序值[36]，分別反映了球流流場的滯留現象、流動的徑向均勻性和流動次序性。利用定量評價指標，可評判不同摩擦系數組合下球流流場各方面的表現。例如，整體流中所有顆粒物料處于一致的運動中，并以均勻的速度向下流動，顆粒流動服從“先進先出”序列，上下層顆粒具備較好的流動次序，各顆粒在流場中經歷的時間相當，這點在球床堆中尤其重要，可保證燃料球燃耗均勻，與反應堆安全性、經濟性直接相關。相比之下，漏斗流的流場徑向差異大，中心區流動明顯快于近壁區，近壁區顆粒流速慢，易滯留。整體流指數(Mass Flow Index, MFI)定義為流場中同高度的貼壁速度與中心速度的比值，當該值大于0.3時，流型被視為整體流，反之，該值小于0.3時，定性為漏斗流。

圖7 球-壁摩擦(a)及球-球摩擦(b)對球流型的影響Fig.7 Cross-sectional distribution of vertical velocity with different particle-wall (a) and particle-particle (b) friction coefficients

比較發現(圖8)，減小球-壁摩擦和球-球摩擦均可改善球流流場的均勻性，而球-壁摩擦是決定流型的主要因素。摩擦效應變化與流型改變之間的關系是非線性的。隨著壁面摩擦逐漸減小，流型會在一段摩擦系數(μ)的區間內集中完成由非均勻漏斗流向均勻整體流的轉變。所以，通過確定該流型轉變區的區間，并將壁面摩擦系數控制在轉變區下限以下，即可維持球床內為均勻流場。通過比較球床不同部位的壁面摩擦對整體流場的影響貢獻，證實基底部分的摩擦影響明顯大于本體部分，建議工程設計上應優先考慮降低球床基底部分的壁面摩擦，而在工程施工中應盡量避免基底斜坡的凹凸，或提高基底斜面的光滑性。

圖8 球-壁摩擦(a)及球-球摩擦(b)對球流速度的影響Fig.8 Radial distribution of vertical velocity with different particle-wall (a) and particle-particle (b) friction coefficients

4) 極緩慢流態

顆粒流理論方面的研究比較滯后，顆粒流流態的定義和劃分仍是模糊和粗糙的。比較受認可的理論是Campbell 20世紀90年代提出的依據顆粒間是否明顯存在接觸力鏈的分類方法[37]：能形成穩定接觸力鏈的是彈性流，無明顯力鏈的為慣性流。彈性流一般是稠密流，又被細分為彈性-準靜態流(接觸應力與剪切速率無關)和彈性-慣性流(接觸應力與剪切速率呈正比)。慣性流較稀疏，顆粒間為瞬間的碰撞行為，難以形成穩定力鏈。根據碰撞是否為嚴格的兩體碰撞，慣性流又被細分為慣性-碰撞流(嚴格兩體碰撞)和慣性-非碰撞流(非嚴格兩體碰撞)。

極緩慢流體的流態極為特殊，可從不同角度對球流流態進行分類。鑒于球流內存在間歇性瞬時速度的局部整體滑移特性(圖9)，可從如下幾方面考慮。

(1) 從能量角度考慮。提出了重力驅動稠密顆粒流的系統動能級數符合的單峰形概率密度分布，并依據統計分布規律定義了能級標準差和能級跨度(圖10)。采用能級跨度和能級標準差將重力驅動稠密顆粒流細分為3個子流態(表1)，分別為：極緩慢流(能級跨度SOE>1，能級標準差σ<0.1)、連續流(SOE<1，σ<0.1)和間歇流(σ>0.1)[38]。

a——發生整體滑移的流場演化過程； b——對應的速度、加速度、法向接觸力的變化圖9 球床內顆粒速度場的整體滑移特性Fig.9 Whole slip characteristic of particle velocity field in pebble bed

圖10 能級跨度和能級標準差Fig.10 Span of energy magnitude and standard deviation

表1 流態劃分采集時間獨立性檢驗Table 1 Independence test on simulation time

圖11 基于參數分析的稠密顆粒流細分流型圖Fig.11 Diagram of parameter analysis for categorizing dense pebble flow

(2) 從流動間歇性特征角度考慮。首先采用無量綱的脈動動能指數(RPKE)描述顆粒流的脈動性。并基于多重分形分析提出了間歇性指數C，表征了相鄰脈動運動的平均時間間隔，在時間尺度上描述連續均勻運動的相對時間。基于C和RPKE，將其分為3種典型流型(圖11)：顯著間歇流(C≥0.5)；連續流(C<0.5，RPKE≤0.1)；過渡脈動流(C<0.5，RPKE>0.1)。顯著間歇流，球床內顆粒球在JP3絕大多數時間內保持接近于靜止的狀態。此時，脈動運動以較大的速度為特征，相鄰脈動運動的平均時間間隔較大。連續流，球床流場呈現較快的流動速度，其對應的顆粒流動呈現較小的脈動特性，流動較為平穩。過渡脈動流，涵蓋了明顯間歇流到連續流的過渡區域[39]。

(3) 從主導流態的機理因素考慮。將動力學量主導的流態稱之為動力學流態，典型的表征量是顆粒受力(或加速度、位矢的二階導數)；對應地將運動學量主導的流態稱為運動學流態，典型的表征量是速度(位矢的一階導數)。在數學上相關系數為正則化的協方差，表示兩個獨立變量之間的關聯性，其值為1或-1表示正相關或負相關，0表示不相關。對球流的平均速度和平均受力進行相關性分析，可發現在高速球流中力的相關系數ρFF無明顯周期性，而低速球流中ρFF有明顯周期性和主頻(圖12a)；反之，在高速球流中速度的相關系數ρVV有明顯周期性和主頻，而低速球流中無明顯周期性和主頻(圖12b)[40]。

圖12 球床受力和速度的相關系數及其頻譜Fig.12 Autocorrelation functions and power spectra of mean force and mean velocity

圖13 短程、長程和微觀表面輻射模型對球床等效導熱系數估計的對比(a) 以及短程輻射模型與南非HTTU實驗結果的比較(b)Fig.13 Effective thermal conductivities predicted by short-range, long-range and microscopic radiation models (a) and comparison of simulation results obtained by short-range radiation model with HTTU experiment result (b)

2.2 球床高溫熱輻射模型

為采用數值方法研究球床高溫輻射，在不同空間尺度下建立了長程輻射模型(LRM)與短程輻射模型(SRM)、微觀表面輻射模型(MSM)，以及亞顆粒輻射模型(SCM，又稱子單元輻射模型)[41-45]。長程輻射模型中考慮了所有存在輻射換熱關系的顆粒對，而短程輻射模型僅考慮距離較近的單層Voronoi單元相鄰顆粒對。因模型中認為材料內導熱系數為無限大，結果表明當顆粒材料導熱系數遠大于輻射有效導熱系數時，長程輻射模型適用于高溫球床輻射計算，但當兩者在同一數量級時，長程模型高估了材料內導熱系數的影響，而短程模型因截斷了長程輻射熱流，從而抵消了對材料導熱系數的高估，顯示出良好的計算精度[41]。微觀表面輻射模型通過對顆粒表面的詳細剖分并詳細計算角系數，呈現較高的計算精度，但其計算量過于龐大，無法滿足存在大量顆粒的球床工程計算的需求。圖13示出短程、長程和微觀表面輻射模型對球床等效導熱系數ke估計的對比以及短程輻射模型與南非HTTU實驗結果的比較。圖中，ZBS為Zehner-Bauer-Schlünder。圖14示出亞顆粒模型預測的等效輻射導熱系數kr與短程輻射模型、經驗模型的比較。

圖14 亞顆粒模型預測的等效輻射導熱系數 與短程輻射模型、經驗模型的比較Fig.14 Effective thermal conductivity for particle radiation predicted by SCM and comparison with SRM and empirical correlation results

進一步，推導出了短程輻射模型和亞顆粒輻射模型下材料有效導熱系數的解析表達式。計算結果表明，短程輻射模型在溫度超過1 200 ℃時略有偏高，而亞顆粒輻射模型是通用的理論性方法，在超高溫度下(大于1 500 ℃)仍適用[44-45]。

3 總結

通過對清華大學核能與新能源技術研究院近幾年對球流和球床等效導熱系數的主要研究成果的回顧和總結，實驗和模擬結果極大地加深了對球床反應堆的內在規律和機理的認識，給出了球床內關鍵工程參數的準確估計。因此，可認為球流實驗和球床等效導熱系數實驗有效揭示了球流的基本規律，而球流模擬及基于元胞建模的顆粒輻射模型能為球流研究、堆芯優化、重要熱工參數的估計和比較提供有效的幫助。

但球流仍存在較多需解決的理論、計算、實驗和工程問題，僅列如下幾方面。

理論方面：1) 如何給出能描述球流在擬連續介質力學框架中的力學本構方程，以能跨越離散介質和連續介質描述的鴻溝，實現理論的統一性；2) 如何給出極緩慢球流分區、緩慢與快速流動的物理本征特征及數學表達。

計算方面：1) 在離散介質力學模型方面如何開發更高效的算法，提高計算效率并保證計算精度；2) 如何開發球流計算與球床輻射導熱計算、球床冷卻劑流動傳熱計算的多重耦合模型及高精度算法。

工程方面：1) 如何給出球流的間歇頻率、滯留特性與球床、球等幾何、物性和運行參數之間的關系，以更好地指導球床堆在不同運行環境及應用中的反應堆設計；2) 球床球流設計與反應堆物理設計的接口及迭代問題。

上述內容是將來還可進一步深入探討的問題，因此下一步將在球流全3D測量、球流及球床等效導熱系數快速計算、球床式反應堆內的氦氣和球流兩相耦合及換熱的測量和建模等方面進行深入研究。