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

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

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

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

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

1 堆芯幾何模型

根據(jù)HPR1000堆型堆芯功率和溫度分布的對(duì)稱性，選取如圖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

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

表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í)間步長為1 s，且將每個(gè)時(shí)間步長細(xì)分為6個(gè)子時(shí)間步。每經(jīng)過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í)，熔融的燃料棒單元將離開原來的位置，在對(duì)應(yīng)的流道中下落。在下落過程中，為計(jì)算高溫熔融物單元的位置，考慮熔融物的重力和與固體燃料棒壁面的摩擦力的綜合影響來建立運(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及流體與壁面的接觸面積呈正比。通過Darcy-Werbach公式可推導(dǎo)得到壁面的摩擦因數(shù)λ=4f，λ的計(jì)算可采用柯列勃洛克(Colebrook)公式[6]：

(3)

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

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

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

圖4 熔融單元下落過程中的速度和雷諾數(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)情況下，首先是不銹鋼等部件熔融，隨后是鋯合金，最后是燃料芯塊。通過耦合熔融物質(zhì)的運(yùn)動(dòng)和傳熱過程，從理論上模擬出熔融堆芯的遷移過程。

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

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

其次是熔融堆芯的堵塞。在流道中，下滑的熔融物與低溫棒面持續(xù)換熱，導(dǎo)致溫度重新降回熔點(diǎn)。對(duì)于進(jìn)入固液兩相區(qū)的熔融物，在流道未被堵塞時(shí)，冷卻劑蒸汽通過自然對(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ū)的熔融物，固液比α取值影響總體遷移過程中的運(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]。

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

3 結(jié)果與討論

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

本文所計(jì)算的HPR1000堆型在失水事故(LOCA)后堆芯熔融份額隨時(shí)間的變化如圖7所示。由圖7可見：在考慮鋯水反應(yīng)后，在2 500 s附近堆芯發(fā)生熔融時(shí)刻與文獻(xiàn)[10]的相近；壓力容器中的水在2 500 s時(shí)已蒸干，堆芯上部喪失冷卻開始熔融，而文獻(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 堆芯熔融遷移過程

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

先從小尺度上觀察相鄰3根燃料棒的熔融物遷移，如圖10所示。整個(gè)遷移過程可分為3個(gè)階段：熔融物生成累積、燃料棒上下部全面熔融和熔融物的穩(wěn)定發(fā)展。由圖10可見，隨時(shí)間發(fā)展，燃料棒中上部發(fā)生熔融，熔融物在中部發(fā)生累積，上部的芯塊下落形成下凹，隨后下部的芯塊開始熔融，加快棒束的熔融進(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開始取8個(gè)代表時(shí)間點(diǎn)觀察熔融物遷移過程。在前3個(gè)時(shí)刻，堆芯熔融物生成并累積形成冷凝殼層；到中間3個(gè)時(shí)刻，隨堆芯下部的進(jìn)一步熔融，冷凝殼從中部遷移到下部與新生成冷凝殼匯聚，并在 4 900 s后形成1.5 m高的穩(wěn)定熔池；最后兩時(shí)刻經(jīng)過2 400 s，熔池充分發(fā)展，不久便會(huì)擊穿堆芯底部，且逃逸出堆芯。

對(duì)于堆芯而言，熔池形態(tài)和組件功率分布緊密相關(guān)，第45排組件功率分布與熔池溫度對(duì)比如圖12所示。由圖12可見，對(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附近，下部堆芯熔融從而加大熔融比曲線斜率；在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所示。熔融前期，來自率先熔融組件的熔融物滑落至堆芯中部，并成片狀分布，造成流道堵塞，驗(yàn)證了熱力學(xué)控制方程中不具備自然對(duì)流條件的假設(shè)。熔融中期，碗狀的冷凝殼層已完全成形，且各組件熔融進(jìn)程與1/4堆芯組件徑向功率分布[11]相吻合。熔融后期，堆芯內(nèi)大部分組件熔融且在堆芯底部形成穩(wěn)定熔池，其上方形成巨大空穴，同時(shí)外圍組件仍未發(fā)生大規(guī)模熔融，再次驗(yàn)證了其能緩解熔融物橫向擴(kuò)散的判斷。

在同一水平面上，熔融物面積占比如圖15所示，其隨高度的降低而減小，可見熔池并非均勻遷移。一方面由于部分熔融物凝固附著在燃料棒上脫節(jié)，另一方面，伴隨著冷凝殼層熔融破裂，其所承載的大量熔融物開始在約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)重工況下壓水堆堆芯熔融坍塌問題，本文建立了精確的三維堆芯運(yùn)動(dòng)模型求解堆芯熔融物的瞬時(shí)位置，采用大空間自然換熱模型和經(jīng)典導(dǎo)熱模型求解堆芯瞬時(shí)溫度，使用了HPR1000堆型的真實(shí)參數(shù)以模擬堆芯升溫及熔融進(jìn)程。研究結(jié)果表明：停堆后約2 400 s開始出現(xiàn)熔融現(xiàn)象，熔融物在堆芯活性區(qū)域內(nèi)下落且發(fā)生多重相變過程；在4 900 s后，熔融物在堆芯底部形成約1.5 m高的穩(wěn)定熔池；由于外圍組件與低溫圍欄裝置換熱，最外圍的組件不會(huì)發(fā)生熔融。

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