A Multi-Physical Modelling Tool For Reverse Electrodialysis: Alessandro - Tamburini@unipa - It
A Multi-Physical Modelling Tool For Reverse Electrodialysis: Alessandro - Tamburini@unipa - It
A Multi-Physical Modelling Tool For Reverse Electrodialysis: Alessandro - Tamburini@unipa - It
Abstract
the energy content of two salt solutions into electrical current by means of a selective controlled
mixing. RED process physics involves the interaction of several phenomena of different nature and
space-time scales. Therefore, mathematical modelling and numerical simulation tools are crucial for
In this work, a multi-physical modelling approach for the simulation of RED units was
developed. A periodic portion of a single cell pair was simulated in 2-D. Fluid dynamics was
simulated by the Navier-Stokes and continuity equations, and ionic mass transfer was simulated by
the Nernst–Planck approach along with the local electroneutrality condition. Moreover, the Donnan
A sensitivity analysis of the process performance was carried out. Different membrane/channel
geometrical configurations were investigated, including flat membranes, either with or without non-
conductive spacers, and profiled membranes. The influence of feeds concentration/velocity was also
evaluated. Results showed that, with respect to the ideal case of plane (empty) channels and planar
membranes, non-conductive spacers always reduce the power produced, while profiled membranes
may or may not perform better, depending on stack features and operating conditions.
1
Keywords: Reverse electrodialysis, multi-physical model, finite element method, power
1 Introduction
In the last decades, the need of an energy diversification and the environmental issues have
prompted research towards renewable resources. Salinity gradient power, i.e. the power extractable
from two salt solutions with different concentrations, can play an important role. Different
technologies have been proposed in order to exploit salinity gradients [1–4]. Reverse electrodialysis
(RED) was the first concept proposed [5] and appears to be the most promising, as it exhibits the
highest power density [6] and can become competitive with the other renewable energy processes
potential of a salinity gradient (e.g. from NaCl solutions) into electrical energy supplied to an external
load (Figure 1). The basic unit of a RED stack is the cell pair, comprising: an anionic exchange
membrane (CEM), and a dilute compartment (DIL, e.g. river water). Flat AEMs and CEMs are
usually separated by net spacers, but the process performance may be enhanced by the use of self-
supporting profiled membranes [8]. In the following, we will denote the three spatial directions as
axial (along the flow), cross-membrane and lateral. The dimensions of the stack in the axial and lateral
directions will be denoted by Lstack and Wstack, respectively. The projected area of the stack will thus
2
Figure 1. Sketch of a reverse electrodialysis stack.
Due to the electrochemical equilibrium, an ion exchange membrane (IEM) immersed between
two solutions at different concentrations is subject to a voltage difference (Donnan exclusion) [9],
and the sum of all the membrane potentials of a stack is the open circuit voltage (OCV). Under closed
circuit conditions, redox reactions arise due to the voltage difference at the electrodes, thus providing
an electrical current to the load. At the same time, a selective ionic transport takes place within the
The theoretical (maximum) electromotive force, or OCV, is proportional to the logarithm of the
ratio between the ion activities in the two solutions. However, the use of highly concentrated solutions
may imply a significant reduction of the IEMs ability to allow the passage of a single ionic species
(permselectivity) [10,11].
The voltage over the external load can be expressed as OCV less the voltage loss due to the
internal resistance of the stack, which takes into account different phenomena.
Ohmic losses (ηΩ) are mainly due to the resistance of the dilute channel when low concentration
solutions, e.g. river water, are used [12–14]. Stacks with net spacers made by non-conductive material
3
may cause an increase of ηΩ by more than 50% [10,12–15], while stacks with profiled membranes
The ion transport across IEMs causes a streamwise concentration variation in the bulk of the
solution and a concentration variation perpendicular to the membranes within the diffusion boundary
layers. These concentration non-uniformities give rise to further voltage drops, referred to as non-
Ohmic losses (η∆c and ηBL) [17]. η∆c is determined by the ion mass balance within the channels. In
stacks fed by river water - seawater solutions, under conditions of maximum net power density, η∆c
Concentration polarisation phenomena cause a change of the Donnan potential with consequent
reduction of the actual voltage over the membrane [18]. With the standard river water – seawater
solutions ηBL may be significant [13,14,16,19]; however, at the flow rates allowing the maximum net
power density it is usually lower than the other contributions to the overall voltage drop when using
spacers [10,14,16], although Vermaas et al. [16] found it to be larger than η∆c and only slightly lower
than ηΩ when using profiled membranes (mixing not favoured). As the feeds concentrations increase,
concentration variations are relatively lower and thus non-Ohmic effects decrease [10,20].
The pumping power consumption may considerably reduce the actual power. Net spacers may
increase by several times hydraulic friction in the channel with respect to the ideal case of a spacer-
less channel [8,20,21], while simple profiled membranes may be preferable [8,16,22]. Several
experimental data showed pressure drop within the manifolds comparable or higher than those
distributed in the channels [8,16,23–25]. However, pressure drop within the manifolds can be
It is clear that RED is a complex process where various phenomena occur at different space-
time scales and interact with one another. In this context, the mathematical description and the
numerical simulation can provide a modelling tool very useful for predictive purposes and for process
4
1.3 State of the art of the small scale simulation
The simulation of fluid dynamics and mass transport phenomena occurring at small scale by
Computational Fluid Dynamics (CFD) can provide important basic data to a process simulator, thus
greatly enhancing its predictive power. In previous papers we performed CFD simulations in order to
investigate spacer-filled and profiled membrane channels for RED applications [20–22,27]. Fluid
flow governing equations (Navier–Stokes and continuity) and a transport equation (diffusion-
convection) for a binary electrolyte (NaCl) were solved by using the 3-D CFD code Ansys CFX™
[28]. Synthetic parameters (dimensionless numbers) quantifying pressure drop and mass transfer,
such as Sherwood numbers and friction coefficients, were calculated. Important insights on the
features that a high performance RED channel should possess were thus obtained.
Pawlowski et al. [29] followed a similar simulation approach by using the open-source Open-
FOAM software [30]. CFD results were then used as input data in a simplified model for the
calculation of the net power density of a RED stack. In particular, the authors proposed a profiled
membrane channel with chevron corrugations, i.e. with staggered herringbone ribs as those used in
micro-mixers.
investigating some specific aspects concerning RED units. For example, Kim et al. [31] studied a
module with channels in serial configuration, Weiner et al. [32] proposed to improve the stack
performance by blending the low salinity feed with a higher salinity stream before the stack entrance,
Moya [33] evaluated the thickness of the diffusion boundary layer and the optimal load resistance.
equation is not sufficient, but ionic migration must be taken into account. The Nernst–Planck equation
is commonly adopted along with the local electroneutrality condition. Several works were carried out
Among these works, Jeong et al. [34] proposed a FEM model of RED implemented in
COMSOL™ [35]. A cell pair of 0.4 m long spacer-less channels was simulated in 2-D (axial + cross-
5
membrane); membranes were not included in the computational domain, but were modelled by
imposing boundary conditions. The maximum net power density was obtained by choosing a channel
Zourmand [36] developed a similar model, based on the COMSOL™ software [35], for ED
devices. In this case, membranes were simulated as domain with Ohmic behaviour, i.e. by neglecting
the effect of concentration variations. Tadimeti et al. [37] extended the model in order to investigate
the effect of obstacles (spacer/corrugations), and validated it against experimental data. Simulation
results showed that, as the number of obstacles per unit length and/or the obstacle height increased,
pressure drop increased. The effect of the geometry on concentration polarization was not
straightforward, because the obstacles generated stagnant zones but also zones with improved
convective motion.
Tedesco et al. [38] conducted 2-D simulations of transport phenomena in ED/RED systems.
The computational domain consisted of two empty half channels and one membrane, i.e. a symmetric
behaviour was assumed for the two ionic species of a 1:1 salt and for the two membranes. The effect
of co-ion transport on the process performance was assessed, and permselectivity and electrical
resistance of membranes were predicted as functions of the fixed charge density and of the position
From the (brief, but sufficiently complete) review conducted in the previous section, it appears
clear that the modelling of electro-membrane processes is a very complex task. The mathematical
models proposed in the literature are based on simplifying assumptions on (i) the geometry of the
6
The aim of this work was to develop a relatively complete modelling tool for RED systems,
able physically to describe the transport phenomena within the cell pair and to compare the
simulates simultaneously fluid dynamics and ionic transport in two dimensions within a periodic
geometrical unit of the cell pair. The only feature of the stack that the model does not simulate is the
stack length with the related axial voltage drop (η∆c). Electrochemical mass transport is simulated by
the Nernst–Planck equation along with the local electroneutrality condition, and double layer
phenomena are simulated by the Donnan exclusion theory. Membranes were explicitly simulated as
parts of the computational domain, and the 2-D spatial distribution of concentration, current density
and electrical potential is predicted within the whole cell pair. Various membrane/channel
geometrical configurations were studied with good spatial resolution. Moreover the influence of
concentration and velocity of feed solutions was assessed. The electrical power delivered by a stack
to an external load was computed, thus giving important information for the optimization of stack
The geometrical domain of a RED cell pair was simulated in 2-D (axial + cross-membrane). As
Needless to say, in a 2-D model obstacles spanning the whole channel thickness cannot be simulated.
Therefore, spacers and profiled membranes were described by non-obstructing shapes which occupy
only a fraction of the thickness of each channel but maintain some of the features of real 3-D
7
Figure 2. Computational domain and boundary conditions.
Electrodic compartments were not explicitly simulated, but their resistance (blank resistance,
Rblank) was taken into account as shown below. By considering the equivalent electrical circuit of the
system (Figure 3), the electrical quantities relevant to a lab-scale stack of n = 10 cell pairs and a
The model was implemented in the FEM commercial software COMSOL Multiphysics® [35].
Rext
n·OCVcp
n·Rcp Rblank
n·Vcp
Vstack
Continuity and momentum transport (Navier-Stokes) equations for a Newtonian fluid are:
8
𝜕𝜌
+ ∇𝜌𝑢⃗ = 0 (1)
𝜕𝑡
𝜕𝜌𝑢⃗
+ 𝑢⃗∇𝜌𝑢⃗ = −∇𝑝 + ∇𝜇∇𝑢⃗ (2)
𝜕𝑡
migration + convection):
𝑁⃗ = −𝐷 ∇𝑐 − 𝑧 𝑢 , 𝐹𝑐 ∇𝜙 + 𝑢⃗𝑐 (3)
where the subscript i indicates the ionic species, Di is the diffusivity, ci is the ion concentration, zi is
the valence, um,i is the mobility, F is Faraday’s constant and ϕ is the electrical potential. The mass
𝜕𝑐
+ ∇𝑁⃗ = 0 (4)
𝜕𝑡
∑𝑧 𝑐 = 0 (5)
In the membrane domains, also the concentration of fixed charges was included in the terms ci.
𝚤⃗ = 𝐹 𝑧 𝑁⃗ (6)
Multiplication of Eq. (4) by ziF and summation over the species taking Eqs. (5) and (6) into
∇𝚤⃗ = 0 (7)
With reference to the equivalent electrical circuit in Figure 3, by linking the distributed values
of current density and electrical potential with the quantities relevant to the external circuit, and
9
adding Ohm’s law for the external circuit, one can obtain the boundary conditions for current density
and electrical potential. The length-average of voltage and current density over one cell pair are
calculated at the boundaries of the computational domain in Figure 2 (arbitrarily located at mid-AEM)
1
〈𝑖〉 = 𝑖 𝑑𝑙 (8)
𝑙
1
〈𝜙〉 = 𝜙 𝑑𝑙 (9)
𝑙
where l denotes the length of the 2-D domain simulated. The voltage over the cell pair is calculated
𝑉 = 〈𝜙 〉 − 〈𝜙 〉 (10)
The average voltage and the average current density over the cell pair are supposed to be uniform
along the flow direction (i.e. the axial voltage drop η∆c is neglected). This is an acceptable assumption
in the present study, since it is aimed at comparing different membrane/channel arrangements rather
than at predicting the stack performance. Thus the electrical current is:
𝐼 = 〈𝑖〉𝑆 (11)
By considering the external circuit, the electrical current and the cell pair voltage are related as
𝑛𝑉
𝐼= (12)
𝑅 +𝑅
where n is the number of cell pairs; Rblank was fixed to 1.0 Ω (from experimental results); and Rext was
treated as a variable input parameter. The voltage over the stack, and thus over the external load, is
𝑉 =𝐼𝑅 (13)
𝑂𝐶𝑉 − 𝑉
𝑅 = (14)
𝐼
10
where OCVcp is the open circuit voltage of the cell pair, computed for Rext → ∞ Ω. It should be kept
in mind that Rcp includes also non-Ohmic resistances and thus depends on the total current I.
The gross power density per cell pair and per unit projected membrane area produced by the
stack is
𝑉 ∙𝐼
𝐺𝑃𝐷 = (15)
𝑛𝑆
∆𝑝 𝑄 + ∆𝑝 𝑄
𝑃𝑃𝐷 = (16)
𝑆
where ∆pSOL and QSOL are the pressure drop and the volume flow rate in a single concentrate (CONC)
∆𝑝 = , ,
(17)
𝑄 =𝑢 𝐻 𝑊 (18)
in which pinflow,SOL and pouflow,SOL are the pressure at the inflow and outflow periodic boundaries,
respectively, l is the length of the periodic domain, uSOL is the superficial velocity of the solution and
HSOL is the channel height (inter-membrane distance). Finally, the net power density per cell pair is
A sketch of the domain simulated was shown in Figure 2 along with the sizes and the boundary
conditions set. The external boundaries were placed in correspondence of the midplane of the AEM.
Periodic conditions were set at the inflow / outflow boundaries, thus simulating fully developed
flow and concentration fields. Appropriate source terms were added at the right side of Eq. (2) and
Eq. (4), taking into account the large-scale gradients of pressure and concentration, respectively.
11
Reducing the computational domain to a small periodic tract of the cell pair allows the spatial
resolution to increase while keeping the total number of elements compatible with an acceptable
memory requirement and computing time (more details on the mesh are given in section 2.7).
However, this implies renouncing to calculating η∆c, but it could be included by coupling the present
model with a higher scale model taking into account mass balances from inlet to outlet. For both
channels, the superficial velocity was made to vary in the range 0.3-5 cm/s. The bulk concentration
(which remains fixed at its initial value in each fluid) was set at 4 M for the concentrate channel and
was made to vary in the range 0.005-0.5 M in the dilute channel. The choice of the concentrations
was made in order to simulate optimal conditions for the power output, with possible applications in
areas where concentrate brines are available or for a closed loop [39,40] where the solutions are
The external boundaries were set as periodic with respect to concentration and current density.
At the IEM-solution interfaces, the Donnan electrochemical equilibrium conditions were set
according to the following formulae [18,41], in which, for simplicity, activity coefficients were set to
1:
𝑅𝑇 𝑐 ,
𝜙 =𝜙 −𝜙 = 𝑙𝑛 (20)
𝑧 𝐹 𝑐 ,
𝑐 , −𝑐 , = 𝑐 + 4𝑐 , 𝑐 , −𝑐 +
(21)
+𝑎𝑐 −𝑐 ,
where 𝜙 is the Donnan potential, and indicate IEM side and solution side values,
respectively, at the IEM-solution interface, R is the universal gas constant, T is temperature, IEM
stands for AEM or CEM, the subscripts co, counter and fix indicate co-ion, counter-ion and fixed
charges, respectively, a is the fraction (0.1%) of fixed charges having the same sign of counter-ions
[41].
12
2.5 Membrane/channel configurations
Different configurations of the cell pair were simulated as shown in Figure 4. In all cases, the
thickness of both concentrate and dilute channels was 270 μm. The side or diameter of the obstacles
was slightly more than one half of the channel thickness (150 µm). Of course, empty channels are
only an ideal case (as those simulated by Jeong et al. [34]), due to the need for mechanically
supporting the membranes. The configurations with obstacles in Figure 4 are necessarily only an
approximation of actual three-dimensional shapes because in 2-D simulations the obstacles cannot
occupy the whole channel height. Similar simplified 2-D shapes have been often simulated in the past
in the context of membrane processes, see e.g. Shakaib et al. [42] or Tadimeti et al. [37]. The distance
between two obstacles on the same side was 600 µm and was chosen as the length of the periodic
computational domain.
Figure 4. Different membrane/channel configurations simulated: (I) flat membranes and empty channels, (II) flat
membranes and non-conductive round spacers placed alternately on both sides of the membranes, (III) flat membranes
and non-conductive square spacers placed alternately on both sides of the membranes, (IV) profiled membranes with
different arrangements.
13
2.6 Materials
All the simulations were run at the temperature of 20 °C. NaCl aqueous solutions were
simulated. Physical properties within each fluid were set as uniform, as computed at the bulk
concentrations. The ions mobility in solution was calculated by the Nernst-Einstein equation
𝐷
𝑢 , = (22)
𝑅𝑇
Membranes were assumed to be homogeneous and isotropic and were simulated as electrolytic
solutions with a null flow field. Diffusivity and mobility within the membranes were obtained by
The concentrations of fixed charges in the membranes were assumed to be uniform and equal
to 4267 mol/m3 for the AEM and 4833 mol/m3 for the CEM on the basis of experimental data on ion
exchange capacity and water uptake (FUJIFILM™ membranes). The diffusion coefficient and the
Hybrid finite element meshes were used (see example in Figure 5). They were composed by
quadrilateral elements near the boundaries and triangular elements elsewhere. Grid independence was
preliminarily addressed, and a mesh with a total number of elements of ~40,000 was chosen. This
mesh exhibits a very good spatial resolution, along with moderate computing times and memory
requirements, while finer grids with an order of magnitude of elements more yield only negligibly
different results, against a very large increase in computing time and, especially, memory
requirements. This feature makes also an accurate simulation of the cell pair from inlet to outlet
prohibitive.
The numerical solution of the equation system was structured in two steps: the first step solved
the fluid dynamics (laminar flow module in COMSOL); the second step solved the electrochemistry
14
(tertiary current distribution module), by using the velocity field from the first step. The fully coupled
Figure 5. FEM discretization of the computational domain in the case of profiled membranes.
A large number of simulations was carried out, but, for the sake of brevity, only the main results
will be shown and discussed in this paper. Results are reported in this section on the distribution of
the main quantities within the cell pair in order to verify the physical consistency of the model
predictions. The effects of the various parameters will be discussed in the following section.
Figure 6 shows maps of the velocity module for various geometries simulated. The superficial
velocity was assumed equal to 1 cm/s in the two channels. At the very low Reynolds numbers typical
of RED channels (< 10) the flow is steady. Within empty channels, the flow is parallel and the velocity
exhibits a parabolic profile. When obstacles are present, velocity components perpendicular to the
membranes arise, i.e. the flow path becomes tortuous. The larger surface causing friction and the
inertial effects cause higher pressure drops; on the other hand, a mixing enhancement is generally
obtained, although stagnant regions in the proximity of the obstacles can negatively affect this aspect.
Finally, some small differences are caused by the obstacle shape (round vs. square), and the presence
of obstacles only on one side of the channel causes larger stagnant regions in the proximity of that
side.
15
Figure 6. Velocity map for different configurations, at a superficial velocity of 1 cm/s in the two channels.
3.2 Concentration
Transport phenomena of ions involve the whole cell pair, occurring both in the solutions and in
the membranes. Figure 7 shows an example of concentration profiles predicted in the case of empty
channels. Double layer phenomena are simulated as sudden jumps of concentration at each IEM-
solution interface. Because of the electroneutrality condition, in the solutions the concentrations of
Na+ and Cl- are equal, while in the membranes they differ by a quantity equal to the concentration of
the fixed charges. The co-ions concentration inside the membranes is not negligible, due to the high
The insets in Figure 7 show also concentration polarization phenomena within the fluid
domains. Under open circuit conditions, since the membranes are not perfectly permselective, a
diffusive transport takes place, thus generating concentration gradients. Under closed circuit
conditions (as in the case shown in Figure 7), a current flows through the stack. In the solutions it is
transported in a similar amount by cations and anions (transport numbers ≈ 0.5), while in the
membranes it is transported almost exclusively by counter-ions (transport number close to 1). The
total flux in the solutions is maintained constant by diffusive fluxes (concentration polarization).
16
AEM CONC CEM DIL AEM
8,000
Na+
7,000 Cl-
fix,IEM
6,000
5,000
ci [mol/m 3]
4,000
520
3,000 4,010
490
2,000 0. 45 0. 6
3,980
0. 19 0. 34
1,000
0
0 0.2 0.4 0.6 0.8
x [mm]
Figure 7. Concentration profiles within the cell pair with empty channels, fed by 4M/0.5M solutions at 1 cm/s in the two
channels. Rext was set equal to OCV/Ishort-circuit. The values are taken along an arbitrary horizontal line.
The Donnan equilibrium for the electrical potential is computed as a sudden jump at each IEM-
solution interface, as shown in Figure 8, which reports the potential profiles for different external
loads. In particular, the IEM-CONC interfaces are associated with a voltage loss, while the IEM-DIL
interfaces are associated with a voltage gain, larger than the IEM-CONC loss.
At open circuit no electrical current circulates, so that there are no voltage losses. The electrical
potential is flat also in the membranes since the diffusion coefficients were set at the same values for
At closed circuit, an electrical current flows through the cell pair and the potential exhibits some
losses. The largest voltage losses occur in short-circuit conditions, i.e. when the maximum current is
circulating. The maximum gross power density is obtained when the stack voltage is approximately
equal to OCV/2, corresponding to an external resistance equal to the internal one (actually, non-
Ohmic phenomena cause a small deviation). The voltage drop is significant within the membranes
and, to a lesser extent, within the 0.5M dilute solution. Conversely, the potential drop within the 4M
17
concentrate solution is negligible, due to the higher electrical conductivity. Note that the large
variation of concentration inside the membrane causes a significant variation of conductivity, thus
-0.01
-0.03
-0.05
-0.07
-0.09
0 0.2 0.4 0.6 0.8
x [mm]
Figure 8. Potential profile within the cell pair with empty channels, fed by 4M/0.5M solutions at 1 cm/s in the two
channels, at different external loads. Rblank was set equal to zero in order to visualize all the internal losses inside the cell
solutions and membranes. In this regard, it is interesting to consider the current density distribution
in complex geometries. Figure 9 shows the current density maps and current lines for a cell pair with
profiled membranes fed by cCONC = 4M and cDIL = 0.01 (a), 0.05 (b) and 0.1M (c). When the profiles
are immersed in a highly conductive solution such as the 4M CONC, they exhibit an almost null
current density and are bypassed by the current lines. Conversely, when the profiles are immersed in
a less conductive solution such as the 0.01M DIL, plot (a), they are crossed by a significant current
density. A concentration of 0.05M, plot (b), is sufficient to confer the dilute solution a conductivity
higher than that of the present membranes, and a concentration of 0.1M, plot (c), makes the current
18
density negligible inside the profiles on the DIL side. Therefore, the features of membranes and
solutions significantly affect the usefulness of the profiles for increasing the active area and reducing
Figure 9. Current density distribution and current lines within the cell pair with profiled membranes, fed by 4M
concentrate solution and various dilute solutions at a superficial velocity of 1 cm/s in the two channels. Rext was set equal
to OCV/Ishort-circuit.
This section will describe a sensitivity analysis carried out by maintaining the concentrate at
4M and letting the dilute concentration and the channel configuration vary. The solution velocity was
also made to vary in order to find the maximum net power density.
In all cases investigated, the Gross Power Density increases asymptotically towards a maximum
as the fluid velocity increases, due to the reduction of the polarization loss ηBL. In particular, the
increase of GPD is larger when the dilute concentration is lower, as expected. By way of example,
Figure 10 reports GPD as a function of the superficial velocity (assumed to be equal in the two
channels) for stacks with different membrane/channel configurations, fed by 0.01M/4M (left) and
0.05M/4M (rigth) solutions. For 0.01M dilute solution, left plot, the highest GPD values are obtained
with the AEM/CEM-DIL configuration of Figure 4. The same finding emerges also for cDIL = 0.005M
19
(not shown for brevity). On the contrary, in the case of cDIL = 0.05M, right plot, and higher the
maximum GPD values were provided by stacks with empty channels. This behaviour highlights that
the use of profiles reduces the overall stack Ohmic resistance only in the presence of a dilute solution
0.01M/4M 0.05M/4M
5.0 5.0
4.5 4.5
4.0 4.0
GPD [W/m 2]
GPD [W/m 2]
3.5 3.5
3.0 3.0
Empty Round Sp ace rs Empty Round Space rs
Square Spacers Pro filed Membran es Squ are Spacers Pro filed Membran es
2.5 Pro filed AE M- DIL Pro filed AE M- CO NC 2.5 Pro filed AE M- DIL Pro filed AE M- CO NC
Pro filed CEM-DIL Pro filed CEM-CONC Pro filed CEM-DIL Pro filed CEM-CONC
Pro filed AE M/CEM-DIL Pro filed AE M/CEM-DIL
2.0 2.0
0 1 2 3 4 5 0 1 2 3 4 5
uSOL [cm/s] uSOL [cm/s]
Figure 10. Gross Power Density for stacks (10 cell pairs) with different membrane/channel configurations, fed by
0.01M/4M (left) and 0.05M/4M (rigth) solutions, as a function of the superficial velocity (assumed to be equal in the two
The channel geometry significantly affects the pressure drop. Figure 11 reports the sum of the
pressure losses in the concentrate and in the dilute channels per unit stack length for various
configurations as a function of the superficial velocity (assumed to be equal in the two channels). One
can observe that the presence of obstacles may significantly increase pressure drop. In particular,
square obstacles on both sides of both membranes cause the highest values of pressure drop. Round
spacers produce a lower hydraulic friction. As the overall number of obstacles per unit length
decreases, pressure drop decreases; for example, moving from the AEM/CEM-DIL configuration to
the AEM-DIL or CEM-DIL configurations, one can observe an almost twofold reduction of pressure
20
drop. Notably, the pressure drop is higher when profiles are present inside the concentrate channel,
100
80
60
40
20
0
0 1 2 3 4 5
uSOL [cm/s]
Figure 11. Pressure drop for different membrane/channel configurations fed by 0.5M/4M solutions, as a function of the
Figure 12 reports the Net Power Density for the same cases as in Figure 10. All curves exhibit
a maximum as a function of the fluid velocity. As for the case of GPD, at cDIL ≤ 0.01M the
AEM/CEM-DIL configuration yields the highest NPD values (~4.4 W/m2 for cDIL = 0.01M), while at
cDIL ≥ 0.05M it is the empty channel configuration that yields the highest NPD values (~4.1 W/m2 for
cDIL = 0.05M).
21
0.01M/4M 0.05M/4M
5.0 5.0
4.5 4.5
4.0 4.0
NPD [W/m2]
NPD [W/m 2]
3.5 3.5
3.0 3.0
Empty Round Sp ace rs Empty Round Sp ace rs
Squ are Spacers Pro filed Membran es Squ are Spacers Pro filed Membran es
2.5 Pro filed AE M- DIL Pro filed AE M- CO NC 2.5 Profiled AE M- DIL Pro filed AE M- CO NC
Pro filed CEM-DIL Pro filed CEM-CONC Profiled CEM-DIL Pro filed CEM-CONC
Pro filed AE M/CEM-DIL Profiled AE M/CEM-DIL
2.0 2.0
0 1 2 3 4 5 0 1 2 3 4 5
uSOL [cm/s] uSOL [cm/s]
Figure 12. Net Power Density for stacks (10 cell pairs) with different membrane/channel configurations, fed by 0.01M/4M
(left) and 0.05M/4M (rigth) solutions, as a function of the superficial velocity (assumed to be equal in the two channels).
This section compares the results obtained for different dilute concentrations and different
geometries at the fluid velocity that, for each configuration, maximizes the Net Power Density. In all
The values of the cell pair electric resistance under conditions of maximum NPD are shown in
Figure 13. This is the total resistance, due to both Ohmic and non-Ohmic phenomena. In order to give
an idea of the significance of boundary layer effects, we mention that by using an empty channel 270
μm thick for the dilute solution and the couple 4M/0.01M, the non-Ohmic resistance associated to
concentration polarization is ~15% of Rcp. The presence of obstacles within the channel enhances
mixing, thus reducing the non-Ohmic resistance. However, the influence of the obstacles on the
Ohmic resistance may be larger; therefore, non-conductive spacers cause larger resistances, especially
square spacers which occupy a larger volume with respect to round spacers. On the contrary, profiles
more conductive than the solution in which they are immersed, e.g. cDIL ≤ 0.01M, reduce Rcp, while
lower boundary layer effect [20]. With a dilute concentration of 0.005M, Rcp is about an order of
Figure 13. Cell pair resistance for different dilute concentrations and membrane/channel configurations, at the superficial
velocities (assumed to be equal in the two channels) that maximize the net power density, and with 4M concentrate
The Gross Power Density GPD (see Figure 14) depends on Rcp and also on OCV, which, in its
turn, is strongly affected by the feeds concentrations. As a consequence, the trend of GPD with cDIL
exhibits a maximum. This is obtained at cDIL = 0.01M-0.05M, depending on the geometry. Most
profiled membranes enhance the process performance with respect to the empty channels for cDIL ≤
0.01M. The differences in GPD among the various membrane/channel configurations decrease as cDIL
increases. Note that the present simulations are for a fully-developed tract of the cell pair and do not
take axial concentration changes into account. If the salt intake into the dilute channel from the
concentrate channel were considered, the average concentration of the diluate would be higher than
the inlet value, especially in the case of long channels. Therefore, the actual optimal value of the inlet
concentration for the dilute stream might be lower than 0.01M. For the same reason, (ideal) stacks
23
with empty channels might provide GPD values higher than stacks with profiled membranes even at
Figure 14. Gross Power Density for different dilute concentrations and membrane/channel configurations, at the
superficial velocities (assumed to be equal in the two channels) that maximize the net power density, and with 4M
The features of GPD are reflected into NPD, which is also influenced by the pumping power
and, thus, again, by the geometry. Figure 15 reports the maximum Net Power Densities obtained in
all the cases simulated. As in the case of GPD, the highest values of NPD are obtained for cDIL =
0.01M-0.05M. Stacks with spacers are characterized by higher Rcp and pressure drops, thus leading
to lower NPD, especially in the case of the square shape. Profiled membranes increase NPD with
respect to the empty channels for cDIL ≤ 0.01M, yielding the highest NPD ≈ 4.4 W/m2 with the
AEM/CEM-DIL configuration. Stacks with empty channels yield higher NPD than other
24
Figure 15. Maximum Net Power Density for different dilute concentrations and membrane/channel configurations, with
5 Conclusions
A multi-physical model of RED units was developed. Fluid dynamics, electrochemical mass
transport and double layer phenomena were modelled in a two-dimensional periodic domain of a cell
pair by the Navier–Stokes equations complemented by the Nernst–Planck approach and the Donnan
exclusion theory. The model is fully predictive and requires only few empirical data as input
parameters.
A sensitivity analysis of the stack performance to the stack features and operating conditions,
such as membrane/channel configuration, solution concentrations and feed velocities, was conducted.
With respect to the ideal case of empty (i.e. spacer-less) channels, cell pairs with non-conductive
spacers were characterized by higher electrical resistance and higher pressure drop, and thus by lower
net electrical power. The membranes considered in this study were more conductive than solutions
with a concentration up to 0.01M, so that profiled membranes reduced the cell electrical resistance
with respect to empty channels for cDIL less than this value. This resulted in an increase of gross and
25
net power density at cDIL ≤ 0.01M, despite the increment in pressure drop. Therefore, the effectiveness
of profiled membranes for the active area increase and the stack resistance reduction depends on the
stack features and operating conditions, and only in some cases the profiles can actually enhance the
stack performance. Of course, manufacturing highly conductive membranes would be important for
the process optimization. On the other hand, empty channels are only an ideal configuration and
cannot be used for real applications, due to the need of a mechanical support for the membranes.
The multi-physical model presented possesses great capabilities for treating the complex
phenomenology of the RED process. The model is a valid tool in order to understand the behaviour
of RED systems and to investigate the effect of the various operating parameters and stack features.
In order to develop a more effective predicting tool, the model should be extended to 3-D geometries
and coupled with a higher scale simulation taking into account mass balance from inlet to outlet.
Acknowledgements
This work has been performed within the RED-Heat-to-Power (Conversion of Low Grade Heat
to Power through closed loop Reverse Electro-Dialysis) project – Horizon 2020 programme, Grant
Nomenclature
GPD Gross power density of the stack per membrane area [W m-2]
NPD Net power density of the stack per membrane area [W m-2]
PPD Pumping power density of the stack per membrane area [W m-2]
p Pressure [Pa]
t Time [s]
𝑢⃗ Velocity [m s-1]
Greek letters
η∆c Voltage loss due to the streamwise concentration change in the bulk of the solutions
[V]
ηBL Voltage loss due to the concentration polarization in the boundary layers [V]
Subscripts/superscripts
CONC Concentrate
27
DIL Dilute
References
[1] M. Turek, B. Bandura, Renewable energy by reverse electrodialysis, Desalination. 205 (2007)
67–74. doi:10.1016/j.desal.2006.04.041.
[2] D. Brogioli, Extracting renewable energy from a salinity difference using a capacitor., Phys.
[3] A. Achilli, A.E. Childress, Pressure retarded osmosis: From the vision of Sidney Loeb to the
doi:10.1016/j.desal.2010.06.017.
[4] F. La Mantia, M. Pasta, H.D. Deshazer, B.E. Logan, Y. Cui, Batteries for efficient energy
doi:10.1021/nl200500s.
[5] R.E. Pattle, Production of Electric Power by mixing Fresh and Salt Water in the Hydroelectric
[6] A. Cipollina, G. Micale, eds., Sustainable Energy from Salinity Gradients, 1st ed., Woodhead
[7] E. Brauns, Salinity gradient power by reverse electrodialysis: effect of model parameters on
[8] D.A. Vermaas, M. Saakes, K. Nijmeijer, Enhanced mixing in the diffusive boundary layer for
doi:10.1016/j.memsci.2013.11.005.
[9] H. Strathmann, Ion-exchange membrane separation processes, first ed., Elsevier, Amsterdam,
28
2004.
[10] A. Daniilidis, D.A. Vermaas, R. Herber, K. Nijmeijer, Experimentally obtainable energy from
mixing river water, seawater or brines with reverse electrodialysis, Renew. Energy. 64 (2014)
123–131. doi:10.1016/j.renene.2013.11.001.
[11] M. Tedesco, A. Cipollina, A. Tamburini, I.D.L. Bogle, G. Micale, A simulation tool for
analysis and design of reverse electrodialysis using concentrated brines, Chem. Eng. Res. Des.
[12] J.W. Post, H.V.M. Hamelers, C.J.N. Buisman, Energy recovery from controlled mixing salt
and fresh water with a reverse electrodialysis system, Environ. Sci. Technol. 42 (2008) 5785–
5790. doi:10.1021/es8004317.
electrodialysis as process for sustainable energy generation, Environ. Sci. Technol. 43 (2009)
6888–6894. doi:10.1021/es9009635.
[14] D.A. Vermaas, M. Saakes, K. Nijmeijer, Doubled power density from salinity gradients at
doi:10.1021/es2012758.
[15] S. Pawlowski, P. Sistat, J.G. Crespo, S. Velizarov, Mass Transfer in Reverse Electrodialysis:
Flow Entrance Effects and Diffusion Boundary Layer Thickness, J. Memb. Sci. 471 (2014)
72–83. doi:10.1016/j.memsci.2014.07.075.
[16] D.A. Vermaas, M. Saakes, K. Nijmeijer, Power generation using profiled membranes in
doi:10.1016/j.memsci.2011.09.043.
[17] D.A. Vermaas, E. Guler, M. Saakes, K. Nijmeijer, Theoretical power density from salinity
doi:10.1016/j.egypro.2012.03.018.
[18] K. Kontturi, M. Lasse, J.A. Manzanares, Ionic Transport Processes in Electrochemistry and
29
Membrane Science, Oxford University Press Inc., New York, 2008.
doi:10.1093/acprof:oso/9780199533817.001.0001.
resistances of membrane, diffusion boundary layer and double layer in ion exchange membrane
polarization phenomena in spacer-filled channels for reverse electrodialysis, J. Memb. Sci. 468
[21] L. Gurreri, A. Tamburini, A. Cipollina, G. Micale, M. Ciofalo, Flow and mass transfer in
spacer-filled channels for reverse electrodialysis: a CFD parametrical study, J. Memb. Sci. 497
[22] L. Gurreri, M. Ciofalo, A. Cipollina, A. Tamburini, W. Van Baak, G. Micale, CFD modelling
20. doi:10.1080/19443994.2014.940651.
[23] J. Veerman, M. Saakes, S.J. Metz, G.J. Harmsen, Reverse electrodialysis: Performance of a
stack with 50 cells on the mixing of sea and river water, J. Memb. Sci. 327 (2009) 136–144.
doi:10.1016/j.memsci.2008.11.015.
[24] J. Veerman, M. Saakes, S.J. Metz, G.J. Harmsen, Electrical power from sea and river water by
reverse electrodialysis: A first step from the laboratory to a real power plant, Environ. Sci.
Experimental and modeling studies for stacks with variable number of cell pairs, J. Memb. Sci.
[26] L. Gurreri, A. Tamburini, A. Cipollina, G. Micale, CFD analysis of the fluid flow behavior in
doi:10.1080/19443994.2012.705966.
30
[27] L. Gurreri, A. Tamburini, A. Cipollina, G. Micale, M. Ciofalo, Pressure drop in woven-spacer-
filled channels for reverse electrodialysis: CFD prediction and experimental validation,
[29] S. Pawlowski, V. Geraldes, J.G. Crespo, S. Velizarov, Computational fluid dynamics (CFD)
[31] K.S. Kim, W. Ryoo, M.S. Chun, G.Y. Chung, Simulation of enhanced power generation by
reverse electrodialysis stack module in serial configuration, Desalination. 318 (2013) 79–87.
doi:doi.org/10.1016/j.desal.2013.03.023.
[32] A.M. Weiner, R.K. McGovern, J.H. Lienhard V, Increasing the power density and reducing
[33] A.A. Moya, A numerical comparison of optimal load and internal resistances in ion-exchange
membrane systems under reverse electrodialysis conditions, Desalination. 392 (2016) 25–33.
doi:10.1016/j.desal.2016.04.016.
[34] H.-I. Jeong, H.J. Kim, D.-K. Kim, Numerical analysis of transport phenomena in reverse
doi:10.1016/j.energy.2014.03.013.
doi:10.1016/j.desal.2014.12.008.
for effective ion transport in electrodialysis, J. Memb. Sci. 499 (2016) 418–428.
31
doi:10.1016/j.memsci.2015.11.001.
[38] M. Tedesco, H.V.M. Hamelers, P.M. Biesheuvel, Nernst-Planck transport theory for (reverse)
electrodialysis: I. Effect of co-ion transport through the membranes, J. Memb. Sci. 510 (2016)
370–381. doi:10.1016/j.memsci.2016.03.012.
[39] M. Reali, Closed cycle osmotic power plants for electric power production, Energy. 5 (1980)
325–329. doi:10.1016/0360-5442(80)90033-X.
[40] W. Huang, W.S. Walker, Y. Kim, Junction potentials in thermolytic reverse electrodialysis,
[41] A.H. Galama, J.W. Post, M.A. Cohen Stuart, P.M. Biesheuvel, Validity of the Boltzmann
[42] M. Shakaib, S.M.F. Hasani, I. Ahmed, R.M. Yunus, A CFD study on the effect of spacer
32