壓水堆堆芯熔融物的運(yùn)動(dòng)及傳熱數(shù)值分析

丁 碩，黃海富，鐘汝浩，張小英,*，展德奎

(1.中山大學(xué) 中法核工程與技術(shù)學(xué)院，廣東 珠海 519082；2.中廣核研究院有限公司，廣東 深圳 518031)

嚴(yán)重事故下堆芯燃料的升溫熔融及遷移過(guò)程一直是一個(gè)重要的研究課題。在核電站嚴(yán)重事故中，堆芯由于缺少冷卻劑將持續(xù)升溫，極可能發(fā)生熔融現(xiàn)象。熔融的燃料棒單元會(huì)沿下方未熔融棒的壁面滴落，高溫熔融物接觸到溫度較低的壁面時(shí)，由于熱交換，溫度降到凝固溫度，會(huì)形成凝固殼堵塞局部通道。當(dāng)流道被完全堵塞時(shí)，堆芯冷卻劑無(wú)法進(jìn)入流道，最終導(dǎo)致堆芯由于缺乏有效冷卻而熔融坍塌。堆芯熔融物的堆內(nèi)運(yùn)動(dòng)與傳熱瞬態(tài)分析對(duì)于預(yù)測(cè)凝固堵塞位置、評(píng)估堆芯熔毀狀況及制定事故緩解措施有著重要意義。

由于受到實(shí)驗(yàn)條件的制約，目前國(guó)際上針對(duì)嚴(yán)重事故下堆芯熔融的實(shí)驗(yàn)很少，研究者主要借助仿真軟件來(lái)模擬嚴(yán)重事故下的熔堆進(jìn)程和現(xiàn)象，并對(duì)其后果進(jìn)行客觀評(píng)價(jià)。MAAP[1]是一款用于模擬嚴(yán)重事故現(xiàn)象及進(jìn)程的一體化、模塊化的輕水堆核電站嚴(yán)重事故分析程序。它可模擬燃料升溫，堆芯材料熔化及遷移，裂變產(chǎn)物的釋放、傳輸和沉積等過(guò)程及現(xiàn)象，但空間節(jié)點(diǎn)較粗糙，且節(jié)點(diǎn)劃分是相對(duì)固定的，對(duì)用戶(hù)而言缺乏可操作性。MELCOR[2]是全集成封裝的工程計(jì)算程序，它利用粗網(wǎng)格大步長(zhǎng)，使用了大量的參數(shù)化建模方法，使得不同計(jì)算階段的計(jì)算數(shù)據(jù)可在程序的不同模塊之間傳遞，但計(jì)算結(jié)果不確定，在各關(guān)鍵環(huán)節(jié)，中間結(jié)果的正確性需要用戶(hù)手動(dòng)檢查[3]。SCDAP/RELAP5程序[4]主要由RELAP5和SCDAP耦合而成，不僅繼承了RELAP5在熱工水力計(jì)算方面的優(yōu)點(diǎn)，且并入了SCDAP的嚴(yán)重事故現(xiàn)象分析模型。這些程序雖能模擬堆芯材料的熔融和遷移，但節(jié)點(diǎn)都很粗糙，無(wú)法精細(xì)計(jì)算熔融物質(zhì)的堆內(nèi)運(yùn)動(dòng)及傳熱過(guò)程。

因此，針對(duì)HPR1000堆型嚴(yán)重事故下反應(yīng)堆堆芯熔融坍塌問(wèn)題，在獲取嚴(yán)重事故工況下堆芯參數(shù)后，本文針對(duì)堆芯建立精細(xì)、符合實(shí)際的三維幾何模型，應(yīng)用Fortran語(yǔ)言通過(guò)數(shù)值模擬方法模擬嚴(yán)重事故工況下壓水堆堆芯升溫及堆芯熔融情況。

1 堆芯幾何模型

根據(jù)HPR1000堆型堆芯功率和溫度分布的對(duì)稱(chēng)性，選取如圖1所示的1/4堆芯作為研究對(duì)象，其中包括52個(gè)組件，289根燃料棒和控制棒。模擬中采用的物性參數(shù)列于表1。由于控制棒主要由金屬材料構(gòu)成，認(rèn)為在燃料棒達(dá)到熔融狀態(tài)前，所有控制棒已完全熔化。選取18號(hào)組件中任一燃料棒，事故發(fā)生后其燃料芯塊軸向中心外表面溫度變化如圖2所示。當(dāng)芯塊外表面溫度達(dá)到鋯包殼熔點(diǎn)時(shí)刻，約800 s后溫度升至芯塊熔點(diǎn)。同時(shí)鋯水反應(yīng)釋放的化學(xué)熱將加速鋯包殼升溫進(jìn)程[5]，故在燃料芯塊發(fā)生熔融前，鋯包殼全部熔融并有充分的時(shí)間掉落至堆芯底部。

圖1 1/4堆芯的組件結(jié)構(gòu)Fig.1 1/4 core cross-sectional shape

將燃料棒考慮為長(zhǎng)3 675 mm、直徑8.19 mm的實(shí)心圓棒，為精細(xì)計(jì)算燃料棒單元的瞬態(tài)運(yùn)動(dòng)與傳熱過(guò)程，將每根燃料棒在軸向上等分為60個(gè)單元，每個(gè)單元在軸向上又被等分為6層，在徑向上被劃分為等徑差的10格，如圖3所示。在計(jì)算過(guò)程中認(rèn)為未熔融的燃料棒單元軸向溫度相同，徑向存在溫差，而熔融的燃料棒單元每層具有各自的溫度。當(dāng)燃料棒的溫度達(dá)到熔點(diǎn)時(shí)，熔融的物質(zhì)會(huì)離開(kāi)原始位置。由流體力學(xué)可知，6層熔融單元將呈環(huán)狀一起沿下方固體燃料棒壁面在相應(yīng)的流道中下落，因此將方形流道進(jìn)行圓形等效處理。同時(shí)，熔融單元位置上方的固體單元由于熔融物離開(kāi)初始位置而缺少支持力，將在重力作用下下落，并與下方未熔融的燃料棒單元相連。

表1 堆芯材料物性參數(shù)Table 1 Material property of core component

圖2 軸向中心外表面溫度的變化Fig.2 Change of external surface temperature of axial center

圖3 燃料棒單元的軸向與徑向劃分模型Fig.3 Fuel rod model in axial and radial directions

獲取嚴(yán)重事故下堆芯的初始溫度和功率分布后，使用時(shí)間推進(jìn)方法構(gòu)建傳熱與運(yùn)動(dòng)的耦合模型。取瞬態(tài)時(shí)間步長(zhǎng)為1 s，且將每個(gè)時(shí)間步長(zhǎng)細(xì)分為6個(gè)子時(shí)間步。每經(jīng)過(guò)1個(gè)子時(shí)間步，程序?qū)⒂?jì)算各燃料棒單元的溫度，并對(duì)每個(gè)運(yùn)動(dòng)的熔融單元求解一次運(yùn)動(dòng)和傳熱方程，以確定每個(gè)時(shí)刻熔融燃料棒單元的位置與溫度信息，同時(shí)每1 s輸出1次計(jì)算結(jié)果。

2 數(shù)值模型

2.1 熔融物瞬態(tài)運(yùn)動(dòng)模型

當(dāng)燃料棒單元的溫度達(dá)到熔點(diǎn)時(shí)，熔融的燃料棒單元將離開(kāi)原來(lái)的位置，在對(duì)應(yīng)的流道中下落。在下落過(guò)程中，為計(jì)算高溫熔融物單元的位置，考慮熔融物的重力和與固體燃料棒壁面的摩擦力的綜合影響來(lái)建立運(yùn)動(dòng)方程。

z′=-u

(1)

(2)

其中：z為熔融單元的位移；u為熔融單元的運(yùn)動(dòng)速度；g為重力加速度；t為計(jì)算時(shí)間；Ff為熔融單元與固體燃料棒壁面的摩擦阻力；dm=ρdV為熔融單元的質(zhì)量，ρ為二氧化鈾的密度，dV為熔融單元的體積。

由流體力學(xué)知識(shí)可知，由流體黏性而產(chǎn)生的摩擦阻力Ff與流體的動(dòng)能因子f及流體與壁面的接觸面積呈正比。通過(guò)Darcy-Werbach公式可推導(dǎo)得到壁面的摩擦因數(shù)λ=4f，λ的計(jì)算可采用柯列勃洛克(Colebrook)公式[6]：

(3)

其中：ε為絕對(duì)粗糙度(壁面凸出部分的平均高度)；d為水力學(xué)直徑。

柯列勃洛克公式是隱式方程，計(jì)算過(guò)程中可使用牛頓法求解。對(duì)于運(yùn)動(dòng)方程的求解，本研究應(yīng)用四階Runge-Kutta迭代方法來(lái)確定每個(gè)單位時(shí)刻內(nèi)熔融物的移動(dòng)速度和相對(duì)位移，從而解出熔融物新時(shí)刻的位置。

在不考慮傳熱的情況下，選取單個(gè)熔融燃料棒單元進(jìn)行計(jì)算，其運(yùn)動(dòng)結(jié)果如圖4所示。由圖4可看出，熔融單元離開(kāi)初始位置后，在重力作用下速度不斷增大，在1.5 s左右到達(dá)堆芯底部。在這一過(guò)程中，由于熔融單元的速度增大，熔融單元與固體燃料棒壁面的摩擦力也隨之增大，導(dǎo)致其下落加速度減小，速度升高趨勢(shì)減緩，并最終趨向平衡狀態(tài)。同時(shí)，Re正比于熔融單元的速度，因此Re變化趨勢(shì)與速度相同。此時(shí)Re=4×103～3×104，滿(mǎn)足柯列勃洛克公式的適用條件。

圖4 熔融單元下落過(guò)程中的速度和雷諾數(shù)Fig.4 Velocity and Re of melting material in dropping process

2.2 熔融物瞬態(tài)熱傳導(dǎo)模型

進(jìn)行熔融物與燃料棒換熱分析時(shí)作如下假設(shè)：1) 相較導(dǎo)熱、對(duì)流換熱兩種熱交換形式，不考慮熱輻射影響；2) 強(qiáng)迫對(duì)流換熱系數(shù)遠(yuǎn)大于導(dǎo)熱系數(shù)，熔融物運(yùn)動(dòng)狀態(tài)下控制方程可直接忽略接觸導(dǎo)熱項(xiàng)；3) 由于流道堵塞，不考慮自然對(duì)流影響，此假設(shè)在之后的結(jié)果分析中也將得到驗(yàn)證。

采用大空間自然對(duì)流模型計(jì)算對(duì)流換熱熱流密度。在時(shí)間Δτ內(nèi)，熔融單元運(yùn)動(dòng)狀態(tài)下與燃料棒發(fā)生熱量交換的控制方程為：

mcΔT=(h(Trod-Tmelt)A+qΦΔV)Δτ

(4)

采用經(jīng)典導(dǎo)熱模型計(jì)算導(dǎo)熱熱流密度。在時(shí)間Δτ內(nèi)，熔融單元靜止?fàn)顟B(tài)下與燃料棒發(fā)生熱量交換的控制方程為：

(5)

其中：m為熔融單元質(zhì)量；c為比熱容；ΔT為熔融單元溫度的變化；h為強(qiáng)迫對(duì)流換熱系數(shù)；Trod為熔融單元所接觸燃料棒表面的對(duì)應(yīng)溫度；Tmelt為熔融單元溫度；A為熔融單元與燃料棒的接觸面積；qΦ為熔融單元單位體積的熱功率；ΔV為熔融單元體積；λ為導(dǎo)熱系數(shù)；δ為熔融單元徑向厚度。

2.3 熔融物運(yùn)動(dòng)和傳熱的耦合模型

在衰變加熱和喪失冷卻作用下，堆芯工況不斷惡化。隨溫度升高，堆芯內(nèi)部具有不同熔點(diǎn)的組件相繼發(fā)生熔融坍塌。在不考慮共晶反應(yīng)情況下，首先是不銹鋼等部件熔融，隨后是鋯合金，最后是燃料芯塊。通過(guò)耦合熔融物質(zhì)的運(yùn)動(dòng)和傳熱過(guò)程，從理論上模擬出熔融堆芯的遷移過(guò)程。

圖5 熔融物和燃料棒的傳熱Fig.5 Heat transfer between melting material and rod

首先當(dāng)堆芯組件熔融。一部分燃料芯塊先達(dá)到熔點(diǎn)并發(fā)生相變，形成黏性熔融氧化物。熔融物和燃料棒的傳熱如圖5所示。由圖5可見(jiàn)：下滑過(guò)程中熔融物發(fā)生了以強(qiáng)迫對(duì)流換熱為主導(dǎo)的傳熱過(guò)程；冷凝后換熱方式主要以接觸導(dǎo)熱為主；另外的輻射傳熱過(guò)程暫不做考慮。在堆芯熔融過(guò)程中，考慮熔融區(qū)上部的芯塊依次掉落疊加，保持原有的豎直狀態(tài)。從俯視角度觀察，堆芯因中心功率較高而形成凹陷的熔池。

其次是熔融堆芯的堵塞。在流道中，下滑的熔融物與低溫棒面持續(xù)換熱，導(dǎo)致溫度重新降回熔點(diǎn)。對(duì)于進(jìn)入固液兩相區(qū)的熔融物，在流道未被堵塞時(shí)，冷卻劑蒸汽通過(guò)自然對(duì)流進(jìn)一步冷卻熔融物；而在冷卻劑蒸干和流道堵塞時(shí)，不考慮自然對(duì)流。根據(jù)Veshchunov等[7]的理論，熔融物外部形成凝固殼，內(nèi)部是液態(tài)熔融物并緩慢下滑。本文中，由于形成的固體阻礙，認(rèn)為熔融物一旦降到熔點(diǎn)便停止運(yùn)動(dòng)，等待進(jìn)一步的衰變加熱。固液比α對(duì)熔融物遷移影響如圖6所示。由圖6可知，對(duì)于處于兩相區(qū)的熔融物，固液比α取值影響總體遷移過(guò)程中的運(yùn)動(dòng)和靜止的交換頻率。

圖6 固液比α對(duì)熔融物遷移的影響Fig.6 Influence of solid-liquid ratio α on melting material migration

對(duì)于“M”型分布的軸向功率因子[8]，各組件的燃料棒按相同數(shù)值分布。因此，各組件均在相似高度發(fā)生熔融，并根據(jù)不同的棒束溫度在不同高度發(fā)生冷凝，棒束間冷凝的熔融物形成碗狀的冷凝殼層，熔融池最終呈包絡(luò)狀態(tài)[9]。鑒于流道被堵塞，后續(xù)產(chǎn)生的熔融物以液態(tài)形式堆積在冷凝殼層上。不考慮熔融物水平方向運(yùn)動(dòng)，可得到如三哩島熔毀堆芯的4層結(jié)構(gòu)：完整下部、冷凝殼層、液態(tài)熔池和碎片床[3]。

最后是冷凝殼的再熔融。通過(guò)衰變釋熱和接觸導(dǎo)熱，?？康睦淠龤雍芸炀瓦_(dá)到熱平衡。接著冷凝殼和所接觸的未熔融芯塊一起衰變同步升溫并重新達(dá)到熔點(diǎn)。待到底層的冷凝殼開(kāi)始液化，它帶動(dòng)整個(gè)熔池向下遷移，并發(fā)生持續(xù)的換熱降溫，直到底層熔融物再次降到熔點(diǎn)冷凝。由圖6可見(jiàn)，對(duì)于3 m多高的堆芯，在300多秒內(nèi)發(fā)生多次熔融-冷凝遷移。因此不斷擴(kuò)大的熔池的遷移速度受傳熱和運(yùn)動(dòng)耦合過(guò)程的影響。

3 結(jié)果與討論

3.1 驗(yàn)證性計(jì)算

本文所計(jì)算的HPR1000堆型在失水事故(LOCA)后堆芯熔融份額隨時(shí)間的變化如圖7所示。由圖7可見(jiàn)：在考慮鋯水反應(yīng)后，在2 500 s附近堆芯發(fā)生熔融時(shí)刻與文獻(xiàn)[10]的相近；壓力容器中的水在2 500 s時(shí)已蒸干，堆芯上部喪失冷卻開(kāi)始熔融，而文獻(xiàn)[10]在4 870 s前冷卻劑尚未蒸干，堆芯呈緩慢融化狀態(tài)；堆芯下部于3 500 s后也發(fā)生熔融，加劇了熔融份額，6 000 s后堆芯熔融份額相對(duì)穩(wěn)定，其熔融進(jìn)程與文獻(xiàn)[10]后期較為吻合，且熔融份額總體相近。

圖7 堆芯熔融份額隨時(shí)間的變化Fig.7 Change of core melting proportion with time

3.2 堆芯溫度分布

LOCA發(fā)生后，堆芯余熱使得堆芯內(nèi)一部分芯塊在2 400 s左右達(dá)到熔點(diǎn)。取1/4堆芯作為研究對(duì)象，等距地取8個(gè)縱截面研究熔融物的軸向遷移，其初始溫度如圖8所示。圖8中縱截面高溫區(qū)在3.25 m處。根據(jù)溫度分布對(duì)橫截面取特征明顯的0.6、1.2、1.8、2.4、3和3.25 m處截面，其橫截面在2 400 s時(shí)的溫度如圖9所示。

3.3 堆芯熔融遷移過(guò)程

1) 堆芯熔融物縱截面分布

先從小尺度上觀察相鄰3根燃料棒的熔融物遷移，如圖10所示。整個(gè)遷移過(guò)程可分為3個(gè)階段：熔融物生成累積、燃料棒上下部全面熔融和熔融物的穩(wěn)定發(fā)展。由圖10可見(jiàn)，隨時(shí)間發(fā)展，燃料棒中上部發(fā)生熔融，熔融物在中部發(fā)生累積，上部的芯塊下落形成下凹，隨后下部的芯塊開(kāi)始熔融，加快棒束的熔融進(jìn)程，3根棒束熔融物進(jìn)入充分發(fā)展階段。

圖8 2 400 s時(shí)8個(gè)特征縱截面的初始溫度Fig.8 Initial temperature of 8 vertical sections at 2 400 s

圖9 2 400 s時(shí)不同橫截面處的初始溫度Fig.9 Initial temperature of different cross sections at 2 400 s

t1～t6為特征時(shí)間圖10 3根相鄰燃料棒的熔融進(jìn)程Fig.10 Melting process of 3 neighbor rods

從大尺度上觀察1/4堆芯的縱截面，第28排1/4堆芯熔融進(jìn)程如圖11所示，從第2 400 s開(kāi)始取8個(gè)代表時(shí)間點(diǎn)觀察熔融物遷移過(guò)程。在前3個(gè)時(shí)刻，堆芯熔融物生成并累積形成冷凝殼層；到中間3個(gè)時(shí)刻，隨堆芯下部的進(jìn)一步熔融，冷凝殼從中部遷移到下部與新生成冷凝殼匯聚，并在 4 900 s后形成1.5 m高的穩(wěn)定熔池；最后兩時(shí)刻經(jīng)過(guò)2 400 s，熔池充分發(fā)展，不久便會(huì)擊穿堆芯底部，且逃逸出堆芯。

對(duì)于堆芯而言，熔池形態(tài)和組件功率分布緊密相關(guān)，第45排組件功率分布與熔池溫度對(duì)比如圖12所示。由圖12可見(jiàn)，對(duì)于第45排縱截面的熔池發(fā)展中期形貌，各組件的熔融進(jìn)程和功率大小吻合。另外堆芯外圍組件由于低功率因子以至難以發(fā)生熔融現(xiàn)象，形成了如屏蔽體一樣的堆芯結(jié)構(gòu)，將高溫熔池盛在內(nèi)部，一定程度上緩解了熔融物的外泄。

圖11 第28排1/4堆芯熔融進(jìn)程Fig.11 Melting process of row 28 in 1/4 core

圖12 第45排組件熔池的溫度(a)與功率分布(b)Fig.12 Molten pool temperature (a) and assembly power distribution (b) of row 45 assembly

圖13 4個(gè)特征縱截面熔融比Fig.13 Molten ratio of 4 vertical sections

4個(gè)代表縱截面的熔融比如圖13所示。在3 500 s附近，下部堆芯熔融從而加大熔融比曲線(xiàn)斜率；在4 900 s附近，堆芯內(nèi)部大部分芯塊熔融，熔池趨于穩(wěn)定，后續(xù)外圍組件緩慢發(fā)生熔融。若考慮熔融物橫向遷移和輻射換熱，堆芯外圍組件會(huì)加快熔融進(jìn)程，破壞熔融物在堆芯的滯留狀態(tài)，使其遷移到下支撐板和下封頭。

2) 堆芯熔融物橫截面分布

從大尺度觀察1/4堆芯的橫截面，取3個(gè)特征時(shí)刻下不同高度橫截面，觀察水平方向上熔融物的分布及熔池的形成發(fā)展進(jìn)程。

圖9中的橫截面初始溫度分布與1/4堆芯組件徑向功率分布[11]決定了各組件發(fā)生熔融現(xiàn)象的先后時(shí)間。4個(gè)特征橫截面的熔融物分布如圖14所示。熔融前期，來(lái)自率先熔融組件的熔融物滑落至堆芯中部，并成片狀分布，造成流道堵塞，驗(yàn)證了熱力學(xué)控制方程中不具備自然對(duì)流條件的假設(shè)。熔融中期，碗狀的冷凝殼層已完全成形，且各組件熔融進(jìn)程與1/4堆芯組件徑向功率分布[11]相吻合。熔融后期，堆芯內(nèi)大部分組件熔融且在堆芯底部形成穩(wěn)定熔池，其上方形成巨大空穴，同時(shí)外圍組件仍未發(fā)生大規(guī)模熔融，再次驗(yàn)證了其能緩解熔融物橫向擴(kuò)散的判斷。

在同一水平面上，熔融物面積占比如圖15所示，其隨高度的降低而減小，可見(jiàn)熔池并非均勻遷移。一方面由于部分熔融物凝固附著在燃料棒上脫節(jié)，另一方面，伴隨著冷凝殼層熔融破裂，其所承載的大量熔融物開(kāi)始在約3 500 s后直接掉落至堆芯底部。

在冷凝殼層所承載的熔池逐漸在堆芯底部沉積的同時(shí)，熔池上方的空穴慢慢深入擴(kuò)大。燃料棒熔融下沉導(dǎo)致的空位數(shù)與1/4堆芯燃料棒總數(shù)之比示于圖16，在3 500 s左右空穴深入到堆芯中下部，隨后由于空穴擴(kuò)展至難以熔融下沉的外圍組件，各水平高度空穴面積變化變緩且逐漸趨同，空穴逐漸由錐體變成貫穿堆芯的柱體。

第1行為熔融前期，2 500 s時(shí)刻；第2行為熔融中期，3 600 s時(shí)刻；第3行為熔融后期，4 900 s時(shí)刻圖14 4個(gè)特征橫截面的熔融物分布Fig.14 Melting material distribution of 4 cross sections

圖15 4個(gè)特征橫截面的熔融物面積占比Fig.15 Melting material area ratio of 4 cross sections

圖16 5個(gè)特征橫截面的空位比Fig.16 Vacancy ratio of 5 cross sections

4 結(jié)論

為研究嚴(yán)重工況下壓水堆堆芯熔融坍塌問(wèn)題，本文建立了精確的三維堆芯運(yùn)動(dòng)模型求解堆芯熔融物的瞬時(shí)位置，采用大空間自然換熱模型和經(jīng)典導(dǎo)熱模型求解堆芯瞬時(shí)溫度，使用了HPR1000堆型的真實(shí)參數(shù)以模擬堆芯升溫及熔融進(jìn)程。研究結(jié)果表明：停堆后約2 400 s開(kāi)始出現(xiàn)熔融現(xiàn)象，熔融物在堆芯活性區(qū)域內(nèi)下落且發(fā)生多重相變過(guò)程；在4 900 s后，熔融物在堆芯底部形成約1.5 m高的穩(wěn)定熔池；由于外圍組件與低溫圍欄裝置換熱，最外圍的組件不會(huì)發(fā)生熔融。

模型改進(jìn)方面，若計(jì)算中考慮熱輻射項(xiàng)，預(yù)計(jì)將加速熔池的橫向擴(kuò)展。另外，希望能找到計(jì)算強(qiáng)迫對(duì)流換熱系數(shù)的有效公式，從常量?jī)?yōu)化為與速度、溫度等參數(shù)相關(guān)的變量，預(yù)計(jì)熔融物運(yùn)動(dòng)和靜止的交換頻率將會(huì)加快。最后，燃料芯塊所采取的物性參數(shù)均為純二氧化鈾，若考慮實(shí)際燃料芯塊，共晶熔融物的物性問(wèn)題有待解決。