Development and application of a multi-physics and multi-scale coupling program for lead-cooled fast reactor

Development and application of a multi-physics and multi-scale coupling program for lead-cooled fast reactor 1 first-author 1 1 1 1 2 corresp Zhao Chen, E-mail address: chen-zhao@cgnpc.com.cn Zhao Chen, E-mail address: chen-zhao@cgnpc.com.cn chen-zhao@cgnpc.com.cn 1 corresp Hong-Li Chen E-mail address: hlchen1@ustc.edu.cn; Hong-Li Chen E-mail address: 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 China Nuclear Power Technology Research Institute Co. Ltd, Shenzhen 518000, China China Nuclear Power Technology Research Institute Co. Ltd, Shenzhen 518000, China In this study, a multi-physics and multi-scale coupling program, Fluent/KMC-sub/NDK, was developed based on the user-defined functions (UDF) of Fluent, in which the KMC-sub code is a sub-channel thermal-hydraulic code and the NDK code is a neutron diffusion code. The coupling program framework adopt the "master-slave" mode, in which Fluent is the master program while NDK and KMC-sub are coupled internally and compiled into the dynamic link library (DLL) as slave codes. The domain decomposition method was adopted, in which the reactor core was simulated by NDK and KMC-sub, while the rest of the primary loop was simulated using Fluent. A simulation of the reactor shutdown process of M 2 LFR-1000 was carried out using the coupling program, and the code-to-code verification was performed with ATHLET, demonstrating a good agreement, with absolute deviation was smaller than 0.2%. The results show an obvious thermal stratification phenomenon during the shutdown process, which occurs 10 s after shutdown, and the change in thermal stratification phenomena is also captured by the coupling program. At the same time, the change in the neutron flux density distribution of the reactor was also obtained. Multi-physics and multi-scale coupling method User-defined functions Dynamic link library Thermal stratification Lead-cooled fast reactor Introduction 1 Safety analysis is a key tool for proving the safety of nuclear power plants during the design phase. Usually, the safety analysis of nuclear power plants is carried out using a system analysis code, which adopts the one-dimensional approximate simplification method to model the real plant system. However, there is three-dimensional flow and thermal stratification in the pool structure of lead-cooled fast reactors (LFRs), which cannot be simulated by system codes, especially in the accident state. On the other hand, CFD tools such as Fluent can capture the complex three-dimensional flow and heat transfer phenomenon in detail, which can carry out accurate numerical simulation analysis of the plenum in LFRs. However, they may require very large computational effort to simulate such important phenomena. In this case, a multiscale coupling method is proposed to meet this demand. Aumille et al. proposed a RELAP-3D/CFD coupling methodology and performed a proof-of-principle calculation based on the coupling code [ 1 1 - 2 2 1-2 ]. Unlike the explicit CFD coupling method, this coupling method is numerically stable, provided the Courant limit is not violated in RELAP5-3D or at the coupling locations. Gesellschaft für Anlagen- und Reaktorsicherheit (GRS) proposed an explicit coupling method and developed a 1D system/3D-CFD coupling program based on the self-developed system analysis code ATHLET [ 3 3 ]. The preliminary comparative calculations are complicated between the coupling program and ATHLET stand-alone, which showed good agreement. Similarly, the Royal Institute of Technology (KTH) also proposed a coupling method between the RELAP5/STAR-CCM+ coupling program and compared its simulation results with GRS directly [ 4 4 ]. Experiment T01 of the TALL-3D experimental facility was chosen as the benchmark, and the simulation results were directly compared with the GRS (ATHLET–ANSYS CFX) coupling code. The French Alternative Energies and Atomic Energy Commission (CEA) coupled the 3D computational fluid dynamics (CFD) code TRIO_U with the thermal-hydraulic system analysis code CATHARE to perform a simulation of the natural circulation test of the French prototype sodium fast reactor (SFR) Phénix, which closely followed the available experimental data and the results of a single CATHARE calculation [ 5 5 ]. In addition, we successfully completed the ATHLET-OpenFOAM coupling code application in our previous work and simulated the unprotected loss of flow accident (ULOF) process for LFR[ 6 6 ]. However, in the above coupled programs, for the reactor core, the lumped parameter method is usually used to simplify the calculation in the system analysis code, which cannot accurately simulate the internal thermal-hydraulic behaviors. The sub-channel code can correctly and quickly simulate thermal-hydraulic phenomena inside the reactor core and can assess the thermal-hydraulic safety characteristics. In previous works, Chi Wang et al. developed a code by coupling FLUENT with the sub-channel thermal-hydraulic code (KMC-sub)[ 7 7 ]. In addition, the point source neutron-kinetic model was used to simplify the reactor core, which could not evaluate its spatial performance. For example, when the reactor is shut down and the control rod is inserted, the nuclear thermal coupling effect is strong and the local effect is obvious, which cannot be correctly captured by using the point source neutron-kinetic model. To solve this problem, a three-dimensional neutron diffusion code should be used in the reactor core to replace the point-source neutron-kinetic model. In our previous work, a three-dimensional neutron diffusion code was developed [ 8 8 ]. In addition, there are also some nuclear thermal coupling codes, such as RMC/COBRA (a Monte Carlo and subchannel coupling code for PWR multi-physics analysis) [ 9 9 , 10 10 ], GENFOAM (a multi-physics solver based on OpenFOAM) [ 11 11 - 13 13 11-13 ], and FLUENT/PK (a CFD/neutron kinetics coupled code for LFR transient safety analysis) [ 14 14 ]. Therefore, it is necessary to develop a multi-physics and multiscale coupling program to simulate the internal physical state and thermal-hydraulic characteristics of pool-type reactors more accurately. In general, a coupling program FLUENT/KMC-sub/NDK was developed in this study and used to simulate the transient accident of the entire primary loop of pool type LFRs. Unlike the previous work mentioned above, a more refined coupling model is used in this work, which adopts a three-dimensional neutron diffusion code (NDK) to handle reactor core neutronic physics. In addition, the CFD program (FLUENT) was used as the master code to capture the three-dimensional thermal-hydraulic phenomenon in the reactor pool, and the sub-channel thermal-hydraulic code (KMC-sub) was used to simulate the reactor core. A simulation of the reactor shutdown process of M 2 LFR-1000 was carried out using the coupling program, and the code-to-code verification was performed with ATHLET, demonstrating good agreement. Codes to be coupled 1 There are many commercial (generic) CFD codes and special CFD codes for nuclear safety applications in the nuclear field. Among them, ANSYS Fluent is widely used in nuclear engineering. To realize detailed three-dimensional flow and heat transfer calculations, using UDF [ 15 15 ] to realize data transfer between Fluent and other codes is a good choice. However, the smaller mesh sizes and larger mesh counts required by CFD codes make it difficult to simulate the entire system, especially for the reactor core, under current computer hardware conditions. The subchannel code can be used to model the entire reactor core, in which the conservation equations of continuity, momentum (both axial and lateral), and energy are solved in a control volume with radial and axial interconnections, and the correlations of flow resistance and heat transfer obtained from experiments or simulations are added. The inter-channel mixing model considering the transverse momentum balance and cross-flow between adjacent sub-channels was used to evaluate the thermal-hydraulic and coolant flow distribution in the assembly. In our previous work, a sub-channel code KMC-sub was developed for the thermal-hydraulic analysis of LFR [ 16 16 ]. In addition, it was coupled with Fluent and applied to M 2 LFR-1000, an LFR design proposed by USTC [ 17 17 ]. According to the neutron transport equation, assuming an isotropic angular distribution of neutron flux density, a typical neutron diffusion kinetic model is obtained, which can be used to calculate the reactor power distribution and reactivity changes. The multidimensional space-time neutron diffusion equation and the delayed neutron pioneer nuclear equilibrium equation are shown in equation (1). <math display="block" id="M1"><mtable columnalign="left"><mtr><mtd><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>g</mi></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mi>g</mi></msub><mo>∇</mo><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>r</mi><mo>,</mo><mi>g</mi></mrow></munder><mrow><msub><mi>ϕ</mi><mi>g</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>≠</mo><mi>g</mi></mrow></munder><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>s</mi><mo>,</mo><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></munder><mrow><msub><mi>ϕ</mi><mrow><mi>g</mi><mo>'</mo></mrow></msub></mrow></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>+</mo></mrow></mstyle><msub><mi>χ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>f</mi><mo>,</mo><mi>g</mi><mo>'</mo></mrow></munder><mrow><msub><mi>ϕ</mi><mrow><mi>g</mi><mo>'</mo></mrow></msub></mrow></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>N</mi><mi>D</mi></mrow></munderover><mrow><msub><mi>χ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle><mo>,</mo><mi>g</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>G</mi></mtd></mtr><mtr><mtd><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>β</mi><mi>i</mi></msub><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mi>ν</mi><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>f</mi><mo>,</mo><mi>g</mi><mo>'</mo></mrow></munder><mrow><msub><mi>ϕ</mi><mrow><mi>g</mi><mo>'</mo></mrow></msub></mrow></mstyle><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow></mstyle><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><mover accent="true"><mi>r</mi><mo stretchy="true">→</mo></mover><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>N</mi><mi>D</mi></mtd></mtr></mtable></math> A detailed explanation of these variables can be found in [ 8 8 ]. In this study, the two-group method is adopted, which means that thermal neutrons are treated as one group (energy < 0.1 eV) and the energy of the other energy group is greater than E c , where E c is the delimitation energy (~ 0.4 eV). This NDK code based on the neutron diffusion equation was developed. Its result was compared with a 2D-TWIGL benchmark and demonstrated good agreement [ 8 8 ]. For the transient diffusion equation, because it is a time-rigid equation, the implicit difference scheme is selected to solve the equation on the time scale. At each time step, it is reduced to a neutron diffusion equation with a fixed source term. Coupling approach 1 In this study, the UDF of Fluent was used to realize the data exchange between Fluent, KMC-sub, and NDK. The “master-slave” mode was selected for coupling, in which FLUENT was used as the "master program,” while the “slave code” KMC-sub and NDK coupling code were compiled into Dynamic Link Library (DLL), as shown in Fig. 1(a) Fig. 1(a) . The "master program” Fluent calls DLL file through UDF to realize data exchange and process control. Once the two programs are coupled under the control of Fluent, data exchange occurs in the "master program" boundary condition. In our previous work, an explicit sequential coupling scheme between Fluent and KMC-sub was implemented, and a steady case and transient case were completed [ 7 7 ]. To expand the scope of coupling subjects, a multiphysics and multi-scale coupling scheme is proposed in this paper, as shown in Fig. 1(b) Fig. 1(b) . The UDF can receive the data on the boundary from KMC-sub by calling the API functions in DLL files and sends the boundary data to the KMC-sub/NDK coupling code in each CFD time step, where the KMC-sub module transmits the temperature of the coolant, cladding, and pellet to the NDK module, while the NDK module transmits the spatial distribution of power to the KMC-sub module. The data transferred between the "master program" Fluent and the DLL is mainly the coolant temperature, pressure, and mass flow at the inlet and outlet of the core. To flatten the power distribution, the core assemblies are usually divided into several zones loaded with fuels of different enrichments. When Fluent transfers the data to KMC-sub at the boundary in a fuel zone, the temperature and pressure are the average values of the Fluent interface, while the mass flow is the sum of the Fluent interface, and KMC-sub calculates the flow mass rate distribution of each assembly in this zone, according to the flow-split model. The data transferred from KMC-sub to Fluent are transferred according to the division of fuel assembly zones, and the pressure and temperature on the Fluent interface of different zones were obtained from the average value of the number of KMC-sub assemblies in the corresponding zone, while the mass flow of the Fluent interface was obtained by the sum of the number of corresponding KMC-sub assemblies. In this study, the time-explicit method was used in the coupling scheme. The three coupling codes have their own time steps and converge residuals, and Fluent can control the time advancement of the KMC-sub/NDK module through the UDF to maintain the same calculation time as Fluent. The KMC-sub module and Fluent have the same time step, while the NDK module has a much smaller time step than Fluent, one-tenth of the Fluent time step. At the beginning of each time step, Fluent starts the iterative calculation. After convergence or reaching the maximum number of iteration steps, the data are transferred to the KMC module, then the KMC sub-module starts the iteration. After the KMC-sub module converges, the temperature parameters are transferred to the NDK module, and the NDK module performs the iterative calculation. When the NDK module converges and reaches the same time as the KMC-sub module, the NDK module transfers power to the KMC-sub module; otherwise, the NDK module continues to increase the number of time iterations. When the overall time of the three modules is the same, the KMC sub-module transfers the data to CFD, and the next iteration begins. In each step of the calculation, an absolute error limit and a specified number of steps were used for each code as their convergence conditions. The overall structure is an explicit sequential calculation, and the data are exchanged at the end of a time step. When the CFD reaches the set time or converges residuals, the overall calculation is terminated. Application of the coupled code Fluent/KMC-sub/NDK 1 M 2 LFR-1000 reactor 2 M 2 LFR-1000 is a typical pool-type lead-cooled fast reactor (LFR) proposed by the University of Science and Technology of China (USTC), which has completed the overall conceptual design and preliminary safety analysis [ 17 17 - 19 19 17-19 ]. The overall design of M 2 LFR-1000 is shown in Fig. 2 Fig. 2 , and the detailed design parameters of M 2 LFR are shown in Table 1 Table 1 . 1 1 Parameters 1 1 Value 1 1 Thermal power (MW th ) 1 1 1000 1 1 Primary coolant 1 1 Lead 1 1 Core inlet temp. (℃) 1 1 400 1 1 Core outlet temp. (℃) 1 1 480 1 1 Cooling method 1 1 Forced circulation 1 1 Number of total fuel assemblies 1 1 355 1 1 Number of pins in each fuel assembly 1 1 169 1 1 Assembly geometry/rod arrangement 1 1 Hexagon/triangle 1 1 Hexagonal wrapper apothem (mm) 1 1 97.5 1 1 Number of primary pumps 1 1 4 1 1 Number of steam generators (SG) 1 1 8 1 1 Secondary medium 1 1 Water 1 1 Inlet temp. of secondary loop (℃) 1 1 335 1 1 Outlet temp. of secondary loop (℃) 1 1 475 1 1 Equivalent diameter of reactor core (m) 1 1 4.048 1 1 Total core height (m) 1 1 2.15 Description of the shutdown process 2 To evaluate the performance and accuracy of the Fluent/KMC-sub/NDK coupling program, a reactor shutdown process simulation was carried out and was compared with the results obtained from ATHLET3.1, which is a system analysis code for nuclear reactors or facilities [ 20 20 ]. The transient boundary conditions for this case are given below. (1) First, the control rods and adjusting rods were not inserted. When shutting down, all control rods were inserted within 2 s. (2) The reactor mass flow rate remains unchanged, and the variations in the coolant temperature and mass flow rate at the outlet of the heat exchanger are input from the ATHLET results through the UDF. This is displayed in detail in Sect. 4.3. (3) The core decay heat power accounts for about 5% of normal power. ATHLET model and its input conditions 2 To evaluate the results of the coupled program more reasonably, ATHLET3.1 was used to calculate the mass flow rate and temperature parameters at the outlet of the heat exchanger during the shutdown transient to provide boundary conditions for the heat exchanger in the coupled program. The ATHLET node diagram for M 2 LFR-1000 is shown in Fig. 3(a) Fig. 3(a) , in which the reactor core is divided into four zones: the inner average fuel zone, inner hottest fuel zone, outer average fuel zone, and outer hottest fuel zone, which correspond to the inner and outer fuel assemblies in the core of M 2 LFR-1000, as shown in Fig. 3(b) Fig. 3(b) . In the ATHLET calculation model, the reactor was in steady-state operation before 0 s. From 0 s, 5100 pcm of negative reactivity is added within 2 s to simulate the reactor shutdown process, as shown in Fig. 4(a) Fig. 4(a) . During the shutdown process, the pump speed remains constant, and the changes in the core power and primary circuit mass flow are shown in Fig. 4(b) Fig. 4(b) . The change in the coolant temperature at the outlet heat exchanger is shown in Fig. 4(c) Fig. 4(c) . Model of coupling program 2 The domains of the three modules and their data exchange are shown in Fig. 5(a) Fig. 5(a) , which is divided into three domains, where the hot plenum and cold plenum are calculated by Fluent to simulate the thermal-hydraulic phenomenon, while the core is calculated by KMC-sub and NDK modules to simulate the thermal-hydraulic phenomenon and physical state. As shown in the C1 and C2 interface data mapping diagram ( Fig. 5(b) Fig. 5(b) ), the geometry of KMC-sub consists of many assemblies, which are divided into inner and outer fuel zones. To solve the mapping problem between the coarse mesh adopted by KMC-sub and the finer mesh in the CFD model, the following method was used in this study. The average value parameter of each assembly simulated by Fluent is transferred to the subchannel code when Fluent transfers data to KMC-sub, while KMC-sub returns lumped parameters of each assembly to Fluent. As for the reactor core, KMC-sub can calculate the spatial temperature distribution of the coolant, cladding, and pellet, and then transmits it to the NDK domain, and then updates the material’s absorption cross section according to the material’s temperature and calculates the distribution of neutron flux density. After renewing the neutron flux density, a new power distribution can be calculated and transmitted to the KMC-sub. As shown in Fig. 5(a) Fig. 5(a) , the cold plenum simulated by Fluent transmits the mass flow rate and temperature to the KMC-sub and NDK modules at the first interface, whereas the hot plenum simulated by Fluent transmits the pressure to the KMC-sub and NDK modules at the second interface. At the same time, the KMC-sub and NDK modules transmit the pressure to Fluent at the first interface and transmit the mass flow and temperature to Fluent at the second interface. CFD model 3 In this study, the entire primary loop model, except for the reactor core, was built in Fluent using an unstructured mesh obtained from ICEM CFD, as shown in Fig. 5(c) Fig. 5(c) . The number of meshes used for this simulation is approximately 3.86 million. In the CFD model, the inlet boundary was the reactor core outlet, where the mass flow rate and temperature were obtained from KMC-sub, and the pressure was transferred from Fluent to KMC-sub. The outlet boundary is the reactor core inlet and adopts a pressure outlet where the pressure is transferred from KMC-sub, and the mass flow and temperature are transferred from Fluent to KMC-sub. In addition, the outlet temperature and mass flow of the heat exchanger (see location A in Fig. 5(b) Fig. 5(b) ) provided by ATHLET is another inlet boundary condition of Fluent. At present, it has been widely used by scholars to simulate the flow and heat transfer of liquid metals based on the modified full flow model and wall function method. The k–ε model was used to calculate the thermal stratification in the cylindrical hot pool of LMFBR, which is a liquid metal fast reactor [ 21 21 ]. The standard k – ε model, k – ε (RNG) model, and Reynolds stress model were used to investigate thermal stratification in the sodium reactor and compared to experimental results, which showed that the RNG model has a milder temperature gradient [ 22 22 ]. The k – ε (RNG) turbulent model was used to calculate the primary cooling system of the lead-bismuth reactor, MYRRHA [ 23 23 ]. In this study, the k – ε (RNG) model and standard wall functions were used to calculate the turbulent flow inside the entire primary loop. Neutronic/thermal-hydraulic model 3 Because the KMC-sub and NDK modules are written based on C++, the classes of KMC-sub and NDK were constructed, respectively, and the two modules realize the transfer and exchange of data through API functions written in DLL. In the iterative calculation of a time step, the temperature of the coolant and fuel rods are calculated by the KMC-sub and transferred to the NDK module through the corresponding relationship of coordinates. The neutron flux density distribution is calculated by the NDK module and transferred to the power density distribution of the KMC-sub module, and the next time step calculation is carried out. The KMC-sub module and NDK module are coupled according to the mapping of the coordinate relation where the NDK code is built by the rectangular coordinate system, while the KMC-sub has a hexagonal mesh in every assembly. It is worth noting that NDK calculates the nuclear properties of materials in different areas according to whether the node coordinates are in the hexagonal assembly. As shown in Fig. 5(d) Fig. 5(d) , only the x - y plane is shown, and the z -direction (axis direction) uses the same method to map mesh data with the corresponding height. Compared with the mesh in KMC-sub representing an assembly, the mesh size of the NDK module is much smaller. Additionally, a KMC-sub mesh contains several NDK meshes. Each NDK mesh and KMC submesh have the coordinates of the center point. Before data mapping, it is determined whether the NDK mesh coordinates are in the KMC-sub control volume. If it is, data mapping will be carried out; otherwise, it will continue to search for an NDK mesh that meets the determination conditions in a KMC-sub mesh. When KMC-sub transmits temperature data to NDK, one KMC-sub mesh cell is mapped to multiple NDK mesh cells and transmits the same value to NDK; for power distribution transfer from NDK to KMC-sub, the average value of multiple meshes is adopted. Each time the mesh data mapping begins from the x-y plane first, then the z-direction is processed, and finally, the data transmission of different meshes is completed. In the core calculation, the temperature distribution of the coolant, cladding, and pellet is simulated in the KMC-sub calculation domain; when the temperature distribution is transferred to the NDK calculation domain, the NDK module updates the cross-section of the material according to the temperature, and then updates the core neutron flux density distribution, and then provides the corresponding power distribution to the KMC-sub calculation domain. In the process of time-step iteration, because of the delayed neutron effect, the time step of the NDK module is much smaller than that of KMC-sub. Therefore, it is necessary to perform multiple NDK module iterations in the time-step iteration of a KMC-sub. Analysis of shutdown transient of M 2 LFR-1000 and validation 1 Pre-simulation 2 A steady-state coupling in M 2 LFR-1000 is performed before the shutdown transient is implemented, in which Fluent first calculates independently to reach steady state, and then transfers coupling data with KMC-sub/NDK through UDF until all three codes reach a steady state. The heat source (reactor core) power shape and its value set into the first-step Fluent-alone calculation are obtained from the design values, which are calculated by the Monte Carlo code during the design stage. Some local and global steady-state physical and thermal-hydraulic parameters are shown in Fig. 6 Fig. 6 and Table 2 Table 2 , respectively. In addition, the design parameters are listed to verify the steady-state results obtained from the coupling program, and the value of the core pressure drop is obtained from the numerical simulation results in the preliminary design process. It can be seen from the results that the steady-state calculation results are in good agreement with the design values, laying the foundation for the following transient calculations. 1 1 Parameters 1 1 Predicted value 1 1 Design value 1 1 Deviation 1 1 Core outlet temperature 1 1 479.92 ℃ 1 1 480 ℃ 1 1 0.0017% 1 1 Total pressure drop (within core) 1 1 3.12×10 5 Pa 1 1 3.116×10 5 Pa 1 1 0.13% Result and discussion 2 The transient shutdown calculation of the coupled program was based on the results of the above steady-state calculation. Figure 7(a) Figure 7(a) shows the power change calculated by the coupling program after shutdown, and because the continuous insertion of control rods and the negative reactivity introduction reactor power drops quickly. When control rods arrive at the reactor bottom, the speed of power dropping becomes slow, and finally, the power reaches the level of decay heat. Figure 7(b) Figure 7(b) shows the average coolant temperature change at the core inlet and outlet calculated by the coupling program. It can be seen from Fig. 7(b) Fig. 7(b) that at 0 s the coolant temperature in the coupling program and ATHLET are 500.8 ℃ and 495.8 ℃, respectively. Then, as the reactor shutdown, the coolant temperature continues to decrease and finally stabilizes at 362.7 ℃ and 352.9 ℃ respectively at 250 s. At the same time, the highest temperature of cladding in the coupling program and ATHLET decrease from 530.6 ℃ and 526.1 ℃ to 363.1 ℃ and 354.6 ℃, respectively. The highest temperature of the pellet decreased from 1654.6 ℃ and 1551.5 ℃ to 442.5 ℃ and 387.0 ℃, respectively. Compared with the curve of power change, the temperature change has a certain lag, especially when the speed of the fuel pellet temperature decreases far exceeding that of the cladding and coolant because it takes time for energy to transfer in the rod bundles. Figure 7(c) and (d) Figure 7(c) and (d) show that the highest temperature of the coolant, pellet, and cladding calculated by the coupling program and ATHLET, respectively, change after shutdown. Comparing the calculation results between Fluent/KMC-sub/NDK and ATHLET during the shutdown process, the highest temperatures of pellet, cladding, and coolant calculated by Fluent/KMC-sub/NDK were all higher than those of ATHLET. This is mainly due to the lump-parameter method adopted in the ATHLET model, which has a certain number of pipes modeled as the reactor zone, and the maximum temperature is only the average value in this zone. However, all assemblies and sub-channels in the Fluent/KMC-sub/NDK model were modeled separately, so the maximum temperature obtained was more accurate. In general, the maximum errors of the pellet, cladding, and coolant in the initial and final states were 14.3%, 2.4%, and 2.8%, respectively, which demonstrates a good agreement with ATHLET, ensuring the correctness of the coupling program. Another advantage of the multi-scale and multi-physical coupling program is that it can capture the three-dimensional phenomena of the reactor core. The contours of temperature and velocity distribution after shutdown are shown in Fig. 8 Fig. 8 , from which we can see that after reactor shutdown, there exists an obvious thermal stratification phenomenon within the upper plenum region. In reactor pools of liquid metal cooled fast reactors (LMFRs), especially in the upper plenum, the magnitudes of the inertial and buoyancy forces are usually of similar order [ 24 24 ]. If the high-temperature liquid layer separates from the low-temperature liquid layer, the phenomenon called thermal stratification occurs. At the beginning of our simulation, thermal stratification occurred between the outer and inner regions of the upper plenum owing to the lack of coolant mixing in the upper core structure. In the early stage, the thermal stratification occurring between the outer and inner sides weakens, and the thermal stratification phenomena occurring between the top and bottom parts of the upper plenum begin to hold a dominant position. The Richardson number is introduced to distinguish the thermal stratification phenomenon in the reactor pool, which is used to judge the strength of the natural convection capacity and forced convection capacity [ 25 25 ]. g represents the acceleration of gravity, β is the thermal expansion coefficient of the lead coolant, D is the diameter of the upper plenum, Δ T is the coolant temperature change, and V is the velocity of the coolant. When the Richardson number is greater than one, thermal stratification occurs. <math display="block" id="M2"><mrow><mtext>R</mtext><mtext>I</mtext><mo>=</mo><mi>g</mi><mi>β</mi><mi>D</mi><mi>Δ</mi><mi>T</mi><mo>/</mo><msup><mi>V</mi><mn>2</mn></msup></mrow></math> As shown in Fig. 8(e-f) Fig. 8(e-f) , the coolant temperature at the reactor outlet decreases quickly, and because of the constant mass flow, the coolant with a lower temperature goes to the upper plenum at 10 s. At that time, the coolant temperature difference was 50 °C, the value of RI was 1.194, and the thermal stratification phenomena began to appear. Subsequently, the coolant temperature drops continuously, and at the same time, owing to the influence of buoyancy, the coolant with a higher temperature at the top of the upper plenum is gradually cooled by the heat exchanger, and then the reactor reaches a new steady state. At the same time, although the coolant mass flow remains constant, there are different velocity distributions owing to the coolant temperature at the outlet of the reactor decreasing continuously and the influence of buoyancy. In addition, from the perspective of engineering safety, this type of thermal stratification may cause the temperature boundary and stress effects of the surrounding structural boundaries to change, and further analysis of its influence is required by coupling the thermal-structural mechanics code. Compared to the point reactor neutron dynamics model in ATHLET, the diffusion neutron dynamics model can capture the spatial variation of the neutron flux density in the reactor and evaluate the neutronic and thermal characteristics of the reactor more accurately. Figure 9(a-d) Figure 9(a-d) show the contour of the neutron flux density changes with time at the horizontal section of the middle height of the core active zone. Initially, the control rod and safety rod were not inserted. During the shutdown process, all the control rods were inserted completely within 2 s, and because of the introduction of a large amount of negative reactivity during the control insertion, the neutron flux density of the reactor core decreased rapidly. In particular, around the control rod, the absorption of the control rod caused the nearby neutron flux density to tend to zero. When reaching a steady state, the distribution of the neutron flux density is relatively uniform. After the control rod is inserted, because of the role of the control rod, the neutron flux density presents a spatial change closely related to the position of the control rod. Further analysis of the change in neutron flux density with time and space shows that the overall neutron flux density decreases with time, but the neutron flux density in the core center is higher than that in other regions. The Fluent/KMC-sub/NDK coupling program can consider physical and thermal-hydraulic interactions. Both the overall and local parameters, that is, the coolant temperature, fuel temperature, pellet temperature, and other thermal parameters of each assembly can be captured. During the shutdown process, the changes in the coolant temperature at the core outlet are shown in Fig. 9(e-h) Fig. 9(e-h) . As time goes on, the temperature of the coolant at the outlet of the core decreases with time, which is the same as the change in the neutron flux density. Because the neutron flux density is relative to the reactor power, under the condition of constant inlet mass flow, the coolant outlet temperature gradually decreases. Further analysis of the spatial distribution of temperature shows that in the steady state, the temperature distribution is relatively uniform because the neutron flux density distribution of the reactor is also relatively uniform. With the insertion of control rods, the spatial distribution of power becomes higher at the center, which leads to the same distribution of coolant temperature. However, at 250 s, the minimum temperature at the outlet of the core coolant appeared in the external fuel area because the neutron flux density of the reactor was very small and the heat generated by fission tended to zero, and the energy of the coolant was mainly provided by the decay heat of the radioactive products in the reactor. At this time, the dominant factor of the coolant temperature changes from power to inlet flow, and because the flow rate of the center is relatively higher than that of the corner in the assembly, the temperature of the external fuel area is higher than that of the central fuel area. By analyzing the simulation results of the Fluent/KMC-sub/NDK coupling program, it was found that the numerical simulation results were consistent with the actual reactor shutdown process, which shows that the calculation of the coupled program is reasonable and accurate, and can be used for the thermal safety performance evaluation of the reactor. Conclusion 1 Based on the previously developed sub-channel code KMC-sub, neutron diffusion kinetic program NDK, a multi-physics and multiscale program Fluent/KMC-sub/NDK was developed using “Master-slave” mode. The coupling scheme is explicit order, and the reactor core domain is simulated by KMC-sub and NDK, in which the power and temperature distributions are transferred, while the rest of the domain of the primary loop is simulated by Fluent, where the temperature, mass flow, and pressure are transferred between KMC-sub and Fluent at the boundary. Based on the steady-state calculation, a shutdown process was carried out and compared with the ATHLET results for verification. The following conclusions can be drawn based on detailed analyses. (1) A multi-physics and multi-scale program, Fluent/KMC-sub/NDK, was developed using UDF and an explicit coupling method, and a shutdown process was simulated and compared with ATHLET, finally demonstrating good agreement. At the same time, the neutron flux density distribution in the reactor core changes with time is observed successfully. (2) In the process of shutdown, as the control rods are inserted within 2 s, the neutron flux density decreases rapidly, which causes the power to fall promptly. However, the maximum temperature of the coolant, fuel pellet, and cladding have a certain lag compared to power, especially the speed of fuel pellet temperature decrease far exceeding that of cladding and coolant. This proves that energy conversion is much faster than the transfer of energy between rod bundles. (3) During the shutdown process, thermal stratification occurs in the inner and outer parts of the upper plenum at the beginning, and when the RI number exceeds 1.0, the thermal stratification phenomena occurring between the upper and lower parts of the upper plenum begin to hold a dominant position. Finally, the thermal stratification phenomena generally disappear, but there are some mild thermal stratification phenomena between the inner and outer parts of the upper plenum due to the lack of mixture in the upper core structure. Future work will focus on two aspects: (1) thorough and deep research on the interplay between time steps, spatial mesh size, and spatial stencils on the convergence and precision of the multi-scale and multi-physics coupling simulation, and (2) expanding the capability of the multiscale and multiphysics coupling program by implementing more sub-tools (i.e., secondary simulation programs, thermomechanical codes, etc.) References D.L. Aumiller , E.T. Tomlinson , R.C. Bauer , A coupled RELAP5-3D/CFD methodology with a proof-of-principle calculation . Nucl. Eng. Des. 205 ( 1-2 ), 83 - 90 ( 2001 ). doi: 10.1016/S0029-5493(00)00370-8 http://doi.org/10.1016/S0029-5493(00)00370-8 A coupled RELAP5-3D/CFD methodology with a proof-of-principle calculation D.L. Aumiller , E.T. Tomlinson , W.L. Weaver , An integrated relap5-3d and multiphase cfd code system utilizing a semi-implicit coupling technique . Nucl. Eng. Des. 216 ( 1-3 ), 77 - 87 ( 2002 ). doi: 10.1016/S0029-5493(01)00522-2 http://doi.org/10.1016/S0029-5493(01)00522-2 An integrated relap5-3d and multiphase cfd code system utilizing a semi-implicit coupling technique A. Papukchiev , G. Lerchl , C. Waata et al. , Extension of the simulation capabilities of the 1D system code ATHLET by coupling with the 3D CFD software package ANSYS CFX . Paper presented at the proceedings of the 13th international topical meeting on nuclear reactor thermal-hydraulics (NURETH-13) , Kanazawa City, Japan , September 27-October 2, 2009 . Extension of the simulation capabilities of the 1D system code ATHLET by coupling with the 3D CFD software package ANSYS CFX A. Papukchiev , M. Jeltsov , K. Kööp et al. , Comparison of different coupling CFD–STH approaches for pre-test analysis of a TALL-3D experiment . Nucl. Eng. Des. 290 , 135 - 143 ( 2015 ). doi: 10.1016/j.nucengdes.2014.11.008 http://doi.org/10.1016/j.nucengdes.2014.11.008 Comparison of different coupling CFD–STH approaches for pre-test analysis of a TALL-3D experiment R. Baviere , N. Tauveron , F. Perdu et al. , A first system/CFD coupled simulation of a complete nuclear reactor transient using CATHARE2 and TRIO_U. Preliminary validation on the Phénix Reactor Natural Circulation Test . Nucl. Eng. Des. 277 , 124 - 137 ( 2014 ). doi: 10.1016/j.nucengdes.2014.05.031 http://doi.org/10.1016/j.nucengdes.2014.05.031 A first system/CFD coupled simulation of a complete nuclear reactor transient using CATHARE2 and TRIO_U. Preliminary validation on the Phénix Reactor Natural Circulation Test Z. R. Zou , C. Shen , X.L. Zhang et al. , 3D thermal hydraulic characteristics analysis of pool-type upper plenum for lead-cooled fast reactor with multi-scale coupling program . Nucl. Eng. Des. 370 ( 2020 ): 110892 . 3D thermal hydraulic characteristics analysis of pool-type upper plenum for lead-cooled fast reactor with multi-scale coupling program C. Wang , L.K. Cao , X.L. Zhang et al. , Development and application of CFD and subchannel coupling analysis code for lead‐cooled fast reactor . Int. J. Energy Res. 43 ( 14 ), 8447 - 8462 ( 2019 ). doi: 10.1002/er.4845 http://doi.org/10.1002/er.4845 Development and application of CFD and subchannel coupling analysis code for lead‐cooled fast reactor X.B. Zhang , Z. Qin , H.L. Chen , Development and validation of a coupled neutron diffusion-thermal-hydraulic calculation procedure for fast reactor applications . Ann. Nucl. Energy. 139 , 107234 ( 2020 ). doi: 10.1016/j.anucene.2019.107243 http://doi.org/10.1016/j.anucene.2019.107243 Development and validation of a coupled neutron diffusion-thermal-hydraulic calculation procedure for fast reactor applications K. Wang , S.C. Liu , Z.G. Li et al. , Analysis of BEAVRS two-cycle benchmark using RMC based on full core detailed model . Prog. Nucl. Energ. 98 , 301 - 312 ( 2017 ). doi: 10.1016/j.pnucene.2017.04.009 http://doi.org/10.1016/j.pnucene.2017.04.009 Analysis of BEAVRS two-cycle benchmark using RMC based on full core detailed model S. C. Liu , J. G. Liang , Q. Wu et al. , BEAVRS full core burnup calculation in hot full power condition by RMC code . Ann. Nucl. Energy. 101 ( 2017 ): 434 - 446 . doi: 10.1016/j.pnucene.2017.04.009 http://doi.org/10.1016/j.pnucene.2017.04.009 BEAVRS full core burnup calculation in hot full power condition by RMC code C. Fiorina , I. Clifford , M. Aufiero et al. , GeN-Foam: a novel OpenFOAM® based multi-physics solver for 2D/3D transient analysis of nuclear reactors . Nucl. Eng. Des. 294 , 24 - 37 ( 2015 ). doi: 10.1016/j.nucengdes.2015.05.035 http://doi.org/10.1016/j.nucengdes.2015.05.035 GeN-Foam: a novel OpenFOAM® based multi-physics solver for 2D/3D transient analysis of nuclear reactors C. Fiorina , N. Kerkar , K. Mikityuk et al. , Development and verification of the neutron diffusion solver for the GeN-Foam multi-physics platform . Ann. Nucl. Energy 96 , 212 - 222 ( 2016 ). doi: 10.1016/j.anucene.2016.05.023 http://doi.org/10.1016/j.anucene.2016.05.023 Development and verification of the neutron diffusion solver for the GeN-Foam multi-physics platform M.R. Altahhan , S. Bhaskar , D. Ziyad et al. , Preliminary design and analysis of Liquid Fuel Molten Salt Reactor using multi-physics code GeN-Foam . Nucl. Eng. Des. 369 , 110826 ( 2020 ). doi: 10.1016/j.nucengdes.2020.110826 http://doi.org/10.1016/j.nucengdes.2020.110826 Preliminary design and analysis of Liquid Fuel Molten Salt Reactor using multi-physics code GeN-Foam Z. Chen , X.N. Chen , A. Rineiski et al. , Coupling a CFD code with neutron kinetics and pin thermal models for nuclear reactor safety analyses . Ann. Nucl. Energy 83 , 41 - 49 ( 2015 ). doi: 10.1016/j.anucene.2015.03.023 http://doi.org/10.1016/j.anucene.2015.03.023 Coupling a CFD code with neutron kinetics and pin thermal models for nuclear reactor safety analyses A.R. ANSYS® , Release 15.0, ANSYS FLUENT UDF Manual , ANSYS, Inc , 2013 . Release 15.0, ANSYS FLUENT UDF Manual S.Z. Li , L.K. Cao , M.S. Khan et al. , Development of a sub-channel thermal-hydraulic analysis code and its application to lead cooled fast reactor . Appl. Therm. Eng. 117 , 443 - 451 ( 2017 ). doi: 10.1016/j.applthermaleng.2017.02.044 http://doi.org/10.1016/j.applthermaleng.2017.02.044 Development of a sub-channel thermal-hydraulic analysis code and its application to lead cooled fast reactor H.L. Chen , X.L. Zhang , Y.S. Zhao et al. , Preliminary design of a medium‐power modular lead‐cooled fast reactor with the application of optimization methods . Int. J. Energy Res. 42 , 3643 - 3657 ( 2018 ). doi: 10.1002/er.4112 http://doi.org/10.1002/er.4112 Preliminary design of a medium‐power modular lead‐cooled fast reactor with the application of optimization methods C. Shen , X.L. Zhang , C. Wang et al. , Transient safety analysis of M2LFR-1000 reactor using ATHLET . Nucl. Eng. Technol. 51 , 116 - 124 ( 2019 ). doi: 10.1016/j.net.2018.08.011 http://doi.org/10.1016/j.net.2018.08.011 Transient safety analysis of M2LFR-1000 reactor using ATHLET X. Luo , X.L. Zhang , S. Wang et al. , Automated core design code development for a lead‐cooled fast reactor and its core optimization . Int. J. Energy Res. 45 , 11721 - 11734 ( 2021 ). doi: 10.1002/er.5553 http://doi.org/10.1002/er.5553 Automated core design code development for a lead‐cooled fast reactor and its core optimization A. Wielenberg , C. Bals , Nuclear thermal hydraulics with the AC² system code package . Nuclear Power Plant Design and Analysis Codes. 12 , 277 - 311 ( 2021 ). doi: 10.1016/B978-0-12-818190-4.00012-7 http://doi.org/10.1016/B978-0-12-818190-4.00012-7 Nuclear thermal hydraulics with the AC² system code package S.D. Shyam , K.S. Akshaya , G. Padmakumar et al. , CFD analysis of thermal stratification and sensitivity study of model parameters for k– model in a cylindrical hot plenum . Nucl. Eng. Des. 250 , 417 - 435 ( 2012 ). doi: 10.1016/j.nucengdes.2012.04.008 http://doi.org/10.1016/j.nucengdes.2012.04.008 CFD analysis of thermal stratification and sensitivity study of model parameters for k– model in a cylindrical hot plenum S. Ohno , H. Ohshima , H. Ohki et al. , Validation of numerical simulation method for thermal stratification in reactor vessel upper plenum of fast reactor . Paper presented at the proceeding of 6 th Japan-Korea symposium on nuclear thermal-hydraulics and safety(NTHAS6), Nov. 24-27 , Okinawa, Japan , 2008 . Validation of numerical simulation method for thermal stratification in reactor vessel upper plenum of fast reactor M. Vanderhaegen , J. Vierendeels , B. Arien , CFD analysis of the MYRRHA primary cooling system . Nucl. Eng. Des. 241 ( 3 ), 775 - 784 ( 2011 ). doi: 10.1016/j.nucengdes.2010.12.009 http://doi.org/10.1016/j.nucengdes.2010.12.009 CFD analysis of the MYRRHA primary cooling system F. Roelofs . Thermal hydraulics aspects of liquid metal cooled nuclear reactors . M. Woodhead Publishing , 2018 Thermal hydraulics aspects of liquid metal cooled nuclear reactors B. Cai , H.F. Gu , Y. Weng , et al . Numerical investigation on the thermal stratification in a pressurizer surge line . Ann. Nucl. Energy 101 , 293 - 300 ( 2017 ). doi: 10.1016/j.anucene.2016.11.024 http://doi.org/10.1016/j.anucene.2016.11.024 Numerical investigation on the thermal stratification in a pressurizer surge line All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xiao Luo, Chi Wang, Ze-Ren Zou, Lian-Kai Cao, Shuai Wang, Zhao Chen, 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. We thank GRS for providing the ATHLET code. 10.1007/s41365-022-01008-y This work was supported by Science and Technology on Reactor System Design Technology Laboratory, Chengdu, China (No. LRSDT2020106). 2021-07-27 2021-12-21 2022-01-13 2022-02-25 2022-03-18 © 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-01008-y http://dx.doi.org10.1007/s41365-022-01008-y . 2022 Nuclear Science and Techniques 02 33