Characterization of exchange flows within pressure exchanger device by Lattice Boltzmann method

Seawater reverse osmosis (SWRO) with pressure exchanger (PX) been used widely on the earth for desalination to produce fresh water. Mixing due to exchange flow is considered to be dominant. In this paper, full-depth exchange flow is numerically simulated by Palabos, an open source program based on Lattice Boltzmann Method. GIF pictures and data produced by Palabos are analyzed to figure out how exchange flow processes. The velocity of the current front U=0.45 √(g(1 − γ)H is obtained from analysis. The results from the simulation are consistent with the Benjamin’s energy-conserving theory very well for full-depth exchange flow.


Introduction
Seawater reverse osmosis (SWRO) [1,2] process is an economic technology for large-scale seawater desalination processes. PX (pressure exchanger) is an energy recovery system for SWRO. High-pressure brine and low-pressure seawater contact in PX to achieve energy recovery.
Dispersion is mentioned to basic concept of mixing in PX. There are two laws in the three fundamental transport processed. Diffusion is a basic molecular motion which is a random motion in all directions. Fick's Law [3] can be used as a theory to explain diffusion process. The second law is the law of conservation of mass within a control volume with a control surface boundary. The1-D dispersion problem with intermittent oscillatory pressure forcing in PX has predicted the time for the mixing zone of seawater and brine reject to reach the ends of the duct [4]. The speed of the front of the current was observed to be nearly constant and the interface shape kept typical with a raised head and nearly horizontal interface behind the head [5]. A model has been formed [6] for steadily gravity currents in infinite ambient fluid with two deductions based on inviscid Bernoulli equation that are an angle 60°between the front shape and the bottom boundary and constant velocity in the light fluid. So for currents propagate at less than one half depths, there were energy losses [7] by energy dissipation and the dimensionless energy flux entering from upstream [8]. It was pointed out that the energy dissipation could be described by an equivalent head loss independent of height which are also unrealistic [9]. Additionally, it also be thought h/H > 0.347 was not able to be achieved with the Froude number F H in infinite ambient fluid as √2 [9]. Two different cases of lock exchange flow in high Reynolds numbers have been carried out with ignoring viscous effects [8]. After the lock is removed, lock exchange flow immediately forms with the dense fluid propagating. The motion of the following current can be obviously observed to have three phases [8,10,11]. For this study, exchange flow is greatly more important than the mixing and is dominant in PX devices. Therefore, a full-depth release exchange flow will be simulated and researched in this paper.

The Lattice Boltzmann Method
Palabos, an open-source program, is the numerical simulation tool for CFD (Computational fluid dynamics). Lock exchange flow is simulated by Palabos in this project. Palabos provides a library of a framework in C++ for computational fluid dynamics based on lattice Boltzmann method (LBM), which has an appropriate balance between generality, ease of use, and numerical efficiency.
The aim of Palabos is to set up lattice Boltzmann models for fluid flow simulation. The LBM is the solution of nonlinear partial differential equations such as Navier-stokes equation, which can be applied to simulate single and multi-phase flows of single and multi-component fluids [12]. The lattice Boltzmann method is a kinetic-based approach for fluid flow simulations with microscopic models [13], which is different from previous conventional numerical ways of discretizing macroscopic continuum equations. Macroscopic properties which comply with the expected macroscopic equations can be obtained by simplified kinetic models that follow the essential mechanics of microscopic processes. It is introduced that the macroscopic dynamics of fluid in these simplified kinetic-type models for fluid flows results from the collective behavior of many microscopic particles [14], which is required to avoid solving complicated kinetic equations and following each particle as in molecular dynamics simulations to get simplification.
LBM is different from some other computational fluid methods in many aspects. Firstly, LBM is derived from the discrete velocity model with a group of first-order partial differential equations rather than second-order. Secondly, the convection operator or streaming process of the LBM is linear. Thirdly, the LBM is kinetic-based. Fourthly, the incompressible Navier-Stokes equations can be obtained in the nearly incompressible limit of the LBM. The data communication of LBM is always local and the pressure of the LBM is obtained from an equation of state. Meanwhile, in other numerical simulations of the incompressible NS equations, the pressure satisfies a Poisson equation with required velocity strains, which often causes numerical difficulties.

The Lattice Boltzmann Equation
The LBM is derived from lattice gas automata (LGA) [15], a discrete particle kinetics with a discrete lattice and discrete time. The LBM can also be viewed as a special finite difference scheme for the kinetic equation of the discrete-velocity distribution function. The lattice Boltzmann method adopts a simplified, fictitious molecular dynamic in which space, time, and particle velocities are all discrete. A discrete kinetic equation for the particle distribution function can be expressed as [16], Where fi is the particle velocity distribution function along the ith direction; , is the collision operator which represents the rate of change of fi resulting from collision and depends only on the local distribution function. ∆t and ∆x are time and space increments, respectively.
The density and momentum density u is defined as particle velocity moments of the distribution function, fi, = ∑ , u= ∑ . Ω satisfies conservation of total mass and total momentum at each lattice: ∑ 0 ,: ∑ 0 . If only the physics in the long-wave-length and low-frequency limit are of interest, the lattice spacing ∆x and the time increment ∆t in Equation 1 can be regarded as small parameters of the same order . Carrying out a Taylor expansion in time and space, the continuum form of the kinetic equation accurate to second order ε can be obtained: 2 2 (0.5) : : A particular simple linearized version of the collision operator makes use of a relaxation time towards the local equilibrium using a single time relaxation [13]. The relaxation term is the Bhatnagar-Gross-Krook (BGK) collision operator. The lattice BGK model makes the computations more efficient and allows flexibility of the transport coefficients. So LBGK equation is: There are two steps in the LBM in LBGK model: Collision step : Streaming step 1) : (, Where represents the post-collision state. Here, make ∆t as 1 and ∆x as 1.

D2Q9 Model
D2Q9 model is a kind of LBGK models, which represents two dimensions and nine velocities. A particle can move to directions, shown in Figure 1(0 means staying at original point) [12]. From the LBGK equation, Collision of the fluid particles is considered as relaxation towards local equilibrium and the equilibrium distribution for this model are defined as, Where i with 4/9 for i=0, 1/9 for i=1, 2, 3, 4 and 1/36 for i=5, 6, 7, 8.

Shan & Chen Model
It has been proved that Shan & Chen model is powerful for simulation of multi-phase and multicomponent flows [13], in which microscopic interactions are used to modify the surface-tension-related collision operator so that the interface operation can be automatically processed as above said to conduct collision and stream step. In this model, the collision operator does not satisfy local momentum conservation, which can be reasonable because interactions between phases change the point-wise local momenta at the positions. However, it is also a limitation of this model. For simplicity, only the nearestneighbor interactions were involved in the Shan & Chen model.

No-slip Boundary Conditions
No-slip boundary conditions are applied to the viscous flow, which means that there is no shear velocity on the boundary for flows. In the LBM, a particle distribution bounce-back scheme for a solid wall is used to realize no-slip boundary conditions. With some nodes near the boundary, the neighboring nodes outside the boundary cannot communicate with inside nodes. When particles strike into the boundary wall, the wall will reflect particles and just reverse the distribution direction without changing the value of velocity distribution instead of collision step in bounce-back boundary condition.

Simulation of Exchange Flow with Densities of 0 and 1
At the beginning of simulation, it is easy to simulate full-depth lock exchange flow with densities of 0 and 1. Therefore, simulation of exchange flow with densities of 0 and 1 is the first step. The initial condition of simulation model is set up with heavy liquid (density 1) on the right side as a red block and light liquid (density 0) on the left side as a blue block, as shown in Figure 2.

Simulation of Exchange Flow with Practical Densities
To improve the practical significance, the next step to do is to change the densities into the more real ones that are 1046 kg/m 3 and 1028 kg/m 3 . So, with the program of 0 and 1, the program has the same structure and the same implementation with initial condition. Time interval of simulation is 100 in the program. Choice of lattice dimensions is 1000 of x axis and 200 of y axis. The shapes of flows are almost same with ones with densities 0 and1. But the difference between them is the velocity of the leading part of the current, which is that the exchange flow with density 1.046 and 1.028 travels faster. From Figure  4, the location of front of density 0 and 1 is 318 at time step 1400 while for practical density is 311. The front of dense liquid reach to the left end at time step 4000 for density of 0 and 1 on the other hand it is at time step 3900 for density of 1.046 and 1.028. As a result, it can be deduced that the density difference can affect the process of exchange flow, especially the velocity of the front. And the higher the differential of densities is, the faster the current front.  Parts of the results of the practical simulation are shown in Figure 5.  These pictures show that the constant-velocity phase of exchange flow. The speed of the front is constant. And the height of the front is seen to be equivalent. From Benjamin's energy-conserving theory, the energy-conserving depth of full-depth is the line where h=1/2H (the black horizontal line). The dense flow and the light flow occupy nearly one half of the channel, which means that the energy of the exchange flow is conserved.
The exchange flow is almost symmetrical about the black centerline so the characteristic of the two counter-flow currents should be the same. Besides the most important is to research the time when brine reject reaches the end. Therefore, the dense flow will be analyzed.

Front Location
With the GIF pictures produced by numerical models, numerical analysis can be done on the front of exchange flow. At first, the units of data from the simulation are different from the physical meaning because in Palabos the physical model is discretized into single microscopic particles with lattice Boltzmann method. So, it is required to know how to convert the units between LBM models and physical models. A method including two steps [17]. First, a data conversion between a physical model and a dimensionless model is made. It is to make units independent on the physical scale. Second, a data conversion between a physical model and a discrete model is made. So, the dimensionless model is the link of the physical and discrete one. For Palabos, the dimensionless parameter is Reynolds number (Re) and Froude number (F H ), which are the same in the three models. The conversion is shown in Figure 6. The unit conversion in this simulation by Palabos is listed in Table1.  After conversion of units, the next step is to analyze the front speed. As above results, references and location of the front is marked, with lattice units, as in the figure. Choose the time interval 300, and t and l (the distance the front travels) will be non-dimensionalized by the channel depth into dimensionless parameters with: Where γ is density ratio of 1028/1046 and H=200. These data are listed in Table 2 as followed. With the handling data, the curve of dimensionless parameter l/H and t 1 / can be plotted as Figure 7, The curve shows relationship between the dimensionless front location and time. The front location l/H is linear to the time t 1 / and so that the velocity of the front of dense current is constant. This characteristic, the front with constant velocity, is complied with Benjamin's energy-conserving theory.
In the previous research, Shin have done a same data analysis of relationship between front location and time on the experiment of full-depth lock exchange flow and obtained the curve of full-depth exchange flow. Compared with the curve of Shin's, the relationship is the same complied with the approximate linear relation. Additionally, it is easy to observe that the slope of Figure 7 with γ density ratio of 0.983 (1028/1046) and the slope of Shin's with γ density ratio of 0.993 are almost the same. So, with the same slope, the velocity of the front of exchange flow depends on H the channel depth for Boussinesq flow (the density ratio γ nearly to 1). It also indicates that for Boussinesq flow the differential of densities has little effect to the speed of the current front. However, the exchange flow with density 0 and 1 is not Boussinesq flow so the density difference can affect the flow. Back to Figure 7, the slope is approximate 0.45. So Then with transposition and division out √H, get: Because the speed of the front is considered to be constant in the simulation, l/t can become U, which means the front speed. So, the expression of the front speed U can be obtained: On the other hand, Froude number F H is Here the g will be converted to g': ' (1 ) Where = ρ 1 /ρ 2 ≈ 1 for Boussinesq flow, ρ 1 <ρ 2 . So, FH is 0.45 in this simulation by Palabos with density 1046 kg/m 3 and 1028 kg/m 3 . In Benjamin's energy-conserving theory, the FH is predicted to be 0.5 for full-depth lock exchange flow presented by Shin et al. [8] and the simulation is consistent with Benjamin's energy-conserving theory. Therefore, Benjamin's energy-conserving theory can predict fulldepth lock exchange flow very well.
For exchange flow in PX devices, setting a duct with length (L) of 1m and radius (r) of 0.01m rotating around a centre axis with r/R=0.1 (R is rotation radius), the time for brine reject reaching the end can be calculated. Densities are assumed to be 1046 kg/m 3 and 1028 kg/m 3 and density ratio is 1028/1046= 0.983. Then centrifugal acceleration is obtained a=r 2 = 100 m 2 /s. So according to the simulation, the front speed is: Replace g with Where a is 100, is 0.983, H is 0.01, so U is 0.0587 m/s 2 . So = / = 1 /0 .0 5 8 7 = 1 7 s T L U (15) Therefore, the reach time for brine reject in 1m duct is 17s from the simulation.

Local Particle Velocity
A series of local particle velocity vector can be obtained from the particle velocity value in the results of simulation. Although the velocity of the front is constant in the macroscopic level, the process of local particle can be observed in the microscopic level from results of Figure 8. The velocities of particles of the front can be observed to decrease and has a trend towards bottom boundary. The decrease is caused by the counter-force in the contact interface from the light liquid. The trend of particles towards bottom boundary is complied with the macro behavior of dense flow that travels along the bottom boundary.
At the gate location, the horizontal velocity of particles at the same point is observed to keep constant with time changing. But the direction of vertical velocity changes, which is considered to be due to the bounce-back boundary condition and not a characteristic of exchange flow.
Velocities of particles of the part between the front and the gate location, a rectangular area, keep constant at the same point as well. As a consequence, the flow between the front and the gate location acts like uniform flow. In addition, the velocities of particles near the boundary are much smaller than those away from the boundary, which indicates the no-slip boundary conditions.

Conclusion
In this paper the full-depth exchange flow by numerical simulation has been studied and the velocity of the current front is predicted so the efficiency of PX can be improved by reducing the dense current velocity. The simulation with Palabos adopts lattice Boltzmann method to build a model through discretizing flows into discrete lattices with discrete time and space rather than solve the Navier-Stokes equation. The model is about Boussinesq flow. The following points can be outlined as outcomes of the study: (1) From the results of the numerical simulation, full-depth lock exchange flow is symmetrical about the centerline and each current (dense and light current) occupies half of channel along its own boundary. Furthermore, a curve of the front locations and time, drawn through dimensionless variables from the simulation, is consistent with the curve that Shin et al. [8] has obtained.
(2) Froude number is 0.45 while the Benjamin's energy-conserving theory can get this number of 0.5 close to 0.45 for full-depth lock exchange flow. The results from the simulation are consistent with the Benjamin's energy-conserving theory very well. Therefore, the Benjamin's energy-conserving theory can explain full-depth lock exchange flow with low density difference. (4) The time needs to be more for dense current reaching the end to improve the efficiency of PX so that the length of duct needs to be bigger or U needs to be smaller. The gyration radius, the diameter of ducts and the angular velocity need to small. On the other hand, if values of these parameters are reduced the energy recovery efficiency can be reduced. So, adjustment of these parameters is required to balance the efficiency to improve the PX devices and achieve better performance of PX.