Large eddy simulation of unsteady flow in gas–liquid separator applied in thorium molten salt reactor

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Large eddy simulation of unsteady flow in gas–liquid separator applied in thorium molten salt reactor

Jing-Jing Li，

Ya-Lan Qian，

Jun-Lian Yin，

Hua Li，

Wei Liu，

De-Zhong Wang

Nuclear Science and Techniques

Vol.29, No.5

Article number 62

Published in print 01 May 2018

Available online 29 Mar 2018

DOI：10.1007/s41365-018-0411-3

49700

Axial gas–liquid separators have been adopted in fission gas removal systems for the development of thorium molten salt reactors. In our previous study, we observed an unsteady flow phenomenon in which the flow pattern is directly dependent on the backpressure in a gas–liquid separator; however, the underlying flow mechanism is still unknown. In order to move a step further in clarifying how the flow pattern evolves with a variation in backpressure, a large eddy simulation (LES) was adopted to study the flow field evolution. In the simulation, an artificial boundary was applied at the separator outlet under the assumption that the backpressure increases linearly. The numerical results indicate that the unsteady flow feature is captured by the LES approach, and the flow transition is mainly due to the axial velocity profile redistribution induced by the backpressure variation. With the increase in backpressure, the axial velocity near the downstream orifice transits from negative to positive. This change in the axial velocity sign forces the unstable spiral vortex to become a stable rectilinear vortex.

Swirl flowThorium molten salt reactorComputational fluid dynamicsLarge eddy simulation

1. Introduction

In the development of the fission gas removal system applied in thorium molten salt reactors (TMSRs), an axial-type swirl tube is used as a gas–liquid separator (see Fig. 1). The separator consists of the swirl vane, swirl chamber, and recovery vane. When a bubbling flow moves though the separator from the right-hand side to the left-hand side, bubbles concentrate in the center of the swirl chamber and coalesce into an air core; the concentrated gas–liquid mixture is removed from the upstream and downstream orifices. In previous experiments [1, 2], during the air core establishment process, the authors observed that the two-phase flow regime illustrated in Fig. 2 exhibits a three-stage phase as the separator outlet backpressure increases. The flow development is similar to the well-known procession vortex core (PVC) phenomenon [3], which is encountered in various swirl flow applications such as the swirl combustor [4], cyclone [5], vortex amplifier, and vortex diode [6]. In the PVC, the center of the swirl vortex deviates from the geometric center of the swirl device concerned, and rotates like a helix. Because the PVC plays a key role in swirl flow devices, the underlying mechanism dominating the flow transition has attracted the interest of several researchers. A comprehensive review of oscillation mechanisms and the role of the PVC in swirl combustion systems was reported in a study by Nicholas [3]; it stated that the occurrence of the PVC is related to the swirl number (S) and the presence of a central recirculation zone, as well as the mode of the fuel entry, combustor configuration, and equivalence ratio. Pisarev et al.[5] studied the "end of vortex" phenomenon in a reversed flow centrifugal separator, and concluded that the instability observed could be related to Reynolds number and swirl chamber length. Experimental studies conducted by Shtork et al. [7], using a high-speed camera and laser Doppler velocimetry on swirl flows in a lean premixed swirl-stabilized combustor, and a gas–liquid cylindrical cyclone separator observed by Hreiz et al. [8], confirmed that the occurrence of PVC is strongly related to a variation in the axial velocity profile. Furthermore, Alekseenko et al. [8] pointed out that the commonly used swirl flow parameters (Reynolds and swirl numbers) do not uniquely characterize the flow structure and the axial velocity profile is another key parameter affecting swirl flow pattern evolution. Moreover, Ragab and Sreedhar [10] used a numerical simulation to demonstrate that large-scale helical sheets of vorticity will be produced when a swirl vortex that has axial velocity deficits is affected by little disturbances.

Fig. 1

2D geometry of axial gas–liquid separator

Fig. 2

(Color online) Two-phase flow evolution with S = 1.17 and Re = 119 891

Based on the above review, a common qualitative conclusion for explaining the PVC phenomenon is the variation in the axial velocity component [11]. Inspired by this concept, the authors felt that by conducting a numerical study on the evolution of flow field in the gas–liquid separator with increasing backpressure, they could produce results that will be significant in interpreting flow unsteadiness. Due to this hypothesis, an unsteady flow simulation based on the large eddy simulation (LES) method was carried out.

2. Numerical modeling

2.1. Geometry, grid, and boundary conditions

Computational fluid dynamics (CFD) models were constructed in order to conform to the configuration of swirl tubes used in the experiments. All geometrical and operational parameters of the swirl tube, including the shapes of swirl and recovery vanes, were carefully copied in the CFD models. The 3D models were built using a commercial program called "universal graphics (UG) Graphix," and discretization was carried out using integrated computer engineering and manufacturing (ICEM) computational fluid dynamics (CFD). Fig. 3 provides an overview of the grid geometry.

Fig. 3

(Color online) Grid geometry and boundary conditions in LES computation

The grid density was controlled in such a manner that the y+ value of the first grid point of all the walls concerned was less than 1, while the swirl chamber’s core zone was refined in order to resolve the vortex motion more accurately. The total number of cells was approximately 10,000,000, including approximately 2,000,000 cells in the swirl vane domain, 2,000,000 in the recovery vane domain, and 6,000,000 in the swirl chamber domain. Pure water was selected as the working fluid for the simulations. The temperature was taken as 293 K, and the properties of the water were as follows: density = 997 kg/m³; dynamic viscosity = 8.899 × 10 ^–4 Pa.s; and molecular mass = 18.02 kg/kmol. By considering a state of an unsteady simulation, the time step was set to 0.0005 s; this guarantees that the Courant number is around 1. A simulation of 0.25 s of real time is sufficient for obtaining a regular unsteady flow. In terms of the boundary conditions, a velocity inlet with V_in = 2.41 m/s (Re = 119, 891) was set at the separator inlet, while the opening boundary conditions with atmospheric static pressures were set at the outlets of the upstream and downstream orifices. At the separator outlet, a linear pressure outlet assumed by Eq. 1 was adopted:

P_{out} = 1013250 \times t [P a]

(1)

2.2 Turbulence modeling

Appropriate turbulence model selection is a crucial factor in numerical simulations of confined swirling turbulent flows. Based on comprehensive CFD studies [12-17] on gas-solid cyclones, it has been established that standard eddy viscosity-based models cannot predict the velocity profiles for swirl flows; this is as a result of their isotropic modeling of Reynolds stresses. Second-order closure models like the Reynolds stress model (RSM) can predict the averaged velocity profiles, which conform to some reasonable degree with experiments. Regarding the capability of predicting unsteady characteristics, no models can achieve better performance than LES, which was validated by Pisarev et al. [5] during their unsteady numerical simulations at the end of vortex evolution in cyclones. LES involves a 3D time-dependent computation of the large-scale turbulent motions that are mainly responsible for turbulent mixing, while those with scales smaller than the computational grid are parameterized; that is, LES solves in the filtered velocity field where the filtering operation allows for separation of the fluid motion scales at the grid level, with small motion scales being considered by a subgrid-scale model. The filtered velocity field is computed as the solution of the filtered Naiver-Stokes equations:

{\begin{cases} \frac{\partial {\bar{u}}_{i}}{\partial x_{i}} = 0, \\ \frac{\partial {\bar{u}}_{i}}{\partial t} + \frac{\partial {\bar{u}}_{i} {\bar{u}}_{j}}{\partial x_{i}} = - \frac{1}{ρ} \frac{\partial \bar{p}}{\partial x_{i}} + \frac{\partial}{\partial x_{j}} ((v + v_{s g s}) \frac{\partial {\bar{u}}_{i}}{\partial x_{j}}) \end{cases},

(2)

where the additional eddy viscosity ν_sgs must be modeled in order to close the system. The localized dynamic Smagorinsky model was used in this study, in which the eddy viscosity is calculated as follows:

v_{s g s} = {(C Δ)}^{2} | \bar{S} (x, t) |,

(3)

where $| \bar{S} (x, t) |$ denotes the norm of the filtered strain rate tensor. The grid size is defined as $Δ = （ V o l)^{1 / 3}$ ; "Vol" is the cell volume, and the coefficient, C, is dynamically computed following Germano’s definition [18].

2.3 Computation

A three-stage calculation procedure was applied to implement the unsteady simulation. In order to obtain appropriate initial flow conditions, first, a steady CFD result was computed by the RSM turbulence model, and taken as the initial condition for the unsteady LES simulation with a constant backpressure at t = 0. A total time of Δt₁ = 0.15 s was used for the first LES simulation in order to obtain a permanent flow regime. Then, the linear varying pressure outlet was activated and the unsteady flow simulation was initiated, with a time step of 0.0001 s, and a total time of 0.25 s. The computation running in parallel mode in a cluster with 20 CPUs (2.4 GHz) and 96 GB RAM typically required 168 h for the entire process.

3. Results and discussions

Post-processing of the transient results obtained during the total computed time indicates that the time scale used to capture the flow transition is T = 0.15 s. Therefore, four points with an equivalent time step of 0.05 s were defined to describe the unsteady process, and the non-dimensional time was defined as t^* = t/T (where t^* = 0; t^* = 0.33; t^* = 0.67; t^* = 1). The following results are discussed from part to whole. The swirl flow regime was characterized by the vortex region detected by the Q-criterion, which is defined as the second invariant of the velocity gradient tensor and has been widely adopted to illustrate spatial vortices [19]. In order to validate the numerical model, a comparison between the calculated vortex and experimental results obtained by means of visualization [1], is illustrated in Fig. 4. A qualitative agreement can be observed in terms of the vortex core shape. Fig. 5 presents the vortex region evolution with increasing backpressure, from which we can observe a significant change in the vortex core, transiting from a cylindrical shape (t^* = 0) to a double and single helix. Another observation of Fig. 2 reveals that the basic flow phenomenon in terms of the air core shape representing the vortex core shape is captured; again, this validates that backpressure plays a key role in vortex dynamics. From the perspective of a gas–liquid separator’s separation function, we can conclude that unstable vortex patterns such as the cylindrical shape, double helix, and single helix, prevent the bubbles from being separated. In order to bring more clarity on the mode in which the flow field dominates vortex dynamics, the flow field downstream of the vortex region: near the recovery vane, is illustrated by the velocity vector field, streamlines, and a zero axial velocity isosurface in Fig. 6. It can be observed that the velocity field experiences complex flow unsteadiness. With the increase in backpressure, the axial velocity sign changes from negative to positive, which is also depicted by the change in area of the zero axial velocity isosurface. As noted in our previous study, when the initial backpressure approaches the atmospheric pressure level, air entrainment occurs at both the upstream and downstream orifices. In the simulation, the water entrainment results in a negative axial velocity, which initiates prominent secondary vortices near the zero axial velocity surface (see Fig. 6a). The axial velocity magnitude is reduced with increasing backpressure, as indicated in Fig. 6b. The zero axial velocity isosurface transforms into a spiral style, which in turn induces the secondary vortices into an oscillatory mode. The existence of the oscillating secondary vortices explains the periodic double helix. The further increase in backpressure shrinks the negative axial velocity zone from a continuous spiral to dispersed pieces close to the downstream orifice where a single helix is formed (see Fig. 6c). The vanishing of the zero axial velocity surface in Fig. 6d simultaneously drives the periodic single helix into a rectilinear air core. It can be concluded from the above analysis that the axial velocity sign dominates the vortex dynamics. The vortex shape evolution represented by the variation in the air core shape in Fig. 2 is strongly dependent on the change in axial velocity.

Fig. 4

(Color online) Comparison of CFD and experimental result flow patterns

Fig. 5

(Color online) Vortex region evolution illustrated by the Q-criterion

Fig. 6

Velocity distribution evolution near downstream orifice

In order to elucidate the details of the variation in the velocity vector, the velocity distribution extracted from four axial positions (see Fig. 7) with an equivalent axial distance of 0.5 D near the downstream orifice is illustrated in Fig. 8. Qualitatively, with increasing backpressure, the prominent change in velocity is in the axial velocity component, the sign of which varies from negative to positive. In order to obtain a quantitative comparison of the variation in velocity, both the axial and tangential velocity components normalized by the average velocity at the separator inlet are shown in Figs. 9 and 10. It can be clearly identified from Fig. 9 that the initial zone (t^* = 0) with negative axial velocity occupies the range of –0.15 < r* < 0.15 and gradually decreases and diminishes when t^*= 1. The transition of the axial velocity sign results in the vortex motion shown in Fig. 10. When the axial velocity is negative, the vortex center wraps around the centerline; when the axial velocity is positive, the vortex center is confined to the geometrical center of the swirl chamber.

Fig. 7

(Color online) Four axial positions used to extract velocity vectors

Fig. 8

(Color online) Variation in velocity vector at four axial locations

Fig. 9

Variation in non-dimensional axial velocity profile with time

Fig. 10

Variation in non-dimensional tangential velocity profile with time

4. Conclusion

In this paper, we presented a numerical study to investigate the unsteady flow phenomenon in a gas–liquid separator. During the experiment, we observed that the air core in the swirl chamber of the separator exhibited a cylindrical shape, a double helix, a single helix, and a rectilinear shape when the backpressure at the separator outlet gradually increased. In order to explain the underlying mechanism dominating the flow regime transition, a single-phase LES simulation with an artificial boundary condition was carried out to represent the varying backpressure. The vortex region evolution indicated that the change in the flow pattern observed in the experiment can be captured via a numerical approach. The variation in the vortex region was analyzed using the velocity vector distribution. We found that the vortex dynamics are dependent on the change in the axial velocity distribution, which can induce unsteady secondary vortices and an inner counter-rotating vortex. When the axial velocity is positive, the vortex approaches a steady state. It should be noted that, in the numerical study, the assumed boundary condition for backpressure is not based on experimental data because a transient experiment with a fixed Reynolds number and increasing backpressure is not available. Moreover, a comparison between the vortex shape evolution and air core was not presented because a single-phase medium was used in the simulation. A two-phase simulation, which can predict the dispersed flow corresponding to the separation process and gas–liquid interface simultaneously, would be more appropriate and provide high fidelity; however, this simulation is currently very difficult to carry out.

References:

J. Yin, J. Li, Y. Ma, et al.,

Study on the Air Core Formation of a Gas-Liquid Separator

. J. Fluid. Eng.-T. ASME. 137, 91301 (2015). doi: 10.1115/1.4030198