壓水堆堆芯熔融物的運動及傳熱數值分析

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

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

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

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

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

1 堆芯幾何模型

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

圖1 1/4堆芯的組件結構Fig.1 1/4 core cross-sectional shape

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

表1 堆芯材料物性參數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

獲取嚴重事故下堆芯的初始溫度和功率分布后，使用時間推進方法構建傳熱與運動的耦合模型。取瞬態時間步長為1 s，且將每個時間步長細分為6個子時間步。每經過1個子時間步，程序將計算各燃料棒單元的溫度，并對每個運動的熔融單元求解一次運動和傳熱方程，以確定每個時刻熔融燃料棒單元的位置與溫度信息，同時每1 s輸出1次計算結果。

2 數值模型

2.1 熔融物瞬態運動模型

當燃料棒單元的溫度達到熔點時，熔融的燃料棒單元將離開原來的位置，在對應的流道中下落。在下落過程中，為計算高溫熔融物單元的位置，考慮熔融物的重力和與固體燃料棒壁面的摩擦力的綜合影響來建立運動方程。

z′=-u

(1)

(2)

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

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

(3)

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

柯列勃洛克公式是隱式方程，計算過程中可使用牛頓法求解。對于運動方程的求解，本研究應用四階Runge-Kutta迭代方法來確定每個單位時刻內熔融物的移動速度和相對位移，從而解出熔融物新時刻的位置。

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

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

2.2 熔融物瞬態熱傳導模型

進行熔融物與燃料棒換熱分析時作如下假設：1) 相較導熱、對流換熱兩種熱交換形式，不考慮熱輻射影響；2) 強迫對流換熱系數遠大于導熱系數，熔融物運動狀態下控制方程可直接忽略接觸導熱項；3) 由于流道堵塞，不考慮自然對流影響，此假設在之后的結果分析中也將得到驗證。

采用大空間自然對流模型計算對流換熱熱流密度。在時間Δτ內，熔融單元運動狀態下與燃料棒發生熱量交換的控制方程為：

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

(4)

采用經典導熱模型計算導熱熱流密度。在時間Δτ內，熔融單元靜止狀態下與燃料棒發生熱量交換的控制方程為：

(5)

其中：m為熔融單元質量；c為比熱容；ΔT為熔融單元溫度的變化；h為強迫對流換熱系數；Trod為熔融單元所接觸燃料棒表面的對應溫度；Tmelt為熔融單元溫度；A為熔融單元與燃料棒的接觸面積；qΦ為熔融單元單位體積的熱功率；ΔV為熔融單元體積；λ為導熱系數；δ為熔融單元徑向厚度。

2.3 熔融物運動和傳熱的耦合模型

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

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

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

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

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

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

最后是冷凝殼的再熔融。通過衰變釋熱和接觸導熱，停靠的冷凝殼層很快就達到熱平衡。接著冷凝殼和所接觸的未熔融芯塊一起衰變同步升溫并重新達到熔點。待到底層的冷凝殼開始液化，它帶動整個熔池向下遷移，并發生持續的換熱降溫，直到底層熔融物再次降到熔點冷凝。由圖6可見，對于3 m多高的堆芯，在300多秒內發生多次熔融-冷凝遷移。因此不斷擴大的熔池的遷移速度受傳熱和運動耦合過程的影響。

3 結果與討論

3.1 驗證性計算

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

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

3.2 堆芯溫度分布

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

3.3 堆芯熔融遷移過程

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

先從小尺度上觀察相鄰3根燃料棒的熔融物遷移，如圖10所示。整個遷移過程可分為3個階段：熔融物生成累積、燃料棒上下部全面熔融和熔融物的穩定發展。由圖10可見，隨時間發展，燃料棒中上部發生熔融，熔融物在中部發生累積，上部的芯塊下落形成下凹，隨后下部的芯塊開始熔融，加快棒束的熔融進程，3根棒束熔融物進入充分發展階段。

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

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

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

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

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

圖11 第28排1/4堆芯熔融進程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個特征縱截面熔融比Fig.13 Molten ratio of 4 vertical sections

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

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

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

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

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

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

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

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

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

4 結論

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

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