高溫氣冷球床堆有效導熱系數模型的改進

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

(重慶大學 低品位能源利用技術及系統教育部重點實驗室，重慶 400044)

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

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

在前期工作的基礎上[6,11-12,16]，本文數值模擬不同球床結構的有效導熱系數，通過多元線性回歸獲得接觸面積系數φ的計算表達式。通過與SANA實驗結果及前期實驗結果的對比，評估改進后的ZBS模型的預測能力。

1 物理模型和數值方法

1.1 ZBS模型

Bauer等[7]基于圓柱單位元提出并聯的3條傳熱路徑，并采用分布系數加權疊加接觸導熱、氣固導熱和輻射換熱來預測球床結構的有效導熱系數，這種理論模型被稱為ZBS模型，由于該模型具有良好的預測性而被廣泛關注。ZBS模型如下：

(1)

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

φ=aNcm·γn

(2)

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

1.2 數值計算模型

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

結構單元：a——SC；b——BCC；c——FCC；d——無序

表1 不同幾何模型的參數

1.3 數值方法

本文計算中流體無流動，球床內無內熱源，只需求解固體區域和流體區域的能量方程。流體和固體區域的能量方程為：

(3)

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

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

(4)

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

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

圖2 模型1計算區域和邊界條件

1.4 網格獨立性驗證

采用ICEM軟件對幾何模型劃分了3套非結構化網格，對模型1進行網格無關性考核，結果列于表2。和第3套網格的計算結果相比，第1和第2套網格計算得到的有效導熱系數相對偏差分別為1.43%和0.29%，因此可認為第2套網格獲得網格無關解。以無序球床結構單元為例，模型10網格劃分示意圖如圖3所示，接觸區域和顆粒表面均進行了適當的加密。

表2 模型1網格獨立性驗證

圖3 模型10網格劃分示意圖

2 接觸面積系數的多元線性回歸

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

圖4 模型1模擬結果與ZBS模型預測值的對比

圖5 相對誤差的絕對值分布

表3 不同配位數Nc和接觸直徑比γ對應下的φ

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

lnφ=lna+mlnNc+nlnγ

(5)

然后采用最小二乘法對表3中的數據進行擬合，獲得φ的表達式為：

φ=0.21Ncγ

(6)

對于真實的無序球床結構，Nc可通過以下表達式[17]計算：

Nc=25.952ε3-62.364ε2+

39.724ε-2.023 3 0.26≤ε≤0.48

(7)

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

(8)

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

3 改進后的ZBS模型預測能力評估

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

圖6 改進后的ZBS模型與SANA實驗結果對比

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

圖7 實驗中球床示意圖

圖8 改進后的ZBS模型和前期實驗結果對比

4 結論

為提高ZBS模型對球床結構有效導熱系數的預測能力，本文數值模擬了不同結構球床的有效導熱系數，通過多元線性回歸獲得了ZBS模型中經驗參數接觸面積系數φ的計算表達式。與SANA實驗結果對比表明，改進后的ZBS模型的預測能力優于其他球床有效導熱系數模型。和前期的實驗結果對比，改進后的ZBS模型的計算結果在球床主體區域吻合很好，但在近壁面區域吻合較差。因此，改進后的ZBS模型可用于孔隙率為0.26～0.48、溫度在273～1 100 K之間的球床結構主體區域有效導熱系數預測。