Simulation Investigation on the Flow and Mixing in Ducts of the Rotary Energy Recovery Device

The rotary energy recovery device (RERD) is widely equipped in desalination to reduce the system energy consumption. In this study, the fluid dynamics and mixing performance of a typical structure RERD and a visualization apparatus of a RERD (V-RERD) had been compared by simulation. The effects of rotating components on fluid dynamics and mixing had been researched. Simulation results indicated that a swirling flow can be observed from flow fields in the device duct. In the RERD case, the swirling flow changed its rotating direction in the center of the duct, while in the V-RERD case, its rotating direction was unchanged. Then, a swirling number Sn was applied to characterize the degree of swirl intensity, and its formation mechanism in RERD had been discussed. In addition, the Q criterion was adopted to visualize the three-dimensional flow structures and identify vortex structures in the duct. The evolution of vortices in the working process had been investigated. It was found that vortices significantly affected the mixing performance, and the detached vortex could lead to high turbulence and mixing in the duct. Suppressing the vortex shedding may reduce the flow turbulence and gain a lower volumetric mixing rate.


Introduction
With the rapid growth of the population and the development of industry, an increase in global water usage is expected. Desalination has been extensively employed in the last few decades to ensure the efficiency and purity of seawater treatment, alleviating the stresses of the freshwater supply. Nowadays, Reverse Osmosis (RO) is regarded as the leading and the most optimized desalination process, accounting for 61% of the global share [1,2]. Despite the widespread use, RO is an expensive way to supply freshwater. The production cost of freshwater using the RO is 1.25 to 2.5 times as high as that directly from a freshwater source [3]. To push water through semipermeable membranes in the RO desalination process, the seawater is pressurized by the high-pressure pump to typical pressures ranging about 5.5~8.0 MPa, which requires a large amount of electrical energy. Therefore, the energy recovery device (ERD) was employed in the desalination to reduce energy needs by recovering the remaining pressure energy from the concentrate flow. Nowadays, with the implementation of ERDs, the energy consumption of the desalination process based on RO has reached as low as about 2 kWh/m 3 [4,5].
The rotary energy recovery device (RERD), having a high work efficiency, is one kind of ERDs. The RERD, as shown in Figure 1, has two key components: one is the rotor with numbers of axial ducts for passing through fluids and energy transfer; the other is the end cover for distribution of the low-pressure (LP) or high-pressure (HP) fluid. The sole rotating component in the RERD is the rotor, rotating around the shaft. With the rotor rotation, the RERD experienced three working stages as shown in Figure 1. When the RERD works at the seawater supercharging stage, the seawater stream would firstly experience an acceleration stage (the fluid velocity increased significantly) and then a stabilization stage (the fluid kept a steady velocity). In this seawater supercharging stage, the seawater stream with low pressure will meet the high-pressure brine stream in the duct, and it will be pushed out by the brine stream and recovered pressure energy from the brine stream. Hence, the RERD can decrease the electric power consumption in the desalination process. More details on the working principle of the RERD can be found in References [6,7].
Due to no physical piston in the duct of the RERD [8], it can achieve the highest efficiency and lowest lifecycle cost, compared to other types of ERDs [9,10]. But it was unable to avoid the mixing of brine and seawater streams. As a consequence, the salinity of the supercharged seawater stream would increase, which would cause a high energy consumption of the RO system [7]. Therefore, several researchers have investigated the mixing process in the RERD to design a high-efficiency RERD. Liu et al. [11] and Zhou et al. [12] investigated the mixing zone formed in the RERD with simulations. Their results showed that the mixing zone can significantly affect RERD performance. Xu et al. [13] performed a numerical study to research salinity concentration distribution and investigate the effect of operational conditions on the RERD performance. Yin et al. [14] discussed the mixing process of a RERD integrated with a pressure boost device and optimized its structure parameters and working conditions. Cao et al. [15] researched the salinity concentration distribution in a RERD. They found that the flow at a low turbulence intensity in the duct was beneficial to reduce the seawater salinity.
The mechanism of the mixing process and the salinity concentration distribution in the RERD has attracted much interest in recent years, but it is still not fully clear. Moreover, those studies [11][12][13][14][15] have difficulty in validating their simulation results, especially the flow pattern and mixing process, with comparative experimental data. In our previous research [6,16,17], we proposed a visualization apparatus of a RERD (V-RERD) to research. Compared with the typical RERD, rotary end covers and stationary ducts were the only two different components in the V-RERD. Figure 2 shows   2 Geofluids the flow schematic of the V-RERD. Besides the two different components, its structure was the same as that of the typical RERD. Meanwhile, it has the same working principle and working stages as the typical RERD. In those studies [6,16,17]    3 Geofluids results in our experimental research. Moreover, a further study is required to research the effect of the rotating duct of RERD on the mixing process by comparing with the experimental research on the static duct of V-RERD.
To obtain and analyze the mixing process and fluid behavior in the RERD and V-RERD, computational fluid dynamics (CFD) were adopted in this work. The effect of the rotating component on fluid behavior had been discussed. A dimensionless parameter of a swirling number Sn was applied to characterize the degree of swirl intensity. Moreover, the Q criterion was adopted to visualize and identify vortex structures. At last, the influence of the vortex formation on the flow pattern and the salinity concentration distribution had been researched in one working cycle. Figure 3 shows geometric models of the RERD and V-RERD for CFD calculations. Each device consisted of two end covers and 12 ducts. The rotor including 12 ducts was the sole rotating component in the RERD case, while the two end covers were the rotating component in the V-RERD case. Each duct consisted of a square transverse with 30 × 30 mm 2 and a length of 265 mm. Those parameters were from our previous studies [6,16]. Simulation models in two device cases had the same structural parameters. Hence, the simulation models in two device cases could share the same mesh model, as shown in Figure 4. The ICEM CFD v14.5 as a preprocessor of Fluent was used to generate three-dimensional hexahedral grids. The grid independence study had been performed, and the solutions of mesh design of 4:1 × 10 6 were chosen for simulation. Figure 5 shows 12 phases of one duct at the specific duct positions in one working cycle. Phases 2# to 6# represented the seawater supercharging stage. Phases 1# and 7# represented the sealed stage. Phase 8#~phase 12# represented the brine discharging stage.

Governing Equations.
The fluid flow in ERDs was assumed as incompressible and unsteady flow, and computations were carried out as a three-dimensional RANS solution. Considering the slight discrepancy in temperature between the two streams, the heat transfer was ignored in computations. The equations for the continuity and balance of momentum are defined as where p means the static pressure, ρ means the liquid density, ρg ! means the gravitational body force, u ! means the velocity

Geofluids
vector, and τ means the stress tensor. The τ can be described as where I means the unit tensor and μ means the molecular viscosity.
To calculate the mixing process in turbulent flows, the species transport equation is given by where Y i means the mass fraction and J i ! means the diffusion flux. The J i ! is described as where D i,m means the mass diffusion coefficient for species i in the mixture, Sc t means the turbulent Schmidt number, and μ t means the turbulent viscosity.

Numerical Details.
According to our previous work, the Realizable k-ε (RKE) turbulence model showed better agreement with experimental data, capturing the main feature of flow characteristics in the ERD [16]. Thus, the RKE turbulence model was chosen to study. The inlets of devices were set as the velocity boundary condition, and the outlets were set as the pressure boundary condition. Other boundaries were set as no-slip wall boundary conditions. The inflow rate was Q in = 3:6 m 3 /h for the two device cases in this study. The mass fraction of NaCl was 3.5% in the HP inlet and 1.8% in the LP inlet, respectively. The pressure at the HP outlet for the supercharged seawater stream was 6.0 MPa, while the pressure at the LP outlet for the discharged brine stream was 0.2 MPa. As shown in Figure 5, the rotational speed of ducts in the RERD was 120 rpm in the anticlockwise direction, while the rotational speed of end covers in the V-RERD was 120 rpm in the clockwise direction. The sliding

Geofluids
mesh approach was performed to model the movement of ducts or end covers. Sliding mesh interfaces were formed between the duct and end cover domains. The fluid flow was modeled using the standard wall function in the nearwall region.

CFD Solution Methods.
In this study, all the numerical simulations were carried out using Fluent v14.5. The SIMPLE pressure-velocity coupling scheme and the second-order implicit scheme were used in this study. More details of solver settings can be found in Table 1.

Results and Discussion
3.1. Working Performance. The mixing degree in ERDs can be characterized by the volumetric mixing rate M recommended by Stover [7], according to the following definition: As shown in Table 2, the volumetric mixing rate of the two devices showed acceptable values compared to other numerical and experimental results [7,13,15], validating the accuracy of the CFD model in this study. Figure 6 shows comparisons of salinity concentration variations in the HP outlet of the two devices. It can be found that the two curves showed a similar downward trend. As the device started to rotate, the salinity concentration decreased rapidly. Then, a relevant steady salinity concentration was formed at about 6.0 s in V-RERD and 4.0 s in RERD, showing that the mixing process had reached the equilibrium. But there was a discrepancy in Figure 6 that the salinity concentration in the V-RERD case exhibited an approximately periodical variation with a lower fluctuation amplitude. Although these two devices had the same structural parameters and operational parameters, their working performance was a little different as shown in Table 2 and Figure 6. The rotating component, the sole difference factor between the two device cases, may influence the working performance. Therefore, the effect of the rotating component on the fluid dynamics will be discussed in a later section.

Swirling Flow and Vortex Structure.
A swirling flow had been observed in the duct by our previous numerical study [16]. For further research of this phenomenon in the duct of the two devices, flow structures at five transverse planes z 1 = 1/6 L, z 2 = 2/6 L, z 3 = 3/6 L, z 4 = 4/6 L, and z 5 = 5/6 L, as shown in Figure 7, were investigated. Figures 8 and 9 show instantaneous vorticity contours, ω z , in the cross-section of the duct in the V-RERD and RERD cases, respectively. Due to the duct rotation, the relative velocity was used to calculate the relative vorticity and characterize the instantaneous flow field in RERD cases. In Figure 8, a large-scale vortex was clearly observed at z 1 , z 2 , z 4 , and z 5 in the V-RERD cases. It should be pointed out that the direction of the vortex was clockwise (negative vorticity), the same as that of the rotating end cover, showing that the rotating end cover had a significant effect on the flow development at the ends of the duct [6]. But in the center of the duct (z 3 = 3/6 L), the vortex motion diminished markedly. When the work process began into the seawater stream supercharging stage (phase 2#~phase 6#), the seawater stream flowed along the positive direction of the z-axis into the duct. From Figure 8, the large-scale vortex still existed at the duct entrance and exited in the seawater stream supercharging stage, due to the effect of the end cover rotation. Furthermore, it must be noted that the vortex direction was unchanged in the two working stages (phase 1#~phase 6#). These results confirmed the formation of the swirling flow in the duct [16]. Figure 9 shows the different and complex flow structures in instantaneous vorticity contours of RERD. A large-scale vortex with positive vorticity was clearly observed only at z 1 = 1/6 L in the sealed stage (phase 1#). This vortex gradually diminished in phase 2# to phase 4#; then, a vortex with negative vorticity formed in phase 5# and enlarged in phase 6#. It should be pointed out that a large-scale vortex, with an opposite direction of the duct rotation direction, existed at z 1 = 3/6 L. In comparison with stationary ducts in the V-RERD case, the main difference for rotating ducts in  7 Geofluids the RERD case was that there existed the Coriolis force in the duct, affecting the flow structure and motion [18]. As is known to all, the Coriolis force acts in a direction perpendicular to the flow velocity and to the rotation axis. Therefore, the Coriolis force in the RERD case acted the flow rotating in the clockwise direction, which was against the duct rotation direction. As a consequence, the large-scale vortex with a clockwise direction was formed and could be observed at z 1 = 3/6 L in the two working stages as shown in phase 1# to phase 6#.
To characterize the degree of swirl intensity, the dimensionless swirl number, Sn, recommended by Gupta et al. [19] was used, which is defined as where G φ means the axial flux of tangential momentum, G z means the axial flux of axial momentum, and R means  Figure 12: Q criterion isosurface (Q = 0:03) colored by salinity magnitude in RERD. 8 Geofluids the duct radius. The parameters, G φ and G z , can be expressed as where u z and u φ are the axial and tangential components of the velocity, respectively. Figure 10 shows variations of the swirling number Sn in the duct of the RERD and V-RERD. When the duct of the V-RERD started to connect the end cover in phase 2#, the swirl number sharply increased at the duct entrance. That was because fluid had a large value of tangential velocity induced by the end cover rotation. According to the previous study, the value of axial velocity increased from phase 2# to phase 4# [20]. Therefore, the Sn gradually decreased as axial velocity increased at the duct entrance from phase 2# to phase 4#. In the whole seawater stream supercharging stage (phase 2#~phase 6#), the Sn in V-RERD cases kept a positive value in the duct, except at the duct exit in phases 2# and 3# (caused by the backflow). These results confirmed the swirling flow, caused by the end cover rotation, in the duct of the V-RERD. Meanwhile, there was a minimum swirl number in the center of the duct, showing that the rotating end cover had little effect on this area. In the comparison of high swirl intensity at the duct entrance and exit in V-RERD cases, a lower swirl intensity with a negative value was obtained at most part of the duct in RERD cases. In addition, a maximum swirl number with a positive value induced by the Coriolis force can be observed in the center of the duct.
The swirling flow would cause the formation of a vortex, leading to flow instabilities in the duct. The Q criterion [21,22] as a vortex identification method to reveal the complex flow structures was used in this paper, expressed as where Ω is the rotation-rate tensor and S is the strain-rate tensor. Vortex structures were identified and existed in the region with positive values (Q > 0), namely, the rotating force prevailing over the strain force.  Figures 11 and 12 show the visualization of threedimensional vortex structures identified by the Q criterion (Q = 0:03) and colored by salinity magnitude in the V-RERD and RERD, respectively. The vortex formation, expansion, detachment, and diminishment could be observed in the duct of the two devices. In two device cases, the vortex formed and expanded at the duct entrance in the acceleration stage (at the left side of the duct in phase 2#~phase 4#); then, a small-scale vortex detached from it and moved forward in the stabilization stage (in phase 4#~phase 6#) and sealed stage (in phase 7#). Finally, vortices at the duct exit gradually diminished in the brine discharging stage (at the right side of the duct in phase 8#~phase 12#). Although the repeating sequence of the vortex can be observed both in the two devices, there was a discrepancy in the vortex structure. A longitudinal vortex (denoted by the orange arrow in Figure 11) caused by the swirling flow was observed in the V-RERD case, while a transverse vortex (denoted by the orange arrow in Figure 12) was observed in the RERD case. Additionally, the detached vortex moved much farther in the duct of the RERD, leading to a high turbulence level between the two streams. Hence, a high mixing degree could be induced in the working process of the RERD. Figures 13-15 show the salinity distributions and variations in the duct of two devices in one working cycle. The blue and red parts in Figure 13 represented the seawater and brine streams, respectively. During the working process, the mixing layer (the interface of the seawater and brine streams) moved with the inflow movement in the V-RERD case, as shown in Figure 13. Furthermore, its shape was almost unaffected by the detached vortex. However, compared with salinity distributions in V-RERD cases as shown in Figure 13, the mixing layer was twisted and a larger mixing zone can be found in RERD cases. Combined with the results of vortex distributions in Figure 12, the thicker mixing layer in RERD cases in Figure 13 may be induced by the detached vortex, which moved much farther from the duct entrance and stayed a longer time in the center of the duct.

Salinity Distributions.
Owing to the formation of large-scale vortices with a high turbulence level shown in Figure 12, an obvious mixing process can be observed in the RERD working process (see Figures 13 and 15). It should be pointed out that large-scale vortices closed to the duct exit in V-RERD cases (see Figure 12) could induce a higher mass transfer rate and an obvious salinity variation as shown in Figure 15. This phenomenon can be used to explain why the salinity variation in RERD cases had a high fluctuation amplitude in the HP outlet, as shown in Figure 6. Compared with the results in Figure 15, curves of the average salinity in Figure 14 clearly moved with the phase change (inflow movement). Moreover, the curve of the average salinity in V-RERD cases demonstrated a smooth distribution, showing a lower mass transfer rate along the x-axis. It could infer that the vortex formation in RERD cases had a more serious impact on the mixing process than that in V-RERD cases. Therefore, flow control techniques [23,24] should be adopted to suppress the formation of large-scale vortices, finally decreasing the seawater salinity.

Conclusions
In this study, the flow fields in the ducts of RERDs, a typical RERD and a V-RERD, had been investigated by threedimensional simulations. From simulation results, a swirling flow was confirmed in the two devices' dust, which had a significant effect on the fluid dynamics. In the V-RERD case, the swirling flow formation was caused by the end cover rotation, and its direction remained unchanged in the working process. While in the RERD case, due to the effect of the Coriolis force, the swirling flow changed its direction in the center of the duct. To characterize the degree of swirl intensity, the dimensionless parameter of a swirling number, Sn, was adopted. Comparing to the swirling number in the RERD case, the swirling number in the V-RERD case showed a higher magnitude at the ends of the duct, while the lowest magnitude was located in the center of the duct. The threedimensional vortex structures, identified by the Q criterion, were clearly observed in the duct of the two devices. A longitudinal vortex structure was in the V-RERD case, while a transverse vortex structure was in the RERD case. The evolution of vortices showed that the detached vortex had a significant effect on the mixing process. Suppressing the detached vortex formation may be considered an effective way to decrease the volumetric mixing rate. These findings will help the understanding of the mixing process in a RERD and enhance the theoretical basis for the RERD performance improvement.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.