球床反應(yīng)堆體心棱柱單元局部高溫?zé)狳c(diǎn)與渦特性分析

周銳韜，郭雪巖

（上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093）

高溫氣冷堆(HTR)是第四代核反應(yīng)堆，以石墨為慢化劑包覆天然鈾組成球形燃料元件，氦氣流通球床進(jìn)行冷卻，從而實(shí)現(xiàn)反應(yīng)堆的固有安全性。固有安全性可以保證燃料不熔化，但局部高溫?zé)狳c(diǎn)依然可能破壞燃料球。因而，為了保證燃料球的完整性，尋找球床反應(yīng)堆內(nèi)部高溫區(qū)域是十分必要的。傳統(tǒng)實(shí)驗(yàn)方法在尋找堆積球內(nèi)部高溫點(diǎn)時(shí)比較困難，所以越來越多研究者開始采用數(shù)值分析方法探究球床內(nèi)部溫度問題。

實(shí)際球床反應(yīng)堆中燃料球隨機(jī)堆積，整體結(jié)構(gòu)非常復(fù)雜，因而很多研究都對(duì)球床問題做了簡化，并給出了四種規(guī)則排列方式：簡單立方(SC)、體心立方(BCC)、面心立方(FCC)和體心棱柱(BCP)。Hassan于2007年用大渦模擬(LES)方法對(duì)BCC單層結(jié)構(gòu)建立了非結(jié)構(gòu)網(wǎng)格并進(jìn)行了分析。Ferng等于2013年采用雷諾應(yīng)力模型(RSM)對(duì)BCC與FCC單層結(jié)構(gòu)進(jìn)行了研究。Shams等于2013年采用近似直接數(shù)值求解(q-DNS)方法對(duì)FCC單層結(jié)構(gòu)的流場溫度場進(jìn)行求解。蔣旭等于2017年在研究BCC結(jié)構(gòu)時(shí)比較了-模型與LES方法的差異。姚強(qiáng)等于2019年對(duì)8層排列FCC球陣建立結(jié)構(gòu)化網(wǎng)格，并采用LES方法進(jìn)行數(shù)值求解。吳浩等于2019年采用DEM模型對(duì)HTR-10球床氣固傳熱過程進(jìn)行了模擬計(jì)算，其結(jié)果基本與經(jīng)驗(yàn)?zāi)Ｐ鸵恢隆堧p寶等于2020年建立耦合模型模擬穩(wěn)態(tài)流動(dòng)換熱過程，并將CFD結(jié)果與Thermix計(jì)算結(jié)果作比較，事實(shí)證明CFD模型精準(zhǔn)可靠。孫世妍等于2020年研究了不同工況下HTR-10堆芯溫度場的分布，以確定最高溫度變化與冷卻劑的關(guān)系并尋找高溫極限。結(jié)合眾多研究成果可以發(fā)現(xiàn)，規(guī)則排列球陣中，有關(guān)BCP結(jié)構(gòu)的研究較少；RANS湍流模型由于捕捉瞬態(tài)流場的不足和對(duì)流動(dòng)分離的低敏感性，在研究多曲面流道的球床問題上并不適用；使用結(jié)構(gòu)網(wǎng)格在流道狹窄處可以得到更精細(xì)的結(jié)果，適合復(fù)雜的球床問題；單層結(jié)構(gòu)由于整體尺度小，受出入口效應(yīng)影響較大，可能最終會(huì)影響溫度分布，使用多層結(jié)構(gòu)并選取中間層進(jìn)行研究可以較好地避免這一問題。

綜上所述，本文選擇四種規(guī)則排列方式中空隙率最小、排列最緊密的體心棱柱(BCP)結(jié)構(gòu)，燃料球間接觸位置以面接觸方式處理，湍流模型選擇LES方法，對(duì)8層燃料球陣建立結(jié)構(gòu)網(wǎng)格求解溫度場與速度場，尋找流場內(nèi)高溫區(qū)域，分析高溫區(qū)域產(chǎn)生的原因。同時(shí)，充分發(fā)揮LES方法的優(yōu)越性，尋找并分析球床內(nèi)部不同位置渦對(duì)流動(dòng)換熱的影響。

1 計(jì)算模型及方法

1.1 計(jì)算模型

在規(guī)則球陣數(shù)值計(jì)算研究中，出入口效應(yīng)的影響不可忽略，單層結(jié)構(gòu)若不經(jīng)過特殊設(shè)置其結(jié)果必然存在較大偏差。建立多層結(jié)構(gòu)，然后選取其中間層的燃料球進(jìn)行研究是一個(gè)比較好的選擇。入口效應(yīng)最多影響到第3層，出口效應(yīng)的影響相對(duì)不明顯，第4層結(jié)構(gòu)受出、入口效應(yīng)的影響已經(jīng)不大，而采取8層結(jié)構(gòu)的球陣并研究第4、5層球則能完全避免出、入口效應(yīng)。本文采用的模型為面接觸體心棱柱(BCP)排列，如圖1所示，圖中由左向右分為8層，左端為速度入口，右端為壓力出口，氦氣作為冷卻氣體從左端通入，不考慮重力。燃料球直徑為60 mm，接觸面設(shè)為直徑的1%。燃料球面材料為石墨，且球面設(shè)為熱源8 317 W·m。燃料球陣四周設(shè)為對(duì)稱面，不考慮壁面效應(yīng)帶來的影響。物性參數(shù)參考文獻(xiàn)［14］中的數(shù)據(jù)設(shè)置，具體如表1所示。

表1 氦氣與石墨的物性參數(shù)Tab. 1 Physical parameters of helium and graphite

圖1 面接觸體心棱柱(BCP)排列模型Fig. 1 Model of body-centered prism (BCP) by surface contact

在選擇湍流模型時(shí)，為了準(zhǔn)確捕捉流場中的渦結(jié)構(gòu)，本文采用大渦模擬(LES)方法。LES方法通過一個(gè)濾波函數(shù)對(duì)大渦直接求解，對(duì)小渦過濾，并建立亞格子應(yīng)力(SGS)模型。在眾多SGS模型中，本文選擇了自適應(yīng)更好的動(dòng)態(tài)Smagorinsky-Lilly模型。網(wǎng)格方面，本文使用Icem軟件，先尋找體心棱柱結(jié)構(gòu)的最小幾何單元(1/32)，然后將最小幾何單元?jiǎng)澐殖?1個(gè)六面體，通過點(diǎn)、線、面關(guān)聯(lián)構(gòu)建最小幾何單元的結(jié)構(gòu)化網(wǎng)格，并最終完成完整的面接觸體心棱柱球陣的結(jié)構(gòu)網(wǎng)格。BCP面接觸結(jié)構(gòu)網(wǎng)格如圖2所示。最終達(dá)成單層體心棱柱球陣131萬網(wǎng)格、8層1 050萬網(wǎng)格以滿足大渦模擬LES方法的計(jì)算要求。

圖2 BCP面接觸結(jié)構(gòu)網(wǎng)格Fig. 2 Structured grid of BCP surface contact

1.2 計(jì)算驗(yàn)證

為了驗(yàn)證網(wǎng)格與計(jì)算方法的適用性，本文選取與文獻(xiàn)[16]中同樣的中心球之間截線位置，并對(duì)距離與時(shí)均速度進(jìn)行了無量綱處理，對(duì)比結(jié)果如圖3所示。由于模型和計(jì)算方法的差異，截線中心區(qū)域的結(jié)果差別明顯，而且LES速度結(jié)果相比q-DNS速度結(jié)果更加對(duì)稱。除此之外，兩結(jié)果時(shí)均速度趨勢一致，在高速區(qū)域和近接觸面(點(diǎn))低速區(qū)域結(jié)果基本一致，由此可以說明本文所采用的計(jì)算方法與網(wǎng)格數(shù)量可以滿足計(jì)算要求。

圖3 兩種方法計(jì)算結(jié)果對(duì)比Fig. 3 Comparison of the calculation results by two methods

2 結(jié)果與討論

考慮到出、入口效應(yīng)的影響，本文選擇了中心位置第4層燃料球及其接觸燃料球作為主要的研究對(duì)象。同時(shí)，為了便于觀察各接觸面附近區(qū)域的速度場溫度場等細(xì)節(jié)，在體心棱柱球陣內(nèi)部流域的3個(gè)位置截取截面1～3，并在截面1上取線1～5，具體位置如圖4(a)、(b)所示。

圖4 各位置示意圖Fig. 4 Schematic diagram of each position

圖5為體心棱柱單層結(jié)構(gòu)的示意圖。根據(jù)接觸面位置的不同可以將接觸面分為：體心棱柱頂點(diǎn)燃料球相互接觸形成的接觸面1、體心棱柱同層中心燃料球相互接觸形成的接觸面2、體心棱柱頂部燃料球與中心燃料球相互接觸形成的接觸面3以及體心棱柱底部燃料球與中心燃料球相互接觸形成的接觸面4。

圖5 四種接觸面位置Fig. 5 Four locations of surface contact

2.1 BCP結(jié)構(gòu)面接觸與人工間隙處理方式對(duì)比

在處理規(guī)則排列球陣接觸面的問題上，人為拉開燃料球產(chǎn)生間隙來簡化問題的處理方式比較常見。本文將對(duì)面接觸與人工間隙處理方式作簡單對(duì)比分析。為了便于觀察，本文采用球坐標(biāo)系，代表天頂角，代表方位角。在球面取=0°和=45°兩處特殊位置曲線，如圖6所示。

圖6 球面取線位置Fig. 6 Line locations on the sphere surface

圖7為球面溫度對(duì)比。人工間隙處理方式由于不存在接觸，整體溫度比面接觸方式的低5～10 K左右。面接觸結(jié)構(gòu)中由于接觸面的存在，曲線并不連續(xù)，其中=0°曲線經(jīng)過接觸面3，=45°曲線經(jīng)過接觸面2、4。隨著氦氣的流動(dòng)，值由正變負(fù)，在接觸面曲線不連續(xù)的位置，值較大的位置稱為接觸面前沿，值較小的位置稱為接觸面后沿。不管是面接觸還是人工間隙處理方式，原接觸面位置必然會(huì)產(chǎn)生溫度升高的現(xiàn)象，而面接觸方式的溫差更為明顯。接觸面前沿過渡到接觸面后沿時(shí)溫度驟升，且接觸面后沿的溫度遠(yuǎn)大于同位置人工間隙處理方式的溫度。=0°這條曲線顯示氦氣流經(jīng)接觸面3時(shí)，對(duì)應(yīng)的高溫區(qū)域基本一致；而=45°的曲線則顯示出接觸面2、4處采用人工間隙處理方式時(shí)的溫度趨勢相比于面接觸的出現(xiàn)了延遲。當(dāng)面接觸方式的溫度由最高點(diǎn)下降后，人工間隙處理方式的溫度剛剛到達(dá)最高溫度，由此可見不同的接觸面影響各不相同。

圖7 球面溫度對(duì)比Fig. 7 Comparison of temperature on the sphere surface

總體來看，采用人工間隙處理方式的溫度分布更為均勻，但其高溫區(qū)出現(xiàn)的位置與實(shí)際高溫區(qū)域并不相符；面接觸處理方式是兩球真實(shí)接觸得出的結(jié)果，溫度的不均勻代表流場的復(fù)雜性。因此，面接觸方式更為合理。

2.2 BCP結(jié)構(gòu)流動(dòng)傳熱特性

BCP結(jié)構(gòu)內(nèi)部流場由于各種接觸面的存在十分復(fù)雜，本小節(jié)將從流域進(jìn)行分析。

圖8、9分別為流域內(nèi)各處截面時(shí)均與瞬時(shí)溫度和速度分布。觀察溫度結(jié)果可以發(fā)現(xiàn)，接觸面附近存在溫度高于其他流域溫度的區(qū)域，這些高溫區(qū)域集中于接觸面后沿。與時(shí)均結(jié)果相比，瞬時(shí)溫度分布更加混亂，呈現(xiàn)出一種不對(duì)稱性，局部位置瞬時(shí)溫度明顯高于時(shí)均值。觀察速度結(jié)果可以發(fā)現(xiàn)，在中心球頂部與底部存在低速區(qū)域，頂部區(qū)域?yàn)闇箙^(qū)，底部區(qū)域?yàn)槲擦鲄^(qū)。滯止區(qū)實(shí)際上流道逐漸變寬，氦氣沿著球面繼續(xù)流動(dòng)，頂部只流進(jìn)少量氦氣并附著停留；而尾流區(qū)流道逐漸收窄，同時(shí)因?yàn)榻咏佑|面后沿，流動(dòng)出現(xiàn)了分離，分離出來的氦氣與下一中心燃料球頂部滯止區(qū)留存的氦氣相互作用。這種相互作用導(dǎo)致氦氣的再附著，因此尾流區(qū)低速面積比滯止區(qū)的更大。對(duì)比速度場與溫度場可以看到，接觸面位置處速度低且始終對(duì)應(yīng)高溫，可見流道狹窄小的接觸面流動(dòng)換熱欠佳；而滯止區(qū)并未出現(xiàn)明顯的高溫區(qū)，可見速度分布不是影響溫度分布的唯一因素。

圖8 各截面時(shí)均與瞬時(shí)溫度分布Fig. 8 Time-averaged and instantaneous temperature distribution in each cross section

圖10為截面1取線位置時(shí)均速度。由圖中可以觀察到，在接觸面附近存在非常明顯的低速區(qū)域，而兩中心燃料球之間區(qū)域接近尾流區(qū)和滯止區(qū)，速度也較?。桓咚賲^(qū)域則代表未受接觸面影響的氦氣脫離接觸面、尾流區(qū)、滯止區(qū)以較快的速度進(jìn)入下一層BCP結(jié)構(gòu)。線1～5位置處的時(shí)均速度變化表明在從接觸面前沿向后沿過渡時(shí)氦氣速度不斷下降，此時(shí)流動(dòng)逐漸產(chǎn)生分離，這也是速度變小的根本原因。從前沿向后沿過渡時(shí)，接觸面附近會(huì)產(chǎn)生一個(gè)較小的速度波動(dòng)且隨著來流方向越來越明顯，這個(gè)速度波動(dòng)代表接觸面后沿產(chǎn)生了漩渦。

圖9 各截面時(shí)均與瞬時(shí)速度分布Fig. 9 Time-averaged and instantaneous velocity distribution in each cross section

圖10 截面1取線位置時(shí)均速度Fig. 10 Time-averaged velocity of the lines on cross section 1

2.3 BCP結(jié)構(gòu)中心球渦分析

為了判斷渦的存在，本文將采用Q準(zhǔn)則和渦量進(jìn)行對(duì)比分析，并選取更合適的判據(jù)。Q準(zhǔn)則由Hunt等于1988年提出，定義為

當(dāng)＞0時(shí)，表示流體微團(tuán)的旋轉(zhuǎn)克服了剪切應(yīng)力并在流動(dòng)中占主導(dǎo)作用，可以判斷流場中產(chǎn)生了漩渦。渦量通常用來描述渦旋強(qiáng)度，定義為流體速度矢量的旋度。

圖11為值與渦量的比較。從圖中可以看到，瞬時(shí)渦量圖和值圖都顯示了在多個(gè)接觸面互相影響下產(chǎn)生了眾多大小不一的渦，在不同接觸面之間互相影響的區(qū)域使用Q準(zhǔn)則能夠判斷識(shí)別出更多的渦，而使用渦量判斷則效果不明顯；渦量圖中頂部流經(jīng)中心接觸面、中心接觸面流經(jīng)底部的區(qū)域顯示有大量渦存在，而該區(qū)域?qū)嶋H上是氦氣沿著球面流動(dòng)，流動(dòng)分離并不劇烈，只會(huì)出現(xiàn)少量旋渦，可見渦量判據(jù)出現(xiàn)了誤判。誤判的根本原因還是在于渦量對(duì)于旋轉(zhuǎn)與剪切兩種流動(dòng)一視同仁，在多曲面問題上無法準(zhǔn)確判斷旋渦，這也是渦量的局限性所在。因此，在后續(xù)分析中本文將采用Q準(zhǔn)則判斷渦區(qū)域。

圖11 Q值與渦量比較Fig. 11 Comparison of Q value and vorticity

圖12、13為中心燃料球表面速度場與溫度場的時(shí)均結(jié)果與瞬時(shí)結(jié)果。在燃料球頂部和底部可見兩塊低速區(qū)域，該低速區(qū)域的形成是因接觸面1的存在擠壓了流道，阻礙了來流氦氣的流動(dòng)。所有接觸面后沿都形成局部高溫?zé)岚撸佑|面前沿則存在低溫區(qū)域，這是因?yàn)榻佑|面后沿流道較為狹小，氦氣流經(jīng)此處無法充分換熱，而前沿處于流道較寬闊地帶，熱量能夠被來流氦氣帶走。頂部接觸面3與底部接觸面4之間的垂直區(qū)域是低速高溫區(qū)且變化不大，而接觸面2與接觸面3之間的三角區(qū)域、接觸面2與接觸面4之間的三角區(qū)域有較大的溫度和速度變化。因此，基本可以判斷中心接觸面2的存在直接改變了氦氣流向，從而影響了對(duì)流換熱。

圖12 中心燃料球時(shí)均溫度與燃料球外0.6 mm處時(shí)均速度分布Fig. 12 Time-averaged temperature of central fuel ball and time-averaged velocity distribution at 0.6 mm outside of fuel ball

圖13 中心燃料球瞬時(shí)溫度與燃料球外0.6 mm處瞬時(shí)速度分布Fig. 13 Instantaneous temperature of central fuel ball and instantaneous velocity distribution at 0.6 mm outside of fuel ball

圖14為燃料球外0.6 mm處時(shí)均值分布。由圖中可以看到，渦強(qiáng)度較高的區(qū)域集中在接觸面2和接觸面3之間的三角區(qū)域，而燃料球頂部接觸面2與底部接觸面4處產(chǎn)生的渦強(qiáng)度較小。結(jié)合圖12時(shí)均速度和時(shí)均溫度分布可以發(fā)現(xiàn)，在高速區(qū)域渦的強(qiáng)度很大，而接觸面后沿幾乎無渦，速度幾乎為零，低溫區(qū)域則完全涵蓋了存在渦的區(qū)域。由此說明：低速區(qū)域基本不會(huì)產(chǎn)生渦，高速區(qū)域必然會(huì)產(chǎn)生大量的渦；渦的存在有利于氦氣的流動(dòng)，高速氦氣帶走了熱量從而降低了球面溫度。

從能量方程對(duì)流項(xiàng)可知，除了速度大小，速度矢量與溫度梯度夾角也決定了速度與溫度梯度點(diǎn)乘的大小，且夾角越小對(duì)流換熱越好。圖15為燃料球外0.6 mm處時(shí)均速度與溫度梯度夾角分布。結(jié)合圖14、15可以發(fā)現(xiàn)，在頂部與底部渦分布較少的區(qū)域，時(shí)均夾角普遍大于90°，部分區(qū)域甚至達(dá)到120°，過大的時(shí)均夾角對(duì)換熱不利；而有漩渦產(chǎn)生的區(qū)域時(shí)均夾角普遍小于等于45°，此時(shí)漩渦加強(qiáng)了對(duì)流，增強(qiáng)了換熱。在接觸面2與接觸面3、接觸面2與接觸面4共同作用的區(qū)域，大量渦的存在使得時(shí)均夾角分布不均，此處漩渦不如頂部與底部漩渦的作用明顯。

圖14 燃料球外0.6 mm處時(shí)均Q值分布Fig. 14 Time-averaged Q value distribution at 0.6 mm outside of fuel ball

圖15 燃料球外0.6 mm處時(shí)均速度與溫度梯度夾角分布Fig. 15 Distribution of angles between time-averaged velocity and temperature gradient at 0.6 mm outside of fuel ball

2.4 BCP結(jié)構(gòu)與BCC、FCC結(jié)構(gòu)對(duì)比

BCP結(jié)構(gòu)與常見的BCC結(jié)構(gòu)都為體心結(jié)構(gòu)，但是前者接觸位置更多，結(jié)構(gòu)更為緊湊；BCP結(jié)構(gòu)與FCC結(jié)構(gòu)相比，兩者均具有相同的空隙率，但幾何構(gòu)成卻大不相同。參考文獻(xiàn)［9］、［10］，本文對(duì)3種結(jié)構(gòu)做一簡單對(duì)比。

圖16為3種規(guī)則排列球陣在相同工況、相同位置溫度場的差異。觀察軸向截面結(jié)果，BCC結(jié)構(gòu)體心球與球之間距離近，從而導(dǎo)致球間流域溫度高于BCP和FCC結(jié)構(gòu)；FCC結(jié)構(gòu)由于其面心排布，面心球之間距離和BCP結(jié)構(gòu)體心球之間距離相近，因而在球間流域有相似的結(jié)果，氦氣流動(dòng)換熱效果較好；BCP結(jié)構(gòu)頂部燃料球之間互相接觸產(chǎn)生的局部高溫?zé)狳c(diǎn)使整體溫度分布更不均勻，而這種接觸在另外兩種結(jié)構(gòu)中并不存在，這也是BCP結(jié)構(gòu)最大的特點(diǎn)。再來觀察徑向截面結(jié)果，BCC結(jié)構(gòu)中心球與其他中心球之間并無接觸，間隙處存在4塊低溫區(qū)域；而FCC結(jié)構(gòu)與BCP結(jié)構(gòu)則存在接觸，并產(chǎn)生了高溫?zé)狳c(diǎn)，熱點(diǎn)后方則有8小塊低溫區(qū)域，這些低溫區(qū)域代表氦氣沿球面流動(dòng)而未發(fā)生流動(dòng)分離。接觸面的存在打斷了FCC和BCP結(jié)構(gòu)球面的連續(xù)性，故而使流動(dòng)更加復(fù)雜。

圖16 3種結(jié)構(gòu)的時(shí)均溫度對(duì)比Fig. 16 Comparison of time-averaged temperature for three models

3 結(jié) 論

本文采用大渦模擬方法，對(duì)比了面接觸與人工間隙處理方式的差異，對(duì)面接觸BCP結(jié)構(gòu)的速度場與溫度場進(jìn)行求解，得到流場中的瞬時(shí)和時(shí)均速度和溫度分布，并了解了其流動(dòng)換熱的特點(diǎn)。同時(shí)充分發(fā)揮LES方法的優(yōu)越性，對(duì)流場中的渦進(jìn)行了分析。最后得到的主要結(jié)論為：

(1) 人工間隙處理方式對(duì)于高溫區(qū)域存在誤判，采用面接觸方式更為合理。在接觸面后沿出現(xiàn)顯著高溫?zé)狳c(diǎn)，且隨著來流方向接觸面越多，溫度越高。在燃料球頂部與底部滯止區(qū)與尾流區(qū)氦氣分離再附著造成局部高溫區(qū)域，且滯止區(qū)高溫面積小于尾流區(qū)。

(2) BCP排列球陣四個(gè)不同位置的接觸面構(gòu)成的三角形區(qū)域存在大量附著與分離流動(dòng)共存的現(xiàn)象。

(3) BCP結(jié)構(gòu)頂部球之間的互相接觸使得流域被多種接觸面分割，與BCC和FCC結(jié)構(gòu)相比，其溫度分布更為復(fù)雜。

(4) 渦在接觸面附近、接觸面三角區(qū)域、接觸面垂直區(qū)域會(huì)影響溫度分布，且從整體來看渦的存在對(duì)流動(dòng)換熱是有利的。