Modeling and Simulation of Flow and Uranium Isotopes Separation in Gas Centrifuges Using Implicit Coupled Density-Based Solver in OpenFOAM

Valiyollah Ghazanfari¹, Ali Akbar Salehi², Ali Reza Keshtkar^1,*, Mohammad Mahdi Shadman³ and Mohammad Hossein Askari³

¹Materials and Nuclear Fuel Research School, Nuclear Science and Technology Research Institute, AEOI, P.O. Box: 11365-8486, Tehran-Iran

²Department of Energy Engineering, Sharif University of Technology, P.O. Box: 14565-1114, Tehran-Iran

³Advanced Technology Company of Iran, AEOI, P.O. Box: 14399-55431, Tehran-Iran

Email: akeshtkar@aeoi.org.ir

* Corresponding Author

Received 05 March 2020; Accepted 03 May 2020; Publication 08 July 2020

Abstract

The performance of a gas centrifuge that is used for isotopes separation is dependent on the gas flow inside it. In this study, for modeling the UF₆ gas flow, an Implicit Coupled Density-Based (ICDB) solver, was developed in OpenFOAM. To validate the ICDB solver, the gas flow within the rotor in total reflux state was compared with the analytical solution obtained by Onsager model and the numerical solution obtained by the Fluent software. The results showed that the ICDB solver had acceptable accuracy and validity. Also the computational efficiency of Roe, AUSM (Advection Upstream Splitting Method) and AUSM⁺ up schemes were compared. The results showed AUSM⁺ up scheme is efficient. Then, the uranium isotopes separation in a gas centrifuge was simulated. It was revealed that all gas flow characteristics including velocity, pressure, temperature and axial mass flux, as well as uranium isotope separation parameters including separation power and separation coefficients could well be predicted.

Keywords: Gas flow modeling, uranium isotopes separation, gas centrifuge, ICDB solver, OpenFOAM.

1 Introduction

In a gas centrifuge, the strong centrifugal acceleration existing under operating conditions tends to collect the heavier isotopes toward the wall, while the lighter isotopes tend to concentrate toward the axis. In this way, radial separation occurs between the isotopes. However, the separation resulting from the action of the centrifugal field is very limited. In order to increase the separation process, the axial mass flux in the axial direction is generated by combining two mechanisms: mechanical drives and thermal drives [1]. By forming the axial flow, the axial separation takes place inside the rotor, and thus, the rotor’s separation performance increases. Therefore, it is necessary to simulate the behavior of the dynamics of gas flow. The gas flow simulation of gas centrifuges is difficult due to the very high rotational speed of the rotor, which causes the pressure at the wall to be typically 5–6 orders of magnitude larger than that at the center of the rotor [2]. In addition, the other difficulty is due to the shock caused by the collision of the gas flow with the scoop heads [3]. The theoretical investigation of gas flow in the rotor has been made by two different approaches: analytical and numerical. Analytical methods such as long rotor solution [4] and boundary layer (Onsager) solution [5–8] are applied to simulate the gas flow in the rotor. It is to be noted that these problems have been solved with various hypotheses and simplifications.

Nowadays, due to the upgrading of computers in terms of hardware and software, numerical methods are applied to simulate the gas flow in the rotor. In 1976, Kai studied compressible flow in axisymmetric state with numerical integral procedure through Navier-Stokes equations. He used the wall temperature gradient and feed to generate axial mass flux [9]. In 1977, Soubbaramayer et al. using finite element method, studied the velocity, temperature, pressure and concentration of isotopes inside the rotor of centrifuge in an axisymmetric state. The calculations were performed with the FORTRAN code [10]. In 1984, Ribando examined the gas flow inside the centrifuge in axisymmetric state using the finite difference method, in which linear Navier-Stokes equations were solved by time-marching method [11]. In 2000, Borisevich et al. examined the numerical solution of Navier-Stokes equations in axisymmetric state. The calculations showed that the dependence of the machine’s separation power on geometric parameters and the position of the bellows are as important as the scoop drag force and feed flow. They used finite difference method in their study [12]. Zeng et al. (2006) examined the simulation of gas flow inside the centrifuge using a computational fluid dynamics (CFD) method. They used finite volume method and upwind implicit second-order equation to solve Navier-Stokes equations [13]. In 2013, Bogovalov et al. used a numerical method to study the gas flow and the isotopes separation in axisymmetric state. The applied method’s accuracy was compared with that of semi-analytical methods [2].

In all of these studies, the performance of centrifuge has been investigated by applying numerical methods. Therefore, the use of advanced methods, codes and computing tools, as well as their development is essential for simulating the gas flow inside the rotor. One of the computing tools is OpenFOAM software. OpenFOAM is an open source and free software that has the ability to modify and change codes through adding new equations and models [14–17]. In this software, there is the ability to solve the molecular region using the Direct Simulation Monte Carlo (DSMC) method and also develop the solver for continuous and molecular region couplings [18–20]. It is worth noting that the codes of this software are written in C++. In 2013, the capabilities of CFD solvers in OpenFOAM software for supersonic air flows were verified by Chun et al. [21]. Comparison of the obtained results with experimental data confirmed the high accuracy of these solvers. In 2016, Farber et al. developed Navier-Stokes and DSMC equations for simulating diluted gas flows for production of organic OLED using OpenFOAM software. Finally, they compared their results with experimental data, which yielded good results [22]. In 2018, White et al. performed air flow simulations on a variety of items, including supersonic air flow on a flat surface, free molecular air flow on a cylinder using the DSMC method in OpenFOAM software. The results indicated the good accuracy of the DSMC solver in OpenFOAM software comparing to other methods [21].

The purpose of the present study is modeling the uranium hexafluoride (UF₆) gas flow inside the rotor of centrifuge in axisymmetric state using OpenFOAM software. The simulation is performed using a finite-volume numerical method by developing Implicit Coupled Density-Based (ICDB) solver to solve Navier-Stokes equations in OpenFOAM framework. First, the ICDB solver is verified with analytical solution (Olander’s study) and numerical solution (Fluent) for UF₆ gas flow within a rotor in total reflux state. Also the computational efficiency of AUSM⁺ up scheme is investigated. Then the simulation of gas flow inside a rotor is done by considering the feed flow and thermal and mechanical drives, such as baffle and scoop. The purpose of the numerical simulation of the gas flow inside the rotor is to calculate the velocity, pressure and temperature distribution in order to determine the distribution of uranium isotopes’ concentration within the rotor. By calculating the distribution of uranium isotopes’ concentration, the separation coefficients and separation power, which are amongst the important parameters for studying the rotor’s separation performance, will be determined. In this study, UF₆ gas flow into the centrifuge’s rotor is considered axisymmetric, steady and laminar.

2 Governing Equations

2.1 Hydrodynamic Equations

The governing equations on the UF6 gas flow inside the centrifuge rotor are the Navier-Stokes equations. These equations in axisymmetric and steady states are as follows [23]:

\frac{1}{r} \frac{\partial}{\partial r} (r u ρ) + \frac{\partial}{\partial z} (w ρ) = 0

(1)

\begin{aligned} ρ (u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} - \frac{v^{2}}{r}) & = - \frac{\partial p}{\partial r} + \frac{\partial}{\partial r} [μ (- \frac{2}{3} \nabla \cdot \bar{V} + 2 \frac{\partial u}{\partial r})] \\ + \frac{\partial}{\partial z} [μ (\frac{\partial u}{\partial z} + \frac{\partial w}{\partial r})] + \frac{2 μ}{r} (\frac{\partial u}{\partial r} - \frac{u}{r}) \end{aligned}

(2)

\begin{aligned} ρ (u \frac{\partial v}{\partial r} + w \frac{\partial v}{\partial z} + \frac{u v}{r}) & = \frac{\partial}{\partial z} [μ (\frac{\partial v}{\partial z}) + \frac{\partial}{\partial r} [μ (\frac{\partial v}{\partial r} - \frac{v}{r})] \\ + \frac{2 μ}{r} (\frac{\partial v}{\partial r} - \frac{v}{r}) \end{aligned}

(3)

\begin{aligned} ρ (u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z}) & = - \frac{\partial p}{\partial z} + \frac{\partial}{\partial z} [μ (- \frac{2}{3} \nabla \cdot \bar{V} + 2 \frac{\partial w}{\partial z})] \\ + \frac{\partial}{\partial r} [μ (\frac{\partial w}{\partial r} + \frac{\partial w}{\partial z})] + \frac{μ}{r} [\frac{\partial u}{\partial z} + \frac{\partial w}{\partial r}] \end{aligned}

(4)

\begin{aligned} ρ c_{p} (u \frac{\partial T}{\partial r} + w \frac{\partial T}{\partial z}) & = \frac{1}{r} \frac{\partial}{\partial r} (k r \frac{\partial T}{\partial r}) + \frac{\partial}{\partial z} (k \frac{\partial T}{\partial z}) \\ + β T [u \frac{\partial p}{\partial r} + w \frac{\partial p}{\partial z}] + \emptyset \end{aligned}

(5)

where,

\begin{aligned} \emptyset & = 2 μ [{(\frac{\partial u}{\partial r})}^{2} + {(\frac{u}{r})}^{2} + {(\frac{\partial w}{\partial z})}^{2}] \\ + μ [{(\frac{\partial v}{\partial r} - \frac{v}{r})}^{2} + {(\frac{\partial v}{\partial z})}^{2} + {(\frac{\partial w}{\partial r} + \frac{\partial u}{\partial z})}^{2}] \\ - \frac{2}{3} μ {(\frac{1}{r} \frac{\partial (r u)}{\partial r} + \frac{\partial w}{\partial z})}^{2} \end{aligned}

The equation of state for ideal gas is as follows:

p = (γ - 1) ρ I, γ = \frac{c_{p}}{c_{v}}, I = c_{p} T

(6)

where, ρ is the density, u, v and w are the velocity in radial, azimuthal and axial directions, respectively, p is the pressure, k is the conductivity, and μ is the viscosity, β is the thermal expansion coefficient, T is the temperature, and c_p and c_v are the specific heat capacity at constant pressure and volume. In this study, the UF₆ gas flow into the centrifuge rotor has been considered steady.

Note that in the conventional manner, the analysis of flow in a gas centrifuge assumes that the motion of gas is axisymmetric. Thus, the azimuthal derivatives vanish in the equations, the viscosity, thermal conductivity and specific heat of the UF₆ gas are constant throughout a cylinder, the cylinder surfaces are perfect conductors of heat, the gravitational acceleration is negligible compared to the centrifugal one, and the process gas is an ideal gas.

2.2 Mass Transport Equation

In order to obtain the distribution of uranium isotopes’ concentration in the rotor of centrifuge, the mass transport equation for isotopes is used as follows [24]:

\nabla \cdot {\vec{j}}_{A} + \nabla \cdot (C_{A} \vec{V}) = 0

(7)

where, the molecular concentration of isotope A in the mixture is equal to:

C_{A} = x_{A} ρ

(8)

In which, x_A is the molar fraction of isotope A, and ρ is the density.

Diffusion flux ${\vec{J}}_{A}$ will be obtained based on the combination of the pressure diffusion and the concentration diffusion flux. Thermal and force diffusions are not discussed here.

Pressure diffusion:

{\vec{j}}_{A, P} = - \frac{D_{A B}}{R T_{0}} (x_{A} (1 - x_{A}) \frac{△ M}{M_{A}}) \nabla P

(9)

Concentration diffusion:

{\vec{j}}_{A, C} = - ρ D_{A B} \nabla x_{A}

(10)

Here, D_AB is the diffusion coefficient of isotope A in isotope B, ΔM is the difference between the molecular mass of the light isotope (M_A) and the molecular mass of the heavy isotope (M_B), R the universal gas constant, and T₀ is the average gas temperature.

For separation of uranium isotopes, the chemical composition of a mixture may be expressed in terms of mass fraction or molar fraction corresponding to each isotope. The z_F sign represents the molar fraction of the desired isotope in the input feed flow, y_P is the molar fraction of the desired isotope in the enriched flow, and x_W is the molar fraction of the desired isotope in the tail flow that are removed from the separation unit.

The separation coefficients are defined as follows [25]:

α = \frac{(\frac{y_{P}}{1 - y_{P}})}{(\frac{x_{W}}{1 - x_{W}})}

(11)

β = \frac{(\frac{y_{P}}{1 - y_{P}})}{(\frac{z_{F}}{1 - z_{F}})}

(12)

γ = \frac{(\frac{z_{F}}{1 - z_{F}})}{(\frac{x_{W}}{1 - x_{W}})}

(13)

where, α, β and γ are separation, enrichment and stripping factors, respectively.

Also the ratio of product flow rate to feed flow rate is called “cut coefficient” and defined as follows:

θ = P / F

(14)

The separation power can be calculated as follows [25]:

\begin{aligned} δ U & = P (2 y_{P} - 1) l n \frac{y_{P}}{1 - y_{P}} + W (2 x_{W} - 1) l n \frac{x_{W}}{1 - x_{W}} \\ - F (2 z_{F} - 1) l n \frac{z_{F}}{1 - z_{F}} \end{aligned}

(15)

Where, F, P and W are the mass flow rate of the feed, enriched, and tailed flows, respectively.

3 Numerical Solution

To solve the governing equations and simulate the gas flow in the centrifuge rotor, we use a solver with special capability that has been developed in OpenFOAM by the present authors [26]. This solver is called ICDB (Implicit Coupled Density-Based) solver.

ICDB solver is based on the AUSM⁺ up scheme using LU-SGS (Lower-Upper Symmetric Gauss–Seidel) algorithm and GMRES (Generalized Minimal Residual). LU-SGS time scheme is used to perform the time march of Navier-Stokes equations. By using the AUSM⁺ up scheme, a wide range of Mach number present in the centrifuge rotor (such that it is low near the axis and high near the wall of the rotor) is simulated. ICDB solver is uniformly valid for all speed regimes. This solver has been validated for supersonic and hypersonic flows for the three cases of airflow over a flat plate, a 2-D cylinder, and a prism [26].

To our knowledge, there is no modeling and simulation of the uranium gas flow using found in the literature the ICDB solver in OpenFOAM.

4 Simulation Rotor

The main purpose of the present work is the simulation of the UF₆ gas flow within a rotor with radius “a” and height “H” that rotates with a constant velocity. To simulate the gas flow, the full Navier-Stokes equations are applied, which are solved using finite volume method in OpenFOAM.

The specifications of UF₆ as working gas are described in Table 1.

Table 1 Basic parameters of the rotor [2]

Parameter	Value	Parameter	Value
Molecular weight (kg/mol)	0.352	μ (Pa.s)	1.83*10⁻⁵
ρD (Pa.s)	2.306*10⁻⁵	$γ = \frac{C_{p}}{C_{v}}$ (−)	1.067
Thermal conductivity (W/m.K)	0.0061	Cp(J/kg·K)	385
Diameter of UF₆ molecule (m)	6*10⁻¹⁰	Molar fraction of feed for UF₆ 235 (−)	0.00711
Radius of rotor (mm)	100	Height of rotor (mm)	1000

4.1 Verification of ICDB Solver for Gas Flow in Gas Centrifuge

In order to evaluate the accuracy and validity of the solver to simulate high speed rotary flows, simulation of a rotor in a total reflux state is performed. For this purpose, a rotor with the length of 3.35 m and the radius of 0.0915 m, the wall pressure of 13400 Pa, the wall linear velocity of 700 m/s, and the average UF₆ gas temperature of 300 K is considered [27].

The simulation UF₆ gas flow in the rotor is done by analytical solution (Onsager equation) as implemented by Olander [27]. The other simulation is numerical solution that is obtained by Fluent in present time. The numerical solution of Fluent is done by using the ROE scheme. The difference between the numerical study (Fluent and OpenFOAM) and the analytical study is that the Onsager equation is solved by considering simplifications in terms of equations, while in the numerical solution, the full Navier-Stokes equations are completely solved.

In Figure 1, the distribution of the axial mass flux in total reflux state simulated by OpenFOAM and Fluent is compared with the simulated results by Olander in the rotor’s mid-section (Z/H = 0.5). Note that in plotting of graphs, the scale height $ξ = A^{2} (1 - \frac{r^{2}}{a^{2}})$ ) is used in which A² is the velocity parameter $A^{2} = \frac{Ω^{2} a^{2}}{2 \bar{R} T_{0}}$ , a is the radius of the rotor’s wall, r is radius, T₀ is the average gas temperature, Ω is the angular velocity, and $\bar{R}$ is the mass-specific gas constant [27]. By comparing the results, the accuracy of the simulation performed with the ICDB solver is acceptable and the accuracy of the simulation is visible.

Figure 1 The distribution of the countercurrent of UF₆ gas flow in total reflux state in a centrifuge (Z/H = 0.5) [32].

It was found that the ICDB solver is capable of simulating high speed rotary flows and can be used to simulate the flow inside the rotor.

4.2 Modeling a Rotor Using Roe, AUSM and AUSM⁺ Up Schemes

Computational fluid dynamics codes are becoming a promising tool for regular and routine use in engineering. Hence, more careful attention must be paid for developing a numerical scheme that is reliable for a wide range of applications like physical modeling [28]. One of the important extensions is to allow the application of the existing compressible flow codes to reliably predict low and high speed flows.

Flux difference splitting (FDS) schemes such as the Roe scheme [29] (The Roe scheme, devised by Phil Roe, is an approximate Riemann scheme based on the Godunov scheme) have very high resolution for both contact discontinuity and boundary layer in compressible flows. Flux vector splitting (FVS) schemes have much better robustness in capturing strong discontinuities; however, they have a large numerical dissipation on contact discontinuities and in boundary layers. The AUSM (Advection Upstream Splitting Method) family schemes enjoy the advantages of both FDS and FVS schemes like high resolution for contact discontinuity, low numerical dissipation, and high computational efficiency [30].

A new version of the AUSM-family schemes, called AUSM⁺ up based on the low Mach number asymptotic analysis, is described in Liou’s study [28]. It has been demonstrated to be reliable and effective not only for low Mach numbers, but also for all speed regimes, as well as for a wide variety of flow problems over different geometries and grids [28]. In this scheme, central and upwind interpolations are used for subsonic and supersonic flows, respectively [31]. The equations related to Roe and AUSM⁺ up schemes are found in Refs. [28, 29].

Figure 2 The distributions of axial mass flux obtained from Roe, AUSM and AUSM⁺ UP schemes.

In this section, the simulation of the UF₆ gas flux inside the rotor with radius “a” and height “H” using Roe, AUSM and AUSM⁺ up schemes is performed in OpenFOAM framework. The rotor is considered in total reflux state with thermal drives. The wall velocity is 6000 rad/s and thermal derive is 20 K.

In Figure 2 the distributions of axial mass flux obtained from Roe, AUSM and AUSM⁺ up schemes are compared; as shown, the difference between the schemes is small.

The number of computational cells in this model is 75000 and the width of the smallest cell near the wall is 2 micrometers. The computational efficiency of the Roe, AUSM and AUSM⁺ up schemes is calculated. The features of the computing system for performing the above mentioned solutions are considered to be identical and with 30 cores with a specific processor (Intel 2.2 GHz).

The computational efficiency of different schemes is presented in Table 2.

Table 2 Computational efficiency of different schemes

Scheme	Roe	AUSM	AUSM⁺	UP
Time (Hour)	52	45	40

It is shown that the best computational efficiency belongs to AUSM⁺ up scheme (The computational time of AUSM⁺ up scheme is 40 hour).

Based on the above results, ICDB solver based on AUSM⁺ up scheme can be used as reliable solver to simulate the gas flow in the rotor of centrifuge.

4.3 Rotor Including Feed, Scoops and Baffle

In this section, the simulation of the UF₆ gas flow inside the rotor with radius “a” and height “H” in feed flow state, and thermal and mechanical drives is considered simultaneously.

The schematic diagram of the rotor geometry and the type of its boundaries including the feed, product scoop, waste scoop, and baffle is shown in Figure 3. The grid is also denser toward the rotor walls, where the flow variables have large gradients. It is worth mentioning that all geometric dimensions of the rotor are dimensionless (height is Z = z/H and radius is R = r/a).

The present simulation is considered in axisymmetric state, in which changes in angular position are ignored, so the derivative in the direction of the angle is equal to zero: ∂/∂θ = 0. In total reflux, due to the perfect symmetry of geometry, axisymmetric assumption can be well applied; however, in the case of rotor with scoop, the axisymmetric assumption should be modified in relation with the scoops. The scoops are the only 3-D parts of the centrifuge. For modeling the scoops, it is necessary to retain the two major roles they play in the separation process. First, gas extraction is allowed from the rotor and then it induces the mechanical drives from creating the axial gradient of the concentration of the different isotopes. Thus, the modeling of the scoops needs to reflect these two roles. In the present study, the scoop is modeled with a disk instead of a pipe. To modify the modeling of scoop, the drag force (F_D) created by the scoop in a 3-D model is considered to be equivalent with the shear force (F_{τ_w}) applied to the disk surface in axisymmetric state:

F_{D} = F_{τ_{w}}

(16)

Figure 3 (a) Computational domain of the rotor, and (b) Structured mesh in the fluid domain of the rotor.

The drag force is:

F_{D} = \frac{1}{2} C_{D} ρ v^{2}

(17)

The coefficient C_D may itself depend on the obstacle geometry and on the normal Mach number of the gas. For this purpose, we average this force over the 2π radians of the centrifuge, in the case of scoop with normal incidence.

Here, the drag force applied as a drive in a 3-D state is equal to the angular velocity variation (ΔΩ), so that by changing the angular velocity, the scoop drive can be modified in axisymmetric state.

Table 3 Boundary conditions of the rotor

Boundary	Velocity	Temperature
Wall rotor	ω = 6000 (rad/s)	$T (z) = [300 + (\frac{z}{H} * 20)]$ (K)
Axis	(0 0 0) (m/s)	310 (K)
Top cap	ω = 6000 (rad/s)	320 (K)
Bottom cap	ω = 6000 (rad/s)	300 (K)
Baffle	ω = 6000 (rad/s)	300 (K)
Waste scoop	ω = 4000 (rad/s)	$\frac{\partial T}{\partial n} = 0$ (K/m)
Product scoop	ω = 4500 (rad/s)	$\frac{\partial T}{\partial n} = 0$ (K/m)

In the present study, the boundary conditions are determined for the velocity, pressure and temperature variables.

The Flow Rate Outlet Velocity boundary condition is considered for the exhaust gas from the product and the waste scoops, and the values of the flow rates are set using the cut coefficient (θ).

To apply the fixed inlet flow in the rotor through the feed entrance boundary, pressure, temperature and velocity are used as fixed value conditions, where the boundary conditions for the feed are equal to 86.65 m/s (in the direction x), 300 K and 300 Pa for velocity, temperature and pressure, respectively. The other boundary conditions applied are according to Table 3.

4.3.1 Grid independence study

There are rotating and stationary inserts in the centrifuge rotor for the inflow and outflow of gas. These inserts have to be incorporated as solid boundaries in a simulation, resulting in a complex simulation domain, which in turn increases the complexity of the grid, and the number of cells at which the hydrodynamic fields are evaluated. To study the grid independence, four cases with a different number of cells in the radial and axial directions are selected; the number of cells from small to large is equal to 25913, 109288, 255636, and 547990, respectively. The characteristics of these grids and the size of the cells near the wall rotor are shown in Table 4.

In this part, the flow characteristics including temperature, pressure, and axial mass flux within the rotor are studied.

In Figure 4, the flow characteristics are shown at a section near baffle. Pressure, temperature and axial mass flux distributions are obtained in the rotor for different grids.

Table 4 Characteristics of the grids

Case	Number of cells	Size of the cells near the wall rotor (μm)
1	25913	2
2	109288	2
3	255636	1
4	547990	0.5

It was found that by changing the cell size, no effect was observed in the pressure distribution. Also by increasing the number of cells, the difference in results related to temperature and axial mass flux decreased. It is clear that the obtained results from 25913 and 109288 cells for temperature and axial mass flux are about 5% and 30%, respectively. While there is no difference between the obtained results of 255636 and 547990 cells, the results are quite consistent with each other.

By reviewing the studies on flow properties such as pressure, temperature and axial mass flux, it can be concluded that the grid with 255636 cells has a good accuracy for simulating the behavior of the gas inside the rotor.

4.3.2 Characteristics of UF₆ gas flow in the rotor

According to the grid independence study, the grid (255636 cells) with 789 cells in radial direction and 324 cells in axial direction (size of the cells near the wall rotor is 1 μm) have an acceptable accuracy and the results obtained are independent of the size of cells. The specifications of the computational system for present simulation are 30 cores with a specific processor (Intel 2.2 GHz). The computational time on convergence of the problem by this system is about 120 hours.

After the convergence of the numerical solution and reaching the maximum residual 10⁻¹⁰, the results are extracted. In this part, the velocity, temperature, pressure and axial mass flux distributions within the rotor are shown at various sections.

Figure 5 illustrates the pressure contour and the pressure distributions in logarithmic scale in different sections of the rotor as well as near the baffle. As seen, due to the presence of centrifugal force, the gas has been accumulated near the wall and has led to an increase in pressure variations in this area. It is noteworthy that the maximum pressure on the wall top of the baffle is visible and its ratio is about 4 times larger than in the bottom of the baffle. The value of wall pressure is 2200 Pa and the value of wall pressure in the bottom of the baffle is 539 Pa. In a centrifuge, there is at least one baffle in each machine to remove the effect of one of the scoops on the axial flow, which is separated by a baffle from the main separation chamber. The presence of baffle is essential for the formation of desired axial flow.

Figure 4 Distribution of pressure (a), temperature (b) and axial mass flux (c) in different grids.

Figure 5 Comparison of the pressure distributions along the rotor.

The temperature, azimuthal velocity, and axial mass flux distributions along the rotor are shown in Figure 6.

Figure 6 Comparison of the temperature (a), azimuthal velocity (b), and axial mass flux (c) distributions along the rotor.

The temperature distributions are shown in Figure 6(a). As observed, due to the viscous stress in the gas layers, mechanical energy is converted into thermal energy, which increases the temperature of the gas inside the rotor. The axial temperature is fixed at 310 K and the temperature gradient rotor is in the range of 300–320 K.

Figure 6(b) shows the distribution of the swirl velocity (υ_θ) in three sections of the rotor. As it is seen, the UF₆ gas velocity along the wall is equal to the its velocity in the wall, and as the rotor wall moves away, the velocity values decrease in the radial direction. It should be noted that in the dense and continuous regions (near the rotor wall), the velocity changes are linear and proportional to the radius. While in rarified regions, variations in velocity are different with those in the continuous region. In addition, by observing the graph at Z = 0.4 and Z = 0, the effect of scoop waste and feed on the swirl velocity is visible.

The distribution of axial mass flux is also shown in three sections in radial direction in Figure 6(c). Due to the use of feed flow, thermal and mechanical drives to generate the axial flow, it is expected that the axial mass flux in the vicinity of the drives would be larger. As can be seen, the values of axial mass flux in the vicinity of scoop waste and in the middle region are more than in the other regions.

To investigate the effect of driving force in the rotor, the flow function is used, which is defined as follows:

F (r \cdot z) = 2 π \int_{0}^{r} ρ w (\overset{´}{r} \cdot \overset{´}{z}) d \overset{´}{r}

(18)

where, r and z denote the radial and axial positions in the rotor, respectively, and ρw is axial mass flux.

The distribution of flow function in three sections in radial direction is shown in Figure 7. As can be seen, the maximum value of the flow function is 431 kg/s, so the maximum driving force is at the center of the rotor (Z = 0). The second maximum driving force is near the waste scoop (Z = −0.4) and the maximum value of the flow function is 300 kg/s. Hence, it can be said the driving forces including feed flow, thermal and mechanical drives are applied simultaneously. Furthermore, they affect and reinforce each other.

4.3.3 Characteristics of gas flow around the feed inlet

Due to the importance of the feed and its effect on the behavior of the gas inside the rotor, radial velocity, temperature and pressure distributions in the feed region are investigated.

Figure 7 Comparison of the flow function distributions along the rotor.

As explained, the UF₆ gas expands in front of the feed region because there is high pressure difference in this region. With the expansion of the gas, the temperature drops to 285 K and the radial velocity increases to 120 m/s. The pressure, velocity and temperature distributions in the axial direction and in the area around the feed are shown in Figure 8. It is to be noted that since the Navier-Stokes equations are valid under the assumption of flow continuity, while in the area around the feed, the gas flow regime is diluted; therefore, for a more precise examination of this region, it is suggested that the feed region be simulated by the molecular methods such as DSMC.

4.3.4 Concentration distribution

In order to obtain the concentration values inside the rotor, the distribution of flow characteristics is first obtained. Then, by solving the mass transport equation, the distribution of the concentration is computed. The concentration distribution in the axial direction of the rotor is shown in Figure 9. The concentration of the feed for UF₆ (235) is equal to 0.00711.

Given the distribution of the obtained concentration, the separation performance of the rotor is presented in Table 5.

Figure 8 Distributions of the pressure (a), radial velocity (b), and temperature (c) near the feed.

Figure 9 Distribution of the concentration UF₆ (235) in rotor.

Table 5 Separation data of the rotor

Case	Number of cells	δU (K_g UF₆ SWU/yr)	α(−)	β(−)	γ(−)
1	255636	9.580	1.462	1.215	1.204

5 Conclusion

In the present study, using OpenFOAM software, the simulation of the flow of uranium hexafluoride was done by feed the flow, thermal and mechanical drives in axisymmetric state. By applying the drives, an axial flow was created inside the rotor, which is an effective parameter in evaluating the centrifuge’s separation performance.

Due to the high compressibility of gas flow inside the rotor, the results of the solution are highly dependent on the size of the cells; therefore, grid independence study was necessary.

In the case of pressure distribution, it was observed that due to centrifugal force, a large mass of gas was collected along the wall. The maximum value of the wall pressure was 2200 Pa and the minimum value of the pressure near the axis was 0.006 Pa. Also the effect of baffle to decrease the pressure was observed. In the case of temperature distribution, the gas temperature was increased due to the viscous dissipation between the gas layers. Regarding the distribution of angular velocity, its changes in the dense region were found to be linear, which deviated from the linear state in the diluted region. There was also a variation in the rate of changes in the areas around the feed and scoop. The results obtained for axial mass flux revealed that due to the application of the feed flow, thermal and mechanical derives, the axial flow was generated. It was also observed that the value of the axial mass flux in the middle section of the rotor and near the waste scoop was greater than in the other regions. The same trend was shown in the distribution of the flow function. Comparing the variation in radial velocity and temperature in the areas around the feed, it was determined that due to sudden pressure change, the gas expanded, the radial velocity increased to 120 m/s and the temperature decreased to 285 K. After the flow analysis, the separation performance of the rotor was examined and it was found out that the separation power was 9.58 (kg UF₆ SWU/yr.) and its separation factor was 1.462. By investigation of the results, it can be concluded that the ICDB solver has high capability to model and simulate the gas flow in the centrifuge rotor as well. Also it was shown that the best computational efficiency belongs to AUSM⁺ up scheme in ICDB solver.

References

[1] Soubbaramayer, “Centrifugation,” Applied Physics, vol. 35, pp. 183–244, 1979.

[2] S. Bogovalov, V. kislov and I. Tronin, “Verification of numerical codes for modeling of the flow and isotope separation in gas centrifuges,” Computers & Fluids, vol. 86, pp. 177–184, 2013.

[3] V. Borman, S. Bogovalov, V. Borisevich and I. Tronin, “The computer simulation of 3d gas dynamics in a gas,” in Journal of Physics: Conference Series 751, 2016.

[4] J. Hu, C. Ying and S. Zeng, “Overall Separation Factor in a Gas Centrifuge Using a Purely Axial Flow Model,” Separation Science and Technology, vol. 40, pp. 2139–2152, 2005.

[5] D. R. Olander, “The Theory of Uranium Enrichment by the Gas Centrifuge,” Progress in Nuclear Energy, vol. 8, pp. 1–33, 1981.

[6] J. Brouwers, “On Compressible Flow in a Gas Centrifuge and its Effect on the Maximum Separative Power,” Nuclear Technology, vol. 39, 1978.

[7] P. Migliorini, “Modeling and Simulation of Gas Centrifuge Cascade for Enhancing the Efficiency of IAEA Safeguards,” pp. 1–158, 2013.

[8] R. Bourn, T. Peterson and H. Wood, “Solution of the pancake model for flow in a gas centrifuge by means of a temperature potential,” Computer Methods in Applied Mechanics and Engineering, vol. 178, pp. 183–197, 1999.

[9] T. Kai, “Basic Characteristics of Centrifuges, (III) Analysis of Fluid Flow in Centrifuges,” Journal of Nuclear Science and Technology, vol. 14, pp. 267–281, 1976.

[10] Soubbaramayer and J. Lahargue, “A numerical model for the investigation of the flow and isotope concentration field in an ultracentrifuge,” Computer Method in Applied Mechanic and Engineering, vol. 15, pp. 259–273, 1978.

[11] R. J. Ribando, “A finite-difference solution of onsager’s model for flow in a gas centrifuge,” Computers & Fluid, pp. 235–252, 1984.

[12] V. Borisevich, O. Morozov and O. Godisov, “Numerical simulation of bellows effect on flow and separation,” Nuclear Instruments and Methods in Physics Research A, vol. 455, pp. 487–494, 2000.

[13] D. Jiang and S. Zeng, “CFD simulation of 3D flowfield in a gas centrifuge,” International Conference on Nuclear Engineering July 17–20, Miami, Florida, USA, 2006.

[14] L. Silva and P. Lage, “Development and implementation of a polydispersed multiphase flow model in OpenFOAM,” Comput. Chem. Eng., vol. 35, pp. 2653–2666, 2011.

[15] OpenFOAM, “The Open Source CFD Toolbox, User Guide, ESI-OpenCFD Ltd,” 2019.

[16] H. Kassem, X. Liu and J. Banerjee, “Transonic flutter analysis using a fully coupled density based solver for inviscid flow,” Advances in Engineering Software, vol. 95, pp. 1–6, 2016.

[17] S. Damián, J. M. Giménez and M. Nigro Norberto, “gdbOF: A debugging tool for OpenFOAM,” Advances in Engineering Software, vol. 47, pp. 17–23, 2012.

[18] L. Marcantoni, J. Tamagno and S. Elaskar, “High spped flow simulation using OpenFOAM,” Mecánica Computacional, pp. 2939–2959, 2012.

[19] F. Moukalled, L. Mangani and M. Darwish, The Finite Volume Method in Computational Fluid Dynamics (An Advanced Introduction with OpenFOAM^® and Matlab^®), vol. 113, Springer International Publishing Switzerland, 2016.

[20] M. Darbandi and E. Roohi, “A hybrid DSMC/Navier–Stokes frame to solve mixed rarefied/nonrarefied hypersonic flows over nano-plate and micro-cylinder,” International Journal for Numerical Methods in Fluids, vol. 72, pp. 937–966, 2013.

[21] S. Chun, S. Fengxian and X. Xinlin, “Analysis on capabilities of density-based solvers within OpenFOAM to distinguish aerothermal variables in diffusion boundary layer,” Chinese Journal of Aeronautics, vol. 26, pp. 1370–1379, 2013.

[22] K. Farber, P. Farber, J. Gräbel, S. Krick, J. Reitz and P. Ueberholz, “Development and validation of a coupled Navier–Stokes/DSMC simulation for rarefied gas flow in the production process for OLEDs,” Applied Mathematics and Computation, vol. 272, pp. 648–656, 2016.

[23] F. Moukalled, L. Mangani and M. Darwish, The finite volume method in computational fluid dynamics (an advanced introduction with OpenFOAM^® and Matlab^®), vol. 113, Springer International Publishing Switzerland, 2016.

[24] M. Benedict, “Nuclear Chemical Engineering,” Mcgraw-Hill Book Company, vol. chapter 14, 1981.

[25] K. Cohen, The Theory of Isotope Separation, McGraw Book Company, 1951.

[26] V. Ghazanfari, A. Salehi, A. Keshtkar, M. Shadman and M. Askari, “Numerical simulation using a modified solver within OpenFOAM for compressible viscous flows,” European Journal of Computational Mechanics, vol. 28, pp. 541–572, 2020.

[27] H. G. Wood and J. B. Morton, “Onsager’s Pancake approximation for the fluid dynamics of a gas centrifuge,” Fluid Mech, vol. 101, pp. 1–31, 1981.

[28] M. Liou, “A sequel to AUSM, Part II: AUSM+− up for all speeds,” Journal of Computational Physics, vol. 214, pp. 137–170, 2006.

[29] P. Roe, “Approximate Riemann solvers schemes,” Computer Physic, vol. 43, pp. 357–372, 1981.

[30] S. Chun, S. Fengxian and X. Xinlin, “Analysis on capabilities of density-based solvers within OpenFOAM to distinguish aerothermal variables in diffusion boundary layer,” Chinese Journal of Aeronautics, vol. 26, pp. 1370–1379, 2013.

[31] K. Kitamura and A. Hashimoto, “Reduced dissipation AUSM-family fluxes: HR-SLAU2 and HR-AUSM+− up for high resolution unsteady flow simulations,” Computers & Fluids, vol. 126, pp. 41–57, 2016.

Biographies

images

Valiyollah Ghazanfari earned a Ph.D. in nuclear engineering from Nuclear Science and Technology Institute in 2020, Iran. His research focuses on thermodynamic, fluid mechanic, CFD and numerical simulation using Open-FOAM and Fluent. He is currently working in Advanced Technology Company of Iran and Materials and Nuclear Fuel Research School.

images

Ali Akbar Salehi received a Ph.D. in nuclear engineering from Massachusetts Institute of Technology in 1977. He is full professor and was chancellor of Sharif University of Technology. His research focuses on Theoretical Physics. He is currently head of the Atomic Energy Organization of Iran.

images

Ali Reza Keshtkar earned a Ph.D. in chemical engineering from Tehran University, Iran. He is full professor and his research interests include design and analysis of separation processes. He is currently head of the Material and Nuclear Fuel Research School.

images

Mohammad Mahdi Shadman received a Ph.D. in chemical engineering from Tarbiat Modarres University, Iran. His research focuses on thermo-kinetic and fluid-mechanic in chemical engineering. He is currently working in Advanced Technology Company of Iran.

images

Mohammad Hossein Askari earned a Ph.D. in mechanical engineering from Tehran University, Iran. His research focuses on fluid mechanic, numerical simulation and CFD method using Fluent software. He is currently working in Advanced Technology Company of Iran.