Next Article in Journal
Impacts of Climate Change and Land-Use Change on Hydrological Extremes in the Jinsha River Basin
Previous Article in Journal
Dissolved Carbon Transport and Processing in North America’s Largest Swamp River Entering the Northern Gulf of Mexico
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Investigation of the Effect of Two-Dimensional Surface Waviness on the Current Density of Ion-Selective Membranes for Electrodialysis

1
Mechanical and Aerospace Engineering Department, New Mexico State University, Las Cruces, NM 88003, USA
2
Civil Engineering Department, New Mexico State University, Las Cruces, NM 88003, USA
*
Author to whom correspondence should be addressed.
Water 2019, 11(7), 1397; https://doi.org/10.3390/w11071397
Submission received: 13 June 2019 / Revised: 4 July 2019 / Accepted: 5 July 2019 / Published: 7 July 2019
(This article belongs to the Section Wastewater Treatment and Reuse)

Abstract

:
Ion-selective membranes are an important component of electrodialysis stacks for desalination. Manufacturing imperfections or slight inhomogeneity of the material can lead to minute membrane surface imperfections. Two-dimensional solutions of the coupled Poisson–Nernst–Planck and Navier–Stokes equations were sought for a perfectly smooth membrane and for membranes with well-defined small-amplitude harmonic surface roughness. The simulations were carried out with the validated rheoEFoam solver by Pimenta and Alves. In the overlimiting regime, the electric field is strong enough for an electrokinetic instability to occur. The instability leads to disturbance growth and the formation of electro-convection cells, which strongly increase the current density. The present simulations show that with an increasing ion concentration and applied voltage, the instability becomes stronger and the overlimiting regime is reached earlier. The limiting current density shows a noticeable dependence on the wavelength of the surface roughness. When the wavelength of the surface roughness is incommensurate with the wavelength of the naturally occurring instability, the limiting current density is increased. Since production membranes will always have some degree of surface roughness, this suggests that membrane surface treatments which favor certain wavelengths may have an effect on the overall membrane performance.

1. Introduction

Desalination of alternative waters, such as brackish water, seawater, municipal, and industrial wastewater, has become a critical strategy to expand traditional water supplies and to alleviate water shortages [1]. Electrodialysis (ED) is a membrane desalination technology that uses semi-permeable ion-exchange membranes (IEMs) to selectively separate salt ions in water under the influence of an electric field [2]. An ED stack consists of pairs of anion-exchange membranes (AEMs) and cation-exchange membranes (CEMs) arranged alternatingly between an anode and a cathode (Figure 1). The positively charged cations migrate toward the cathode, pass the negatively charged CEM, and are rejected by the positively charged AEM. The opposite occurs when the negatively charged anions migrate toward the anode. This results in an alternatively increasing ion concentration in one compartment (concentrate) and decreasing salt concentration in the other (diluate). The process can be viewed as a dialysis process, where the ion diffusion across the membranes is amplified and orientated by an electric field. Electrodialysis reversal (EDR) is a modified ED process, where the electrical polarity of the electrodes is reversed periodically. This results in a reversal in the direction of the ion transport so that the concentrated stream becomes the diluate stream and vice versa. This self-flushing of scale forming ions and fouling matter in concentrate compartments allows EDR to operate at a higher water recovery than ED.
Electrodialysis and EDR have been utilized for decades to desalinate brackish water, treat municipal and industrial wastewater, and produce NaCl from seawater [1]. It has also found applications in chemical processes, as well as the food, beverage, and drug industries. Despite the broad application of ED/EDR, the number of experiments and simulations concerned with the thermodynamics and fluid dynamics of the ED process is very limited [3]. For example, a key parameter for the design and operation of ED is the limiting current density (LCD). It is typically identified as the diffusion-limited current density, corresponding to the complete solute depletion in the layer adjacent to the membrane surface [2]. The limiting regime in ED is associated with the appearance of electroconvective structures near the membranes [3]. The structures are caused by an electrokinetic instability (EKI). The phenomena involved in the limiting region are not yet fully understood. Any improvement of the effectiveness of the ED process based on a reduction of the energy losses (e.g., Chehayeb and Lienhard [3]) and methods that overcome inherent operating limitations have strong implications for environmental sustainability and the water, food, and chemical industries.
This paper is concerned with the factors governing the onset of the overlimiting regime in ED. Of particular interest are the effects of the applied voltage, salt concentration, and membrane surface geometry on the onset of the EKI. Towards that end, simulations of the phenomena near a CEM were carried out for a binary electrolyte solution. Because the hydrodynamics and ion transport for CEMs and AEMs are similar, the present findings transfer over to AEMs. The paper starts out with a discussion of the onset of the overlimiting regime and how it is related to the EKI. This section is followed by a discussion of the governing equations, the setup of the simulations, and the defining parameters. Instantaneous visualizations and time histories of the current density are employed to investigate the onset of the EKI. The paper concludes with a brief discussion of the results and their relevance for ED.

2. Overlimiting Regime and Electrokinetic Instability

When the voltage is low enough, the cations in the diffusion layer are drawn to the CEM surface by both migration and diffusion. For the anions, the diffusive and electromigration fluxes cancel each other out. The solution remains locally charge neutral with the exception of a thin electric double layer (EDL) with a thickness of tens of micrometers that forms on the selective membrane surface. As the voltage and hence the concentration gradients are increased, the electroneutrality condition demands very small anion concentrations near the CEM surface. For some limiting voltage, the electroneutrality condition in the diffusion layer is violated and another layer emerges between the EDL and the diffusion layer, which is known as extended space charge layer. Beyond a critical voltage, an instability sets in (here referred to as electrokinetic instability, EKI) and disturbance amplification results in the formation of organized flow structures and electroconvection. The flow structures eject fluid and charge density away from the membrane and substantially increase the current density. For mono-valent binary solutions, the magnitude of the limiting current density is inversely proportional to the boundary layer thickness [4]. The electroconvective structures are similar to those arising from the Rayleigh–Bénard instability, which is driven by opposing gravity and buoyancy forces.
Based on experiments, Rubinstein et al. [5] made a connection between the onset of the overlimiting regime, which is associated with an unsteadiness of the current density, and the occurrence of the EKI. Later, simulations and stability investigations followed (e.g., Demekhin et al. [6], Rubinstein and Zaltzman [7,8]). A review on electroconvective micro-structures and related phenomena was provided by Chang et al. [9]. Despite the small length-scales associated with electroconvection, the instability can lead to broad spectra and an energy cascade similar to those observed in turbulent flows. Dye visualizations by Yossifon et al. [10] revealed vortices with an initial wavelength of 10 to 100   μ m (close to the thickness of the boundary layer). As the vortices grow in size they merge and the wavelength is reduced. For narrow microchannels, the channel height determines the final vortex size [10]. When fully developed, the size of the structures is of the order of the distance between the membranes, which is typically less than 1   mm . For very narrow channels, the vortices never form, which, of course, suggests that the instability is suppressed.
The simulation of reverse osmosis (RO) and even more so ED processes is challenging because of the complex flow physics and the large spread of length-scales from nanometer-sized pores to millimeter-sized spacers. For liquids, the continuum assumption breaks down for length scales around 10   nm . Although some researchers argue that lattice Boltzmann methods should be employed to describe the flow in ED stacks (because of the small scales), the general consensus is that the continuous flow model (Navier–Stokes equations) is accurate enough. Aside from the Navier–Stokes (N–S) equations, for ED processes additional conservation equations for the molar concentrations, c i , of the ions, the Nernst–Planck (N–P) equations, have to be solved. Diffusion and electromigration dominate the ionic transport in the boundary layer. Since the water flux through the IEMs in ED is typically insignificant, the convection terms are often neglected. For example, Kim et al. [11] developed a mathematical model for the steady-state transport of three ion species through an ED cell based on diffusion and migration.
Research aimed at understanding the physics of the boundary layer near the IEMs has to consider the complete N–P equations. Druzgalski et al. [12] solved the coupled two-dimensional (2-D) N–P and N–S equations to investigate the EKI for a binary electrolyte between an IEM and a reservoir (water flow channel). They also presented a framework for the development of ensemble-averaged models (similar to Reynolds-averaged Navier–Stokes models) for the inclusion of the resulting electroconvective-driven mixing processes. Similar simulations were performed by Zourmand et al. [13] and Enciso et al. [14]. Enciso et al. [14] performed experiments and found a 16% difference between the computed and measured mass transfer coefficients. Pimenta and Alves [15] developed an OpenFOAM solver, rheoEFoam, which is part of the rheoTool suite, for the simulation of electrically-driven flows. Importance was placed on the conservation of the ionic species, numerical robustness, and an overall second-order-accuracy of the discretization. RheoEFoam was validated by comparison with the binary electrolyte problem by Druzgalski et al. [12].
Demekhin et al. [16] carried out three-dimensional (3-D) simulations of a reservoir-membrane configuration and found three typical flow structures arising from the EKI: Two-dimensional (2-D) counter-rotating rolls, regular hexagonal (“bee hive”) structures, and more chaotic 3-D structures. The onset of the instability and the type of the resulting flow structures were found to depend on the applied voltage and a coupling coefficient, κ = ε ψ 2 / ( η D ) , between the hydrodynamics and the electrostatics. Here, ε is the permittivity, ψ is the electric potential, η is the dynamic viscosity, and D is the ionic diffusivity. When the voltage exceeded a first critical value, counter-rotating rolls formed. When the voltage was increased further beyond a second critical value, hexagonal structures arose. Finally, for very large voltages, the flow became chaotic and the spectra filled up. Druzgalski and Mani [17] employed 3-D simulations to investigate the fully chaotic regime. They found short-lived high-current-density spots that appeared randomly on the membrane surface and contributed significantly to the mean current density.
Rubinstein and Zaltzman [8] showed that the EKI leads to the growth of disturbances with a preferred natural wavelength whose initial amplitude can be increased by imparting a minute waviness on the membrane surface. This paper is concerned with the dependence of the current density on the wavelength of the waviness of the membrane surface. The present effort has been motivated by research in fluid dynamics that has shown that the growth of instability waves can be influenced. Saric et al. [18] investigated crossflow transition on swept wings. By adding micro-sized roughness elements that were spaced at a wavelength below the most unstable (natural) wavelength of the instability, the growth rate of the steady crossflow mode could be reduced and transition to turbulence could be delayed. Embacher and Fasel [19] stabilized a secondary instability by forcing the primary instability at a wavelength that was incommensurate with the most amplified natural wavelength. Finally, streamwise riblets with a spanwise wavelength close to the wavelength of the naturally occurring near-wall streaks in turbulent boundary layers (the wavelength can be as small as 10   µ m ) were found to reduce the turbulent skin-friction drag [20].
This paper shows that by transferring these ideas to ED applications, the onset of the EKI and the overlimiting regime can be delayed. Of particular relevance for the present investigations are the results for wavy membranes by Rubinstein and Zaltzman [8] and the crossflow instability research by Saric et al. [18]. As will be shown later, the required dimensions of the membrane surface structures will be much smaller than what has already been considered in the literature (e.g., Güler et al. [21]). All of the present simulations were carried out with the rheoEFoam code, which was developed by Pimenta and Alves [15]. The majority of the simulations are 2-D. The concentration and the strength of the electric field were varied. For one case, a 3-D simulation was carried out to demonstrate the relevance of the 2-D simulations.

3. Methodology

3.1. Governing Equations

The equations governing the unsteady, incompressible, isothermal, single-phase, laminar flow of a Newtonian fluid under the action of an electric field were solved with rheoEFoam by Pimenta and Alves [15]. Within rheoEFoam, after the initialization of the field variables, a Poisson equation for the electric potential, ψ , with unit V = N m / C , is solved:
× ( ε ψ ) =   ρ E ,
where ε is the electric permittivity in C 2 / ( Nm 2 ) and ρ E = F z i c i is the charge density in C / m 3 . Here, F = 96 , 485   C / mol is Faraday’s constant, and z i and c i are the charge valences or charge numbers and molar concentrations (in mol / m 3 ) of species i . Updated molar concentrations, c i , are obtained from the Nernst–Planck equation:
c i t + v × c i =   × ( D i c i ) +   × [ ( μ i ψ ) c i ]   ,
where v is the velocity vector in m / s , D i are the diffusion coefficients in m 2 / s , and μ i   = D i e z i / k T are the electrical mobilities in m 2 / ( Vs ) . Here, e = 1.602   ×   10 19   C is the elementary charge, k = 8.617   ×   10 5   eV / K is the Boltzmann constant, and T is the temperature in K . RheoEFoam converges Equations (1) and (2) before solving the Navier–Stokes equations:
× v = 0   ,
ρ ( v t + v × v ) =   p +   η 2 v + f E   ,
to obtain updates for the velocities, v , and pressure, p (in N / m 2 ). The mixture density, ρ = c i M i in kg / m 3 , is computed from the molar masses, M i , in kg / mol . η is the dynamic viscosity in k g / ( m s ) . For the computation of the electric force per unit volume, f E = ρ E E in N / m 3 , the electric permittivity is assumed as constant. Finally, E = ψ is the electric field in V / m . Once Equations (3) and (4) are solved, rheoEFoam goes back to computing Equations (1) and (2) until the entire system of Equations (1)–(4) converges.

3.2. Boundary Conditions

The present simulations were set up according to the 2-D problem by Druzgalski et al. [12]. This two-species problem is defined by a reservoir at the top surface and a CEM at the bottom surface that allows cationic species to pass (Figure 2). An electric field ( E   =   Δ V / H ) is applied in the direction perpendicular to the reservoir and the membrane.
The CEM was placed at y = 0. The membrane was modeled by no-slip and no-penetration conditions for the velocities, v = 0 , a zero normal pressure gradient, a zero potential, ψ = 0 , a fixed cation concentration, c + , and a no-flux condition for the anions. No-slip and no-penetration conditions for the velocities and a zero normal gradient for the pressure were also applied at the reservoir boundary at y = H . The potential, ψ , at the reservoir boundary was set according to the applied voltage differential, Δ V . The cation and anion concentrations, c + and c , at the reservoir boundary were held constant. Periodic (i.e., cyclic) boundary conditions were employed at all other boundaries.

3.3. Defining Parameters

A case with very low NaCl concentration ( 1   mol / m 3 ) matching that of an earlier simulation by Pimenta and Alves [15] was chosen as the baseline. In addition, cases with higher salt concentrations matching slightly and moderately brackish water (low and medium), as well as seawater (high concentration), were considered. The concentrations for the different cases are listed in Table 1. The temperature was set to room temperature, T = 300.7   K , the diffusion coefficients were set to D = 10 9   m 2 / s . According to Vitagliano and Lyons [22], the diffusion coefficient for NaCl solutions has the same order of magnitude and varies slightly with concentration. The mixture density was ρ = 1000   kg / m 3 (i.e., water), the dynamic viscosity was taken as η = 0.001   kg / ( ms ) , and the kinematic viscosity was ν = 10 6   m 2 / s .

3.4. Computational Grid and Timestep

Structured computational grids with equidistant line spacing in the direction parallel to the membrane were created with OpenFoam blockMesh (an example is shown in Figure 3). Grid line clustering was employed near the membrane to properly resolve the boundary layer. The domain dimensions (width × height × depth) are denoted as L   ×   H   ×   Z . It was found that the grid resolution and the timestep had to be adjusted depending on the concentrations and electric field strength to resolve all relevant scales and maintain numerical stability. The domain dimensions as well as the number of cells and the computational timestep are listed in Table 2. The computational timestep is dictated by numerical stability and was found through trial and error. All simulations were carried out on a 24-core Dell Precision T5600 desktop computer. As in the work by Pimenta and Alves [15], a one-dimensional (1-D) solution was computed first. The 1-D solution was then extended in the second and third dimension and a 1% random disturbance was added to the concentrations to raise the initial disturbance amplitudes.
For the baseline case, 2-D simulations with a wavy wall were carried out to investigate the effect of the membrane surface geometry on the current density (Figure 4). The 2-D roughness was periodic along the width of the domain. The average domain height, H , width, and depth of the domain as well as the number of cells in all three directions were left unchanged. The wavy wall geometry was defined with splines in blockMesh. The amplitude of the roughness, A , was either 0.01   μ m (for the majority of the cases; 0.02   μ m peak to peak) or 0.5   μ m   ( 1   μ m peak to peak).

3.5. Non-Dimensionalization

In the following, tildes identify dimensionless quantities. The voltage:
V ˜ = Δ V V T   ,
was made dimensionless with the thermal voltage, V T = k T / ( e z ) . Time was made dimensionless with the diffusion coefficient and the domain height:
t ˜ = t D H 2 .
Finally, the dimensionless current density is:
J ˜ = J H D c 0 z F .

4. Results

4.1. Effect of Dimensionless Potential

As two representative cases for the baseline ( 1   mol / m 3 NaCl concentration) 2-D simulations, in Figure 5 and Figure 6, instantaneous visualizations of the anion concentration, vertical velocity, and electric potential are shown for dimensionless voltages of V ˜ = 30 and 40. Initially, at t ˜ = 0.0008 for V ˜ = 30 and t ˜ = 0.0004 for V ˜ = 40, the solutions are very close to the 1-D initial condition. However, the vertical velocity iso-contours reveal that disturbances are growing as a result of the EKI. The instability sets in earlier for the higher voltage. Once the disturbance amplitude is above a certain threshold, periodic structures can be discerned above the membrane. The structures are also seen in the anion concentration but less so in the electrical potential. As time progresses, the disturbances are growing and for t ˜ = 0.01 ( V ˜   = 30) and 0.0008 ( V ˜ = 40) they are strong enough to result in a significant distortion of the electrical potential. The merging of neighboring structures results in the formation of larger structures such as seen for the last two time instances for V ˜ = 40.
When the voltage is increased to V ˜ = 120, the instability is stronger, the disturbances grow faster, and the structures show up much earlier (compare anion concentrations in Figure 5, Figure 6 and Figure 7). For t ˜ = 0.001, smaller structures have coalesced into bigger structures in a process similar to vortex merging in fluid mechanics. For t ˜ = 0.001 and t ˜ = 0.1, the solution appears chaotic and much of the ion transfer across the height of the domain is accomplished by convection instead of diffusion.
In Figure 8, the current density, averaged over the surface area of the membrane, for six different voltages is plotted versus time. For voltages of V ˜ = 10 and 20, the current density remains steady around J ˜ 2 , which suggests that no instability occurs. As the voltage is increased to V ˜ = 30, the current density shoots up to J ˜ 5 around t ˜ = 0.004. At this time instant, the near-membrane structures are growing strongly as a result of the primary instability. The current density then drops slowly to an average value of J ˜ 3 and remains essentially constant apart from a weak unsteadiness (Figure 8). Figure 5 reveals nearly steady structures ( t ˜   = 0.05 and 0.1). The structures increase the current density. An observed slight “wobble” of the structures leads to the minute oscillation of the current density. For voltages of V ˜ = 40 and above, the current density also exhibits an initial overshoot and then settles on a lower value. The oscillations are, however, much stronger and more stochastic than for V ˜ = 30. Based on a comparison of Figure 5 and Figure 6, it can be concluded that for V ˜ 30 , the structures become unsteady (in a fluid dynamics context, this would be referred to as the onset of secondary instability). In the simulations, the unsteadiness grows from truncation and machine roundoff errors. In experiments, slight imperfections of the experimental setup or environmental factors can provide the initial disturbances. The unsteadiness appears random (randomness would typically be associated with turbulence) and the amplitude of the unsteadiness grows when the electric field strength is increased.
For the baseline case, a 3-D simulation was also performed for a voltage of V ˜ = 120. Instantaneous visualizations for t ˜ = 0.00056 are shown in Figure 9. As for the 2-D simulation, the EKI leads to the appearance of large-scale structures that effectively transport ions from the membrane to the reservoir. The wavelengths of the dominant structures for the 2-D and 3-D simulation are similar (Figure 7 and Figure 9) and the initial rise of the current density is almost identical (Figure 10). Because of this close similarity and the high cost of the 3-D simulation, it was decided to exclusively employ 2-D simulations for the following investigations.
The dimensionless current density, J ˜ , obtained from the 2-D simulations was averaged for 0.02 < t ˜ < 0.1 (Figure 11). Similar to experimental J-V curves (e.g., [9,23,24]), three distinct regimes can be observed for the smooth membrane. (i) For the Ohmic regime, the current density scales linearly with the applied voltage; (ii) in the limiting regime, the rate of increase of the current density decreases and the current density approaches an asymptotic value that is determined by the diffusion-limited transport of ions; and (iii) for the overlimiting regime, the current density increases again with the applied voltage. The existence of the overlimiting regime has been commonly attributed to electroconvection, which refers to the generation of vertical flow structures near the membrane. These structures, which are also observed in the present simulations for V ˜ 30 , facilitate a very effective convective mixing, which overcomes the diffusion-limited transport of ions. As a result, the current density increases dramatically. In the present simulations, the transition from the limiting regime to the overlimiting regime occurred at V ˜ 20 . A fine mesh simulation with a smooth membrane (black “x” in Figure 11) was carried out for V ˜ = 20 . Since the current density for the coarse and fine mesh were almost identical, the conclusion was drawn that the present results are grid-converged.

4.2. Effect of Salt Concentration

The NaCl concentration was raised to 10   mol / m 3 and 100   mol / m 3 (simulating slightly and moderately brackish water), and finally 400   mol / m 3 (seawater). For the medium and high concentration cases, the grid extent was decreased by a factor of 10 in all directions. Since the number of cells was kept constant, this resulted in a factor of 10 increase of the grid resolution and the electric field strength. As the concentration was increased, the structures appeared earlier and their size was reduced. The same mesh dimensions ( 60   μ m   ×   10   μ m ) were employed for the 1 and 10   mol / m 3 simulations (Figure 5, Figure 6 and Figure 7 and Figure 12a). For the latter, the initial structures were much smaller and appeared earlier. Similar observations can be made for the 100 and 400   mol / m 3 concentration cases, which were both computed on a 6   μ m   ×   1   μ m mesh (Figure 12b,c). This suggests that the onset of the instability is also (in addition to the voltage) dependent on the ion concentration and that the wavelength of the most unstable disturbances scales with the concentration.

4.3. Wavy Membrane

The main emphasis of the present investigation is on the effect of surface roughness on the current density. According to Rubinstein and Zaltzman [8], by making the membrane surface wavy, the initial wavelength of the disturbances can be controlled. All simulations with wavy membrane are for the baseline case with a 1   mol / m 3 concentration. It was decided to first try a membrane waviness with an amplitude of A = 0.5   μ m and a wavelength of λ = 10   μ m . Iso-contours of the anion concentration for a dimensionless voltage of V ˜ = 120 are shown in Figure 13. Compared to the smooth membrane results in Figure 7, the structures form much earlier and up to t ˜ = 0.00016 have a pronounced periodicity that is dictated by the wall waviness. Especially for t ˜ = 0.00012 and 0.00016, the structures are mushroom-shaped and similar in appearance to the structures resulting from buoyancy-driven instability. For t ˜ 0.001 , the initial periodicity is lost and the flow is more chaotic.
The wall waviness raises the initial amplitudes of the λ = 10   μ m disturbances, which dominate at the beginning. As the λ = 10   μ m mode saturates and stops growing, other modes are still growing and reach comparable amplitudes, which makes the anion concentration field appear chaotic. The earlier onset of the EKI suggests that the waviness increased the initial amplitude of an unstable mode. The decision was made to lower the amplitude of the membrane waviness to A = 0.01   μ m and to investigate a broader range of wavelengths.
The wavelength of the waviness was varied from λ = 1.5   μ m (40 periods per domain width; 12 grid intervals per period) to 20 μ m (3 periods per domain width; 120 grid intervals per period). Instantaneous visualizations of the anion concentration for five cases are provided in Figure 14 and Figure 15. Despite the very small amplitude of the membrane waviness, for λ = 3.75   μ m , 5   μ m , and 7.5   μ m , the wavelength of the primary structures resulting from the EKI is identical to the wavelength of the waviness (Figure 14). For λ = 3.75   μ m and λ = 5   μ m , the structures maintain the same wavelength up to t ˜ = 0.00016, after which the flow becomes chaotic. For λ = 7.5   μ m , the structures sustain their initial wavelength only up to t ˜ = 0.0008, after which they split up in two (halving of wavelength). The halving of the wavelength implies that the higher harmonic, λ = 3.75   μ m , grows more strongly.
For the larger wavelengths of λ = 10   μ m and 15   μ m , the visualizations for t ˜ = 0.00004 provide clear evidence that, although the disturbance that is directly related to the forcing is present, higher harmonics grow more strongly (Figure 15). For λ = 10   μ m , the wavelength is reduced to λ = 3.75   μ m (and interestingly not 5   μ m ) for t ˜ = 0.00016. This provides supporting evidence that the wavelength of the most amplified disturbance is close to λ = 3.75   μ m . For λ = 15   μ m , each structure splits up into two structures with a wavelength of approximately λ = 7.5   μ m (subharmonic of λ = 3.75   μ m ).
In Figure 16, the dimensionless current density for λ = 15   μ m is plotted versus time. Compared to the smooth membrane (Figure 8), the initial overshoots occur earlier, which implies that the structures form sooner. For V ˜ = 20, the current density is increased for t ˜ > 0.02 , suggesting that the wavy membrane leads to the formation of structures for a voltage where no structures form for the smooth membrane. This means that the overlimiting regime starts at a lower voltage as also seen in Figure 11. Irrespective of the applied voltage, the average current density for t ˜ > 0.02 is noticeably higher, which, in accordance with Rubinstein and Zaltzman [8], suggests that the waviness enhances the electroconvection. This is also seen in Figure 11 (compare the value for λ = 15   μ m with the value for the smooth membrane).
In Figure 17, the time- and area-averaged current density for the baseline 2-D case is plotted versus the wavelength, λ , of the membrane waviness. For very short wavelengths, λ 3   μ m , the current density matches that of the smooth membrane, J ˜ 2.1 , and the waviness (i.e., surface roughness) has no effect. For a minimal increase to λ = 3.75   μ m (wavelength of the most amplified structures), the current density is strongly increased ( J ˜ 2.7 ). Even stronger maxima of the current density are observed for integer multiples of this wavelength, λ = 2   ×   3.75   μ m = 7.5   μ m and λ = 4   ×   3.75   μ m = 15   μ m . By introducing wavelengths that are incommensurate with this natural disturbance wavelength (similar to Saric et al. [18]), the growth of the disturbance waves can be reduced. As a result, the current density is reduced. Figure 17 displays a strong reduction of the current density for λ = 5   μ m and 10   μ m .
Fine grid simulations were also carried out for two cases with the wavy membrane ( λ = 10   μ m and 15   μ m ). The coarse and fine grid solutions are in good agreement for the former and in reasonable agreement for the latter (Figure 11 and Figure 17).

5. Discussion

The overlimiting regime in electrodialysis (ED) is associated with the appearance of electroconvective structures near the membranes. The structures are caused by an electrokinetic instability (EKI). The present two-dimensional (2-D) simulations of a generic two-species problem with a smooth ion-selective membrane on the bottom and a water flow channel (reservoir) on the top indicate a dependence of the onset of the EKI on the applied voltage and concentration. With an increasing applied voltage and concentration, the instability becomes stronger and the structures form earlier. The initial perturbations were found to have a well-defined wavelength that is inversely related to the concentration. For very low ion concentrations ( 1   mol / m 3 ), the wavelength was of the order of 5   μ m . For ion concentrations of 100   mol / m 3 (moderately brackish water) and 400   mol / m 3 (seawater), the wavelength was in the order of 0.5   μ m .
For ED applications, it is desirable to delay the onset of the overlimiting regime. Similar to hydrodynamic stability problems in fluid mechanics, it was speculated that a favored natural most-amplified wavelength exists. Related research by Saric et al. [18] of the crossflow instability on swept wings has shown that by forcing a wavelength that is incommensurate with the wavelength of the most amplified natural disturbances, the appearance of nonlinear flow structures can be delayed. Whether the same idea can be transferred to ED applications was investigated. In the present simulations, disturbances were introduced through a 2-D periodic waviness of the membrane surface. It was found that a very minute waviness amplitude had a large effect on the initial growth of the disturbances and the time- and area-averaged current density. When the wavelength was 3.75   μ m , 2   ×   3.75 = 7.5   μ m , or 4   ×   3.75 = 15   μ m , in agreement with Rubinstein and Zaltzman [8], the current density increased and the overlimiting regime set in at a lower voltage. It was speculated that for these cases, the initial amplitude of the naturally occurring unstable mode was increased. Forcing with a wavelength of 5   μ m and 2   ×   5 = 10   μ m , on the other hand, delayed the formation of the structures and the onset of the overlimiting regime. When the wavelength of the membrane waviness was below approximately 3   μ m , it had no effect on the current density and the smooth membrane result was recovered.
Membranes are never perfectly smooth. The proposed simulations show a clear dependence of the minimum voltage for the onset of the EKI and thus the overlimiting regime on the geometric surface properties of the membrane. The membrane fine structure (order of 1   μ m for low concentrations and 0.1   μ m for higher concentrations) has implications for the onset of the overlimiting regime and the current density. It was found that by adding fine-scale structures to the membrane surface, the onset of the instability can be controlled. Compared to Güler et al. [21], the wavelength is several orders of magnitude smaller and the control is based on a different physical mechanism. Whether such membrane surface modifications are feasible remains to be seen. Finally, the 2-D assumption was made for the present parameter study. The extension to 3-D will be explored in the future.

Author Contributions

Conceptualization and methodology, A.G., X.X., and P.X.; investigation, software, formal analysis, and visualization, A.M. and P.C.G.; writing—original draft, A.G. and P.X.; writing—review and editing, A.G., X.X., and P.X.; funding acquisition, project administration, and supervision, A.G. and P.X.

Funding

This work was supported by a New Mexico State University College of Engineering mini-grant, and the National Science Foundation Engineering Research Center ReNUWIt (EEC-1028968).

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Xu, P.; Cath, T.Y.; Robertson, A.P.; Reinhard, M.; Leckie, J.O.; Drewes, J.E. Critical Review of Desalination Concentrate Management, Treatment and Beneficial Use. Environ. Eng. Sci. 2013, 30, 502–514. [Google Scholar] [CrossRef] [Green Version]
  2. Xu, P.; Capito, M.; Cath, T.Y. Selective removal of arsenic and monovalent ions from brackish water reverse osmosis concentrate. J. Hazard. Mater. 2013, 260, 885–891. [Google Scholar] [CrossRef] [PubMed]
  3. Chehayeb, K.M.; Lienhard, J.H. Entropy Generation Analysis of Electrodialysis. Desalination 2017, 413, 184–198. [Google Scholar] [CrossRef]
  4. Campione, A.; Gurreri, L.; Ciofalo, M.; Micale, G.; Tamburini, A.; Cipollina, A. Electrodialysis for water desalination: A critical assessment of recent developments on process fundamentals, models and applications. Desalination 2018, 434, 121–160. [Google Scholar] [CrossRef]
  5. Rubinstein, I.; Staude, E.; Kedem, O. Role of the membrane surface in concentration polarization at ion-exchange membrane. Desalination 1988, 69, 101–114. [Google Scholar] [CrossRef]
  6. Demekhin, E.A.; Shapar, E.M.; Lapchenko, V.V. Initiation of electroconvection in semipermeable electric membranes. Dokl. Phys. 2008, 53, 450–453. [Google Scholar] [CrossRef]
  7. Rubinstein, I.; Zaltzman, B. Electro-convective versus electroosmotic instability in concentration polarization. Adv. Colloid Interface Sci. 2007, 134, 190–200. [Google Scholar] [CrossRef] [PubMed]
  8. Rubinstein, I.; Zaltzman, B. Electro-osmotically induced convection at a permselective membrane. Phys. Rev. E 2000, 62, 2238. [Google Scholar] [CrossRef]
  9. Chang, H.-C.; Yossifon, G.; Demekhin, E.A. Nanoscale Electrokinetics and Microvortices: How Microhydrodynamics Affects Nanofluidic Ion Flux. Ann. Rev. Fluid Mech. 2012, 44, 401–426. [Google Scholar] [CrossRef] [Green Version]
  10. Yossifon, G.; Mushenheim, P.; Chang, H.-C. Controlling the nanoslot overlimiting current by varying the micro-chamber depth. Europhys. Lett. 2010, 90, 64004. [Google Scholar] [CrossRef]
  11. Kim, Y.; Walker, W.S.; Lawler, D.F. Electrodialysis with spacers: Effects of variation and correlation of boundary layer thickness. Desalination 2011, 274, 54–63. [Google Scholar] [CrossRef]
  12. Druzgalski, C.L.; Andersen, M.B.; Mani, A. Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface. Phys. Fluids 2013, 25, 110804. [Google Scholar] [CrossRef]
  13. Zourmand, Z.; Faridirad, F.; Kasiri, N.; Mohammadi, T. Mass transfer modeling of desalination through an electrodialysis cell. Desalination 2015, 359, 41–51. [Google Scholar] [CrossRef]
  14. Enciso, R.; Delgadillo, J.A.; Dominguez, O.; Rodríguez-Torres, I. Analysis and validation of the hydrodynamics of an electrodialysis cell using computational fluid dynamics. Desalination 2017, 408, 127–132. [Google Scholar] [CrossRef]
  15. Pimenta, F.; Alves, M.A. Numerical simulation of electrically-driven flows using OpenFOAM. arXiv 2018, arXiv:Physics/1802.02843v2. [Google Scholar]
  16. Demekhin, E.A.; Nikitin, N.V.; Shelistov, V.S. Three-dimensional coherent structures of electrokinetic instability. Phys. Rev. E 2014, 90, 013031. [Google Scholar] [Green Version]
  17. Druzgalski, C.L.; Mani, A. Statistical analysis of electroconvection near an ion-selective membrane in the highly chaotic regime. Phys. Rev. Fluids 2016, 1, 073601. [Google Scholar] [CrossRef]
  18. Saric, W.S.; Carpenter, A.L.; Reed, H.L. Passive control of transition in three-dimensional boundary layers, with emphasis on discrete roughness elements. Philos. Trans. R. Soc. A 2011, 369, 1352–1364. [Google Scholar] [CrossRef]
  19. Embacher, M.; Fasel, H.F. Direct numerical simulations of laminar separation bubbles: Investigation of absolute instability and active flow control of transition to turbulence. J. Fluid Mech. 2014, 747, 141–185. [Google Scholar] [CrossRef]
  20. Walsh, M.J. Riblets as a viscous drag reduction technique. AIAA J. 1983, 21, 485–486. [Google Scholar] [CrossRef]
  21. Güler, E.; Elizen, R.; Saakes, M.; Nijmeijer, K. Micro-structured membranes for electricity generation by reverse electrodialysis. J. Membr. Sci. 2014, 458, 136–148. [Google Scholar] [CrossRef]
  22. Vitagliano, V.; Lyons, P.A. Diffusion Coefficients for Aqueous Solutions of Sodium Chloride and Barium Chloride. J. Am. Chem. Soc. 1956, 78, 1549–1552. [Google Scholar] [CrossRef]
  23. Rubinstein, I.; Shtilman, L. Voltage against current curves of cation exchange membrane. J. Chem. Soc. Faraday Trans. 2 Mol. Chem. Phys. 1979, 75, 231–246. [Google Scholar] [CrossRef]
  24. Krol, J.J.; Wessling, M.; Strathmann, H. Concentration polarization with monopolar ion exchange membranes: Current-voltage curves and water dissociation. J. Membr. Sci. 1999, 162, 145–154. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of electrodialysis stack.
Figure 1. Schematic diagram of electrodialysis stack.
Water 11 01397 g001
Figure 2. Generic problem of a 2-D reservoir with a cation-selective membrane (CEM).
Figure 2. Generic problem of a 2-D reservoir with a cation-selective membrane (CEM).
Water 11 01397 g002
Figure 3. Zoom-in of computational grid near membrane.
Figure 3. Zoom-in of computational grid near membrane.
Water 11 01397 g003
Figure 4. Outline of computational domain with wavy wall.
Figure 4. Outline of computational domain with wavy wall.
Water 11 01397 g004
Figure 5. Baseline 2-D case, V ˜ = 30. Smooth membrane. From left to right: Anion concentration, vertical velocity in m / s , electric potential in V .
Figure 5. Baseline 2-D case, V ˜ = 30. Smooth membrane. From left to right: Anion concentration, vertical velocity in m / s , electric potential in V .
Water 11 01397 g005
Figure 6. Baseline 2-D case, V ˜   = 40. Smooth membrane. From left to right: Anion concentration, vertical velocity in m / s , electric potential in V .
Figure 6. Baseline 2-D case, V ˜   = 40. Smooth membrane. From left to right: Anion concentration, vertical velocity in m / s , electric potential in V .
Water 11 01397 g006
Figure 7. Baseline 2-D case ( 1   mol / m 3 concentration), V ˜ = 120. Smooth membrane. Anion concentration.
Figure 7. Baseline 2-D case ( 1   mol / m 3 concentration), V ˜ = 120. Smooth membrane. Anion concentration.
Water 11 01397 g007
Figure 8. Baseline 2-D case. Area-averaged current density for the smooth membrane.
Figure 8. Baseline 2-D case. Area-averaged current density for the smooth membrane.
Water 11 01397 g008
Figure 9. Baseline 3-D case, V ˜ = 120. Smooth membrane. Iso-contours of anion concentration and iso-surfaces of anion concentration of 0.2 for t ˜ = 0.00056 .
Figure 9. Baseline 3-D case, V ˜ = 120. Smooth membrane. Iso-contours of anion concentration and iso-surfaces of anion concentration of 0.2 for t ˜ = 0.00056 .
Water 11 01397 g009
Figure 10. Baseline case, V ˜ = 120. Area-averaged current density for the smooth membrane.
Figure 10. Baseline case, V ˜ = 120. Area-averaged current density for the smooth membrane.
Water 11 01397 g010
Figure 11. Baseline 2-D case. Time- and area-averaged current density.
Figure 11. Baseline 2-D case. Time- and area-averaged current density.
Water 11 01397 g011
Figure 12. Two-dimensional simulations for V ˜   = 120 and (a) 10   mol / m 3 , (b) 100   mol / m 3 , and (c) 400   mol / m 3 . Smooth membrane. Anion concentration.
Figure 12. Two-dimensional simulations for V ˜   = 120 and (a) 10   mol / m 3 , (b) 100   mol / m 3 , and (c) 400   mol / m 3 . Smooth membrane. Anion concentration.
Water 11 01397 g012
Figure 13. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.5   μ m and λ = 10   μ m . Anion concentration.
Figure 13. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.5   μ m and λ = 10   μ m . Anion concentration.
Water 11 01397 g013
Figure 14. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.01   μ m and (a) λ = 3.75   μ m , (b) λ = 5   μ m and (c) λ = 7.5   μ m . Anion concentration.
Figure 14. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.01   μ m and (a) λ = 3.75   μ m , (b) λ = 5   μ m and (c) λ = 7.5   μ m . Anion concentration.
Water 11 01397 g014
Figure 15. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.01   μ m and (a) λ = 10   μ m and (b) λ = 15   μ m . Anion concentration.
Figure 15. Baseline 2-D case, V ˜ = 120. Wavy membrane with A = 0.01   μ m and (a) λ = 10   μ m and (b) λ = 15   μ m . Anion concentration.
Water 11 01397 g015
Figure 16. Baseline 2-D case. Area-averaged current density for the wavy membrane with A = 0.01   μ m and λ = 15   μ m .
Figure 16. Baseline 2-D case. Area-averaged current density for the wavy membrane with A = 0.01   μ m and λ = 15   μ m .
Water 11 01397 g016
Figure 17. Baseline 2-D case, V ˜ = 20 (transition from limiting to overlimiting regime). Time- and area-averaged current density.
Figure 17. Baseline 2-D case, V ˜ = 20 (transition from limiting to overlimiting regime). Time- and area-averaged current density.
Water 11 01397 g017
Table 1. Simulated cases.
Table 1. Simulated cases.
Case Reservoir   Species   Concentrations   in   m o l / m 3 Reservoir   Species   Concentrations   in   m g / L Membrane   Cation   Concentration   in   m o l / m 3
Baseline178.52
Low10157020
Medium1007850200
High40031,400800
Table 2. Details of computational grids and timestep.
Table 2. Details of computational grids and timestep.
Case Domain   Dimensions ,   L   ×   H   ×   Z   in   μ m Number of Cells Timestep   for   1 D   Simulations   in   s Timestep   for   2 D / 3 D   Simulations   in   s
Baseline, 2-D60 × 10 × 1480 × 90 × 1 2   ×   10 6 1   ×   10 7
Baseline, 2-D, fine mesh60 × 10 × 1960 × 180 × 1 1   ×   10 6 0.5   ×   10 7
Baseline, 3-D30 × 10 × 30240 × 90 × 240 2   ×   10 6 1   ×   10 7
Low, 2-D60 × 10 × 1480 × 90 × 1 1   ×   10 7 2   ×   10 9
Medium, 2-D6 × 1 × 0.1480 × 90 × 1 2   ×   10 8 4   ×   10 10
High, 2-D6 × 1 × 0.1480 × 90 × 1 1   ×   10 10 1   ×   10 10

Share and Cite

MDPI and ACS Style

Gross, A.; Morvezen, A.; Castillo Gomez, P.; Xu, X.; Xu, P. Numerical Investigation of the Effect of Two-Dimensional Surface Waviness on the Current Density of Ion-Selective Membranes for Electrodialysis. Water 2019, 11, 1397. https://doi.org/10.3390/w11071397

AMA Style

Gross A, Morvezen A, Castillo Gomez P, Xu X, Xu P. Numerical Investigation of the Effect of Two-Dimensional Surface Waviness on the Current Density of Ion-Selective Membranes for Electrodialysis. Water. 2019; 11(7):1397. https://doi.org/10.3390/w11071397

Chicago/Turabian Style

Gross, Andreas, Arthur Morvezen, Pedro Castillo Gomez, Xuesong Xu, and Pei Xu. 2019. "Numerical Investigation of the Effect of Two-Dimensional Surface Waviness on the Current Density of Ion-Selective Membranes for Electrodialysis" Water 11, no. 7: 1397. https://doi.org/10.3390/w11071397

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop