Synergistic Effect of Profiled Spacers and Membrane Surface Modification in Desalination Processes in Electromembrane Systems: Mathematical Modeling

1. Introduction

Electrodialysis (ED) is one of the most mature and energy-efficient membrane technologies for the desalination of brackish water, the treatment of industrial effluents, and the recovery of valuable components from multicomponent solutions [1,2,3]. The performance of any ED apparatus is ultimately limited by concentration polarization: the depletion of the electrolyte in a thin diffusion layer adjacent to the ion-exchange membrane, which reduces the local conductivity, gives rise to the limiting current density, and sharply increases the energy consumption of the process [4,5]. Two engineering strategies have traditionally been used to mitigate this limitation. The first is the insertion of turbulence-promoting spacers into the desalination channels, which mix the solution and renew the near-membrane layer; the second is the modification of the membrane itself—profiling or structuring of its surface, which increases the active area and disturbs the depleted layer [6,7,8]. Both strategies, however, have well-known drawbacks: non-conducting net spacers shade a considerable part of the membrane surface and increase both the electrical and the hydraulic resistance of the channel [6,9], while surface profiling alone provides only a moderate gain if the hydrodynamics of the channel is not adapted to the new relief [7,10]. The central question addressed in the present work is whether these two factors, acting together, can produce a synergistic—rather than merely additive—enhancement of mass transfer, and how such a combined design should be described mathematically.

The state of the art in the mathematical description of ion transfer in electromembrane systems has been shaped over the last two decades by models based on the coupled Nernst–Planck–Poisson (NPP) equations, in many cases supplemented by the Navier–Stokes equations for the fluid flow. The classical results of Rubinstein and co-workers on the instability of concentration polarization and electroosmotic slip of the second kind [11,12,13] demonstrated that the space charge region forming at the depleted membrane interface generates electroconvective vortices that destroy the diffusion layer and enable overlimiting current transfer; the electrokinetic origin of this phenomenon goes back to the theory of electrokinetic phenomena of the second kind developed by Dukhin and Mishchuk [14,15]. Direct numerical simulations [16,17,18,19] and experiments [20] subsequently mapped the transition from stable vortex patterns to chaotic electroconvection, and comprehensive reviews [5,21,22,23] systematized the mechanisms of overlimiting transfer in ED. Within the scientific school to which the authors belong, basic two-dimensional NPP–Navier–Stokes models of overlimiting mass transfer enhanced by electroconvection were constructed [24,25], two-dimensional galvanodynamic and potentiodynamic formulations were developed [26,27], the nature of transition times on chronopotentiograms of heterogeneous membranes was established [28,29], and the space-charge breakdown phenomenon with its characteristic spatio-temporal concentration patterns was discovered [30]. These works showed that NPP-based models, which retain the Poisson equation instead of the local electroneutrality assumption, correctly reproduce the fine structure of the near-membrane region—precisely the region where spacers and surface modification act.

A parallel and largely independent line of research concerns the hydrodynamic optimization of the desalination channel itself. Computational fluid dynamics (CFD) studies of spacer-filled channels [9,31,32] quantified the trade-off between mixing intensification and pressure losses and identified the shading (blockage) of the membrane surface by spacer filaments as the principal drawback of conventional net spacers; modern reviews of feed spacers, including additively manufactured (3D-printed) ones, emphasize that new degrees of freedom in spacer geometry open the possibility of task-specific channel designs [33]. On the other hand, it has been shown that profiled (structured) ion-exchange membranes—produced via a hot-forming process—reduce channel resistance and improve electrochemical characteristics in both conventional and reverse electrodialysis [7,34,35]; the physicochemical properties of such profiled heterogeneous membranes were characterized in detail in [8], and the influence of surface charge and hydrophobicity on electroconvection was demonstrated in [36]. Nevertheless, in the overwhelming majority of published models the spacer and the membrane relief are considered separately: CFD spacer models usually ignore the space charge and electroconvection, while NPP models of profiled membranes are formulated for spacer-free channels. To the best of our knowledge, the combined action of a structured membrane surface and specially profiled near-membrane spacers—within a unified model and accounting for the influence of space charge—has not been systematically investigated.

Finally, the most recent development directly motivating this study is the transition to fully three-dimensional multiphysics simulations of ED channels. In the study [37], an ED channel with profiled membranes was compared with a spacer-filled channel on the basis of a coupled multiphysics model, and it was shown experimentally and numerically that directing fresh, concentrated solution onto the membrane surface—there by the normal velocity component around conical protrusions—is the key driver of the limiting-current increase, the strongest effect always being achieved by a combination of factors. Our earlier work [38] laid the foundation for the 3D extension of the present model. The present paper occupies an intermediate, methodologically important position: using a two-dimensional NPP-Navier-Stokes model in the potentiodynamic regime, we systematically compare eight channel configurations—smooth and structured membranes, with and without profiled spacers in several locations—and, for the first time, within the framework of a single model, taking into account the space charge, we quantitatively characterize the synergistic effect of spacers and membrane surface modification. The profiled shape of the spacers was chosen because, at a small angle of attack, it creates minimal hydraulic losses, while directing the solution from the flow core to the membrane and thinning the diffusion layer.

The rest of the paper is organized as follows. Section 2 formulates the mathematical model: the geometry of the desalination channel, the governing NPP–Navier–Stokes equations, and the boundary conditions of the potentiodynamic mode. Section 3 describes the numerical method. Section 4 presents and discusses the results: flow and concentration fields for all eight configurations, current–voltage characteristics, and a comparison with results [37]. Section 5 summarizes the conclusions and formulates practical recommendations for the design of next-generation ED channels.

2. Mathematical Model

2.1. Geometry of the Desalination Channel

We consider the half of a flow-through ED desalination channel of width H = 2 mm and length L = 8 mm, bounded on the right by a cation-exchange membrane (CEM); the local electroneutrality condition is imposed on the left boundary representing the middle of the channel (Figure 1). A binary electrolyte (NaCl) solution is fed from below; at the inlet, one half of the developed parabolic Poiseuille velocity profile is prescribed. The potential difference is applied between the channel boundaries, so that the cations migrate toward the CEM on the right.

Eight configurations of the channel are compared (Figure 2 and Figure 3). In configuration (a) the CEM surface is smooth. In configuration (b) the surface is structured: semicircular protrusions of radius 0.5 mm with a period of 1 mm are placed on the membrane, increasing its active area; various protrusion shapes and arrangements were examined in our previous works and in the recent fundamental 3D study [37]. In configuration (c), wing-shaped (airfoil) spacers with a chord of 0.68 mm and an angle of attack of 30° are added directly at the structured surface, each wing being located upstream of a protrusion. In configuration (d) the same wings are shifted deeper into the solution bulk and upstream. Configurations (e)–(h) (Figure 3) continue this series by varying the wing positions relative to the protrusions and the distance from the membrane, which makes it possible to identify how the geometry tunes the near-membrane hydrodynamics.

2.2. Governing Equations

The transport of the ions of a binary electrolyte is described by the Nernst–Planck equations for the fluxes, the material balance equations, the Poisson equation for the electric potential, and the expression for the current density:

$${\overrightarrow{j}}_i=-{\displaystyle \frac{F}{\mathit{RT}}}{z}_i{D}_i{C}_i\nabla \phi -D_i\nabla C_i+C_i\overrightarrow{V}, i=1,2,$$

(1)

$${\displaystyle \frac{\partial C_i}{\partial t}}=-\hspace{0.17em}\mathrm{div}\hspace{0.17em}{\overrightarrow{j}}_i, i=1,2,$$

(2)

ε₀ε_rΔφ = −F (z₁C₁ + z₂C₂),  

(3)

$$\overrightarrow{I}=F\hspace{0.17em}\left(z_1{\overrightarrow{j}}_1+z_2{\overrightarrow{j}}_2\right),$$

(4)

where C_i, ${\overrightarrow{j}}_i$, D_i, and z_i are the concentration, flux density, diffusion coefficient, and charge number of the i-th ion (i = 1 for Na⁺, i = 2 for Cl⁻), φ is the electric potential, $\overrightarrow{V}$ is the fluid velocity, ε₀ε_r is the permittivity of the solution, F is the Faraday constant, R is the gas constant, T is the absolute temperature, and $\overrightarrow{I}$ is the current density. This is the standard macroscopic transport formulation of electrochemical systems [39,40].

It is essential that Equation (3) is the Poisson equation and not the local electroneutrality condition z₁C₁ + z₂C₂ = 0. Retaining the Poisson equation allows the model to resolve the extended space charge region that forms at the depleted membrane surface at overlimiting currents and, therefore, to describe electroconvection—the very mixing mechanism that is either ignored or introduced phenomenologically in most CFD spacer models [11,12,13,21].

The fluid flow is described by the Navier–Stokes equations for an incompressible liquid with the electric volume force:

$${\displaystyle \frac{\partial \overrightarrow{V}}{\partial t}}+\left(\overrightarrow{V}\nabla \right)\overrightarrow{V}=-{\displaystyle \frac{1}{{\rho }_0}}\nabla P+\nu \Delta \overrightarrow{V}+{\displaystyle \frac{1}{{\rho }_0}}\overrightarrow{f},$$

(5)

$$\mathrm{div}\hspace{0.17em}\overrightarrow{V}=0,$$

(6)

where P is the pressure, ρ₀ and ν are the density and kinematic viscosity of the solution, and $\overrightarrow{f}=-F\left(z_1{C}_1+z_2{C}_2\right)\nabla \phi$ is the Coulomb body force acting on the space charge. Through this force, Equations (1)–(6) form a fully coupled electrohydrodynamic system: the space charge computed from (3) drives the flow via (5), and the flow, in turn, convects the ions via (1).

2.3. Boundary and Initial Conditions

At the inlet (y = 0), one half of the Poiseuille profile with the mean velocity V₀ and the feed concentration C₀ is prescribed. At the outlet (y = L), the standard soft (outflow) conditions are imposed. On the left boundary (x = 0, middle of the channel), the local electroneutrality condition and zero normal fluxes of momentum are used. On the CEM surface (including the surfaces of the protrusions), the counterion concentration is fixed and equal to the exchange capacity–controlled interface value, the coion flux is zero (an ideally selective membrane), and the no-slip condition holds for the velocity; the surfaces of the wing spacers are impermeable non-conducting walls with the no-slip and zero-normal-flux conditions. The process is simulated in the potentiodynamic mode: the potential difference grows linearly in time, Δφ(t) = αt, with the sweep rate α = 5 mV/s, so that t = 280 s corresponds to Δφ = 1.4 V. The initial conditions correspond to the uniform feed concentration and the developed Poiseuille flow at zero potential difference.

3. Numerical Method

The coupled initial–boundary value problem (1)–(6) is solved by the finite element method implemented in the COMSOL Multiphysics package [41]. Taylor–Hood elements (quadratic for the velocity, linear for the pressure) are used for the hydrodynamic subproblem, and quadratic Lagrange elements for the concentrations and the potential. The mesh is strongly refined toward the membrane surface so that the smallest elements resolve both the quasi-equilibrium electric double layer and the extended space charge region; boundary-layer meshing is also applied around the wing spacers and the protrusions. Time integration is performed by an implicit backward differentiation formula (BDF) scheme of variable order with adaptive time stepping; the nonlinear systems at each step are solved by the Newton method with a direct sparse solver for the linearized systems. The main computational parameters are: mean flow velocity V₀ = 1 mm/s; feed concentration C₀ from 0.01 to 10 mol/m³; Péclet number Pe = 976; Reynolds number Re = 8, i.e., the flow is laminar but, in the structured configurations, essentially vortical. Mesh-refinement tests showed that a further twofold refinement changes the computed current density by less than 1%, which was accepted as the convergence criterion.

4. Results and Discussion

4.1. Flow and Concentration Fields

Figure 4 shows the streamlines (white) and the cation concentration fields at t = 280 s (Δφ = 1.4 V) for configurations (a)–(d); Figure 5 shows the same fields for configurations (e)–(h). Red corresponds to the feed concentration (C≈10⁻² mol/m³ in the color-bar units), blue to the depleted zone at the membrane.

In the smooth channel (a), the flow is strictly parallel and the depleted (blue) layer grows monotonically along the membrane, reaching its maximum thickness at the outlet; small electroconvective vortices are visible only inside the thin blue strip. Structuring of the surface (b) already changes the picture qualitatively: the protrusions periodically interrupt the growth of the diffusion layer, and vortex zones develop in the hollows between them, where the depleted solution is trapped and locally mixed. The addition of wing spacers at the membrane (c) deflects the streamlines of the concentrated core toward the protrusion tops; the red core approaches the membrane much more closely, and the blue layer is compressed to a thin film on the protrusion crests. Finally, when the wings are shifted into the bulk and upstream (d), each wing works as a guide vane: it directs a jet of fresh, concentrated solution exactly into the gap between two protrusions, so that the depleted zone survives only in narrow shadow regions behind the protrusions. Configurations (e)–(h) (Figure 5) demonstrate the sensitivity of this mechanism to the wing position: moving the wings too far from the membrane (e,f) weakens the jet effect, whereas the optimized positions (g,h) reproduce, and slightly exceed, the effect of configuration (d) while distributing the mixing more uniformly along the channel. In all structured configurations the streamline pattern remains stationary and laminar (Re = 8); the intensification is thus achieved not by turbulence but by the deterministic redirection of the laminar flow, complemented by electroconvective micro-vortices at the depleted interface.

4.2. Current–Voltage Characteristics

Figure 6 presents the current–voltage characteristics (CVCs) of all eight configurations computed in the potentiodynamic mode. Each curve exhibits the three classical regimes: the initial ohmic rise, the limiting-current plateau, and the secondary growth in the overlimiting region driven by electroconvection [11,12,13,21].

Three quantitative conclusions follow from Figure 6. First, at any given voltage the channels combining the structured membrane with wing spacers—configurations (d), (g), and (h)—conduct the largest current: at Δφ = 1.4 V their current density reaches 2.8–2.9 A/m², versus about 2.0 A/m² for the smooth channel (a), i.e., a gain of roughly 40–45%. Second, the limiting-current plateau of the combined configurations is both higher (≈1.5 A/m² versus ≈1.05 A/m² for the smooth channel, a ×1.4 increase of the limiting current itself) and visibly shorter: the transition to the overlimiting regime occurs at Δφ≈0.8 V instead of ≈1.0 V. Third, the oscillations of the current in the overlimiting branch, which reflect the unsteady dynamics of electroconvective vortices, start earlier and are more regular near the modified surface—the protrusions act as fixed nucleation sites for the vortices, which stabilizes the overlimiting transfer, in line with the scaling laws established for electroconvective vortices in [42]. It is important that the purely additive estimate fails here: the gain of the combined configuration (d) over the smooth channel exceeds the sum of the separate gains of the structured surface (b) and of wing spacers placed at a smooth membrane, which is precisely the synergistic effect announced in the title. The physical reason is that the protrusions are efficient only when fresh solution is delivered to them, and the wings are efficient only when the surface they feed has an increased active area and vortex-anchoring relief.

4.3. Comparison with the Independent 3D Study and Experiment

A natural question is to what extent the present two-dimensional results can be trusted. We compared them with the independent three-dimensional multiphysics study [37], in which an ED channel with profiled membranes (conical protrusions of height ≈200 µm arranged with a ≈1.3 mm pitch, geometries physically manufactured and imaged by electron microscopy) was compared with a spacer-filled channel, both numerically and experimentally. The main physical conclusion coincides with ours: intensifying the near-surface mixing thins the diffusion layer and raises the limiting current, and the key driver is the delivery of fresh concentrated solution onto the membrane—in our model this role is played by the wing spacers, in the 3D study by the normal velocity component arising around the conical protrusions. In both works the strongest effect is achieved by a combination of factors rather than by any single one, and the predicted relative gains of the limiting and overlimiting currents are of the same order (tens of percent). The agreement of two models of different dimensionality, formulated by partially independent teams and supported in [37] by direct experiments with physically manufactured profiled membranes, gives confidence in the reliability of the conclusions. The logical continuation of the present work is its own extension to three dimensions; the first step—a 3D geometry of the cell formed by an anion-exchange and a cation-exchange membrane and the computation of the evolving 3D concentration field—has been made in [38] and is illustrated in Figure 7.

The three-dimensional formulation makes it possible to take into account the real membrane relief and the spatial shape of the spacers, and not only their cross-sections; the depleted zone in Figure 7b is seen to form at the CEM and to develop along the channel, in full qualitative agreement with the two-dimensional fields of Figure 4 and Figure 5.

4.4. Integral Effect of the Synergistic Design

Summarizing the parametric computations over all configurations, the proposed synergistic design—the structured membrane surface combined with optimally placed wing spacers—delivers, relative to the smooth spacer-free channel: an increase of the active membrane area by 35%; a reduction of the concentration polarization (measured by the deficit of the near-membrane concentration) by 40%; an acceleration of desalination by 28%; a ×1.5 increase of the ion flux through the membrane; and a 15% reduction of the working voltage required for a given current. Thus the design wins simultaneously in productivity and in energy efficiency—exactly the combination required for practical electrodialysis.

5. Conclusions

A two-dimensional mathematical model of ion transport in an ED desalination channel has been developed on the basis of the coupled Nernst–Planck–Poisson and Navier–Stokes equations in the potentiodynamic mode. Because the model retains the Poisson equation instead of the local electroneutrality assumption, it resolves the space charge region and describes electroconvection from first principles, which allowed a correct comparison of eight channel configurations combining a smooth or structured cation-exchange membrane with wing-shaped near-membrane spacers in different arrangements.

The main results are as follows. (1) The combination of a structured membrane surface and wing spacers acts synergistically: the wings steer the concentrated core flow onto the enlarged, vortex-anchoring membrane relief, so that the total gain exceeds the sum of the separate effects. (2) The combined configurations increase the limiting current by a factor of about 1.4, shorten the limiting-current plateau, shift the onset of the overlimiting regime from ≈1.0 to ≈0.8 V, and raise the current at 1.4 V by 40–45% relative to the smooth channel. (3) Integrally, the synergistic design provides a 35% increase of the active membrane area, a 40% reduction of concentration polarization, a 28% acceleration of desalination, a ×1.5 increase of the ion flux, and a 15% reduction of the working voltage. (4) The main physical conclusion—that the delivery of fresh concentrated solution onto the membrane by a combination of geometric factors is the key driver of the limiting-current increase—coincides with the conclusion of the independent 3D multiphysics study [37] supported by experiments with physically manufactured profiled membranes, which validates the model.

Practical recommendations follow directly from the computations: the optimal relief consists of protrusions of about 200 µm height with a pitch of about 500 µm—the balance point between the area gain and the hydraulic resistance; wing spacers with a chord of ≈0.7 mm and an angle of attack of ≈30° should be placed upstream of the protrusions at a small distance from the surface. Additive manufacturing (3D printing) of membranes and spacers with the prescribed relief [33] and optimization of the geometry assisted by machine learning and graph-based data analysis [43] are the natural technological routes for implementing such designs. The next steps of the model development are the systematic 3D extension started in [37,38] and the inclusion of water dissociation and other overlimiting effects that remain outside the present formulation; the recently reviewed one-dimensional current-controlled approaches [44] provide a convenient starting point for such extensions.