Development of a subchannel code for blockage accidents of LMFRs based on the 3D fuel rod model

Development of a subchannel code for blockage accidents of LMFRs based on the 3D fuel rod model 1 first-author 1 1 1 corresp hlchen1@ustc.edu.cn hlchen1@ustc.edu.cn hlchen1@ustc.edu.cn https://orcid.org/0000-0001-6855-1069 School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230026, China School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230026, China To predict the thermal-hydraulic (T/H) parameters of the reactor core for liquid-metal-cooled fast reactors (LMFRs), especially under flow blockage accidents, we developed a subchannel code called KMC-FB. This code uses a time-dependent, four-equation, single-phase flow model together with a 3D heat conduction model for the fuel rods, which is solved by numerical methods based on the finite difference method with a staggered mesh. Owing to the local effect of the blockage on the flow field, low axial flow, increased forced cross flow, and backflow occur. To more accurately simulate this problem, we implemented a robust and novel solution method. We verified the code with a low-flow (~0.01 m/s) and large-scale blockage case. For the preliminary validation, we compared our results with the experimental data of the NACIE-UP BFPS blockage test and the KIT 19ROD blockage test. The results revealed that KMC-FB has sufficient solution accuracy and can be used in future flow blockage analyses for LMFRs. Subchannel method Code development Blockage accident Liquid metal cooled fast reactor Introduction 1 Fourth-generation liquid-metal-cooled fast reactors (LMFRs), including the lead-cooled fast reactor (LFR) and sodium-cooled fast reactor (SFR), have application prospects in the future nuclear power sector because of their high level of inherent safety, advantages in new fissile fuel breeding, and minimization of nuclear waste. In particular, LMFRs are expected to provide a safe and self-protected response to the Fukushima scenario because of the above-mentioned features providing high thermal inertia. Although LMFRs have the above advantages, some safety issues have been noted by researchers during the design and operation of the above-mentioned prototypes or experimental reactors, such as the sodium boiling problem [ 1 1 ], flow blockage accidents [ 2 2 ], erosion problems [ 3 3 ], the production of polonium-210 under neutron irradiation of Bi [ 4 4 ], the temperature control problem, and the in-reactor detection of opaque liquid metal [ 5 5 ]. Therefore, in recent years, efforts have been directed toward improve the safety of LMFRs. One method that can avoid the high cost of experiments and quickly predict the operating parameters under the design parameters of the reactor is the development and application of various numerical simulation programs. Flow blockage is one of these postulated accidents and events affecting LMFRs; these blockages can form for the following reasons [ 6 6 ]: (1) accumulation of corrosion products, (2) swelling and curving of the cladding, (3) shedding of the wire-wrap spacer, and (4) accumulation of other impurities. The detailed accident process diagram is shown in Figure 1 Figure 1 , which we divided into three parts: (1) the formation of the blockage accident, (2) the dynamic development of the core parameters after the blockage, and (3) the accident process after the pellets melted. At present, research on reactor blockage accidents can be roughly categorized into two aspects: experimental research and numerical simulation. In terms of experimental research, the United States, Germany, Japan, Italy, and other countries have successively conducted a series of experimental studies on the blockage accidents of single-fuel assemblies since the 1970s, considering liquid sodium, air, water, and lead–bismuth alloys (LBE) as coolants [ 6 6 ]-[ 11 11 ]. For example, KIT conducted a heat transfer experiment with LBE flowing through a partially blocked 19-rod bundle fixed by a wire spacer [ 10 10 ]. ENEA constructed an experimental platform called a blocked fuel pin bundle simulator (BFPS) installed in the NACIE-UP facility to study the thermal-hydraulic features under different blockage situations [ 11 11 ]. However, owing to factors such as the high cost and requirements of the experimental facility, conducting such blockage experiments is difficult. Most of the research has focused on numerical simulations and the study of related theoretical models for blockage accidents, which can be divided into five main categories: · Computational fluid dynamics (CFD) can be used to carry out three-dimensional (3D) thermal-hydraulic analysis of blockage accidents [ 12 12 , 13 13 ]. With a refined model and an appropriate turbulence model, the local flow field can be accurately obtained, but the simulation cost is much higher than that of the following tools. · One-dimensional system analysis codes can be used to obtain the thermal hydraulic parameters under flow blockage conditions. Examples of these methods include RELAP and ATHLET [ 14 14 - 16 16 14-16 ]. These codes can quickly calculate the core outlet temperature in the case of a flow blockage, but the local parameters inside the fuel assembly cannot be accurately obtained owing to the application of the lumped parameter method. · A 1.5-dimensional subchannel code can be used to predict the maximum coolant temperature, maximum cladding temperature, and pellet temperature in the fuel assembly under steady-state or transient flow blockage accidents. Examples include COBRA-IV, SABRE4, and MATRA-LMR [ 17 17 - 19 19 17-19 ]. However, most of the existing subchannel codes for LMFRs adopt a relatively simple blockage model, and the details still remain to be resolved. · Serious accident analysis codes can be used for accidents that occur after a large-scale blockage and may result in core melting, such as SIMMER, ATHLET-CD, and RELAP5/SCDAP [ 20 20 - 22 22 20-22 ]. These codes have the advantages of performing dynamic system analysis under severe accidents of the reactor core but perform poorly for fine local analysis. · Multi-physics coupling programs can be used to analyze the comprehensive impact on the neutron physics, T/H, and structural safety [ 23 23 - 25 25 23-25 ]. However, the basis of coupling is the mature development of each subtool. The subchannel thermal hydraulic analysis method has been used in nuclear engineering for many years. In this method, the fuel assembly is artificially divided into several subchannels, and the conservation equation of the flow field and the heat conduction equation of the fuel rod are solved simultaneously to obtain the distribution of temperature, flow rate, and pressure in the assembly. As such, in this study, we applied the subchannel method to the reactor core flow blockage accident. We mainly focused on the second part of Fig. 1 Fig. 1 : calculating the consequences of the accident without considering the fictitious accident caused by the blockage and melting. In our previous work, we developed a subchannel code KMC-SUBtra for transient thermal hydraulic analysis of LMFRs [ 26 26 ]. In this study, our objective was to expand the capability of KMC-SUBtra in dealing with this type of accident. Based on this goal, we constructed a new numerical solution method, adopted a more refined three-dimensional heat conduction model, and performed an in-depth study of blockage-related models (i.e., axial resistance model, heat transfer model, etc.), which laid the foundation for our upgraded code (referred to as KMC-FB). In addition, to ensure the correctness and accuracy of KMF-FB, a low-flow (~0.01 m/s) and large-scale blockage case was used to test the code, and we preliminary validated the code by comparing its results with the experimental data of the NACIE-UP BFPS blockage test and the KIT 19ROD blockage test. Basic conservation equations and physical models 1 Basic conservation equations 2 The subchannel code KMC-FB adopts the finite difference method (FDM) with a staggered mesh, which adopts a scalar field grid for the conservation continuity equation and energy equation, and uses a vector grid for the momentum equation (axial and lateral). The height of the control cell is half-staggered between the two grids, as shown in Fig. 2 Fig. 2 . Owing to the assumption that the liquid metal does not boil, a four-equation, single-phase flow model including the continuity equation, axial and lateral momentum equation, and energy conservation equation is used, which is expressed in the form of differential equations, as shown in Eqs. (1)–(4). The meanings of the variables used in the following equations are provided in the nomenclature. Continuity equation: <math display="block" id="M1"><mrow><mi>A</mi><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>w</mi><mi>A</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>n</mtext><mtext>u</mtext><mtext>m</mtext></mrow></munder><mrow><mfrac><mrow><mo>∂</mo><msub><mrow><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>v</mi><msub><mi>A</mi><mtext>T</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mrow></mstyle><mo>=</mo><mtext>0</mtext></mrow></math> Axial momentum equation: <math display="block" id="M2"><mtable columnalign="left"><mtr><mtd><mi>A</mi><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>w</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>A</mi><mi>ρ</mi><msup><mi>w</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>w</mi><mi>v</mi><msub><mi>A</mi><mtext>T</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mrow></mstyle></mtd></mtr><mtr><mtd><mo>=</mo><mo>−</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>P</mi><mi>A</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><msub><mi>f</mi><mtext>T</mtext></msub><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><msup><mi>V</mi><mo>'</mo></msup></mrow></mstyle><mi>Δ</mi><mi>w</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mfrac><mi>f</mi><mrow><msub><mi>D</mi><mtext>e</mtext></msub></mrow></mfrac><mo>+</mo><mfrac><mi>K</mi><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mi>ρ</mi><mi>w</mi><mo>|</mo><mi>w</mi><mo>|</mo><mi>A</mi></mtd></mtr><mtr><mtd><mo>−</mo><mi>A</mi><mi>ρ</mi><mi>g</mi><mtext>cos</mtext><mi>θ</mi></mtd></mtr></mtable></math> Lateral momentum equation: <math display="block" id="M3"><mtable><mtr><mtd><msub><mi>A</mi><mtext>T</mtext></msub><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>v</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>w</mi><mi>V</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mi>S</mi><mrow><msub><mi>l</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>P</mi><mi>i</mi></msub><mo>−</mo><msub><mi>P</mi><mi>k</mi></msub></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mrow><mn>2</mn><msub><mi>l</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mfrac><msub><mi>K</mi><mtext>G</mtext></msub><mfrac><mrow><mrow><mo>|</mo><mi>v</mi><mo>|</mo></mrow><mi>v</mi></mrow><mrow><mi>ρ</mi><mi>S</mi></mrow></mfrac></mtd></mtr><mtr><mtd><mo>−</mo><mi>S</mi><mi>ρ</mi><mi>g</mi><mtext>sin</mtext><mi>θ</mi></mtd></mtr></mtable></math> Energy conservation equation: <math display="block" id="M4"><mtable><mtr><mtd><mi>A</mi><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>h</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>m</mi><mi>h</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><mi>V</mi><mi>h</mi></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><mi>λ</mi><mfrac><mi>S</mi><mrow><msub><mi>l</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>T</mi><mi>i</mi></msub><mo>-</mo><msub><mi>T</mi><mi>k</mi></msub></mrow><mo>)</mo></mrow></mrow></mstyle></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>n</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><msub><mi>P</mi><mi>w</mi></msub></mrow></mstyle><mrow><mo>〈</mo><mrow><msup><mi>q</mi><mrow><mo>'</mo><mo>'</mo></mrow></msup></mrow><mo>〉</mo></mrow></mtd></mtr><mtr><mtd><mo>−</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>k</mi><mo>∈</mo><mtext>num</mtext></mrow></munder><mrow><msub><mrow><mo stretchy="false">(</mo><mi>Δ</mi><mi>h</mi><msup><mi>V</mi><mo>'</mo></msup><mo stretchy="false">)</mo></mrow><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mstyle><mo>.</mo></mtd></mtr></mtable></math> Heat conduction equation 2 Currently, most subchannel codes, such as COBRA and MATRA, adopt a one-dimensional model to simulate the heat conduction of fuel rods. Under normal conditions, the temperature difference between each subchannel and the circumferential unevenness of the fuel rod are small, which indicates that this method is acceptable. However, for some reactor types whose coolant physical properties vary widely with temperature (i.e., supercritical water-cooled reactors [SCWRs]) or where serious asymmetry of the fuel rod circumference is caused by an accident (i.e., flow blockage accident), the limitations of this method are obvious. In this study, we used a three-dimensional heat conduction model tightly coupled with fluid conservation equations. For comparison, we also implanted a one-dimensional heat conduction equation and a circumferential subnode model in KMC-FB. The one-dimensional heat conduction equation usually does not consider the heat conduction of the circumferential nodes, and the axial temperature depends on the boundary conditions of the convective fluid, as shown in Fig. 3(a) Fig. 3(a) . The circumferential subnode model simply considers that the fuel rods are split according to the surrounding subchannels, as shown in Fig.3(b) Fig.3(b) . The three-dimensional heat conduction model considers the heat conduction of the fuel rod in the axial, radial, and circumferential directions, as shown in Eq. (5) and Fig. 3(c) Fig. 3(c) : <math display="block" id="M5"><mrow><mi>ρ</mi><msub><mi>c</mi><mtext>p</mtext></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mn>1</mn><mi>r</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mi>λ</mi><mi>r</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><mi>r</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>θ</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mi>λ</mi><mi>r</mi></mfrac><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>θ</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mi>λ</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>q</mi><mtext>v</mtext></msub><mo>.</mo></mrow></math> Hybrid difference scheme 2 Blockage accidents can lead to local low flow rates and backflow phenomena, which are similar to the flow around a cylinder and change the flow from Stokes flow to turbulent flow through a laminar region [ 27 27 ]. If the upwind scheme is used in all flow regions, a larger numerical diffusion will occur in the local low-velocity flow and the circulating flow region after the blockage. To avoid this numerical dissipation, the hybrid difference method [ 28 28 ] is adopted, which combines the advantages of the central difference scheme and the upwind difference scheme. The central difference scheme is used when Re is smaller than 2 or Pe is smaller than 2, and the upwind difference scheme is used in other cases. Specifically used in the subchannel conservation equation, the convection term of the momentum and energy equations can be expressed as Eqs. (6)-(8): <math display="block" id="M6"><mrow><mi>w</mi><mfrac><mrow><mo>∂</mo><mi>m</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mi>w</mi><mo>¯</mo></mover><msup><mi>m</mi><mo>*</mo></msup></mrow><mo>)</mo></mrow></mrow><mi>j</mi></msub><mo>−</mo><msub><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mi>w</mi><mo>¯</mo></mover><msup><mi>m</mi><mo>*</mo></msup></mrow><mo>)</mo></mrow></mrow><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>R</mi><mi>e</mi><mo>≥</mo><mn>2</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><msubsup><mi>m</mi><mi>j</mi><mo>*</mo></msubsup><mo>+</mo><msubsup><mi>m</mi><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow><mo>*</mo></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>w</mi><mo>¯</mo></mover><mi>j</mi></msub><mo>−</mo><msub><mover accent="true"><mi>w</mi><mo>¯</mo></mover><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>R</mi><mi>e</mi><mo><</mo><mn>2</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> <math display="block" id="M7"><mrow><mi>V</mi><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mi>V</mi><mrow><mi>j</mi><mo>,</mo><mi>i</mi><mi>k</mi></mrow></msub><msup><mrow/><mo>*</mo></msup><msub><mi>w</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>i</mi><mi>k</mi></mrow></msub><mo>−</mo><msubsup><mi>V</mi><mrow><mi>j</mi><mo>−</mo><mn>1</mn><mo>,</mo><mi>i</mi><mi>k</mi></mrow><mo>*</mo></msubsup><msub><mi>w</mi><mrow><mi>j</mi><mo>,</mo><mi>i</mi><mi>k</mi></mrow></msub></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>R</mi><mi>e</mi><mo>≥</mo><mn>2</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><msubsup><mi>V</mi><mrow><mi>j</mi><mo>,</mo><mi>i</mi><mi>k</mi></mrow><mo>*</mo></msubsup><mo>+</mo><msubsup><mi>V</mi><mrow><mi>j</mi><mo>−</mo><mn>1</mn><mo>,</mo><mi>i</mi><mi>k</mi></mrow><mo>*</mo></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>w</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn><mo>,</mo><mi>i</mi><mi>k</mi></mrow></msub><mo>−</mo><msub><mi>w</mi><mrow><mi>j</mi><mo>,</mo><mi>i</mi><mi>k</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>R</mi><mi>e</mi><mo><</mo><mn>2</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> <math display="block" id="M8"><mrow><mi>m</mi><mfrac><mrow><mo>∂</mo><mi>h</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mi>m</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><msubsup><mi>h</mi><mi>j</mi><mo>*</mo></msubsup><mo>−</mo><msub><mi>m</mi><mi>j</mi></msub><msubsup><mi>h</mi><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow><mo>*</mo></msubsup></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>P</mi><mi>e</mi><mo>≥</mo><mn>2</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><msubsup><mi>h</mi><mi>j</mi><mo>*</mo></msubsup><mo>+</mo><msubsup><mi>h</mi><mrow><mi>j</mi><mo>−</mo><mn>1</mn></mrow><mo>*</mo></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>m</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>m</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><mtext>if </mtext><mi>P</mi><mi>e</mi><mo><</mo><mn>2</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> Flow and heat transfer model 2 Similar to the previous subchannel code KMC-SUBtra, various heat transfer and friction factor correlations are used in the updated code (KMC-FB), and a detailed presentation was provided by Cao et al. [ 26 26 ]. In addition, the heat conduction in the blockage is important for the prediction of the maximum temperature of the cladding and the temperature of the coolant in the case of a blockage. The total heat transfer coefficient ( k ) between the subchannel and the cladding in the blockage flow area can be expressed as the Karlsruhe blockage equation (KBE) [ 10 10 ] in Eq. (9), where R bl is the blockage thermal conductivity and R f is the convection in the fluid. An empirical coefficient r 1 > 0 indicates that thermal conduction dominates in the solid blockage. <math display="block" id="M9"><mrow><mi>k</mi><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>R</mi><mrow><mtext>bl</mtext></mrow></msub><mo>+</mo><msub><mi>R</mi><mtext>f</mtext></msub></mrow></mfrac><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>r</mi><mn>1</mn></msub><mfrac><mrow><msub><mi>D</mi><mtext>e</mtext></msub></mrow><mrow><msub><mi>λ</mi><mrow><mtext>bl</mtext></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>D</mi><mtext>e</mtext></msub></mrow><mrow><msub><mi>λ</mi><mtext>f</mtext></msub><mi>N</mi><msub><mi>u</mi><mrow><mtext>local</mtext></mrow></msub></mrow></mfrac></mrow></mfrac></mrow></math> Flow resistance model 2 The local flow resistance coefficient is used to describe the magnitude of the pressure drop caused by the blockage. Figure 4 Figure 4 shows a schematic of this case. In this study, we used the empirical relationship proposed by IdelChik [ 29 29 ] to calculate the flow resistance coefficient, which is expressed as Eq. (10). <math display="block" id="M10"><mrow><mi>K</mi><mo>=</mo><mrow><mo>[</mo><mrow><mi>C</mi><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>F</mi><mn>0</mn></msub></mrow><mrow><msub><mi>F</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0.75</mn></mrow></msup><mo>+</mo><mi>τ</mi><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>F</mi><mn>0</mn></msub></mrow><mrow><msub><mi>F</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>1.3</mn><mo>/</mo><mn>5</mn></mrow></msup><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>F</mi><mn>0</mn></msub></mrow><mrow><msub><mi>F</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><mi>f</mi><mfrac><mi>L</mi><mrow><msub><mi>D</mi><mtext>h</mtext></msub></mrow></mfrac></mrow><mo>]</mo></mrow><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>F</mi><mn>1</mn></msub></mrow><mrow><msub><mi>F</mi><mn>0</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>,</mo></mrow></math> where K is the flow resistance coefficient; F 0 and F are the flow area at the blockage and the area of the total flow channel (m 2 ), respectively; f is the friction factor; C is a user-defined coefficient, which is determined by the inclination angle of the blocking block relative to the flow direction of the coolant; L and D h are the length of the blockage region and the hydraulic diameter of the flow channel, respectively; and <math id="M11"><mi>τ</mi></math> is a correction coefficient that characterizes the wall thickness and shape of the flow channel, which can be expressed by Eq. (11): <math display="block" id="M12"><mrow><mi>τ</mi><mo>=</mo><mrow><mo>(</mo><mrow><mn>2.4</mn><mo>−</mo><mfrac><mi>L</mi><mrow><msub><mi>D</mi><mtext>h</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mrow><mo>[</mo><mrow><mn>0.25</mn><mo>+</mo><mfrac><mrow><mn>0.535</mn><msup><mrow><mrow><mo>(</mo><mrow><mi>L</mi><mo>/</mo><msub><mi>D</mi><mtext>h</mtext></msub></mrow><mo>)</mo></mrow></mrow><mn>8</mn></msup></mrow><mrow><mn>0.05</mn><msup><mrow><mrow><mo>(</mo><mrow><mi>L</mi><mo>/</mo><msub><mi>D</mi><mtext>h</mtext></msub></mrow><mo>)</mo></mrow></mrow><mn>7</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow></mrow></msup><mo>.</mo></mrow></math> Property models 2 We obtained the property correlations of sodium [ 30 30 ] and integrated them into KMC-FB. We also obtained the property correlations of Pb and LBE [ 31 31 ]. Both are widely used in the thermal hydraulic analysis of LMFRs. Numerical solution method 1 Solution method for flow field 2 KMC-FB used an implicit iterative method as the numerical solution method, which refers to the ideas of the SIMPLE and Newton–Raphson iterative methods but is different enough that its process needs to be described: 1) Initialize a reasonable pressure field, flow velocity field, and temperature field according to the user-input parameters, which favors the convergence speed. 2) Solve the axial and lateral momentum equations and obtain the tentative flow field that satisfies the momentum equation. 3) Calculate the residual error of the continuity equation considering the influence of the upper, lower, and laterally adjacent nodes for a given node; construct the pressure Poisson equation; and obtain the pressure correction in the specified fuel assembly at one time. 4) Update the pressure and flow rate based on the momentum equations. 5) Solve the flow field energy equation and the fuel rod heat conduction equation and update the thermophysical properties. 6) Verify whether the thermal hydraulic parameters converge or exceed the specified number of iteration steps. If yes, exist the iteration process and output the result; if not, return to step (1) to continue the iterative calculation. Solution method for solid 2 To solve the three-dimensional heat conduction equation of a solid, the numerical solution method used in KMC-FB is Brain-alternating direction implicit (ADI), which is unconditionally stable when solving three-dimensional problems [ 32 32 ]. Equation (13) provides a brief introduction to the Brain-ADI process, where <math id="M13"><mrow><msubsup><mi>δ</mi><mi>x</mi><mn>2</mn></msubsup></mrow></math> represents the concise expression of <math id="M14"><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>T</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac></mrow></math> ; U and V represent the intermediate variables obtained during the x-direction and y-direction sweeping, respectively; T is the temperature; and the subscripts n and n +1 represent the previous time-step value and current value, respectively. In general, the ADI method transforms a three-dimensional problem into three implicit one-dimensional problems. For each one-dimensional problem, the tridiagonal matrix algorithm (TDMA, a simplified form of Gaussian elimination) is used to solve it directly: <math display="block" id="M15"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mi>U</mi><mo>−</mo><msub><mi>T</mi><mi>n</mi></msub></mrow><mrow><mi>Δ</mi><mi>t</mi><mo>/</mo><mn>2</mn></mrow></mfrac><mo>=</mo><msubsup><mi>δ</mi><mi>x</mi><mn>2</mn></msubsup><mi>U</mi><mo>+</mo><msubsup><mi>δ</mi><mi>y</mi><mn>2</mn></msubsup><msub><mi>T</mi><mi>n</mi></msub><mo>+</mo><msubsup><mi>δ</mi><mi>z</mi><mn>2</mn></msubsup><msub><mi>T</mi><mi>n</mi></msub></mtd></mtr><mtr><mtd><mfrac><mrow><mi>V</mi><mo>−</mo><msub><mi>T</mi><mi>n</mi></msub></mrow><mrow><mi>Δ</mi><mi>t</mi><mo>/</mo><mn>2</mn></mrow></mfrac><mo>=</mo><msubsup><mi>δ</mi><mi>x</mi><mn>2</mn></msubsup><mi>U</mi><mo>+</mo><msubsup><mi>δ</mi><mi>y</mi><mn>2</mn></msubsup><mi>V</mi><mo>+</mo><msubsup><mi>δ</mi><mi>z</mi><mn>2</mn></msubsup><msub><mi>T</mi><mi>n</mi></msub></mtd></mtr><mtr><mtd><mfrac><mrow><msub><mi>T</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mi>V</mi></mrow><mrow><mi>Δ</mi><mi>t</mi><mo>/</mo><mn>2</mn></mrow></mfrac><mo>=</mo><msubsup><mi>δ</mi><mi>x</mi><mn>2</mn></msubsup><mi>U</mi><mo>+</mo><msubsup><mi>δ</mi><mi>y</mi><mn>2</mn></msubsup><mi>V</mi><mo>+</mo><msubsup><mi>δ</mi><mi>z</mi><mn>2</mn></msubsup><msub><mi>T</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msub></mtd></mtr></mtable></math> Code verification and validation 1 To ensure the correctness and applicability of the above theoretical model and to test the robustness of KMC-FB, we performed preliminary validation and verification (V & V). The verification process was used to ensure that the program functions are consistent with the expectations. Hence, we chose two cases to test the code’s ability to handle low flow rates and circulating flows. In addition, we constructed a validation matrix based on existing liquid metal bundle blockage experiments, the NACIE-UP and KIT THEADS blockage experiments, as shown in Table 1 Table 1 . 2 1 Experiment 1 1 NACIE-UP 1 1   1 1 THEADS 1 1   1 1   2 1 Case 1 1 BFPS-4-0 1 1 BFPS-4-1 1 1 Block-C1 1 1 Block-E1 1 1 Block-center-C6 2 1 Institution/country/year 1 1 ENEA/Italy/2018 1 1   1 1 KIT-KALLA/Germany/2017 2 1 Bundle parameters 1 1   1 1 Rod arrangement and duct type 1 1 Triangular/Hexagon 1 1   1 1 Triangular/Hexagon 1 1   1 1 Number of rods 1 1 19 1 1   1 1 19 1 1   1 1 Rod diameter 1 1 10 1 1   1 1 8.2 1 1   1 1 Pitch 1 1 14 1 1   1 1 10.49 1 1   1 1 Wire diameter (mm) 1 1 /grid/ 1 1   1 1 2.2 1 1   1 1 Wire pitch (mm) 1 1 /grid/ 1 1   1 1 328 1 1   1 1 Width of duct (mm) 1 1 64 1 1   1 1 45.931 1 1   1 1 Inlet zone length (mm) 1 1 100 1 1   1 1 955 1 1   1 1 Heated zone length (mm) 1 1 600 1 1   1 1 870 1 1   1 1 Outlet zone length (mm) 1 1 500 1 1   1 1 404 2 1 Experimental conditions 1 1   1 1 Inlet temperature (℃) 1 1 222.4 1 1 222.4 1 1 200.0 1 1 200.0 1 1 200.0 1 1   1 1 Thermal power (kW) 1 1 22.8 1 1 22.8 1 1 394 1 1 394 1 1 98.5 1 1   1 1 Operation pressure (MPa) 1 1 0.101 1 1 0.101 1 1 0.101 1 1 0.101 1 1 0.101 1 1   1 1 Mass flow rate (kg/s) 1 1 4.05 1 1 4.05 1 1 18.7 1 1 18.7 1 1 18.75 2 1 Blockage characteristics 1 1   1 1 Blockage type 1 1 planar 1 1 planar 1 1 lumpy 1 1 lumpy 1 1 lumpy 1 1   1 1 Length (mm) 1 1 - 1 1 - 1 1 54.6 1 1 54.6 1 1 54.6 1 1   1 1 Size (subchannels) 1 1 Sector 1 1 Sector 1 1 one central subchannel 1 1 one edge subchannel 1 1 Six central subchannels 1 1   1 1 Location 1 1 Inlet of heated zone 1 1 Inlet of heated zone 1 1 492-546.6 1 1 710.7-765.3 1 1 492-546.6 1 1   1 1 Porosity 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1   1 1 Material 1 1 Steel 1 1   1 1 Steel + ceramic 1 1   1 1   Code verification 2 Very low flow (~0.01 m/s) 3 In an unprotected loss of flow accident (ULOF), the liquid metal in the reactor is driven only by natural circulation, which is always followed by a low flow phenomenon. Owing to the low flow rate, the numerical solution method must be sufficiently stable. To test the ability of KMC-FB to handle this situation, we artificially constructed a very-low-flow case: a 61-pin test with wire-wrapper spacers cooled by an LBE. All fuel rods were 13.1 mm in diameter and arranged in a triangular pitch within a hexagonal duct. The wire spacers had a diameter of 1.4 mm and a pitch of 210 mm. All the liquid LBE flowed from the bottom to the top and passed through a 350 mm inlet zone and a 1000 mm heating zone. The inlet flow was 0.5797 kg/s at the mass flow rate and 283.5 ℃. All fuel rods were heated uniformly with a total thermal power of 11.39 kW. The test section was modeled with 100 axial nodes in the heated zone and 30 axial nodes in the inlet zone (total axial nodes = 130, each axial level = 10 mm). Figure 5(a) Figure 5(a) shows the convergence residuals in the calculation process, in which the red dotted line is the convergence residual limit set by the user. As shown in Fig. 5(a) Fig. 5(a) , with the improved novel solution method, the residual error of each physical quantity decreased steadily as the iteration progressed and finally reached the limit set by the user. The subchannel number and simulation results by KMC-FB for the 61-pin low-flow test are depicted in Fig. 5(b) and (c) Fig. 5(b) and (c) . The inlet flow rate was approximately 0.01 m/s and the minimum flow rate was approximately 0.008 m/s located at the edge channels (#1 and #18) owing to their smaller volume power density and higher frictional resistance. Multiple blockages 3 Based on the worst-case hypothesis, we constructed multi-axial level and multi-subchannel flow blockage accident with a 19-pin assembly, and performed a simulation. Figure 6(a) Figure 6(a) shows a schematic diagram of the blockage position, which had four axial planar blockages and six central subchannel blockages. The six blockage subchannels were distributed in four axial positions with an axial interval of 5 cm. The rod bundle was cooled by LBE, which flowed from the bottom to the top and passed through a 100 mm inlet zone and a 600 mm heating zone. All fuel rods were 10 mm in diameter and arranged in a triangular pitch within a hexagonal duct. The simulation results obtained using KMC-FB for the 19-pin multi-blockage test are shown in Figs. 6(b) and (c) Figs. 6(b) and (c) . The obvious backflow phenomenon was captured by KMC-FB, and the maximum downward flow rate was approximately 0.1 m/s. A temperature peak occurred near each blockage, and the temperature gradually increased axially with the increase in the blockage numbers, which caused the overall temperature profile to present a trapezoid in the flow direction, resulting in the outlet temperature being noticeably uneven. Preliminary validation 2 NACIE-UP blockage experiment 3 In the framework of the research activities planned to support the development of the ALFRED reactor, ENEA performed experimental tests to simulate the thermal-hydraulic behavior of a pin fuel bundle that is LBE-cooled in a flow blockage accident scenario, which was performed at the Natural Circulation Experiment-Upgraded) (NACIE-UP) facility loop with a heat source (a BFPS) consisting of 19 electrical pins. The experimental campaign focused on different stationary tests characterized by three different mass flow rates (4, 8, and 12 kg/s) and five fundamental degrees of blockage. In this study, as validation cases, we chose a normal condition (without blockage, referred to as BFPS-4-0) and a sector blockage case (BFPS-4-1). The instrumentation in the NACIE-UP facility allowed us to directly measure the nominal power of the LBE mass flow rate of the LBE temperatures in the loop, localized wall, and subchannel temperatures in the bundle. In this study, we mainly compared the differences in the values of temperatures of the subchannel and the wall surface of the electric heating rod produced by the code and the experiment. Figure 7 Figure 7 represents a top-view section of the BFPS in the case of sector blockage, and the position of the thermocouples. The key geometric parameters and boundary conditions of the experiments are listed in Table 1 Table 1 . The simulation results produced by KMC-FB for the BFPS tests are depicted in Fig. 8(a-e) Fig. 8(a-e) , where the x-coordinate represents the axial height from the entrance of the heated zone (mm), the y-coordinate represents the temperature and 1.5D, 2.5D, and 3D in the legend represent the one-dimensional heat conduction model, circumferential sub-node model, and three-dimensional heat conduction model used in KMC-FB, respectively, which were previously mentioned in Sect. 2.2. The experimental data are displayed found in Fig. 4 Fig. 4 , where SB2 represents the center of subchannel B2. Figure 8(a-d) Figure 8(a-d) , which represents the experiment without a blockage, shows that the calculation results of the code closely followed the experimental data, and the maximum relative error did not exceed 1%. Owing to the uniform heating power, the wall temperature of the rod bundle increased linearly with height. The results in Fig. 8(a)–(d) Fig. 8(a)–(d) show that the calculation results of the 1.5D, 2.5D, and 3D models were basically the same: because the nonuniformity of the fluid temperature around these rods was small, the circumferential subnode model and the three-dimensional heat conduction model showed no obvious effect. Figure 8(e) Figure 8(e) shows the distribution of the temperature of subchannel SB2 along the height of the flow path. The figure shows that the profile shape and the magnitude of temperature described by KMC-FB are in good agreement with the experimental results, and the maximum relative errors of all three models do not exceed 1.1%: their results are basically the same. The BFPS-4-1 case corresponds to approximately 1/6 of the flow area being fully blocked, and the blockage solid being a thin plate located at the beginning of the heated zone. Figure 9(a)–(d) Figure 9(a)–(d) shows the temperature distributions of the #1, #2, #5, and #15 pins along the axial direction. Because pins #1, #5, and #15 were all located in the blocked flow area, the axial temperature distribution showed a similar trend, but the specific value and change range were different, which was closely related to the pin position and the degree of blockage of the circumferential subchannels. Generally, the linear profile was maintained for pin 2, which faced an unblocked sector E, while all the other pins exhibited the behavior of a maximum, a local minimum, followed by a linear trend. From the perspective of the fuel rod models in the blockage area, the average effect of the one-dimensional fuel rod model severely underestimated the temperature peak of the pins behind the blockage. The three-dimensional fuel rod model was more accurate for predicting the local peak temperature under severe asymmetric conditions. However, at 30–100 mm, as shown in Fig. 9(a) Fig. 9(a) , the temperature predicted by KMC-FB showed a certain deviation from the experimental value, which may be due to the weak prediction effect of the corner channel. For the coolant temperature ( Fig. 9(e) Fig. 9(e) ), the calculation results produced by KMC-FB were generally in good agreement with the experimental results, and the maximum deviation did not exceed 3.1%. For the three different fuel rods, the three-dimensional fuel rod model was the closest to the experimental results, which proved that the circumferential heat flux density of the fuel rod was no longer uniformly distributed in this case, and was closer to the real experimental environment. However, we observed a certain deviation in the calculation results in the recirculation zone, mainly because (1) the blockage model was not accurate, resulting in a deviated hydraulic behavior; and (2) the lateral velocity cross-term <math id="M16"><mrow><msub><mi>c</mi><mtext>s</mtext></msub><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>n</mi><mi>i</mi></mrow></munderover><mfrac><mi>S</mi><mrow><msub><mi>l</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mfrac><mrow><mo>〈</mo><mrow><mi>ρ</mi><msup><mi>v</mi><mn>2</mn></msup></mrow><mo>〉</mo></mrow><mtext>cos</mtext><mi>Δ</mi><msup><mi>β</mi><mo>′</mo></msup></mrow></math> is ignored in the subchannel approach, so the calculation of the recirculation zone behind the blockage was not sufficiently accurate. KIT-THEADES blockage experiment 3 Within the European project MAXSIMA, the thermal-hydraulic effects of internal flow blockages in a heated rod bundle cooled by lead–bismuth eutectic (LBE) were investigated using a test section installed in a vertical test port of the existing THEADES loop facility at KIT-KALLA. In this study, we considered two blockage scenarios: a small blockage covering one subchannel, and a large blockage covering the six central subchannels, as shown in 10(a) 10(a) , in which the six gray subchannels represent the Block-center-C6 case, the black central subchannel represents the Block-C1 case, and the black edge subchannel represents the Block-E1 case. In this part of the study, we used KMC-FB to simulate these two blockage scenarios. Because there were too few thermocouples in the coolant subchannel (only three), we mainly focused on the influence of different heat conduction models on the wall temperature of the rod bundles; the simulation results of the coolant temperature were produced by KMC-FB using the three-dimensional heat conduction model. The key geometric parameters and boundary conditions of these cases are listed in Table 1 Table 1 . Figure 10(b) Figure 10(b) compares the temperature of the five thermocouples at the axial center in the blockage zone and the values predicted by the code under the first blockage scenario. We found that the temperature of points A-E calculated by the one-dimensional fuel rod model was basically the same, indicating the model could not reflect the circumferential unevenness caused by the blockage around the fuel rod. In addition, the circumferential subnode model adopts the subnode assumption and does not consider the circular conduction heat; therefore, the circumferential nonuniformity was too large, and the prediction of node B was obviously higher. Generally, the 3D fuel rod model more closely followed the experimental data, and the maximum deviation did not exceed 9%. However, the code did not accurately predict the temperature at points A and C, which may have been due to the local effect of the wire spacer, but the code adopts a simplified C-T correlation and can only consider the average effect of a wire spacer. Figure 10(c) Figure 10(c) compares the temperature in the blocked subchannel, where the distance is expressed as a fraction of the wire pitch ( H = 328 mm). The predicted value of the code was accurate, with a maximum deviation of 2%. Fig. 10(d),(e) Fig. 10(d),(e) shows the simulation results of the second small blockage. Similar to the first case, the code predicted an accurate value with a maximum deviation of 4% at points C and F, with a maximum deviation of 6% in the subchannel temperature. The simulation results produced by KMC-FB for the second scenario (large blockage) are depicted in Figure 10 (f) and (g) Figure 10 (f) and (g) . Because the subchannels around fuel rod #1 were completely blocked, their circumferential temperature was basically the same. Therefore, the values predicted by the three-dimensional fuel rod heat conduction model and the one-dimensional heat conduction model were similar. However, the temperature of the circumferential subnode model at points B and C was obviously higher, mainly due to ignoring the circumferential heat conduction. As for the coolant temperature, the calculation results of the code were in good agreement with the experiment, and the deviation did not exceed 5%. Conclusion 1 The analysis of reactor blockage accidents is an indispensable part of research on reactor safety analysis, especially for liquid metal fast reactors that are featured by a compact reactor core layout and high power density. In this study, we developed a subchannel code called KMC-FB to analyze blockage accidents in a reactor core. The code uses a time-dependent, four-equation, single-phase flow model together with a 3D heat conduction model for the fuel rod, which is solved by numerical methods based on the FDM with a staggered mesh. In addition, we proposed and implemented a novel and robust numerical solution method for the flow field solution in KMC-FB. We verified the code to demonstrate its strong ability to handle low flow rates and circulating flow cases. The first case considered was a very low flow case (~0.01 m/s) with a 61-pin bundle, and KMC-FB successfully calculated its flow velocity distribution. The second case was a multi-level and -subchannel flow blockage accident with a 19-pin assembly, and KMC-FB successfully simulated the locally low flow rate and reverse-flow phenomenon. We performed a preliminary validation by comparing the results produced by KMC-FB with the experimental data of the NACIE-UP BFPS blockage test and the KIT 19ROD blockage test, which included normal working conditions (no blockage), planar blockage, small solid blockage, and large solid blockage. We obtained agreement with the experimental data, with a relative deviation interval of [1%, 9%], which indicated that the code is capable of performing flow blockage analysis for LMFRs. In addition, we used a three-dimensional heat conduction model together with the one-dimensional heat conduction model and the circumferential sub-node model for comparison. The verification results showed that (1) under normal operating conditions, the accuracy of the one-dimensional heat conduction model is sufficient for subchannel thermal-hydraulic analysis. (2) In flow blockage accidents or similar conditions that can cause large circumferential and axial temperature gradients of the fuel rods, the three-dimensional heat conduction model is more accurate. (3) The fuel rod circumferential subnode model can sometimes predict accurately, but at other times, it overestimates the fuel rod temperature, which depends on how much it is blocked in the surrounding subchannels. Future work will focus on several aspects to improve simulation accuracy, which will involve (1) searching for more accurate hydraulic and heat transfer correlations when the flow is blocked or the crossflow is dominant to improve the calculation accuracy of the subchannel code in the recirculation zone; and (2) constructing a more accurate wire-wrapper spacer model. References T-T. Haileyesus , P-M Sara , P. Werner et al. , A review of models for the sodium boiling phenomena in sodium-cooled fast reactor subassemblies . ASME J. Nucl. Rad. Sci. 8 ( 1 ), 011305 ( 2021 ). doi: 10.1115/1.4051066 http://doi.org/10.1115/1.4051066 A review of models for the sodium boiling phenomena in sodium-cooled fast reactor subassemblies W.P. Chang , K.S. Ha , S.D. Suk et al. , A comparative study of the MATRA-LMR-FB calculation with the SABRE result for the flow blockage accident in the sodium cooled fast reactor . Nucl. Eng. Des. 241 , 5225 - 5237 ( 2011 ). doi: 10.1016/j.nucengdes.2011.08.039 http://doi.org/10.1016/j.nucengdes.2011.08.039 A comparative study of the MATRA-LMR-FB calculation with the SABRE result for the flow blockage accident in the sodium cooled fast reactor W.P. Feng , X. Zhang , L.K. Cao et al. , Development of oxygen/corrosion product mass transfer model and oxidation-corrosion model in the lead-alloy cooled reactor core . Corros. Sci. 190 , 109708 ( 2021 ). doi: 10.1016/j.corsci.2021.109708 http://doi.org/10.1016/j.corsci.2021.109708 Development of oxygen/corrosion product mass transfer model and oxidation-corrosion model in the lead-alloy cooled reactor core B. Jacopo , P.L. Eric , C. Kenneth et al. , Studies of polonium removal from molten lead-bismuth for lead-alloy-cooled reactor applications . Nucl. Technol. 147 , 406 - 417 ( 2004 ). doi: 10.13182/NT04-A3539 http://doi.org/10.13182/NT04-A3539 Studies of polonium removal from molten lead-bismuth for lead-alloy-cooled reactor applications S.Y. Wei ; C.L. Wang , W.X. Tian et al. , Research progress in key thermal-hydraulic issue of lead-based fast reactor . Atomic Energy Science and Technology. 53 , 326 ( 2019 ). doi: 10.7538/yzk.2018.youxian.0335 http://doi.org/10.7538/yzk.2018.youxian.0335 (in Chinese) Research progress in key thermal-hydraulic issue of lead-based fast reactor M.H. Fontana , T.S. Kress , T.S. Parsly et al. , Effect of partial blockages in simulated LMFBR fuel assemblies. No. ORNL-TM-4324 . Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United States) , 1973 . URL https://www.osti.gov/biblio/4357059 https://www.osti.gov/biblio/4357059 Effect of partial blockages in simulated LMFBR fuel assemblies. No. ORNL-TM-4324 J.T. Han , M.H. Fontana , Blockages in LMFBR fuel assemblies: A review . ( 1977 ). URL https://www.osti.gov/biblio/5232755 https://www.osti.gov/biblio/5232755 Blockages in LMFBR fuel assemblies: A review J. Kobayashi , T. Isozaki , M. Tanaka et al. , Study on temperature field in porous blockage in a fuel subassembly. 37-pin bundle sodium experimen t . No. JNC-TN--9400-2000-025 . Japan Nuclear Cycle Development Inst. , 1999 . URL https://www.osti.gov/etdeweb/biblio/20239902 https://www.osti.gov/etdeweb/biblio/20239902 Study on temperature field in porous blockage in a fuel subassembly. 37-pin bundle sodium experiment. No. JNC-TN--9400-2000-025 J.M. Creer , D.S. Rowe , J.M. Bates et al. , Effects of sleeve blockages on axial velocity and intensity of turbulence in an unheated 7 × 7 rod bundle.[PWR]. No. BNWL-1965 . Battelle Pacific Northwest Labs., Richland, Wash. (USA) , 1976 . URL https://www.osti.gov/biblio/7290924 https://www.osti.gov/biblio/7290924 Effects of sleeve blockages on axial velocity and intensity of turbulence in an unheated 7 × 7 rod bundle.[PWR]. No. BNWL-1965 J. Pacio , M. Daubner , F. Fellmoser et al. , Heat transfer experiment in a partially (internally) blocked 19-rod bundle with wire spacers cooled by LBE . Nucl. Eng. Des. 330 , 225 - 240 ( 2018 ). doi: 10.1016/j.nucengdes.2018.01.034 http://doi.org/10.1016/j.nucengdes.2018.01.034 Heat transfer experiment in a partially (internally) blocked 19-rod bundle with wire spacers cooled by LBE R. Marinari , I. Di Piazza , M. Tarantino et al. , Blockage fuel pin simulator experiments and simulation . Nucl. Eng. Des. 353 , 110215 ( 2019 ). doi: 10.1016/j.nucengdes.2019.110215 http://doi.org/10.1016/j.nucengdes.2019.110215 Blockage fuel pin simulator experiments and simulation I.D. Piazza , F. Magugliani , M. Tarantino et al. , A CFD analysis of flow blockage phenomena in ALFRED LFR demo fuel assembly . Nucl. Eng. Des. 276 , 202 - 215 ( 2014 ). doi: 10.1016/j.nucengdes.2014.05.033 http://doi.org/10.1016/j.nucengdes.2014.05.033 A CFD analysis of flow blockage phenomena in ALFRED LFR demo fuel assembly H.G. Kim , B. Hwange , J.K. Yung et al. , Investigations of single-phase flow mixing characteristics in a wire-wrapped 37-pin bundle for a sodium-cooled fast reactor . Ann. Nucl. Energy 87 , 541 - 546 ( 2016 ). doi: 10.1016/j.anucene.2015.09.022 http://doi.org/10.1016/j.anucene.2015.09.022 Investigations of single-phase flow mixing characteristics in a wire-wrapped 37-pin bundle for a sodium-cooled fast reactor A. Kaliatka , E. Uspuras , Development and testing of RBMK-1500 model for BDBA analysis employing RELAP/SCDAPSIM code . Ann. Nucl. Energy 35 , 977 - 992 ( 2008 ). doi: 10.1016/j.anucene.2007.11.015 http://doi.org/10.1016/j.anucene.2007.11.015 Development and testing of RBMK-1500 model for BDBA analysis employing RELAP/SCDAPSIM code H.T. Mehdi , T. Saman , M.S. Seyyed et al. , Forced reflood modeling in a 2× 2 rod bundle with a 90% partially blocked region . Ann. Nucl. Energy 131 , 425 - 432 ( 2019 ). doi: 10.1016/j.anucene.2019.04.019 http://doi.org/10.1016/j.anucene.2019.04.019 Forced reflood modeling in a 2× 2 rod bundle with a 90% partially blocked region Z.R. Zou , Study on instability characteristics of natural circulation lead-cooled fast reactor under startup and transient operation. (Master Thesis) . University of Science and Technology of China , Hefei, China , 2021 . (in Chinese) Study on instability characteristics of natural circulation lead-cooled fast reactor under startup and transient operation. (Master Thesis) W.S. Kim , Y.-G. Kim , Y.J. Kim , A subchannel analysis code MATRA-LMR for wire wrapped LMR subassembly . Ann. Nucl. Energy 29 , 303 - 321 ( 2002 ). doi: 10.1016/S0306-4549(01)00041-X http://doi.org/10.1016/S0306-4549(01)00041-X A subchannel analysis code MATRA-LMR for wire wrapped LMR subassembly J.D. Macdougall , J.N. Lillington , The SABRE code for fuel rod cluster thermohydraulics . Nucl. Eng. Des. 82 , 171 - 190 ( 1984 ). doi: 10.1016/0029-5493(84)90210-3 http://doi.org/10.1016/0029-5493(84)90210-3 The SABRE code for fuel rod cluster thermohydraulics C.W. Stewart , C.L. Wheeler , R.J. Cena et al. , COBRA-IV: The model and the method. No. BNWL-2214 . Pacific Northwest Lab., Richland, WA (USA) , 1977 . URL https://www.osti.gov/biblio/5358588 https://www.osti.gov/biblio/5358588 COBRA-IV: The model and the method. No. BNWL-2214 V. Kriventsev , A. Rineiski , W. Maschek , Application of safety analysis code SIMMER-IV to blockage accidents in FASTEF subcritical core . Ann. Nucl. Energy 64 , 114 - 121 ( 2014 ). doi: 10.1016/j.anucene.2013.09.045 http://doi.org/10.1016/j.anucene.2013.09.045 Application of safety analysis code SIMMER-IV to blockage accidents in FASTEF subcritical core M. Di , Valentino , I. Uwe , S. Victor , Validation and application of the system code ATHLET-CD for BWR severe accident analyses . Nucl. Eng. Des. 307 , 284 - 298 ( 2016 ). doi: 10.1016/j.nucengdes.2016.07.013 http://doi.org/10.1016/j.nucengdes.2016.07.013 Validation and application of the system code ATHLET-CD for BWR severe accident analyses A.S. Ekariansyah , P.H. Endiah , S. Sudarmono , Relap5 simulation for severe accident analysis of rsg-gas reactor . Jurnal Teknologi Reaktor Nuklir Tri Dasa Mega. 20 , 23 - 34 ( 2018 ). doi: 10.17146/tdm.2018.20.1.4040 http://doi.org/10.17146/tdm.2018.20.1.4040 Relap5 simulation for severe accident analysis of rsg-gas reactor S.U. Khan , U.K. Shahab , M.J. Peng , Flow blockage accident or loss of flow accident by using comparative approach of NK/TH coupling codes and RELAP5 code . Ann. Nucl. Energy 64 , 311 - 319 ( 2014 ). doi: 10.1016/j.anucene.2013.10.010 http://doi.org/10.1016/j.anucene.2013.10.010 Flow blockage accident or loss of flow accident by using comparative approach of NK/TH coupling codes and RELAP5 code L.Z. Wang , M. Jin , G.Y. Zhang et al. , Transient thermal-hydraulic analysis for fuel assembly blockage accident by coupling CFD code with the neutron kinetics model . Ann. Nucl. Energy 137 , 107138 ( 2020 ). doi: 10.1016/j.anucene.2019.107138 http://doi.org/10.1016/j.anucene.2019.107138 Transient thermal-hydraulic analysis for fuel assembly blockage accident by coupling CFD code with the neutron kinetics model Z.X. Gu ., Q.X. Zhang , Y. Gu et al. , Verification of a self-developed CFD-based multi-physics coupled code MPC-LBE for LBE-cooled reactor . Nucl. Sci. Tech. 32 , 52 ( 2021 ). doi: 10.1007/s41365-021-00887-x http://doi.org/10.1007/s41365-021-00887-x . Verification of a self-developed CFD-based multi-physics coupled code MPC-LBE for LBE-cooled reactor L.K. Cao , G.L. Yang , and H.L. Chen . Transient sub-channel code development for lead-cooled fast reactor using the second-order upwind scheme . Prog. Nucl. Energy 110 , 199 - 212 ( 2019 ). doi: 10.1016/j.pnucene.2018.09.022 http://doi.org/10.1016/j.pnucene.2018.09.022 Transient sub-channel code development for lead-cooled fast reactor using the second-order upwind scheme K.S. Ha , H.Y. Jeong , W.P. Chang et al. , Development of the MATRA-LMR-FB for flow blockage analysis in a LMR . Nucl. Eng. Des. 41 , 797 - 806 ( 2009 ). doi: 10.5516/NET.2009.41.6.797 http://doi.org/10.5516/NET.2009.41.6.797 Development of the MATRA-LMR-FB for flow blockage analysis in a LMR D.B. Spalding , A novel finite difference formulation for differential expressions involving both first and second derivatives . Int. J. Numer. Meth. Eng. 4 , 551 - 559 ( 1972 ). doi: 10.1002/nme.1620040409 http://doi.org/10.1002/nme.1620040409 A novel finite difference formulation for differential expressions involving both first and second derivatives I. E. Idelchik . Handbook of Hydraulic Resistance, Revised and Augmented . ( 2008 ). ISBN Online: 978-1-56700-340-6 Handbook of Hydraulic Resistance, Revised and Augmented J. K. Fink and L. Leibowitz , Thermodynamic and transport properties of sodium liquid and vapor. No. ANL--RE-95/2 . Argonne National Lab. , 1995 . URL https://www.osti.gov/biblio/94649-thermodynamic-transport-properties-sodium-liquid-vapor https://www.osti.gov/biblio/94649-thermodynamic-transport-properties-sodium-liquid-vapor Thermodynamic and transport properties of sodium liquid and vapor. No. ANL--RE-95/2 C. Fazio , V.P. Sobolev , A. Aerts et al. , Handbook on lead-bismuth eutectic alloy and lead properties, materials compatibility, thermal-hydraulics and technologies-2015 edition. No. NEA--7268 . Organisation for Economic Co-Operation and Development , 2015 . URL https://inis.iaea.org/search/search.aspx?orig_q=RN:46133907 https://inis.iaea.org/search/search.aspx?orig_q=RN:46133907 Handbook on lead-bismuth eutectic alloy and lead properties, materials compatibility, thermal-hydraulics and technologies-2015 edition. No. NEA--7268 P.L.T. Brian , A finite‐difference method of high‐order accuracy for the solution of three‐dimensional transient heat conduction problems . AIChE J. 7 , 367 - 370 ( 1961 ). doi: 10.1002/aic.690070305 http://doi.org/10.1002/aic.690070305 A finite‐difference method of high‐order accuracy for the solution of three‐dimensional transient heat conduction problems All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xiao Luo, Lian-Kai Cao, Wen-Pei Feng, and Hong-Li Chen. The first draft of the manuscript was written by Xiao Luo and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01010-4 2021-09-22 2022-01-04 2022-01-12 2022-03-09 2022-05-16 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01010-4 http://dx.doi.org/10.1007/s41365-022-01010-4 . 2022 Abbreviations <math id="M17"><mi>A</mi></math> : Flow area of axial control cell, m 2 <math id="M18"><mrow><msub><mi>A</mi><mtext>T</mtext></msub></mrow></math> : Flow area of the lateral control cell, m 2 <math id="M19"><mrow><msub><mi>c</mi><mtext>s</mtext></msub></mrow></math> : Incomplete coupling coefficient <math id="M20"><mrow><msub><mi>D</mi><mtext>e</mtext></msub></mrow></math> : Hydraulic diameter of flow channel, m <math id="M21"><mi>f</mi></math> : Hydraulic friction factor <math id="M22"><mrow><msub><mi>f</mi><mtext>T</mtext></msub></mrow></math> : Turbulent momentum factor (dimensionless empirical parameter) <math id="M23"><mi>g</mi></math> : Acceleration of gravity <math id="M24"><mi>h</mi></math> : Specific enthalpy of coolant, <math id="M25"><mrow><mtext>J</mtext><mo>⋅</mo><msup><mrow><mtext>kg</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> H : Wire pitch, m <math id="M26"><mi>K</mi></math> : Shape resistance coefficient of axial flow (dimensionless empirical parameter) <math id="M27"><mrow><msub><mi>K</mi><mtext>G</mtext></msub></mrow></math> : Crossflow resistance coefficient (dimensionless empirical parameter) <math id="M28"><mrow><msub><mi>l</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></math> : Centroid distance between adjacent subchannels, m <math id="M29"><mrow><mi>N</mi><mi>u</mi></mrow></math> : Nusselt number <math id="M30"><mi>P</mi></math> : Hydrostatic pressure, Pa Pe : Peclet number <math id="M31"><mrow><msup><mi>q</mi><mrow><mo>'</mo><mo>'</mo></mrow></msup></mrow></math> : Heat flux density on fuel rod surfaces, <math id="M32"><mrow><mtext>W</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>2</mn></mrow></msup></mrow></math> <math id="M33"><mrow><msub><mi>q</mi><mtext>v</mtext></msub></mrow></math> : Power density of unit volume, <math id="M34"><mrow><mtext>W</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>2</mn></mrow></msup></mrow></math> r 1 : Dimensionless empirical coefficient R f : Equivalent heat resistance of convective heat transfer between the fluid and fuel rods, K/W Re : Reynold number R bl : Heat resistance of the blockage, <math id="M35"><mrow><mtext>K</mtext><mo>⋅</mo><msup><mtext>W</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M36"><mi>S</mi></math> : Length of the width gap between adjacent subchannels, m <math id="M37"><mi>t</mi></math> : Time, s <math id="M38"><mi>v</mi></math> : Lateral flow velocity of coolant, <math id="M39"><mrow><mtext>m</mtext><mo>⋅</mo><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M40"><mi>V</mi></math> : Crossflow mass flow rate per unit length, <math id="M41"><mrow><mtext>kg</mtext><mo>⋅</mo><mtext>m</mtext><mo>⋅</mo><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M42"><mrow><msup><mi>V</mi><mo>'</mo></msup></mrow></math> : Crossflow mass flow rate per unit length caused by turbulent mixing between adjacent subchannels, <math id="M43"><mrow><mtext>kg</mtext><mo>⋅</mo><mtext>m</mtext><mo>⋅</mo><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M44"><mi>w</mi></math> : Axial flow velocity of coolant, <math id="M45"><mrow><mtext>m</mtext><mo>⋅</mo><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M46"><mi>x</mi></math> : Lateral position, m <math id="M47"><mi>z</mi></math> : Axial height, m <math id="M48"><mi>ρ</mi></math> : Density of coolant, <math id="M49"><mrow><mtext>kg</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>3</mn></mrow></msup></mrow></math> <math id="M50"><mi>θ</mi></math> : Angle of flow direction relative to vertical, rad <math id="M51"><mrow><msub><mi>λ</mi><mrow><mtext>bl</mtext></mrow></msub></mrow></math> : Thermal conductivity of the blockage material, <math id="M52"><mrow><mtext>W</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>⋅</mo><msup><mtext>K</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M53"><mrow><msub><mi>λ</mi><mtext>f</mtext></msub></mrow></math> : Thermal conductivity of the coolant, <math id="M54"><mrow><mtext>W</mtext><mo>⋅</mo><msup><mtext>m</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>⋅</mo><msup><mtext>K</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> <math id="M55"><mrow><mi>Δ</mi><mi>t</mi></mrow></math> : Time step, s <math id="M56"><mrow><mi>Δ</mi><msup><mi>β</mi><mo>′</mo></msup></mrow></math> : Angle between the gap under investigation and the direction of the adjacent gap and radius Subscripts i : Subchannel number i j : Subchannel number j x : x -direction in the Cartesian coordinate system y : y -direction in the Cartesian coordinate system z : z -direction in the Cartesian coordinate system Superscripts bar: Average value of the physical properties n : Last time step *: Donor cell quantity <math id="M57"><mrow><mo>'</mo><mo>'</mo></mrow></math> : Per unit cell Nuclear Science and Techniques 03 33