陳 波，步珊珊，陳德奇，馬在勇，張盧騰

(重慶大學(xué) 低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，重慶 400044)

高溫氣冷堆是國(guó)際領(lǐng)先的第4代核能技術(shù)，具有經(jīng)濟(jì)性和固有安全性等特點(diǎn)[1]。高溫氣冷球床堆采用球形全陶瓷包覆顆粒作為燃料元件，其安全設(shè)計(jì)準(zhǔn)則要求：在任何事故情況下，球床堆芯的衰變余熱依靠輻射換熱和導(dǎo)熱等非能動(dòng)傳熱機(jī)制導(dǎo)出，實(shí)現(xiàn)事故停堆后的冷卻，以保證燃料包殼不發(fā)生破裂，將放射性物質(zhì)限制在內(nèi)[2]。高溫氣冷球床堆堆芯內(nèi)的熱量傳遞機(jī)制十分復(fù)雜，通常將球床內(nèi)的導(dǎo)熱、輻射換熱等復(fù)雜傳熱機(jī)制等效為簡(jiǎn)單的熱傳導(dǎo)過程，采用一個(gè)綜合的有效導(dǎo)熱系數(shù)來表征。因此有效導(dǎo)熱系數(shù)是表征球床堆芯宏觀傳熱能力的關(guān)鍵參數(shù)[3]，也是高溫氣冷球床堆熱工設(shè)計(jì)與安全分析程序中的基本參數(shù)。

事故工況下，當(dāng)冷卻劑失去強(qiáng)制流動(dòng)后，對(duì)流傳熱不考慮，高溫氣冷球床堆堆芯的有效導(dǎo)熱系數(shù)主要由3部分組成[4]：1) 固體表面間的輻射換熱；2) 球體顆粒的導(dǎo)熱和流體的導(dǎo)熱，即氣-固導(dǎo)熱；3) 相鄰顆粒接觸點(diǎn)之間的接觸導(dǎo)熱。Zehner等[5]基于圓柱單元方法提出了計(jì)算球床有效導(dǎo)熱系數(shù)的ZS模型，這一模型可很好地計(jì)算球床的氣-固導(dǎo)熱[6]，但未考慮接觸導(dǎo)熱和輻射換熱。Bauer等[7]采用了接觸面積系數(shù)φ來考慮接觸導(dǎo)熱，同時(shí)使用Breitbach等[8]模型來計(jì)算輻射換熱，這一改進(jìn)后的模型被稱為ZBS模型。ZBS模型被廣泛應(yīng)用于預(yù)測(cè)球床有效導(dǎo)熱系數(shù)的計(jì)算。當(dāng)溫度低于1 000 ℃ 時(shí)，Wu等[9]發(fā)現(xiàn)ZBS模型可很好地預(yù)測(cè)德國(guó)SANA實(shí)驗(yàn)和南非HTTU實(shí)驗(yàn)的測(cè)量結(jié)果。de Beer等[10]發(fā)現(xiàn)在球床主體區(qū)域，ZBS模型預(yù)測(cè)值略低于實(shí)驗(yàn)值。Bu等[11]的實(shí)驗(yàn)結(jié)果也表明ZBS模型對(duì)球床有效導(dǎo)熱系數(shù)的預(yù)測(cè)性能要優(yōu)于其他模型。然而，ZBS模型預(yù)測(cè)值對(duì)接觸面積系數(shù)φ十分敏感[12]，一般通過實(shí)驗(yàn)值來調(diào)整φ的取值，如Antwerpen等[13]發(fā)現(xiàn)當(dāng)φ=0.01時(shí)，ZBS模型與實(shí)驗(yàn)結(jié)果吻合較好。Bu等[12]發(fā)現(xiàn)，當(dāng)φ=0.07時(shí)ZBS模型預(yù)測(cè)值和簡(jiǎn)單立方球床的實(shí)驗(yàn)測(cè)量結(jié)果的平均相對(duì)誤差小于5%。目前關(guān)于接觸面積系數(shù)φ的計(jì)算表達(dá)式較少，只有Chen等[14]和You等[15]分別給出了φ的計(jì)算表達(dá)式。

在前期工作的基礎(chǔ)上[6,11-12,16]，本文數(shù)值模擬不同球床結(jié)構(gòu)的有效導(dǎo)熱系數(shù)，通過多元線性回歸獲得接觸面積系數(shù)φ的計(jì)算表達(dá)式。通過與SANA實(shí)驗(yàn)結(jié)果及前期實(shí)驗(yàn)結(jié)果的對(duì)比，評(píng)估改進(jìn)后的ZBS模型的預(yù)測(cè)能力。

1 物理模型和數(shù)值方法

1.1 ZBS模型

Bauer等[7]基于圓柱單位元提出并聯(lián)的3條傳熱路徑，并采用分布系數(shù)加權(quán)疊加接觸導(dǎo)熱、氣固導(dǎo)熱和輻射換熱來預(yù)測(cè)球床結(jié)構(gòu)的有效導(dǎo)熱系數(shù)，這種理論模型被稱為ZBS模型，由于該模型具有良好的預(yù)測(cè)性而被廣泛關(guān)注。ZBS模型如下：

(1)

其中：keff為有效導(dǎo)熱系數(shù)，W/(m·K)；kf為流體導(dǎo)熱系數(shù)，W/(m·K)；ε為孔隙率；κ為氣-固導(dǎo)熱系數(shù)比；κr為輻射換熱系數(shù)比；κG為Knudsen狀態(tài)下的氣體導(dǎo)熱系數(shù)比，這里取1；κc為綜合導(dǎo)熱系數(shù)比，計(jì)算方法見文獻(xiàn)[4]；φ為接觸面積系數(shù)，是一個(gè)表征接觸導(dǎo)熱份額的無量綱系數(shù)，可表示為：

φ=aNcm·γn

(2)

其中：Nc為配位數(shù)；γ為接觸直徑比，γ=dc/dp，dc和dp分別為接觸直徑和球徑。系數(shù)a和指數(shù)m、n通過多元線性回歸方法求解，因此需要獲得不同配位數(shù)Nc及接觸直徑比γ對(duì)應(yīng)的接觸面積系數(shù)φ的值。

1.2 數(shù)值計(jì)算模型

為獲得多元線性回歸需要的數(shù)據(jù)，本文針對(duì)不同幾何參數(shù)的球床有效導(dǎo)熱系數(shù)進(jìn)行數(shù)值分析。如圖1所示，采用簡(jiǎn)單立方(SC)、體心立方(BCC)、面心立方(FCC)和無序4種堆積結(jié)構(gòu)單元，這4種結(jié)構(gòu)的配位數(shù)Nc分別為6、8、12和5.5。無序堆積結(jié)構(gòu)單元的生成方法為：基于PFC3D軟件和膨脹法生成直徑為15dp、高度為15dp的大球床，然后從大球床的中心區(qū)域中抽取出來2.5dp×2.4dp×2.6dp的結(jié)構(gòu)單元。每種堆積結(jié)構(gòu)單元分別構(gòu)建了0.035、0.050和0.100 3種不同的接觸直徑比γ，因此有12組不同配位數(shù)Nc及接觸直徑比γ對(duì)應(yīng)的球床結(jié)構(gòu)模型(表1)。

結(jié)構(gòu)單元：a——SC；b——BCC；c——FCC；d——無序

表1 不同幾何模型的參數(shù)

1.3 數(shù)值方法

本文計(jì)算中流體無流動(dòng)，球床內(nèi)無內(nèi)熱源，只需求解固體區(qū)域和流體區(qū)域的能量方程。流體和固體區(qū)域的能量方程為：

(3)

其中：ki為固體(石墨)和流體(氮?dú)?氦氣)的導(dǎo)熱系數(shù)，W/(m·K)，并隨溫度變化；i=s、f，分別代表固相和氣相；T為溫度；x、y、z表示空間坐標(biāo)。

獲得球床溫度分布后，通過逆求解一維傅里葉導(dǎo)熱定律獲得有效導(dǎo)熱系數(shù)：

(4)

其中：q為平均熱流密度，W/m2；ΔT為高溫壁面和低溫壁面的溫差，10 K；Δz為高溫壁面和低溫壁面的距離，m。

模型1的計(jì)算區(qū)域和邊界條件如圖2所示。z方向前后的高溫壁面和低溫壁面是定溫邊界條件；堆積單元的平均溫度范圍為373～1 073 K；垂直于xy平面的另外4個(gè)表面是絕熱的，固體顆粒和氮?dú)?氦氣的交界面無滑移。將以上模型導(dǎo)入ANSYS Fluent 15.0中計(jì)算，求解器基于壓力求解，求解算法選用SIMPLE算法，離散格式采用二階迎風(fēng)格式，采用殘差小于10-12作為收斂判據(jù)。本文數(shù)值計(jì)算模型的驗(yàn)證詳見文獻(xiàn)[6]。

圖2 模型1計(jì)算區(qū)域和邊界條件

1.4 網(wǎng)格獨(dú)立性驗(yàn)證

采用ICEM軟件對(duì)幾何模型劃分了3套非結(jié)構(gòu)化網(wǎng)格，對(duì)模型1進(jìn)行網(wǎng)格無關(guān)性考核，結(jié)果列于表2。和第3套網(wǎng)格的計(jì)算結(jié)果相比，第1和第2套網(wǎng)格計(jì)算得到的有效導(dǎo)熱系數(shù)相對(duì)偏差分別為1.43%和0.29%，因此可認(rèn)為第2套網(wǎng)格獲得網(wǎng)格無關(guān)解。以無序球床結(jié)構(gòu)單元為例，模型10網(wǎng)格劃分示意圖如圖3所示，接觸區(qū)域和顆粒表面均進(jìn)行了適當(dāng)?shù)募用堋?/p>

表2 模型1網(wǎng)格獨(dú)立性驗(yàn)證

圖3 模型10網(wǎng)格劃分示意圖

2 接觸面積系數(shù)的多元線性回歸

通過調(diào)整接觸面積系數(shù)φ的值，使得ZBS模型的預(yù)測(cè)值與數(shù)值計(jì)算結(jié)果的平均相對(duì)誤差最小，來獲得不同配位數(shù)Nc及接觸直徑比γ對(duì)應(yīng)的接觸面積系數(shù)φ的值。以模型1為例，如圖4所示，ZBS模型預(yù)測(cè)值對(duì)φ十分敏感。隨著φ的取值減小，ZBS的預(yù)測(cè)值也減小。如圖5所示，不同φ對(duì)應(yīng)的預(yù)測(cè)值和數(shù)值計(jì)算結(jié)果的平均相對(duì)誤差先減小后增大。平均相對(duì)誤差的絕對(duì)值最小時(shí)對(duì)應(yīng)的φ就是ZBS模型中φ的最優(yōu)取值，因此可獲得模型1對(duì)應(yīng)的φ為0.105。通過同樣的分析方法獲得了12組不同配位數(shù)Nc及接觸直徑比γ對(duì)應(yīng)的φ，結(jié)果列于表3。

圖4 模型1模擬結(jié)果與ZBS模型預(yù)測(cè)值的對(duì)比

圖5 相對(duì)誤差的絕對(duì)值分布

表3 不同配位數(shù)Nc和接觸直徑比γ對(duì)應(yīng)下的φ

為進(jìn)行多元線性回歸，式(1)變換為下面的形式：

lnφ=lna+mlnNc+nlnγ

(5)

然后采用最小二乘法對(duì)表3中的數(shù)據(jù)進(jìn)行擬合，獲得φ的表達(dá)式為：

φ=0.21Ncγ

(6)

對(duì)于真實(shí)的無序球床結(jié)構(gòu)，Nc可通過以下表達(dá)式[17]計(jì)算：

Nc=25.952ε3-62.364ε2+

39.724ε-2.023 3 0.26≤ε≤0.48

(7)

接觸直徑比γ=dc/dp，dc可根據(jù)下式[18]計(jì)算：

(8)

其中：μp為泊松比；Ep為楊氏彈性模量，Pa；F為外加作用力，N；rp為球的半徑。因此接觸直徑比γ也可求解。確定了球床結(jié)構(gòu)的Nc及γ，ZBS模型中的φ可根據(jù)式(6)求解獲得。式(6)適用于孔隙率為0.260～0.48、溫度范圍為273～1 100 K的球床結(jié)構(gòu)。

3 改進(jìn)后的ZBS模型預(yù)測(cè)能力評(píng)估

采用SANA實(shí)驗(yàn)結(jié)果[19]評(píng)估改進(jìn)后的ZBS模型預(yù)測(cè)能力。SANA實(shí)驗(yàn)裝置針對(duì)氮?dú)?氦氣氛圍下的球床等效導(dǎo)熱系數(shù)進(jìn)行了測(cè)量，石墨球床的顆粒直徑dp為60 mm，圓筒堆積球床直徑為25dp，高度為1 m，孔隙率為0.39，測(cè)試溫度范圍為273～1 223 K。改進(jìn)后的ZBS模型和SANA實(shí)驗(yàn)結(jié)果對(duì)比如圖6所示。由圖6可見，改進(jìn)后的ZBS模型預(yù)測(cè)能力在中低溫工況下(小于1 000 K)要優(yōu)于其他模型(包括IAEA ZS模型[19]、MSUC模型[13]、Chen模型[14]和You模型[15])。在高溫工況下(大于1 000 K)，改進(jìn)后的ZBS模型預(yù)測(cè)能力和IAEA ZS模型[19]和MSUC模型[13]相當(dāng)。當(dāng)溫度大于1 000 K左右，球床主導(dǎo)的傳熱機(jī)制從接觸導(dǎo)熱轉(zhuǎn)變?yōu)檩椛鋫鳠釞C(jī)制[20]。Chen模型和You模型均是基于ZBS模型的改進(jìn)，其中Chen模型可預(yù)測(cè)兩種不同直徑顆粒堆積成的二元球床的有效導(dǎo)熱系數(shù)，You模型改進(jìn)了ZBS模型中的輻射傳熱模型。由于Chen模型和You模型中提出的接觸面積系數(shù)φ預(yù)測(cè)值(分別為0.019和0.003)要比本文的計(jì)算值(0.033)小很多，低估了接觸導(dǎo)熱，因此在中低溫工況下對(duì)球床有效導(dǎo)熱系數(shù)的預(yù)測(cè)值偏小。IAEA ZS模型中將氣-固導(dǎo)熱、接觸導(dǎo)熱和輻射傳熱機(jī)制對(duì)應(yīng)的傳熱系數(shù)直接疊加，其中接觸導(dǎo)熱部分采用Chen-Tien模型[21]計(jì)算。MSUC模型中接觸導(dǎo)熱系數(shù)的計(jì)算方法與IAEA ZS模型相同。綜合以上分析，改進(jìn)后的ZBS模型對(duì)于接觸導(dǎo)熱的貢獻(xiàn)有較好地預(yù)測(cè)能力，這說明本文提出的φ的計(jì)算表達(dá)式有較高的可靠性。值得注意的是，目前所有模型的預(yù)測(cè)值均比SANA實(shí)驗(yàn)結(jié)果小。原因可能是SANA實(shí)驗(yàn)中自然對(duì)流效應(yīng)比較顯著，而以上所有模型都未考慮自然對(duì)流的影響。

圖6 改進(jìn)后的ZBS模型與SANA實(shí)驗(yàn)結(jié)果對(duì)比

本課題組前期開展了球床有效導(dǎo)熱系數(shù)實(shí)驗(yàn)測(cè)量[11，16]，實(shí)驗(yàn)系統(tǒng)主要包括長(zhǎng)方體球床堆芯、氮?dú)夤?yīng)系統(tǒng)、加熱電源與控制系統(tǒng)、冷卻水系統(tǒng)和溫度數(shù)據(jù)采集系統(tǒng)5部分。如圖7所示，直流電加熱板與石墨球床的高溫壁面接觸，球床的4個(gè)側(cè)壁面均包裹保溫層，盡量減少堆芯散熱，使熱量主要向另外一側(cè)的低溫壁面?zhèn)鬟f，然后被低溫壁面?zhèn)鹊睦鋮s水帶走熱量。在實(shí)驗(yàn)測(cè)量中，壓力控制在40 kPa，自然對(duì)流效應(yīng)被大大抑制。石墨球床的顆粒直徑dp為15 mm，長(zhǎng)方體堆積球床為6dp×6dp×14dp，孔隙率為0.40，實(shí)驗(yàn)中高溫壁面溫度范圍為380～800 K。如圖8所示，在球床主體區(qū)域(距壁面大于4個(gè)顆粒直徑)，改進(jìn)后的ZBS模型和實(shí)驗(yàn)結(jié)果吻合得很好，最大相對(duì)誤差只有12.6%。但在近壁面區(qū)域(距壁面4個(gè)顆粒直徑以內(nèi))和實(shí)驗(yàn)值誤差較大，這說明改進(jìn)后的ZBS模型在球床近壁面區(qū)域的預(yù)測(cè)能力較差。由于壁面限制，球床在近壁面區(qū)域的幾何結(jié)構(gòu)與主體區(qū)域有顯著不同，如孔隙率在近壁面區(qū)域大同時(shí)波動(dòng)強(qiáng)烈(壁面效應(yīng))，因此近壁面區(qū)域有效導(dǎo)熱系數(shù)會(huì)明顯小于主體區(qū)域的有效導(dǎo)熱系數(shù)，但ZBS模型(包括基于ZBS模型的改進(jìn)模型)的推導(dǎo)未考慮到這種壁面效應(yīng)，因此ZBS模型在球床近壁面區(qū)域的預(yù)測(cè)能力較差。

圖7 實(shí)驗(yàn)中球床示意圖

圖8 改進(jìn)后的ZBS模型和前期實(shí)驗(yàn)結(jié)果對(duì)比

4 結(jié)論

為提高ZBS模型對(duì)球床結(jié)構(gòu)有效導(dǎo)熱系數(shù)的預(yù)測(cè)能力，本文數(shù)值模擬了不同結(jié)構(gòu)球床的有效導(dǎo)熱系數(shù)，通過多元線性回歸獲得了ZBS模型中經(jīng)驗(yàn)參數(shù)接觸面積系數(shù)φ的計(jì)算表達(dá)式。與SANA實(shí)驗(yàn)結(jié)果對(duì)比表明，改進(jìn)后的ZBS模型的預(yù)測(cè)能力優(yōu)于其他球床有效導(dǎo)熱系數(shù)模型。和前期的實(shí)驗(yàn)結(jié)果對(duì)比，改進(jìn)后的ZBS模型的計(jì)算結(jié)果在球床主體區(qū)域吻合很好，但在近壁面區(qū)域吻合較差。因此，改進(jìn)后的ZBS模型可用于孔隙率為0.26～0.48、溫度在273～1 100 K之間的球床結(jié)構(gòu)主體區(qū)域有效導(dǎo)熱系數(shù)預(yù)測(cè)。