A Water/Ion Separation Device: Theoretical and Numerical Investigation

Featured Application: The proposed application could be exploited for the design of a desalination device at various scales, depending on the desired ﬂow rate of clean water. Abstract: An array of ion separation cells is presented in this work, to propose a novel desalination device. Molecular Dynamics simulations have been incorporated to establish the theoretical background and calculate all parameters that could lead the manufacturing step. The main system component is an ion separation cell, in which water/NaCl solution ﬂows due to an external pressure difference and ions are directed towards the non-permeable walls under the effect of an electric ﬁeld, with direction perpendicular to the ﬂow. Clean water is gathered from the output, while the remaining, high-concentration water/ion solution is re-cycled in the cells. The strength of the electric ﬁeld, cell dimensions, and wall/ﬂuid interactions are investigated over a wide range, and shear viscosity and the volumetric ﬂow rate are calculated for each case.


Introduction
Water flow at the nanoscale has gained increased research interest in nanofluidics, due to its key role in environmental, biological, and chemical processes that can lead the development of novel devices in the neighboring fields of physics, chemistry, and materials science [1]. As the classical Navier theory of slip at the macroscale breaks down, predicting fluid behavior at the nanoscale and understanding the link between molecular structure and macroscopic behavior requires deep investigation over a number of parameters that can affect its flow and transport properties [2]. Only after having established deep understanding at these scales, we will be able to design devices with desired properties for practical applications [3].
The experimental investigation of devices at the nanoscale is uniquely important, as it would not only unravel phenomena taking place there, but also, would offer fabrication guidelines for commercial use. Nevertheless, the whole experimental processes are complex and challenging. Going down to less than 10 nm in pore sizes has only recently been accessible experimentally [4], and the best way to overcome these barriers emerges from numerical approaches [5]. Molecular dynamics (MD) simulations have given a boost in establishing water properties at the nanoscale through numerical calculations, bridging fluid mechanics, statistical mechanics, and condensed matter physics and setting the basis for the design and fabrication of nanofluidic devices [6]. Among these, novel desalination and water purification techniques have been proposed, covering membrane-based applications, reverse osmosis, Capacitive De-Ionization (CDI) and flow-electrode CDI [7][8][9].
Similar to the CDI philosophy, the recently proposed electric ion drift (EID) method does not use membranes or special electrodes [10,11]. The EID method is of the same energy consumption as the capacitive deionization method. This method is approached by the Poisson-Nerst-Planck (PNP) theory, which can be less computationally demanding than MD but, on the other hand, it succeeds only at the bulk and cannot incorporate wall effects 2 of 9 in regions near the channel walls. At the macroscale, the PNP scheme is a popular choice to explain the relation between the electrostatic problem and the mass transport [12,13].
System design, either in real applications or simulation models, has the main role in achieving an efficient desalination process. First, the choice of materials is primarily important. Towards minimization of energy consumption during the desalination process, novel materials have been employed, although sometimes they have been found to have only a slight impact on increasing energy efficiency [14]. On the other hand, carbonbased materials, such as graphene, are among the most popular choices and have been incorporated for the construction of permeable membranes due to their porous structure and wetting characteristics [15]. Recently, synthetic nano-conduits have been introduced, such as artificial water channels (AWCs), to develop highly selective membranes for water desalination [16]. In a neighboring research field with biological interest, bacterial membrane potential dynamics have been investigated for exploitation in microfluidicbased platforms [17].
Another significant parameter that should be borne in mind is the choice of the driving method that removes unwanted substances from the solution. The application of an external electric field, E, parallel or perpendicular to the walls, is common choice. Care has to be taken as the presence of the electric field may alter the hydrogen bond network structure, relaxation times, and rotational dynamics of water [18]. Depending on E direction and strength, water molecules could be subject to phase change [19] but, in general, higher E value leads to higher ion rejection rate [20]. In salt solutions, the electric field along with ion concentration have also a considerable impact on the mutual diffusion coefficients [21][22][23]. In contrast, the water flux is insensitive to the applied voltage and depends on pore size, pressure, thickness, ion concentration, water salinity, and pore geometry [24,25]. Hydrophobic/hydrophilic behavior of the walls is also a matter of investigation, since it could diminish the high flow resistance owing to surface and confinement effects [26]. It has been shown that cotton beads acting as pseudosuperhydrophilic surfaces could enhance clean water flow in a novel, condensation-based desalination system [27,28].
In this work, an array of desalination cells has been incorporated to suggest a novel desalination device. Each cell lies in the order of some nanometers, which is also the typical size-scale for membrane separation techniques for nanofiltration [29], and when the cells are connected in parallel, they comprise a device that could even reach at the macroscale. It is of importance to highlight water/ion behavior in each cell to guide possible fabrication. Towards this direction, in the next sections we present the structure of the proposed model, calculate the volumetric flow rate, and the shear viscosity (in direction parallel and perpendicular to the flow) versus the cell size, wettability, and the strength of the applied electric field. Figure 1 presents the architecture of the proposed desalination device, built by an array of N cells. Each cell resembles a nanochannel, where water molecules and Na + and Cl − ions flow between two planar, periodic, carbon surfaces, theoretically analogous to Poiseuille flow (see the inset of Figure 1). In a previous work, we have investigated the effect of the electric field strength, E, and the distance between the two wall surfaces, h, on desalination performance in a stand-alone nanochannel [30]. In this work, we consider a parallel network of such nanochannels and extend the investigation over the whole network. walls, which is expected to lead Na + near the upper wall, Cl − ions near the lower wall, and clean water in the region around the cell centerline. During the simulations, the solution density is constant, ρ = 1078 kg/m 3 , and temperature remains constant at T = 300 K with the application of Nose-Hoover thermostats.

System Model
Cells of various heights (i.e., the distance between the walls) are simulated and, along with the number of cells, N, inside the device, the flow rate of clean water is estimated in Section 3. Simulation parameters for all atoms in the system are described by the Lennard-Jones (LJ) potential. For two i and j species, the potential is given by [31] where parameters ε and σ are the energy and size parameters, respectively, and rc = 9 Å is the cut-off radius. Moreover, coulombic interaction between hydrogen, oxygen, and ions is given by = 0 , where C the energy conversion constant, qi and qj the charges of interacting atoms, and ε0 is the dielectric constant. For two different species i and j we consider = ( + ) 2 ⁄ , = √( ). When i and j correspond to wall (w) and fluid (f) atoms, the ratio ⁄ represents whether the wall is treated as hydrophobic ( ⁄ = 0.1 and 0.2) or hydrophilic ( ⁄ = 0.5 and 1.0). Water molecules comply to the SPC/E (extended simple point charge) pair potential, which has been found to reproduce adequately the structural and dynamic properties of water [32], especially under confinement [33]. Molecular Dynamics is a common method for interpretation and investigation of nanoscale phenomena. This does not exclude investigation with Navier-Stokes equations or quasi-continuum approaches, provided that one considers the appropriate scale factors and coefficients [34]. For example, shear viscosity across a nanochannel is not constant and this has to be addressed when applying the classical N-S theory [35]. The breakdown Green and red circles are Na + and Cl − ions, respectively. The external electric field E z is applied perpendicular to the piston-driven flow. Clean water is gathered from each cell to the common outlet, while the remaining water/ion solution (brine) is recirculated to further separate water from ions.
In each cell, water/ion solution is inserted on the left, while water and ions are separated at the outlet, as the solution flows from left to right, driven by a piston. This is achieved by the application of an external electric field, E, acting perpendicular to the walls, which is expected to lead Na + near the upper wall, Cl − ions near the lower wall, and clean water in the region around the cell centerline. During the simulations, the solution density is constant, ρ = 1078 kg/m 3 , and temperature remains constant at T = 300 K with the application of Nose-Hoover thermostats.
Cells of various heights (i.e., the distance between the walls) are simulated and, along with the number of cells, N, inside the device, the flow rate of clean water is estimated in Section 3.
Simulation parameters for all atoms in the system are described by the Lennard-Jones (LJ) potential. For two i and j species, the potential is given by [31] where parameters ε and σ are the energy and size parameters, respectively, and r c = 9 Å is the cut-off radius. Moreover, coulombic interaction between hydrogen, oxygen, and ions is given by where C the energy conversion constant, q i and q j the charges of interacting atoms, and ε 0 is the dielectric constant. For two different species i and j we consider σ ij = σ i + σ j /2, ε ij = ε i ε j . When i and j correspond to wall (w) and fluid (f) atoms, the ratio ε w f /ε f f represents whether the wall is treated as hydrophobic (ε w f /ε f f = 0.1 and 0.2) or hydrophilic (ε w f /ε f f = 0.5 and 1.0). Water molecules comply to the SPC/E (extended simple point charge) pair potential, which has been found to reproduce adequately the structural and dynamic properties of water [32], especially under confinement [33].
Molecular Dynamics is a common method for interpretation and investigation of nanoscale phenomena. This does not exclude investigation with Navier-Stokes equations or quasi-continuum approaches, provided that one considers the appropriate scale factors and coefficients [34]. For example, shear viscosity across a nanochannel is not constant and this has to be addressed when applying the classical N-S theory [35]. The breakdown of the continuum hypothesis of the no-slip condition is believed as the main reason that continuum relations have to alter. Moreover, the occurrence of EDLs and specific adsorption effects are to be expected, while possible wall roughness would have a pronounced effect on transport properties [36].

Mathematical Formulation
Shear viscosity, µ, is given by the Stokes-Einstein's relation [37,38] as where constant α = 1.7 Å is the effective hydrodynamic diameter of the water molecule, and k B the Boltzmann constant. The channel diffusion coefficient, D ch , is obtained from the time-averaging of the mean square displacement of N fluid particles [39] as and where r j is the position vector of the jth atom. The brackets indicate time average. We distinguish a parallel, D P , and a perpendicular to the wall direction, D T , component, as where, r xy j is the position vector in xand y-dimensions, and r z j the position vector in z-dimension of the jth atom.
To calculate shear viscosity in direction parallel to the cell walls, µ P , or perpendicular, µ T , we, respectively, substitute in Equation (2) the diffusion coefficients D P and D T . Diffusion coefficient calculations from the above equations, characterizing the mechanism of mass transfer at the lower scales, improve our understanding on the microscopic effects taking place at these scales. This would be of interest for problems that have to do with electric fields, since the system behavior is investigated at the lower possible level and all hidden effects are revealed [40]. Equation (2) provides calculations on shear viscosity, bypassing a significant amount of time-dependent statistical mechanics simulations needed to extract shear viscosity by other methods, such as the Green-Kubo formalism, and has been widely considered in nanoscale systems [37,38]. It shows that the size effect of D ch is directly influenced by the size effect of the shear viscosity µ. As D ch is decreasing with the channel width, µ is increasing. Shear viscosity is connected to friction, i.e., a macroscopic property, and it becomes clear that it also affects the ion and clean water flow rate for the proposed device.
In order to argue on the achievable performance of the proposed desalination method, the volumetric flow rate in every channel is calculated, as where A is the cross-sectional channel area, h o the outlet height, υ the average fluid velocity and L y is the channel y-dimension. To disjoin calculations from L y , since the simulation Appl. Sci. 2021, 11, 8548 5 of 9 system is periodic in the y-direction, we consider the more general volumetric flow rate per unit length to be

Shear Viscosity Calculations
Shear viscosity is extracted from the diffusion coefficients according to Equation (2). By separating the calculations in a parallel and a perpendicular direction, µ P and µ T , respectively, we obtain a close-up view of the internal mechanism that affects ion movement inside each cell. For the narrowest cell under investigation (h = 3 nm), in Figure 2a, we observe that there exists increased shear viscosity perpendicular to the wall, µ T , significantly higher that the parallel component, µ P . As the strength of E increases from E = 0.0 to 1.0 V/Å, µ T decreases monotonically. In contrast, µ P remains unaffected by the electric field value. The channel shear viscosity is, though, unaffected, remaining around the water bulk value for the range of E investigated here. In a weaker E range, the viscosity component parallel to the electric field has been found to increase monotonically with the electric field strength, while the perpendicular component first decreases and then increases with E [41]. We also note that experimental results do not always adhere to simulation, and this is attributed to the lower E strength applied, compared to simulated nanoscale systems [42].
where A is the cross-sectional channel area, ℎ the outlet height, ̄ the average fluid velocity and Ly is the channel y-dimension. To disjoin calculations from Ly, since the simulation system is periodic in the y-direction, we consider the more general volumetric flow rate per unit length to be = =̄× ℎ (8)

Shear Viscosity Calculations
Shear viscosity is extracted from the diffusion coefficients according to Equation (2). By separating the calculations in a parallel and a perpendicular direction, μP and μT, respectively, we obtain a close-up view of the internal mechanism that affects ion movement inside each cell. For the narrowest cell under investigation (h = 3 nm), in Figure 2a, we observe that there exists increased shear viscosity perpendicular to the wall, μT, significantly higher that the parallel component, μP. As the strength of E increases from E = 0.0 to 1.0 V/Å, μT decreases monotonically. In contrast, μP remains unaffected by the electric field value. The channel shear viscosity is, though, unaffected, remaining around the water bulk value for the range of E investigated here. In a weaker E range, the viscosity component parallel to the electric field has been found to increase monotonically with the electric field strength, while the perpendicular component first decreases and then increases with E [41]. We also note that experimental results do not always adhere to simulation, and this is attributed to the lower E strength applied, compared to simulated nanoscale systems [42].
For the h = 6 nm cell, where the effect of wall confinement is weaker on the transport properties of fluids [43], we observe that shear viscosity anisotropy is smaller compared to h = 3 nm, and its behavior is almost isotropic at E = 1.0 V/Å. Shear viscosity is also examined on its dependence on wall/fluid interaction ⁄ in Figure 3, for cells of height h = 3-15 nm. For h = 3 nm (Figure 3a), all shear viscosity components acquire their maximum value for the strongly hydrophilic wall ( ⁄ = 1.0). However, this effect is weak and does not have a strong impact on μ values, and this is also the case for h = 6, 9 and 15 nm in Figure 3b-d, respectively. For the h = 6 nm cell, where the effect of wall confinement is weaker on the transport properties of fluids [43], we observe that shear viscosity anisotropy is smaller compared to h = 3 nm, and its behavior is almost isotropic at E = 1.0 V/Å.
Shear viscosity is also examined on its dependence on wall/fluid interaction ε w f /ε f f in Figure 3, for cells of height h = 3-15 nm. For h = 3 nm (Figure 3a), all shear viscosity components acquire their maximum value for the strongly hydrophilic wall (ε w f /ε f f = 1.0). However, this effect is weak and does not have a strong impact on µ values, and this is also the case for h = 6, 9 and 15 nm in Figure 3b-d, respectively.
From a macroscopic point of view, studied under the Poisson-Boltzmann (PB) theory, it has been shown that the response of an electrolyte to an applied voltage leads to sharp ion concentration near a surface [21,22]. This has major implications for the system dynamics, as changes in diffusion and viscosity have been observed. In fact, the viscosity, which increases in high-concentration solutions, dominates the concentration dependence of D for some systems. Nevertheless, PB theory has many shortcomings because it neglects ion-ion interactions and steric effects and it is proposed that the electrolyte dynamics are better described by a modified form of the PNP equations [12,23]. From a macroscopic point of view, studied under the Poisson-Boltzmann (PB) theory, it has been shown that the response of an electrolyte to an applied voltage leads to sharp ion concentration near a surface [21,22]. This has major implications for the system dynamics, as changes in diffusion and viscosity have been observed. In fact, the viscosity, which increases in high-concentration solutions, dominates the concentration dependence of D for some systems. Nevertheless, PB theory has many shortcomings because it neglects ion-ion interactions and steric effects and it is proposed that the electrolyte dynamics are better described by a modified form of the PNP equations [12,23].

Volumetric Flow Rate
The volumetric flow rate per unit length, q, for each cell is shown in Figure 4 ⁄ and for all cell heights investigated here, q obtains its maximum value. This is attributed to the existence of slip near the solid boundaries, which facilitates flow and, in the case of narrower channels, is more pronounced [26]. As the cell height increases, e.g., for h = 21 nm, the flow rate is stabilized to a constant value within statistical accuracy, without being further affected by the wall/fluid interaction.

Volumetric Flow Rate
The volumetric flow rate per unit length, q, for each cell is shown in Figure 4, calculated over the range 0.1 ≤ ε w f /ε f f ≤ 1.0. For the h = 3 nm cell, q presents a monotonic, decreasing behavior as the walls vary from being strongly hydrophobic (ε w f /ε f f = 0.1) to strongly hydrophilic (ε w f /ε f f = 1.0). Moreover, for ε w f /ε f f = 0.1 and for all cell heights investigated here, q obtains its maximum value. This is attributed to the existence of slip near the solid boundaries, which facilitates flow and, in the case of narrower channels, is more pronounced [26]. As the cell height increases, e.g., for h = 21 nm, the flow rate is stabilized to a constant value within statistical accuracy, without being further affected by the wall/fluid interaction.

Configuration Issues
We now turn our attention to configuration issues and how the previous theoretical analysis can be transformed to real applications, such as the ion separation device presented in Figure 1. Apart from piston and water/ion feeding issues, which are not examined in this work, characteristics of the cell array have to be pointed out, such as the cell dimensions, the number of cells, and, consequently, the device dimensions, to be assessed over the induced flow rate of pure water from the output channel. All these parameters

Configuration Issues
We now turn our attention to configuration issues and how the previous theoretical analysis can be transformed to real applications, such as the ion separation device presented Appl. Sci. 2021, 11, 8548 7 of 9 in Figure 1. Apart from piston and water/ion feeding issues, which are not examined in this work, characteristics of the cell array have to be pointed out, such as the cell dimensions, the number of cells, and, consequently, the device dimensions, to be assessed over the induced flow rate of pure water from the output channel. All these parameters are presented in Table 1. At first, the minimum channel length L x for a cell of height h that could yield satisfying ion separation results has been estimated in a previous work and is given by L x = 320h − 27.46 [30]. For the device dimension on z-direction equal to L z = 1 cm, the number of cells, N, in the array can be calculated, as well as the resulting volumetric flow rate q (per unit length). Depending on the configuration capability, one may choose to construct an array of 395,257 cells in a (L x , L z ) = (1 cm, 7.33×10 −4 cm) where q = 1.54×10 −7 lt/day or an array of only 80 cells with macroscale dimensions (L x , L z ) = (1 cm, 3.64 cm) where q = 3.76 lt/day.
The proposed ion separation device could be incorporated to construct small/mediumsized systems and cannot be considered as a replacement for well-established, industrial water desalination plants. However, this approach seems to have advantages over classical methods, and further manufacturing issues are to be explored in the future.

Conclusions
A novel desalination procedure has been presented in this work, consisting of an array of cells that can function in parallel, based on simulation data derived from molecular dynamics simulations. The investigation has been based on nanoscale simulations to extract each cell property and project the findings on a wider scale.
Each cell comprises a nanochannel filled with water/ion solution. A perpendicular to the walls electric field is applied to drive ions towards the walls, while an external force with direction parallel to the walls is applied to induce flow. Clean water is gathered at the outlet. Nanoscale surface effects are expected to affect water/ion properties inside the cell. In terms of shear viscosity, calculations have shown that there exists increased shear viscosity perpendicular to the wall, which decreases monotonically as the strength of the external electric field E increases from E = 0.0 to 1.0 V/Å. On the other hand, the shear viscosity component parallel to the flow direction remains unaffected. Wall hydrophobicity/hydrophilicity does not seem to affect shear viscosity significantly, at least at the range investigated here.
Moreover, calculations on the volumetric flow rate per unit length for each cell has revealed a monotonic, decreasing effect behavior as the walls turn from hydrophobic to hydrophilic. This effect is more pronounced in the smaller cell dimensions studied, for h = 3 nm. However, for every cell dimension shown here, the flow rate obtains its maximum value for the strongly hydrophobic channel.
By establishing a parallel array of cells investigated in this work, it is possible to achieve practical water desalination rates, without the need to experiment on expensive membrane materials and/or electrodes. Some indicative cases have been tabulated and we believe that this will assist on establishing a practical desalination device in the future.