張寶銳，夏兆陽，周志偉

(清華大學 核能與新能源技術研究院，北京 100084)

聚變堆包層作為聚變堆的關鍵部件，承擔著氚增殖、能量轉換及輻射屏蔽等功能。包層氚增殖劑按其形態(tài)可分為固態(tài)包層和液態(tài)包層，其中固態(tài)包層采用鋰陶瓷顆粒作為氚增殖材料，球床采用低壓氦氣作為吹掃氣體將氚載出，并摻雜一定量的氫氣以加快氚的載帶[1-2]。球床的內部結構會對吹掃氣體的流動特性及球床內的傳熱傳質行為產(chǎn)生影響，因此獲取球床氣體的詳細流動特性及傳熱傳質特性對包層氚輸運分析及氚提取系統(tǒng)設計具有重要意義。

在球床堆積方面，Gong等[3]采用不同尺寸的鋁合金顆粒進行了一元及二元尺寸球床的堆積實驗，填充率均可達到設計目標。他們同時采用離散單元法(DEM)計算程序進行了數(shù)值模擬，填充率和孔隙率分布與實驗結果符合得較好。在球床流動阻力特性研究方面，Abou-Sena等[4]基于其設計的壓降測量裝置，測量了不同填充因子、流速及粒徑下的氦氣流動壓降，驗證了Ergun方程的適用性；伍振興等[5]發(fā)現(xiàn)在顆粒粒徑較小(低于1 mm)時，Ergun方程預測值低于實驗值，通過理論推導得出了考慮氦氣可壓縮性影響的修正Ergun方程，與實驗值符合良好；汪衛(wèi)華等[6]首先采用計算流體力學(CFD)方法對規(guī)則堆積的氚增殖球床進行了數(shù)值模擬，獲得了球床內的流速及壓力分布；Zhang等[7]、陳有華等[8]分別對一元及二元球床進行了DEM-CFD耦合計算，獲得了隨機堆積球床內的吹掃氣體速度及壓力分布，發(fā)現(xiàn)在壁面處出現(xiàn)流體加速，而在堆積較密實的區(qū)域存在流動遲滯區(qū)。在球床等效導熱系數(shù)研究方面，Donne等[9]、Liu等[10]采用穩(wěn)態(tài)熱流法分別對不同堆積率下的硅酸鋰球床進行了等效導熱系數(shù)的測量實驗，獲得了等效導熱系數(shù)與平均溫度的線性擬合關系；駱貝貝等[11]在中國先進研究堆(CARR)中開展了輻照環(huán)境下硅酸鋰球床等效導熱系數(shù)的測量，并與理論值和堆外實驗擬合值進行了比較。

球床內部結構復雜，為減小計算量，目前針對固態(tài)包層模塊的熱工水力及氚輸運分析將氚增殖球床等效為固體域或多孔介質域，因此需獲得球床結構下流體流動特性及球床等效導熱系數(shù)，并對氚在吹掃氣體內的擴散系數(shù)進行修正。本文以核工業(yè)西南物理研究院提出的一種氦冷固態(tài)包層概念為研究對象[12]，采用離散單元法生成滿足填充率要求的球床隨機堆積結構，并對該精細模型開展吹掃氣體流動特性研究與球床導熱模擬，為后續(xù)包層的傳熱傳質多物理場耦合模擬提供宏觀輸入?yún)?shù)。

1 隨機堆積模型與網(wǎng)格劃分

1.1 球床隨機堆積模型的DEM計算

氦冷固態(tài)包層采用粒徑為1 mm的硅酸鋰顆粒作為氚增殖劑，為模擬真實的球床堆積結構，本文采用DEM對硅酸鋰球床的堆積結構進行數(shù)值模擬。DEM的基本思想是將每個顆?？醋骶哂匈|量、形狀、速度等參數(shù)的獨立單元，基于不同的接觸模型，在每個時間步內確定顆粒的接觸力，并對顆粒的加速度、速度、位移等參數(shù)進行更新，直至球床達到應力平衡狀態(tài)。計算流程如圖1所示。

圖1 DEM計算流程Fig.1 Simulation procedure of DEM

DEM計算程序采用PFC3D，所需的硅酸鋰顆粒物性參數(shù)列于表1[3]。圓柱形球床工業(yè)應用廣泛，其填充結構也有較成熟的研究，在球床直徑與顆粒直徑比小于10時，壁面效應會對球床填充率有較大影響[13-15]，因此本研究中設定計算模型為半徑6.5 mm、高度9.0 mm的圓柱體，初始時刻在圓柱體頂部生成位置隨機的硅酸鋰顆粒，隨后硅酸鋰顆粒在重力作用下落入圓柱體內部，顆粒-顆粒與顆粒-壁面間接觸力的確定基于Hertz接觸模型[16]。

表1 用于球床堆積模擬的硅酸鋰物性參數(shù)Table 1 Main parameter of Li4SiO4for DEM simulation

經(jīng)多次迭代后球床達到應力平衡狀態(tài)，得到氚增殖球床隨機堆積結構，如圖2所示。

圖2 硅酸鋰球床離散堆積模型Fig.2 Structure model of Li4SiO4 pebble bed

球床模型共模擬了1 348個硅酸鋰顆粒，總體堆積率為59%，球床內徑向孔隙率驗證多采用Klerk提出的徑向孔隙率分布經(jīng)驗關系式[17]：

(1)

其中：εr為徑向孔隙率；R=(r-rc)/d為無量綱距離，r為測點與圓柱體中軸距離(0≤r≤rc)，rc為圓柱體半徑，d為顆粒粒徑；εb為球床中心平均孔隙率。

模型徑向孔隙率分布與Klerk關系式的對比如圖3所示。由圖3可知，模型徑向孔隙率分布與經(jīng)驗值符合良好，可認為DEM計算得到的球床堆積結構真實合理。

圖3 徑向孔隙率分布Fig.3 Radial porosity distribution

1.2 CFD模型重建與網(wǎng)格劃分

1) CFD模型重建

采用有限元計算軟件COMSOL Multiphysics進行氚在球床結構內的等效擴散系數(shù)及吹掃氣體流動的數(shù)值模擬。COMSOL除了可視化建模，其提供的方法功能也可運行COMSOL自己開發(fā)的類JAVA語法命令流實現(xiàn)參數(shù)化建模，同時COMSOL提供了與MATLAB的LiveLink接口[18]，可實現(xiàn)數(shù)據(jù)及命令流的實時交互。因此本文基于MATLAB，讀取DEM計算得到的穩(wěn)定堆積狀態(tài)下硅酸鋰顆粒的位置坐標，并通過LiveLink接口實現(xiàn)在COMSOL內顆粒幾何結構的重構。模型重構流程如圖4所示。

圖4 COMSOL幾何重建流程Fig.4 Pebble bed structure rebuilding procedure in COMSOL

重建的CFD計算模型如圖5a所示，球床在軸向設置了相應的進出口段以減小氣體流動的出入口效應。通過布爾運算抽取球床內的復雜流道進行網(wǎng)格劃分，在網(wǎng)格劃分時顆粒接觸位置網(wǎng)格變形較大，因此需對模型內顆粒間接觸進行處理。常用于球床顆粒間接觸的處理方法有縮徑法[19]、擴徑法[20]、剖切法[21]、搭橋法[22]等。縮徑法的思路是通過減小顆粒半徑以避免顆粒間的點接觸，擴徑法、搭橋法的思路是將點接觸擴大為面接觸以降低網(wǎng)格劃分難度。其中搭橋法最大程度地保留了球床顆粒間的接觸信息，且用于搭橋的圓柱體半徑一般小于顆粒半徑的1/10，對球床的孔隙率分布影響極小，因此本文采用搭橋法處理顆粒間接觸，處理后的顆粒間接觸如圖5b所示。

a——重建的COMSOL模型；b——搭橋法處理的顆粒接觸圖5 重建的COMSOL模型與處理后的顆粒接觸Fig.5 Rebuild model in COMSOL and particle contact after treatment

2) 網(wǎng)格無關性驗證

氦冷固體包層采用氦氣作吹掃氣體，運行壓力為101 kPa。以流動場計算為例，計算模型采用速度入口、壓力出口，內部顆粒壁面及圓柱容器壁面為無滑移邊界，在吹掃氣體流速范圍(0.1～2.0 m/s)內Re均處于層流模型適用范圍，因此計算模型采用層流穩(wěn)態(tài)計算。模型采用四面體網(wǎng)格，分別劃分了224萬、292萬、326萬與401萬4套網(wǎng)格進行網(wǎng)格無關性驗證，入口流速0.1 m/s下的球床壓降計算結果列于表2。由表2可知，當網(wǎng)格量大于326萬后壓降計算結果差別很小(相對誤差小于1%)，因此本文采用網(wǎng)格量為326萬的計算結果作為分析對象，其網(wǎng)格劃分示意圖如圖6所示。

表2 網(wǎng)格無關性驗證Table 2 Mesh independence verification

2 計算結果分析與討論

2.1 球床結構下氚在吹掃氣體內的等效擴散系數(shù)

氚在氦氣中的自由擴散系數(shù)由Chapman-Enskog公式計算：

(2)

其中：k為玻爾茲曼常數(shù)；Tg為氣體溫度；M為等效摩爾質量；ni、nj分別為組分i、j的原子數(shù)密度；ε0為修正因子；Ωi,j為擴散碰撞積分因子，詳細推導過程參見文獻[23]。593 K下，氚在氦氣中的擴散系數(shù)為1.78×10-4m2/s。

在球床結構下，復雜曲折的內部區(qū)域使氣體分子平均輸運路徑增加，因此若使用多孔介質模型對氚增殖球床進行氚的輸運行為模擬，需對氚在氦氣中的擴散系數(shù)進行修正。具體計算流程如下：1) 建立帶有入口段和出口段的球床隨機堆積模型；2) 物理場采用稀物質傳遞場，在入口段和出口段設定濃度邊界；3) 對稀物質傳遞場進行穩(wěn)態(tài)計算；4) 由通量連續(xù)條件，得到球床段等效擴散系數(shù)。

穩(wěn)態(tài)下氚濃度分布如圖7所示，計算模型入口段、球床段、出口段通量連續(xù)，即：

圖7 穩(wěn)態(tài)下氚濃度分布Fig.7 Concentration of tritium in pebble bed under steady state

J入口段=J球床段=J出口段

(3)

由Fick定律，可得球床段氚的等效擴散系數(shù)：

(4)

其中：D為氚在氦氣中的擴散系數(shù)；Deff為球床結構下氚的等效擴散系數(shù)；ΔcT為各段氚濃度差值；Δx為各段長度。因此，若采用多孔介質模型對氚增殖球床進行氚輸運模擬，球床結構影響下氚在吹掃氣中的等效擴散系數(shù)Deff=0.4D。

2.2 吹掃氣體流動壓降分析

球床通道內流體的流動壓降用Ergun方程描述：

(5)

其中：ε為球床孔隙率；μ為氦氣動力黏度；dc為球床直徑；ρ為氦氣密度；C1與C2分別為黏性阻力與Forchheimer阻力系數(shù)。層流下，式(5)可簡化為Blake-Kozeny壓降方程：

(6)

不同入口流速下球床通道內流體的流動壓降示于圖8，擬合所得Blake-Kozeny方程系數(shù)C1=87。

圖8 不同入口流速下的流動壓降Fig.8 Pressure drop under different inlet velocities

2.3 吹掃氣體流動特性分析

吹掃氣體軸向速度分布如圖9所示?？煽吹剑虼矃^(qū)域氣體流速變化劇烈，吹掃氣體不斷經(jīng)歷加速、減速過程，在孔隙率較大的位置由于流動阻力較低，出現(xiàn)氣體流動加速現(xiàn)象。在球床內部填充較緊實的區(qū)域顆粒背部氣體流速緩慢，出現(xiàn)較多流動遲滯區(qū)(流速小于10-2m/s)，在近壁面區(qū)域由于球床壁面效應，孔隙率較大，流動阻力降低，氣體流速較大。

圖9 y=0截面速度云圖Fig.9 Velocity distribution of axial central cross section

2.4 吹掃氣體流動特性多孔介質模型驗證

為對球床精細模型得到的流動阻力特性進行驗證，采用自定義徑向孔隙率的二維軸對稱模型進行數(shù)值模擬，如圖10a所示，模型尺寸與球床精細模型保持一致，采用結構化網(wǎng)格進行網(wǎng)格剖分，如圖10b所示。流動壓降與多孔介質模型計算結果示于圖11，可見二者符合較好。徑向流速分布如圖12所示，呈現(xiàn)與孔隙率分布(圖3)相反的波動現(xiàn)象，即在孔隙率較低處流動阻力減小，流動加速。

圖10 二維軸對稱多孔介質模型Fig.10 Two-dimensional axisymmetric porous media model

圖11 壓降計算結果對比Fig.11 Comparison of pressure drop calculation results

圖12 徑向流速分布Fig.12 Radial flow velocity distribution

2.5 球床等效導熱系數(shù)模擬

駱貝貝等[11]針對硅酸鋰球床開展了堆內輻照環(huán)境下的等效導熱系數(shù)的測量實驗，并與理論值和堆外實驗結果進行了對比，結果表明輻照實驗下硅酸鋰球床的等效導熱系數(shù)略低于理論值和堆外實驗值。由于ZBS模型及堆外實驗均基于外加熱流的邊界條件，因此本研究進行了外加熱流與內熱源兩種工況下的模擬分析。由于提氚氣體(氦氣)流速較低，從保守角度假設氦氣為靜止狀態(tài)，為減小模型計算量及降低壁面效應的影響，提取球床內部徑向孔隙率變化較低的部分進行模擬計算。模型如圖13所示，藍色區(qū)域為硅酸鋰顆粒，灰色區(qū)域為氦氣，采用四面體網(wǎng)格進行網(wǎng)格劃分，在顆粒-氦氣及顆粒-顆粒接觸處進行局部加密，總網(wǎng)格量達1 122萬，網(wǎng)格劃分如圖14所示。

圖13 球床導熱模型Fig.13 Thermal conductivity model of pebble bed

圖14 球床網(wǎng)格劃分示意Fig.14 Schematic view of mesh

1) 外加熱流工況

該工況下模型頂部為加熱面，熱流密度為1 000 W/m2，模型底部為定溫壁面，其余邊界為絕熱邊界。改變模型底部溫度，分別計算底部溫度為473、573、673、773、873 K時球床平均溫度對應的等效導熱系數(shù)。由于熱流連續(xù)，等效導熱系數(shù)可直接由傅里葉定律計算：

keff=q″/(?T/?z)

(7)

其中：keff為球床等效導熱系數(shù)，W/(m·K)；q″為模型z方向熱流密度，W/m2；?T/?z為z方向溫度梯度，K/m。對5種平均溫度對應的等效導熱系數(shù)進行線性擬合，得到:

keff=0.951+2.778×10-4T

(8)

其中，T為球床平均溫度，K。

2) 內熱源工況

設定該工況下硅酸鋰顆粒為體熱源，功率密度為1 MW/m3，模型圓周面為定溫邊界，上下表面為絕熱邊界。同樣進行了定溫邊界溫度為473、573、673、773、873 K時球床平均溫度對應的等效導熱系數(shù)計算。由于有內熱源工況下球床內傳熱路徑復雜，不滿足熱流連續(xù)條件，因此從導熱控制方程出發(fā)進行等效導熱系數(shù)推導。

對于硅酸鋰顆粒、氦氣，有：

(9)

其中：ks、kf分別為硅酸鋰顆粒與氦氣的導熱系數(shù)，W/(m·K)；qv為硅酸鋰顆粒的功率密度，W/m3。

將球床等效為均勻多孔介質，則球床整體滿足有內熱源下的圓柱導熱方程：

(10)

(11)

其中，tw為定溫邊界溫度。對模型不同半徑處溫度取平均值，可得到平均溫度下的等效導熱系數(shù)，經(jīng)線性擬合，有：

keff=0.834+2.818×10-4T

(12)

3) 球床等效導熱系數(shù)分析

Liu等[10]對采用穩(wěn)態(tài)熱流法對堆積率為59.7%的硅酸鋰球床的等效導熱系數(shù)進行了實驗測量，線性擬合曲線為：

keff=0.753+5.17×10-4T

(13)

駱貝貝等[11]基于CARR堆內輻照實驗得到的球床導熱系數(shù)與溫度的關系為：

keff=0.647+4.074×10-4T

(14)

外加熱流工況、內熱源工況、CARR堆內輻照實驗、Liu等[10]的穩(wěn)態(tài)熱流實驗得到的球床等效導熱系數(shù)對比示于圖15。由圖15可看出，內熱源工況下的球床等效導熱系數(shù)均略低于外加熱流工況下的等效導熱系數(shù)，CARR堆內輻照實驗的球床等效導熱系數(shù)最低，可能是由于CARR實驗中硅酸鋰球床堆積率較低(56%)，導致球床接觸熱阻相對較大；而本研究中兩種工況采用相同的球床結構，有內熱源工況下的球床等效導熱系數(shù)仍低于外加熱流工況。劉子平等[24]對高溫氣冷堆采用的包覆型顆粒燃料的等效導熱系數(shù)進行了理論推導，同樣得到在有內熱源情況下燃料球等效導熱系數(shù)低于無內熱源情況下等效導熱系數(shù)。相較于外加熱流工況，在有內熱源工況下，球床內傳熱工況更加復雜，外加熱流下的熱阻簡單并聯(lián)模型并不適用，內熱源的存在降低了熱流方向上的溫度差，即增加了熱流方向上的熱阻，導致球床整體等效導熱系數(shù)降低。有內熱源工況更符合硅酸鋰球床實際運行狀態(tài)，從保守分析角度，建議球床導熱系數(shù)采用有內熱源工況下的球床等效導熱系數(shù)。

圖15 硅酸鋰球床等效導熱系數(shù)模擬值與實驗值對比Fig.15 Comparison of simulated and experimental results of Li4SiO4 pebble bed

3 結論

通過離散單元法生成了硅酸鋰顆粒隨機堆積模型，并對球床結構下氚在吹掃氣體內的等效擴散系數(shù)及吹掃氣體流動及球床導熱進行了CFD計算。計算結果表明：球床復雜曲折的內部結構使氣體分子平均輸運路徑增加，氚在吹掃氣體內的等效擴散系數(shù)Deff=0.4D；吹掃氣體在球床內的流動壓降可由Blake-Kozeny方程預測，經(jīng)擬合得到Blake-Kozeny方程系數(shù)C1=87；內熱源的存在降低了熱流方向上的溫度差，導致球床整體等效導熱系數(shù)低于外加熱流工況下的球床等效導熱系數(shù)。

本文通過對硅酸鋰球床開展DEM-CFD耦合計算，獲得了氚在吹掃氣體內的等效擴散系數(shù)、吹掃氣體流動壓降及球床等效導熱系數(shù)等宏觀參數(shù)，可在此基礎上基于多孔介質模型對固態(tài)包層進一步開展球床傳熱及氚輸運行為模擬。此外，本文結果對包層氚提取系統(tǒng)設計也具有一定的參考價值。