Development and verification of CFD model of carbothermal synthesis furnace for production of mixed nitride uranium-plutonium fuel✩

Rinat N. Shamsutdinov; Sergey V. Pavlov; Anton Yu. Leshchenko

doi:10.3897/nucet.6.51253

Development and verification of CFD model of carbothermal synthesis furnace for production of mixed nitride uranium-plutonium fuel^✩

Rinat N. Shamsutdinov^‡§, Sergey V. Pavlov^‡§, Anton Yu. Leshchenko^‡

‡ Sosny R&D Company, Dimitrovgrad, Russia

§ DETI MEPhI, Dimitrovgrad, Russia

Corresponding author: Rinat N. Shamsutdinov ( srn@sosny.ru )

Academic editor: Georgy Tikhomirov

This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Citation: Shamsutdinov RN, Pavlov SV, Leshchenko AYu (2020) Development and verification of CFD model of carbothermal synthesis furnace for production of mixed nitride uranium-plutonium fuel. Nuclear Energy and Technology 6(1): 37-42. https://doi.org/10.3897/nucet.6.51253

Abstract

The retort batch furnace for the carbothermal synthesis of uranium and plutonium nitrides is a component of the mixed nitride uranium-plutonium (MNUP) fuel Fabrication/Refabrication Module (FRM), a part of the Pilot Demonstration Energy Complex (PDEC). A CFD model of the furnace was built in the SolidWorks Flow Simulation code to check the feasibility of its thermophysical operating modes (heating and cooling rates, temperature of the product loaded into the furnace). The experimental data obtained from the performance tests was used to verify the developed CFD model and to confirm its adequacy. The relative deviation of the calculated temperature of the loaded product from the experimental data in the process of isothermal annealing does not exceed 0.7%.

The temperatures of the loaded product predicted by the CFD model were used to justify engineering solutions for the uranium and plutonium nitrides carbothermal synthesis furnace. The CFD model can be used to define the furnace operation mode by selecting the gas flow rates inside and outside the retort, the heater temperature, and heating and cooling rates for the product loaded in the furnace.

Keywords

CFD model, MNUP fuel, SolidWorks Flow Simulationv, furnace, carbothermal synthesis, mathematical model, grid model

Introduction

For the past three decades, there has been an increasing tendency to use the Computational Fluid Dynamics (CFD) methods for computational justifications of engineering solutions in designing various industrial facilities (Spalart and Venkatakrishnan 2016, Raynal et al. 2016, Mizzi et al. 2015, Joshi and Nayak 2019). Using modern CFD codes enables modeling heat and mass transport, hydro- and gas dynamics hard to reproduce experimentally. The results of the CFD simulation modeling allow choosing acceptable engineering and process parameters for the equipment (liquid/gas flow rate and concentration, temperature, pressure, etc.) that may improve the reliability and simplify the design of the facility modeled (Besagni et al. 2015, Wang and Li 2015, Thompson et al. 2017).

Unlike other industries, the nuclear industry took the report by the joint technical meeting of the International Atomic Energy Agency (IAEA) and the Nuclear Energy Agency (NEA) held in Pisa, Italy, on 11 through 14 November 2002 as the starting point for a more extensive application of the CFD codes for computational studies and safety justifications of technologies (IAEA-TECDOC-1379 20028). Since that time, the IAEA has been engaged in comprehensive activities to encourage the use of the CFD technique in the nuclear industry, with the CFD code verification making part of it (Mahaffy 2015, Smith 2015, Bestion 2014). The primary goal of the verification is to confirm adequacy of the computational model and its numerical implementation by comparing the modeling results against experimental data. As shown by the latest tendency in the CFD code development in Russia and noted by many experts, the verification of the CFD codes to the extent required for application in the nuclear industry has not yet been completed (Bykov and Kurnosov 2018). Primarily, this due to the lack of experimental data for a broad range of CFD problems. In this connection, a topical task is to confirm the applicability of the CFD codes for verifying the engineering solutions for development of nuclear power equipment.

As part of the retort batch furnace design for the carbothermal synthesis of uranium and plutonium nitrides, requirements were defined for the thermophysical modes of the furnace operation, such as heating and cooling rates, the temperature of the product loaded into the furnace, flow rates of process gases overpressure, and others. To check the feasibility of the thermophysical modes set for the retort batch furnace operation, a need arose for building a computational thermophysical model. A CFD code, SolidWorks Flow Simulation (Alyamovsky et al. 2005), was used to solve the task. The purpose of the study is to build a CFD furnace model to calculate the temperature distribution inside the furnace retort for the given values of the heater temperatures and the flow rates of the process gases fed into the retort and into the space behind it. The paper presents a description of the CFD model designed to simulate the thermophysical operating modes of the furnace for the carbothermal synthesis of uranium and plutonium nitrides. A comparison is provided between the calculation results and the experimental data for the temperature distribution on the product loaded into the furnace in the working mode of operation.

Uranium and plutonium nitrides carbothermal synthesis furnace

The retort batch furnace for the carbothermal synthesis of uranium and plutonium nitrides, a 3D model of which is shown in Fig. 1, is a component of the equipment chain in a mixed uranium-plutonium nitride (MNUP) fuel fabrication/refabrication module (FRM), a part of the Pilot Demonstration Energy Complex (PDEC).

Figure 1.

3D model of the carbothermal synthesis furnace: 1, 6 – inlet (outlet) gate valves; 2 – furnace body water cooling circuit; 3 – furnace gas discharge pipe; 4 – electric heaters; 5 – furnace body; 7 – retort; 8 – boats with loaded product; 9 – lining.

The furnace is designed to produce uranium and plutonium nitrides via carbon-based carbothermal reduction of uranium and plutonium dioxides in nitrogen and nitrogen-hydrogen atmosphere at temperature of 1650 ± 50 °С. The retort batch furnace design and the carbothermal synthesis process are described in (Smirnov et al. 2019). Compacted powders of uranium and plutonium dioxides mixed with carbon and a binding agent (zinc stearate) are placed in rectangular open top containers (boats) inside the furnace retort. The furnace is equipped with a three-compartment heater with individual current control in each compartment, this providing for the uniform heating throughout the retort length. For a simpler and more reliable furnace design, the temperature of the loaded product in the carbothermal synthesis process is determined indirectly, that is, from the heater temperature.

Physical statement of the problem and a CFD model of the uranium and plutonium nitrides carbothermal synthesis furnace

Physical problem statement

In a retort batch furnace of the height H, the length L, and the width S, nitrogen with the temperature T_ni is fed into the retort space, and argon with the temperature T_ar is fed into the space behind the retort with the respective flow rates Q_r and Q_b. In the process of the heat transfer from the heaters installed inside the space behind the retort, the simulators of product loaded into the boats, which are placed inside the furnace retort, are heated to the given temperature Т_h followed by their isothermal annealing and cooling. A lining of the thickness H_l, on the outer surface of which there is a system of water cooling tubes with the overall flow rate Q_w, is used as the thermal insulation for the furnace hearth. Key assumptions: liquid and gases are considered to be single-phase incompressible Newtonian continuums, and the inner and outer walls of the wetted parts are considered to be hydraulically smooth surfaces. We need to find the temperature distribution within the furnace during non-steady-state and steady-state modes of operation.

The SolidWorks Flow Simulation code was selected to solve the problem at hand. In a broad range of commercial software products for the CFD problem solving, the SolidWorks Flow Simulation code has a number of advantages the most prominent of which are

application of a rectangular structured computational grid (Sobachkin and Dumnov 2014);
capability to simulate complex geometries of the calculated assemblies;
user friendliness and a self-explanatory interface.

The developed CFD model of the furnace consists of a mathematical model and a grid model.

Mathematical model

Reynolds-averaged Navier-Stokes equations are used in the SolidWorks Flow Simulation as the mathematical model’s base equations to describe the motion and heat exchange in liquid and gas in a 3D case. The mass, impulse and energy conservation laws for these fluids (Alyamovsky et al. 2005) are

$\frac{\partial ρ}{\partial t} + \frac{\partial (ρ u_{i})}{\partial x_{i}} = 0$

$\frac{\partial (ρ u_{i})}{\partial t} + \frac{\partial (ρ u_{i} u_{j})}{\partial x_{j}} + \frac{\partial P}{\partial x_{i}} = \frac{\partial (τ_{i j} + τ_{i j}^{R})}{\partial x_{j}} + S_{i}$

$\frac{\partial ρ H}{\partial t} + \frac{\partial ρ u_{i} H}{\partial x_{i}} = \frac{\partial (u_{j} (τ_{i j} + τ_{i j}^{R}) + q_{i})}{\partial x_{i}} + \frac{\partial p}{\partial t} - τ_{i j}^{R} \frac{\partial u_{i}}{\partial x_{j}} + ρ ε + S_{i} u_{i} + Q_{H}$

H = h + u²/2,

where u is the fluid velocity; ρ is the fluid density; P is the fluid pressure; S_i is the external mass forces acting on the fluid unit mass; H is the total energy of the fluid unit mass; Q_H is the heat emitted by the source in the fluid unit volume; τ_ij is the tensor of elastic shear stresses; and q_i is the diffusion heat flux.

These equations are supplemented with a gas or liquid state equation and with dependences of density, viscosity, heat capacity and thermal conductivity coefficient on temperature. A two-parameter turbulence model, k–ε (Wilcox 1998), is used to determine the turbulent viscosity and kinetic energy of turbulent pulsations which make a part of the equation for the viscous shear stress tensor.

Heat transfer in solid-state bodies is simulated using the equation

$\frac{\partial ρ c T}{\partial t} = \frac{\partial}{\partial x_{i}} (λ_{i} \frac{\partial T}{\partial x_{i}}) + Q_{H}$

where c is the specific heat capacity; T is the temperature; λ is the thermal conductivity coefficient; and Q_H is the specific (volumetric) heat generation by the heat source.

Heat exchange by radiation is simulated in accordance with the Stefan-Boltzman law:

Q_H = εσ₀(T_w⁴ – T_s⁴),

where ε is the surface emissivity factor; σ₀ is the Stefan-Boltzman constant; T_w is the radiation-emitting surface temperature; and T_s is the environment temperature. The thermophysical properties of the furnace materials, gas and water (density, heat capacity, thermal conductivity coefficient, surface emissivity factor) were defined according to references (Chirkin 1968). The properties of the furnace lining and heater materials were defined based on the manufacturer’s technical documentation.

In the case of coupled heat transfer, into account is taken of the heat flux between the solid-state body and the liquid (gas), the solid-state body wall temperature, the near-wall layer characteristics, and the heat exchange by radiation (where required).

The initial and boundary conditions are defined to close the system of differential equations. The furnace material, gas and water temperature: T (x, y, z, 0) = 20 °C, and the gas and water pressure and velocity: P (x, y, z, 0) = 101325 Pa and V (x, y, z, 0) = 0 m/s, where x, y, z are their coordinates, were preset as the boundary conditions.

The Derichlet’s boundary conditions: T(x_h, y_h, z_h, t) = T_i = const, where x_h, y_h, z_h are the heater surface coordinates, and T_i is the value of the heater temperature et each time i, are defined for the heater surface (Fig. 2). The Derichlet’s boundary conditions: V (x_w, y_w, z_w, t) = const and T (x_w, y_w, z_w, t) = const, where x_w, y_w, z_w are the coordinates of the gas supply and water cooling surfaces, are defined for the surfaces of the nozzles for the gas supply into the retort space and into the space behind the retort, and for the surfaces of the water cooling tubes. The Newmann’s boundary condition are defined for the outlet of the tubes for the gas supply into the retort area and into the space behind the retort and for the outlet of the water cooling tubes

∂V (x_out, y_out, z_out, t) / ∂n = 0,

∂T (x_out, y_out, z_out, t) / ∂n = 0,

∂P (x_out, y_out,z_out, t) / ∂n = 0,

where n is the vector of the normal to the boundary; and x_out, y_out, z_вout are the coordinates of the outlet tubes. The Derichlet’s boundary conditions for the pressure P (x_b, y_b, z_b, t) = 101325 Pa, the temperature T (x_b, y_b, z_b, t) = 20 °C, and the air velocity V (x_b, y_b, z_b, t) = 0 m/s, where x_b, y_b, z_b are the coordinates of the computational domain boundary, and t is the time, are defined for the computational domain boundary (see Fig. 2), except the boundary formed by the symmetry plane perpendicular to the longitudinal axis of the furnace retort. The symmetry boundary means the equality to zero of the gas velocity component normal to the boundary, and of the gradients for all other values (Φ – temperature, pressure):

V (x_s, y_s, z_s, t) = 0 m/s, ∂Φ / ∂n = 0,

where n is the vector of the normal to the symmetry boundary; and x_s, y_s, z_s are the coordinates of the symmetry boundary.

Finite volume method is used to discretize the resultant system of differential equations in the SolidWorks Flow Simulation code. The test item is presented as a grid model the cells in which are parallelepiped-shaped (the values of independent variables are calculated at the cell centers, and the mass, impulse, and energy flows are calculated on the cell faces). Linearized equations describing the conservation law for the investigated scalar physical quantity are written for each cell. Spatial derivatives are approximated using implicit difference operators of the second-order accuracy (modified implicit Leonard QUICK-approximation (Roache 1998) and TVD total variation minimization method (Hirsch 1991)). With a time discretization (SIMPLE-type method (Patankar 1980)), the maximum allowable time step, depending on both the values of physical quantities and the discretization step over the space in that cell, is determined for each grid cell using the Courant condition. Generalized conjugate gradients method (Saad 1996) is used to solve the system of fully implicit linear discrete equations involving incomplete LU-factorization method (Alyamovsky et al. 2005). The convergence of the solution was estimated from the nature of the changes in the values of the target parameters (temperature on the loaded product surface, balance of the nitrogen and cooling water flow rates in the retort area and in the space behind the retort) with a fixed time step. The target parameter is considered to be converged if the calculated value of the target parameter fluctuation amplitude range (delta) becomes smaller than the established convergence criterion (for the temperature of 1 °C, for the nitrogen and cooling water flow rates of 0.000001 kg/h).

Grid model

To save the computational resources, a furnace half (see Fig. 2) was used as the computational domain for building the grid model. The computational domain represents a parallelepiped symmetrical to the parallelepiped’s median plane (passes through the central boat). The symmetry plane is perpendicular to the longitudinal axis of the furnace retort.

The grid convergence was studied to estimate the dependence of the calculated temperature values for the furnace’s local regions on the size of the grid model. Fig. 3 presents the dependence of the calculated steady-state temperature on the loaded product surface in the boat at one of the extreme positions on the computational grid size. A computational grid with the size of 7.65 million cells was chosen as the result of the grid convergence analysis. The explanation for the choice is that the relative deviation of the temperature on the loaded product surface is 0.17% for the grid dimensionality of 5.36 million to 3.68 million, and 0.14% for the dimensionality of 7.65 million to 5.36 million. These deviations are 2.7 and 2.3 °С in the absolute temperature values. Further fragmentation and an increase in the grid dimensionality will lead to a smaller deviation. For instance, the value of the experimental error for this temperature range is ± 8 °C. Therefore, no increase in the grid size is reasonable, including for the sake of saving the computational resources.

Local fragmentation of the grid cells in the solid-state body and fluid contact region is used to resolve relatively small geometrical features of the grid model. Regions with temperature gradients in solid-state bodies of the computational model (boats with the loaded product, the retort, nitrogen supply tubes, etc.) are solved in the same way.

Figure 2.

Computation region of the CFD furnace model.

Figure 3.

Grid convergence investigation.

CFD furnace model verification results

The results of temperature testing the retort batch furnace for the carbothermal synthesis of uranium and plutonium nitrides were used to check the adequacy of the developed CFD furnace model. The purpose of the experimental study is to estimate the temperature distribution on the loaded product in the boats during different temperature modes of the furnace operation (heating, annealing, and cooling). Five tungsten boats were placed inside the furnace retort, into which simulators of product (silica brick) with similar thermophysical properties were loaded instead of compacted uranium and plutonium dioxide powders. A mode close to the furnace commercial operation conditions was chosen for the experiment. Temperature values on the heaters and the loaded product, values of the nitrogen flow rates fed into the retort and into the space behind the retort, and the water flow rate in the furnace body water cooling circuit were recorded during the test. The loaded product temperatures were measured using three tungsten-rhenium thermocouples with a temperature measurement accuracy of ± 5 °C in a range of 0 to 1000 °C and ± 0.005×t for the range of 1000 to 1950 °C (t is the measured temperature value). The working ends of the two thermocouples (hot junctions) were contained inside the boats at the extreme positions (boat 1 and boat 5), resting on the loaded product surface and at a distance of ~15 mm from the boat outer wall – area 1 and area 3 respectively (Fig. 4). Thermocouple 3 was placed on the loaded product surface (area 2) at the center of boat 3.

Figure 4.

Boats 1 through 5 shown inside the furnace retort.

The dependences of the loaded product temperature in the furnace boats at the thermocouple locations on time (Figs 5, 6) were determined based on the experiment and calculation results. The values of the loaded product temperature recorded during the experiment in the boats at the extreme positions (1 and 5) differ insignificantly, this indicating that the resultant temperature profile is symmetrical (see Fig. 5).

Figure 5.

Temperature of the loaded product in the boats at the extreme positions inside the retort (areas 1 and 3).

Figure 6.

Temperature of the loaded product in the central boat inside the retort (area 2).

It can be seen in Figs 5 and 6 that the calculated and experimental values of the loaded product temperature in the boat at an extreme position and in the central boat agree well both during the isothermal annealing, heating and cooling. There is also a major deviation observed between the calculated and experimental temperature values during the loaded product heating and cooling in the boats at the extreme positions as compared with the central boat.

The maximum relative deviation of the calculated average rate of the loaded product heating in the central boat and in the boats at the extreme positions from the experimental value for the time span of 0.15 to 0.33 rel. units is 7.6 and 3.6% respectively.

The maximum relative deviation of the calculated average rate of the loaded product cooling in the central boat and in the boats at the extreme positions from the experimental value after 0.45 rel. time units is 10 and 6% respectively.

Fig. 7 shows a screenshot with the calculated and experimental values of the loaded product surface temperature in the boats in the isothermal annealing process (a steady-state mode). The relative deviation of the calculated temperature values from the experimental values is

0.7% in the boats at the extreme positions;
and 0.1% in the central boat.

The obtained values of the relative deviation between the calculated and experimental data indicate that the accuracy of the calculations based on a CFD furnace model is sufficiently high.

Figure 7.

Loaded product surface temperature profile in the isothermal annealing process: –––– – calculation; ♦ – experiment (a screenshot).

Conclusion

1. A CFD model of a furnace for the carbothermal synthesis of uranium and plutonium nitrides has been developed in the SolidWorks Flow Simulation code, which makes it possible to simulate thermophysical modes of the furnace operation (heating, annealing and cooling of the loaded product) without regard for the physicochemical transformations taking place in the loaded fuel carbothermal synthesis process.

2. The developed CFD furnace model has been verified and its adequacy has been confirmed based on the experimental temperature test data for the uranium and plutonium nitrides carbothermal synthesis furnace. The relative deviation between the values of the calculated loaded product temperature in the boats at the extreme positions and in the central boat and the experimental data in the isothermal annealing process is 0.7% and 0.1% respectively. The estimated values of the loaded product temperature obtained using the developed CFD furnace model were used to justify the designs adopted for the uranium and plutonium nitrides carbothermal synthesis furnace. These results can be also used to simulate the kinetics of the MNUP fuel carbothermal synthesis reaction in the nitrogen (hydrogen-nitrogen mixture) flow in this furnace.

3. The developed CFD model of the uranium and plutonium nitrides carbothermal synthesis furnace can be used to choose the furnace operating mode by selecting among the values of the process gas flow rates fed into the furnace’s retort area and into the space behind the retort, the heater temperature, and the loaded product heating and cooling rates. The CFD model was used in the process of developing and justifying the serviceability of the commercial furnace design. The results obtained using the CFD model are a part of the detailed design documentation (DDD) and the technical reports for the considered furnace.

References

Alyamovsky AA, Sobachkin AA, Odintsov YeV, Kharitonovich AI, Ponomarev NB (2005) SolidWorks. Computer Simulation in Engineering Practice. Saint Petersburg. BKhV-Petersburg Publ., 800 pp. [in Russian]

Besagni G, Mereu R, Inzoli F (2015) CFD study of ejector flow behavior in a blast furnace gas galvanizing plant. Journal of Thermal Science 24(1): 58–66. https://doi.org/10.1007/s11630-015-0756-4

Bestion D (2014) Extension of CFD Codes Application to Two-Phase Flow Safety Problems. NEA-CSNI-R-2014-13, Nuclear Energy Agency of the OECD (NEA).

Bykov MA, Kurnosov MM (2018) Designs of JSC OKB “Gidropress” for the Implementation of the Research and Development Investment Project “Introduction of CFD Methods to Support the Design and Justification of Projects». Tasks, Composition and Progress of Implementation. Proc. of the Scientific and Technical Workshop “The Problem of Verification and Application of CFD Codes in Nuclear Power Engineering”. Nizhniy Novgorod. OAO “OKBM Afrikantova” Publ., 5–16. [in Russian]

Chirkin VS (1968) Thermophysical Properties of Nuclear Equipment Materials. Moscow. Atomizdat Publ., 485 pp. https://doi.org/10.1016/0029-5493(68)90086-1 [in Russian]

Hirsch C (1991) Numerical Computation of Internal and External Flows. Vol. 1. Chichester. John Wiley&Sons Ltd, 538 pp.

IAEA-TECDOC-1379 (2002) Use of computational fluid dynamics codes for safety analysis of nuclear reactor systems: 1–50.

Joshi J, Nayak A (2019) Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment. Duxford, United Kingdom: Woodhead Publ., 888 pp.

Mahaffy J (2015) Best Practice Guidelines for the Use of CFD in Nuclear Reactor Safety Applications. NEA-CSNI-R-2014-11, Nuclear Energy Agency of the OECD (NEA).

Mizzi K, Kellett P, Demirel YK, Martin R, Turan O (2015) HPC and CFD in the marine industry: past, present and future. Proc. of the III^rd Intern. Conf. on Exascale Applications and Software. Edinburgh: 92–97.

Patankar SV (1980) Numerical Heat Transfer and Fluid Flow. Washington. Hemisphere, 197 pp.

Raynal L, Augier F, Bazer-Bachi F, Haroun Y, Pereira da Fonte C (2016) CFD applied to process development in the oil and gas industry. Oil & Gas Science and Technology 71(3): 1–24. https://doi.org/10.2516/ogst/2015019

Roache PJ (1998) Fundamentals of Computational Fluid Dynamics. Albuquerque, New Mexico. Hermosa Publ., 648 pp.

Saad Y (1996) Iterative Methods for Sparse Linear Systems. Boston. PWS Publ., 184 pp.

Smirnov VP, Pavlov SV, Ivanov DV, Shamsutdinov RN, Inyutin NV, Sidorenko DM, Solopeko AV, Sayfutdinov SYu, Davydov AV, Chamovskikh YuV, Sergeyev NG, Zozulya DV (2019) A Carbothermal Facility for the Synthesis of Nitrides for Mixed Uranium-Plutonium Nitride Fuel Production. Atomnaya energiya 125(5): 322–325. https://doi.org/10.1007/s10512-019-00488-3 [in Russian]

Smith BL (2015) Assessment of CFD Codes for Nuclear Reactor Safety Problems. NEA-CSNI-R-2014-12, Nuclear Energy Agency of the OECD (NEA).

Sobachkin A, Dumnov G (2014) Numerical Basis of CAD-Embedded CFD. Dassault systems: 1–19.

Spalart PR, Venkatakrishnan V (2016) On the role and challenges of CFD in the aerospace industry. The Aeronautical Journal, 120(1223): 209–232. https://doi.org/10.1017/aer.2015.10

Thompson JS, Walton S, Hassan S, Rolland SA, Sienz J (2017) The use of CFD and multi-objective optimization techniques to customise an industrial pre-mixer. Structural and Multisiciplinary Optimization 55(6): 2339–2351. https://doi.org/10.1007/s00158-016-1643-7

Wang X, Li Y (2015) A comprehensive 3D mathematical model of the electroslag remelting process. Metallurgical and Materials Transactions B 46(4): 1837–1849. https://doi.org/10.1007/s11663-015-0342-x

Wilcox DC (1998) Turbulence Modeling for CFD. La Canada, California: DCW Publ., 460.

* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2019, n. 4, pp. 130–141.